libMesh
Loading...
Searching...
No Matches
h1_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/fe_interface.h"
19#include "libmesh/h1_fe_transformation.h"
20#include "libmesh/tensor_value.h"
21#include "libmesh/elem.h"
22#include "libmesh/int_range.h"
23
24namespace libMesh
25{
26template<typename OutputShape>
27void H1FETransformation<OutputShape>::init_map_phi(const FEGenericBase<OutputShape> & /* fe */ ) const
28{
29 // Nothing needed
30}
31
32
33
34template<typename OutputShape>
39
40
41
42template<typename OutputShape>
43void H1FETransformation<OutputShape>::init_map_d2phi(const FEGenericBase<OutputShape> & fe) const
44{
45 fe.get_fe_map().get_dxidx();
46#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
47 fe.get_fe_map().get_d2xidxyz2();
48#endif
49}
50
51template<typename OutputShape>
52void H1FETransformation<OutputShape>::map_phi( const unsigned int dim,
53 const Elem * const elem,
54 const std::vector<Point> & qp,
55 const FEGenericBase<OutputShape> & fe,
56 std::vector<std::vector<OutputShape>> & phi,
57 const bool add_p_level) const
58{
59 FEInterface::all_shapes<OutputShape>(dim, fe.get_fe_type(), elem, qp, phi, add_p_level);
60}
61
62
63template<typename OutputShape>
64void H1FETransformation<OutputShape>::map_dphi( const unsigned int dim,
65 const Elem * const,
66 const std::vector<Point> &,
67 const FEGenericBase<OutputShape> & fe,
68 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> & dphi,
69 std::vector<std::vector<OutputShape>> & dphidx,
70 std::vector<std::vector<OutputShape>> & dphidy,
71 std::vector<std::vector<OutputShape>> & dphidz ) const
72{
73 switch(dim)
74 {
75 case 0: // No derivatives in 0D
76 {
77 for (auto i : index_range(dphi))
78 for (auto p : index_range(dphi[i]))
79 dphi[i][p] = 0.;
80 break;
81 }
82
83 case 1:
84 {
85 const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
86
87 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
88#if LIBMESH_DIM>1
89 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
90#endif
91#if LIBMESH_DIM>2
92 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
93#endif
94
95 for (auto i : index_range(dphi))
96 for (auto p : index_range(dphi[i]))
97 {
98 // dphi/dx = (dphi/dxi)*(dxi/dx)
99 dphi[i][p].slice(0) = dphidx[i][p] = dphidxi[i][p]*dxidx_map[p];
100
101#if LIBMESH_DIM>1
102 dphi[i][p].slice(1) = dphidy[i][p] = dphidxi[i][p]*dxidy_map[p];
103#endif
104#if LIBMESH_DIM>2
105 dphi[i][p].slice(2) = dphidz[i][p] = dphidxi[i][p]*dxidz_map[p];
106#endif
107 }
108
109 break;
110 }
111
112 case 2:
113 {
114 const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
115 const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
116
117 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
118 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
119#if LIBMESH_DIM > 2
120 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
121#endif
122
123 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
124 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
125#if LIBMESH_DIM > 2
126 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
127#endif
128
129 for (auto i : index_range(dphi))
130 for (auto p : index_range(dphi[i]))
131 {
132 // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx)
133 dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
134 dphideta[i][p]*detadx_map[p]);
135
136 // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy)
137 dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
138 dphideta[i][p]*detady_map[p]);
139
140#if LIBMESH_DIM > 2
141 // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz)
142 dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
143 dphideta[i][p]*detadz_map[p]);
144#endif
145 }
146
147 break;
148 }
149
150 case 3:
151 {
152 const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
153 const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
154 const std::vector<std::vector<OutputShape>> & dphidzeta = fe.get_dphidzeta();
155
156 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
157 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
158 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
159
160 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
161 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
162 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
163
164 const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
165 const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
166 const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
167
168 for (auto i : index_range(dphi))
169 for (auto p : index_range(dphi[i]))
170 {
171 // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
172 dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
173 dphideta[i][p]*detadx_map[p] +
174 dphidzeta[i][p]*dzetadx_map[p]);
175
176 // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
177 dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
178 dphideta[i][p]*detady_map[p] +
179 dphidzeta[i][p]*dzetady_map[p]);
180
181 // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
182 dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
183 dphideta[i][p]*detadz_map[p] +
184 dphidzeta[i][p]*dzetadz_map[p]);
185 }
186 break;
187 }
188
189 default:
190 libmesh_error_msg("Invalid dim = " << dim);
191 } // switch(dim)
192}
193
194#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
195template<typename OutputShape>
196void H1FETransformation<OutputShape>::map_d2phi(const unsigned int dim,
197 const std::vector<Point> &,
198 const FEGenericBase<OutputShape> & fe,
199 std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> & d2phi,
200 std::vector<std::vector<OutputShape>> & d2phidx2,
201 std::vector<std::vector<OutputShape>> & d2phidxdy,
202 std::vector<std::vector<OutputShape>> & d2phidxdz,
203 std::vector<std::vector<OutputShape>> & d2phidy2,
204 std::vector<std::vector<OutputShape>> & d2phidydz,
205 std::vector<std::vector<OutputShape>> & d2phidz2) const
206{
207 switch (dim)
208 {
209 case 0: // No derivatives in 0D
210 {
211 for (auto i : index_range(d2phi))
212 for (auto p : index_range(d2phi[i]))
213 d2phi[i][p] = 0.;
214
215 break;
216 }
217
218 case 1:
219 {
220 const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.get_d2phidxi2();
221
222 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
223#if LIBMESH_DIM>1
224 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
225#endif
226#if LIBMESH_DIM>2
227 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
228#endif
229
230 // Shape function derivatives in reference space
231 const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
232
233 // Inverse map second derivatives
234 const std::vector<std::vector<Real>> & d2xidxyz2 = fe.get_fe_map().get_d2xidxyz2();
235
236 for (auto i : index_range(d2phi))
237 for (auto p : index_range(d2phi[i]))
238 {
239 // phi_{x x}
240 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
241 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + // (xi_x)^2 * phi_{xi xi}
242 d2xidxyz2[p][0]*dphidxi[i][p]; // xi_{x x} * phi_{xi}
243
244#if LIBMESH_DIM>1
245 // phi_{x y}
246 d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
247 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + // xi_x * xi_y * phi_{xi xi}
248 d2xidxyz2[p][1]*dphidxi[i][p]; // xi_{x y} * phi_{xi}
249#endif
250
251#if LIBMESH_DIM>2
252 // phi_{x z}
253 d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
254 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + // xi_x * xi_z * phi_{xi xi}
255 d2xidxyz2[p][2]*dphidxi[i][p]; // xi_{x z} * phi_{xi}
256#endif
257
258
259#if LIBMESH_DIM>1
260 // phi_{y y}
261 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
262 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + // (xi_y)^2 * phi_{xi xi}
263 d2xidxyz2[p][3]*dphidxi[i][p]; // xi_{y y} * phi_{xi}
264#endif
265
266#if LIBMESH_DIM>2
267 // phi_{y z}
268 d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
269 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + // xi_y * xi_z * phi_{xi xi}
270 d2xidxyz2[p][4]*dphidxi[i][p]; // xi_{y z} * phi_{xi}
271
272 // phi_{z z}
273 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
274 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + // (xi_z)^2 * phi_{xi xi}
275 d2xidxyz2[p][5]*dphidxi[i][p]; // xi_{z z} * phi_{xi}
276#endif
277 }
278 break;
279 }
280
281 case 2:
282 {
283 const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.get_d2phidxi2();
284 const std::vector<std::vector<OutputShape>> & d2phidxideta = fe.get_d2phidxideta();
285 const std::vector<std::vector<OutputShape>> & d2phideta2 = fe.get_d2phideta2();
286
287 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
288 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
289#if LIBMESH_DIM > 2
290 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
291#endif
292
293 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
294 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
295#if LIBMESH_DIM > 2
296 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
297#endif
298
299 // Shape function derivatives in reference space
300 const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
301 const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
302
303 // Inverse map second derivatives
304 const std::vector<std::vector<Real>> & d2xidxyz2 = fe.get_fe_map().get_d2xidxyz2();
305 const std::vector<std::vector<Real>> & d2etadxyz2 = fe.get_fe_map().get_d2etadxyz2();
306
307 for (auto i : index_range(d2phi))
308 for (auto p : index_range(d2phi[i]))
309 {
310 // phi_{x x}
311 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
312 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + // (xi_x)^2 * phi_{xi xi}
313 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] + // (eta_x)^2 * phi_{eta eta}
314 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] + // 2 * xi_x * eta_x * phi_{xi eta}
315 d2xidxyz2[p][0]*dphidxi[i][p] + // xi_{x x} * phi_{xi}
316 d2etadxyz2[p][0]*dphideta[i][p]; // eta_{x x} * phi_{eta}
317
318 // phi_{x y}
319 d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
320 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + // xi_x * xi_y * phi_{xi xi}
321 d2phideta2[i][p]*detadx_map[p]*detady_map[p] + // eta_x * eta_y * phi_{eta eta}
322 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) + // (xi_x*eta_y + eta_x*xi_y) * phi_{xi eta}
323 d2xidxyz2[p][1]*dphidxi[i][p] + // xi_{x y} * phi_{xi}
324 d2etadxyz2[p][1]*dphideta[i][p]; // eta_{x y} * phi_{eta}
325
326#if LIBMESH_DIM > 2
327 // phi_{x z}
328 d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
329 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + // xi_x * xi_z * phi_{xi xi}
330 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] + // eta_x * eta_z * phi_{eta eta}
331 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) + // (xi_x*eta_z + eta_x*xi_z) * phi_{xi eta}
332 d2xidxyz2[p][2]*dphidxi[i][p] + // xi_{x z} * phi_{xi}
333 d2etadxyz2[p][2]*dphideta[i][p]; // eta_{x z} * phi_{eta}
334#endif
335
336 // phi_{y y}
337 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
338 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + // (xi_y)^2 * phi_{xi xi}
339 d2phideta2[i][p]*detady_map[p]*detady_map[p] + // (eta_y)^2 * phi_{eta eta}
340 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] + // 2 * xi_y * eta_y * phi_{xi eta}
341 d2xidxyz2[p][3]*dphidxi[i][p] + // xi_{y y} * phi_{xi}
342 d2etadxyz2[p][3]*dphideta[i][p]; // eta_{y y} * phi_{eta}
343
344#if LIBMESH_DIM > 2
345 // phi_{y z}
346 d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
347 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + // xi_y * xi_z * phi_{xi xi}
348 d2phideta2[i][p]*detady_map[p]*detadz_map[p] + // eta_y * eta_z * phi_{eta eta}
349 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) + // (xi_y*eta_z + eta_y*xi_z) * phi_{xi eta}
350 d2xidxyz2[p][4]*dphidxi[i][p] + // xi_{y z} * phi_{xi}
351 d2etadxyz2[p][4]*dphideta[i][p]; // eta_{y z} * phi_{eta}
352
353 // phi_{z z}
354 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
355 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + // (xi_z)^2 * phi_{xi xi}
356 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] + // (eta_z)^2 * phi_{eta eta}
357 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] + // 2 * xi_z * eta_z * phi_{xi eta}
358 d2xidxyz2[p][5]*dphidxi[i][p] + // xi_{z z} * phi_{xi}
359 d2etadxyz2[p][5]*dphideta[i][p]; // eta_{z z} * phi_{eta}
360#endif
361 }
362
363 break;
364 }
365
366 case 3:
367 {
368 const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.get_d2phidxi2();
369 const std::vector<std::vector<OutputShape>> & d2phidxideta = fe.get_d2phidxideta();
370 const std::vector<std::vector<OutputShape>> & d2phideta2 = fe.get_d2phideta2();
371 const std::vector<std::vector<OutputShape>> & d2phidxidzeta = fe.get_d2phidxidzeta();
372 const std::vector<std::vector<OutputShape>> & d2phidetadzeta = fe.get_d2phidetadzeta();
373 const std::vector<std::vector<OutputShape>> & d2phidzeta2 = fe.get_d2phidzeta2();
374
375 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
376 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
377 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
378
379 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
380 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
381 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
382
383 const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
384 const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
385 const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
386
387 // Shape function derivatives in reference space
388 const std::vector<std::vector<OutputShape>> & dphidxi = fe.get_dphidxi();
389 const std::vector<std::vector<OutputShape>> & dphideta = fe.get_dphideta();
390 const std::vector<std::vector<OutputShape>> & dphidzeta = fe.get_dphidzeta();
391
392 // Inverse map second derivatives
393 const std::vector<std::vector<Real>> & d2xidxyz2 = fe.get_fe_map().get_d2xidxyz2();
394 const std::vector<std::vector<Real>> & d2etadxyz2 = fe.get_fe_map().get_d2etadxyz2();
395 const std::vector<std::vector<Real>> & d2zetadxyz2 = fe.get_fe_map().get_d2zetadxyz2();
396
397 for (auto i : index_range(d2phi))
398 for (auto p : index_range(d2phi[i]))
399 {
400 // phi_{x x}
401 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
402 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] + // (xi_x)^2 * phi_{xi xi}
403 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] + // (eta_x)^2 * phi_{eta eta}
404 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p] + // (zeta_x)^2 * phi_{zeta zeta}
405 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] + // 2 * xi_x * eta_x * phi_{xi eta}
406 2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] + // 2 * xi_x * zeta_x * phi_{xi zeta}
407 2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] + // 2 * eta_x * zeta_x * phi_{eta zeta}
408 d2xidxyz2[p][0]*dphidxi[i][p] + // xi_{x x} * phi_{xi}
409 d2etadxyz2[p][0]*dphideta[i][p] + // eta_{x x} * phi_{eta}
410 d2zetadxyz2[p][0]*dphidzeta[i][p]; // zeta_{x x} * phi_{zeta}
411
412 // phi_{x y}
413 d2phi[i][p].slice(0).slice(1) = d2phi[i][p].slice(1).slice(0) = d2phidxdy[i][p] =
414 d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] + // xi_x * xi_y * phi_{xi xi}
415 d2phideta2[i][p]*detadx_map[p]*detady_map[p] + // eta_x * eta_y * phi_{eta eta}
416 d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] + // zeta_x * zeta_y * phi_{zeta zeta}
417 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) + // (xi_x*eta_y + eta_x*xi_y) * phi_{xi eta}
418 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetady_map[p] + dzetadx_map[p]*dxidy_map[p]) + // (zeta_x*xi_y + xi_x*zeta_y) * phi_{xi zeta}
419 d2phidetadzeta[i][p]*(detadx_map[p]*dzetady_map[p] + dzetadx_map[p]*detady_map[p]) + // (zeta_x*eta_y + eta_x*zeta_y) * phi_{eta zeta}
420 d2xidxyz2[p][1]*dphidxi[i][p] + // xi_{x y} * phi_{xi}
421 d2etadxyz2[p][1]*dphideta[i][p] + // eta_{x y} * phi_{eta}
422 d2zetadxyz2[p][1]*dphidzeta[i][p]; // zeta_{x y} * phi_{zeta}
423
424 // phi_{x z}
425 d2phi[i][p].slice(0).slice(2) = d2phi[i][p].slice(2).slice(0) = d2phidxdz[i][p] =
426 d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] + // xi_x * xi_z * phi_{xi xi}
427 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] + // eta_x * eta_z * phi_{eta eta}
428 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] + // zeta_x * zeta_z * phi_{zeta zeta}
429 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) + // (xi_x*eta_z + eta_x*xi_z) * phi_{xi eta}
430 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetadz_map[p] + dzetadx_map[p]*dxidz_map[p]) + // (zeta_x*xi_z + xi_x*zeta_z) * phi_{xi zeta}
431 d2phidetadzeta[i][p]*(detadx_map[p]*dzetadz_map[p] + dzetadx_map[p]*detadz_map[p]) + // (zeta_x*eta_z + eta_x*zeta_z) * phi_{eta zeta}
432 d2xidxyz2[p][2]*dphidxi[i][p] + // xi_{x z} * phi_{xi}
433 d2etadxyz2[p][2]*dphideta[i][p] + // eta_{x z} * phi_{eta}
434 d2zetadxyz2[p][2]*dphidzeta[i][p]; // zeta_{x z} * phi_{zeta}
435
436 // phi_{y y}
437 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
438 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] + // (xi_y)^2 * phi_{xi xi}
439 d2phideta2[i][p]*detady_map[p]*detady_map[p] + // (eta_y)^2 * phi_{eta eta}
440 d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p] + // (zeta_y)^2 * phi_{zeta zeta}
441 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] + // 2 * xi_y * eta_y * phi_{xi eta}
442 2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] + // 2 * xi_y * zeta_y * phi_{xi zeta}
443 2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] + // 2 * eta_y * zeta_y * phi_{eta zeta}
444 d2xidxyz2[p][3]*dphidxi[i][p] + // xi_{y y} * phi_{xi}
445 d2etadxyz2[p][3]*dphideta[i][p] + // eta_{y y} * phi_{eta}
446 d2zetadxyz2[p][3]*dphidzeta[i][p]; // zeta_{y y} * phi_{zeta}
447
448 // phi_{y z}
449 d2phi[i][p].slice(1).slice(2) = d2phi[i][p].slice(2).slice(1) = d2phidydz[i][p] =
450 d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] + // xi_y * xi_z * phi_{xi xi}
451 d2phideta2[i][p]*detady_map[p]*detadz_map[p] + // eta_y * eta_z * phi_{eta eta}
452 d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] + // zeta_y * zeta_z * phi_{zeta zeta}
453 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) + // (xi_y*eta_z + eta_y*xi_z) * phi_{xi eta}
454 d2phidxidzeta[i][p]*(dxidy_map[p]*dzetadz_map[p] + dzetady_map[p]*dxidz_map[p]) + // (zeta_y*xi_z + xi_y*zeta_z) * phi_{xi zeta}
455 d2phidetadzeta[i][p]*(detady_map[p]*dzetadz_map[p] + dzetady_map[p]*detadz_map[p]) + // (zeta_y*eta_z + eta_y*zeta_z) * phi_{eta zeta}
456 d2xidxyz2[p][4]*dphidxi[i][p] + // xi_{y z} * phi_{xi}
457 d2etadxyz2[p][4]*dphideta[i][p] + // eta_{y z} * phi_{eta}
458 d2zetadxyz2[p][4]*dphidzeta[i][p]; // zeta_{y z} * phi_{zeta}
459
460 // phi_{z z}
461 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
462 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] + // (xi_z)^2 * phi_{xi xi}
463 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] + // (eta_z)^2 * phi_{eta eta}
464 d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p] + // (zeta_z)^2 * phi_{zeta zeta}
465 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] + // 2 * xi_z * eta_z * phi_{xi eta}
466 2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] + // 2 * xi_z * zeta_z * phi_{xi zeta}
467 2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] + // 2 * eta_z * zeta_z * phi_{eta zeta}
468 d2xidxyz2[p][5]*dphidxi[i][p] + // xi_{z z} * phi_{xi}
469 d2etadxyz2[p][5]*dphideta[i][p] + // eta_{z z} * phi_{eta}
470 d2zetadxyz2[p][5]*dphidzeta[i][p]; // zeta_{z z} * phi_{zeta}
471 }
472
473 break;
474 }
475
476 default:
477 libmesh_error_msg("Invalid dim = " << dim);
478 } // switch(dim)
479}
480#endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
481
482template<>
483void H1FETransformation<Real>::map_curl(const unsigned int,
484 const Elem * const,
485 const std::vector<Point> &,
486 const FEGenericBase<Real> &,
487 std::vector<std::vector<Real>> &) const
488{
489 libmesh_error_msg("Computing the curl of a shape function only \nmakes sense for vector-valued elements.");
490}
491
492template<>
493void H1FETransformation<RealGradient>::map_curl(const unsigned int dim,
494 const Elem * const,
495 const std::vector<Point> &,
496 const FEGenericBase<RealGradient> & fe,
497 std::vector<std::vector<RealGradient>> & curl_phi ) const
498{
499 switch (dim)
500 {
501 case 0:
502 case 1:
503 libmesh_error_msg("The curl only make sense in 2D and 3D");
504
505 case 2:
506 {
507 const std::vector<std::vector<RealGradient>> & dphidxi = fe.get_dphidxi();
508 const std::vector<std::vector<RealGradient>> & dphideta = fe.get_dphideta();
509
510 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
511 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
512#if LIBMESH_DIM > 2
513 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
514#endif
515
516 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
517 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
518#if LIBMESH_DIM > 2
519 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
520#endif
521
522 // For 2D elements in 3D space:
523 // curl = ( -dphi_y/dz, dphi_x/dz, dphi_y/dx - dphi_x/dy )
524 for (auto i : index_range(curl_phi))
525 for (auto p : index_range(curl_phi[i]))
526 {
527
528 Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
529 + (dphideta[i][p].slice(1))*detadx_map[p];
530
531 Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
532 + (dphideta[i][p].slice(0))*detady_map[p];
533
534 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
535
536#if LIBMESH_DIM > 2
537 Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
538 + (dphideta[i][p].slice(1))*detadz_map[p];
539
540 Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
541 + (dphideta[i][p].slice(0))*detadz_map[p];
542
543 curl_phi[i][p].slice(0) = -dphiy_dz;
544 curl_phi[i][p].slice(1) = dphix_dz;
545#endif
546 }
547
548 break;
549 }
550 case 3:
551 {
552 const std::vector<std::vector<RealGradient>> & dphidxi = fe.get_dphidxi();
553 const std::vector<std::vector<RealGradient>> & dphideta = fe.get_dphideta();
554 const std::vector<std::vector<RealGradient>> & dphidzeta = fe.get_dphidzeta();
555
556 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
557 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
558 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
559
560 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
561 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
562 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
563
564 const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
565 const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
566 const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
567
568 // In 3D: curl = ( dphi_z/dy - dphi_y/dz, dphi_x/dz - dphi_z/dx, dphi_y/dx - dphi_x/dy )
569 for (auto i : index_range(curl_phi))
570 for (auto p : index_range(curl_phi[i]))
571 {
572 Real dphiz_dy = (dphidxi[i][p].slice(2))*dxidy_map[p]
573 + (dphideta[i][p].slice(2))*detady_map[p]
574 + (dphidzeta[i][p].slice(2))*dzetady_map[p];
575
576 Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
577 + (dphideta[i][p].slice(1))*detadz_map[p]
578 + (dphidzeta[i][p].slice(1))*dzetadz_map[p];
579
580 Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
581 + (dphideta[i][p].slice(0))*detadz_map[p]
582 + (dphidzeta[i][p].slice(0))*dzetadz_map[p];
583
584 Real dphiz_dx = (dphidxi[i][p].slice(2))*dxidx_map[p]
585 + (dphideta[i][p].slice(2))*detadx_map[p]
586 + (dphidzeta[i][p].slice(2))*dzetadx_map[p];
587
588 Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
589 + (dphideta[i][p].slice(1))*detadx_map[p]
590 + (dphidzeta[i][p].slice(1))*dzetadx_map[p];
591
592 Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
593 + (dphideta[i][p].slice(0))*detady_map[p]
594 + (dphidzeta[i][p].slice(0))*dzetady_map[p];
595
596 curl_phi[i][p].slice(0) = dphiz_dy - dphiy_dz;
597
598 curl_phi[i][p].slice(1) = dphix_dz - dphiz_dx;
599
600 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
601 }
602
603 break;
604 }
605 default:
606 libmesh_error_msg("Invalid dim = " << dim);
607 }
608}
609
610
611template<>
612void H1FETransformation<Real>::map_div(const unsigned int,
613 const Elem * const,
614 const std::vector<Point> &,
615 const FEGenericBase<Real> &,
616 std::vector<std::vector<FEGenericBase<Real>::OutputDivergence>> & ) const
617{
618 libmesh_error_msg("Computing the divergence of a shape function only \nmakes sense for vector-valued elements.");
619}
620
621
622template<>
623void H1FETransformation<RealGradient>::map_div(const unsigned int dim,
624 const Elem * const,
625 const std::vector<Point> &,
626 const FEGenericBase<RealGradient> & fe,
627 std::vector<std::vector<FEGenericBase<RealGradient>::OutputDivergence>> & div_phi) const
628{
629 switch (dim)
630 {
631 case 0:
632 case 1:
633 libmesh_error_msg("The divergence only make sense in 2D and 3D");
634
635 case 2:
636 {
637 const std::vector<std::vector<RealGradient>> & dphidxi = fe.get_dphidxi();
638 const std::vector<std::vector<RealGradient>> & dphideta = fe.get_dphideta();
639
640 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
641 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
642
643 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
644 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
645
646 // In 2D: div = dphi_x/dx + dphi_y/dy
647 for (auto i : index_range(div_phi))
648 for (auto p : index_range(div_phi[i]))
649 {
650 Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
651 + (dphideta[i][p].slice(0))*detadx_map[p];
652
653 Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
654 + (dphideta[i][p].slice(1))*detady_map[p];
655
656 div_phi[i][p] = dphix_dx + dphiy_dy;
657 }
658 break;
659 }
660 case 3:
661 {
662 const std::vector<std::vector<RealGradient>> & dphidxi = fe.get_dphidxi();
663 const std::vector<std::vector<RealGradient>> & dphideta = fe.get_dphideta();
664 const std::vector<std::vector<RealGradient>> & dphidzeta = fe.get_dphidzeta();
665
666 const std::vector<Real> & dxidx_map = fe.get_fe_map().get_dxidx();
667 const std::vector<Real> & dxidy_map = fe.get_fe_map().get_dxidy();
668 const std::vector<Real> & dxidz_map = fe.get_fe_map().get_dxidz();
669
670 const std::vector<Real> & detadx_map = fe.get_fe_map().get_detadx();
671 const std::vector<Real> & detady_map = fe.get_fe_map().get_detady();
672 const std::vector<Real> & detadz_map = fe.get_fe_map().get_detadz();
673
674 const std::vector<Real> & dzetadx_map = fe.get_fe_map().get_dzetadx();
675 const std::vector<Real> & dzetady_map = fe.get_fe_map().get_dzetady();
676 const std::vector<Real> & dzetadz_map = fe.get_fe_map().get_dzetadz();
677
678 // In 3D: div = dphi_x/dx + dphi_y/dy + dphi_z/dz
679 for (auto i : index_range(div_phi))
680 for (auto p : index_range(div_phi[i]))
681 {
682 Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
683 + (dphideta[i][p].slice(0))*detadx_map[p]
684 + (dphidzeta[i][p].slice(0))*dzetadx_map[p];
685
686 Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
687 + (dphideta[i][p].slice(1))*detady_map[p]
688 + (dphidzeta[i][p].slice(1))*dzetady_map[p];
689
690 Real dphiz_dz = (dphidxi[i][p].slice(2))*dxidz_map[p]
691 + (dphideta[i][p].slice(2))*detadz_map[p]
692 + (dphidzeta[i][p].slice(2))*dzetadz_map[p];
693
694 div_phi[i][p] = dphix_dx + dphiy_dy + dphiz_dz;
695 }
696
697 break;
698 }
699
700 default:
701 libmesh_error_msg("Invalid dim = " << dim);
702 } // switch(dim)
703
704 return;
705}
706
707// Explicit Instantiations
708template class LIBMESH_EXPORT H1FETransformation<Real>;
709template class LIBMESH_EXPORT H1FETransformation<RealGradient>;
710
711} //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::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_d2phidxidzeta
const std::vector< std::vector< OutputShape > > & get_d2phidxidzeta() const
Definition fe_base.h:395
libMesh::FEGenericBase::get_dphidxi
const std::vector< std::vector< OutputShape > > & get_dphidxi() const
Definition fe_base.h:293
libMesh::FEGenericBase::get_dphidzeta
const std::vector< std::vector< OutputShape > > & get_dphidzeta() const
Definition fe_base.h:309
libMesh::FEGenericBase::get_d2phidxi2
const std::vector< std::vector< OutputShape > > & get_d2phidxi2() const
Definition fe_base.h:379
libMesh::FEGenericBase::get_d2phideta2
const std::vector< std::vector< OutputShape > > & get_d2phideta2() const
Definition fe_base.h:403
libMesh::FEGenericBase::get_d2phidzeta2
const std::vector< std::vector< OutputShape > > & get_d2phidzeta2() const
Definition fe_base.h:419
libMesh::FEGenericBase::get_d2phidetadzeta
const std::vector< std::vector< OutputShape > > & get_d2phidetadzeta() const
Definition fe_base.h:411
libMesh::FEGenericBase::OutputDivergence
TensorTools::DecrementRank< OutputShape >::type OutputDivergence
Definition fe_base.h:122
libMesh::FEGenericBase::get_d2phidxideta
const std::vector< std::vector< OutputShape > > & get_d2phidxideta() const
Definition fe_base.h:387
libMesh::FEGenericBase::get_dphideta
const std::vector< std::vector< OutputShape > > & get_dphideta() const
Definition fe_base.h:301
libMesh::FEMap::get_dxidz
const std::vector< Real > & get_dxidz() const
Definition fe_map.h:325
libMesh::FEMap::get_detadx
const std::vector< Real > & get_detadx() const
Definition fe_map.h:333
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_d2zetadxyz2
const std::vector< std::vector< Real > > & get_d2zetadxyz2() const
Second derivatives of "zeta" reference coordinate wrt physical coordinates.
Definition fe_map.h:395
libMesh::FEMap::get_dxidx
const std::vector< Real > & get_dxidx() const
Definition fe_map.h:309
libMesh::FEMap::get_d2etadxyz2
const std::vector< std::vector< Real > > & get_d2etadxyz2() const
Second derivatives of "eta" reference coordinate wrt physical coordinates.
Definition fe_map.h:388
libMesh::FEMap::get_dzetadx
const std::vector< Real > & get_dzetadx() const
Definition fe_map.h:357
libMesh::FEMap::get_d2xidxyz2
const std::vector< std::vector< Real > > & get_d2xidxyz2() const
Second derivatives of "xi" reference coordinate wrt physical coordinates.
Definition fe_map.h:381
libMesh::FEMap::get_dxidy
const std::vector< Real > & get_dxidy() const
Definition fe_map.h:317
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::H1FETransformation
This class handles the computation of the shape functions in the physical domain for H1 conforming el...
Definition h1_fe_transformation.h:41
libMesh::H1FETransformation::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 H1 conforming elements.
Definition h1_fe_transformation.C:64
libMesh::H1FETransformation::map_d2phi
virtual void map_d2phi(const unsigned int dim, const std::vector< Point > &qp, const FEGenericBase< OutputShape > &fe, std::vector< std::vector< typename FEGenericBase< OutputShape >::OutputTensor > > &d2phi, std::vector< std::vector< OutputShape > > &d2phidx2, std::vector< std::vector< OutputShape > > &d2phidxdy, std::vector< std::vector< OutputShape > > &d2phidxdz, std::vector< std::vector< OutputShape > > &d2phidy2, std::vector< std::vector< OutputShape > > &d2phidydz, std::vector< std::vector< OutputShape > > &d2phidz2) const override
Evaluates shape function Hessians in physical coordinates based on H1 conforming finite element trans...
Definition h1_fe_transformation.C:196
libMesh::H1FETransformation::init_map_phi
virtual void init_map_phi(const FEGenericBase< OutputShape > &fe) const override
Pre-requests any necessary data from FEMap.
Definition h1_fe_transformation.C:27
libMesh::H1FETransformation::map_curl
virtual void map_curl(const unsigned int dim, const Elem *const elem, const std::vector< Point > &qp, const FEGenericBase< OutputShape > &fe, std::vector< std::vector< OutputShape > > &curl_phi) const override
Evaluates the shape function curl in physical coordinates based on H1 conforming finite element trans...
libMesh::H1FETransformation::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 H1 conforming finite element...
libMesh::H1FETransformation::map_phi
virtual void map_phi(const unsigned int, const Elem *const, const std::vector< Point > &, const FEGenericBase< OutputShape > &, std::vector< std::vector< OutputShape > > &, bool add_p_level=true) const override
Evaluates shape functions in physical coordinates for H1 conforming elements.
Definition h1_fe_transformation.C:52
libMesh::H1FETransformation::init_map_dphi
virtual void init_map_dphi(const FEGenericBase< OutputShape > &fe) const override
Pre-requests any necessary data from FEMap.
Definition h1_fe_transformation.C:35
libMesh::H1FETransformation::init_map_d2phi
virtual void init_map_d2phi(const FEGenericBase< OutputShape > &fe) const override
Pre-requests any necessary data from FEMap.
Definition h1_fe_transformation.C:43
libMesh::VectorValue
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
Definition vector_value.h:48
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::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Fri Aug 21 2026 13:32:29 for libMesh by doxygen 1.9.8