Branch data 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 : 118 : GaussLegendreQuadrature(const ElementType& ref_element,
100 : 0 : const size_t& max_newton_itersations = 10,
101 [ + - + - : 90 : const size_t& min_newton_iterations = 1)
+ - + - +
- + - + -
+ - ][ + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - ]
102 : : : GaussJacobiQuadrature<T, NumNodes1D, ElementType>(
103 [ + - ]: 118 : 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 : 28 : 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 [ + - ]: 28 : 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
|
