Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #ifndef PALACE_FEM_INTEGRATOR_HPP
5 : #define PALACE_FEM_INTEGRATOR_HPP
6 :
7 : #include <mfem.hpp>
8 : #include "fem/libceed/ceed.hpp"
9 :
10 : namespace palace
11 : {
12 :
13 : class MaterialPropertyCoefficient;
14 :
15 : //
16 : // Classes which implement or extend bilinear and linear form integrators. In doc strings u
17 : // refers to the trial function, and v the test function.
18 : //
19 :
20 : namespace fem
21 : {
22 :
23 : // Helper functions for creating an integration rule to exactly integrate polynomials of
24 : // order 2 * p_trial + order(|J|) + q_extra.
25 : struct DefaultIntegrationOrder
26 : {
27 : inline static int p_trial = 1;
28 : inline static bool q_order_jac = true;
29 : inline static int q_order_extra_pk = 0;
30 : inline static int q_order_extra_qk = 0;
31 : static int Get(const mfem::IsoparametricTransformation &T);
32 : static int Get(const mfem::ElementTransformation &T);
33 : static int Get(const mfem::Mesh &mesh, mfem::Geometry::Type geom);
34 : };
35 :
36 : } // namespace fem
37 :
38 : // Base class for libCEED-based bilinear form integrators.
39 : class BilinearFormIntegrator
40 : {
41 : protected:
42 : const MaterialPropertyCoefficient *Q;
43 : bool assemble_q_data;
44 : bool transpose;
45 :
46 : public:
47 : BilinearFormIntegrator(const MaterialPropertyCoefficient *Q = nullptr,
48 : const bool transpose = false)
49 2238 : : Q(Q), assemble_q_data(false), transpose(transpose)
50 : {
51 : }
52 : BilinearFormIntegrator(const MaterialPropertyCoefficient &Q, const bool transpose = false)
53 4002 : : Q(&Q), assemble_q_data(false), transpose(transpose)
54 : {
55 : }
56 : virtual ~BilinearFormIntegrator() = default;
57 :
58 : virtual void Assemble(Ceed ceed, CeedElemRestriction trial_restr,
59 : CeedElemRestriction test_restr, CeedBasis trial_basis,
60 : CeedBasis test_basis, CeedVector geom_data,
61 : CeedElemRestriction geom_data_restr, CeedOperator *op) const = 0;
62 :
63 7665 : virtual void SetMapTypes(int trial_type, int test_type) {}
64 :
65 0 : void AssembleQuadratureData() { assemble_q_data = true; }
66 : };
67 :
68 : // Integrator for a(u, v) = (Q u, v) for H1 elements (also for vector (H1)ᵈ spaces).
69 444 : class MassIntegrator : public BilinearFormIntegrator
70 : {
71 : public:
72 654 : using BilinearFormIntegrator::BilinearFormIntegrator;
73 :
74 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
75 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
76 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
77 : };
78 :
79 : // Integrator for a(u, v) = (Q u, v) for vector finite elements.
80 552 : class VectorFEMassIntegrator : public BilinearFormIntegrator
81 : {
82 : protected:
83 : int trial_map_type = mfem::FiniteElement::UNKNOWN_MAP_TYPE;
84 : int test_map_type = mfem::FiniteElement::UNKNOWN_MAP_TYPE;
85 :
86 : public:
87 1119 : using BilinearFormIntegrator::BilinearFormIntegrator;
88 :
89 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
90 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
91 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
92 :
93 4127 : void SetMapTypes(int trial_type, int test_type) override
94 : {
95 4127 : trial_map_type = trial_type;
96 4127 : test_map_type = test_type;
97 4127 : }
98 : };
99 :
100 : // Integrator for a(u, v) = (Q grad u, grad v) for H1 elements.
101 216 : class DiffusionIntegrator : public BilinearFormIntegrator
102 : {
103 : public:
104 441 : using BilinearFormIntegrator::BilinearFormIntegrator;
105 :
106 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
107 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
108 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
109 : };
110 :
111 : // Integrator for a(u, v) = (Q curl u, curl v) for Nedelec elements.
112 162 : class CurlCurlIntegrator : public BilinearFormIntegrator
113 : {
114 : public:
115 234 : using BilinearFormIntegrator::BilinearFormIntegrator;
116 :
117 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
118 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
119 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
120 : };
121 :
122 : // Integrator for a(u, v) = (Q div u, div v) for Raviart-Thomas elements.
123 108 : class DivDivIntegrator : public BilinearFormIntegrator
124 : {
125 : public:
126 114 : using BilinearFormIntegrator::BilinearFormIntegrator;
127 :
128 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
129 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
130 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
131 : };
132 :
133 : // Integrator for a(u, v) = (Qd grad u, grad v) + (Qm u, v) for H1 elements.
134 : class DiffusionMassIntegrator : public BilinearFormIntegrator
135 : {
136 : protected:
137 : const MaterialPropertyCoefficient *Q_mass;
138 : bool transpose_mass;
139 :
140 : public:
141 : using BilinearFormIntegrator::BilinearFormIntegrator;
142 : DiffusionMassIntegrator(const MaterialPropertyCoefficient &Q,
143 : const MaterialPropertyCoefficient &Q_mass,
144 : const bool transpose = false, const bool transpose_mass = false)
145 12 : : BilinearFormIntegrator(Q, transpose), Q_mass(&Q_mass), transpose_mass(transpose_mass)
146 : {
147 : }
148 :
149 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
150 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
151 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
152 : };
153 :
154 : // Integrator for a(u, v) = (Qc curl u, curl v) + (Qm u, v) for Nedelec elements.
155 : class CurlCurlMassIntegrator : public BilinearFormIntegrator
156 : {
157 : protected:
158 : const MaterialPropertyCoefficient *Q_mass;
159 : bool transpose_mass;
160 :
161 : public:
162 : using BilinearFormIntegrator::BilinearFormIntegrator;
163 : CurlCurlMassIntegrator(const MaterialPropertyCoefficient &Q,
164 : const MaterialPropertyCoefficient &Q_mass,
165 : const bool transpose = false, const bool transpose_mass = false)
166 12 : : BilinearFormIntegrator(Q, transpose), Q_mass(&Q_mass), transpose_mass(transpose_mass)
167 : {
168 : }
169 :
170 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
171 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
172 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
173 : };
174 :
175 : // Integrator for a(u, v) = (Qd div u, div v) + (Qm u, v) for Raviart-Thomas elements.
176 : class DivDivMassIntegrator : public BilinearFormIntegrator
177 : {
178 : protected:
179 : const MaterialPropertyCoefficient *Q_mass;
180 : bool transpose_mass;
181 :
182 : public:
183 : using BilinearFormIntegrator::BilinearFormIntegrator;
184 : DivDivMassIntegrator(const MaterialPropertyCoefficient &Q,
185 : const MaterialPropertyCoefficient &Q_mass,
186 : const bool transpose = false, const bool transpose_mass = false)
187 12 : : BilinearFormIntegrator(Q, transpose), Q_mass(&Q_mass), transpose_mass(transpose_mass)
188 : {
189 : }
190 :
191 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
192 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
193 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
194 : };
195 :
196 : // Integrator for a(u, v) = (Q grad u, v) for u in H1 and v in H(curl) or H(div).
197 270 : class MixedVectorGradientIntegrator : public BilinearFormIntegrator
198 : {
199 : protected:
200 : int trial_map_type = mfem::FiniteElement::UNKNOWN_MAP_TYPE;
201 : int test_map_type = mfem::FiniteElement::UNKNOWN_MAP_TYPE;
202 :
203 : public:
204 540 : using BilinearFormIntegrator::BilinearFormIntegrator;
205 :
206 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
207 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
208 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
209 :
210 2067 : void SetMapTypes(int trial_type, int test_type) override
211 : {
212 2067 : trial_map_type = trial_type;
213 2067 : test_map_type = test_type;
214 2067 : }
215 : };
216 :
217 : // Integrator for a(u, v) = -(Q u, grad v) for u in H(curl) and v in H1.
218 162 : class MixedVectorWeakDivergenceIntegrator : public BilinearFormIntegrator
219 : {
220 : public:
221 324 : using BilinearFormIntegrator::BilinearFormIntegrator;
222 :
223 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
224 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
225 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
226 : };
227 :
228 : // Integrator for a(u, v) = (Q curl u, v) for u in H(curl) and v in H(div) or H(curl).
229 108 : class MixedVectorCurlIntegrator : public BilinearFormIntegrator
230 : {
231 : protected:
232 : int trial_map_type = mfem::FiniteElement::UNKNOWN_MAP_TYPE;
233 : int test_map_type = mfem::FiniteElement::UNKNOWN_MAP_TYPE;
234 :
235 : public:
236 216 : using BilinearFormIntegrator::BilinearFormIntegrator;
237 :
238 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
239 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
240 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
241 :
242 774 : void SetMapTypes(int trial_type, int test_type) override
243 : {
244 774 : trial_map_type = trial_type;
245 774 : test_map_type = test_type;
246 774 : }
247 : };
248 :
249 : // Integrator for a(u, v) = -(Q u, curl v) for u in H(div) or H(curl) and v in H(curl).
250 108 : class MixedVectorWeakCurlIntegrator : public BilinearFormIntegrator
251 : {
252 : protected:
253 : int trial_map_type = mfem::FiniteElement::UNKNOWN_MAP_TYPE;
254 : int test_map_type = mfem::FiniteElement::UNKNOWN_MAP_TYPE;
255 :
256 : public:
257 216 : using BilinearFormIntegrator::BilinearFormIntegrator;
258 :
259 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
260 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
261 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
262 :
263 774 : void SetMapTypes(int trial_type, int test_type) override
264 : {
265 774 : trial_map_type = trial_type;
266 774 : test_map_type = test_type;
267 774 : }
268 : };
269 :
270 : // Integrator for a(u, v) = (Q grad u, v) for u in H1 and v in (H1)ᵈ.
271 108 : class GradientIntegrator : public BilinearFormIntegrator
272 : {
273 : public:
274 108 : using BilinearFormIntegrator::BilinearFormIntegrator;
275 :
276 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
277 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector geom_data,
278 : CeedElemRestriction geom_data_restr, CeedOperator *op) const override;
279 : };
280 :
281 : // Base class for all discrete interpolators.
282 : class DiscreteInterpolator
283 : {
284 : public:
285 : void Assemble(Ceed ceed, CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
286 : CeedBasis interp_basis, CeedOperator *op, CeedOperator *op_t);
287 : };
288 :
289 : // Interpolator for the identity map, where the domain space is a subspace of the range
290 : // space (discrete embedding matrix).
291 : using IdentityInterpolator = DiscreteInterpolator;
292 :
293 : // Interpolator for the discrete gradient map from an H1 space to an H(curl) space.
294 : using GradientInterpolator = DiscreteInterpolator;
295 :
296 : // Interpolator for the discrete curl map from an H(curl) space to an H(div) space.
297 : using CurlInterpolator = DiscreteInterpolator;
298 :
299 : // Interpolator for the discrete divergence map from an H(div) space to an L2 space.
300 : using DivergenceInterpolator = DiscreteInterpolator;
301 :
302 : // Similar to MFEM's VectorFEBoundaryTangentLFIntegrator for ND spaces, but instead of
303 : // computing (n x f, v), this just computes (f, v). Also eliminates the a and b quadrature
304 : // parameters and uses fem::DefaultIntegrationOrder instead.
305 : class VectorFEBoundaryLFIntegrator : public mfem::LinearFormIntegrator
306 : {
307 : private:
308 : mfem::VectorCoefficient &Q;
309 : mfem::DenseMatrix vshape;
310 : mfem::Vector f_loc, f_hat;
311 :
312 : public:
313 21 : VectorFEBoundaryLFIntegrator(mfem::VectorCoefficient &QG) : Q(QG) {}
314 :
315 : void AssembleRHSElementVect(const mfem::FiniteElement &fe, mfem::ElementTransformation &T,
316 : mfem::Vector &elvect) override;
317 : };
318 :
319 : // Similar to MFEM's BoundaryLFIntegrator for H1 spaces, but eliminates the a and b
320 : // quadrature parameters and uses fem::DefaultIntegrationOrder instead.
321 : class BoundaryLFIntegrator : public mfem::LinearFormIntegrator
322 : {
323 : private:
324 : mfem::Coefficient &Q;
325 : mfem::Vector shape;
326 :
327 : public:
328 0 : BoundaryLFIntegrator(mfem::Coefficient &QG) : Q(QG) {}
329 :
330 : void AssembleRHSElementVect(const mfem::FiniteElement &fe, mfem::ElementTransformation &T,
331 : mfem::Vector &elvect) override;
332 : };
333 :
334 : using VectorFEDomainLFIntegrator = VectorFEBoundaryLFIntegrator;
335 : using DomainLFIntegrator = BoundaryLFIntegrator;
336 :
337 : } // namespace palace
338 :
339 : #endif // PALACE_FEM_INTEGRATOR_HPP
|
