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           36 :     Quadrature<T, NumNodes1D, ElementType>::getWeightsForRefElement() const {
      20           36 :         Vector<T, NumNodes1D> w = this->getWeights1D(0.0, 1.0);
      21              : 
      22           36 :         Vector<T, std::remove_reference_t<decltype(*this)>::numElementNodes> tensor_prod_w;
      23              : 
      24           36 :         Vector<unsigned, ElementType::dim> nd_index(0);
      25          470 :         for (unsigned i = 0; i < std::remove_reference_t<decltype(*this)>::numElementNodes; ++i) {
      26          434 :             tensor_prod_w[i] = 1.0;
      27         1602 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      28         1168 :                 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          604 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      35          568 :                 if (++nd_index[d] < NumNodes1D)
      36          398 :                     break;
      37          170 :                 nd_index[d] = 0;
      38              :             }
      39              :         }
      40              : 
      41           36 :         return tensor_prod_w;
      42           36 :     }
      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           54 :     Quadrature<T, NumNodes1D, ElementType>::getIntegrationNodesForRefElement() const {
      48           54 :         Vector<T, NumNodes1D> q = this->getIntegrationNodes1D(0.0, 1.0);
      49              : 
      50           54 :         Vector<Vector<T, ElementType::dim>, std::remove_reference_t<decltype(*this)>::numElementNodes> tensor_prod_q;
      51              : 
      52           54 :         Vector<unsigned, ElementType::dim> nd_index(0);
      53          600 :         for (unsigned i = 0; i < std::remove_reference_t<decltype(*this)>::numElementNodes; ++i) {
      54         1998 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      55         1452 :                 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          786 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      62          732 :                 if (++nd_index[d] < NumNodes1D)
      63          492 :                     break;
      64          240 :                 nd_index[d] = 0;
      65              :             }
      66              :         }
      67              : 
      68           54 :         return tensor_prod_q;
      69           54 :     }
      70              : 
      71              :     template <typename T, unsigned NumNodes1D, typename ElementType>
      72          118 :     Vector<T, NumNodes1D> Quadrature<T, NumNodes1D, ElementType>::getIntegrationNodes1D(
      73              :         const T& a, const T& b) const {
      74          118 :         assert(b > a);
      75              :         // scale the integration nodes from the local domain [a_m, b_m] to the given one [a, b]
      76              : 
      77          118 :         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           82 :     Vector<T, NumNodes1D> Quadrature<T, NumNodes1D, ElementType>::getWeights1D(const T& a,
      82              :                                                                                const T& b) const {
      83           82 :         assert(b > a);
      84           82 :         return this->weights_m * (b - a) / (this->b_m - this->a_m);
      85              :     }
      86              : 
      87              : }  // namespace ippl