Line data Source code
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/fe_interface.h"
20 : #include "libmesh/int_range.h"
21 :
22 : namespace {
23 :
24 : using namespace libMesh;
25 :
26 : // The contravariant Piola map phi = J^{-1} * (dx/dxi) * phihat, shared by
27 : // map_phi and map_dphi (which also needs the physical shape function value
28 : // for the term coming from the derivative of J^{-1}).
29 : template<typename OutputShape>
30 32384506 : OutputShape hdiv_piola_map(const RealGradient & dxyz_dxi,
31 : const RealGradient & dxyz_deta,
32 : Real J,
33 : const OutputShape & phi_ref)
34 : {
35 38081978 : return (dxyz_dxi*phi_ref(0) + dxyz_deta*phi_ref(1))/J;
36 : }
37 :
38 : template<typename OutputShape>
39 104280598 : OutputShape hdiv_piola_map(const RealGradient & dxyz_dxi,
40 : const RealGradient & dxyz_deta,
41 : const RealGradient & dxyz_dzeta,
42 : Real J,
43 : const OutputShape & phi_ref)
44 : {
45 124102426 : return (dxyz_dxi*phi_ref(0) + dxyz_deta*phi_ref(1) + dxyz_dzeta*phi_ref(2))/J;
46 : }
47 :
48 : } // anonymous namespace
49 :
50 : namespace libMesh
51 : {
52 :
53 : template<typename OutputShape>
54 4455198 : void HDivFETransformation<OutputShape>::init_map_phi(const FEGenericBase<OutputShape> & fe) const
55 : {
56 : // We only need to pre-request one first-order derivative map piece because there is only a
57 : // single calculate_dxyz boolean flag that toggles first derivative computations in the FEMap
58 : // and in addition to covering all components it also covers both forward and inverse mapping.
59 : // We choose to document pre-requesting the forward (reference -> physical) map here because
60 : // that is what is used in map_phi (there is no inverse map usage).
61 319974 : fe.get_fe_map().get_dxyzdxi();
62 4455198 : }
63 :
64 :
65 :
66 : template<typename OutputShape>
67 1695140 : void HDivFETransformation<OutputShape>::init_map_dphi(const FEGenericBase<OutputShape> & fe) const
68 : {
69 : // See above comment in init_map_phi. In map_dphi we actually use both the forward and inverse
70 : // first derivative maps so the choice here is a little more arbitrary. We choose to be
71 : // consistent with init_map_phi.
72 111482 : fe.get_fe_map().get_dxyzdxi();
73 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
74 : // As for first order derivatives, there is only a single boolean flag (calculate_d2xyz)
75 : // controlling second order derivative computations in FEMap so we only bother prerequesting
76 : // one piece. We document pre-requesting the forward map because those are the second
77 : // derivatives used in map_dphi.
78 111482 : fe.get_fe_map().get_d2xyzdxi2();
79 : #endif
80 1695140 : }
81 :
82 :
83 :
84 : template<typename OutputShape>
85 0 : void HDivFETransformation<OutputShape>::init_map_d2phi(const FEGenericBase<OutputShape> & /*fe*/) const
86 : {
87 : // We choose not to pre-request any computations here since we have not yet implemented
88 : // map_d2phi. We don't need to give the map unnecessary work.
89 0 : }
90 :
91 :
92 :
93 : template<typename OutputShape>
94 1726781 : void HDivFETransformation<OutputShape>::map_phi(const unsigned int dim,
95 : const Elem * const elem,
96 : const std::vector<Point> & qp,
97 : const FEGenericBase<OutputShape> & fe,
98 : std::vector<std::vector<OutputShape>> & phi,
99 : const bool /*add_p_level*/) const
100 : {
101 1726781 : switch (dim)
102 : {
103 0 : case 0:
104 : case 1:
105 0 : libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
106 :
107 498352 : case 2:
108 : {
109 38662 : const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
110 38662 : const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
111 :
112 38662 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
113 :
114 : // phi = J^{-1} * (dx/dxi) * \hat{phi}
115 3425468 : for (auto i : index_range(phi))
116 27874510 : for (auto p : index_range(phi[i]))
117 : {
118 : // Need to temporarily cache reference shape functions
119 : // We are computing mapping basis functions, so we explicitly ignore
120 : // any non-zero p_level() the Elem might have.
121 2016562 : OutputShape phi_ref;
122 28980518 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
123 :
124 35030204 : phi[i][p] = hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], J[p], phi_ref);
125 : }
126 :
127 38662 : break;
128 : }
129 1228429 : case 3:
130 : {
131 81444 : const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
132 81444 : const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
133 81444 : const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
134 :
135 81444 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
136 :
137 : // phi = J^{-1} * (dx/dxi) * \hat{phi}
138 6917339 : for (auto i : index_range(phi))
139 81869520 : for (auto p : index_range(phi[i]))
140 : {
141 : // Need to temporarily cache reference shape functions
142 : // We are computing mapping basis functions, so we explicitly ignore
143 : // any non-zero p_level() the Elem might have.
144 4879292 : OutputShape phi_ref;
145 85939194 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
146 :
147 95697778 : phi[i][p] = hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], dxyz_dzeta[p], J[p], phi_ref);
148 : }
149 :
150 81444 : break;
151 : }
152 :
153 0 : default:
154 0 : libmesh_error_msg("Invalid dim = " << dim);
155 : } // switch(dim)
156 1726781 : }
157 :
158 : template<typename OutputShape>
159 402024 : void HDivFETransformation<OutputShape>::map_dphi(const unsigned int dim,
160 : const Elem * const elem,
161 : const std::vector<Point> & qp,
162 : const FEGenericBase<OutputShape> & fe,
163 : std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi,
164 : std::vector<std::vector<OutputShape>> & dphidx,
165 : std::vector<std::vector<OutputShape>> & dphidy,
166 : std::vector<std::vector<OutputShape>> & dphidz) const
167 : {
168 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
169 402024 : switch (dim)
170 : {
171 0 : case 0:
172 : case 1:
173 0 : libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
174 :
175 99354 : case 2:
176 : {
177 7374 : const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
178 7374 : const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
179 :
180 7374 : const std::vector<RealGradient> & d2xyz_dxi2 = fe.get_fe_map().get_d2xyzdxi2();
181 7374 : const std::vector<RealGradient> & d2xyz_deta2 = fe.get_fe_map().get_d2xyzdeta2();
182 7374 : const std::vector<RealGradient> & d2xyz_dxideta = fe.get_fe_map().get_d2xyzdxideta();
183 :
184 7374 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
185 :
186 7374 : const std::vector<Real> & dxi_dx = fe.get_fe_map().get_dxidx();
187 7374 : const std::vector<Real> & dxi_dy = fe.get_fe_map().get_dxidy();
188 7374 : const std::vector<Real> & deta_dx = fe.get_fe_map().get_detadx();
189 7374 : const std::vector<Real> & deta_dy = fe.get_fe_map().get_detady();
190 : #if LIBMESH_DIM > 2
191 7374 : const std::vector<Real> & dxi_dz = fe.get_fe_map().get_dxidz();
192 7374 : const std::vector<Real> & deta_dz = fe.get_fe_map().get_detadz();
193 : #endif
194 :
195 7374 : const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
196 7374 : const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
197 :
198 799978 : for (auto i : index_range(dphi))
199 10986472 : for (auto p : index_range(dphi[i]))
200 : {
201 : // Need to temporarily cache reference shape functions
202 : // We are computing mapping basis functions, so we explicitly ignore
203 : // any non-zero p_level() the Elem might have.
204 832174 : OutputShape phi_ref;
205 11950196 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
206 :
207 : const OutputShape phi_val =
208 11950196 : hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], J[p], phi_ref);
209 :
210 : // Jacobi's formula: dJ/dxi_n = J * tr(F^{-1} * dF/dxi_n)
211 10285848 : Real dJ_dxi =
212 12782370 : J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) +
213 11118022 : dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1));
214 10285848 : Real dJ_deta =
215 10285848 : J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) +
216 11118022 : dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1));
217 : #if LIBMESH_DIM > 2
218 11118022 : dJ_dxi += J[p] * (dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2));
219 10285848 : dJ_deta += J[p] * (dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2));
220 : #endif
221 :
222 41143392 : for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
223 : {
224 : // dphi_k/dx_l = A + B + C, where (n, m summed over reference directions):
225 : // A = -phi_k(x) * SUM_n (dxi_n/dx_l) * (1/J) * (dJ/dxi_n)
226 : // B = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m
227 : // C = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n
228 :
229 : // Term A: -phi_k(x) * SUM_n (dxi_n/dx_l)*(1/J)*(dJ/dxi_n)
230 30857544 : const Real phik = phi_val(k);
231 :
232 : // Term B: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m
233 30857544 : const Real d2xk_dxi2 = d2xyz_dxi2[p](k);
234 30857544 : const Real d2xk_deta2 = d2xyz_deta2[p](k);
235 30857544 : const Real d2xk_dxideta = d2xyz_dxideta[p](k);
236 :
237 30857544 : const Real dphik_dxi_B = (d2xk_dxi2*phi_ref(0) + d2xk_dxideta*phi_ref(1))/J[p];
238 30857544 : const Real dphik_deta_B = (d2xk_dxideta*phi_ref(0) + d2xk_deta2*phi_ref(1))/J[p];
239 :
240 : // Term C: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n
241 30857544 : const Real Fk_xi = dxyz_dxi[p](k);
242 30857544 : const Real Fk_eta = dxyz_deta[p](k);
243 :
244 33354066 : const Real dphik_dxi_C = (Fk_xi*dphi_dxi[i][p](0) + Fk_eta*dphi_dxi[i][p](1))/J[p];
245 33354066 : const Real dphik_deta_C = (Fk_xi*dphi_deta[i][p](0) + Fk_eta*dphi_deta[i][p](1))/J[p];
246 :
247 30857544 : const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C;
248 30857544 : const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C;
249 :
250 33354066 : dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total;
251 33354066 : dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total;
252 : #if LIBMESH_DIM > 2
253 35850588 : dphidz[i][p](k) = dxi_dz[p]*dphik_dxi_total + deta_dz[p]*dphik_deta_total;
254 : #endif
255 : }
256 :
257 13614544 : dphi[i][p].slice(0) = dphidx[i][p];
258 11118022 : dphi[i][p].slice(1) = dphidy[i][p];
259 : #if LIBMESH_DIM > 2
260 11118022 : dphi[i][p].slice(2) = dphidz[i][p];
261 : #endif
262 : }
263 :
264 7374 : break;
265 : }
266 302670 : case 3:
267 : {
268 18670 : const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
269 18670 : const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
270 18670 : const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
271 :
272 18670 : const std::vector<RealGradient> & d2xyz_dxi2 = fe.get_fe_map().get_d2xyzdxi2();
273 18670 : const std::vector<RealGradient> & d2xyz_deta2 = fe.get_fe_map().get_d2xyzdeta2();
274 18670 : const std::vector<RealGradient> & d2xyz_dzeta2 = fe.get_fe_map().get_d2xyzdzeta2();
275 18670 : const std::vector<RealGradient> & d2xyz_dxideta = fe.get_fe_map().get_d2xyzdxideta();
276 18670 : const std::vector<RealGradient> & d2xyz_dxidzeta = fe.get_fe_map().get_d2xyzdxidzeta();
277 18670 : const std::vector<RealGradient> & d2xyz_detadzeta = fe.get_fe_map().get_d2xyzdetadzeta();
278 :
279 18670 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
280 :
281 18670 : const std::vector<Real> & dxi_dx = fe.get_fe_map().get_dxidx();
282 18670 : const std::vector<Real> & dxi_dy = fe.get_fe_map().get_dxidy();
283 18670 : const std::vector<Real> & dxi_dz = fe.get_fe_map().get_dxidz();
284 18670 : const std::vector<Real> & deta_dx = fe.get_fe_map().get_detadx();
285 18670 : const std::vector<Real> & deta_dy = fe.get_fe_map().get_detady();
286 18670 : const std::vector<Real> & deta_dz = fe.get_fe_map().get_detadz();
287 18670 : const std::vector<Real> & dzeta_dx = fe.get_fe_map().get_dzetadx();
288 18670 : const std::vector<Real> & dzeta_dy = fe.get_fe_map().get_dzetady();
289 18670 : const std::vector<Real> & dzeta_dz = fe.get_fe_map().get_dzetadz();
290 :
291 18670 : const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
292 18670 : const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
293 18670 : const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
294 :
295 1731274 : for (auto i : index_range(dphi))
296 29528592 : for (auto p : index_range(dphi[i]))
297 : {
298 : // Need to temporarily cache reference shape functions
299 : // We are computing mapping basis functions, so we explicitly ignore
300 : // any non-zero p_level() the Elem might have.
301 1727984 : OutputShape phi_ref;
302 31555956 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
303 :
304 : const OutputShape phi_val =
305 35011924 : hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], dxyz_dzeta[p], J[p], phi_ref);
306 :
307 : // Jacobi's formula: dJ/dxi_n = J * tr(F^{-1} * dF/dxi_n)
308 28099988 : const Real dJ_dxi =
309 38467892 : J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) + dzeta_dx[p]*d2xyz_dxidzeta[p](0) +
310 31555956 : dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1) + dzeta_dy[p]*d2xyz_dxidzeta[p](1) +
311 31555956 : dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2) + dzeta_dz[p]*d2xyz_dxidzeta[p](2));
312 28099988 : const Real dJ_deta =
313 29827972 : J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) + dzeta_dx[p]*d2xyz_detadzeta[p](0) +
314 31555956 : dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1) + dzeta_dy[p]*d2xyz_detadzeta[p](1) +
315 29827972 : dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2) + dzeta_dz[p]*d2xyz_detadzeta[p](2));
316 28099988 : const Real dJ_dzeta =
317 28099988 : J[p] * (dxi_dx[p]*d2xyz_dxidzeta[p](0) + deta_dx[p]*d2xyz_detadzeta[p](0) + dzeta_dx[p]*d2xyz_dzeta2[p](0) +
318 31555956 : dxi_dy[p]*d2xyz_dxidzeta[p](1) + deta_dy[p]*d2xyz_detadzeta[p](1) + dzeta_dy[p]*d2xyz_dzeta2[p](1) +
319 28099988 : dxi_dz[p]*d2xyz_dxidzeta[p](2) + deta_dz[p]*d2xyz_detadzeta[p](2) + dzeta_dz[p]*d2xyz_dzeta2[p](2));
320 :
321 112399952 : for (unsigned int k = 0; k < 3; ++k)
322 : {
323 : // dphi_k/dx_l = A + B + C, where (n, m summed over reference directions):
324 : // A = -phi_k(x) * SUM_n (dxi_n/dx_l) * (1/J) * (dJ/dxi_n)
325 : // B = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m
326 : // C = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n
327 :
328 : // Term A: -phi_k(x) * SUM_n (dxi_n/dx_l)*(1/J)*(dJ/dxi_n)
329 84299964 : const Real phik = phi_val(k);
330 :
331 : // Term B: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m
332 84299964 : const Real d2xk_dxi2 = d2xyz_dxi2[p](k);
333 84299964 : const Real d2xk_deta2 = d2xyz_deta2[p](k);
334 84299964 : const Real d2xk_dzeta2 = d2xyz_dzeta2[p](k);
335 84299964 : const Real d2xk_dxideta = d2xyz_dxideta[p](k);
336 84299964 : const Real d2xk_dxidzeta = d2xyz_dxidzeta[p](k);
337 84299964 : const Real d2xk_detadzeta = d2xyz_detadzeta[p](k);
338 :
339 84299964 : const Real dphik_dxi_B = (d2xk_dxi2*phi_ref(0) + d2xk_dxideta*phi_ref(1) + d2xk_dxidzeta*phi_ref(2))/J[p];
340 84299964 : const Real dphik_deta_B = (d2xk_dxideta*phi_ref(0) + d2xk_deta2*phi_ref(1) + d2xk_detadzeta*phi_ref(2))/J[p];
341 84299964 : const Real dphik_dzeta_B = (d2xk_dxidzeta*phi_ref(0) + d2xk_detadzeta*phi_ref(1) + d2xk_dzeta2*phi_ref(2))/J[p];
342 :
343 : // Term C: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n
344 84299964 : const Real Fk_xi = dxyz_dxi[p](k);
345 84299964 : const Real Fk_eta = dxyz_deta[p](k);
346 84299964 : const Real Fk_zeta = dxyz_dzeta[p](k);
347 :
348 89483916 : 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];
349 89483916 : 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];
350 89483916 : 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];
351 :
352 84299964 : const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C;
353 84299964 : const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C;
354 84299964 : const Real dphik_dzeta_total = -phik*dJ_dzeta/J[p] + dphik_dzeta_B + dphik_dzeta_C;
355 :
356 89483916 : dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total + dzeta_dx[p]*dphik_dzeta_total;
357 89483916 : dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total + dzeta_dy[p]*dphik_dzeta_total;
358 94667868 : dphidz[i][p](k) = dxi_dz[p]*dphik_dxi_total + deta_dz[p]*dphik_deta_total + dzeta_dz[p]*dphik_dzeta_total;
359 : }
360 :
361 35011924 : dphi[i][p].slice(0) = dphidx[i][p];
362 29827972 : dphi[i][p].slice(1) = dphidy[i][p];
363 29827972 : dphi[i][p].slice(2) = dphidz[i][p];
364 : }
365 :
366 18670 : break;
367 : }
368 :
369 0 : default:
370 0 : libmesh_error_msg("Invalid dim = " << dim);
371 : } // switch(dim)
372 : #else
373 : libmesh_ignore(dim, elem, qp, fe, dphi, dphidx, dphidy, dphidz);
374 : libmesh_error_msg("HDiv shape function gradients require the library to be configured "
375 : "with --enable-second-derivatives (LIBMESH_ENABLE_SECOND_DERIVATIVES).");
376 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
377 402024 : }
378 :
379 : template<typename OutputShape>
380 824650 : void HDivFETransformation<OutputShape>::map_div(const unsigned int dim,
381 : const Elem * const,
382 : const std::vector<Point> &,
383 : const FEGenericBase<OutputShape> & fe,
384 : std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence>> & div_phi) const
385 : {
386 824650 : switch (dim)
387 : {
388 0 : case 0:
389 : case 1:
390 0 : libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
391 :
392 15756 : case 2:
393 : {
394 15756 : const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
395 15756 : const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
396 :
397 15756 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
398 :
399 : // div(phi) = J^{-1} * div(\hat{phi})
400 1652046 : for (auto i : index_range(div_phi))
401 18442044 : for (auto p : index_range(div_phi[i]))
402 : {
403 25271740 : div_phi[i][p] = (dphi_dxi[i][p](0) + dphi_deta[i][p](1))/J[p];
404 : }
405 :
406 15756 : break;
407 : }
408 38492 : case 3:
409 : {
410 38492 : const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
411 38492 : const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
412 38492 : const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
413 :
414 38492 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
415 :
416 : // div(phi) = J^{-1} * div(\hat{phi})
417 3501096 : for (auto i : index_range(div_phi))
418 41817204 : for (auto p : index_range(div_phi[i]))
419 : {
420 58424340 : div_phi[i][p] = (dphi_dxi[i][p](0) + dphi_deta[i][p](1) + dphi_dzeta[i][p](2))/J[p];
421 : }
422 :
423 38492 : break;
424 : }
425 :
426 0 : default:
427 0 : libmesh_error_msg("Invalid dim = " << dim);
428 : } // switch(dim)
429 824650 : }
430 :
431 : template class LIBMESH_EXPORT HDivFETransformation<RealGradient>;
432 :
433 : template<>
434 0 : void HDivFETransformation<Real>::init_map_phi(const FEGenericBase<Real> & ) const
435 : {
436 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
437 : }
438 :
439 : template<>
440 0 : void HDivFETransformation<Real>::init_map_dphi(const FEGenericBase<Real> & ) const
441 : {
442 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
443 : }
444 :
445 : template<>
446 0 : void HDivFETransformation<Real>::init_map_d2phi(const FEGenericBase<Real> & ) const
447 : {
448 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
449 : }
450 :
451 : template<>
452 0 : void HDivFETransformation<Real>::map_phi(const unsigned int,
453 : const Elem * const,
454 : const std::vector<Point> &,
455 : const FEGenericBase<Real> &,
456 : std::vector<std::vector<Real>> &,
457 : bool) const
458 : {
459 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
460 : }
461 :
462 : template<>
463 0 : void HDivFETransformation<Real>::map_dphi(const unsigned int,
464 : const Elem * const,
465 : const std::vector<Point> &,
466 : const FEGenericBase<Real> &,
467 : std::vector<std::vector<FEGenericBase<Real>::OutputGradient>> &,
468 : std::vector<std::vector<Real>> &,
469 : std::vector<std::vector<Real>> &,
470 : std::vector<std::vector<Real>> &) const
471 : {
472 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
473 : }
474 :
475 : template<>
476 0 : void HDivFETransformation<Real>::map_div(const unsigned int,
477 : const Elem * const,
478 : const std::vector<Point> &,
479 : const FEGenericBase<Real> &,
480 : std::vector<std::vector<Real>> &) const
481 : {
482 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
483 : }
484 :
485 :
486 : } // namespace libMesh
|
