Branch data     Line data    Source code

       1                 :             : 
       2                 :             : namespace ippl {
       3                 :             :     template <typename T, unsigned NumNodes1D, typename ElementType>
       4                 :         272 :     Quadrature<T, NumNodes1D, ElementType>::Quadrature(const ElementType& ref_element)
       5         [ +  - ]:         272 :         : ref_element_m(ref_element) {}
       6                 :             : 
       7                 :             :     template <typename T, unsigned NumNodes1D, typename ElementType>
       8                 :           0 :     size_t Quadrature<T, NumNodes1D, ElementType>::getOrder() const {
       9                 :           0 :         return this->degree_m + 1;
      10                 :             :     }
      11                 :             : 
      12                 :             :     template <typename T, unsigned NumNodes1D, typename ElementType>
      13                 :             :     size_t Quadrature<T, NumNodes1D, ElementType>::getDegree() const {
      14                 :             :         return this->degree_m;
      15                 :             :     }
      16                 :             : 
      17                 :             :     template <typename T, unsigned NumNodes1D, typename ElementType>
      18                 :             :     Vector<T, Quadrature<T, NumNodes1D, ElementType>::numElementNodes>
      19                 :          30 :     Quadrature<T, NumNodes1D, ElementType>::getWeightsForRefElement() const {
      20         [ +  - ]:          30 :         Vector<T, NumNodes1D> w = this->getWeights1D(0.0, 1.0);
      21                 :             : 
      22         [ +  - ]:          30 :         Vector<T, this->numElementNodes> tensor_prod_w;
      23                 :             : 
      24         [ +  - ]:          30 :         Vector<unsigned, ElementType::dim> nd_index(0);
      25         [ +  + ]:         162 :         for (unsigned i = 0; i < this->numElementNodes; ++i) {
      26                 :         132 :             tensor_prod_w[i] = 1.0;
      27         [ +  + ]:         436 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      28                 :         304 :                 tensor_prod_w[i] *= w[nd_index[d]];
      29                 :             :             }
      30                 :             : 
      31                 :             :             // Update nd_index for next iteration
      32                 :             :             // Increment the nd_index variable in the first dimension, or if it
      33                 :             :             // is already at the maximum value reset it and, go to the higher dimension
      34         [ +  + ]:         214 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      35         [ +  + ]:         184 :                 if (++nd_index[d] < NumNodes1D)
      36                 :         102 :                     break;
      37                 :          82 :                 nd_index[d] = 0;
      38                 :             :             }
      39                 :             :         }
      40                 :             : 
      41                 :          30 :         return tensor_prod_w;
      42                 :          30 :     }
      43                 :             : 
      44                 :             :     template <typename T, unsigned NumNodes1D, typename ElementType>
      45                 :             :     Vector<Vector<T, Quadrature<T, NumNodes1D, ElementType>::dim>,
      46                 :             :            Quadrature<T, NumNodes1D, ElementType>::numElementNodes>
      47                 :          48 :     Quadrature<T, NumNodes1D, ElementType>::getIntegrationNodesForRefElement() const {
      48         [ +  - ]:          48 :         Vector<T, NumNodes1D> q = this->getIntegrationNodes1D(0.0, 1.0);
      49                 :             : 
      50         [ +  - ]:          48 :         Vector<Vector<T, ElementType::dim>, this->numElementNodes> tensor_prod_q;
      51                 :             : 
      52         [ +  - ]:          48 :         Vector<unsigned, ElementType::dim> nd_index(0);
      53         [ +  + ]:         292 :         for (unsigned i = 0; i < this->numElementNodes; ++i) {
      54         [ +  + ]:         832 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      55                 :         588 :                 tensor_prod_q[i][d] = q[nd_index[d]];
      56                 :             :             }
      57                 :             : 
      58                 :             :             // Update nd_index for next iteration
      59                 :             :             // Increment the nd_index variable in the first dimension, or if it
      60                 :             :             // is already at the maximum value reset it and, go to the higher dimension
      61         [ +  + ]:         396 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      62         [ +  + ]:         348 :                 if (++nd_index[d] < NumNodes1D)
      63                 :         196 :                     break;
      64                 :         152 :                 nd_index[d] = 0;
      65                 :             :             }
      66                 :             :         }
      67                 :             : 
      68                 :          48 :         return tensor_prod_q;
      69                 :          48 :     }
      70                 :             : 
      71                 :             :     template <typename T, unsigned NumNodes1D, typename ElementType>
      72                 :         112 :     Vector<T, NumNodes1D> Quadrature<T, NumNodes1D, ElementType>::getIntegrationNodes1D(
      73                 :             :         const T& a, const T& b) const {
      74         [ -  + ]:         112 :         assert(b > a);
      75                 :             :         // scale the integration nodes from the local domain [a_m, b_m] to the given one [a, b]
      76                 :             : 
      77   [ #  #  #  #  :         112 :         return (this->integration_nodes_m - this->a_m) / (this->b_m - this->a_m) * (b - a) + a;
             #  #  #  # ]
           [ +  -  +  -  
          +  -  +  -  +  
                      - ]
      78                 :             :     }
      79                 :             : 
      80                 :             :     template <typename T, unsigned NumNodes1D, typename ElementType>
      81                 :          76 :     Vector<T, NumNodes1D> Quadrature<T, NumNodes1D, ElementType>::getWeights1D(const T& a,
      82                 :             :                                                                                const T& b) const {
      83         [ -  + ]:          76 :         assert(b > a);
      84   [ #  #  #  # ]:          76 :         return this->weights_m * (b - a) / (this->b_m - this->a_m);
           [ +  -  +  -  
                   +  - ]
      85                 :             :     }
      86                 :             : 
      87                 :             : }  // namespace ippl