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" 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,
57 const bool add_p_level)
const 59 FEInterface::all_shapes<OutputShape>(
dim, fe.
get_fe_type(), elem, qp, phi, add_p_level);
63 template<
typename OutputShape>
66 const std::vector<Point> &,
69 std::vector<std::vector<OutputShape>> & dphidx,
70 std::vector<std::vector<OutputShape>> & dphidy,
71 std::vector<std::vector<OutputShape>> & dphidz )
const 85 const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
99 dphi[i][p].slice(0) = dphidx[i][p] = dphidxi[i][p]*dxidx_map[p];
102 dphi[i][p].slice(1) = dphidy[i][p] = dphidxi[i][p]*dxidy_map[p];
105 dphi[i][p].slice(2) = dphidz[i][p] = dphidxi[i][p]*dxidz_map[p];
114 const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
115 const std::vector<std::vector<OutputShape>> & dphideta = fe.
get_dphideta();
133 dphi[i][p].slice(0) = dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
134 dphideta[i][p]*detadx_map[p]);
137 dphi[i][p].slice(1) = dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
138 dphideta[i][p]*detady_map[p]);
142 dphi[i][p].slice(2) = dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
143 dphideta[i][p]*detadz_map[p]);
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();
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]);
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]);
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]);
190 libmesh_error_msg(
"Invalid dim = " <<
dim);
194 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 195 template<
typename OutputShape>
197 const std::vector<Point> &,
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 220 const std::vector<std::vector<OutputShape>> & d2phidxi2 = fe.
get_d2phidxi2();
231 const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
240 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
241 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
242 d2xidxyz2[p][0]*dphidxi[i][p];
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] +
248 d2xidxyz2[p][1]*dphidxi[i][p];
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] +
255 d2xidxyz2[p][2]*dphidxi[i][p];
261 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
262 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
263 d2xidxyz2[p][3]*dphidxi[i][p];
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] +
270 d2xidxyz2[p][4]*dphidxi[i][p];
273 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
274 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
275 d2xidxyz2[p][5]*dphidxi[i][p];
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();
300 const std::vector<std::vector<OutputShape>> & dphidxi = fe.
get_dphidxi();
301 const std::vector<std::vector<OutputShape>> & dphideta = fe.
get_dphideta();
311 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
312 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
313 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +
314 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +
315 d2xidxyz2[p][0]*dphidxi[i][p] +
316 d2etadxyz2[p][0]*dphideta[i][p];
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] +
321 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +
322 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) +
323 d2xidxyz2[p][1]*dphidxi[i][p] +
324 d2etadxyz2[p][1]*dphideta[i][p];
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] +
330 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +
331 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) +
332 d2xidxyz2[p][2]*dphidxi[i][p] +
333 d2etadxyz2[p][2]*dphideta[i][p];
337 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
338 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
339 d2phideta2[i][p]*detady_map[p]*detady_map[p] +
340 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +
341 d2xidxyz2[p][3]*dphidxi[i][p] +
342 d2etadxyz2[p][3]*dphideta[i][p];
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] +
348 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +
349 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) +
350 d2xidxyz2[p][4]*dphidxi[i][p] +
351 d2etadxyz2[p][4]*dphideta[i][p];
354 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
355 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
356 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +
357 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +
358 d2xidxyz2[p][5]*dphidxi[i][p] +
359 d2etadxyz2[p][5]*dphideta[i][p];
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();
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();
401 d2phi[i][p].slice(0).slice(0) = d2phidx2[i][p] =
402 d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
403 d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +
404 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p] +
405 2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +
406 2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] +
407 2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] +
408 d2xidxyz2[p][0]*dphidxi[i][p] +
409 d2etadxyz2[p][0]*dphideta[i][p] +
410 d2zetadxyz2[p][0]*dphidzeta[i][p];
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] +
415 d2phideta2[i][p]*detadx_map[p]*detady_map[p] +
416 d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] +
417 d2phidxideta[i][p]*(dxidx_map[p]*detady_map[p] + detadx_map[p]*dxidy_map[p]) +
418 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetady_map[p] + dzetadx_map[p]*dxidy_map[p]) +
419 d2phidetadzeta[i][p]*(detadx_map[p]*dzetady_map[p] + dzetadx_map[p]*detady_map[p]) +
420 d2xidxyz2[p][1]*dphidxi[i][p] +
421 d2etadxyz2[p][1]*dphideta[i][p] +
422 d2zetadxyz2[p][1]*dphidzeta[i][p];
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] +
427 d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +
428 d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] +
429 d2phidxideta[i][p]*(dxidx_map[p]*detadz_map[p] + detadx_map[p]*dxidz_map[p]) +
430 d2phidxidzeta[i][p]*(dxidx_map[p]*dzetadz_map[p] + dzetadx_map[p]*dxidz_map[p]) +
431 d2phidetadzeta[i][p]*(detadx_map[p]*dzetadz_map[p] + dzetadx_map[p]*detadz_map[p]) +
432 d2xidxyz2[p][2]*dphidxi[i][p] +
433 d2etadxyz2[p][2]*dphideta[i][p] +
434 d2zetadxyz2[p][2]*dphidzeta[i][p];
437 d2phi[i][p].slice(1).slice(1) = d2phidy2[i][p] =
438 d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
439 d2phideta2[i][p]*detady_map[p]*detady_map[p] +
440 d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p] +
441 2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +
442 2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] +
443 2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] +
444 d2xidxyz2[p][3]*dphidxi[i][p] +
445 d2etadxyz2[p][3]*dphideta[i][p] +
446 d2zetadxyz2[p][3]*dphidzeta[i][p];
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] +
451 d2phideta2[i][p]*detady_map[p]*detadz_map[p] +
452 d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] +
453 d2phidxideta[i][p]*(dxidy_map[p]*detadz_map[p] + detady_map[p]*dxidz_map[p]) +
454 d2phidxidzeta[i][p]*(dxidy_map[p]*dzetadz_map[p] + dzetady_map[p]*dxidz_map[p]) +
455 d2phidetadzeta[i][p]*(detady_map[p]*dzetadz_map[p] + dzetady_map[p]*detadz_map[p]) +
456 d2xidxyz2[p][4]*dphidxi[i][p] +
457 d2etadxyz2[p][4]*dphideta[i][p] +
458 d2zetadxyz2[p][4]*dphidzeta[i][p];
461 d2phi[i][p].slice(2).slice(2) = d2phidz2[i][p] =
462 d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
463 d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +
464 d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p] +
465 2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +
466 2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] +
467 2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] +
468 d2xidxyz2[p][5]*dphidxi[i][p] +
469 d2etadxyz2[p][5]*dphideta[i][p] +
470 d2zetadxyz2[p][5]*dphidzeta[i][p];
477 libmesh_error_msg(
"Invalid dim = " <<
dim);
480 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES 485 const std::vector<Point> &,
487 std::vector<std::vector<Real>> &)
const 489 libmesh_error_msg(
"Computing the curl of a shape function only \nmakes sense for vector-valued elements.");
495 const std::vector<Point> &,
497 std::vector<std::vector<RealGradient>> & curl_phi )
const 503 libmesh_error_msg(
"The curl only make sense in 2D and 3D");
507 const std::vector<std::vector<RealGradient>> & dphidxi = fe.
get_dphidxi();
508 const std::vector<std::vector<RealGradient>> & dphideta = fe.
get_dphideta();
528 Real dphiy_dx = (dphidxi[i][p].slice(1))*dxidx_map[p]
529 + (dphideta[i][p].slice(1))*detadx_map[p];
531 Real dphix_dy = (dphidxi[i][p].slice(0))*dxidy_map[p]
532 + (dphideta[i][p].slice(0))*detady_map[p];
534 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
537 Real dphiy_dz = (dphidxi[i][p].slice(1))*dxidz_map[p]
538 + (dphideta[i][p].slice(1))*detadz_map[p];
540 Real dphix_dz = (dphidxi[i][p].slice(0))*dxidz_map[p]
541 + (dphideta[i][p].slice(0))*detadz_map[p];
543 curl_phi[i][p].slice(0) = -dphiy_dz;
544 curl_phi[i][p].slice(1) = dphix_dz;
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();
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];
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];
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];
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];
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];
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];
596 curl_phi[i][p].slice(0) = dphiz_dy - dphiy_dz;
598 curl_phi[i][p].slice(1) = dphix_dz - dphiz_dx;
600 curl_phi[i][p].slice(2) = dphiy_dx - dphix_dy;
606 libmesh_error_msg(
"Invalid dim = " <<
dim);
614 const std::vector<Point> &,
618 libmesh_error_msg(
"Computing the divergence of a shape function only \nmakes sense for vector-valued elements.");
625 const std::vector<Point> &,
633 libmesh_error_msg(
"The divergence only make sense in 2D and 3D");
637 const std::vector<std::vector<RealGradient>> & dphidxi = fe.
get_dphidxi();
638 const std::vector<std::vector<RealGradient>> & dphideta = fe.
get_dphideta();
650 Real dphix_dx = (dphidxi[i][p].slice(0))*dxidx_map[p]
651 + (dphideta[i][p].slice(0))*detadx_map[p];
653 Real dphiy_dy = (dphidxi[i][p].slice(1))*dxidy_map[p]
654 + (dphideta[i][p].slice(1))*detady_map[p];
656 div_phi[i][p] = dphix_dx + dphiy_dy;
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();
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];
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];
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];
694 div_phi[i][p] = dphix_dx + dphiy_dy + dphiz_dz;
701 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
This is the base class from which all geometric element types are derived.
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
TensorTools::DecrementRank< OutputShape >::type OutputDivergence
const std::vector< std::vector< Real > > & get_d2etadxyz2() const
Second derivatives of "eta" reference coordinate wrt physical coordinates.
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
Second derivatives of "zeta" reference coordinate wrt physical coordinates.
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
Second derivatives of "xi" reference coordinate wrt physical coordinates.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
This class forms the foundation from which generic finite elements may be derived.
const std::vector< Real > & get_detadx() const