libMesh
Loading...
Searching...
No Matches
hdiv_fe_transformation.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4// This library is free software; you can redistribute it and/or
5// modify it under the terms of the GNU Lesser General Public
6// License as published by the Free Software Foundation; either
7// version 2.1 of the License, or (at your option) any later version.
8
9// This library is distributed in the hope that it will be useful,
10// but WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12// Lesser General Public License for more details.
13
14// You should have received a copy of the GNU Lesser General Public
15// License along with this library; if not, write to the Free Software
16// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18#include "libmesh/hdiv_fe_transformation.h"
19#include "libmesh/elem.h"
20#include "libmesh/fe_interface.h"
21#include "libmesh/int_range.h"
22#include "libmesh/tensor_value.h"
23
24namespace {
25
26using namespace libMesh;
27
28// The contravariant Piola map phi = J^{-1} * (dx/dxi) * phihat, shared by
29// map_phi and map_dphi (which also needs the physical shape function value
30// for the term coming from the derivative of J^{-1}).
31template<typename OutputShape>
32OutputShape hdiv_piola_map(const RealGradient & dxyz_dxi,
33 const RealGradient & dxyz_deta,
34 Real J,
35 const OutputShape & phi_ref)
36{
37 return (dxyz_dxi*phi_ref(0) + dxyz_deta*phi_ref(1))/J;
38}
39
40template<typename OutputShape>
41OutputShape hdiv_piola_map(const RealGradient & dxyz_dxi,
42 const RealGradient & dxyz_deta,
43 const RealGradient & dxyz_dzeta,
44 Real J,
45 const OutputShape & phi_ref)
46{
47 return (dxyz_dxi*phi_ref(0) + dxyz_deta*phi_ref(1) + dxyz_dzeta*phi_ref(2))/J;
48}
49
50} // anonymous namespace
51
52namespace libMesh
53{
54
55template<typename OutputShape>
56void HDivFETransformation<OutputShape>::init_map_phi(const FEGenericBase<OutputShape> & fe) const
57{
58 // We only need to pre-request one first-order derivative map piece because there is only a
59 // single calculate_dxyz boolean flag that toggles first derivative computations in the FEMap
60 // and in addition to covering all components it also covers both forward and inverse mapping.
61 // We choose to document pre-requesting the forward (reference -> physical) map here because
62 // that is what is used in map_phi (there is no inverse map usage).
63 fe.get_fe_map().get_dxyzdxi();
64}
65
66
67
68template<typename OutputShape>
69void HDivFETransformation<OutputShape>::init_map_dphi(const FEGenericBase<OutputShape> & fe) const
70{
71 // See above comment in init_map_phi. In map_dphi we actually use both the forward and inverse
72 // first derivative maps so the choice here is a little more arbitrary. We choose to be
73 // consistent with init_map_phi.
74 fe.get_fe_map().get_dxyzdxi();
75#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
76 // As for first order derivatives, there is only a single boolean flag (calculate_d2xyz)
77 // controlling second order derivative computations in FEMap so we only bother prerequesting
78 // one piece. We document pre-requesting the forward map because those are the second
79 // derivatives used in map_dphi.
80 fe.get_fe_map().get_d2xyzdxi2();
81#endif
82}
83
84
85
86template<typename OutputShape>
87void HDivFETransformation<OutputShape>::init_map_d2phi(const FEGenericBase<OutputShape> & /*fe*/) const
88{
89 // We choose not to pre-request any computations here since we have not yet implemented
90 // map_d2phi. We don't need to give the map unnecessary work.
91}
92
93
94
95template<typename OutputShape>
96void HDivFETransformation<OutputShape>::map_phi(const unsigned int dim,
97 const Elem * const elem,
98 const std::vector<Point> & qp,
99 const FEGenericBase<OutputShape> & fe,
100 std::vector<std::vector<OutputShape>> & phi,
101 const bool /*add_p_level*/) const
102{
103 switch (dim)
104 {
105 case 0:
106 case 1:
107 libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
108
109 case 2:
110 {
111 const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
112 const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
113
114 const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
115
116 // phi = J^{-1} * (dx/dxi) * \hat{phi}
117 for (auto i : index_range(phi))
118 for (auto p : index_range(phi[i]))
119 {
120 // Need to temporarily cache reference shape functions
121 // We are computing mapping basis functions, so we explicitly ignore
122 // any non-zero p_level() the Elem might have.
123 OutputShape phi_ref;
124 FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
125
126 phi[i][p] = hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], J[p], phi_ref);
127 }
128
129 break;
130 }
131 case 3:
132 {
133 const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
134 const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
135 const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
136
137 const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
138
139 // phi = J^{-1} * (dx/dxi) * \hat{phi}
140 for (auto i : index_range(phi))
141 for (auto p : index_range(phi[i]))
142 {
143 // Need to temporarily cache reference shape functions
144 // We are computing mapping basis functions, so we explicitly ignore
145 // any non-zero p_level() the Elem might have.
146 OutputShape phi_ref;
147 FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
148
149 phi[i][p] = hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], dxyz_dzeta[p], J[p], phi_ref);
150 }
151
152 break;
153 }
154
155 default:
156 libmesh_error_msg("Invalid dim = " << dim);
157 } // switch(dim)
158}
159
160template<typename OutputShape>
161void HDivFETransformation<OutputShape>::map_dphi(const unsigned int dim,
162 const Elem * const elem,
163 const std::vector<Point> & qp,
164 const FEGenericBase<OutputShape> & fe,
165 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi,
166 std::vector<std::vector<OutputShape>> & dphidx,
167 std::vector<std::vector<OutputShape>> & dphidy,
168 std::vector<std::vector<OutputShape>> & dphidz) const
169{
170 switch (dim)
171 {
172 case 0:
173 case 1:
174 libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
175
176 case 2:
177 {
178 const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
179 const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
180
181#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
182 const std::vector<RealGradient> & d2xyz_dxi2 = fe.get_fe_map().get_d2xyzdxi2();
183 const std::vector<RealGradient> & d2xyz_deta2 = fe.get_fe_map().get_d2xyzdeta2();
184 const std::vector<RealGradient> & d2xyz_dxideta = fe.get_fe_map().get_d2xyzdxideta();
185#else
186 // Anything but a triangle mesh or simple grid may need --enable-second
187 libmesh_error_msg_if
188 (!elem->has_affine_map(),
189 "HDiv FE gradients on non-affine-mapped elements require second-derivative support");
190 auto n_qp = dxyz_dxi.size();
191
192 const std::vector<RealGradient> d2xyz_dxi2(n_qp, RealGradient());
193 const std::vector<RealGradient> d2xyz_deta2(n_qp, RealGradient());
194 const std::vector<RealGradient> d2xyz_dxideta(n_qp, RealGradient());
195#endif
196
197 const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
198
199 const std::vector<Real> & dxi_dx = fe.get_fe_map().get_dxidx();
200 const std::vector<Real> & dxi_dy = fe.get_fe_map().get_dxidy();
201 const std::vector<Real> & deta_dx = fe.get_fe_map().get_detadx();
202 const std::vector<Real> & deta_dy = fe.get_fe_map().get_detady();
203#if LIBMESH_DIM > 2
204 const std::vector<Real> & dxi_dz = fe.get_fe_map().get_dxidz();
205 const std::vector<Real> & deta_dz = fe.get_fe_map().get_detadz();
206#endif
207
208 const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
209 const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
210
211 for (auto i : index_range(dphi))
212 for (auto p : index_range(dphi[i]))
213 {
214 // Need to temporarily cache reference shape functions
215 // We are computing mapping basis functions, so we explicitly ignore
216 // any non-zero p_level() the Elem might have.
217 OutputShape phi_ref;
218 FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
219
220 const OutputShape phi_val =
221 hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], J[p], phi_ref);
222
223 // Jacobi's formula: dJ/dxi_n = J * tr(F^{-1} * dF/dxi_n)
224 Real dJ_dxi =
225 J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) +
226 dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1));
227 Real dJ_deta =
228 J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) +
229 dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1));
230#if LIBMESH_DIM > 2
231 dJ_dxi += J[p] * (dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2));
232 dJ_deta += J[p] * (dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2));
233#endif
234
235 for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
236 {
237 // dphi_k/dx_l = A + B + C, where (n, m summed over reference directions):
238 // A = -phi_k(x) * SUM_n (dxi_n/dx_l) * (1/J) * (dJ/dxi_n)
239 // B = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m
240 // C = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n
241
242 // Term A: -phi_k(x) * SUM_n (dxi_n/dx_l)*(1/J)*(dJ/dxi_n)
243 const Real phik = phi_val(k);
244
245 // Term B: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m
246 const Real d2xk_dxi2 = d2xyz_dxi2[p](k);
247 const Real d2xk_deta2 = d2xyz_deta2[p](k);
248 const Real d2xk_dxideta = d2xyz_dxideta[p](k);
249
250 const Real dphik_dxi_B = (d2xk_dxi2*phi_ref(0) + d2xk_dxideta*phi_ref(1))/J[p];
251 const Real dphik_deta_B = (d2xk_dxideta*phi_ref(0) + d2xk_deta2*phi_ref(1))/J[p];
252
253 // Term C: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n
254 const Real Fk_xi = dxyz_dxi[p](k);
255 const Real Fk_eta = dxyz_deta[p](k);
256
257 const Real dphik_dxi_C = (Fk_xi*dphi_dxi[i][p](0) + Fk_eta*dphi_dxi[i][p](1))/J[p];
258 const Real dphik_deta_C = (Fk_xi*dphi_deta[i][p](0) + Fk_eta*dphi_deta[i][p](1))/J[p];
259
260 const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C;
261 const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C;
262
263 dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total;
264 dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total;
265#if LIBMESH_DIM > 2
266 dphidz[i][p](k) = dxi_dz[p]*dphik_dxi_total + deta_dz[p]*dphik_deta_total;
267#endif
268 }
269
270 dphi[i][p].slice(0) = dphidx[i][p];
271 dphi[i][p].slice(1) = dphidy[i][p];
272#if LIBMESH_DIM > 2
273 dphi[i][p].slice(2) = dphidz[i][p];
274#endif
275 }
276
277 break;
278 }
279 case 3:
280 {
281 const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
282 const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
283 const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
284
285#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
286 const std::vector<RealGradient> & d2xyz_dxi2 = fe.get_fe_map().get_d2xyzdxi2();
287 const std::vector<RealGradient> & d2xyz_deta2 = fe.get_fe_map().get_d2xyzdeta2();
288 const std::vector<RealGradient> & d2xyz_dzeta2 = fe.get_fe_map().get_d2xyzdzeta2();
289 const std::vector<RealGradient> & d2xyz_dxideta = fe.get_fe_map().get_d2xyzdxideta();
290 const std::vector<RealGradient> & d2xyz_dxidzeta = fe.get_fe_map().get_d2xyzdxidzeta();
291 const std::vector<RealGradient> & d2xyz_detadzeta = fe.get_fe_map().get_d2xyzdetadzeta();
292#else
293 // Anything but a tet mesh or simple grid may need --enable-second
294 libmesh_error_msg_if
295 (!elem->has_affine_map(),
296 "HDiv FE gradients on non-affine-mapped elements require second-derivative support");
297 auto n_qp = dxyz_dxi.size();
298
299 const std::vector<RealGradient> d2xyz_dxi2(n_qp, RealGradient());
300 const std::vector<RealGradient> d2xyz_deta2(n_qp, RealGradient());
301 const std::vector<RealGradient> d2xyz_dzeta2(n_qp, RealGradient());
302 const std::vector<RealGradient> d2xyz_dxideta(n_qp, RealGradient());
303 const std::vector<RealGradient> d2xyz_dxidzeta(n_qp, RealGradient());
304 const std::vector<RealGradient> d2xyz_detadzeta(n_qp, RealGradient());
305#endif
306
307 const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
308
309 const std::vector<Real> & dxi_dx = fe.get_fe_map().get_dxidx();
310 const std::vector<Real> & dxi_dy = fe.get_fe_map().get_dxidy();
311 const std::vector<Real> & dxi_dz = fe.get_fe_map().get_dxidz();
312 const std::vector<Real> & deta_dx = fe.get_fe_map().get_detadx();
313 const std::vector<Real> & deta_dy = fe.get_fe_map().get_detady();
314 const std::vector<Real> & deta_dz = fe.get_fe_map().get_detadz();
315 const std::vector<Real> & dzeta_dx = fe.get_fe_map().get_dzetadx();
316 const std::vector<Real> & dzeta_dy = fe.get_fe_map().get_dzetady();
317 const std::vector<Real> & dzeta_dz = fe.get_fe_map().get_dzetadz();
318
319 const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
320 const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
321 const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
322
323 for (auto i : index_range(dphi))
324 for (auto p : index_range(dphi[i]))
325 {
326 // Need to temporarily cache reference shape functions
327 // We are computing mapping basis functions, so we explicitly ignore
328 // any non-zero p_level() the Elem might have.
329 OutputShape phi_ref;
330 FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
331
332 const OutputShape phi_val =
333 hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], dxyz_dzeta[p], J[p], phi_ref);
334
335 // Jacobi's formula: dJ/dxi_n = J * tr(F^{-1} * dF/dxi_n)
336 const Real dJ_dxi =
337 J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) + dzeta_dx[p]*d2xyz_dxidzeta[p](0) +
338 dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1) + dzeta_dy[p]*d2xyz_dxidzeta[p](1) +
339 dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2) + dzeta_dz[p]*d2xyz_dxidzeta[p](2));
340 const Real dJ_deta =
341 J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) + dzeta_dx[p]*d2xyz_detadzeta[p](0) +
342 dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1) + dzeta_dy[p]*d2xyz_detadzeta[p](1) +
343 dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2) + dzeta_dz[p]*d2xyz_detadzeta[p](2));
344 const Real dJ_dzeta =
345 J[p] * (dxi_dx[p]*d2xyz_dxidzeta[p](0) + deta_dx[p]*d2xyz_detadzeta[p](0) + dzeta_dx[p]*d2xyz_dzeta2[p](0) +
346 dxi_dy[p]*d2xyz_dxidzeta[p](1) + deta_dy[p]*d2xyz_detadzeta[p](1) + dzeta_dy[p]*d2xyz_dzeta2[p](1) +
347 dxi_dz[p]*d2xyz_dxidzeta[p](2) + deta_dz[p]*d2xyz_detadzeta[p](2) + dzeta_dz[p]*d2xyz_dzeta2[p](2));
348
349 for (unsigned int k = 0; k < 3; ++k)
350 {
351 // dphi_k/dx_l = A + B + C, where (n, m summed over reference directions):
352 // A = -phi_k(x) * SUM_n (dxi_n/dx_l) * (1/J) * (dJ/dxi_n)
353 // B = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m
354 // C = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n
355
356 // Term A: -phi_k(x) * SUM_n (dxi_n/dx_l)*(1/J)*(dJ/dxi_n)
357 const Real phik = phi_val(k);
358
359 // Term B: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m
360 const Real d2xk_dxi2 = d2xyz_dxi2[p](k);
361 const Real d2xk_deta2 = d2xyz_deta2[p](k);
362 const Real d2xk_dzeta2 = d2xyz_dzeta2[p](k);
363 const Real d2xk_dxideta = d2xyz_dxideta[p](k);
364 const Real d2xk_dxidzeta = d2xyz_dxidzeta[p](k);
365 const Real d2xk_detadzeta = d2xyz_detadzeta[p](k);
366
367 const Real dphik_dxi_B = (d2xk_dxi2*phi_ref(0) + d2xk_dxideta*phi_ref(1) + d2xk_dxidzeta*phi_ref(2))/J[p];
368 const Real dphik_deta_B = (d2xk_dxideta*phi_ref(0) + d2xk_deta2*phi_ref(1) + d2xk_detadzeta*phi_ref(2))/J[p];
369 const Real dphik_dzeta_B = (d2xk_dxidzeta*phi_ref(0) + d2xk_detadzeta*phi_ref(1) + d2xk_dzeta2*phi_ref(2))/J[p];
370
371 // Term C: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n
372 const Real Fk_xi = dxyz_dxi[p](k);
373 const Real Fk_eta = dxyz_deta[p](k);
374 const Real Fk_zeta = dxyz_dzeta[p](k);
375
376 const Real dphik_dxi_C = (Fk_xi*dphi_dxi[i][p](0) + Fk_eta*dphi_dxi[i][p](1) + Fk_zeta*dphi_dxi[i][p](2))/J[p];
377 const Real dphik_deta_C = (Fk_xi*dphi_deta[i][p](0) + Fk_eta*dphi_deta[i][p](1) + Fk_zeta*dphi_deta[i][p](2))/J[p];
378 const Real dphik_dzeta_C = (Fk_xi*dphi_dzeta[i][p](0) + Fk_eta*dphi_dzeta[i][p](1) + Fk_zeta*dphi_dzeta[i][p](2))/J[p];
379
380 const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C;
381 const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C;
382 const Real dphik_dzeta_total = -phik*dJ_dzeta/J[p] + dphik_dzeta_B + dphik_dzeta_C;
383
384 dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total + dzeta_dx[p]*dphik_dzeta_total;
385 dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total + dzeta_dy[p]*dphik_dzeta_total;
386 dphidz[i][p](k) = dxi_dz[p]*dphik_dxi_total + deta_dz[p]*dphik_deta_total + dzeta_dz[p]*dphik_dzeta_total;
387 }
388
389 dphi[i][p].slice(0) = dphidx[i][p];
390 dphi[i][p].slice(1) = dphidy[i][p];
391 dphi[i][p].slice(2) = dphidz[i][p];
392 }
393
394 break;
395 }
396
397 default:
398 libmesh_error_msg("Invalid dim = " << dim);
399 } // switch(dim)
400}
401
402template<typename OutputShape>
403void HDivFETransformation<OutputShape>::map_div(const unsigned int dim,
404 const Elem * const,
405 const std::vector<Point> &,
406 const FEGenericBase<OutputShape> & fe,
407 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence>> & div_phi) const
408{
409 switch (dim)
410 {
411 case 0:
412 case 1:
413 libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
414
415 case 2:
416 {
417 const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
418 const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
419
420 const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
421
422 // div(phi) = J^{-1} * div(\hat{phi})
423 for (auto i : index_range(div_phi))
424 for (auto p : index_range(div_phi[i]))
425 {
426 div_phi[i][p] = (dphi_dxi[i][p](0) + dphi_deta[i][p](1))/J[p];
427 }
428
429 break;
430 }
431 case 3:
432 {
433 const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
434 const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
435 const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
436
437 const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
438
439 // div(phi) = J^{-1} * div(\hat{phi})
440 for (auto i : index_range(div_phi))
441 for (auto p : index_range(div_phi[i]))
442 {
443 div_phi[i][p] = (dphi_dxi[i][p](0) + dphi_deta[i][p](1) + dphi_dzeta[i][p](2))/J[p];
444 }
445
446 break;
447 }
448
449 default:
450 libmesh_error_msg("Invalid dim = " << dim);
451 } // switch(dim)
452}
453
454template class LIBMESH_EXPORT HDivFETransformation<RealGradient>;
455
456template<>
457void HDivFETransformation<Real>::init_map_phi(const FEGenericBase<Real> & ) const
458{
459 libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
460}
461
462template<>
463void HDivFETransformation<Real>::init_map_dphi(const FEGenericBase<Real> & ) const
464{
465 libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
466}
467
468template<>
469void HDivFETransformation<Real>::init_map_d2phi(const FEGenericBase<Real> & ) const
470{
471 libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
472}
473
474template<>
475void HDivFETransformation<Real>::map_phi(const unsigned int,
476 const Elem * const,
477 const std::vector<Point> &,
478 const FEGenericBase<Real> &,
479 std::vector<std::vector<Real>> &,
480 bool) const
481{
482 libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
483}
484
485template<>
486void HDivFETransformation<Real>::map_dphi(const unsigned int,
487 const Elem * const,
488 const std::vector<Point> &,
489 const FEGenericBase<Real> &,
490 std::vector<std::vector<FEGenericBase<Real>::OutputGradient>> &,
491 std::vector<std::vector<Real>> &,
492 std::vector<std::vector<Real>> &,
493 std::vector<std::vector<Real>> &) const
494{
495 libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
496}
497
498template<>
499void HDivFETransformation<Real>::map_div(const unsigned int,
500 const Elem * const,
501 const std::vector<Point> &,
502 const FEGenericBase<Real> &,
503 std::vector<std::vector<Real>> &) const
504{
505 libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
506}
507
508
509} // namespace libMesh
dim
unsigned int dim
Definition adaptivity_ex3.C:114
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::has_affine_map
virtual bool has_affine_map() const
Definition elem.h:1312
libMesh::FEAbstract::get_fe_type
FEType get_fe_type() const
Definition fe_abstract.h:520
libMesh::FEAbstract::get_fe_map
const FEMap & get_fe_map() const
Definition fe_abstract.h:554
libMesh::FEGenericBase
This class forms the foundation from which generic finite elements may be derived.
Definition fe_base.h:86
libMesh::FEGenericBase::get_dphidxi
const std::vector< std::vector< OutputShape > > & get_dphidxi() const
Definition fe_base.h:325
libMesh::FEGenericBase::get_dphidzeta
const std::vector< std::vector< OutputShape > > & get_dphidzeta() const
Definition fe_base.h:341
libMesh::FEGenericBase::OutputDivergence
TensorTools::DecrementRank< OutputShape >::type OutputDivergence
Definition fe_base.h:122
libMesh::FEGenericBase::get_dphideta
const std::vector< std::vector< OutputShape > > & get_dphideta() const
Definition fe_base.h:333
libMesh::FEInterface::shape
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition fe_interface.C:764
libMesh::FEMap::get_dxidz
const std::vector< Real > & get_dxidz() const
Definition fe_map.h:325
libMesh::FEMap::get_dxyzdeta
const std::vector< RealGradient > & get_dxyzdeta() const
Definition fe_map.h:247
libMesh::FEMap::get_detadx
const std::vector< Real > & get_detadx() const
Definition fe_map.h:333
libMesh::FEMap::get_jacobian
const std::vector< Real > & get_jacobian() const
Definition fe_map.h:223
libMesh::FEMap::get_dzetady
const std::vector< Real > & get_dzetady() const
Definition fe_map.h:365
libMesh::FEMap::get_detady
const std::vector< Real > & get_detady() const
Definition fe_map.h:341
libMesh::FEMap::get_d2xyzdxidzeta
const std::vector< RealGradient > & get_d2xyzdxidzeta() const
Definition fe_map.h:292
libMesh::FEMap::get_d2xyzdetadzeta
const std::vector< RealGradient > & get_d2xyzdetadzeta() const
Definition fe_map.h:299
libMesh::FEMap::get_dxidx
const std::vector< Real > & get_dxidx() const
Definition fe_map.h:309
libMesh::FEMap::get_d2xyzdeta2
const std::vector< RealGradient > & get_d2xyzdeta2() const
Definition fe_map.h:271
libMesh::FEMap::get_dzetadx
const std::vector< Real > & get_dzetadx() const
Definition fe_map.h:357
libMesh::FEMap::get_dxyzdzeta
const std::vector< RealGradient > & get_dxyzdzeta() const
Definition fe_map.h:255
libMesh::FEMap::get_dxidy
const std::vector< Real > & get_dxidy() const
Definition fe_map.h:317
libMesh::FEMap::get_d2xyzdzeta2
const std::vector< RealGradient > & get_d2xyzdzeta2() const
Definition fe_map.h:278
libMesh::FEMap::get_d2xyzdxi2
const std::vector< RealGradient > & get_d2xyzdxi2() const
Definition fe_map.h:264
libMesh::FEMap::get_dxyzdxi
const std::vector< RealGradient > & get_dxyzdxi() const
Definition fe_map.h:239
libMesh::FEMap::get_dzetadz
const std::vector< Real > & get_dzetadz() const
Definition fe_map.h:373
libMesh::FEMap::get_detadz
const std::vector< Real > & get_detadz() const
Definition fe_map.h:349
libMesh::FEMap::get_d2xyzdxideta
const std::vector< RealGradient > & get_d2xyzdxideta() const
Definition fe_map.h:285
libMesh::HDivFETransformation
This class handles the computation of the shape functions in the physical domain for HDiv conforming ...
Definition hdiv_fe_transformation.h:37
libMesh::HDivFETransformation::map_dphi
virtual void map_dphi(const unsigned int dim, const Elem *const elem, const std::vector< Point > &qp, const FEGenericBase< OutputShape > &fe, std::vector< std::vector< typename FEGenericBase< OutputShape >::OutputGradient > > &dphi, std::vector< std::vector< OutputShape > > &dphidx, std::vector< std::vector< OutputShape > > &dphidy, std::vector< std::vector< OutputShape > > &dphidz) const override
Evaluates shape function gradients in physical coordinates for conforming elements,...
Definition hdiv_fe_transformation.C:161
libMesh::HDivFETransformation::init_map_d2phi
virtual void init_map_d2phi(const FEGenericBase< OutputShape > &fe) const override
Pre-requests any necessary data from FEMap.
Definition hdiv_fe_transformation.C:87
libMesh::HDivFETransformation::init_map_dphi
virtual void init_map_dphi(const FEGenericBase< OutputShape > &fe) const override
Pre-requests any necessary data from FEMap.
Definition hdiv_fe_transformation.C:69
libMesh::HDivFETransformation::init_map_phi
virtual void init_map_phi(const FEGenericBase< OutputShape > &fe) const override
Pre-requests any necessary data from FEMap.
Definition hdiv_fe_transformation.C:56
libMesh::HDivFETransformation::map_div
virtual void map_div(const unsigned int dim, const Elem *const elem, const std::vector< Point > &qp, const FEGenericBase< OutputShape > &fe, std::vector< std::vector< typename FEGenericBase< OutputShape >::OutputDivergence > > &div_phi) const override
Evaluates the shape function divergence in physical coordinates based on conforming finite element t...
Definition hdiv_fe_transformation.C:403
libMesh::HDivFETransformation::map_phi
virtual void map_phi(const unsigned int dim, const Elem *const elem, const std::vector< Point > &qp, const FEGenericBase< OutputShape > &fe, std::vector< std::vector< OutputShape > > &phi, bool add_p_level=true) const override
Evaluates shape functions in physical coordinates for conforming elements.
Definition hdiv_fe_transformation.C:96
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition int_range.h:153
libMesh::RealGradient
RealVectorValue RealGradient
Definition hp_coarsentest.h:49
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:148
Generated on Thu Nov 12 2026 20:15:22 for libMesh by doxygen 1.9.8