libMesh
h1_fe_transformation.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 
24 namespace libMesh
25 {
26 template<typename OutputShape>
28 {
29  // Nothing needed
30 }
31 
32 
33 
34 template<typename OutputShape>
36 {
37  fe.get_fe_map().get_dxidx();
38 }
39 
40 
41 
42 template<typename OutputShape>
44 {
45  fe.get_fe_map().get_dxidx();
46 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
48 #endif
49 }
50 
51 template<typename OutputShape>
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 
63 template<typename OutputShape>
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
195 template<typename OutputShape>
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 
482 template<>
483 void 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 
492 template<>
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 
611 template<>
612 void 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 
622 template<>
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
708 template class LIBMESH_EXPORT H1FETransformation<Real>;
709 template class LIBMESH_EXPORT H1FETransformation<RealGradient>;
710 
711 } //namespace libMesh
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...
const std::vector< std::vector< OutputShape > > & get_d2phideta2() const
Definition: fe_base.h:403
unsigned int dim
const std::vector< std::vector< OutputShape > > & get_d2phidxideta() const
Definition: fe_base.h:387
const std::vector< Real > & get_detadz() const
Definition: fe_map.h:349
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
virtual void init_map_d2phi(const FEGenericBase< OutputShape > &fe) const override
Pre-requests any necessary data from FEMap.
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The libMesh namespace provides an interface to certain functionality in the library.
const std::vector< std::vector< OutputShape > > & get_dphideta() const
Definition: fe_base.h:301
TensorTools::DecrementRank< OutputShape >::type OutputDivergence
Definition: fe_base.h:122
const std::vector< std::vector< Real > > & get_d2etadxyz2() const
Second derivatives of "eta" reference coordinate wrt physical coordinates.
Definition: fe_map.h:388
This class handles the computation of the shape functions in the physical domain for H1 conforming el...
const std::vector< Real > & get_dxidz() const
Definition: fe_map.h:325
const std::vector< std::vector< OutputShape > > & get_dphidzeta() const
Definition: fe_base.h:309
FEType get_fe_type() const
Definition: fe_abstract.h:499
const std::vector< std::vector< Real > > & get_d2zetadxyz2() const
Second derivatives of "zeta" reference coordinate wrt physical coordinates.
Definition: fe_map.h:395
virtual void init_map_phi(const FEGenericBase< OutputShape > &fe) const override
Pre-requests any necessary data from FEMap.
const std::vector< Real > & get_dzetadx() const
Definition: fe_map.h:357
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.
const std::vector< Real > & get_dxidx() const
Definition: fe_map.h:309
const std::vector< std::vector< OutputShape > > & get_d2phidxi2() const
Definition: fe_base.h:379
const std::vector< Real > & get_dzetady() const
Definition: fe_map.h:365
const std::vector< Real > & get_dxidy() const
Definition: fe_map.h:317
const std::vector< std::vector< OutputShape > > & get_d2phidxidzeta() const
Definition: fe_base.h:395
const std::vector< std::vector< OutputShape > > & get_d2phidetadzeta() const
Definition: fe_base.h:411
const std::vector< Real > & get_dzetadz() const
Definition: fe_map.h:373
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...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< Real > & get_detady() const
Definition: fe_map.h:341
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...
const std::vector< std::vector< OutputShape > > & get_d2phidzeta2() const
Definition: fe_base.h:419
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.
const std::vector< std::vector< OutputShape > > & get_dphidxi() const
Definition: fe_base.h:293
const FEMap & get_fe_map() const
Definition: fe_abstract.h:533
virtual void init_map_dphi(const FEGenericBase< OutputShape > &fe) const override
Pre-requests any necessary data from FEMap.
const std::vector< std::vector< Real > > & get_d2xidxyz2() const
Second derivatives of "xi" reference coordinate wrt physical coordinates.
Definition: fe_map.h:381
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:111
This class forms the foundation from which generic finite elements may be derived.
const std::vector< Real > & get_detadx() const
Definition: fe_map.h:333