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 8661238 : OutputShape hdiv_piola_map(const RealGradient & dx_dxi,
31 : const RealGradient & dx_deta,
32 : Real J,
33 : const OutputShape & phi_ref)
34 : {
35 2848736 : OutputShape phi_val;
36 11509974 : phi_val(0) = (dx_dxi(0)*phi_ref(0) + dx_deta(0)*phi_ref(1))/J;
37 11509974 : phi_val(1) = (dx_dxi(1)*phi_ref(0) + dx_deta(1)*phi_ref(1))/J;
38 : #if LIBMESH_DIM > 2
39 11509974 : phi_val(2) = (dx_dxi(2)*phi_ref(0) + dx_deta(2)*phi_ref(1))/J;
40 : #endif
41 8661238 : return phi_val;
42 : }
43 :
44 : template<typename OutputShape>
45 29182982 : OutputShape hdiv_piola_map(const RealGradient & dx_dxi,
46 : const RealGradient & dx_deta,
47 : const RealGradient & dx_dzeta,
48 : Real J,
49 : const OutputShape & phi_ref)
50 : {
51 6607276 : OutputShape phi_val;
52 29182982 : phi_val(0) = (dx_dxi(0)*phi_ref(0) + dx_deta(0)*phi_ref(1) + dx_dzeta(0)*phi_ref(2))/J;
53 29182982 : phi_val(1) = (dx_dxi(1)*phi_ref(0) + dx_deta(1)*phi_ref(1) + dx_dzeta(1)*phi_ref(2))/J;
54 29182982 : phi_val(2) = (dx_dxi(2)*phi_ref(0) + dx_deta(2)*phi_ref(1) + dx_dzeta(2)*phi_ref(2))/J;
55 29182982 : return phi_val;
56 : }
57 :
58 : } // anonymous namespace
59 :
60 : namespace libMesh
61 : {
62 :
63 : template<typename OutputShape>
64 1339905 : void HDivFETransformation<OutputShape>::init_map_phi(const FEGenericBase<OutputShape> & fe) const
65 : {
66 319974 : fe.get_fe_map().get_dxidx();
67 1339905 : }
68 :
69 :
70 :
71 : template<typename OutputShape>
72 485047 : void HDivFETransformation<OutputShape>::init_map_dphi(const FEGenericBase<OutputShape> & fe) const
73 : {
74 111482 : fe.get_fe_map().get_dxidx();
75 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
76 111482 : fe.get_fe_map().get_d2xyzdxi2();
77 : #endif
78 485047 : }
79 :
80 :
81 :
82 : template<typename OutputShape>
83 0 : void HDivFETransformation<OutputShape>::init_map_d2phi(const FEGenericBase<OutputShape> & fe) const
84 : {
85 0 : fe.get_fe_map().get_dxidx();
86 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
87 0 : fe.get_fe_map().get_d2xidxyz2();
88 : #endif
89 0 : }
90 :
91 :
92 :
93 : template<typename OutputShape>
94 509996 : 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 509996 : 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 157032 : 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 1101900 : for (auto i : index_range(phi))
116 9086222 : 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 12174478 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
123 :
124 18224164 : phi[i][p] = hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], J[p], phi_ref);
125 : }
126 :
127 38662 : break;
128 : }
129 352964 : 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 1977784 : for (auto i : index_range(phi))
139 23079828 : 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 31213592 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
146 :
147 40972176 : 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 509996 : }
157 :
158 : template<typename OutputShape>
159 113819 : 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 113819 : 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 30770 : 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 256670 : for (auto i : index_range(dphi))
199 3594520 : 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 5032968 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
206 :
207 : const OutputShape phi_val =
208 5032968 : 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 3368620 : Real dJ_dxi =
212 5865142 : J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) +
213 4200794 : dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1));
214 3368620 : Real dJ_deta =
215 3368620 : J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) +
216 4200794 : dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1));
217 : #if LIBMESH_DIM > 2
218 4200794 : dJ_dxi += J[p] * (dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2));
219 3368620 : dJ_deta += J[p] * (dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2));
220 : #endif
221 :
222 13474480 : 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 10105860 : const Real phik = (k == 0) ? phi_val(0) : (k == 1) ? phi_val(1) : phi_val(2);
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 10105860 : const Real d2x_k_dxi2 = (k == 0) ? d2xyz_dxi2[p](0) : (k == 1) ? d2xyz_dxi2[p](1) : d2xyz_dxi2[p](2);
234 10105860 : const Real d2x_k_deta2 = (k == 0) ? d2xyz_deta2[p](0) : (k == 1) ? d2xyz_deta2[p](1) : d2xyz_deta2[p](2);
235 10105860 : const Real d2x_k_dxideta = (k == 0) ? d2xyz_dxideta[p](0) : (k == 1) ? d2xyz_dxideta[p](1) : d2xyz_dxideta[p](2);
236 :
237 10105860 : const Real dphik_dxi_B = (d2x_k_dxi2*phi_ref(0) + d2x_k_dxideta*phi_ref(1))/J[p];
238 10105860 : const Real dphik_deta_B = (d2x_k_dxideta*phi_ref(0) + d2x_k_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 10105860 : const Real Fk_xi = (k == 0) ? dxyz_dxi[p](0) : (k == 1) ? dxyz_dxi[p](1) : dxyz_dxi[p](2);
242 10105860 : const Real Fk_eta = (k == 0) ? dxyz_deta[p](0) : (k == 1) ? dxyz_deta[p](1) : dxyz_deta[p](2);
243 :
244 12602382 : const Real dphik_dxi_C = (Fk_xi*dphi_dxi[i][p](0) + Fk_eta*dphi_dxi[i][p](1))/J[p];
245 12602382 : const Real dphik_deta_C = (Fk_xi*dphi_deta[i][p](0) + Fk_eta*dphi_deta[i][p](1))/J[p];
246 :
247 10105860 : const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C;
248 10105860 : const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C;
249 :
250 12602382 : dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total;
251 12602382 : dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total;
252 : #if LIBMESH_DIM > 2
253 15098904 : dphidz[i][p](k) = dxi_dz[p]*dphik_dxi_total + deta_dz[p]*dphik_deta_total;
254 : #endif
255 : }
256 :
257 6697316 : dphi[i][p].slice(0) = dphidx[i][p];
258 4200794 : dphi[i][p].slice(1) = dphidy[i][p];
259 : #if LIBMESH_DIM > 2
260 4200794 : dphi[i][p].slice(2) = dphidz[i][p];
261 : #endif
262 : }
263 :
264 7374 : break;
265 : }
266 83049 : 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 474875 : for (auto i : index_range(dphi))
296 8119800 : 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 11183942 : FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref);
303 :
304 : const OutputShape phi_val =
305 14639910 : 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 7727974 : const Real dJ_dxi =
309 18095878 : J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) + dzeta_dx[p]*d2xyz_dxidzeta[p](0) +
310 11183942 : dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1) + dzeta_dy[p]*d2xyz_dxidzeta[p](1) +
311 11183942 : dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2) + dzeta_dz[p]*d2xyz_dxidzeta[p](2));
312 7727974 : const Real dJ_deta =
313 9455958 : J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) + dzeta_dx[p]*d2xyz_detadzeta[p](0) +
314 11183942 : dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1) + dzeta_dy[p]*d2xyz_detadzeta[p](1) +
315 9455958 : dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2) + dzeta_dz[p]*d2xyz_detadzeta[p](2));
316 7727974 : const Real dJ_dzeta =
317 7727974 : J[p] * (dxi_dx[p]*d2xyz_dxidzeta[p](0) + deta_dx[p]*d2xyz_detadzeta[p](0) + dzeta_dx[p]*d2xyz_dzeta2[p](0) +
318 11183942 : dxi_dy[p]*d2xyz_dxidzeta[p](1) + deta_dy[p]*d2xyz_detadzeta[p](1) + dzeta_dy[p]*d2xyz_dzeta2[p](1) +
319 7727974 : dxi_dz[p]*d2xyz_dxidzeta[p](2) + deta_dz[p]*d2xyz_detadzeta[p](2) + dzeta_dz[p]*d2xyz_dzeta2[p](2));
320 :
321 30911896 : 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 23183922 : const Real phik = (k == 0) ? phi_val(0) : (k == 1) ? phi_val(1) : phi_val(2);
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 23183922 : const Real d2x_k_dxi2 = (k == 0) ? d2xyz_dxi2[p](0) : (k == 1) ? d2xyz_dxi2[p](1) : d2xyz_dxi2[p](2);
333 23183922 : const Real d2x_k_deta2 = (k == 0) ? d2xyz_deta2[p](0) : (k == 1) ? d2xyz_deta2[p](1) : d2xyz_deta2[p](2);
334 23183922 : const Real d2x_k_dzeta2 = (k == 0) ? d2xyz_dzeta2[p](0) : (k == 1) ? d2xyz_dzeta2[p](1) : d2xyz_dzeta2[p](2);
335 23183922 : const Real d2x_k_dxideta = (k == 0) ? d2xyz_dxideta[p](0) : (k == 1) ? d2xyz_dxideta[p](1) : d2xyz_dxideta[p](2);
336 23183922 : const Real d2x_k_dxidzeta = (k == 0) ? d2xyz_dxidzeta[p](0) : (k == 1) ? d2xyz_dxidzeta[p](1) : d2xyz_dxidzeta[p](2);
337 23183922 : const Real d2x_k_detadzeta = (k == 0) ? d2xyz_detadzeta[p](0) : (k == 1) ? d2xyz_detadzeta[p](1) : d2xyz_detadzeta[p](2);
338 :
339 23183922 : const Real dphik_dxi_B = (d2x_k_dxi2*phi_ref(0) + d2x_k_dxideta*phi_ref(1) + d2x_k_dxidzeta*phi_ref(2))/J[p];
340 23183922 : const Real dphik_deta_B = (d2x_k_dxideta*phi_ref(0) + d2x_k_deta2*phi_ref(1) + d2x_k_detadzeta*phi_ref(2))/J[p];
341 23183922 : const Real dphik_dzeta_B = (d2x_k_dxidzeta*phi_ref(0) + d2x_k_detadzeta*phi_ref(1) + d2x_k_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 23183922 : const Real Fk_xi = (k == 0) ? dxyz_dxi[p](0) : (k == 1) ? dxyz_dxi[p](1) : dxyz_dxi[p](2);
345 23183922 : const Real Fk_eta = (k == 0) ? dxyz_deta[p](0) : (k == 1) ? dxyz_deta[p](1) : dxyz_deta[p](2);
346 23183922 : const Real Fk_zeta = (k == 0) ? dxyz_dzeta[p](0) : (k == 1) ? dxyz_dzeta[p](1) : dxyz_dzeta[p](2);
347 :
348 28367874 : 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 28367874 : 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 28367874 : 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 23183922 : const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C;
353 23183922 : const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C;
354 23183922 : const Real dphik_dzeta_total = -phik*dJ_dzeta/J[p] + dphik_dzeta_B + dphik_dzeta_C;
355 :
356 28367874 : dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total + dzeta_dx[p]*dphik_dzeta_total;
357 28367874 : dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total + dzeta_dy[p]*dphik_dzeta_total;
358 33551826 : 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 14639910 : dphi[i][p].slice(0) = dphidx[i][p];
362 9455958 : dphi[i][p].slice(1) = dphidy[i][p];
363 9455958 : 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 113819 : }
378 :
379 : template<typename OutputShape>
380 236485 : 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 236485 : 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 530838 : for (auto i : index_range(div_phi))
401 6041358 : for (auto p : index_range(div_phi[i]))
402 : {
403 13846972 : 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 972315 : for (auto i : index_range(div_phi))
418 11624328 : for (auto p : index_range(div_phi[i]))
419 : {
420 30317370 : 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 236485 : }
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
|
