26 template<
typename OutputShape>
34 template<
typename OutputShape>
42 template<
typename OutputShape>
46 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 51 template<
typename OutputShape>
53 const Elem *
const elem,
54 const std::vector<Point> & qp,
56 std::vector<std::vector<OutputShape>> & phi )
const 64 libmesh_assert_equal_to ( qp.size(), phi[i].size() );
66 FEInterface::shape<OutputShape>(0, fe.
get_fe_type(), elem, i, qp[p], phi[i][p]);
74 libmesh_assert_equal_to ( qp.size(), phi[i].size() );
76 FEInterface::shape<OutputShape>(1, fe.
get_fe_type(), elem, i, qp[p], phi[i][p]);
84 libmesh_assert_equal_to ( qp.size(), phi[i].size() );
86 FEInterface::shape<OutputShape>(2, fe.
get_fe_type(), elem, i, qp[p], phi[i][p]);
94 libmesh_assert_equal_to ( qp.size(), phi[i].size() );
96 FEInterface::shape<OutputShape>(3, fe.
get_fe_type(), elem, i, qp[p], phi[i][p]);
102 libmesh_error_msg(
"Invalid dim = " << dim);
107 template<
typename OutputShape>
110 const std::vector<Point> &,
113 std::vector<std::vector<OutputShape>> & dphidx,
114 std::vector<std::vector<OutputShape>> & dphidy,
115 std::vector<std::vector<OutputShape>> & dphidz )
const 129 const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
143 dphi[i][p].slice(0) = dphidx[i][p] = dphidxi[i][p]*dxidx_map[p];
146 dphi[i][p].slice(1) = dphidy[i][p] = dphidxi[i][p]*dxidy_map[p];
149 dphi[i][p].slice(2) = dphidz[i][p] = dphidxi[i][p]*dxidz_map[p];
158 const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
159 const std::vector<std::vector<OutputShape>> & dphideta = fe.
get_dphideta();
177 dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
178 dphideta[i][p]*detadx_map[p]);
181 dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
182 dphideta[i][p]*detady_map[p]);
186 dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
187 dphideta[i][p]*detadz_map[p]);
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();
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]);
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]);
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]);
234 libmesh_error_msg(
"Invalid dim = " << dim);
238 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 239 template<
typename OutputShape>
241 const std::vector<Point> &,
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 264 const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.
get_d2phidxi2();
275 const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
284 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
285 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
286 d2xidxyz2[p][0]*dphidxi[i][p];
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] +
292 d2xidxyz2[p][1]*dphidxi[i][p];
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] +
299 d2xidxyz2[p][2]*dphidxi[i][p];
305 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
306 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
307 d2xidxyz2[p][3]*dphidxi[i][p];
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] +
314 d2xidxyz2[p][4]*dphidxi[i][p];
317 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
318 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
319 d2xidxyz2[p][5]*dphidxi[i][p];
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();
344 const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
345 const std::vector<std::vector<OutputShape>> & dphideta = fe.
get_dphideta();
355 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
356 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
357 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +
358 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +
359 d2xidxyz2[p][0]*dphidxi[i][p] +
360 d2etadxyz2[p][0]*dphideta[i][p];
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] +
365 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +
366 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) +
367 d2xidxyz2[p][1]*dphidxi[i][p] +
368 d2etadxyz2[p][1]*dphideta[i][p];
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] +
374 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +
375 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) +
376 d2xidxyz2[p][2]*dphidxi[i][p] +
377 d2etadxyz2[p][2]*dphideta[i][p];
381 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
382 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
383 d2phideta2[i][p]*detady_map[p]*detady_map[p] +
384 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +
385 d2xidxyz2[p][3]*dphidxi[i][p] +
386 d2etadxyz2[p][3]*dphideta[i][p];
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] +
392 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +
393 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) +
394 d2xidxyz2[p][4]*dphidxi[i][p] +
395 d2etadxyz2[p][4]*dphideta[i][p];
398 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
399 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
400 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +
401 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +
402 d2xidxyz2[p][5]*dphidxi[i][p] +
403 d2etadxyz2[p][5]*dphideta[i][p];
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();
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();
445 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
446 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
447 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +
448 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p] +
449 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +
450 2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] +
451 2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] +
452 d2xidxyz2[p][0]*dphidxi[i][p] +
453 d2etadxyz2[p][0]*dphideta[i][p] +
454 d2zetadxyz2[p][0]*dphidzeta[i][p];
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] +
459 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +
460 d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] +
461 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) +
462 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetady_map[p] + dzetadx_map[p]*dxidy_map[p]) +
463 d2phidetadzeta[i][p]*(detadx_map[p]*dzetady_map[p] + dzetadx_map[p]*detady_map[p]) +
464 d2xidxyz2[p][1]*dphidxi[i][p] +
465 d2etadxyz2[p][1]*dphideta[i][p] +
466 d2zetadxyz2[p][1]*dphidzeta[i][p];
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] +
471 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +
472 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] +
473 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) +
474 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetadz_map[p] + dzetadx_map[p]*dxidz_map[p]) +
475 d2phidetadzeta[i][p]*(detadx_map[p]*dzetadz_map[p] + dzetadx_map[p]*detadz_map[p]) +
476 d2xidxyz2[p][2]*dphidxi[i][p] +
477 d2etadxyz2[p][2]*dphideta[i][p] +
478 d2zetadxyz2[p][2]*dphidzeta[i][p];
481 d2phi[i][p].slice(1).slice(1) = d2phidxdz[i][p] =
482 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
483 d2phideta2[i][p]*detady_map[p]*detady_map[p] +
484 d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p] +
485 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +
486 2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] +
487 2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] +
488 d2xidxyz2[p][3]*dphidxi[i][p] +
489 d2etadxyz2[p][3]*dphideta[i][p] +
490 d2zetadxyz2[p][3]*dphidzeta[i][p];
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] +
495 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +
496 d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] +
497 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) +
498 d2phidxidzeta[i][p]*(dxidy_map[p]*dzetadz_map[p] + dzetady_map[p]*dxidz_map[p]) +
499 d2phidetadzeta[i][p]*(detady_map[p]*dzetadz_map[p] + dzetady_map[p]*detadz_map[p]) +
500 d2xidxyz2[p][4]*dphidxi[i][p] +
501 d2etadxyz2[p][4]*dphideta[i][p] +
502 d2zetadxyz2[p][4]*dphidzeta[i][p];
505 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
506 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
507 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +
508 d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p] +
509 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +
510 2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] +
511 2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] +
512 d2xidxyz2[p][5]*dphidxi[i][p] +
513 d2etadxyz2[p][5]*dphideta[i][p] +
514 d2zetadxyz2[p][5]*dphidzeta[i][p];
521 libmesh_error_msg(
"Invalid dim = " << dim);
524 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES 529 const std::vector<Point> &,
531 std::vector<std::vector<Real>> &)
const 533 libmesh_error_msg(
"Computing the curl of a shape function only \nmakes sense for vector-valued elements.");
539 const std::vector<Point> &,
541 std::vector<std::vector<RealGradient>> & curl_phi )
const 547 libmesh_error_msg(
"The curl only make sense in 2D and 3D");
551 const std::vector<std::vector<RealGradient>> & dphidxi = fe.
get_dphidxi();
552 const std::vector<std::vector<RealGradient>> & dphideta = fe.
get_dphideta();
572 Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
573 + (dphideta[i][p].slice(1))*detadx_map[p];
575 Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
576 + (dphideta[i][p].slice(0))*detady_map[p];
578 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
581 Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
582 + (dphideta[i][p].slice(1))*detadz_map[p];
584 Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
585 + (dphideta[i][p].slice(0))*detadz_map[p];
587 curl_phi[i][p].slice(0) = -dphiy_dz;
588 curl_phi[i][p].slice(1) = dphix_dz;
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();
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];
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];
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];
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];
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];
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];
640 curl_phi[i][p].slice(0) = dphiz_dy - dphiy_dz;
642 curl_phi[i][p].slice(1) = dphix_dz - dphiz_dx;
644 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
650 libmesh_error_msg(
"Invalid dim = " << dim);
658 const std::vector<Point> &,
662 libmesh_error_msg(
"Computing the divergence of a shape function only \nmakes sense for vector-valued elements.");
669 const std::vector<Point> &,
677 libmesh_error_msg(
"The divergence only make sense in 2D and 3D");
681 const std::vector<std::vector<RealGradient>> & dphidxi = fe.
get_dphidxi();
682 const std::vector<std::vector<RealGradient>> & dphideta = fe.
get_dphideta();
694 Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
695 + (dphideta[i][p].slice(0))*detadx_map[p];
697 Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
698 + (dphideta[i][p].slice(1))*detady_map[p];
700 div_phi[i][p] = dphix_dx + dphiy_dy;
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();
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];
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];
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];
738 div_phi[i][p] = dphix_dx + dphiy_dy + dphiz_dz;
745 libmesh_error_msg(
"Invalid dim = " << dim);
const std::vector< std::vector< OutputShape > > & get_d2phideta2() const
const std::vector< std::vector< OutputShape > > & get_d2phidxideta() const
const std::vector< Real > & get_detadz() const
The base class for all geometric element types.
IntRange< std::size_t > index_range(const std::vector< T > &vec)
const std::vector< std::vector< OutputShape > > & get_dphideta() const
TensorTools::DecrementRank< OutputShape >::type OutputDivergence
const std::vector< std::vector< Real > > & get_d2etadxyz2() const
const std::vector< Real > & get_dxidz() const
const std::vector< std::vector< OutputShape > > & get_dphidzeta() const
FEType get_fe_type() const
const std::vector< std::vector< Real > > & get_d2zetadxyz2() const
const std::vector< Real > & get_dzetadx() const
const std::vector< Real > & get_dxidx() const
const std::vector< std::vector< OutputShape > > & get_d2phidxi2() const
const std::vector< Real > & get_dzetady() const
const std::vector< Real > & get_dxidy() const
const std::vector< std::vector< OutputShape > > & get_d2phidxidzeta() const
const std::vector< std::vector< OutputShape > > & get_d2phidetadzeta() const
const std::vector< Real > & get_dzetadz() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< Real > & get_detady() const
const std::vector< std::vector< OutputShape > > & get_d2phidzeta2() const
const std::vector< std::vector< OutputShape > > & get_dphidxi() const
const FEMap & get_fe_map() const
const std::vector< std::vector< Real > > & get_d2xidxyz2() const
const std::vector< Real > & get_detadx() const