Line data    Source code

       1              : 
       2              : namespace ippl {
       3              :     template <typename T, unsigned NumNodes1D, typename ElementType>
       4          544 :     Quadrature<T, NumNodes1D, ElementType>::Quadrature(const ElementType& ref_element)
       5          544 :         : 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           64 :     Quadrature<T, NumNodes1D, ElementType>::getWeightsForRefElement() const {
      20           64 :         Vector<T, NumNodes1D> w = this->getWeights1D(0.0, 1.0);
      21              : 
      22           64 :         Vector<T, std::remove_reference_t<decltype(*this)>::numElementNodes> tensor_prod_w;
      23              : 
      24           64 :         Vector<unsigned, ElementType::dim> nd_index(0);
      25          924 :         for (unsigned i = 0; i < std::remove_reference_t<decltype(*this)>::numElementNodes; ++i) {
      26          860 :             tensor_prod_w[i] = 1.0;
      27         3180 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      28         2320 :                 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         1184 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      35         1120 :                 if (++nd_index[d] < NumNodes1D)
      36          796 :                     break;
      37          324 :                 nd_index[d] = 0;
      38              :             }
      39              :         }
      40              : 
      41           64 :         return tensor_prod_w;
      42           64 :     }
      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          100 :     Quadrature<T, NumNodes1D, ElementType>::getIntegrationNodesForRefElement() const {
      48          100 :         Vector<T, NumNodes1D> q = this->getIntegrationNodes1D(0.0, 1.0);
      49              : 
      50          100 :         Vector<Vector<T, ElementType::dim>, std::remove_reference_t<decltype(*this)>::numElementNodes> tensor_prod_q;
      51              : 
      52          100 :         Vector<unsigned, ElementType::dim> nd_index(0);
      53         1184 :         for (unsigned i = 0; i < std::remove_reference_t<decltype(*this)>::numElementNodes; ++i) {
      54         3972 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      55         2888 :                 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         1548 :             for (unsigned d = 0; d < ElementType::dim; ++d) {
      62         1448 :                 if (++nd_index[d] < NumNodes1D)
      63          984 :                     break;
      64          464 :                 nd_index[d] = 0;
      65              :             }
      66              :         }
      67              : 
      68          100 :         return tensor_prod_q;
      69          100 :     }
      70              : 
      71              :     template <typename T, unsigned NumNodes1D, typename ElementType>
      72          228 :     Vector<T, NumNodes1D> Quadrature<T, NumNodes1D, ElementType>::getIntegrationNodes1D(
      73              :         const T& a, const T& b) const {
      74          228 :         assert(b > a);
      75              :         // scale the integration nodes from the local domain [a_m, b_m] to the given one [a, b]
      76              : 
      77          228 :         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          156 :     Vector<T, NumNodes1D> Quadrature<T, NumNodes1D, ElementType>::getWeights1D(const T& a,
      82              :                                                                                const T& b) const {
      83          156 :         assert(b > a);
      84          156 :         return this->weights_m * (b - a) / (this->b_m - this->a_m);
      85              :     }
      86              : 
      87              : }  // namespace ippl