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