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_COEFFICIENT_HPP
5 : #define PALACE_FEM_COEFFICIENT_HPP
6 :
7 : #include <complex>
8 : #include <memory>
9 : #include <utility>
10 : #include <vector>
11 : #include <mfem.hpp>
12 : #include "fem/gridfunction.hpp"
13 : #include "linalg/vector.hpp"
14 : #include "models/materialoperator.hpp"
15 : #include "utils/geodata.hpp"
16 : #include "utils/labels.hpp"
17 :
18 : // XX TODO: Add bulk element Eval() overrides to speed up postprocessing (also needed in
19 : // mfem::DataCollection classes.
20 :
21 : namespace palace
22 : {
23 :
24 : //
25 : // Derived coefficients which compute single values on internal boundaries where a possibly
26 : // discontinuous function is given as an input grid function. These are all cheap to
27 : // construct by design. All methods assume the provided grid function is ready for parallel
28 : // comm on shared faces after a call to ExchangeFaceNbrData.
29 : //
30 :
31 : // Base class for coefficients which need to evaluate a GridFunction in a domain element
32 : // attached to a boundary element, or both domain elements on either side for internal
33 : // boundaries.
34 : class BdrGridFunctionCoefficient
35 : {
36 : protected:
37 : // XX TODO: For thread-safety (multiple threads evaluating a coefficient simultaneously),
38 : // the FET, FET.Elem1, and FET.Elem2 objects cannot be shared.
39 : const mfem::ParMesh &mesh;
40 : mfem::FaceElementTransformations FET;
41 : mfem::IsoparametricTransformation T1, T2;
42 :
43 : bool GetBdrElementNeighborTransformations(int i, const mfem::IntegrationPoint &ip)
44 : {
45 : // Get the element transformations neighboring the element, and optionally set the
46 : // integration point too.
47 24300 : return GetBdrElementNeighborTransformations(i, mesh, FET, T1, T2, &ip);
48 : }
49 :
50 : public:
51 120 : BdrGridFunctionCoefficient(const mfem::ParMesh &mesh) : mesh(mesh) {}
52 :
53 : // For a boundary element, return the element transformation objects for the neighboring
54 : // domain elements. FET.Elem2 may be nullptr if the boundary is a true one-sided boundary,
55 : // but if it is shared with another subdomain then it will be populated. Expects
56 : // ParMesh::ExchangeFaceNbrData has been called already.
57 : static bool GetBdrElementNeighborTransformations(
58 : int i, const mfem::ParMesh &mesh, mfem::FaceElementTransformations &FET,
59 : mfem::IsoparametricTransformation &T1, mfem::IsoparametricTransformation &T2,
60 : const mfem::IntegrationPoint *ip = nullptr);
61 :
62 : // Return normal vector to the boundary element at an integration point. For a face
63 : // element, the normal points out of the element (from element 1 into element 2, if it
64 : // exists). This convention can be flipped with the optional parameter. It is assumed
65 : // that the element transformation has already been configured at the integration point
66 : // of interest.
67 36900 : static void GetNormal(mfem::ElementTransformation &T, mfem::Vector &normal,
68 : bool invert = false)
69 : {
70 : MFEM_ASSERT(normal.Size() == T.GetSpaceDim(),
71 : "Size mismatch for normal vector (space dimension = " << T.GetSpaceDim()
72 : << ")!");
73 36900 : mfem::CalcOrtho(T.Jacobian(), normal);
74 36900 : normal /= invert ? -normal.Norml2() : normal.Norml2();
75 36900 : }
76 : };
77 :
78 : // Computes surface current Jₛ = n x H = n x μ⁻¹ B on boundaries from B as a vector grid
79 : // function where n is an inward normal (computes -n x H for outward normal n). For a
80 : // two-sided internal boundary, the contributions from both sides add.
81 0 : class BdrSurfaceCurrentVectorCoefficient : public mfem::VectorCoefficient,
82 : public BdrGridFunctionCoefficient
83 : {
84 : private:
85 : const mfem::ParGridFunction &B;
86 : const MaterialOperator &mat_op;
87 :
88 : public:
89 27 : BdrSurfaceCurrentVectorCoefficient(const mfem::ParGridFunction &B,
90 : const MaterialOperator &mat_op)
91 27 : : mfem::VectorCoefficient(B.VectorDim()),
92 27 : BdrGridFunctionCoefficient(*B.ParFESpace()->GetParMesh()), B(B), mat_op(mat_op)
93 : {
94 27 : }
95 :
96 : using mfem::VectorCoefficient::Eval;
97 4212 : void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
98 : const mfem::IntegrationPoint &ip) override
99 : {
100 : // Get neighboring elements.
101 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
102 : "Unexpected element type in BdrSurfaceCurrentVectorCoefficient!");
103 4212 : bool ori = GetBdrElementNeighborTransformations(T.ElementNo, ip);
104 :
105 : // For interior faces, compute Jₛ = n x H = n x μ⁻¹ (B1 - B2), where B1 (B2) is B in
106 : // element 1 (element 2) and n points into element 1.
107 : double W_data[3], VU_data[3];
108 8424 : mfem::Vector W(W_data, vdim), VU(VU_data, vdim);
109 4212 : B.GetVectorValue(*FET.Elem1, FET.Elem1->GetIntPoint(), W);
110 4212 : mat_op.GetInvPermeability(FET.Elem1->Attribute).Mult(W, VU);
111 4212 : if (FET.Elem2)
112 : {
113 : // Double-sided, not a true boundary. Add result with opposite normal.
114 : double VL_data[3];
115 0 : mfem::Vector VL(VL_data, vdim);
116 0 : B.GetVectorValue(*FET.Elem2, FET.Elem2->GetIntPoint(), W);
117 0 : mat_op.GetInvPermeability(FET.Elem2->Attribute).Mult(W, VL);
118 0 : VU -= VL;
119 : }
120 :
121 : // Orient with normal pointing into element 1.
122 : double normal_data[3];
123 4212 : mfem::Vector normal(normal_data, vdim);
124 4212 : GetNormal(T, normal, ori);
125 4212 : V.SetSize(vdim);
126 4212 : linalg::Cross3(normal, VU, V);
127 4212 : }
128 : };
129 :
130 : // Computes the flux Φₛ = F ⋅ n with F = B or ε D on interior boundary elements using B or
131 : // E given as a vector grid function. For a two-sided internal boundary, the contributions
132 : // from both sides can either add or be averaged.
133 : template <SurfaceFlux Type>
134 0 : class BdrSurfaceFluxCoefficient : public mfem::Coefficient,
135 : public BdrGridFunctionCoefficient
136 : {
137 : private:
138 : const mfem::ParGridFunction *E, *B;
139 : const MaterialOperator &mat_op;
140 : bool two_sided;
141 : const mfem::Vector &x0;
142 :
143 : void GetLocalFlux(mfem::ElementTransformation &T, mfem::Vector &V) const;
144 :
145 : public:
146 18 : BdrSurfaceFluxCoefficient(const mfem::ParGridFunction *E, const mfem::ParGridFunction *B,
147 : const MaterialOperator &mat_op, bool two_sided,
148 : const mfem::Vector &x0)
149 : : mfem::Coefficient(), BdrGridFunctionCoefficient(E ? *E->ParFESpace()->GetParMesh()
150 : : *B->ParFESpace()->GetParMesh()),
151 36 : E(E), B(B), mat_op(mat_op), two_sided(two_sided), x0(x0)
152 : {
153 18 : MFEM_VERIFY((E || (Type != SurfaceFlux::ELECTRIC && Type != SurfaceFlux::POWER)) &&
154 : (B || (Type != SurfaceFlux::MAGNETIC && Type != SurfaceFlux::POWER)),
155 : "Missing E or B field grid function for surface flux coefficient!");
156 18 : }
157 :
158 3888 : double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override
159 : {
160 : // Get neighboring elements.
161 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
162 : "Unexpected element type in BdrSurfaceFluxCoefficient!");
163 3888 : bool ori = GetBdrElementNeighborTransformations(T.ElementNo, ip);
164 :
165 : // For interior faces, compute either F ⋅ n as the average or by adding the
166 : // contributions from opposite sides with opposite normals.
167 3888 : const int vdim = T.GetSpaceDim();
168 : double VU_data[3];
169 : mfem::Vector VU(VU_data, vdim);
170 3888 : GetLocalFlux(*FET.Elem1, VU);
171 3888 : if (FET.Elem2)
172 : {
173 : // Double-sided, not a true boundary.
174 : double VL_data[3];
175 : mfem::Vector VL(VL_data, vdim);
176 0 : GetLocalFlux(*FET.Elem2, VL);
177 0 : if (two_sided)
178 : {
179 : // Add result with opposite normal. This only happens when crack_bdr_elements =
180 : // false (two_sided = true doesn't make sense for an internal boundary without an
181 : // associated BC).
182 0 : VU -= VL;
183 : }
184 : else
185 : {
186 : // Take the average of the values on both sides.
187 0 : add(0.5, VU, VL, VU);
188 : }
189 : }
190 :
191 : // Dot with normal direction and assign appropriate sign. The normal is oriented to
192 : // point into element 1.
193 : double normal_data[3];
194 : mfem::Vector normal(normal_data, vdim);
195 3888 : GetNormal(T, normal, ori);
196 3888 : double flux = VU * normal;
197 3888 : if (two_sided)
198 : {
199 : return flux;
200 : }
201 : else
202 : {
203 : // Orient outward from the surface with the given center.
204 : double x_data[3];
205 : mfem::Vector x(x_data, vdim);
206 0 : T.Transform(ip, x);
207 0 : x -= x0;
208 0 : return (x * normal < 0.0) ? -flux : flux;
209 : }
210 : }
211 : };
212 :
213 : template <>
214 3888 : inline void BdrSurfaceFluxCoefficient<SurfaceFlux::ELECTRIC>::GetLocalFlux(
215 : mfem::ElementTransformation &T, mfem::Vector &V) const
216 : {
217 : // Flux D.
218 : double W_data[3];
219 3888 : mfem::Vector W(W_data, T.GetSpaceDim());
220 3888 : E->GetVectorValue(T, T.GetIntPoint(), W);
221 3888 : mat_op.GetPermittivityReal(T.Attribute).Mult(W, V);
222 3888 : }
223 :
224 : template <>
225 : inline void BdrSurfaceFluxCoefficient<SurfaceFlux::MAGNETIC>::GetLocalFlux(
226 : mfem::ElementTransformation &T, mfem::Vector &V) const
227 : {
228 : // Flux B.
229 0 : B->GetVectorValue(T, T.GetIntPoint(), V);
230 0 : }
231 :
232 : template <>
233 : inline void
234 0 : BdrSurfaceFluxCoefficient<SurfaceFlux::POWER>::GetLocalFlux(mfem::ElementTransformation &T,
235 : mfem::Vector &V) const
236 : {
237 : // Flux E x H = E x μ⁻¹ B.
238 : double W1_data[3], W2_data[3];
239 0 : mfem::Vector W1(W1_data, T.GetSpaceDim()), W2(W2_data, T.GetSpaceDim());
240 0 : B->GetVectorValue(T, T.GetIntPoint(), W1);
241 0 : mat_op.GetInvPermeability(T.Attribute).Mult(W1, W2);
242 0 : E->GetVectorValue(T, T.GetIntPoint(), W1);
243 0 : V.SetSize(W1.Size());
244 0 : linalg::Cross3(W1, W2, V);
245 0 : }
246 :
247 : // Computes a single-valued α Eᵀ E on boundaries from E given as a vector grid function.
248 : // Uses the neighbor element on a user specified side to compute a single-sided value for
249 : // potentially discontinuous solutions for an interior boundary element. The four cases
250 : // correspond to a generic interface vs. specializations for metal-air, metal-substrate,
251 : // and substrate-air interfaces following:
252 : // J. Wenner et al., Surface loss simulations of superconducting coplanar waveguide
253 : // resonators, Appl. Phys. Lett. (2011).
254 : template <InterfaceDielectric Type>
255 0 : class InterfaceDielectricCoefficient : public mfem::Coefficient,
256 : public BdrGridFunctionCoefficient
257 : {
258 : private:
259 : const GridFunction &E;
260 : const MaterialOperator &mat_op;
261 : const double t_i, epsilon_i;
262 :
263 0 : void Initialize(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip,
264 : mfem::Vector *normal)
265 : {
266 : // Get neighboring elements and the normal vector, oriented to point into element 1.
267 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
268 : "Unexpected element type in InterfaceDielectricCoefficient!");
269 0 : bool ori = GetBdrElementNeighborTransformations(T.ElementNo, ip);
270 0 : if (normal)
271 : {
272 0 : GetNormal(T, *normal, ori);
273 : }
274 0 : }
275 :
276 0 : int GetLocalVectorValue(const mfem::ParGridFunction &U, mfem::Vector &V,
277 : bool vacuum_side) const
278 : {
279 : constexpr double threshold = 1.0 - 1.0e-6;
280 : const bool use_elem1 =
281 0 : ((vacuum_side && mat_op.GetLightSpeedMax(FET.Elem1->Attribute) >= threshold) ||
282 0 : (!vacuum_side && mat_op.GetLightSpeedMax(FET.Elem1->Attribute) < threshold));
283 : const bool use_elem2 =
284 0 : (FET.Elem2 &&
285 0 : ((vacuum_side && mat_op.GetLightSpeedMax(FET.Elem2->Attribute) >= threshold) ||
286 0 : (!vacuum_side && mat_op.GetLightSpeedMax(FET.Elem2->Attribute) < threshold)));
287 0 : if (use_elem1)
288 : {
289 0 : U.GetVectorValue(*FET.Elem1, FET.Elem1->GetIntPoint(), V);
290 0 : if (use_elem2)
291 : {
292 : // Double-sided, not a true boundary. Just average the solution from both sides.
293 : double W_data[3];
294 : mfem::Vector W(W_data, V.Size());
295 0 : U.GetVectorValue(*FET.Elem2, FET.Elem2->GetIntPoint(), W);
296 0 : add(0.5, V, W, V);
297 : }
298 0 : return FET.Elem1->Attribute;
299 : }
300 0 : else if (use_elem2)
301 : {
302 0 : U.GetVectorValue(*FET.Elem2, FET.Elem2->GetIntPoint(), V);
303 0 : return FET.Elem2->Attribute;
304 : }
305 : else
306 : {
307 : return 0;
308 : }
309 : }
310 :
311 : public:
312 0 : InterfaceDielectricCoefficient(const GridFunction &E, const MaterialOperator &mat_op,
313 : double t_i, double epsilon_i)
314 0 : : mfem::Coefficient(), BdrGridFunctionCoefficient(*E.ParFESpace()->GetParMesh()), E(E),
315 0 : mat_op(mat_op), t_i(t_i), epsilon_i(epsilon_i)
316 : {
317 : }
318 :
319 : double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override;
320 : };
321 :
322 : template <>
323 0 : inline double InterfaceDielectricCoefficient<InterfaceDielectric::DEFAULT>::Eval(
324 : mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
325 : {
326 : // Get single-sided solution. Don't use lightspeed detection for differentiating side.
327 0 : auto GetLocalVectorValueDefault = [this](const mfem::ParGridFunction &U, mfem::Vector &V)
328 : {
329 0 : U.GetVectorValue(*FET.Elem1, FET.Elem1->GetIntPoint(), V);
330 0 : if (FET.Elem2)
331 : {
332 : // Double-sided, not a true boundary. Just average the field solution from both sides.
333 : double W_data[3];
334 : mfem::Vector W(W_data, V.Size());
335 0 : U.GetVectorValue(*FET.Elem2, FET.Elem2->GetIntPoint(), W);
336 0 : add(0.5, V, W, V);
337 : }
338 0 : };
339 : double V_data[3];
340 0 : mfem::Vector V(V_data, T.GetSpaceDim());
341 : Initialize(T, ip, nullptr);
342 0 : GetLocalVectorValueDefault(E.Real(), V);
343 0 : double V2 = V * V;
344 0 : if (E.HasImag())
345 : {
346 0 : GetLocalVectorValueDefault(E.Imag(), V);
347 0 : V2 += V * V;
348 : }
349 :
350 : // No specific interface, use full field evaluation: 0.5 * t * ε * |E|² .
351 0 : return 0.5 * t_i * epsilon_i * V2;
352 : }
353 :
354 : template <>
355 0 : inline double InterfaceDielectricCoefficient<InterfaceDielectric::MA>::Eval(
356 : mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
357 : {
358 : // Get single-sided solution on air (vacuum) side and neighboring element attribute.
359 : double V_data[3], normal_data[3];
360 0 : mfem::Vector V(V_data, T.GetSpaceDim()), normal(normal_data, T.GetSpaceDim());
361 0 : Initialize(T, ip, &normal);
362 0 : int attr = GetLocalVectorValue(E.Real(), V, true);
363 0 : if (attr <= 0)
364 : {
365 : return 0.0;
366 : }
367 0 : double Vn = V * normal;
368 0 : double Vn2 = Vn * Vn;
369 0 : if (E.HasImag())
370 : {
371 0 : GetLocalVectorValue(E.Imag(), V, true);
372 0 : Vn = V * normal;
373 0 : Vn2 += Vn * Vn;
374 : }
375 :
376 : // Metal-air interface: 0.5 * t / ε_MA * |E_n|² .
377 0 : return 0.5 * (t_i / epsilon_i) * Vn2;
378 : }
379 :
380 : template <>
381 0 : inline double InterfaceDielectricCoefficient<InterfaceDielectric::MS>::Eval(
382 : mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
383 : {
384 : // Get single-sided solution on substrate side and neighboring element attribute.
385 : double V_data[3], W_data[3], normal_data[3];
386 0 : mfem::Vector V(V_data, T.GetSpaceDim()), W(W_data, T.GetSpaceDim()),
387 0 : normal(normal_data, T.GetSpaceDim());
388 0 : Initialize(T, ip, &normal);
389 0 : int attr = GetLocalVectorValue(E.Real(), V, false);
390 0 : if (attr <= 0)
391 : {
392 : return 0.0;
393 : }
394 0 : mat_op.GetPermittivityReal(attr).Mult(V, W);
395 0 : double Vn = W * normal;
396 0 : double Vn2 = Vn * Vn;
397 0 : if (E.HasImag())
398 : {
399 0 : GetLocalVectorValue(E.Imag(), V, false);
400 0 : mat_op.GetPermittivityReal(attr).Mult(V, W);
401 0 : Vn = W * normal;
402 0 : Vn2 += Vn * Vn;
403 : }
404 :
405 : // Metal-substrate interface: 0.5 * t / ε_MS * |(ε_S E)_n|² .
406 0 : return 0.5 * (t_i / epsilon_i) * Vn2;
407 : }
408 :
409 : template <>
410 0 : inline double InterfaceDielectricCoefficient<InterfaceDielectric::SA>::Eval(
411 : mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
412 : {
413 : // Get single-sided solution on air side and neighboring element attribute.
414 : double V_data[3], normal_data[3];
415 0 : mfem::Vector V(V_data, T.GetSpaceDim()), normal(normal_data, T.GetSpaceDim());
416 0 : Initialize(T, ip, &normal);
417 0 : int attr = GetLocalVectorValue(E.Real(), V, true);
418 0 : if (attr <= 0)
419 : {
420 : return 0.0;
421 : }
422 0 : double Vn = V * normal;
423 0 : V.Add(-Vn, normal);
424 0 : double Vn2 = Vn * Vn;
425 0 : double Vt2 = V * V;
426 0 : if (E.HasImag())
427 : {
428 0 : GetLocalVectorValue(E.Imag(), V, true);
429 0 : Vn = V * normal;
430 0 : V.Add(-Vn, normal);
431 0 : Vn2 += Vn * Vn;
432 0 : Vt2 += V * V;
433 : }
434 :
435 : // Substrate-air interface: 0.5 * t * (ε_SA * |E_t|² + 1 / ε_SA * |E_n|²) .
436 0 : return 0.5 * t_i * ((epsilon_i * Vt2) + (Vn2 / epsilon_i));
437 : }
438 :
439 : // Helper for EnergyDensityCoefficient.
440 : enum class EnergyDensityType
441 : {
442 : ELECTRIC,
443 : MAGNETIC
444 : };
445 :
446 : // Returns the local energy density evaluated as 1/2 Dᴴ E or 1/2 Hᴴ B for real-valued
447 : // material coefficients. For internal boundary elements, the solution is averaged across
448 : // the interface.
449 : template <EnergyDensityType Type>
450 : class EnergyDensityCoefficient : public mfem::Coefficient, public BdrGridFunctionCoefficient
451 : {
452 : private:
453 : const GridFunction &U;
454 : const MaterialOperator &mat_op;
455 :
456 : double GetLocalEnergyDensity(mfem::ElementTransformation &T) const;
457 :
458 : public:
459 24 : EnergyDensityCoefficient(const GridFunction &U, const MaterialOperator &mat_op)
460 24 : : mfem::Coefficient(), BdrGridFunctionCoefficient(*U.ParFESpace()->GetParMesh()), U(U),
461 24 : mat_op(mat_op)
462 : {
463 : }
464 :
465 25920 : double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override
466 : {
467 25920 : if (T.ElementType == mfem::ElementTransformation::ELEMENT)
468 : {
469 20736 : return GetLocalEnergyDensity(T);
470 : }
471 5184 : else if (T.ElementType == mfem::ElementTransformation::BDR_ELEMENT)
472 : {
473 : // Get neighboring elements.
474 5184 : GetBdrElementNeighborTransformations(T.ElementNo, ip);
475 :
476 : // For interior faces, compute the average value.
477 5184 : if (FET.Elem2)
478 : {
479 : return 0.5 *
480 0 : (GetLocalEnergyDensity(*FET.Elem1) + GetLocalEnergyDensity(*FET.Elem2));
481 : }
482 : else
483 : {
484 5184 : return GetLocalEnergyDensity(*FET.Elem1);
485 : }
486 : }
487 0 : MFEM_ABORT("Unsupported element type in EnergyDensityCoefficient!");
488 : return 0.0;
489 : }
490 : };
491 :
492 : template <>
493 12960 : inline double EnergyDensityCoefficient<EnergyDensityType::ELECTRIC>::GetLocalEnergyDensity(
494 : mfem::ElementTransformation &T) const
495 : {
496 : // Only the real part of the permittivity contributes to the energy (imaginary part
497 : // cancels out in the inner product due to symmetry).
498 : double V_data[3];
499 12960 : mfem::Vector V(V_data, T.GetSpaceDim());
500 12960 : U.Real().GetVectorValue(T, T.GetIntPoint(), V);
501 25920 : double dot = mat_op.GetPermittivityReal(T.Attribute).InnerProduct(V, V);
502 12960 : if (U.HasImag())
503 : {
504 6480 : U.Imag().GetVectorValue(T, T.GetIntPoint(), V);
505 12960 : dot += mat_op.GetPermittivityReal(T.Attribute).InnerProduct(V, V);
506 : }
507 12960 : return 0.5 * dot;
508 : }
509 :
510 : template <>
511 12960 : inline double EnergyDensityCoefficient<EnergyDensityType::MAGNETIC>::GetLocalEnergyDensity(
512 : mfem::ElementTransformation &T) const
513 : {
514 : double V_data[3];
515 12960 : mfem::Vector V(V_data, T.GetSpaceDim());
516 12960 : U.Real().GetVectorValue(T, T.GetIntPoint(), V);
517 25920 : double dot = mat_op.GetInvPermeability(T.Attribute).InnerProduct(V, V);
518 12960 : if (U.HasImag())
519 : {
520 6480 : U.Imag().GetVectorValue(T, T.GetIntPoint(), V);
521 12960 : dot += mat_op.GetInvPermeability(T.Attribute).InnerProduct(V, V);
522 : }
523 12960 : return 0.5 * dot;
524 : }
525 :
526 : // Compute time-averaged Poynting vector Re{E x H⋆}, without the typical factor of 1/2. For
527 : // internal boundary elements, the solution is taken as the average.
528 : class PoyntingVectorCoefficient : public mfem::VectorCoefficient,
529 : public BdrGridFunctionCoefficient
530 : {
531 : private:
532 : const GridFunction &E, &B;
533 : const MaterialOperator &mat_op;
534 :
535 11664 : void GetLocalPower(mfem::ElementTransformation &T, mfem::Vector &V) const
536 : {
537 : double W1_data[3], W2_data[3];
538 23328 : mfem::Vector W1(W1_data, T.GetSpaceDim()), W2(W2_data, T.GetSpaceDim());
539 11664 : B.Real().GetVectorValue(T, T.GetIntPoint(), W1);
540 11664 : mat_op.GetInvPermeability(T.Attribute).Mult(W1, W2);
541 11664 : E.Real().GetVectorValue(T, T.GetIntPoint(), W1);
542 11664 : V.SetSize(vdim);
543 11664 : linalg::Cross3(W1, W2, V);
544 11664 : if (E.HasImag())
545 : {
546 7776 : B.Imag().GetVectorValue(T, T.GetIntPoint(), W1);
547 7776 : mat_op.GetInvPermeability(T.Attribute).Mult(W1, W2);
548 7776 : E.Imag().GetVectorValue(T, T.GetIntPoint(), W1);
549 7776 : linalg::Cross3(W1, W2, V, true);
550 : }
551 11664 : }
552 :
553 : public:
554 9 : PoyntingVectorCoefficient(const GridFunction &E, const GridFunction &B,
555 : const MaterialOperator &mat_op)
556 9 : : mfem::VectorCoefficient(E.VectorDim()),
557 9 : BdrGridFunctionCoefficient(*E.ParFESpace()->GetParMesh()), E(E), B(B), mat_op(mat_op)
558 : {
559 9 : }
560 :
561 : using mfem::VectorCoefficient::Eval;
562 11664 : void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
563 : const mfem::IntegrationPoint &ip) override
564 : {
565 11664 : if (T.ElementType == mfem::ElementTransformation::ELEMENT)
566 : {
567 9720 : GetLocalPower(T, V);
568 9720 : return;
569 : }
570 1944 : else if (T.ElementType == mfem::ElementTransformation::BDR_ELEMENT)
571 : {
572 : // Get neighboring elements.
573 1944 : GetBdrElementNeighborTransformations(T.ElementNo, ip);
574 :
575 : // For interior faces, compute the value on the desired side.
576 1944 : GetLocalPower(*FET.Elem1, V);
577 1944 : if (FET.Elem2)
578 : {
579 : double W_data[3];
580 : mfem::Vector W(W_data, V.Size());
581 0 : GetLocalPower(*FET.Elem2, W);
582 0 : add(0.5, V, W, V);
583 : }
584 1944 : return;
585 : }
586 0 : MFEM_ABORT("Unsupported element type in PoyntingVectorCoefficient!");
587 : }
588 : };
589 :
590 : // Returns the local vector field evaluated on a boundary element. For internal boundary
591 : // elements the solution is the average.
592 : class BdrFieldVectorCoefficient : public mfem::VectorCoefficient,
593 : public BdrGridFunctionCoefficient
594 : {
595 : private:
596 : const mfem::ParGridFunction &U;
597 :
598 : public:
599 39 : BdrFieldVectorCoefficient(const mfem::ParGridFunction &U)
600 39 : : mfem::VectorCoefficient(U.VectorDim()),
601 39 : BdrGridFunctionCoefficient(*U.ParFESpace()->GetParMesh()), U(U)
602 : {
603 39 : }
604 :
605 : using mfem::VectorCoefficient::Eval;
606 8424 : void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
607 : const mfem::IntegrationPoint &ip) override
608 : {
609 : // Get neighboring elements.
610 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
611 : "Unexpected element type in BdrFieldVectorCoefficient!");
612 8424 : GetBdrElementNeighborTransformations(T.ElementNo, ip);
613 :
614 : // For interior faces, compute the average.
615 8424 : U.GetVectorValue(*FET.Elem1, FET.Elem1->GetIntPoint(), V);
616 8424 : if (FET.Elem2)
617 : {
618 : double W_data[3];
619 : mfem::Vector W(W_data, V.Size());
620 0 : U.GetVectorValue(*FET.Elem2, FET.Elem2->GetIntPoint(), W);
621 0 : add(0.5, V, W, V);
622 : }
623 8424 : }
624 : };
625 :
626 : // Returns the local scalar field evaluated on a boundary element. For internal boundary
627 : // elements the solution is the average.
628 : class BdrFieldCoefficient : public mfem::Coefficient, public BdrGridFunctionCoefficient
629 : {
630 : private:
631 : const mfem::ParGridFunction &U;
632 :
633 : public:
634 : BdrFieldCoefficient(const mfem::ParGridFunction &U)
635 3 : : mfem::Coefficient(), BdrGridFunctionCoefficient(*U.ParFESpace()->GetParMesh()), U(U)
636 : {
637 : }
638 :
639 648 : double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override
640 : {
641 : // Get neighboring elements.
642 : MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
643 : "Unexpected element type in BdrFieldCoefficient!");
644 648 : GetBdrElementNeighborTransformations(T.ElementNo, ip);
645 :
646 : // For interior faces, compute the average.
647 648 : if (FET.Elem2)
648 : {
649 0 : return 0.5 * (U.GetValue(*FET.Elem1, FET.Elem1->GetIntPoint()),
650 0 : U.GetValue(*FET.Elem2, FET.Elem2->GetIntPoint()));
651 : }
652 : else
653 : {
654 648 : return U.GetValue(*FET.Elem1, FET.Elem1->GetIntPoint());
655 : }
656 : }
657 : };
658 :
659 : //
660 : // More helpful coefficient types. Wrapper coefficients allow additions of scalar and vector
661 : // or matrix coefficients. Restricted coefficients only compute the coefficient if for the
662 : // given list of attributes. Sum coefficients own a list of coefficients to add.
663 : //
664 :
665 : class VectorWrappedCoefficient : public mfem::VectorCoefficient
666 : {
667 : private:
668 : std::unique_ptr<mfem::Coefficient> coeff;
669 :
670 : public:
671 : VectorWrappedCoefficient(int dim, std::unique_ptr<mfem::Coefficient> &&coeff)
672 : : mfem::VectorCoefficient(dim), coeff(std::move(coeff))
673 : {
674 : }
675 :
676 : using mfem::VectorCoefficient::Eval;
677 0 : void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
678 : const mfem::IntegrationPoint &ip) override
679 : {
680 0 : V.SetSize(vdim);
681 0 : V = coeff->Eval(T, ip);
682 0 : }
683 : };
684 :
685 : class MatrixWrappedCoefficient : public mfem::MatrixCoefficient
686 : {
687 : private:
688 : std::unique_ptr<mfem::Coefficient> coeff;
689 :
690 : public:
691 : MatrixWrappedCoefficient(int dim, std::unique_ptr<mfem::Coefficient> &&coeff)
692 : : mfem::MatrixCoefficient(dim), coeff(std::move(coeff))
693 : {
694 : }
695 :
696 : void Eval(mfem::DenseMatrix &K, mfem::ElementTransformation &T,
697 : const mfem::IntegrationPoint &ip) override
698 : {
699 : K.Diag(coeff->Eval(T, ip), height);
700 : }
701 : };
702 :
703 : template <typename Coefficient>
704 : class RestrictedCoefficient : public Coefficient
705 : {
706 : private:
707 : mfem::Array<int> attr_marker;
708 :
709 : public:
710 : template <typename... T>
711 0 : RestrictedCoefficient(const mfem::Array<int> &attr_list, T &&...args)
712 : : Coefficient(std::forward<T>(args)...),
713 0 : attr_marker(mesh::AttrToMarker(attr_list.Size() ? attr_list.Max() : 0, attr_list))
714 : {
715 0 : }
716 :
717 0 : double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override
718 : {
719 0 : return (T.Attribute > attr_marker.Size() || !attr_marker[T.Attribute - 1])
720 0 : ? 0.0
721 0 : : Coefficient::Eval(T, ip);
722 : }
723 : };
724 :
725 : template <typename Coefficient>
726 : class RestrictedVectorCoefficient : public Coefficient
727 : {
728 : private:
729 : mfem::Array<int> attr_marker;
730 :
731 : public:
732 : template <typename... T>
733 21 : RestrictedVectorCoefficient(const mfem::Array<int> &attr_list, T &&...args)
734 : : Coefficient(std::forward<T>(args)...),
735 21 : attr_marker(mesh::AttrToMarker(attr_list.Size() ? attr_list.Max() : 0, attr_list))
736 : {
737 21 : }
738 :
739 756 : void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
740 : const mfem::IntegrationPoint &ip) override
741 : {
742 756 : if (T.Attribute > attr_marker.Size() || !attr_marker[T.Attribute - 1])
743 : {
744 0 : V.SetSize(this->vdim);
745 0 : V = 0.0;
746 : }
747 : else
748 : {
749 324 : Coefficient::Eval(V, T, ip);
750 : }
751 756 : }
752 : };
753 :
754 : template <typename Coefficient>
755 : class RestrictedMatrixCoefficient : public Coefficient
756 : {
757 : private:
758 : mfem::Array<int> attr_marker;
759 :
760 : public:
761 : template <typename... T>
762 : RestrictedMatrixCoefficient(const mfem::Array<int> &attr_list, T &&...args)
763 : : Coefficient(std::forward<T>(args)...),
764 : attr_marker(mesh::AttrToMarker(attr_list.Size() ? attr_list.Max() : 0, attr_list))
765 : {
766 : }
767 :
768 : void Eval(mfem::DenseMatrix &K, mfem::ElementTransformation &T,
769 : const mfem::IntegrationPoint &ip) override
770 : {
771 : if (T.Attribute > attr_marker.Size() || !attr_marker[T.Attribute - 1])
772 : {
773 : K.SetSize(this->height, this->width);
774 : K = 0.0;
775 : }
776 : else
777 : {
778 : Coefficient::Eval(K, T, ip);
779 : }
780 : }
781 : };
782 :
783 : class SumCoefficient : public mfem::Coefficient
784 : {
785 : private:
786 : std::vector<std::pair<std::unique_ptr<mfem::Coefficient>, double>> c;
787 :
788 : public:
789 : SumCoefficient() : mfem::Coefficient() {}
790 :
791 : bool empty() const { return c.empty(); }
792 :
793 : void AddCoefficient(std::unique_ptr<mfem::Coefficient> &&coeff, double a = 1.0)
794 : {
795 : c.emplace_back(std::move(coeff), a);
796 : }
797 :
798 0 : double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override
799 : {
800 : double val = 0.0;
801 0 : for (auto &[coeff, a] : c)
802 : {
803 0 : val += a * coeff->Eval(T, ip);
804 : }
805 0 : return val;
806 : }
807 : };
808 :
809 18 : class SumVectorCoefficient : public mfem::VectorCoefficient
810 : {
811 : private:
812 : std::vector<std::pair<std::unique_ptr<mfem::VectorCoefficient>, double>> c;
813 :
814 : public:
815 18 : SumVectorCoefficient(int d) : mfem::VectorCoefficient(d) {}
816 :
817 : bool empty() const { return c.empty(); }
818 :
819 21 : void AddCoefficient(std::unique_ptr<mfem::VectorCoefficient> &&coeff, double a = 1.0)
820 : {
821 21 : MFEM_VERIFY(coeff->GetVDim() == vdim,
822 : "Invalid VectorCoefficient dimensions for SumVectorCoefficient!");
823 21 : c.emplace_back(std::move(coeff), a);
824 21 : }
825 :
826 : void AddCoefficient(std::unique_ptr<mfem::Coefficient> &&coeff, double a = 1.0)
827 : {
828 : c.emplace_back(std::make_unique<VectorWrappedCoefficient>(vdim, std::move(coeff)), a);
829 : }
830 :
831 : using mfem::VectorCoefficient::Eval;
832 756 : void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
833 : const mfem::IntegrationPoint &ip) override
834 : {
835 : double U_data[3];
836 756 : mfem::Vector U(U_data, vdim);
837 756 : V.SetSize(vdim);
838 756 : V = 0.0;
839 1512 : for (auto &[coeff, a] : c)
840 : {
841 756 : coeff->Eval(U, T, ip);
842 756 : V.Add(a, U);
843 : }
844 756 : }
845 : };
846 :
847 : class SumMatrixCoefficient : public mfem::MatrixCoefficient
848 : {
849 : private:
850 : std::vector<std::pair<std::unique_ptr<mfem::MatrixCoefficient>, double>> c;
851 :
852 : public:
853 : SumMatrixCoefficient(int d) : mfem::MatrixCoefficient(d) {}
854 : SumMatrixCoefficient(int h, int w) : mfem::MatrixCoefficient(h, w) {}
855 :
856 : bool empty() const { return c.empty(); }
857 :
858 : void AddCoefficient(std::unique_ptr<mfem::MatrixCoefficient> &&coeff, double a)
859 : {
860 : MFEM_VERIFY(coeff->GetHeight() == height && coeff->GetWidth() == width,
861 : "Invalid MatrixCoefficient dimensions for SumMatrixCoefficient!");
862 : c.emplace_back(std::move(coeff), a);
863 : }
864 :
865 : void AddCoefficient(std::unique_ptr<mfem::Coefficient> &&coeff, double a)
866 : {
867 : MFEM_VERIFY(width == height, "MatrixWrappedCoefficient can only be constructed for "
868 : "square MatrixCoefficient objects!");
869 : c.emplace_back(std::make_unique<MatrixWrappedCoefficient>(height, std::move(coeff)), a);
870 : }
871 :
872 : void Eval(mfem::DenseMatrix &K, mfem::ElementTransformation &T,
873 : const mfem::IntegrationPoint &ip) override
874 : {
875 : double M_data[9];
876 : mfem::DenseMatrix M(M_data, height, width);
877 : K.SetSize(height, width);
878 : K = 0.0;
879 : for (auto &[coeff, a] : c)
880 : {
881 : coeff->Eval(M, T, ip);
882 : K.Add(a, M);
883 : }
884 : }
885 : };
886 :
887 : } // namespace palace
888 :
889 : #endif // PALACE_FEM_COEFFICIENT_HPP
|
