Line data Source code
1 : // Class GaussJacobiQuadrature
2 : // The GaussJacobiQuadrature class. This is a class representing a Gauss-Jacobi quadrature
3 : // The algoirthm in computeNodesAndWeights is based on the LehrFEM++ library.
4 : // https://craffael.github.io/lehrfempp/gauss__quadrature_8cc_source.html
5 :
6 : #ifndef IPPL_GAUSSJACOBIQUADRATURE_H
7 : #define IPPL_GAUSSJACOBIQUADRATURE_H
8 :
9 : #include "FEM/Quadrature/Quadrature.h"
10 :
11 : namespace ippl {
12 :
13 : enum InitialGuessType {
14 : Chebyshev,
15 : LehrFEM,
16 : };
17 :
18 : /**
19 : * @brief This is class represents the Gauss-Jacobi quadrature rule
20 : * on a reference element.
21 : *
22 : * @tparam T floating point number type of the quadrature nodes and weights
23 : * @tparam NumNodes1D number of quadrature nodes for one dimension
24 : * @tparam ElementType element type for which the quadrature rule is defined
25 : */
26 : template <typename T, unsigned NumNodes1D, typename ElementType>
27 : class GaussJacobiQuadrature : public Quadrature<T, NumNodes1D, ElementType> {
28 : public:
29 : // a higher precision floating point number type used for the computation
30 : // of the quadrature nodes and weights
31 : using scalar_t = long double; // might be equivalant to double, depending on compiler
32 :
33 : /**
34 : * @brief Construct a new Gauss Jacobi Quadrature rule object
35 : *
36 : * @param ref_element reference element to compute the quadrature nodes on
37 : * @param alpha first Jacobi parameter alpha
38 : * @param beta second Jacobi parameter beta
39 : * @param max_newton_itersations maximum number of Newton iterations (default 10)
40 : * @param min_newton_iterations minimum number of Newton iterations (default 1)
41 : */
42 : GaussJacobiQuadrature(const ElementType& ref_element, const T& alpha, const T& beta,
43 0 : const size_t& max_newton_itersations = 10,
44 0 : const size_t& min_newton_iterations = 1);
45 :
46 : /**
47 : * Computes the quadrature nodes and weights and stores them in the
48 : * quadrature nodes and weights arrays.
49 : */
50 : void computeNodesAndWeights() override;
51 :
52 : /**
53 : * @brief Returns the i-th Chebyshev node, used as initial guess for the Newton iterations.
54 : *
55 : * @param i index of the Chebyshev node
56 : *
57 : * @return scalar_t - i-th Chebyshev node
58 : */
59 : scalar_t getChebyshevNodes(const size_t& i) const; // FIXME maybe move somewhere else?
60 :
61 : private:
62 : /**
63 : * @brief Computes the initial guess for the Newton iterations, the way they are computed in
64 : * the implementation from LehrFEM++.
65 : *
66 : * @param i index of the initial guess (corresponding to the i-th quadrature node)
67 : * @param integration_nodes the integration nodes
68 : *
69 : * @return scalar_t - initial guess
70 : */
71 : scalar_t getLehrFEMInitialGuess(
72 : const size_t& i, const Vector<scalar_t, NumNodes1D>& integration_nodes) const;
73 :
74 : const T alpha_m;
75 : const T beta_m;
76 :
77 : const size_t max_newton_iterations_m;
78 : const size_t min_newton_iterations_m;
79 : };
80 :
81 : /**
82 : * @brief This is class represents the Gauss-Legendre quadrature rule.
83 : * It is a special case of the Gauss-Jacobi quadrature rule with alpha = beta = 0.0.
84 : *
85 : * @tparam T floating point number type of the quadrature nodes and weights
86 : * @tparam NumNodes1D number of quadrature nodes for one dimension
87 : * @tparam ElementType element type for which the quadrature rule is defined
88 : */
89 : template <typename T, unsigned NumNodes1D, typename ElementType>
90 : class GaussLegendreQuadrature : public GaussJacobiQuadrature<T, NumNodes1D, ElementType> {
91 : public:
92 : /**
93 : * @brief Construct a new Gauss Legendre Quadrature rule object
94 : *
95 : * @param ref_element reference element to compute the quadrature nodes on
96 : * @param max_newton_itersations maximum number of Newton iterations (default 10)
97 : * @param min_newton_iterations minimum number of Newton iterations (default 1)
98 : */
99 236 : GaussLegendreQuadrature(const ElementType& ref_element,
100 0 : const size_t& max_newton_itersations = 10,
101 180 : const size_t& min_newton_iterations = 1)
102 : : GaussJacobiQuadrature<T, NumNodes1D, ElementType>(
103 236 : ref_element, 0.0, 0.0, max_newton_itersations, min_newton_iterations) {}
104 : };
105 :
106 : /**
107 : * @brief This is class represents the Chebyshev-Gauss quadrature rule.
108 : * It is a special case of the Gauss-Jacobi quadrature rule with alpha = beta = -0.5.
109 : *
110 : * @tparam T floating point number type of the quadrature nodes and weights
111 : * @tparam NumNodes1D number of quadrature nodes for one dimension
112 : * @tparam ElementType element type for which the quadrature rule is defined
113 : */
114 : template <typename T, unsigned NumNodes1D, typename ElementType>
115 : class ChebyshevGaussQuadrature : public GaussJacobiQuadrature<T, NumNodes1D, ElementType> {
116 : public:
117 : /**
118 : * @brief Construct a new Chebyshev Gauss Quadrature rule object
119 : *
120 : * @param ref_element reference element to compute the quadrature nodes on
121 : * @param max_newton_itersations maximum number of Newton iterations (default 10)
122 : * @param min_newton_iterations minimum number of Newton iterations (default 1)
123 : */
124 56 : ChebyshevGaussQuadrature(const ElementType& ref_element,
125 : const size_t& max_newton_itersations = 10,
126 : const size_t& min_newton_iterations = 1)
127 : : GaussJacobiQuadrature<T, NumNodes1D, ElementType>(
128 56 : ref_element, -0.5, -0.5, max_newton_itersations, min_newton_iterations) {}
129 : };
130 :
131 : } // namespace ippl
132 :
133 : #include "FEM/Quadrature/GaussJacobiQuadrature.hpp"
134 :
135 : #endif
|
