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