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/elem.h"
20 : #include "libmesh/fe_interface.h"
21 : #include "libmesh/int_range.h"
22 : #include "libmesh/tensor_value.h"
23 :
24 : namespace {
25 :
26 : using 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}).
31 : template<typename OutputShape>
32 32384506 : OutputShape hdiv_piola_map(const RealGradient & dxyz_dxi,
33 : const RealGradient & dxyz_deta,
34 : Real J,
35 : const OutputShape & phi_ref)
36 : {
37 38081978 : return (dxyz_dxi*phi_ref(0) + dxyz_deta*phi_ref(1))/J;
38 : }
39 :
40 : template<typename OutputShape>
41 104280598 : OutputShape 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 124102426 : return (dxyz_dxi*phi_ref(0) + dxyz_deta*phi_ref(1) + dxyz_dzeta*phi_ref(2))/J;
48 : }
49 :
50 : } // anonymous namespace
51 :
52 : namespace libMesh
53 : {
54 :
55 : template<typename OutputShape>
56 4455198 : void 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 319974 : fe.get_fe_map().get_dxyzdxi();
64 4455198 : }
65 :
66 :
67 :
68 : template<typename OutputShape>
69 1695140 : void 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 111482 : 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 111482 : fe.get_fe_map().get_d2xyzdxi2();
81 : #endif
82 1695140 : }
83 :
84 :
85 :
86 : template<typename OutputShape>
87 0 : void 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 0 : }
92 :
93 :
94 :
95 : template<typename OutputShape>
96 1726781 : void 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 1726781 : switch (dim)
104 : {
105 0 : case 0:
106 : case 1:
107 0 : libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
108 :
109 498352 : case 2:
110 : {
111 38662 : const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
112 38662 : const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
113 :
114 38662 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
115 :
116 : // phi = J^{-1} * (dx/dxi) * \hat{phi}
117 3425468 : for (auto i : index_range(phi))
118 27874510 : 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 2016562 : OutputShape phi_ref;
124 28980518 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
125 :
126 35030204 : phi[i][p] = hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], J[p], phi_ref);
127 : }
128 :
129 38662 : break;
130 : }
131 1228429 : case 3:
132 : {
133 81444 : const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
134 81444 : const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
135 81444 : const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
136 :
137 81444 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
138 :
139 : // phi = J^{-1} * (dx/dxi) * \hat{phi}
140 6917339 : for (auto i : index_range(phi))
141 81869520 : 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 4879292 : OutputShape phi_ref;
147 85939194 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
148 :
149 95697778 : phi[i][p] = hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], dxyz_dzeta[p], J[p], phi_ref);
150 : }
151 :
152 81444 : break;
153 : }
154 :
155 0 : default:
156 0 : libmesh_error_msg("Invalid dim = " << dim);
157 : } // switch(dim)
158 1726781 : }
159 :
160 : template<typename OutputShape>
161 402024 : void 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 402024 : switch (dim)
171 : {
172 0 : case 0:
173 : case 1:
174 0 : libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
175 :
176 99354 : case 2:
177 : {
178 7374 : const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
179 7374 : const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
180 :
181 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
182 7374 : const std::vector<RealGradient> & d2xyz_dxi2 = fe.get_fe_map().get_d2xyzdxi2();
183 7374 : const std::vector<RealGradient> & d2xyz_deta2 = fe.get_fe_map().get_d2xyzdeta2();
184 7374 : 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 7374 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
198 :
199 7374 : const std::vector<Real> & dxi_dx = fe.get_fe_map().get_dxidx();
200 7374 : const std::vector<Real> & dxi_dy = fe.get_fe_map().get_dxidy();
201 7374 : const std::vector<Real> & deta_dx = fe.get_fe_map().get_detadx();
202 7374 : const std::vector<Real> & deta_dy = fe.get_fe_map().get_detady();
203 : #if LIBMESH_DIM > 2
204 7374 : const std::vector<Real> & dxi_dz = fe.get_fe_map().get_dxidz();
205 7374 : const std::vector<Real> & deta_dz = fe.get_fe_map().get_detadz();
206 : #endif
207 :
208 7374 : const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
209 7374 : const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
210 :
211 799978 : for (auto i : index_range(dphi))
212 10986472 : 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 832174 : OutputShape phi_ref;
218 11950196 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
219 :
220 : const OutputShape phi_val =
221 11950196 : 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 10285848 : Real dJ_dxi =
225 12782370 : J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) +
226 11118022 : dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1));
227 10285848 : Real dJ_deta =
228 10285848 : J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) +
229 11118022 : dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1));
230 : #if LIBMESH_DIM > 2
231 11118022 : dJ_dxi += J[p] * (dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2));
232 10285848 : dJ_deta += J[p] * (dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2));
233 : #endif
234 :
235 41143392 : 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 30857544 : 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 30857544 : const Real d2xk_dxi2 = d2xyz_dxi2[p](k);
247 30857544 : const Real d2xk_deta2 = d2xyz_deta2[p](k);
248 30857544 : const Real d2xk_dxideta = d2xyz_dxideta[p](k);
249 :
250 30857544 : const Real dphik_dxi_B = (d2xk_dxi2*phi_ref(0) + d2xk_dxideta*phi_ref(1))/J[p];
251 30857544 : 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 30857544 : const Real Fk_xi = dxyz_dxi[p](k);
255 30857544 : const Real Fk_eta = dxyz_deta[p](k);
256 :
257 33354066 : const Real dphik_dxi_C = (Fk_xi*dphi_dxi[i][p](0) + Fk_eta*dphi_dxi[i][p](1))/J[p];
258 33354066 : const Real dphik_deta_C = (Fk_xi*dphi_deta[i][p](0) + Fk_eta*dphi_deta[i][p](1))/J[p];
259 :
260 30857544 : const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C;
261 30857544 : const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C;
262 :
263 33354066 : dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total;
264 33354066 : dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total;
265 : #if LIBMESH_DIM > 2
266 35850588 : dphidz[i][p](k) = dxi_dz[p]*dphik_dxi_total + deta_dz[p]*dphik_deta_total;
267 : #endif
268 : }
269 :
270 13614544 : dphi[i][p].slice(0) = dphidx[i][p];
271 11118022 : dphi[i][p].slice(1) = dphidy[i][p];
272 : #if LIBMESH_DIM > 2
273 11118022 : dphi[i][p].slice(2) = dphidz[i][p];
274 : #endif
275 : }
276 :
277 7374 : break;
278 : }
279 302670 : case 3:
280 : {
281 18670 : const std::vector<RealGradient> & dxyz_dxi = fe.get_fe_map().get_dxyzdxi();
282 18670 : const std::vector<RealGradient> & dxyz_deta = fe.get_fe_map().get_dxyzdeta();
283 18670 : const std::vector<RealGradient> & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta();
284 :
285 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
286 18670 : const std::vector<RealGradient> & d2xyz_dxi2 = fe.get_fe_map().get_d2xyzdxi2();
287 18670 : const std::vector<RealGradient> & d2xyz_deta2 = fe.get_fe_map().get_d2xyzdeta2();
288 18670 : const std::vector<RealGradient> & d2xyz_dzeta2 = fe.get_fe_map().get_d2xyzdzeta2();
289 18670 : const std::vector<RealGradient> & d2xyz_dxideta = fe.get_fe_map().get_d2xyzdxideta();
290 18670 : const std::vector<RealGradient> & d2xyz_dxidzeta = fe.get_fe_map().get_d2xyzdxidzeta();
291 18670 : 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 18670 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
308 :
309 18670 : const std::vector<Real> & dxi_dx = fe.get_fe_map().get_dxidx();
310 18670 : const std::vector<Real> & dxi_dy = fe.get_fe_map().get_dxidy();
311 18670 : const std::vector<Real> & dxi_dz = fe.get_fe_map().get_dxidz();
312 18670 : const std::vector<Real> & deta_dx = fe.get_fe_map().get_detadx();
313 18670 : const std::vector<Real> & deta_dy = fe.get_fe_map().get_detady();
314 18670 : const std::vector<Real> & deta_dz = fe.get_fe_map().get_detadz();
315 18670 : const std::vector<Real> & dzeta_dx = fe.get_fe_map().get_dzetadx();
316 18670 : const std::vector<Real> & dzeta_dy = fe.get_fe_map().get_dzetady();
317 18670 : const std::vector<Real> & dzeta_dz = fe.get_fe_map().get_dzetadz();
318 :
319 18670 : const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
320 18670 : const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
321 18670 : const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
322 :
323 1731274 : for (auto i : index_range(dphi))
324 29528592 : 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 1727984 : OutputShape phi_ref;
330 31555956 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
331 :
332 : const OutputShape phi_val =
333 35011924 : 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 28099988 : const Real dJ_dxi =
337 38467892 : J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) + dzeta_dx[p]*d2xyz_dxidzeta[p](0) +
338 31555956 : dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1) + dzeta_dy[p]*d2xyz_dxidzeta[p](1) +
339 31555956 : dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2) + dzeta_dz[p]*d2xyz_dxidzeta[p](2));
340 28099988 : const Real dJ_deta =
341 29827972 : J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) + dzeta_dx[p]*d2xyz_detadzeta[p](0) +
342 31555956 : dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1) + dzeta_dy[p]*d2xyz_detadzeta[p](1) +
343 29827972 : dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2) + dzeta_dz[p]*d2xyz_detadzeta[p](2));
344 28099988 : const Real dJ_dzeta =
345 28099988 : J[p] * (dxi_dx[p]*d2xyz_dxidzeta[p](0) + deta_dx[p]*d2xyz_detadzeta[p](0) + dzeta_dx[p]*d2xyz_dzeta2[p](0) +
346 31555956 : dxi_dy[p]*d2xyz_dxidzeta[p](1) + deta_dy[p]*d2xyz_detadzeta[p](1) + dzeta_dy[p]*d2xyz_dzeta2[p](1) +
347 28099988 : dxi_dz[p]*d2xyz_dxidzeta[p](2) + deta_dz[p]*d2xyz_detadzeta[p](2) + dzeta_dz[p]*d2xyz_dzeta2[p](2));
348 :
349 112399952 : 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 84299964 : 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 84299964 : const Real d2xk_dxi2 = d2xyz_dxi2[p](k);
361 84299964 : const Real d2xk_deta2 = d2xyz_deta2[p](k);
362 84299964 : const Real d2xk_dzeta2 = d2xyz_dzeta2[p](k);
363 84299964 : const Real d2xk_dxideta = d2xyz_dxideta[p](k);
364 84299964 : const Real d2xk_dxidzeta = d2xyz_dxidzeta[p](k);
365 84299964 : const Real d2xk_detadzeta = d2xyz_detadzeta[p](k);
366 :
367 84299964 : const Real dphik_dxi_B = (d2xk_dxi2*phi_ref(0) + d2xk_dxideta*phi_ref(1) + d2xk_dxidzeta*phi_ref(2))/J[p];
368 84299964 : const Real dphik_deta_B = (d2xk_dxideta*phi_ref(0) + d2xk_deta2*phi_ref(1) + d2xk_detadzeta*phi_ref(2))/J[p];
369 84299964 : 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 84299964 : const Real Fk_xi = dxyz_dxi[p](k);
373 84299964 : const Real Fk_eta = dxyz_deta[p](k);
374 84299964 : const Real Fk_zeta = dxyz_dzeta[p](k);
375 :
376 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];
377 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];
378 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];
379 :
380 84299964 : const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C;
381 84299964 : const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C;
382 84299964 : const Real dphik_dzeta_total = -phik*dJ_dzeta/J[p] + dphik_dzeta_B + dphik_dzeta_C;
383 :
384 89483916 : dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total + dzeta_dx[p]*dphik_dzeta_total;
385 89483916 : dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total + dzeta_dy[p]*dphik_dzeta_total;
386 94667868 : 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 35011924 : dphi[i][p].slice(0) = dphidx[i][p];
390 29827972 : dphi[i][p].slice(1) = dphidy[i][p];
391 29827972 : dphi[i][p].slice(2) = dphidz[i][p];
392 : }
393 :
394 18670 : break;
395 : }
396 :
397 0 : default:
398 0 : libmesh_error_msg("Invalid dim = " << dim);
399 : } // switch(dim)
400 402024 : }
401 :
402 : template<typename OutputShape>
403 824650 : void 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 824650 : switch (dim)
410 : {
411 0 : case 0:
412 : case 1:
413 0 : libmesh_error_msg("These element transformations only make sense in 2D and 3D.");
414 :
415 15756 : case 2:
416 : {
417 15756 : const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
418 15756 : const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
419 :
420 15756 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
421 :
422 : // div(phi) = J^{-1} * div(\hat{phi})
423 1652046 : for (auto i : index_range(div_phi))
424 18442044 : for (auto p : index_range(div_phi[i]))
425 : {
426 25271740 : div_phi[i][p] = (dphi_dxi[i][p](0) + dphi_deta[i][p](1))/J[p];
427 : }
428 :
429 15756 : break;
430 : }
431 38492 : case 3:
432 : {
433 38492 : const std::vector<std::vector<OutputShape>> & dphi_dxi = fe.get_dphidxi();
434 38492 : const std::vector<std::vector<OutputShape>> & dphi_deta = fe.get_dphideta();
435 38492 : const std::vector<std::vector<OutputShape>> & dphi_dzeta = fe.get_dphidzeta();
436 :
437 38492 : const std::vector<Real> & J = fe.get_fe_map().get_jacobian();
438 :
439 : // div(phi) = J^{-1} * div(\hat{phi})
440 3501096 : for (auto i : index_range(div_phi))
441 41817204 : for (auto p : index_range(div_phi[i]))
442 : {
443 58424340 : div_phi[i][p] = (dphi_dxi[i][p](0) + dphi_deta[i][p](1) + dphi_dzeta[i][p](2))/J[p];
444 : }
445 :
446 38492 : break;
447 : }
448 :
449 0 : default:
450 0 : libmesh_error_msg("Invalid dim = " << dim);
451 : } // switch(dim)
452 824650 : }
453 :
454 : template class LIBMESH_EXPORT HDivFETransformation<RealGradient>;
455 :
456 : template<>
457 0 : void HDivFETransformation<Real>::init_map_phi(const FEGenericBase<Real> & ) 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>::init_map_dphi(const FEGenericBase<Real> & ) const
464 : {
465 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
466 : }
467 :
468 : template<>
469 0 : void HDivFETransformation<Real>::init_map_d2phi(const FEGenericBase<Real> & ) const
470 : {
471 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
472 : }
473 :
474 : template<>
475 0 : void 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 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
483 : }
484 :
485 : template<>
486 0 : void 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 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
496 : }
497 :
498 : template<>
499 0 : void 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 0 : libmesh_error_msg("HDiv transformations only make sense for vector-valued elements.");
506 : }
507 :
508 :
509 : } // namespace libMesh
|
