LehrFEM++ 1.0.0
A simple Finite Element Library for teaching
make_quad_rule.cc
1
9#include "make_quad_rule.h"
10
11#include <unsupported/Eigen/KroneckerProduct>
12
13#include "gauss_quadrature.h"
14
15namespace lf::quad {
16namespace detail {
17template <base::RefElType REF_EL, int Degree>
18QuadRule HardcodedQuadRule();
19}
20
21QuadRule make_QuadRule(base::RefEl ref_el, unsigned degree) {
22 if (ref_el == base::RefEl::kSegment()) {
23 quadDegree_t n = degree / 2 + 1;
24 auto [points, weights] = GaussLegendre(n);
25 return QuadRule(base::RefEl::kSegment(), points.transpose(),
26 std::move(weights), 2 * n - 1);
27 }
28 if (ref_el == base::RefEl::kQuad()) {
29 quadDegree_t n = degree / 2 + 1;
30 auto [points1d, weights1d] = GaussLegendre(n);
31 Eigen::MatrixXd points2d(2, n * n);
32 points2d.row(0) = Eigen::kroneckerProduct(points1d.transpose(),
33 Eigen::MatrixXd::Ones(1, n));
34 points2d.row(1) = points1d.transpose().replicate(1, n);
35 return QuadRule(base::RefEl::kQuad(), std::move(points2d),
36 Eigen::kroneckerProduct(weights1d, weights1d), 2 * n - 1);
37 }
38 if (ref_el == base::RefEl::kTria()) {
39 switch (degree) {
40 case 1:
41 return detail::HardcodedQuadRule<base::RefEl::kTria(), 1>();
42 case 2:
43 return detail::HardcodedQuadRule<base::RefEl::kTria(), 2>();
44 case 3: // use degree 4 rule instead
45 case 4:
46 return detail::HardcodedQuadRule<base::RefEl::kTria(), 4>();
47 case 5:
48 return detail::HardcodedQuadRule<base::RefEl::kTria(), 5>();
49 case 6:
50 return detail::HardcodedQuadRule<base::RefEl::kTria(), 6>();
51 case 7:
52 return detail::HardcodedQuadRule<base::RefEl::kTria(), 7>();
53 case 8:
54 return detail::HardcodedQuadRule<base::RefEl::kTria(), 8>();
55 case 9:
56 return detail::HardcodedQuadRule<base::RefEl::kTria(), 9>();
57 case 10:
58 return detail::HardcodedQuadRule<base::RefEl::kTria(), 10>();
59 case 11:
60 return detail::HardcodedQuadRule<base::RefEl::kTria(), 11>();
61 case 12:
62 return detail::HardcodedQuadRule<base::RefEl::kTria(), 12>();
63 case 13:
64 return detail::HardcodedQuadRule<base::RefEl::kTria(), 13>();
65 case 14:
66 return detail::HardcodedQuadRule<base::RefEl::kTria(), 14>();
67 case 15:
68 return detail::HardcodedQuadRule<base::RefEl::kTria(), 15>();
69 case 16:
70 return detail::HardcodedQuadRule<base::RefEl::kTria(), 16>();
71 case 17:
72 return detail::HardcodedQuadRule<base::RefEl::kTria(), 17>();
73 case 18:
74 return detail::HardcodedQuadRule<base::RefEl::kTria(), 18>();
75 case 19:
76 return detail::HardcodedQuadRule<base::RefEl::kTria(), 19>();
77 case 20:
78 return detail::HardcodedQuadRule<base::RefEl::kTria(), 20>();
79 case 21:
80 return detail::HardcodedQuadRule<base::RefEl::kTria(), 21>();
81 case 22:
82 return detail::HardcodedQuadRule<base::RefEl::kTria(), 22>();
83 case 23:
84 return detail::HardcodedQuadRule<base::RefEl::kTria(), 23>();
85 case 24:
86 return detail::HardcodedQuadRule<base::RefEl::kTria(), 24>();
87 case 25:
88 return detail::HardcodedQuadRule<base::RefEl::kTria(), 25>();
89 case 26:
90 return detail::HardcodedQuadRule<base::RefEl::kTria(), 26>();
91 case 27:
92 return detail::HardcodedQuadRule<base::RefEl::kTria(), 27>();
93 case 28:
94 return detail::HardcodedQuadRule<base::RefEl::kTria(), 28>();
95 case 29:
96 return detail::HardcodedQuadRule<base::RefEl::kTria(), 29>();
97 case 30:
98 return detail::HardcodedQuadRule<base::RefEl::kTria(), 30>();
99 case 31:
100 return detail::HardcodedQuadRule<base::RefEl::kTria(), 31>();
101 case 32:
102 return detail::HardcodedQuadRule<base::RefEl::kTria(), 32>();
103 case 33:
104 return detail::HardcodedQuadRule<base::RefEl::kTria(), 33>();
105 case 34:
106 return detail::HardcodedQuadRule<base::RefEl::kTria(), 34>();
107 case 35:
108 return detail::HardcodedQuadRule<base::RefEl::kTria(), 35>();
109 case 36:
110 return detail::HardcodedQuadRule<base::RefEl::kTria(), 36>();
111 case 37:
112 return detail::HardcodedQuadRule<base::RefEl::kTria(), 37>();
113 case 38:
114 return detail::HardcodedQuadRule<base::RefEl::kTria(), 38>();
115 case 39:
116 return detail::HardcodedQuadRule<base::RefEl::kTria(), 39>();
117 case 40:
118 return detail::HardcodedQuadRule<base::RefEl::kTria(), 40>();
119 case 41:
120 return detail::HardcodedQuadRule<base::RefEl::kTria(), 41>();
121 case 42:
122 return detail::HardcodedQuadRule<base::RefEl::kTria(), 42>();
123 case 43:
124 return detail::HardcodedQuadRule<base::RefEl::kTria(), 43>();
125 case 44:
126 return detail::HardcodedQuadRule<base::RefEl::kTria(), 44>();
127 case 45:
128 return detail::HardcodedQuadRule<base::RefEl::kTria(), 45>();
129 case 46:
130 return detail::HardcodedQuadRule<base::RefEl::kTria(), 46>();
131 case 47:
132 return detail::HardcodedQuadRule<base::RefEl::kTria(), 47>();
133 case 48:
134 return detail::HardcodedQuadRule<base::RefEl::kTria(), 48>();
135 case 49:
136 return detail::HardcodedQuadRule<base::RefEl::kTria(), 49>();
137 case 50:
138 return detail::HardcodedQuadRule<base::RefEl::kTria(), 50>();
139 default:
140 // Create a quadrule using tensor product quadrature rule + duffy
141 // transform
142 quadDegree_t n = degree / 2 + 1;
143 auto [leg_p, leg_w] = GaussLegendre(n);
144 auto [jac_p, jac_w] = GaussJacobi(n, 1, 0);
145 jac_p.array() = (jac_p.array() + 1) / 2.; // rescale to [0,1]
146 jac_w.array() *= 0.25;
147 Eigen::MatrixXd points2d(2, n * n);
148 points2d.row(0) = Eigen::kroneckerProduct(
149 leg_p.transpose(), (1 - jac_p.transpose().array()).matrix());
150 points2d.row(1) = jac_p.transpose().replicate(1, n);
151 return QuadRule(base::RefEl::kTria(), std::move(points2d),
152 Eigen::kroneckerProduct(leg_w, jac_w), 2 * n - 1);
153 }
154 }
155 LF_VERIFY_MSG(
156 false, "No Quadrature rules implemented for this reference element yet.");
157}
158
159QuadRule make_TriaQR_MidpointRule() {
160 Eigen::MatrixXd points(2, 1);
161 Eigen::VectorXd weights(1);
162
163 points(0, 0) = 1.0 / 3.0;
164 points(1, 0) = 1.0 / 3.0;
165
166 weights(0) = 0.5;
167
168 return QuadRule(base::RefEl::kTria(), std::move(points), std::move(weights),
169 1);
170}
171
172QuadRule make_TriaQR_EdgeMidpointRule() {
173 Eigen::MatrixXd points(2, 3);
174 Eigen::VectorXd weights(3);
175
176 points(0, 0) = 0.5;
177 points(1, 0) = 0.0;
178 points(0, 1) = 0.5;
179 points(1, 1) = 0.5;
180 points(0, 2) = 0.0;
181 points(1, 2) = 0.5;
182
183 weights(0) = 0.16666666666666665741;
184 weights(1) = 0.16666666666666665741;
185 weights(2) = 0.16666666666666665741;
186
187 return QuadRule(base::RefEl::kTria(), std::move(points), std::move(weights),
188 2);
189}
190
191QuadRule make_QuadQR_EdgeMidpointRule() {
192 Eigen::MatrixXd points(2, 4);
193 Eigen::VectorXd weights(4);
194
195 points(0, 0) = 0.5;
196 points(1, 0) = 0.0;
197 points(0, 1) = 1.0;
198 points(1, 1) = 0.5;
199 points(0, 2) = 0.5;
200 points(1, 2) = 1.0;
201 points(0, 3) = 0.0;
202 points(1, 3) = 0.5;
203
204 weights(0) = 0.25;
205 weights(1) = 0.25;
206 weights(2) = 0.25;
207 weights(3) = 0.25;
208
209 return QuadRule(base::RefEl::kQuad(), std::move(points), std::move(weights),
210 1);
211}
212
213QuadRule make_TriaQR_P6O4() {
214 Eigen::MatrixXd points(2, 6);
215 Eigen::VectorXd weights(6);
216
217 points(0, 0) = 0.5;
218 points(1, 0) = 0.0;
219 points(0, 1) = 0.5;
220 points(1, 1) = 0.5;
221 points(0, 2) = 0.0;
222 points(1, 2) = 0.5;
223 points(0, 3) = 1.0 / 6.0;
224 points(1, 3) = 1.0 / 6.0;
225 points(0, 4) = 1.0 / 6.0;
226 points(1, 4) = 2.0 / 3.0;
227 points(0, 5) = 2.0 / 3.0;
228 points(1, 5) = 1.0 / 6.0;
229
230 weights << 1.0, 1.0, 1.0, 9.0, 9.0, 9.0;
231 weights /= 60.0;
232
233 return QuadRule(base::RefEl::kTria(), std::move(points), std::move(weights),
234 3);
235}
236
237QuadRule make_TriaQR_P7O6() {
238 return detail::HardcodedQuadRule<base::RefEl::kTria(), 5>();
239}
240
241} // namespace lf::quad
lf::base::RefEl
Represents a reference element with all its properties.
Definition: ref_el.h:106
lf::base::RefEl::kSegment
static constexpr RefEl kSegment()
Returns the (1-dimensional) reference segment.
Definition: ref_el.h:150
lf::base::RefEl::kTria
static constexpr RefEl kTria()
Returns the reference triangle.
Definition: ref_el.h:158
lf::base::RefEl::kQuad
static constexpr RefEl kQuad()
Returns the reference quadrilateral.
Definition: ref_el.h:166
lf::quad::QuadRule
Represents a Quadrature Rule over one of the Reference Elements.
Definition: quad_rule.h:58
lf::quad::make_TriaQR_MidpointRule
QuadRule make_TriaQR_MidpointRule()
midpoint quadrature rule for triangles
Definition: make_quad_rule.cc:159
lf::quad::make_TriaQR_P7O6
QuadRule make_TriaQR_P7O6()
Seven point triangular quadrature rule of order 6.
Definition: make_quad_rule.cc:237
lf::quad::make_QuadQR_EdgeMidpointRule
QuadRule make_QuadQR_EdgeMidpointRule()
edge midpoint quadrature rule for unit square (= reference quad)
Definition: make_quad_rule.cc:191
lf::quad::make_TriaQR_P6O4
QuadRule make_TriaQR_P6O4()
Six point triangular quadrature rule of order 4.
Definition: make_quad_rule.cc:213
lf::quad::make_TriaQR_EdgeMidpointRule
QuadRule make_TriaQR_EdgeMidpointRule()
edge midpoint quadrature rule for reference triangles
Definition: make_quad_rule.cc:172
lf::quad
Rules for numerical quadrature on reference entity shapes.
Definition: gauss_quadrature.cc:15
lf::quad::GaussJacobi
std::tuple< Eigen::VectorXd, Eigen::VectorXd > GaussJacobi(quadDegree_t num_points, double alpha, double beta)
Computes the quadrature points and weights for the interval [-1,1] of a Gauss-Jacobi quadrature rule ...
Definition: gauss_quadrature.cc:71
lf::quad::GaussLegendre
std::tuple< Eigen::VectorXd, Eigen::VectorXd > GaussLegendre(unsigned num_points)
Definition: gauss_quadrature.cc:16
lf::quad::make_QuadRule
QuadRule make_QuadRule(base::RefEl ref_el, unsigned degree)
Returns a QuadRule object for the given Reference Element and Degree.
Definition: make_quad_rule.cc:21
lf::quad::quadDegree_t
unsigned int quadDegree_t
Definition: quad_rule.h:19
Generated by doxygen 1.9.3