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