libMesh
Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | Friends | List of all members
libMesh::FEMap Class Reference

Class contained in FE that encapsulates mapping (i.e. More...

#include <fe_map.h>

Inheritance diagram for libMesh::FEMap:
[legend]

Public Member Functions

 FEMap (Real jtol=0)
 
virtual ~FEMap ()
 
template<unsigned int Dim>
void init_reference_to_physical_map (const std::vector< Point > &qp, const Elem *elem)
 
void compute_single_point_map (const unsigned int dim, const std::vector< Real > &qw, const Elem *elem, unsigned int p, const std::vector< const Node * > &elem_nodes, bool compute_second_derivatives)
 Compute the jacobian and some other additional data fields at the single point with index p. More...
 
virtual void compute_affine_map (const unsigned int dim, const std::vector< Real > &qw, const Elem *elem)
 Compute the jacobian and some other additional data fields. More...
 
virtual void compute_null_map (const unsigned int dim, const std::vector< Real > &qw)
 Assign a fake jacobian and some other additional data fields. More...
 
virtual void compute_map (const unsigned int dim, const std::vector< Real > &qw, const Elem *elem, bool calculate_d2phi)
 Compute the jacobian and some other additional data fields. More...
 
virtual void compute_face_map (int dim, const std::vector< Real > &qw, const Elem *side)
 Same as compute_map, but for a side. More...
 
void compute_edge_map (int dim, const std::vector< Real > &qw, const Elem *side)
 Same as before, but for an edge. More...
 
template<unsigned int Dim>
void init_face_shape_functions (const std::vector< Point > &qp, const Elem *side)
 Initializes the reference to physical element map for a side. More...
 
template<unsigned int Dim>
void init_edge_shape_functions (const std::vector< Point > &qp, const Elem *edge)
 Same as before, but for an edge. More...
 
const std::vector< Point > & get_xyz () const
 
const std::vector< Real > & get_jacobian () const
 
const std::vector< Real > & get_JxW () const
 
const std::vector< RealGradient > & get_dxyzdxi () const
 
const std::vector< RealGradient > & get_dxyzdeta () const
 
const std::vector< RealGradient > & get_dxyzdzeta () const
 
const std::vector< RealGradient > & get_d2xyzdxi2 () const
 
const std::vector< RealGradient > & get_d2xyzdeta2 () const
 
const std::vector< RealGradient > & get_d2xyzdzeta2 () const
 
const std::vector< RealGradient > & get_d2xyzdxideta () const
 
const std::vector< RealGradient > & get_d2xyzdxidzeta () const
 
const std::vector< RealGradient > & get_d2xyzdetadzeta () const
 
const std::vector< Real > & get_dxidx () const
 
const std::vector< Real > & get_dxidy () const
 
const std::vector< Real > & get_dxidz () const
 
const std::vector< Real > & get_detadx () const
 
const std::vector< Real > & get_detady () const
 
const std::vector< Real > & get_detadz () const
 
const std::vector< Real > & get_dzetadx () const
 
const std::vector< Real > & get_dzetady () const
 
const std::vector< Real > & get_dzetadz () const
 
const std::vector< std::vector< Real > > & get_d2xidxyz2 () const
 Second derivatives of "xi" reference coordinate wrt physical coordinates. More...
 
const std::vector< std::vector< Real > > & get_d2etadxyz2 () const
 Second derivatives of "eta" reference coordinate wrt physical coordinates. More...
 
const std::vector< std::vector< Real > > & get_d2zetadxyz2 () const
 Second derivatives of "zeta" reference coordinate wrt physical coordinates. More...
 
const std::vector< std::vector< Real > > & get_psi () const
 
const std::vector< std::vector< Real > > & get_phi_map () const
 
const std::vector< std::vector< Real > > & get_dphidxi_map () const
 
const std::vector< std::vector< Real > > & get_dphideta_map () const
 
const std::vector< std::vector< Real > > & get_dphidzeta_map () const
 
const std::vector< std::vector< Point > > & get_tangents () const
 
const std::vector< Point > & get_normals () const
 
const std::vector< Real > & get_curvatures () const
 
void print_JxW (std::ostream &os) const
 Prints the Jacobian times the weight for each quadrature point. More...
 
void print_xyz (std::ostream &os) const
 Prints the spatial location of each quadrature point (on the physical element). More...
 
std::vector< std::vector< Real > > & get_psi ()
 
std::vector< std::vector< Real > > & get_dpsidxi ()
 
const std::vector< std::vector< Real > > & get_dpsidxi () const
 
std::vector< std::vector< Real > > & get_dpsideta ()
 
const std::vector< std::vector< Real > > & get_dpsideta () const
 
std::vector< std::vector< Real > > & get_d2psidxi2 ()
 
const std::vector< std::vector< Real > > & get_d2psidxi2 () const
 
std::vector< std::vector< Real > > & get_d2psidxideta ()
 
const std::vector< std::vector< Real > > & get_d2psidxideta () const
 
std::vector< std::vector< Real > > & get_d2psideta2 ()
 
const std::vector< std::vector< Real > > & get_d2psideta2 () const
 
std::vector< std::vector< Real > > & get_phi_map ()
 
std::vector< std::vector< Real > > & get_dphidxi_map ()
 
std::vector< std::vector< Real > > & get_dphideta_map ()
 
std::vector< std::vector< Real > > & get_dphidzeta_map ()
 
std::vector< std::vector< Real > > & get_d2phidxi2_map ()
 
std::vector< std::vector< Real > > & get_d2phidxideta_map ()
 
std::vector< std::vector< Real > > & get_d2phidxidzeta_map ()
 
std::vector< std::vector< Real > > & get_d2phideta2_map ()
 
std::vector< std::vector< Real > > & get_d2phidetadzeta_map ()
 
std::vector< std::vector< Real > > & get_d2phidzeta2_map ()
 
std::vector< Real > & get_JxW ()
 
void set_jacobian_tolerance (Real tol)
 Set the Jacobian tolerance used for determining when the mapping fails. More...
 

Static Public Member Functions

static std::unique_ptr< FEMapbuild (FEType fe_type)
 
static FEFamily map_fe_type (const Elem &elem)
 
static Point map (const unsigned int dim, const Elem *elem, const Point &reference_point)
 
static Point map_deriv (const unsigned int dim, const Elem *elem, const unsigned int j, const Point &reference_point)
 
static Point inverse_map (const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
 
static void inverse_map (unsigned int dim, const Elem *elem, const std::vector< Point > &physical_points, std::vector< Point > &reference_points, const Real tolerance=TOLERANCE, const bool secure=true)
 Takes a number points in physical space (in the physical_points vector) and finds their location on the reference element for the input element elem. More...
 

Protected Member Functions

void determine_calculations ()
 Determine which values are to be calculated. More...
 
void resize_quadrature_map_vectors (const unsigned int dim, unsigned int n_qp)
 A utility function for use by compute_*_map. More...
 
Real dxdxi_map (const unsigned int p) const
 Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected. More...
 
Real dydxi_map (const unsigned int p) const
 Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected. More...
 
Real dzdxi_map (const unsigned int p) const
 Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected. More...
 
Real dxdeta_map (const unsigned int p) const
 Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected. More...
 
Real dydeta_map (const unsigned int p) const
 Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected. More...
 
Real dzdeta_map (const unsigned int p) const
 Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected. More...
 
Real dxdzeta_map (const unsigned int p) const
 Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected. More...
 
Real dydzeta_map (const unsigned int p) const
 Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected. More...
 
Real dzdzeta_map (const unsigned int p) const
 Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected. More...
 

Protected Attributes

std::vector< Pointxyz
 The spatial locations of the quadrature points. More...
 
std::vector< RealGradientdxyzdxi_map
 Vector of partial derivatives: d(x)/d(xi), d(y)/d(xi), d(z)/d(xi) More...
 
std::vector< RealGradientdxyzdeta_map
 Vector of partial derivatives: d(x)/d(eta), d(y)/d(eta), d(z)/d(eta) More...
 
std::vector< RealGradientdxyzdzeta_map
 Vector of partial derivatives: d(x)/d(zeta), d(y)/d(zeta), d(z)/d(zeta) More...
 
std::vector< RealGradientd2xyzdxi2_map
 Vector of second partial derivatives in xi: d^2(x)/d(xi)^2, d^2(y)/d(xi)^2, d^2(z)/d(xi)^2. More...
 
std::vector< RealGradientd2xyzdxideta_map
 Vector of mixed second partial derivatives in xi-eta: d^2(x)/d(xi)d(eta) d^2(y)/d(xi)d(eta) d^2(z)/d(xi)d(eta) More...
 
std::vector< RealGradientd2xyzdeta2_map
 Vector of second partial derivatives in eta: d^2(x)/d(eta)^2. More...
 
std::vector< RealGradientd2xyzdxidzeta_map
 Vector of second partial derivatives in xi-zeta: d^2(x)/d(xi)d(zeta), d^2(y)/d(xi)d(zeta), d^2(z)/d(xi)d(zeta) More...
 
std::vector< RealGradientd2xyzdetadzeta_map
 Vector of mixed second partial derivatives in eta-zeta: d^2(x)/d(eta)d(zeta) d^2(y)/d(eta)d(zeta) d^2(z)/d(eta)d(zeta) More...
 
std::vector< RealGradientd2xyzdzeta2_map
 Vector of second partial derivatives in zeta: d^2(x)/d(zeta)^2. More...
 
std::vector< Realdxidx_map
 Map for partial derivatives: d(xi)/d(x). More...
 
std::vector< Realdxidy_map
 Map for partial derivatives: d(xi)/d(y). More...
 
std::vector< Realdxidz_map
 Map for partial derivatives: d(xi)/d(z). More...
 
std::vector< Realdetadx_map
 Map for partial derivatives: d(eta)/d(x). More...
 
std::vector< Realdetady_map
 Map for partial derivatives: d(eta)/d(y). More...
 
std::vector< Realdetadz_map
 Map for partial derivatives: d(eta)/d(z). More...
 
std::vector< Realdzetadx_map
 Map for partial derivatives: d(zeta)/d(x). More...
 
std::vector< Realdzetady_map
 Map for partial derivatives: d(zeta)/d(y). More...
 
std::vector< Realdzetadz_map
 Map for partial derivatives: d(zeta)/d(z). More...
 
std::vector< std::vector< Real > > d2xidxyz2_map
 Second derivatives of "xi" reference coordinate wrt physical coordinates. More...
 
std::vector< std::vector< Real > > d2etadxyz2_map
 Second derivatives of "eta" reference coordinate wrt physical coordinates. More...
 
std::vector< std::vector< Real > > d2zetadxyz2_map
 Second derivatives of "zeta" reference coordinate wrt physical coordinates. More...
 
std::vector< std::vector< Real > > phi_map
 Map for the shape function phi. More...
 
std::vector< std::vector< Real > > dphidxi_map
 Map for the derivative, d(phi)/d(xi). More...
 
std::vector< std::vector< Real > > dphideta_map
 Map for the derivative, d(phi)/d(eta). More...
 
std::vector< std::vector< Real > > dphidzeta_map
 Map for the derivative, d(phi)/d(zeta). More...
 
std::vector< std::vector< Real > > d2phidxi2_map
 Map for the second derivative, d^2(phi)/d(xi)^2. More...
 
std::vector< std::vector< Real > > d2phidxideta_map
 Map for the second derivative, d^2(phi)/d(xi)d(eta). More...
 
std::vector< std::vector< Real > > d2phidxidzeta_map
 Map for the second derivative, d^2(phi)/d(xi)d(zeta). More...
 
std::vector< std::vector< Real > > d2phideta2_map
 Map for the second derivative, d^2(phi)/d(eta)^2. More...
 
std::vector< std::vector< Real > > d2phidetadzeta_map
 Map for the second derivative, d^2(phi)/d(eta)d(zeta). More...
 
std::vector< std::vector< Real > > d2phidzeta2_map
 Map for the second derivative, d^2(phi)/d(zeta)^2. More...
 
std::vector< std::vector< Real > > psi_map
 Map for the side shape functions, psi. More...
 
std::vector< std::vector< Real > > dpsidxi_map
 Map for the derivative of the side functions, d(psi)/d(xi). More...
 
std::vector< std::vector< Real > > dpsideta_map
 Map for the derivative of the side function, d(psi)/d(eta). More...
 
std::vector< std::vector< Real > > d2psidxi2_map
 Map for the second derivatives (in xi) of the side shape functions. More...
 
std::vector< std::vector< Real > > d2psidxideta_map
 Map for the second (cross) derivatives in xi, eta of the side shape functions. More...
 
std::vector< std::vector< Real > > d2psideta2_map
 Map for the second derivatives (in eta) of the side shape functions. More...
 
std::vector< std::vector< Point > > tangents
 Tangent vectors on boundary at quadrature points. More...
 
std::vector< Pointnormals
 Normal vectors on boundary at quadrature points. More...
 
std::vector< Realcurvatures
 The mean curvature (= one half the sum of the principal curvatures) on the boundary at the quadrature points. More...
 
std::vector< Realjac
 Jacobian values at quadrature points. More...
 
std::vector< RealJxW
 Jacobian*Weight values at quadrature points. More...
 
bool calculations_started
 Have calculations with this object already been started? Then all get_* functions should already have been called. More...
 
bool calculate_xyz
 Should we calculate physical point locations? More...
 
bool calculate_dxyz
 Should we calculate mapping gradients? More...
 
bool calculate_d2xyz
 Should we calculate mapping hessians? More...
 
Real jacobian_tolerance
 The Jacobian tolerance used for determining when the mapping fails. More...
 

Private Member Functions

void compute_inverse_map_second_derivs (unsigned p)
 A helper function used by FEMap::compute_single_point_map() to compute second derivatives of the inverse map. More...
 

Private Attributes

std::vector< const Node * > _elem_nodes
 Work vector for compute_affine_map() More...
 

Friends

template<unsigned int Dim, FEFamily T>
class FE
 FE classes should be able to reset calculations_started in a few special cases. More...
 

Detailed Description

Class contained in FE that encapsulates mapping (i.e.

from physical space to reference space and vice-versa) quantities and operations.

Author
Paul Bauman
Date
2012

Computes finite element mapping function values, gradients, etc.

Definition at line 47 of file fe_map.h.

Constructor & Destructor Documentation

◆ FEMap()

libMesh::FEMap::FEMap ( Real  jtol = 0)

Definition at line 64 of file fe_map.C.

64  :
65  calculations_started(false),
66  calculate_xyz(false),
67  calculate_dxyz(false),
68 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
69  calculate_d2xyz(false),
70 #endif
71  jacobian_tolerance(jtol)
72 {}

◆ ~FEMap()

virtual libMesh::FEMap::~FEMap ( )
inlinevirtual

Definition at line 52 of file fe_map.h.

52 {}

Member Function Documentation

◆ build()

std::unique_ptr< FEMap > libMesh::FEMap::build ( FEType  fe_type)
static

Definition at line 76 of file fe_map.C.

77 {
78  switch( fe_type.family )
79  {
80  case XYZ:
81  return libmesh_make_unique<FEXYZMap>();
82 
83  default:
84  return libmesh_make_unique<FEMap>();
85  }
86 }

References libMesh::FEType::family, and libMesh::XYZ.

◆ compute_affine_map()

void libMesh::FEMap::compute_affine_map ( const unsigned int  dim,
const std::vector< Real > &  qw,
const Elem elem 
)
virtual

Compute the jacobian and some other additional data fields.

Takes the integration weights as input, along with a pointer to the element. The element is assumed to have a constant Jacobian

Definition at line 1277 of file fe_map.C.

1280 {
1281  // Start logging the map computation.
1282  LOG_SCOPE("compute_affine_map()", "FEMap");
1283 
1284  libmesh_assert(elem);
1285 
1286  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1287 
1288  // Resize the vectors to hold data at the quadrature points
1289  this->resize_quadrature_map_vectors(dim, n_qp);
1290 
1291  // Determine the nodes contributing to element elem
1292  unsigned int n_nodes = elem->n_nodes();
1293  _elem_nodes.resize(elem->n_nodes());
1294  for (unsigned int i=0; i<n_nodes; i++)
1295  _elem_nodes[i] = elem->node_ptr(i);
1296 
1297  // Compute map at quadrature point 0
1298  this->compute_single_point_map(dim, qw, elem, 0, _elem_nodes, /*compute_second_derivatives=*/false);
1299 
1300  // Compute xyz at all other quadrature points
1301  if (calculate_xyz)
1302  for (unsigned int p=1; p<n_qp; p++)
1303  {
1304  xyz[p].zero();
1305  for (auto i : index_range(phi_map)) // sum over the nodes
1306  xyz[p].add_scaled (*_elem_nodes[i], phi_map[i][p]);
1307  }
1308 
1309  // Copy other map data from quadrature point 0
1310  if (calculate_dxyz)
1311  for (unsigned int p=1; p<n_qp; p++) // for each extra quadrature point
1312  {
1313  dxyzdxi_map[p] = dxyzdxi_map[0];
1314  dxidx_map[p] = dxidx_map[0];
1315  dxidy_map[p] = dxidy_map[0];
1316  dxidz_map[p] = dxidz_map[0];
1317 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1318  // The map should be affine, so second derivatives are zero
1319  if (calculate_d2xyz)
1320  d2xyzdxi2_map[p] = 0.;
1321 #endif
1322  if (dim > 1)
1323  {
1324  dxyzdeta_map[p] = dxyzdeta_map[0];
1325  detadx_map[p] = detadx_map[0];
1326  detady_map[p] = detady_map[0];
1327  detadz_map[p] = detadz_map[0];
1328 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1329  if (calculate_d2xyz)
1330  {
1331  d2xyzdxideta_map[p] = 0.;
1332  d2xyzdeta2_map[p] = 0.;
1333  }
1334 #endif
1335  if (dim > 2)
1336  {
1337  dxyzdzeta_map[p] = dxyzdzeta_map[0];
1338  dzetadx_map[p] = dzetadx_map[0];
1339  dzetady_map[p] = dzetady_map[0];
1340  dzetadz_map[p] = dzetadz_map[0];
1341 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1342  if (calculate_d2xyz)
1343  {
1344  d2xyzdxidzeta_map[p] = 0.;
1345  d2xyzdetadzeta_map[p] = 0.;
1346  d2xyzdzeta2_map[p] = 0.;
1347  }
1348 #endif
1349  }
1350  }
1351  jac[p] = jac[0];
1352  JxW[p] = JxW[0] / qw[0] * qw[p];
1353  }
1354 }

References _elem_nodes, calculate_d2xyz, calculate_dxyz, calculate_xyz, compute_single_point_map(), d2xyzdeta2_map, d2xyzdetadzeta_map, d2xyzdxi2_map, d2xyzdxideta_map, d2xyzdxidzeta_map, d2xyzdzeta2_map, detadx_map, detady_map, detadz_map, dim, dxidx_map, dxidy_map, dxidz_map, dxyzdeta_map, dxyzdxi_map, dxyzdzeta_map, dzetadx_map, dzetady_map, dzetadz_map, libMesh::index_range(), jac, JxW, libMesh::libmesh_assert(), n_nodes, libMesh::Elem::n_nodes(), libMesh::Elem::node_ptr(), phi_map, resize_quadrature_map_vectors(), and xyz.

Referenced by compute_map().

◆ compute_edge_map()

void libMesh::FEMap::compute_edge_map ( int  dim,
const std::vector< Real > &  qw,
const Elem side 
)

Same as before, but for an edge.

Useful for some projections.

Definition at line 950 of file fe_boundary.C.

953 {
954  libmesh_assert(edge);
955 
956  if (dim == 2)
957  {
958  // A 2D finite element living in either 2D or 3D space.
959  // The edges here are the sides of the element, so the
960  // (misnamed) compute_face_map function does what we want
961  this->compute_face_map(dim, qw, edge);
962  return;
963  }
964 
965  libmesh_assert_equal_to (dim, 3); // 1D is unnecessary and currently unsupported
966 
967  LOG_SCOPE("compute_edge_map()", "FEMap");
968 
969  // We're calculating now!
970  this->determine_calculations();
971 
972  // The number of quadrature points.
973  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
974 
975  // Resize the vectors to hold data at the quadrature points
976  if (calculate_xyz)
977  this->xyz.resize(n_qp);
978  if (calculate_dxyz)
979  {
980  this->dxyzdxi_map.resize(n_qp);
981  this->dxyzdeta_map.resize(n_qp);
982  this->tangents.resize(n_qp);
983  this->normals.resize(n_qp);
984  this->JxW.resize(n_qp);
985  }
986 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
987  if (calculate_d2xyz)
988  {
989  this->d2xyzdxi2_map.resize(n_qp);
990  this->d2xyzdxideta_map.resize(n_qp);
991  this->d2xyzdeta2_map.resize(n_qp);
992  this->curvatures.resize(n_qp);
993  }
994 #endif
995 
996  // Clear the entities that will be summed
997  for (unsigned int p=0; p<n_qp; p++)
998  {
999  if (calculate_xyz)
1000  this->xyz[p].zero();
1001  if (calculate_dxyz)
1002  {
1003  this->tangents[p].resize(1);
1004  this->dxyzdxi_map[p].zero();
1005  this->dxyzdeta_map[p].zero();
1006  }
1007 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1008  if (calculate_d2xyz)
1009  {
1010  this->d2xyzdxi2_map[p].zero();
1011  this->d2xyzdxideta_map[p].zero();
1012  this->d2xyzdeta2_map[p].zero();
1013  }
1014 #endif
1015  }
1016 
1017  // compute x, dxdxi at the quadrature points
1018  for (unsigned int i=0,
1019  psi_map_size=cast_int<unsigned int>(psi_map.size());
1020  i != psi_map_size; i++) // sum over the nodes
1021  {
1022  const Point & edge_point = edge->point(i);
1023 
1024  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
1025  {
1026  if (calculate_xyz)
1027  this->xyz[p].add_scaled (edge_point, this->psi_map[i][p]);
1028  if (calculate_dxyz)
1029  this->dxyzdxi_map[p].add_scaled (edge_point, this->dpsidxi_map[i][p]);
1030 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1031  if (calculate_d2xyz)
1032  this->d2xyzdxi2_map[p].add_scaled (edge_point, this->d2psidxi2_map[i][p]);
1033 #endif
1034  }
1035  }
1036 
1037  // Compute the tangents at the quadrature point
1038  // FIXME: normals (plural!) and curvatures are uncalculated
1039  if (calculate_dxyz)
1040  for (unsigned int p=0; p<n_qp; p++)
1041  {
1042  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
1043  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
1044 
1045  // compute the jacobian at the quadrature points
1046  const Real the_jac = std::sqrt(this->dxdxi_map(p)*this->dxdxi_map(p) +
1047  this->dydxi_map(p)*this->dydxi_map(p) +
1048  this->dzdxi_map(p)*this->dzdxi_map(p));
1049 
1050  libmesh_assert_greater (the_jac, 0.);
1051 
1052  this->JxW[p] = the_jac*qw[p];
1053  }
1054 }

References calculate_d2xyz, calculate_dxyz, calculate_xyz, compute_face_map(), curvatures, d2psidxi2_map, d2xyzdeta2_map, d2xyzdxi2_map, d2xyzdxideta_map, determine_calculations(), dim, dpsidxi_map, dxdxi_map(), dxyzdeta_map, dxyzdxi_map, dydxi_map(), dzdxi_map(), JxW, libMesh::libmesh_assert(), normals, libMesh::Elem::point(), psi_map, libMesh::Real, std::sqrt(), tangents, and xyz.

◆ compute_face_map()

void libMesh::FEMap::compute_face_map ( int  dim,
const std::vector< Real > &  qw,
const Elem side 
)
virtual

Same as compute_map, but for a side.

Useful for boundary integration.

Reimplemented in libMesh::FEXYZMap.

Definition at line 587 of file fe_boundary.C.

589 {
590  libmesh_assert(side);
591 
592  // We're calculating now!
593  this->determine_calculations();
594 
595  LOG_SCOPE("compute_face_map()", "FEMap");
596 
597  // The number of quadrature points.
598  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
599 
600  const FEFamily mapping_family = FEMap::map_fe_type(*side);
601  const Order mapping_order (side->default_order());
602  const FEType map_fe_type(mapping_order, mapping_family);
603  const ElemType mapping_elem_type (side->type());
604  const unsigned int n_mapping_shape_functions =
605  FEInterface::n_shape_functions(dim, map_fe_type, mapping_elem_type);
606 
607  switch (dim)
608  {
609  case 1:
610  {
611  // A 1D finite element, currently assumed to be in 1D space
612  // This means the boundary is a "0D finite element", a
613  // NODEELEM.
614 
615  // Resize the vectors to hold data at the quadrature points
616  {
617  if (calculate_xyz)
618  this->xyz.resize(n_qp);
619  if (calculate_dxyz)
620  normals.resize(n_qp);
621 
622  if (calculate_dxyz)
623  this->JxW.resize(n_qp);
624  }
625 
626  // If we have no quadrature points, there's nothing else to do
627  if (!n_qp)
628  break;
629 
630  // We need to look back at the full edge to figure out the normal
631  // vector
632  const Elem * elem = side->parent();
633  libmesh_assert (elem);
634  if (calculate_dxyz)
635  {
636  if (side->node_id(0) == elem->node_id(0))
637  normals[0] = Point(-1.);
638  else
639  {
640  libmesh_assert_equal_to (side->node_id(0),
641  elem->node_id(1));
642  normals[0] = Point(1.);
643  }
644  }
645 
646  // Calculate x at the point
647  if (calculate_xyz)
648  libmesh_assert_equal_to (this->psi_map.size(), 1);
649  // In the unlikely event we have multiple quadrature
650  // points, they'll be in the same place
651  for (unsigned int p=0; p<n_qp; p++)
652  {
653  if (calculate_xyz)
654  {
655  this->xyz[p].zero();
656  this->xyz[p].add_scaled (side->point(0), this->psi_map[0][p]);
657  }
658  if (calculate_dxyz)
659  {
660  normals[p] = normals[0];
661  this->JxW[p] = 1.0*qw[p];
662  }
663  }
664 
665  // done computing the map
666  break;
667  }
668 
669  case 2:
670  {
671  // A 2D finite element living in either 2D or 3D space.
672  // This means the boundary is a 1D finite element, i.e.
673  // and EDGE2 or EDGE3.
674  // Resize the vectors to hold data at the quadrature points
675  {
676  if (calculate_xyz)
677  this->xyz.resize(n_qp);
678  if (calculate_dxyz)
679  {
680  this->dxyzdxi_map.resize(n_qp);
681  this->tangents.resize(n_qp);
682  this->normals.resize(n_qp);
683 
684  this->JxW.resize(n_qp);
685  }
686 
687 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
688  if (calculate_d2xyz)
689  {
690  this->d2xyzdxi2_map.resize(n_qp);
691  this->curvatures.resize(n_qp);
692  }
693 #endif
694  }
695 
696  // Clear the entities that will be summed
697  // Compute the tangent & normal at the quadrature point
698  for (unsigned int p=0; p<n_qp; p++)
699  {
700  if (calculate_xyz)
701  this->xyz[p].zero();
702  if (calculate_dxyz)
703  {
704  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
705  this->dxyzdxi_map[p].zero();
706  }
707 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
708  if (calculate_d2xyz)
709  this->d2xyzdxi2_map[p].zero();
710 #endif
711  }
712 
713  // compute x, dxdxi at the quadrature points
714  for (unsigned int i=0; i<n_mapping_shape_functions; i++) // sum over the nodes
715  {
716  const Point & side_point = side->point(i);
717 
718  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
719  {
720  if (calculate_xyz)
721  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
722  if (calculate_dxyz)
723  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
724 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
725  if (calculate_d2xyz)
726  this->d2xyzdxi2_map[p].add_scaled(side_point, this->d2psidxi2_map[i][p]);
727 #endif
728  }
729  }
730 
731  // Compute the tangent & normal at the quadrature point
732  if (calculate_dxyz)
733  {
734  for (unsigned int p=0; p<n_qp; p++)
735  {
736  // The first tangent comes from just the edge's Jacobian
737  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
738 
739 #if LIBMESH_DIM == 2
740  // For a 2D element living in 2D, the normal is given directly
741  // from the entries in the edge Jacobian.
742  this->normals[p] = (Point(this->dxyzdxi_map[p](1), -this->dxyzdxi_map[p](0), 0.)).unit();
743 
744 #elif LIBMESH_DIM == 3
745  // For a 2D element living in 3D, there is a second tangent.
746  // For the second tangent, we need to refer to the full
747  // element's (not just the edge's) Jacobian.
748  const Elem * elem = side->parent();
749  libmesh_assert(elem);
750 
751  // Inverse map xyz[p] to a reference point on the parent...
752  Point reference_point = FEMap::inverse_map(2, elem, this->xyz[p]);
753 
754  // Get dxyz/dxi and dxyz/deta from the parent map.
755  Point dx_dxi = FEMap::map_deriv (2, elem, 0, reference_point);
756  Point dx_deta = FEMap::map_deriv (2, elem, 1, reference_point);
757 
758  // The second tangent vector is formed by crossing these vectors.
759  tangents[p][1] = dx_dxi.cross(dx_deta).unit();
760 
761  // Finally, the normal in this case is given by crossing these
762  // two tangents.
763  normals[p] = tangents[p][0].cross(tangents[p][1]).unit();
764 #endif
765 
766 
767 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
768  // The curvature is computed via the familiar Frenet formula:
769  // curvature = [d^2(x) / d (xi)^2] dot [normal]
770  // For a reference, see:
771  // F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill, p. 310
772  //
773  // Note: The sign convention here is different from the
774  // 3D case. Concave-upward curves (smiles) have a positive
775  // curvature. Concave-downward curves (frowns) have a
776  // negative curvature. Be sure to take that into account!
777  if (calculate_d2xyz)
778  {
779  const Real numerator = this->d2xyzdxi2_map[p] * this->normals[p];
780  const Real denominator = this->dxyzdxi_map[p].norm_sq();
781  libmesh_assert_not_equal_to (denominator, 0);
782  curvatures[p] = numerator / denominator;
783  }
784 #endif
785  }
786 
787  // compute the jacobian at the quadrature points
788  for (unsigned int p=0; p<n_qp; p++)
789  {
790  const Real the_jac = this->dxyzdxi_map[p].norm();
791 
792  libmesh_assert_greater (the_jac, 0.);
793 
794  this->JxW[p] = the_jac*qw[p];
795  }
796  }
797 
798  // done computing the map
799  break;
800  }
801 
802 
803 
804  case 3:
805  {
806  // A 3D finite element living in 3D space.
807  // Resize the vectors to hold data at the quadrature points
808  {
809  if (calculate_xyz)
810  this->xyz.resize(n_qp);
811  if (calculate_dxyz)
812  {
813  this->dxyzdxi_map.resize(n_qp);
814  this->dxyzdeta_map.resize(n_qp);
815  this->tangents.resize(n_qp);
816  this->normals.resize(n_qp);
817  this->JxW.resize(n_qp);
818  }
819 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
820  if (calculate_d2xyz)
821  {
822  this->d2xyzdxi2_map.resize(n_qp);
823  this->d2xyzdxideta_map.resize(n_qp);
824  this->d2xyzdeta2_map.resize(n_qp);
825  this->curvatures.resize(n_qp);
826  }
827 #endif
828  }
829 
830  // Clear the entities that will be summed
831  for (unsigned int p=0; p<n_qp; p++)
832  {
833  if (calculate_xyz)
834  this->xyz[p].zero();
835  if (calculate_dxyz)
836  {
837  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
838  this->dxyzdxi_map[p].zero();
839  this->dxyzdeta_map[p].zero();
840  }
841 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
842  if (calculate_d2xyz)
843  {
844  this->d2xyzdxi2_map[p].zero();
845  this->d2xyzdxideta_map[p].zero();
846  this->d2xyzdeta2_map[p].zero();
847  }
848 #endif
849  }
850 
851  // compute x, dxdxi at the quadrature points
852  for (unsigned int i=0; i<n_mapping_shape_functions; i++) // sum over the nodes
853  {
854  const Point & side_point = side->point(i);
855 
856  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
857  {
858  if (calculate_xyz)
859  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
860  if (calculate_dxyz)
861  {
862  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
863  this->dxyzdeta_map[p].add_scaled(side_point, this->dpsideta_map[i][p]);
864  }
865 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
866  if (calculate_d2xyz)
867  {
868  this->d2xyzdxi2_map[p].add_scaled (side_point, this->d2psidxi2_map[i][p]);
869  this->d2xyzdxideta_map[p].add_scaled(side_point, this->d2psidxideta_map[i][p]);
870  this->d2xyzdeta2_map[p].add_scaled (side_point, this->d2psideta2_map[i][p]);
871  }
872 #endif
873  }
874  }
875 
876  // Compute the tangents, normal, and curvature at the quadrature point
877  if (calculate_dxyz)
878  {
879  for (unsigned int p=0; p<n_qp; p++)
880  {
881  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
882  this->normals[p] = n.unit();
883  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
884  this->tangents[p][1] = n.cross(this->dxyzdxi_map[p]).unit();
885 
886 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
887  if (calculate_d2xyz)
888  {
889  // Compute curvature using the typical nomenclature
890  // of the first and second fundamental forms.
891  // For reference, see:
892  // 1) http://mathworld.wolfram.com/MeanCurvature.html
893  // (note -- they are using inward normal)
894  // 2) F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill
895  const Real L = -this->d2xyzdxi2_map[p] * this->normals[p];
896  const Real M = -this->d2xyzdxideta_map[p] * this->normals[p];
897  const Real N = -this->d2xyzdeta2_map[p] * this->normals[p];
898  const Real E = this->dxyzdxi_map[p].norm_sq();
899  const Real F = this->dxyzdxi_map[p] * this->dxyzdeta_map[p];
900  const Real G = this->dxyzdeta_map[p].norm_sq();
901 
902  const Real numerator = E*N -2.*F*M + G*L;
903  const Real denominator = E*G - F*F;
904  libmesh_assert_not_equal_to (denominator, 0.);
905  curvatures[p] = 0.5*numerator/denominator;
906  }
907 #endif
908  }
909 
910  // compute the jacobian at the quadrature points, see
911  // http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/thtoc.htm
912  for (unsigned int p=0; p<n_qp; p++)
913  {
914  const Real g11 = (dxdxi_map(p)*dxdxi_map(p) +
915  dydxi_map(p)*dydxi_map(p) +
916  dzdxi_map(p)*dzdxi_map(p));
917 
918  const Real g12 = (dxdxi_map(p)*dxdeta_map(p) +
919  dydxi_map(p)*dydeta_map(p) +
920  dzdxi_map(p)*dzdeta_map(p));
921 
922  const Real g21 = g12;
923 
924  const Real g22 = (dxdeta_map(p)*dxdeta_map(p) +
925  dydeta_map(p)*dydeta_map(p) +
926  dzdeta_map(p)*dzdeta_map(p));
927 
928 
929  const Real the_jac = std::sqrt(g11*g22 - g12*g21);
930 
931  libmesh_assert_greater (the_jac, 0.);
932 
933  this->JxW[p] = the_jac*qw[p];
934  }
935  }
936 
937  // done computing the map
938  break;
939  }
940 
941 
942  default:
943  libmesh_error_msg("Invalid dimension dim = " << dim);
944  }
945 }

References calculate_d2xyz, calculate_dxyz, calculate_xyz, libMesh::TypeVector< T >::cross(), curvatures, d2psideta2_map, d2psidxi2_map, d2psidxideta_map, d2xyzdeta2_map, d2xyzdxi2_map, d2xyzdxideta_map, libMesh::Elem::default_order(), determine_calculations(), dim, dpsideta_map, dpsidxi_map, dxdeta_map(), dxdxi_map(), dxyzdeta_map, dxyzdxi_map, dydeta_map(), dydxi_map(), dzdeta_map(), dzdxi_map(), inverse_map(), JxW, libMesh::libmesh_assert(), map_deriv(), map_fe_type(), libMesh::FEInterface::n_shape_functions(), libMesh::Elem::node_id(), normals, libMesh::Elem::parent(), libMesh::Elem::point(), psi_map, libMesh::Real, std::sqrt(), tangents, libMesh::Elem::type(), libMesh::TypeVector< T >::unit(), and xyz.

Referenced by compute_edge_map().

◆ compute_inverse_map_second_derivs()

void libMesh::FEMap::compute_inverse_map_second_derivs ( unsigned  p)
private

A helper function used by FEMap::compute_single_point_map() to compute second derivatives of the inverse map.

Definition at line 1502 of file fe_map.C.

1503 {
1504 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1505  // Only certain second derivatives are valid depending on the
1506  // dimension...
1507  std::set<unsigned> valid_indices;
1508 
1509  // Construct J^{-1}, A, and B matrices (see JWP's notes for details)
1510  // for cases in which the element dimension matches LIBMESH_DIM.
1511 #if LIBMESH_DIM==1
1512  RealTensor
1513  Jinv(dxidx_map[p], 0., 0.,
1514  0., 0., 0.,
1515  0., 0., 0.),
1516 
1517  A(d2xyzdxi2_map[p](0), 0., 0.,
1518  0., 0., 0.,
1519  0., 0., 0.),
1520 
1521  B(0., 0., 0.,
1522  0., 0., 0.,
1523  0., 0., 0.);
1524 
1526  dxi (dxidx_map[p], 0., 0.),
1527  deta (0., 0., 0.),
1528  dzeta(0., 0., 0.);
1529 
1530  // In 1D, we have only the xx second derivative
1531  valid_indices.insert(0);
1532 
1533 #elif LIBMESH_DIM==2
1534  RealTensor
1535  Jinv(dxidx_map[p], dxidy_map[p], 0.,
1536  detadx_map[p], detady_map[p], 0.,
1537  0., 0., 0.),
1538 
1539  A(d2xyzdxi2_map[p](0), d2xyzdeta2_map[p](0), 0.,
1540  d2xyzdxi2_map[p](1), d2xyzdeta2_map[p](1), 0.,
1541  0., 0., 0.),
1542 
1543  B(d2xyzdxideta_map[p](0), 0., 0.,
1544  d2xyzdxideta_map[p](1), 0., 0.,
1545  0., 0., 0.);
1546 
1548  dxi (dxidx_map[p], dxidy_map[p], 0.),
1549  deta (detadx_map[p], detady_map[p], 0.),
1550  dzeta(0., 0., 0.);
1551 
1552  // In 2D, we have xx, xy, and yy second derivatives
1553  const unsigned tmp[3] = {0,1,3};
1554  valid_indices.insert(tmp, tmp+3);
1555 
1556 #elif LIBMESH_DIM==3
1557  RealTensor
1558  Jinv(dxidx_map[p], dxidy_map[p], dxidz_map[p],
1559  detadx_map[p], detady_map[p], detadz_map[p],
1560  dzetadx_map[p], dzetady_map[p], dzetadz_map[p]),
1561 
1562  A(d2xyzdxi2_map[p](0), d2xyzdeta2_map[p](0), d2xyzdzeta2_map[p](0),
1563  d2xyzdxi2_map[p](1), d2xyzdeta2_map[p](1), d2xyzdzeta2_map[p](1),
1564  d2xyzdxi2_map[p](2), d2xyzdeta2_map[p](2), d2xyzdzeta2_map[p](2)),
1565 
1569 
1571  dxi (dxidx_map[p], dxidy_map[p], dxidz_map[p]),
1572  deta (detadx_map[p], detady_map[p], detadz_map[p]),
1573  dzeta(dzetadx_map[p], dzetady_map[p], dzetadz_map[p]);
1574 
1575  // In 3D, we have xx, xy, xz, yy, yz, and zz second derivatives
1576  const unsigned tmp[6] = {0,1,2,3,4,5};
1577  valid_indices.insert(tmp, tmp+6);
1578 
1579 #endif
1580 
1581  // For (s,t) in {(x,x), (x,y), (x,z), (y,y), (y,z), (z,z)}, compute the
1582  // vector of inverse map second derivatives [xi_{s t}, eta_{s t}, zeta_{s t}]
1583  unsigned ctr=0;
1584  for (unsigned s=0; s<3; ++s)
1585  for (unsigned t=s; t<3; ++t)
1586  {
1587  if (valid_indices.count(ctr))
1588  {
1590  v1(dxi(s)*dxi(t),
1591  deta(s)*deta(t),
1592  dzeta(s)*dzeta(t)),
1593 
1594  v2(dxi(s)*deta(t) + deta(s)*dxi(t),
1595  dxi(s)*dzeta(t) + dzeta(s)*dxi(t),
1596  deta(s)*dzeta(t) + dzeta(s)*deta(t));
1597 
1598  // Compute the inverse map second derivatives
1599  RealVectorValue v3 = -Jinv*(A*v1 + B*v2);
1600 
1601  // Store them in the appropriate locations in the class data structures
1602  d2xidxyz2_map[p][ctr] = v3(0);
1603 
1604  if (LIBMESH_DIM > 1)
1605  d2etadxyz2_map[p][ctr] = v3(1);
1606 
1607  if (LIBMESH_DIM > 2)
1608  d2zetadxyz2_map[p][ctr] = v3(2);
1609  }
1610 
1611  // Increment the counter
1612  ctr++;
1613  }
1614 #else
1615  // to avoid compiler warnings:
1616  libmesh_ignore(p);
1617 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1618 }

References A, d2etadxyz2_map, d2xidxyz2_map, d2xyzdeta2_map, d2xyzdetadzeta_map, d2xyzdxi2_map, d2xyzdxideta_map, d2xyzdxidzeta_map, d2xyzdzeta2_map, d2zetadxyz2_map, detadx_map, detady_map, detadz_map, dxidx_map, dxidy_map, dxidz_map, dzetadx_map, dzetady_map, dzetadz_map, and libMesh::libmesh_ignore().

Referenced by compute_single_point_map().

◆ compute_map()

void libMesh::FEMap::compute_map ( const unsigned int  dim,
const std::vector< Real > &  qw,
const Elem elem,
bool  calculate_d2phi 
)
virtual

Compute the jacobian and some other additional data fields.

Takes the integration weights as input, along with a pointer to the element. Also takes a boolean parameter indicating whether second derivatives need to be calculated, allowing us to potentially skip unnecessary, expensive computations.

Definition at line 1433 of file fe_map.C.

1437 {
1438  if (!elem)
1439  {
1440  compute_null_map(dim, qw);
1441  return;
1442  }
1443 
1444  if (elem->has_affine_map())
1445  {
1446  compute_affine_map(dim, qw, elem);
1447  return;
1448  }
1449 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
1450  libmesh_assert(!calculate_d2phi);
1451 #endif
1452 
1453  // Start logging the map computation.
1454  LOG_SCOPE("compute_map()", "FEMap");
1455 
1456  libmesh_assert(elem);
1457 
1458  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1459 
1460  // Resize the vectors to hold data at the quadrature points
1461  this->resize_quadrature_map_vectors(dim, n_qp);
1462 
1463  // Determine the nodes contributing to element elem
1464  if (elem->type() == TRI3SUBDIVISION)
1465  {
1466  // Subdivision surface FE require the 1-ring around elem
1467  libmesh_assert_equal_to (dim, 2);
1468  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
1470  }
1471  else
1472  {
1473  // All other FE use only the nodes of elem itself
1474  _elem_nodes.resize(elem->n_nodes(), nullptr);
1475  for (auto i : elem->node_index_range())
1476  _elem_nodes[i] = elem->node_ptr(i);
1477  }
1478 
1479  // Compute map at all quadrature points
1480  for (unsigned int p=0; p!=n_qp; p++)
1481  this->compute_single_point_map(dim, qw, elem, p, _elem_nodes, calculate_d2phi);
1482 }

References _elem_nodes, compute_affine_map(), compute_null_map(), compute_single_point_map(), dim, libMesh::MeshTools::Subdivision::find_one_ring(), libMesh::Elem::has_affine_map(), libMesh::libmesh_assert(), libMesh::Elem::n_nodes(), libMesh::Elem::node_index_range(), libMesh::Elem::node_ptr(), resize_quadrature_map_vectors(), libMesh::TRI3SUBDIVISION, and libMesh::Elem::type().

◆ compute_null_map()

void libMesh::FEMap::compute_null_map ( const unsigned int  dim,
const std::vector< Real > &  qw 
)
virtual

Assign a fake jacobian and some other additional data fields.

Takes the integration weights as input. For use on non-element evaluations.

Definition at line 1358 of file fe_map.C.

1360 {
1361  // Start logging the map computation.
1362  LOG_SCOPE("compute_null_map()", "FEMap");
1363 
1364  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1365 
1366  // Resize the vectors to hold data at the quadrature points
1367  this->resize_quadrature_map_vectors(dim, n_qp);
1368 
1369  // Compute "fake" xyz
1370  for (unsigned int p=1; p<n_qp; p++)
1371  {
1372  if (calculate_xyz)
1373  xyz[p].zero();
1374 
1375  if (calculate_dxyz)
1376  {
1377  dxyzdxi_map[p] = 0;
1378  dxidx_map[p] = 0;
1379  dxidy_map[p] = 0;
1380  dxidz_map[p] = 0;
1381  }
1382 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1383  if (calculate_d2xyz)
1384  {
1385  d2xyzdxi2_map[p] = 0;
1386  }
1387 #endif
1388  if (dim > 1)
1389  {
1390  if (calculate_dxyz)
1391  {
1392  dxyzdeta_map[p] = 0;
1393  detadx_map[p] = 0;
1394  detady_map[p] = 0;
1395  detadz_map[p] = 0;
1396  }
1397 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1398  if (calculate_d2xyz)
1399  {
1400  d2xyzdxideta_map[p] = 0.;
1401  d2xyzdeta2_map[p] = 0.;
1402  }
1403 #endif
1404  if (dim > 2)
1405  {
1406  if (calculate_dxyz)
1407  {
1408  dxyzdzeta_map[p] = 0;
1409  dzetadx_map[p] = 0;
1410  dzetady_map[p] = 0;
1411  dzetadz_map[p] = 0;
1412  }
1413 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1414  if (calculate_d2xyz)
1415  {
1416  d2xyzdxidzeta_map[p] = 0;
1417  d2xyzdetadzeta_map[p] = 0;
1418  d2xyzdzeta2_map[p] = 0;
1419  }
1420 #endif
1421  }
1422  }
1423  if (calculate_dxyz)
1424  {
1425  jac[p] = 1;
1426  JxW[p] = qw[p];
1427  }
1428  }
1429 }

References calculate_d2xyz, calculate_dxyz, calculate_xyz, d2xyzdeta2_map, d2xyzdetadzeta_map, d2xyzdxi2_map, d2xyzdxideta_map, d2xyzdxidzeta_map, d2xyzdzeta2_map, detadx_map, detady_map, detadz_map, dim, dxidx_map, dxidy_map, dxidz_map, dxyzdeta_map, dxyzdxi_map, dxyzdzeta_map, dzetadx_map, dzetady_map, dzetadz_map, jac, JxW, resize_quadrature_map_vectors(), and xyz.

Referenced by compute_map().

◆ compute_single_point_map()

void libMesh::FEMap::compute_single_point_map ( const unsigned int  dim,
const std::vector< Real > &  qw,
const Elem elem,
unsigned int  p,
const std::vector< const Node * > &  elem_nodes,
bool  compute_second_derivatives 
)

Compute the jacobian and some other additional data fields at the single point with index p.

Takes the integration weights as input, along with a pointer to the element and a list of points that contribute to the element. Also takes a boolean flag telling whether second derivatives should actually be computed.

Definition at line 446 of file fe_map.C.

452 {
453  libmesh_assert(elem);
455 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
456  libmesh_assert(!compute_second_derivatives);
457 #endif
458 
459  if (calculate_xyz)
460  libmesh_assert_equal_to(phi_map.size(), elem_nodes.size());
461 
462  switch (dim)
463  {
464  //--------------------------------------------------------------------
465  // 0D
466  case 0:
467  {
468  libmesh_assert(elem_nodes[0]);
469  if (calculate_xyz)
470  xyz[p] = *elem_nodes[0];
471  if (calculate_dxyz)
472  {
473  jac[p] = 1.0;
474  JxW[p] = qw[p];
475  }
476  break;
477  }
478 
479  //--------------------------------------------------------------------
480  // 1D
481  case 1:
482  {
483  // Clear the entities that will be summed
484  if (calculate_xyz)
485  xyz[p].zero();
486  if (calculate_dxyz)
487  dxyzdxi_map[p].zero();
488 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
489  if (calculate_d2xyz)
490  {
491  d2xyzdxi2_map[p].zero();
492  // Inverse map second derivatives
493  d2xidxyz2_map[p].assign(6, 0.);
494  }
495 #endif
496 
497  // compute x, dx, d2x at the quadrature point
498  for (auto i : index_range(elem_nodes)) // sum over the nodes
499  {
500  // Reference to the point, helps eliminate
501  // excessive temporaries in the inner loop
502  libmesh_assert(elem_nodes[i]);
503  const Point & elem_point = *elem_nodes[i];
504 
505  if (calculate_xyz)
506  xyz[p].add_scaled (elem_point, phi_map[i][p] );
507  if (calculate_dxyz)
508  dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p]);
509 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
510  if (calculate_d2xyz)
511  d2xyzdxi2_map[p].add_scaled(elem_point, d2phidxi2_map[i][p]);
512 #endif
513  }
514 
515  // Compute the jacobian
516  //
517  // 1D elements can live in 2D or 3D space.
518  // The transformation matrix from local->global
519  // coordinates is
520  //
521  // T = | dx/dxi |
522  // | dy/dxi |
523  // | dz/dxi |
524  //
525  // The generalized determinant of T (from the
526  // so-called "normal" eqns.) is
527  // jac = "det(T)" = sqrt(det(T'T))
528  //
529  // where T'= transpose of T, so
530  //
531  // jac = sqrt( (dx/dxi)^2 + (dy/dxi)^2 + (dz/dxi)^2 )
532 
533  if (calculate_dxyz)
534  {
535  jac[p] = dxyzdxi_map[p].norm();
536 
537  if (jac[p] <= jacobian_tolerance)
538  {
539  // Don't call print_info() recursively if we're already
540  // failing. print_info() calls Elem::volume() which may
541  // call FE::reinit() and trigger the same failure again.
542  static bool failing = false;
543  if (!failing)
544  {
545  failing = true;
546  elem->print_info(libMesh::err);
547  failing = false;
548  if (calculate_xyz)
549  {
550  libmesh_error_msg("ERROR: negative Jacobian " \
551  << jac[p] \
552  << " at point " \
553  << xyz[p] \
554  << " in element " \
555  << elem->id());
556  }
557  else
558  {
559  // In this case xyz[p] is not defined, so don't
560  // try to print it out.
561  libmesh_error_msg("ERROR: negative Jacobian " \
562  << jac[p] \
563  << " at point index " \
564  << p \
565  << " in element " \
566  << elem->id());
567  }
568  }
569  else
570  {
571  // We were already failing when we called this, so just
572  // stop the current computation and return with
573  // incomplete results.
574  return;
575  }
576  }
577 
578  // The inverse Jacobian entries also come from the
579  // generalized inverse of T (see also the 2D element
580  // living in 3D code).
581  const Real jacm2 = 1./jac[p]/jac[p];
582  dxidx_map[p] = jacm2*dxdxi_map(p);
583 #if LIBMESH_DIM > 1
584  dxidy_map[p] = jacm2*dydxi_map(p);
585 #endif
586 #if LIBMESH_DIM > 2
587  dxidz_map[p] = jacm2*dzdxi_map(p);
588 #endif
589 
590  JxW[p] = jac[p]*qw[p];
591  }
592 
593 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
594 
595  if (calculate_d2xyz)
596  {
597 #if LIBMESH_DIM == 1
598  // Compute inverse map second derivatives for 1D-element-living-in-1D case
600 #elif LIBMESH_DIM == 2
601  // Compute inverse map second derivatives for 1D-element-living-in-2D case
602  // See JWP notes for details
603 
604  // numer = x_xi*x_{xi xi} + y_xi*y_{xi xi}
605  Real numer =
606  dxyzdxi_map[p](0)*d2xyzdxi2_map[p](0) +
607  dxyzdxi_map[p](1)*d2xyzdxi2_map[p](1);
608 
609  // denom = (x_xi)^2 + (y_xi)^2 must be >= 0.0
610  Real denom =
611  dxyzdxi_map[p](0)*dxyzdxi_map[p](0) +
612  dxyzdxi_map[p](1)*dxyzdxi_map[p](1);
613 
614  if (denom <= 0.0)
615  {
616  // Don't call print_info() recursively if we're already
617  // failing. print_info() calls Elem::volume() which may
618  // call FE::reinit() and trigger the same failure again.
619  static bool failing = false;
620  if (!failing)
621  {
622  failing = true;
623  elem->print_info(libMesh::err);
624  failing = false;
625  libmesh_error_msg("Encountered invalid 1D element!");
626  }
627  else
628  {
629  // We were already failing when we called this, so just
630  // stop the current computation and return with
631  // incomplete results.
632  return;
633  }
634  }
635 
636  // xi_{x x}
637  d2xidxyz2_map[p][0] = -numer * dxidx_map[p]*dxidx_map[p] / denom;
638 
639  // xi_{x y}
640  d2xidxyz2_map[p][1] = -numer * dxidx_map[p]*dxidy_map[p] / denom;
641 
642  // xi_{y y}
643  d2xidxyz2_map[p][3] = -numer * dxidy_map[p]*dxidy_map[p] / denom;
644 
645 #elif LIBMESH_DIM == 3
646  // Compute inverse map second derivatives for 1D-element-living-in-3D case
647  // See JWP notes for details
648 
649  // numer = x_xi*x_{xi xi} + y_xi*y_{xi xi} + z_xi*z_{xi xi}
650  Real numer =
651  dxyzdxi_map[p](0)*d2xyzdxi2_map[p](0) +
652  dxyzdxi_map[p](1)*d2xyzdxi2_map[p](1) +
653  dxyzdxi_map[p](2)*d2xyzdxi2_map[p](2);
654 
655  // denom = (x_xi)^2 + (y_xi)^2 + (z_xi)^2 must be >= 0.0
656  Real denom =
657  dxyzdxi_map[p](0)*dxyzdxi_map[p](0) +
658  dxyzdxi_map[p](1)*dxyzdxi_map[p](1) +
659  dxyzdxi_map[p](2)*dxyzdxi_map[p](2);
660 
661  if (denom <= 0.0)
662  {
663  // Don't call print_info() recursively if we're already
664  // failing. print_info() calls Elem::volume() which may
665  // call FE::reinit() and trigger the same failure again.
666  static bool failing = false;
667  if (!failing)
668  {
669  failing = true;
670  elem->print_info(libMesh::err);
671  failing = false;
672  libmesh_error_msg("Encountered invalid 1D element!");
673  }
674  else
675  {
676  // We were already failing when we called this, so just
677  // stop the current computation and return with
678  // incomplete results.
679  return;
680  }
681  }
682 
683  // xi_{x x}
684  d2xidxyz2_map[p][0] = -numer * dxidx_map[p]*dxidx_map[p] / denom;
685 
686  // xi_{x y}
687  d2xidxyz2_map[p][1] = -numer * dxidx_map[p]*dxidy_map[p] / denom;
688 
689  // xi_{x z}
690  d2xidxyz2_map[p][2] = -numer * dxidx_map[p]*dxidz_map[p] / denom;
691 
692  // xi_{y y}
693  d2xidxyz2_map[p][3] = -numer * dxidy_map[p]*dxidy_map[p] / denom;
694 
695  // xi_{y z}
696  d2xidxyz2_map[p][4] = -numer * dxidy_map[p]*dxidz_map[p] / denom;
697 
698  // xi_{z z}
699  d2xidxyz2_map[p][5] = -numer * dxidz_map[p]*dxidz_map[p] / denom;
700 #endif //LIBMESH_DIM == 3
701  } // calculate_d2xyz
702 
703 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
704 
705  // done computing the map
706  break;
707  }
708 
709 
710  //--------------------------------------------------------------------
711  // 2D
712  case 2:
713  {
714  //------------------------------------------------------------------
715  // Compute the (x,y) values at the quadrature points,
716  // the Jacobian at the quadrature points
717 
718  if (calculate_xyz)
719  xyz[p].zero();
720 
721  if (calculate_dxyz)
722  {
723  dxyzdxi_map[p].zero();
724  dxyzdeta_map[p].zero();
725  }
726 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
727  if (calculate_d2xyz)
728  {
729  d2xyzdxi2_map[p].zero();
730  d2xyzdxideta_map[p].zero();
731  d2xyzdeta2_map[p].zero();
732  // Inverse map second derivatives
733  d2xidxyz2_map[p].assign(6, 0.);
734  d2etadxyz2_map[p].assign(6, 0.);
735  }
736 #endif
737 
738 
739  // compute (x,y) at the quadrature points, derivatives once
740  for (auto i : index_range(elem_nodes)) // sum over the nodes
741  {
742  // Reference to the point, helps eliminate
743  // excessive temporaries in the inner loop
744  libmesh_assert(elem_nodes[i]);
745  const Point & elem_point = *elem_nodes[i];
746 
747  if (calculate_xyz)
748  xyz[p].add_scaled (elem_point, phi_map[i][p] );
749 
750  if (calculate_dxyz)
751  {
752  dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
753  dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p]);
754  }
755 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
756  if (calculate_d2xyz)
757  {
758  d2xyzdxi2_map[p].add_scaled (elem_point, d2phidxi2_map[i][p]);
759  d2xyzdxideta_map[p].add_scaled (elem_point, d2phidxideta_map[i][p]);
760  d2xyzdeta2_map[p].add_scaled (elem_point, d2phideta2_map[i][p]);
761  }
762 #endif
763  }
764 
765  if (calculate_dxyz)
766  {
767  // compute the jacobian once
768  const Real dx_dxi = dxdxi_map(p),
769  dx_deta = dxdeta_map(p),
770  dy_dxi = dydxi_map(p),
771  dy_deta = dydeta_map(p);
772 
773 #if LIBMESH_DIM == 2
774  // Compute the Jacobian. This assumes the 2D face
775  // lives in 2D space
776  //
777  // Symbolically, the matrix determinant is
778  //
779  // | dx/dxi dx/deta |
780  // jac = | dy/dxi dy/deta |
781  //
782  // jac = dx/dxi*dy/deta - dx/deta*dy/dxi
783  jac[p] = (dx_dxi*dy_deta - dx_deta*dy_dxi);
784 
785  if (jac[p] <= jacobian_tolerance)
786  {
787  // Don't call print_info() recursively if we're already
788  // failing. print_info() calls Elem::volume() which may
789  // call FE::reinit() and trigger the same failure again.
790  static bool failing = false;
791  if (!failing)
792  {
793  failing = true;
794  elem->print_info(libMesh::err);
795  failing = false;
796  if (calculate_xyz)
797  {
798  libmesh_error_msg("ERROR: negative Jacobian " \
799  << jac[p] \
800  << " at point " \
801  << xyz[p] \
802  << " in element " \
803  << elem->id());
804  }
805  else
806  {
807  // In this case xyz[p] is not defined, so don't
808  // try to print it out.
809  libmesh_error_msg("ERROR: negative Jacobian " \
810  << jac[p] \
811  << " at point index " \
812  << p \
813  << " in element " \
814  << elem->id());
815  }
816  }
817  else
818  {
819  // We were already failing when we called this, so just
820  // stop the current computation and return with
821  // incomplete results.
822  return;
823  }
824  }
825 
826  JxW[p] = jac[p]*qw[p];
827 
828  // Compute the shape function derivatives wrt x,y at the
829  // quadrature points
830  const Real inv_jac = 1./jac[p];
831 
832  dxidx_map[p] = dy_deta*inv_jac; //dxi/dx = (1/J)*dy/deta
833  dxidy_map[p] = -dx_deta*inv_jac; //dxi/dy = -(1/J)*dx/deta
834  detadx_map[p] = -dy_dxi* inv_jac; //deta/dx = -(1/J)*dy/dxi
835  detady_map[p] = dx_dxi* inv_jac; //deta/dy = (1/J)*dx/dxi
836 
837  dxidz_map[p] = detadz_map[p] = 0.;
838 
839 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
840  if (compute_second_derivatives)
842 #endif
843 #else // LIBMESH_DIM == 3
844 
845  const Real dz_dxi = dzdxi_map(p),
846  dz_deta = dzdeta_map(p);
847 
848  // Compute the Jacobian. This assumes a 2D face in
849  // 3D space.
850  //
851  // The transformation matrix T from local to global
852  // coordinates is
853  //
854  // | dx/dxi dx/deta |
855  // T = | dy/dxi dy/deta |
856  // | dz/dxi dz/deta |
857  // note det(T' T) = det(T')det(T) = det(T)det(T)
858  // so det(T) = std::sqrt(det(T' T))
859  //
860  //----------------------------------------------
861  // Notes:
862  //
863  // dX = R dXi -> R'dX = R'R dXi
864  // (R^-1)dX = dXi [(R'R)^-1 R']dX = dXi
865  //
866  // so R^-1 = (R'R)^-1 R'
867  //
868  // and R^-1 R = (R'R)^-1 R'R = I.
869  //
870  const Real g11 = (dx_dxi*dx_dxi +
871  dy_dxi*dy_dxi +
872  dz_dxi*dz_dxi);
873 
874  const Real g12 = (dx_dxi*dx_deta +
875  dy_dxi*dy_deta +
876  dz_dxi*dz_deta);
877 
878  const Real g21 = g12;
879 
880  const Real g22 = (dx_deta*dx_deta +
881  dy_deta*dy_deta +
882  dz_deta*dz_deta);
883 
884  const Real det = (g11*g22 - g12*g21);
885 
886  if (det <= 0.)
887  {
888  // Don't call print_info() recursively if we're already
889  // failing. print_info() calls Elem::volume() which may
890  // call FE::reinit() and trigger the same failure again.
891  static bool failing = false;
892  if (!failing)
893  {
894  failing = true;
895  elem->print_info(libMesh::err);
896  failing = false;
897  if (calculate_xyz)
898  {
899  libmesh_error_msg("ERROR: negative Jacobian " \
900  << det \
901  << " at point " \
902  << xyz[p] \
903  << " in element " \
904  << elem->id());
905  }
906  else
907  {
908  // In this case xyz[p] is not defined, so don't
909  // try to print it out.
910  libmesh_error_msg("ERROR: negative Jacobian " \
911  << det \
912  << " at point index " \
913  << p \
914  << " in element " \
915  << elem->id());
916  }
917  }
918  else
919  {
920  // We were already failing when we called this, so just
921  // stop the current computation and return with
922  // incomplete results.
923  return;
924  }
925  }
926 
927  const Real inv_det = 1./det;
928  jac[p] = std::sqrt(det);
929 
930  JxW[p] = jac[p]*qw[p];
931 
932  const Real g11inv = g22*inv_det;
933  const Real g12inv = -g12*inv_det;
934  const Real g21inv = -g21*inv_det;
935  const Real g22inv = g11*inv_det;
936 
937  dxidx_map[p] = g11inv*dx_dxi + g12inv*dx_deta;
938  dxidy_map[p] = g11inv*dy_dxi + g12inv*dy_deta;
939  dxidz_map[p] = g11inv*dz_dxi + g12inv*dz_deta;
940 
941  detadx_map[p] = g21inv*dx_dxi + g22inv*dx_deta;
942  detady_map[p] = g21inv*dy_dxi + g22inv*dy_deta;
943  detadz_map[p] = g21inv*dz_dxi + g22inv*dz_deta;
944 
945 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
946 
947  if (calculate_d2xyz)
948  {
949  // Compute inverse map second derivative values for
950  // 2D-element-living-in-3D case. We pursue a least-squares
951  // solution approach for this "non-square" case, see JWP notes
952  // for details.
953 
954  // A = [ x_{xi xi} x_{eta eta} ]
955  // [ y_{xi xi} y_{eta eta} ]
956  // [ z_{xi xi} z_{eta eta} ]
957  DenseMatrix<Real> A(3,2);
958  A(0,0) = d2xyzdxi2_map[p](0); A(0,1) = d2xyzdeta2_map[p](0);
959  A(1,0) = d2xyzdxi2_map[p](1); A(1,1) = d2xyzdeta2_map[p](1);
960  A(2,0) = d2xyzdxi2_map[p](2); A(2,1) = d2xyzdeta2_map[p](2);
961 
962  // J^T, the transpose of the Jacobian matrix
963  DenseMatrix<Real> JT(2,3);
964  JT(0,0) = dx_dxi; JT(0,1) = dy_dxi; JT(0,2) = dz_dxi;
965  JT(1,0) = dx_deta; JT(1,1) = dy_deta; JT(1,2) = dz_deta;
966 
967  // (J^T J)^(-1), this has already been computed for us above...
968  DenseMatrix<Real> JTJinv(2,2);
969  JTJinv(0,0) = g11inv; JTJinv(0,1) = g12inv;
970  JTJinv(1,0) = g21inv; JTJinv(1,1) = g22inv;
971 
972  // Some helper variables
974  dxi (dxidx_map[p], dxidy_map[p], dxidz_map[p]),
975  deta (detadx_map[p], detady_map[p], detadz_map[p]);
976 
977  // To be filled in below
978  DenseVector<Real> tmp1(2);
979  DenseVector<Real> tmp2(3);
980  DenseVector<Real> tmp3(2);
981 
982  // For (s,t) in {(x,x), (x,y), (x,z), (y,y), (y,z), (z,z)}, compute the
983  // vector of inverse map second derivatives [xi_{s t}, eta_{s t}]
984  unsigned ctr=0;
985  for (unsigned s=0; s<3; ++s)
986  for (unsigned t=s; t<3; ++t)
987  {
988  // Construct tmp1 = [xi_s*xi_t, eta_s*eta_t]
989  tmp1(0) = dxi(s)*dxi(t);
990  tmp1(1) = deta(s)*deta(t);
991 
992  // Compute tmp2 = A * tmp1
993  A.vector_mult(tmp2, tmp1);
994 
995  // Compute scalar value "alpha"
996  Real alpha = dxi(s)*deta(t) + deta(s)*dxi(t);
997 
998  // Compute tmp2 <- tmp2 + alpha * x_{xi eta}
999  for (unsigned i=0; i<3; ++i)
1000  tmp2(i) += alpha*d2xyzdxideta_map[p](i);
1001 
1002  // Compute tmp3 = J^T * tmp2
1003  JT.vector_mult(tmp3, tmp2);
1004 
1005  // Compute tmp1 = (J^T J)^(-1) * tmp3. tmp1 is available for us to reuse.
1006  JTJinv.vector_mult(tmp1, tmp3);
1007 
1008  // Fill in appropriate entries, don't forget to multiply by -1!
1009  d2xidxyz2_map[p][ctr] = -tmp1(0);
1010  d2etadxyz2_map[p][ctr] = -tmp1(1);
1011 
1012  // Increment the counter
1013  ctr++;
1014  }
1015  }
1016 
1017 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1018 
1019 #endif // LIBMESH_DIM == 3
1020  }
1021  // done computing the map
1022  break;
1023  }
1024 
1025 
1026 
1027  //--------------------------------------------------------------------
1028  // 3D
1029  case 3:
1030  {
1031  //------------------------------------------------------------------
1032  // Compute the (x,y,z) values at the quadrature points,
1033  // the Jacobian at the quadrature point
1034 
1035  // Clear the entities that will be summed
1036  if (calculate_xyz)
1037  xyz[p].zero ();
1038  if (calculate_dxyz)
1039  {
1040  dxyzdxi_map[p].zero ();
1041  dxyzdeta_map[p].zero ();
1042  dxyzdzeta_map[p].zero ();
1043  }
1044 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1045  if (calculate_d2xyz)
1046  {
1047  d2xyzdxi2_map[p].zero();
1048  d2xyzdxideta_map[p].zero();
1049  d2xyzdxidzeta_map[p].zero();
1050  d2xyzdeta2_map[p].zero();
1051  d2xyzdetadzeta_map[p].zero();
1052  d2xyzdzeta2_map[p].zero();
1053  // Inverse map second derivatives
1054  d2xidxyz2_map[p].assign(6, 0.);
1055  d2etadxyz2_map[p].assign(6, 0.);
1056  d2zetadxyz2_map[p].assign(6, 0.);
1057  }
1058 #endif
1059 
1060 
1061  // compute (x,y,z) at the quadrature points,
1062  // dxdxi, dydxi, dzdxi,
1063  // dxdeta, dydeta, dzdeta,
1064  // dxdzeta, dydzeta, dzdzeta all once
1065  for (auto i : index_range(elem_nodes)) // sum over the nodes
1066  {
1067  // Reference to the point, helps eliminate
1068  // excessive temporaries in the inner loop
1069  libmesh_assert(elem_nodes[i]);
1070  const Point & elem_point = *elem_nodes[i];
1071 
1072  if (calculate_xyz)
1073  xyz[p].add_scaled (elem_point, phi_map[i][p] );
1074  if (calculate_dxyz)
1075  {
1076  dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
1077  dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p] );
1078  dxyzdzeta_map[p].add_scaled (elem_point, dphidzeta_map[i][p]);
1079  }
1080 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1081  if (calculate_d2xyz)
1082  {
1083  d2xyzdxi2_map[p].add_scaled (elem_point,
1084  d2phidxi2_map[i][p]);
1085  d2xyzdxideta_map[p].add_scaled (elem_point,
1086  d2phidxideta_map[i][p]);
1087  d2xyzdxidzeta_map[p].add_scaled (elem_point,
1088  d2phidxidzeta_map[i][p]);
1089  d2xyzdeta2_map[p].add_scaled (elem_point,
1090  d2phideta2_map[i][p]);
1091  d2xyzdetadzeta_map[p].add_scaled (elem_point,
1092  d2phidetadzeta_map[i][p]);
1093  d2xyzdzeta2_map[p].add_scaled (elem_point,
1094  d2phidzeta2_map[i][p]);
1095  }
1096 #endif
1097  }
1098 
1099  if (calculate_dxyz)
1100  {
1101  // compute the jacobian
1102  const Real
1103  dx_dxi = dxdxi_map(p), dy_dxi = dydxi_map(p), dz_dxi = dzdxi_map(p),
1104  dx_deta = dxdeta_map(p), dy_deta = dydeta_map(p), dz_deta = dzdeta_map(p),
1105  dx_dzeta = dxdzeta_map(p), dy_dzeta = dydzeta_map(p), dz_dzeta = dzdzeta_map(p);
1106 
1107  // Symbolically, the matrix determinant is
1108  //
1109  // | dx/dxi dy/dxi dz/dxi |
1110  // jac = | dx/deta dy/deta dz/deta |
1111  // | dx/dzeta dy/dzeta dz/dzeta |
1112  //
1113  // jac = dx/dxi*(dy/deta*dz/dzeta - dz/deta*dy/dzeta) +
1114  // dy/dxi*(dz/deta*dx/dzeta - dx/deta*dz/dzeta) +
1115  // dz/dxi*(dx/deta*dy/dzeta - dy/deta*dx/dzeta)
1116 
1117  jac[p] = (dx_dxi*(dy_deta*dz_dzeta - dz_deta*dy_dzeta) +
1118  dy_dxi*(dz_deta*dx_dzeta - dx_deta*dz_dzeta) +
1119  dz_dxi*(dx_deta*dy_dzeta - dy_deta*dx_dzeta));
1120 
1121  if (jac[p] <= jacobian_tolerance)
1122  {
1123  // Don't call print_info() recursively if we're already
1124  // failing. print_info() calls Elem::volume() which may
1125  // call FE::reinit() and trigger the same failure again.
1126  static bool failing = false;
1127  if (!failing)
1128  {
1129  failing = true;
1130  elem->print_info(libMesh::err);
1131  failing = false;
1132  if (calculate_xyz)
1133  {
1134  libmesh_error_msg("ERROR: negative Jacobian " \
1135  << jac[p] \
1136  << " at point " \
1137  << xyz[p] \
1138  << " in element " \
1139  << elem->id());
1140  }
1141  else
1142  {
1143  // In this case xyz[p] is not defined, so don't
1144  // try to print it out.
1145  libmesh_error_msg("ERROR: negative Jacobian " \
1146  << jac[p] \
1147  << " at point index " \
1148  << p \
1149  << " in element " \
1150  << elem->id());
1151  }
1152  }
1153  else
1154  {
1155  // We were already failing when we called this, so just
1156  // stop the current computation and return with
1157  // incomplete results.
1158  return;
1159  }
1160  }
1161 
1162  JxW[p] = jac[p]*qw[p];
1163 
1164  // Compute the shape function derivatives wrt x,y at the
1165  // quadrature points
1166  const Real inv_jac = 1./jac[p];
1167 
1168  dxidx_map[p] = (dy_deta*dz_dzeta - dz_deta*dy_dzeta)*inv_jac;
1169  dxidy_map[p] = (dz_deta*dx_dzeta - dx_deta*dz_dzeta)*inv_jac;
1170  dxidz_map[p] = (dx_deta*dy_dzeta - dy_deta*dx_dzeta)*inv_jac;
1171 
1172  detadx_map[p] = (dz_dxi*dy_dzeta - dy_dxi*dz_dzeta )*inv_jac;
1173  detady_map[p] = (dx_dxi*dz_dzeta - dz_dxi*dx_dzeta )*inv_jac;
1174  detadz_map[p] = (dy_dxi*dx_dzeta - dx_dxi*dy_dzeta )*inv_jac;
1175 
1176  dzetadx_map[p] = (dy_dxi*dz_deta - dz_dxi*dy_deta )*inv_jac;
1177  dzetady_map[p] = (dz_dxi*dx_deta - dx_dxi*dz_deta )*inv_jac;
1178  dzetadz_map[p] = (dx_dxi*dy_deta - dy_dxi*dx_deta )*inv_jac;
1179  }
1180 
1181 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1182  if (compute_second_derivatives)
1184 #endif
1185  // done computing the map
1186  break;
1187  }
1188 
1189  default:
1190  libmesh_error_msg("Invalid dim = " << dim);
1191  }
1192 }

References A, calculate_d2xyz, calculate_dxyz, calculate_xyz, calculations_started, compute_inverse_map_second_derivs(), d2etadxyz2_map, d2phideta2_map, d2phidetadzeta_map, d2phidxi2_map, d2phidxideta_map, d2phidxidzeta_map, d2phidzeta2_map, d2xidxyz2_map, d2xyzdeta2_map, d2xyzdetadzeta_map, d2xyzdxi2_map, d2xyzdxideta_map, d2xyzdxidzeta_map, d2xyzdzeta2_map, d2zetadxyz2_map, detadx_map, detady_map, detadz_map, dim, dphideta_map, dphidxi_map, dphidzeta_map, dxdeta_map(), dxdxi_map(), dxdzeta_map(), dxidx_map, dxidy_map, dxidz_map, dxyzdeta_map, dxyzdxi_map, dxyzdzeta_map, dydeta_map(), dydxi_map(), dydzeta_map(), dzdeta_map(), dzdxi_map(), dzdzeta_map(), dzetadx_map, dzetady_map, dzetadz_map, libMesh::err, libMesh::DofObject::id(), libMesh::index_range(), jac, jacobian_tolerance, JxW, libMesh::libmesh_assert(), phi_map, libMesh::Elem::print_info(), libMesh::Real, std::sqrt(), libMesh::DenseMatrix< T >::vector_mult(), and xyz.

Referenced by compute_affine_map(), and compute_map().

◆ determine_calculations()

void libMesh::FEMap::determine_calculations ( )
inlineprotected

Determine which values are to be calculated.

Definition at line 621 of file fe_map.h.

621  {
622  calculations_started = true;
623 
624 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
625  // Second derivative calculations currently have first derivative
626  // calculations as a prerequisite
627  if (calculate_d2xyz)
628  calculate_dxyz = true;
629 #endif
630  }

References calculate_d2xyz, calculate_dxyz, and calculations_started.

Referenced by compute_edge_map(), compute_face_map(), init_edge_shape_functions(), init_face_shape_functions(), init_reference_to_physical_map(), and resize_quadrature_map_vectors().

◆ dxdeta_map()

Real libMesh::FEMap::dxdeta_map ( const unsigned int  p) const
inlineprotected

Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.

Returns
The x value of the pth entry of the dxzydeta_map.

Definition at line 667 of file fe_map.h.

667 { return dxyzdeta_map[p](0); }

References dxyzdeta_map.

Referenced by libMesh::FEXYZMap::compute_face_map(), compute_face_map(), and compute_single_point_map().

◆ dxdxi_map()

Real libMesh::FEMap::dxdxi_map ( const unsigned int  p) const
inlineprotected

Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.

Returns
The x value of the pth entry of the dxzydxi_map.

Definition at line 643 of file fe_map.h.

643 { return dxyzdxi_map[p](0); }

References dxyzdxi_map.

Referenced by compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), and compute_single_point_map().

◆ dxdzeta_map()

Real libMesh::FEMap::dxdzeta_map ( const unsigned int  p) const
inlineprotected

Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.

Returns
The x value of the pth entry of the dxzydzeta_map.

Definition at line 691 of file fe_map.h.

691 { return dxyzdzeta_map[p](0); }

References dxyzdzeta_map.

Referenced by compute_single_point_map().

◆ dydeta_map()

Real libMesh::FEMap::dydeta_map ( const unsigned int  p) const
inlineprotected

Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.

Returns
The y value of the pth entry of the dxzydeta_map.

Definition at line 675 of file fe_map.h.

675 { return dxyzdeta_map[p](1); }

References dxyzdeta_map.

Referenced by libMesh::FEXYZMap::compute_face_map(), compute_face_map(), and compute_single_point_map().

◆ dydxi_map()

Real libMesh::FEMap::dydxi_map ( const unsigned int  p) const
inlineprotected

Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.

Returns
The y value of the pth entry of the dxzydxi_map.

Definition at line 651 of file fe_map.h.

651 { return dxyzdxi_map[p](1); }

References dxyzdxi_map.

Referenced by compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), and compute_single_point_map().

◆ dydzeta_map()

Real libMesh::FEMap::dydzeta_map ( const unsigned int  p) const
inlineprotected

Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.

Returns
The y value of the pth entry of the dxzydzeta_map.

Definition at line 699 of file fe_map.h.

699 { return dxyzdzeta_map[p](1); }

References dxyzdzeta_map.

Referenced by compute_single_point_map().

◆ dzdeta_map()

Real libMesh::FEMap::dzdeta_map ( const unsigned int  p) const
inlineprotected

Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.

Returns
The z value of the pth entry of the dxzydeta_map.

Definition at line 683 of file fe_map.h.

683 { return dxyzdeta_map[p](2); }

References dxyzdeta_map.

Referenced by libMesh::FEXYZMap::compute_face_map(), compute_face_map(), and compute_single_point_map().

◆ dzdxi_map()

Real libMesh::FEMap::dzdxi_map ( const unsigned int  p) const
inlineprotected

Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.

Returns
The z value of the pth entry of the dxzydxi_map.

Definition at line 659 of file fe_map.h.

659 { return dxyzdxi_map[p](2); }

References dxyzdxi_map.

Referenced by compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), and compute_single_point_map().

◆ dzdzeta_map()

Real libMesh::FEMap::dzdzeta_map ( const unsigned int  p) const
inlineprotected

Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.

Returns
The z value of the pth entry of the dxzydzeta_map.

Definition at line 707 of file fe_map.h.

707 { return dxyzdzeta_map[p](2); }

References dxyzdzeta_map.

Referenced by compute_single_point_map().

◆ get_curvatures()

const std::vector<Real>& libMesh::FEMap::get_curvatures ( ) const
inline
Returns
The curvatures for use in face integration.

Definition at line 432 of file fe_map.h.

References calculate_d2xyz, calculations_started, curvatures, and libMesh::libmesh_assert().

◆ get_d2etadxyz2()

const std::vector<std::vector<Real> >& libMesh::FEMap::get_d2etadxyz2 ( ) const
inline

Second derivatives of "eta" reference coordinate wrt physical coordinates.

Definition at line 367 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2etadxyz2_map, and libMesh::libmesh_assert().

Referenced by libMesh::H1FETransformation< OutputShape >::map_d2phi().

◆ get_d2phideta2_map()

std::vector<std::vector<Real> >& libMesh::FEMap::get_d2phideta2_map ( )
inline
Returns
The reference to physical map 2nd derivative

Definition at line 581 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2phideta2_map, and libMesh::libmesh_assert().

◆ get_d2phidetadzeta_map()

std::vector<std::vector<Real> >& libMesh::FEMap::get_d2phidetadzeta_map ( )
inline
Returns
The reference to physical map 2nd derivative

Definition at line 588 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2phidetadzeta_map, and libMesh::libmesh_assert().

◆ get_d2phidxi2_map()

std::vector<std::vector<Real> >& libMesh::FEMap::get_d2phidxi2_map ( )
inline
Returns
The reference to physical map 2nd derivative

Definition at line 560 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2phidxi2_map, and libMesh::libmesh_assert().

◆ get_d2phidxideta_map()

std::vector<std::vector<Real> >& libMesh::FEMap::get_d2phidxideta_map ( )
inline
Returns
The reference to physical map 2nd derivative

Definition at line 567 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2phidxideta_map, and libMesh::libmesh_assert().

◆ get_d2phidxidzeta_map()

std::vector<std::vector<Real> >& libMesh::FEMap::get_d2phidxidzeta_map ( )
inline
Returns
The reference to physical map 2nd derivative

Definition at line 574 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2phidxidzeta_map, and libMesh::libmesh_assert().

◆ get_d2phidzeta2_map()

std::vector<std::vector<Real> >& libMesh::FEMap::get_d2phidzeta2_map ( )
inline
Returns
The reference to physical map 2nd derivative

Definition at line 595 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2phidzeta2_map, and libMesh::libmesh_assert().

◆ get_d2psideta2() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_d2psideta2 ( )
inline
Returns
The reference to physical map 2nd derivative for the side/edge

Definition at line 514 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2psideta2_map, and libMesh::libmesh_assert().

◆ get_d2psideta2() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_d2psideta2 ( ) const
inline
Returns
const reference to physical map 2nd derivative for the side/edge

Definition at line 522 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2psideta2_map, and libMesh::libmesh_assert().

◆ get_d2psidxi2() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_d2psidxi2 ( )
inline
Returns
The reference to physical map 2nd derivative for the side/edge

Definition at line 486 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2psidxi2_map, and libMesh::libmesh_assert().

◆ get_d2psidxi2() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_d2psidxi2 ( ) const
inline
Returns
const reference to physical map 2nd derivative for the side/edge

Definition at line 493 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2psidxi2_map, and libMesh::libmesh_assert().

◆ get_d2psidxideta() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_d2psidxideta ( )
inline
Returns
The reference to physical map 2nd derivative for the side/edge

Definition at line 500 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2psidxideta_map, and libMesh::libmesh_assert().

◆ get_d2psidxideta() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_d2psidxideta ( ) const
inline
Returns
const reference to physical map 2nd derivative for the side/edge

Definition at line 507 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2psidxideta_map, and libMesh::libmesh_assert().

◆ get_d2xidxyz2()

const std::vector<std::vector<Real> >& libMesh::FEMap::get_d2xidxyz2 ( ) const
inline

◆ get_d2xyzdeta2()

const std::vector<RealGradient>& libMesh::FEMap::get_d2xyzdeta2 ( ) const
inline
Returns
The second partial derivatives in eta.

Definition at line 250 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2xyzdeta2_map, and libMesh::libmesh_assert().

◆ get_d2xyzdetadzeta()

const std::vector<RealGradient>& libMesh::FEMap::get_d2xyzdetadzeta ( ) const
inline
Returns
The second partial derivatives in eta-zeta.

Definition at line 278 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2xyzdetadzeta_map, and libMesh::libmesh_assert().

◆ get_d2xyzdxi2()

const std::vector<RealGradient>& libMesh::FEMap::get_d2xyzdxi2 ( ) const
inline
Returns
The second partial derivatives in xi.

Definition at line 243 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2xyzdxi2_map, and libMesh::libmesh_assert().

◆ get_d2xyzdxideta()

const std::vector<RealGradient>& libMesh::FEMap::get_d2xyzdxideta ( ) const
inline
Returns
The second partial derivatives in xi-eta.

Definition at line 264 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2xyzdxideta_map, and libMesh::libmesh_assert().

◆ get_d2xyzdxidzeta()

const std::vector<RealGradient>& libMesh::FEMap::get_d2xyzdxidzeta ( ) const
inline
Returns
The second partial derivatives in xi-zeta.

Definition at line 271 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2xyzdxidzeta_map, and libMesh::libmesh_assert().

◆ get_d2xyzdzeta2()

const std::vector<RealGradient>& libMesh::FEMap::get_d2xyzdzeta2 ( ) const
inline
Returns
The second partial derivatives in zeta.

Definition at line 257 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2xyzdzeta2_map, and libMesh::libmesh_assert().

◆ get_d2zetadxyz2()

const std::vector<std::vector<Real> >& libMesh::FEMap::get_d2zetadxyz2 ( ) const
inline

Second derivatives of "zeta" reference coordinate wrt physical coordinates.

Definition at line 374 of file fe_map.h.

References calculate_d2xyz, calculations_started, d2zetadxyz2_map, and libMesh::libmesh_assert().

Referenced by libMesh::H1FETransformation< OutputShape >::map_d2phi().

◆ get_detadx()

const std::vector<Real>& libMesh::FEMap::get_detadx ( ) const
inline

◆ get_detady()

const std::vector<Real>& libMesh::FEMap::get_detady ( ) const
inline

◆ get_detadz()

const std::vector<Real>& libMesh::FEMap::get_detadz ( ) const
inline

◆ get_dphideta_map() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_dphideta_map ( )
inline
Returns
The reference to physical map derivative

Definition at line 545 of file fe_map.h.

References calculate_dxyz, calculations_started, dphideta_map, and libMesh::libmesh_assert().

◆ get_dphideta_map() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_dphideta_map ( ) const
inline
Returns
The reference to physical map derivative

Definition at line 402 of file fe_map.h.

References calculate_dxyz, calculations_started, dphideta_map, and libMesh::libmesh_assert().

◆ get_dphidxi_map() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_dphidxi_map ( )
inline
Returns
The reference to physical map derivative

Definition at line 538 of file fe_map.h.

References calculate_dxyz, calculations_started, dphidxi_map, and libMesh::libmesh_assert().

◆ get_dphidxi_map() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_dphidxi_map ( ) const
inline
Returns
The reference to physical map derivative

Definition at line 395 of file fe_map.h.

References calculate_dxyz, calculations_started, dphidxi_map, and libMesh::libmesh_assert().

◆ get_dphidzeta_map() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_dphidzeta_map ( )
inline
Returns
The reference to physical map derivative

Definition at line 552 of file fe_map.h.

References calculate_dxyz, calculations_started, dphidzeta_map, and libMesh::libmesh_assert().

◆ get_dphidzeta_map() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_dphidzeta_map ( ) const
inline
Returns
The reference to physical map derivative

Definition at line 409 of file fe_map.h.

References calculate_dxyz, calculations_started, dphidzeta_map, and libMesh::libmesh_assert().

◆ get_dpsideta() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_dpsideta ( )
inline
Returns
The reference to physical map derivative for the side/edge

Definition at line 473 of file fe_map.h.

References calculate_dxyz, calculations_started, dpsideta_map, and libMesh::libmesh_assert().

◆ get_dpsideta() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_dpsideta ( ) const
inline

◆ get_dpsidxi() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_dpsidxi ( )
inline
Returns
The reference to physical map derivative for the side/edge

Definition at line 462 of file fe_map.h.

References calculate_dxyz, calculations_started, dpsidxi_map, and libMesh::libmesh_assert().

◆ get_dpsidxi() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_dpsidxi ( ) const
inline

◆ get_dxidx()

const std::vector<Real>& libMesh::FEMap::get_dxidx ( ) const
inline

◆ get_dxidy()

const std::vector<Real>& libMesh::FEMap::get_dxidy ( ) const
inline

◆ get_dxidz()

const std::vector<Real>& libMesh::FEMap::get_dxidz ( ) const
inline

◆ get_dxyzdeta()

const std::vector<RealGradient>& libMesh::FEMap::get_dxyzdeta ( ) const
inline
Returns
The element tangents in eta-direction at the quadrature points.

Definition at line 226 of file fe_map.h.

References calculate_dxyz, calculations_started, dxyzdeta_map, and libMesh::libmesh_assert().

Referenced by libMesh::HCurlFETransformation< OutputShape >::map_curl().

◆ get_dxyzdxi()

const std::vector<RealGradient>& libMesh::FEMap::get_dxyzdxi ( ) const
inline
Returns
The element tangents in xi-direction at the quadrature points.

Definition at line 218 of file fe_map.h.

References calculate_dxyz, calculations_started, dxyzdxi_map, and libMesh::libmesh_assert().

Referenced by libMesh::HCurlFETransformation< OutputShape >::map_curl().

◆ get_dxyzdzeta()

const std::vector<RealGradient>& libMesh::FEMap::get_dxyzdzeta ( ) const
inline
Returns
The element tangents in zeta-direction at the quadrature points.

Definition at line 234 of file fe_map.h.

References calculate_dxyz, calculations_started, dxyzdzeta_map, and libMesh::libmesh_assert().

Referenced by libMesh::HCurlFETransformation< OutputShape >::map_curl().

◆ get_dzetadx()

const std::vector<Real>& libMesh::FEMap::get_dzetadx ( ) const
inline

◆ get_dzetady()

const std::vector<Real>& libMesh::FEMap::get_dzetady ( ) const
inline

◆ get_dzetadz()

const std::vector<Real>& libMesh::FEMap::get_dzetadz ( ) const
inline

◆ get_jacobian()

const std::vector<Real>& libMesh::FEMap::get_jacobian ( ) const
inline
Returns
The element Jacobian for each quadrature point.

Definition at line 202 of file fe_map.h.

204  calculate_dxyz = true; return jac; }

References calculate_dxyz, calculations_started, jac, and libMesh::libmesh_assert().

Referenced by libMesh::HCurlFETransformation< OutputShape >::map_curl().

◆ get_JxW() [1/2]

std::vector<Real>& libMesh::FEMap::get_JxW ( )
inline
Returns
Writable reference to the element Jacobian times the quadrature weight for each quadrature point.

Definition at line 606 of file fe_map.h.

608  calculate_dxyz = true; return JxW; }

References calculate_dxyz, calculations_started, JxW, and libMesh::libmesh_assert().

◆ get_JxW() [2/2]

const std::vector<Real>& libMesh::FEMap::get_JxW ( ) const
inline
Returns
The element Jacobian times the quadrature weight for each quadrature point.

Definition at line 210 of file fe_map.h.

212  calculate_dxyz = true; return JxW; }

References calculate_dxyz, calculations_started, JxW, and libMesh::libmesh_assert().

◆ get_normals()

const std::vector<Point>& libMesh::FEMap::get_normals ( ) const
inline
Returns
The outward pointing normal vectors for face integration.

Definition at line 423 of file fe_map.h.

References calculate_dxyz, calculations_started, libMesh::libmesh_assert(), and normals.

◆ get_phi_map() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_phi_map ( )
inline
Returns
The reference to physical map for the element

Definition at line 531 of file fe_map.h.

References calculate_xyz, calculations_started, libMesh::libmesh_assert(), and phi_map.

◆ get_phi_map() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_phi_map ( ) const
inline
Returns
The reference to physical map for the element

Definition at line 388 of file fe_map.h.

References calculate_xyz, calculations_started, libMesh::libmesh_assert(), and phi_map.

◆ get_psi() [1/2]

std::vector<std::vector<Real> >& libMesh::FEMap::get_psi ( )
inline
Returns
The reference to physical map for the side/edge

Definition at line 456 of file fe_map.h.

457  { return psi_map; }

References psi_map.

◆ get_psi() [2/2]

const std::vector<std::vector<Real> >& libMesh::FEMap::get_psi ( ) const
inline
Returns
The reference to physical map for the side/edge

Definition at line 382 of file fe_map.h.

383  { return psi_map; }

References psi_map.

◆ get_tangents()

const std::vector<std::vector<Point> >& libMesh::FEMap::get_tangents ( ) const
inline
Returns
The tangent vectors for face integration.

Definition at line 416 of file fe_map.h.

References calculate_dxyz, calculations_started, libMesh::libmesh_assert(), and tangents.

◆ get_xyz()

const std::vector<Point>& libMesh::FEMap::get_xyz ( ) const
inline
Returns
The xyz spatial locations of the quadrature points on the element.

Definition at line 195 of file fe_map.h.

197  calculate_xyz = true; return xyz; }

References calculate_xyz, calculations_started, libMesh::libmesh_assert(), and xyz.

◆ init_edge_shape_functions()

template<unsigned int Dim>
void libMesh::FEMap::init_edge_shape_functions ( const std::vector< Point > &  qp,
const Elem edge 
)

Same as before, but for an edge.

This is used for some projection operators.

Definition at line 510 of file fe_boundary.C.

512 {
513  // Start logging the shape function initialization
514  LOG_SCOPE("init_edge_shape_functions()", "FEMap");
515 
516  libmesh_assert(edge);
517 
518  // We're calculating now!
519  this->determine_calculations();
520 
521  // The element type and order to use in
522  // the map
523  const FEFamily mapping_family = FEMap::map_fe_type(*edge);
524  const Order mapping_order (edge->default_order());
525  const ElemType mapping_elem_type (edge->type());
526  const FEType map_fe_type(mapping_order, mapping_family);
527 
528  // The number of quadrature points.
529  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
530 
531  const unsigned int n_mapping_shape_functions =
532  FEInterface::n_shape_functions(Dim, map_fe_type, mapping_elem_type);
533 
534  // resize the vectors to hold current data
535  // Psi are the shape functions used for the FE mapping
536  if (calculate_xyz)
537  this->psi_map.resize (n_mapping_shape_functions);
538  if (calculate_dxyz)
539  this->dpsidxi_map.resize (n_mapping_shape_functions);
540 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
541  if (calculate_d2xyz)
542  this->d2psidxi2_map.resize (n_mapping_shape_functions);
543 #endif
544 
545  FEInterface::shape_ptr shape_ptr =
547 
548  FEInterface::shape_deriv_ptr shape_deriv_ptr =
550 
551 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
552  FEInterface::shape_second_deriv_ptr shape_second_deriv_ptr =
554 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
555 
556  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
557  {
558  // Allocate space to store the values of the shape functions
559  // and their first and second derivatives at the quadrature points.
560  if (calculate_xyz)
561  this->psi_map[i].resize (n_qp);
562  if (calculate_dxyz)
563  this->dpsidxi_map[i].resize (n_qp);
564 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
565  if (calculate_d2xyz)
566  this->d2psidxi2_map[i].resize (n_qp);
567 #endif
568 
569  // Compute the value of mapping shape function i, and its first
570  // and second derivatives at quadrature point p
571  for (unsigned int p=0; p<n_qp; p++)
572  {
573  if (calculate_xyz)
574  this->psi_map[i][p] = shape_ptr (edge, mapping_order, i, qp[p], false);
575  if (calculate_dxyz)
576  this->dpsidxi_map[i][p] = shape_deriv_ptr (edge, mapping_order, i, 0, qp[p], false);
577 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
578  if (calculate_d2xyz)
579  this->d2psidxi2_map[i][p] = shape_second_deriv_ptr(edge, mapping_order, i, 0, qp[p], false);
580 #endif
581  }
582  }
583 }

References calculate_d2xyz, calculate_dxyz, calculate_xyz, d2psidxi2_map, libMesh::Elem::default_order(), determine_calculations(), dpsidxi_map, libMesh::libmesh_assert(), map_fe_type(), libMesh::FEInterface::n_shape_functions(), psi_map, libMesh::FEInterface::shape_deriv_function(), libMesh::FEInterface::shape_function(), libMesh::FEInterface::shape_second_deriv_function(), and libMesh::Elem::type().

◆ init_face_shape_functions()

template<unsigned int Dim>
void libMesh::FEMap::init_face_shape_functions ( const std::vector< Point > &  qp,
const Elem side 
)

Initializes the reference to physical element map for a side.

This is used for boundary integration.

Definition at line 381 of file fe_boundary.C.

383 {
384  // Start logging the shape function initialization
385  LOG_SCOPE("init_face_shape_functions()", "FEMap");
386 
387  libmesh_assert(side);
388 
389  // We're calculating now!
390  this->determine_calculations();
391 
392  // The element type and order to use in
393  // the map
394  const FEFamily mapping_family = FEMap::map_fe_type(*side);
395  const Order mapping_order (side->default_order());
396  const ElemType mapping_elem_type (side->type());
397  const FEType map_fe_type(mapping_order, mapping_family);
398 
399  // The number of quadrature points.
400  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
401 
402  const unsigned int n_mapping_shape_functions =
403  FEInterface::n_shape_functions(Dim, map_fe_type, mapping_elem_type);
404 
405  // resize the vectors to hold current data
406  // Psi are the shape functions used for the FE mapping
407  if (calculate_xyz)
408  this->psi_map.resize (n_mapping_shape_functions);
409 
410  if (Dim > 1)
411  {
412  if (calculate_dxyz)
413  this->dpsidxi_map.resize (n_mapping_shape_functions);
414 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
415  if (calculate_d2xyz)
416  this->d2psidxi2_map.resize (n_mapping_shape_functions);
417 #endif
418  }
419 
420  if (Dim == 3)
421  {
422  if (calculate_dxyz)
423  this->dpsideta_map.resize (n_mapping_shape_functions);
424 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
425  if (calculate_d2xyz)
426  {
427  this->d2psidxideta_map.resize (n_mapping_shape_functions);
428  this->d2psideta2_map.resize (n_mapping_shape_functions);
429  }
430 #endif
431  }
432 
433  FEInterface::shape_ptr shape_ptr =
435 
436  FEInterface::shape_deriv_ptr shape_deriv_ptr =
438 
439 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
440  FEInterface::shape_second_deriv_ptr shape_second_deriv_ptr =
442 #endif
443 
444  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
445  {
446  // Allocate space to store the values of the shape functions
447  // and their first and second derivatives at the quadrature points.
448  if (calculate_xyz)
449  this->psi_map[i].resize (n_qp);
450  if (Dim > 1)
451  {
452  if (calculate_dxyz)
453  this->dpsidxi_map[i].resize (n_qp);
454 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
455  if (calculate_d2xyz)
456  this->d2psidxi2_map[i].resize (n_qp);
457 #endif
458  }
459  if (Dim == 3)
460  {
461  if (calculate_dxyz)
462  this->dpsideta_map[i].resize (n_qp);
463 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
464  if (calculate_d2xyz)
465  {
466  this->d2psidxideta_map[i].resize (n_qp);
467  this->d2psideta2_map[i].resize (n_qp);
468  }
469 #endif
470  }
471 
472 
473  // Compute the value of mapping shape function i, and its first
474  // and second derivatives at quadrature point p
475  for (unsigned int p=0; p<n_qp; p++)
476  {
477  if (calculate_xyz)
478  this->psi_map[i][p] = shape_ptr (side, mapping_order, i, qp[p], false);
479  if (Dim > 1)
480  {
481  if (calculate_dxyz)
482  this->dpsidxi_map[i][p] = shape_deriv_ptr (side, mapping_order, i, 0, qp[p], false);
483 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
484  if (calculate_d2xyz)
485  this->d2psidxi2_map[i][p] = shape_second_deriv_ptr(side, mapping_order, i, 0, qp[p], false);
486 #endif
487  }
488  // libMesh::out << "this->d2psidxi2_map["<<i<<"][p]=" << d2psidxi2_map[i][p] << std::endl;
489 
490  // If we are in 3D, then our sides are 2D faces.
491  // For the second derivatives, we must also compute the cross
492  // derivative d^2() / dxi deta
493  if (Dim == 3)
494  {
495  if (calculate_dxyz)
496  this->dpsideta_map[i][p] = shape_deriv_ptr (side, mapping_order, i, 1, qp[p], false);
497 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
498  if (calculate_d2xyz)
499  {
500  this->d2psidxideta_map[i][p] = shape_second_deriv_ptr(side, mapping_order, i, 1, qp[p], false);
501  this->d2psideta2_map[i][p] = shape_second_deriv_ptr(side, mapping_order, i, 2, qp[p], false);
502  }
503 #endif
504  }
505  }
506  }
507 }

References calculate_d2xyz, calculate_dxyz, calculate_xyz, d2psideta2_map, d2psidxi2_map, d2psidxideta_map, libMesh::Elem::default_order(), determine_calculations(), dpsideta_map, dpsidxi_map, libMesh::libmesh_assert(), map_fe_type(), libMesh::FEInterface::n_shape_functions(), psi_map, libMesh::FEInterface::shape_deriv_function(), libMesh::FEInterface::shape_function(), libMesh::FEInterface::shape_second_deriv_function(), and libMesh::Elem::type().

◆ init_reference_to_physical_map()

template<unsigned int Dim>
template void libMesh::FEMap::init_reference_to_physical_map< 3 > ( const std::vector< Point > &  qp,
const Elem elem 
)

Definition at line 91 of file fe_map.C.

93 {
94  // Start logging the reference->physical map initialization
95  LOG_SCOPE("init_reference_to_physical_map()", "FEMap");
96 
97  // We're calculating now!
98  this->determine_calculations();
99 
100 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
101  if (elem->infinite())
102  {
103  //This mainly requires to change the FE<>-calls
104  // to FEInterface in this function.
105  libmesh_not_implemented();
106  }
107 #endif
108 
109  // The number of quadrature points.
110  const std::size_t n_qp = qp.size();
111 
112  // The element type and order to use in
113  // the map
114  const FEFamily mapping_family = FEMap::map_fe_type(*elem);
115  const Order mapping_order (elem->default_order());
116  const ElemType mapping_elem_type (elem->type());
117 
118  const FEType map_fe_type(mapping_order, mapping_family);
119 
120  // Number of shape functions used to construct the map
121  // (Lagrange shape functions are used for mapping)
122  const unsigned int n_mapping_shape_functions =
124  mapping_elem_type);
125 
126  if (calculate_xyz)
127  this->phi_map.resize (n_mapping_shape_functions);
128  if (Dim > 0)
129  {
130  if (calculate_dxyz)
131  this->dphidxi_map.resize (n_mapping_shape_functions);
132 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
133  if (calculate_d2xyz)
134  this->d2phidxi2_map.resize (n_mapping_shape_functions);
135 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
136  }
137 
138  if (Dim > 1)
139  {
140  if (calculate_dxyz)
141  this->dphideta_map.resize (n_mapping_shape_functions);
142 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
143  if (calculate_d2xyz)
144  {
145  this->d2phidxideta_map.resize (n_mapping_shape_functions);
146  this->d2phideta2_map.resize (n_mapping_shape_functions);
147  }
148 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
149  }
150 
151  if (Dim > 2)
152  {
153  if (calculate_dxyz)
154  this->dphidzeta_map.resize (n_mapping_shape_functions);
155 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
156  if (calculate_d2xyz)
157  {
158  this->d2phidxidzeta_map.resize (n_mapping_shape_functions);
159  this->d2phidetadzeta_map.resize (n_mapping_shape_functions);
160  this->d2phidzeta2_map.resize (n_mapping_shape_functions);
161  }
162 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
163  }
164 
165 
166  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
167  {
168  if (calculate_xyz)
169  this->phi_map[i].resize (n_qp);
170  if (Dim > 0)
171  {
172  if (calculate_dxyz)
173  this->dphidxi_map[i].resize (n_qp);
174 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
175  if (calculate_d2xyz)
176  {
177  this->d2phidxi2_map[i].resize (n_qp);
178  if (Dim > 1)
179  {
180  this->d2phidxideta_map[i].resize (n_qp);
181  this->d2phideta2_map[i].resize (n_qp);
182  }
183  if (Dim > 2)
184  {
185  this->d2phidxidzeta_map[i].resize (n_qp);
186  this->d2phidetadzeta_map[i].resize (n_qp);
187  this->d2phidzeta2_map[i].resize (n_qp);
188  }
189  }
190 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
191 
192  if (Dim > 1 && calculate_dxyz)
193  this->dphideta_map[i].resize (n_qp);
194 
195  if (Dim > 2 && calculate_dxyz)
196  this->dphidzeta_map[i].resize (n_qp);
197  }
198  }
199 
200  // Optimize for the *linear* geometric elements case:
201  bool is_linear = elem->is_linear();
202 
203  FEInterface::shape_ptr shape_ptr =
205 
206  FEInterface::shape_deriv_ptr shape_deriv_ptr =
208 
209 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
210  FEInterface::shape_second_deriv_ptr shape_second_deriv_ptr =
212 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
213 
214  switch (Dim)
215  {
216  //------------------------------------------------------------
217  // 0D
218  case 0:
219  {
220  if (calculate_xyz)
221  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
222  for (std::size_t p=0; p<n_qp; p++)
223  this->phi_map[i][p] =
224  shape_ptr(elem, mapping_order, i, qp[p], false);
225 
226  break;
227  }
228 
229  //------------------------------------------------------------
230  // 1D
231  case 1:
232  {
233  // Compute the value of the mapping shape function i at quadrature point p
234  // (Lagrange shape functions are used for mapping)
235  if (is_linear)
236  {
237  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
238  {
239  if (calculate_xyz)
240  this->phi_map[i][0] =
241  shape_ptr(elem, mapping_order, i, qp[0], false);
242 
243  if (calculate_dxyz)
244  this->dphidxi_map[i][0] =
245  shape_deriv_ptr(elem, mapping_order, i, 0, qp[0], false);
246 
247 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
248  if (calculate_d2xyz)
249  this->d2phidxi2_map[i][0] =
250  shape_second_deriv_ptr(elem, mapping_order, i, 0, qp[0], false);
251 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
252  for (std::size_t p=1; p<n_qp; p++)
253  {
254  if (calculate_xyz)
255  this->phi_map[i][p] =
256  shape_ptr(elem, mapping_order, i, qp[p], false);
257  if (calculate_dxyz)
258  this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
259 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
260  if (calculate_d2xyz)
261  this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
262 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
263  }
264  }
265  }
266  else
267  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
268  for (std::size_t p=0; p<n_qp; p++)
269  {
270  if (calculate_xyz)
271  this->phi_map[i][p] =
272  shape_ptr (elem, mapping_order, i, qp[p], false);
273  if (calculate_dxyz)
274  this->dphidxi_map[i][p] = shape_deriv_ptr (elem, mapping_order, i, 0, qp[p], false);
275 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
276  if (calculate_d2xyz)
277  this->d2phidxi2_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 0, qp[p], false);
278 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
279  }
280 
281  break;
282  }
283  //------------------------------------------------------------
284  // 2D
285  case 2:
286  {
287  // Compute the value of the mapping shape function i at quadrature point p
288  // (Lagrange shape functions are used for mapping)
289  if (is_linear)
290  {
291  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
292  {
293  if (calculate_xyz)
294  this->phi_map[i][0] =
295  shape_ptr (elem, mapping_order, i, qp[0], false);
296  if (calculate_dxyz)
297  {
298  this->dphidxi_map[i][0] = shape_deriv_ptr (elem, mapping_order, i, 0, qp[0], false);
299  this->dphideta_map[i][0] = shape_deriv_ptr (elem, mapping_order, i, 1, qp[0], false);
300  }
301 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
302  if (calculate_d2xyz)
303  {
304  this->d2phidxi2_map[i][0] = shape_second_deriv_ptr (elem, mapping_order, i, 0, qp[0], false);
305  this->d2phidxideta_map[i][0] = shape_second_deriv_ptr (elem, mapping_order, i, 1, qp[0], false);
306  this->d2phideta2_map[i][0] = shape_second_deriv_ptr (elem, mapping_order, i, 2, qp[0], false);
307  }
308 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
309  for (std::size_t p=1; p<n_qp; p++)
310  {
311  if (calculate_xyz)
312  this->phi_map[i][p] =
313  shape_ptr (elem, mapping_order, i, qp[p], false);
314  if (calculate_dxyz)
315  {
316  this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
317  this->dphideta_map[i][p] = this->dphideta_map[i][0];
318  }
319 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
320  if (calculate_d2xyz)
321  {
322  this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
323  this->d2phidxideta_map[i][p] = this->d2phidxideta_map[i][0];
324  this->d2phideta2_map[i][p] = this->d2phideta2_map[i][0];
325  }
326 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
327  }
328  }
329  }
330  else
331  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
332  for (std::size_t p=0; p<n_qp; p++)
333  {
334  if (calculate_xyz)
335  this->phi_map[i][p] =
336  shape_ptr(elem, mapping_order, i, qp[p], false);
337  if (calculate_dxyz)
338  {
339  this->dphidxi_map[i][p] = shape_deriv_ptr (elem, mapping_order, i, 0, qp[p], false);
340  this->dphideta_map[i][p] = shape_deriv_ptr (elem, mapping_order, i, 1, qp[p], false);
341  }
342 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
343  if (calculate_d2xyz)
344  {
345  this->d2phidxi2_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 0, qp[p], false);
346  this->d2phidxideta_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 1, qp[p], false);
347  this->d2phideta2_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 2, qp[p], false);
348  }
349 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
350  }
351 
352  break;
353  }
354 
355  //------------------------------------------------------------
356  // 3D
357  case 3:
358  {
359  // Compute the value of the mapping shape function i at quadrature point p
360  // (Lagrange shape functions are used for mapping)
361  if (is_linear)
362  {
363  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
364  {
365  if (calculate_xyz)
366  this->phi_map[i][0] =
367  shape_ptr (elem, mapping_order, i, qp[0], false);
368  if (calculate_dxyz)
369  {
370  this->dphidxi_map[i][0] = shape_deriv_ptr (elem, mapping_order, i, 0, qp[0], false);
371  this->dphideta_map[i][0] = shape_deriv_ptr (elem, mapping_order, i, 1, qp[0], false);
372  this->dphidzeta_map[i][0] = shape_deriv_ptr (elem, mapping_order, i, 2, qp[0], false);
373  }
374 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
375  if (calculate_d2xyz)
376  {
377  this->d2phidxi2_map[i][0] = shape_second_deriv_ptr (elem, mapping_order, i, 0, qp[0], false);
378  this->d2phidxideta_map[i][0] = shape_second_deriv_ptr (elem, mapping_order, i, 1, qp[0], false);
379  this->d2phideta2_map[i][0] = shape_second_deriv_ptr (elem, mapping_order, i, 2, qp[0], false);
380  this->d2phidxidzeta_map[i][0] = shape_second_deriv_ptr (elem, mapping_order, i, 3, qp[0], false);
381  this->d2phidetadzeta_map[i][0] = shape_second_deriv_ptr (elem, mapping_order, i, 4, qp[0], false);
382  this->d2phidzeta2_map[i][0] = shape_second_deriv_ptr (elem, mapping_order, i, 5, qp[0], false);
383  }
384 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
385  for (std::size_t p=1; p<n_qp; p++)
386  {
387  if (calculate_xyz)
388  this->phi_map[i][p] =
389  shape_ptr (elem, mapping_order, i, qp[p], false);
390  if (calculate_dxyz)
391  {
392  this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
393  this->dphideta_map[i][p] = this->dphideta_map[i][0];
394  this->dphidzeta_map[i][p] = this->dphidzeta_map[i][0];
395  }
396 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
397  if (calculate_d2xyz)
398  {
399  this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
400  this->d2phidxideta_map[i][p] = this->d2phidxideta_map[i][0];
401  this->d2phideta2_map[i][p] = this->d2phideta2_map[i][0];
402  this->d2phidxidzeta_map[i][p] = this->d2phidxidzeta_map[i][0];
403  this->d2phidetadzeta_map[i][p] = this->d2phidetadzeta_map[i][0];
404  this->d2phidzeta2_map[i][p] = this->d2phidzeta2_map[i][0];
405  }
406 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
407  }
408  }
409  }
410  else
411  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
412  for (std::size_t p=0; p<n_qp; p++)
413  {
414  if (calculate_xyz)
415  this->phi_map[i][p] =
416  shape_ptr(elem, mapping_order, i, qp[p], false);
417  if (calculate_dxyz)
418  {
419  this->dphidxi_map[i][p] = shape_deriv_ptr (elem, mapping_order, i, 0, qp[p], false);
420  this->dphideta_map[i][p] = shape_deriv_ptr (elem, mapping_order, i, 1, qp[p], false);
421  this->dphidzeta_map[i][p] = shape_deriv_ptr (elem, mapping_order, i, 2, qp[p], false);
422  }
423 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
424  if (calculate_d2xyz)
425  {
426  this->d2phidxi2_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 0, qp[p], false);
427  this->d2phidxideta_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 1, qp[p], false);
428  this->d2phideta2_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 2, qp[p], false);
429  this->d2phidxidzeta_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 3, qp[p], false);
430  this->d2phidetadzeta_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 4, qp[p], false);
431  this->d2phidzeta2_map[i][p] = shape_second_deriv_ptr (elem, mapping_order, i, 5, qp[p], false);
432  }
433 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
434  }
435 
436  break;
437  }
438 
439  default:
440  libmesh_error_msg("Invalid Dim = " << Dim);
441  }
442 }

References calculate_d2xyz, calculate_dxyz, calculate_xyz, d2phideta2_map, d2phidetadzeta_map, d2phidxi2_map, d2phidxideta_map, d2phidxidzeta_map, d2phidzeta2_map, libMesh::Elem::default_order(), determine_calculations(), dphideta_map, dphidxi_map, dphidzeta_map, libMesh::Elem::infinite(), libMesh::Elem::is_linear(), map_fe_type(), libMesh::FEInterface::n_shape_functions(), phi_map, libMesh::FEInterface::shape_deriv_function(), libMesh::FEInterface::shape_function(), libMesh::FEInterface::shape_second_deriv_function(), and libMesh::Elem::type().

◆ inverse_map() [1/2]

Point libMesh::FEMap::inverse_map ( const unsigned int  dim,
const Elem elem,
const Point p,
const Real  tolerance = TOLERANCE,
const bool  secure = true 
)
static
Returns
The location (on the reference element) of the point p located in physical space. This function requires inverting the (possibly nonlinear) transformation map, so it is not trivial. The optional parameter tolerance defines how close is "good enough." The map inversion iteration computes the sequence \( \{ p_n \} \), and the iteration is terminated when \( \|p - p_n\| < \mbox{\texttt{tolerance}} \) The parameter secure (always assumed false in non-debug mode) switches on integrity-checks on the mapped points.

Definition at line 1622 of file fe_map.C.

1627 {
1628  libmesh_assert(elem);
1629  libmesh_assert_greater_equal (tolerance, 0.);
1630 
1631 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1632  if (elem->infinite())
1633  return InfFEMap::inverse_map(dim, elem, physical_point, tolerance,
1634  secure);
1635 #endif
1636 
1637  // Start logging the map inversion.
1638  LOG_SCOPE("inverse_map()", "FEMap");
1639 
1640  // How much did the point on the reference
1641  // element change by in this Newton step?
1642  Real inverse_map_error = 0.;
1643 
1644  // The point on the reference element. This is
1645  // the "initial guess" for Newton's method. The
1646  // centroid seems like a good idea, but computing
1647  // it is a little more intensive than, say taking
1648  // the zero point.
1649  //
1650  // Convergence should be insensitive of this choice
1651  // for "good" elements.
1652  Point p; // the zero point. No computation required
1653 
1654  // The number of iterations in the map inversion process.
1655  unsigned int cnt = 0;
1656 
1657  // The number of iterations after which we give up and declare
1658  // divergence
1659  const unsigned int max_cnt = 10;
1660 
1661  // The distance (in master element space) beyond which we give up
1662  // and declare divergence. This is no longer used...
1663  // Real max_step_length = 4.;
1664 
1665 
1666 
1667  // Newton iteration loop.
1668  do
1669  {
1670  // Where our current iterate \p p maps to.
1671  const Point physical_guess = map(dim, elem, p);
1672 
1673  // How far our current iterate is from the actual point.
1674  const Point delta = physical_point - physical_guess;
1675 
1676  // Increment in current iterate \p p, will be computed.
1677  Point dp;
1678 
1679 
1680  // The form of the map and how we invert it depends
1681  // on the dimension that we are in.
1682  switch (dim)
1683  {
1684  // ------------------------------------------------------------------
1685  // 0D map inversion is trivial
1686  case 0:
1687  {
1688  break;
1689  }
1690 
1691  // ------------------------------------------------------------------
1692  // 1D map inversion
1693  //
1694  // Here we find the point on a 1D reference element that maps to
1695  // the point \p physical_point in the domain. This is a bit tricky
1696  // since we do not want to assume that the point \p physical_point
1697  // is also in a 1D domain. In particular, this method might get
1698  // called on the edge of a 3D element, in which case
1699  // \p physical_point actually lives in 3D.
1700  case 1:
1701  {
1702  const Point dxi = map_deriv (dim, elem, 0, p);
1703 
1704  // Newton's method in this case looks like
1705  //
1706  // {X} - {X_n} = [J]*dp
1707  //
1708  // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x1 matrix
1709  // d(x,y,z)/dxi, and we seek dp, a scalar. Since the above
1710  // system is either overdetermined or rank-deficient, we will
1711  // solve the normal equations for this system
1712  //
1713  // [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
1714  //
1715  // which involves the trivial inversion of the scalar
1716  // G = [J]^T [J]
1717  const Real G = dxi*dxi;
1718 
1719  if (secure)
1720  libmesh_assert_greater (G, 0.);
1721 
1722  const Real Ginv = 1./G;
1723 
1724  const Real dxidelta = dxi*delta;
1725 
1726  dp(0) = Ginv*dxidelta;
1727 
1728  // No master elements have radius > 4, but sometimes we
1729  // can take a step that big while still converging
1730  // if (secure)
1731  // libmesh_assert_less (dp.size(), max_step_length);
1732 
1733  break;
1734  }
1735 
1736 
1737 
1738  // ------------------------------------------------------------------
1739  // 2D map inversion
1740  //
1741  // Here we find the point on a 2D reference element that maps to
1742  // the point \p physical_point in the domain. This is a bit tricky
1743  // since we do not want to assume that the point \p physical_point
1744  // is also in a 2D domain. In particular, this method might get
1745  // called on the face of a 3D element, in which case
1746  // \p physical_point actually lives in 3D.
1747  case 2:
1748  {
1749  const Point dxi = map_deriv (dim, elem, 0, p);
1750  const Point deta = map_deriv (dim, elem, 1, p);
1751 
1752  // Newton's method in this case looks like
1753  //
1754  // {X} - {X_n} = [J]*{dp}
1755  //
1756  // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x2 matrix
1757  // d(x,y,z)/d(xi,eta), and we seek {dp}, a 2x1 vector. Since
1758  // the above system is either over-determined or rank-deficient,
1759  // we will solve the normal equations for this system
1760  //
1761  // [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
1762  //
1763  // which involves the inversion of the 2x2 matrix
1764  // [G] = [J]^T [J]
1765  const Real
1766  G11 = dxi*dxi, G12 = dxi*deta,
1767  G21 = dxi*deta, G22 = deta*deta;
1768 
1769 
1770  const Real det = (G11*G22 - G12*G21);
1771 
1772  if (secure)
1773  libmesh_assert_not_equal_to (det, 0.);
1774 
1775  const Real inv_det = 1./det;
1776 
1777  const Real
1778  Ginv11 = G22*inv_det,
1779  Ginv12 = -G12*inv_det,
1780 
1781  Ginv21 = -G21*inv_det,
1782  Ginv22 = G11*inv_det;
1783 
1784 
1785  const Real dxidelta = dxi*delta;
1786  const Real detadelta = deta*delta;
1787 
1788  dp(0) = (Ginv11*dxidelta + Ginv12*detadelta);
1789  dp(1) = (Ginv21*dxidelta + Ginv22*detadelta);
1790 
1791  // No master elements have radius > 4, but sometimes we
1792  // can take a step that big while still converging
1793  // if (secure)
1794  // libmesh_assert_less (dp.size(), max_step_length);
1795 
1796  break;
1797  }
1798 
1799 
1800 
1801  // ------------------------------------------------------------------
1802  // 3D map inversion
1803  //
1804  // Here we find the point in a 3D reference element that maps to
1805  // the point \p physical_point in a 3D domain. Nothing special
1806  // has to happen here, since (unless the map is singular because
1807  // you have a BAD element) the map will be invertible and we can
1808  // apply Newton's method directly.
1809  case 3:
1810  {
1811  const Point dxi = map_deriv (dim, elem, 0, p);
1812  const Point deta = map_deriv (dim, elem, 1, p);
1813  const Point dzeta = map_deriv (dim, elem, 2, p);
1814 
1815  // Newton's method in this case looks like
1816  //
1817  // {X} = {X_n} + [J]*{dp}
1818  //
1819  // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x3 matrix
1820  // d(x,y,z)/d(xi,eta,zeta), and we seek {dp}, a 3x1 vector.
1821  // Since the above system is nonsingular for invertible maps
1822  // we will solve
1823  //
1824  // {dp} = [J]^-1 ({X} - {X_n})
1825  //
1826  // which involves the inversion of the 3x3 matrix [J]
1827  libmesh_try
1828  {
1829  RealTensorValue(dxi(0), deta(0), dzeta(0),
1830  dxi(1), deta(1), dzeta(1),
1831  dxi(2), deta(2), dzeta(2)).solve(delta, dp);
1832  }
1833  libmesh_catch (ConvergenceFailure &)
1834  {
1835  // We encountered a singular Jacobian. The value of
1836  // dp is zero, since it was never changed during the
1837  // call to RealTensorValue::solve(). We don't want to
1838  // continue iterating until max_cnt since there is no
1839  // update to the Newton iterate, and we don't want to
1840  // print the inverse_map_error value since it will
1841  // confusingly be 0. Therefore, in the secure case we
1842  // need to throw an error message while in the !secure
1843  // case we can just return a far away point.
1844  if (secure)
1845  {
1846  libMesh::err << "ERROR: Newton scheme encountered a singular Jacobian in element: "
1847  << elem->id()
1848  << std::endl;
1849 
1850  elem->print_info(libMesh::err);
1851 
1852  libmesh_error_msg("Exiting...");
1853  }
1854  else
1855  {
1856  for (unsigned int i=0; i != dim; ++i)
1857  p(i) = 1e6;
1858  return p;
1859  }
1860  }
1861 
1862  // No master elements have radius > 4, but sometimes we
1863  // can take a step that big while still converging
1864  // if (secure)
1865  // libmesh_assert_less (dp.size(), max_step_length);
1866 
1867  break;
1868  }
1869 
1870 
1871  // Some other dimension?
1872  default:
1873  libmesh_error_msg("Invalid dim = " << dim);
1874  } // end switch(Dim), dp now computed
1875 
1876 
1877 
1878  // ||P_n+1 - P_n||
1879  inverse_map_error = dp.norm();
1880 
1881  // P_n+1 = P_n + dp
1882  p.add (dp);
1883 
1884  // Increment the iteration count.
1885  cnt++;
1886 
1887  // Watch for divergence of Newton's
1888  // method. Here's how it goes:
1889  // (1) For good elements, we expect convergence in 10
1890  // iterations, with no too-large steps.
1891  // - If called with (secure == true) and we have not yet converged
1892  // print out a warning message.
1893  // - If called with (secure == true) and we have not converged in
1894  // 20 iterations abort
1895  // (2) This method may be called in cases when the target point is not
1896  // inside the element and we have no business expecting convergence.
1897  // For these cases if we have not converged in 10 iterations forget
1898  // about it.
1899  if (cnt > max_cnt)
1900  {
1901  // Warn about divergence when secure is true - this
1902  // shouldn't happen
1903  if (secure)
1904  {
1905  // Print every time in devel/dbg modes
1906 #ifndef NDEBUG
1907  libmesh_here();
1908  libMesh::err << "WARNING: Newton scheme has not converged in "
1909  << cnt << " iterations:" << std::endl
1910  << " physical_point="
1911  << physical_point
1912  << " physical_guess="
1913  << physical_guess
1914  << " dp="
1915  << dp
1916  << " p="
1917  << p
1918  << " error=" << inverse_map_error
1919  << " in element " << elem->id()
1920  << std::endl;
1921 
1922  elem->print_info(libMesh::err);
1923 #else
1924  // In optimized mode, just print once that an inverse_map() call
1925  // had trouble converging its Newton iteration.
1926  libmesh_do_once(libMesh::err << "WARNING: At least one element took more than "
1927  << max_cnt
1928  << " iterations to converge in inverse_map()...\n"
1929  << "Rerun in devel/dbg mode for more details."
1930  << std::endl;);
1931 
1932 #endif // NDEBUG
1933 
1934  if (cnt > 2*max_cnt)
1935  {
1936  libMesh::err << "ERROR: Newton scheme FAILED to converge in "
1937  << cnt
1938  << " iterations in element "
1939  << elem->id()
1940  << " for physical point = "
1941  << physical_point
1942  << std::endl;
1943 
1944  elem->print_info(libMesh::err);
1945 
1946  libmesh_error_msg("Exiting...");
1947  }
1948  }
1949  // Return a far off point when secure is false - this
1950  // should only happen when we're trying to map a point
1951  // that's outside the element
1952  else
1953  {
1954  for (unsigned int i=0; i != dim; ++i)
1955  p(i) = 1e6;
1956 
1957  return p;
1958  }
1959  }
1960  }
1961  while (inverse_map_error > tolerance);
1962 
1963 
1964 
1965  // If we are in debug mode do two sanity checks.
1966 #ifdef DEBUG
1967 
1968  if (secure)
1969  {
1970  // Make sure the point \p p on the reference element actually
1971  // does map to the point \p physical_point within a tolerance.
1972 
1973  const Point check = map (dim, elem, p);
1974  const Point diff = physical_point - check;
1975 
1976  if (diff.norm() > tolerance)
1977  {
1978  libmesh_here();
1979  libMesh::err << "WARNING: diff is "
1980  << diff.norm()
1981  << std::endl
1982  << " point="
1983  << physical_point;
1984  libMesh::err << " local=" << check;
1985  libMesh::err << " lref= " << p;
1986 
1987  elem->print_info(libMesh::err);
1988  }
1989 
1990  // Make sure the point \p p on the reference element actually
1991  // is
1992 
1993  if (!FEAbstract::on_reference_element(p, elem->type(), 2*tolerance))
1994  {
1995  libmesh_here();
1996  libMesh::err << "WARNING: inverse_map of physical point "
1997  << physical_point
1998  << " is not on element." << '\n';
1999  elem->print_info(libMesh::err);
2000  }
2001  }
2002 
2003 #endif
2004 
2005  return p;
2006 }

References libMesh::TypeVector< T >::add(), dim, libMesh::err, libMesh::DofObject::id(), libMesh::Elem::infinite(), libMesh::InfFEMap::inverse_map(), libMesh::libmesh_assert(), map(), map_deriv(), libMesh::TypeVector< T >::norm(), libMesh::FEAbstract::on_reference_element(), libMesh::Elem::print_info(), libMesh::Real, and libMesh::Elem::type().

Referenced by libMesh::HPCoarsenTest::add_projection(), assemble_ellipticdg(), assemble_poisson(), libMesh::FEMContext::build_new_fe(), libMesh::FEGenericBase< FEOutputType< T >::type >::coarsened_dof_values(), compute_face_map(), compute_jacobian(), libMesh::FEAbstract::compute_node_constraints(), libMesh::FEGenericBase< FEOutputType< T >::type >::compute_periodic_constraints(), libMesh::FEAbstract::compute_periodic_node_constraints(), libMesh::FEGenericBase< FEOutputType< T >::type >::compute_proj_constraints(), compute_residual(), libMesh::MeshFunction::discontinuous_gradient(), libMesh::MeshFunction::discontinuous_value(), libMesh::FE< Dim, LAGRANGE_VEC >::edge_reinit(), libMesh::DTKEvaluator::evaluate(), libMesh::MeshFunction::gradient(), libMesh::MeshFunction::hessian(), libMesh::InfFEMap::inverse_map(), inverse_map(), libMesh::FE< Dim, LAGRANGE_VEC >::inverse_map(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), libMesh::MeshFunction::operator()(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::ProjectEdges::operator()(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::ProjectSides::operator()(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::ProjectInteriors::operator()(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::Elem::point_test(), libMesh::System::point_value(), libMesh::JumpErrorEstimator::reinit_sides(), libMesh::HPCoarsenTest::select_refinement(), NavierSystem::side_constraint(), FETest< order, family, elem_type >::testGradU(), FETest< order, family, elem_type >::testGradUComp(), and FETest< order, family, elem_type >::testU().

◆ inverse_map() [2/2]

void libMesh::FEMap::inverse_map ( unsigned int  dim,
const Elem elem,
const std::vector< Point > &  physical_points,
std::vector< Point > &  reference_points,
const Real  tolerance = TOLERANCE,
const bool  secure = true 
)
static

Takes a number points in physical space (in the physical_points vector) and finds their location on the reference element for the input element elem.

The values on the reference element are returned in the vector reference_points. The optional parameter tolerance defines how close is "good enough." The map inversion iteration computes the sequence \( \{ p_n \} \), and the iteration is terminated when \( \|p - p_n\| < \mbox{\texttt{tolerance}} \) The parameter secure (always assumed false in non-debug mode) switches on integrity-checks on the mapped points.

Definition at line 2010 of file fe_map.C.

2016 {
2017 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2018  if (elem->infinite())
2019  {
2020  InfFEMap::inverse_map(dim, elem, physical_points, reference_points, tolerance, secure);
2021  return;
2022  // libmesh_not_implemented();
2023  }
2024 #endif
2025 
2026  // The number of points to find the
2027  // inverse map of
2028  const std::size_t n_points = physical_points.size();
2029 
2030  // Resize the vector to hold the points
2031  // on the reference element
2032  reference_points.resize(n_points);
2033 
2034  // Find the coordinates on the reference
2035  // element of each point in physical space
2036  for (std::size_t p=0; p<n_points; p++)
2037  reference_points[p] =
2038  inverse_map (dim, elem, physical_points[p], tolerance, secure);
2039 }

References dim, libMesh::Elem::infinite(), libMesh::InfFEMap::inverse_map(), and inverse_map().

◆ map()

Point libMesh::FEMap::map ( const unsigned int  dim,
const Elem elem,
const Point reference_point 
)
static
Returns
The location (in physical space) of the point p located on the reference element.

Definition at line 2043 of file fe_map.C.

2046 {
2047  libmesh_assert(elem);
2048 
2049 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2050  if (elem->infinite())
2051  return InfFEMap::map(dim, elem, reference_point);
2052 #endif
2053 
2054  Point p;
2055 
2056  const FEFamily mapping_family = FEMap::map_fe_type(*elem);
2057  const ElemType type = elem->type();
2058  const Order order = elem->default_order();
2059  const FEType fe_type (order, mapping_family);
2060 
2061  const unsigned int n_sf = FEInterface::n_shape_functions(dim, fe_type, type);
2062 
2063  FEInterface::shape_ptr shape_ptr =
2064  FEInterface::shape_function(dim, fe_type);
2065 
2066  // Lagrange basis functions are used for mapping
2067  for (unsigned int i=0; i<n_sf; i++)
2068  p.add_scaled (elem->point(i),
2069  shape_ptr(elem, order, i, reference_point, false));
2070 
2071  return p;
2072 }

References libMesh::TypeVector< T >::add_scaled(), libMesh::Elem::default_order(), dim, libMesh::Elem::infinite(), libMesh::libmesh_assert(), libMesh::InfFEMap::map(), map_fe_type(), libMesh::FEInterface::n_shape_functions(), libMesh::Elem::point(), libMesh::FEInterface::shape_function(), and libMesh::Elem::type().

Referenced by inverse_map(), libMesh::InfFEMap::map(), libMesh::FE< Dim, LAGRANGE_VEC >::map(), and libMesh::Elem::point_test().

◆ map_deriv()

Point libMesh::FEMap::map_deriv ( const unsigned int  dim,
const Elem elem,
const unsigned int  j,
const Point reference_point 
)
static
Returns
component j of d(xyz)/d(xi eta zeta) (in physical space) of the point p located on the reference element.

Definition at line 2076 of file fe_map.C.

2080 {
2081  libmesh_assert(elem);
2082 
2083 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2084  if (elem->infinite())
2085  libmesh_not_implemented();
2086 #endif
2087 
2088  Point p;
2089 
2090  const FEFamily mapping_family = FEMap::map_fe_type(*elem);
2091  const ElemType type = elem->type();
2092  const Order order = elem->default_order();
2093  const FEType fe_type (order, mapping_family);
2094  const unsigned int n_sf = FEInterface::n_shape_functions(dim, fe_type, type);
2095 
2096  FEInterface::shape_deriv_ptr shape_deriv_ptr =
2098 
2099  // Lagrange basis functions are used for mapping
2100  for (unsigned int i=0; i<n_sf; i++)
2101  p.add_scaled (elem->point(i),
2102  shape_deriv_ptr(elem, order, i, j, reference_point,
2103  false));
2104 
2105  return p;
2106 }

References libMesh::TypeVector< T >::add_scaled(), libMesh::Elem::default_order(), dim, libMesh::Elem::infinite(), libMesh::libmesh_assert(), map_fe_type(), libMesh::FEInterface::n_shape_functions(), libMesh::Elem::point(), libMesh::FEInterface::shape_deriv_function(), and libMesh::Elem::type().

Referenced by compute_face_map(), inverse_map(), libMesh::FE< Dim, LAGRANGE_VEC >::map_eta(), libMesh::FE< Dim, LAGRANGE_VEC >::map_xi(), and libMesh::FE< Dim, LAGRANGE_VEC >::map_zeta().

◆ map_fe_type()

FEFamily libMesh::FEMap::map_fe_type ( const Elem elem)
inlinestatic

◆ print_JxW()

void libMesh::FEMap::print_JxW ( std::ostream &  os) const

Prints the Jacobian times the weight for each quadrature point.

Definition at line 1486 of file fe_map.C.

1487 {
1488  for (auto i : index_range(JxW))
1489  os << " [" << i << "]: " << JxW[i] << std::endl;
1490 }

References libMesh::index_range(), and JxW.

◆ print_xyz()

void libMesh::FEMap::print_xyz ( std::ostream &  os) const

Prints the spatial location of each quadrature point (on the physical element).

Definition at line 1494 of file fe_map.C.

1495 {
1496  for (auto i : index_range(xyz))
1497  os << " [" << i << "]: " << xyz[i];
1498 }

References libMesh::index_range(), and xyz.

◆ resize_quadrature_map_vectors()

void libMesh::FEMap::resize_quadrature_map_vectors ( const unsigned int  dim,
unsigned int  n_qp 
)
protected

A utility function for use by compute_*_map.

Definition at line 1196 of file fe_map.C.

1197 {
1198  // We're calculating now!
1199  this->determine_calculations();
1200 
1201  // Resize the vectors to hold data at the quadrature points
1202  if (calculate_xyz)
1203  xyz.resize(n_qp);
1204  if (calculate_dxyz)
1205  {
1206  dxyzdxi_map.resize(n_qp);
1207  dxidx_map.resize(n_qp);
1208  dxidy_map.resize(n_qp); // 1D element may live in 2D ...
1209  dxidz_map.resize(n_qp); // ... or 3D
1210  }
1211 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1212  if (calculate_d2xyz)
1213  {
1214  d2xyzdxi2_map.resize(n_qp);
1215 
1216  // Inverse map second derivatives
1217  d2xidxyz2_map.resize(n_qp);
1218  for (auto i : index_range(d2xidxyz2_map))
1219  d2xidxyz2_map[i].assign(6, 0.);
1220  }
1221 #endif
1222  if (dim > 1)
1223  {
1224  if (calculate_dxyz)
1225  {
1226  dxyzdeta_map.resize(n_qp);
1227  detadx_map.resize(n_qp);
1228  detady_map.resize(n_qp);
1229  detadz_map.resize(n_qp);
1230  }
1231 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1232  if (calculate_d2xyz)
1233  {
1234  d2xyzdxideta_map.resize(n_qp);
1235  d2xyzdeta2_map.resize(n_qp);
1236 
1237  // Inverse map second derivatives
1238  d2etadxyz2_map.resize(n_qp);
1239  for (auto i : index_range(d2etadxyz2_map))
1240  d2etadxyz2_map[i].assign(6, 0.);
1241  }
1242 #endif
1243  if (dim > 2)
1244  {
1245  if (calculate_dxyz)
1246  {
1247  dxyzdzeta_map.resize (n_qp);
1248  dzetadx_map.resize (n_qp);
1249  dzetady_map.resize (n_qp);
1250  dzetadz_map.resize (n_qp);
1251  }
1252 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1253  if (calculate_d2xyz)
1254  {
1255  d2xyzdxidzeta_map.resize(n_qp);
1256  d2xyzdetadzeta_map.resize(n_qp);
1257  d2xyzdzeta2_map.resize(n_qp);
1258 
1259  // Inverse map second derivatives
1260  d2zetadxyz2_map.resize(n_qp);
1261  for (auto i : index_range(d2zetadxyz2_map))
1262  d2zetadxyz2_map[i].assign(6, 0.);
1263  }
1264 #endif
1265  }
1266  }
1267 
1268  if (calculate_dxyz)
1269  {
1270  jac.resize(n_qp);
1271  JxW.resize(n_qp);
1272  }
1273 }

References calculate_d2xyz, calculate_dxyz, calculate_xyz, d2etadxyz2_map, d2xidxyz2_map, d2xyzdeta2_map, d2xyzdetadzeta_map, d2xyzdxi2_map, d2xyzdxideta_map, d2xyzdxidzeta_map, d2xyzdzeta2_map, d2zetadxyz2_map, detadx_map, detady_map, detadz_map, determine_calculations(), dim, dxidx_map, dxidy_map, dxidz_map, dxyzdeta_map, dxyzdxi_map, dxyzdzeta_map, dzetadx_map, dzetady_map, dzetadz_map, libMesh::index_range(), jac, JxW, and xyz.

Referenced by compute_affine_map(), compute_map(), and compute_null_map().

◆ set_jacobian_tolerance()

void libMesh::FEMap::set_jacobian_tolerance ( Real  tol)
inline

Set the Jacobian tolerance used for determining when the mapping fails.

The mapping is determined to fail if jac <= jacobian_tolerance.

Definition at line 614 of file fe_map.h.

614 { jacobian_tolerance = tol; }

References jacobian_tolerance.

Friends And Related Function Documentation

◆ FE

template<unsigned int Dim, FEFamily T>
friend class FE
friend

FE classes should be able to reset calculations_started in a few special cases.

Definition at line 1004 of file fe_map.h.

Member Data Documentation

◆ _elem_nodes

std::vector<const Node *> libMesh::FEMap::_elem_nodes
private

Work vector for compute_affine_map()

Definition at line 1023 of file fe_map.h.

Referenced by compute_affine_map(), and compute_map().

◆ calculate_d2xyz

bool libMesh::FEMap::calculate_d2xyz
mutableprotected

◆ calculate_dxyz

bool libMesh::FEMap::calculate_dxyz
mutableprotected

◆ calculate_xyz

bool libMesh::FEMap::calculate_xyz
mutableprotected

◆ calculations_started

bool libMesh::FEMap::calculations_started
mutableprotected

◆ curvatures

std::vector<Real> libMesh::FEMap::curvatures
protected

The mean curvature (= one half the sum of the principal curvatures) on the boundary at the quadrature points.

The mean curvature is a scalar value.

Definition at line 960 of file fe_map.h.

Referenced by compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), and get_curvatures().

◆ d2etadxyz2_map

std::vector<std::vector<Real> > libMesh::FEMap::d2etadxyz2_map
protected

Second derivatives of "eta" reference coordinate wrt physical coordinates.

At each qp: (eta_{xx}, eta_{xy}, eta_{xz}, eta_{yy}, eta_{yz}, eta_{zz})

Definition at line 839 of file fe_map.h.

Referenced by compute_inverse_map_second_derivs(), compute_single_point_map(), get_d2etadxyz2(), and resize_quadrature_map_vectors().

◆ d2phideta2_map

std::vector<std::vector<Real> > libMesh::FEMap::d2phideta2_map
protected

Map for the second derivative, d^2(phi)/d(eta)^2.

Definition at line 888 of file fe_map.h.

Referenced by compute_single_point_map(), get_d2phideta2_map(), and init_reference_to_physical_map().

◆ d2phidetadzeta_map

std::vector<std::vector<Real> > libMesh::FEMap::d2phidetadzeta_map
protected

Map for the second derivative, d^2(phi)/d(eta)d(zeta).

Definition at line 893 of file fe_map.h.

Referenced by compute_single_point_map(), get_d2phidetadzeta_map(), and init_reference_to_physical_map().

◆ d2phidxi2_map

std::vector<std::vector<Real> > libMesh::FEMap::d2phidxi2_map
protected

Map for the second derivative, d^2(phi)/d(xi)^2.

Definition at line 873 of file fe_map.h.

Referenced by compute_single_point_map(), get_d2phidxi2_map(), and init_reference_to_physical_map().

◆ d2phidxideta_map

std::vector<std::vector<Real> > libMesh::FEMap::d2phidxideta_map
protected

Map for the second derivative, d^2(phi)/d(xi)d(eta).

Definition at line 878 of file fe_map.h.

Referenced by compute_single_point_map(), get_d2phidxideta_map(), and init_reference_to_physical_map().

◆ d2phidxidzeta_map

std::vector<std::vector<Real> > libMesh::FEMap::d2phidxidzeta_map
protected

Map for the second derivative, d^2(phi)/d(xi)d(zeta).

Definition at line 883 of file fe_map.h.

Referenced by compute_single_point_map(), get_d2phidxidzeta_map(), and init_reference_to_physical_map().

◆ d2phidzeta2_map

std::vector<std::vector<Real> > libMesh::FEMap::d2phidzeta2_map
protected

Map for the second derivative, d^2(phi)/d(zeta)^2.

Definition at line 898 of file fe_map.h.

Referenced by compute_single_point_map(), get_d2phidzeta2_map(), and init_reference_to_physical_map().

◆ d2psideta2_map

std::vector<std::vector<Real> > libMesh::FEMap::d2psideta2_map
protected

Map for the second derivatives (in eta) of the side shape functions.

Useful for computing the curvature at the quadrature points.

Definition at line 940 of file fe_map.h.

Referenced by libMesh::FEXYZMap::compute_face_map(), compute_face_map(), get_d2psideta2(), and init_face_shape_functions().

◆ d2psidxi2_map

std::vector<std::vector<Real> > libMesh::FEMap::d2psidxi2_map
protected

Map for the second derivatives (in xi) of the side shape functions.

Useful for computing the curvature at the quadrature points.

Definition at line 926 of file fe_map.h.

Referenced by compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), get_d2psidxi2(), init_edge_shape_functions(), and init_face_shape_functions().

◆ d2psidxideta_map

std::vector<std::vector<Real> > libMesh::FEMap::d2psidxideta_map
protected

Map for the second (cross) derivatives in xi, eta of the side shape functions.

Useful for computing the curvature at the quadrature points.

Definition at line 933 of file fe_map.h.

Referenced by libMesh::FEXYZMap::compute_face_map(), compute_face_map(), get_d2psidxideta(), and init_face_shape_functions().

◆ d2xidxyz2_map

std::vector<std::vector<Real> > libMesh::FEMap::d2xidxyz2_map
protected

Second derivatives of "xi" reference coordinate wrt physical coordinates.

At each qp: (xi_{xx}, xi_{xy}, xi_{xz}, xi_{yy}, xi_{yz}, xi_{zz})

Definition at line 833 of file fe_map.h.

Referenced by compute_inverse_map_second_derivs(), compute_single_point_map(), get_d2xidxyz2(), and resize_quadrature_map_vectors().

◆ d2xyzdeta2_map

std::vector<RealGradient> libMesh::FEMap::d2xyzdeta2_map
protected

◆ d2xyzdetadzeta_map

std::vector<RealGradient> libMesh::FEMap::d2xyzdetadzeta_map
protected

Vector of mixed second partial derivatives in eta-zeta: d^2(x)/d(eta)d(zeta) d^2(y)/d(eta)d(zeta) d^2(z)/d(eta)d(zeta)

Definition at line 762 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_d2xyzdetadzeta(), and resize_quadrature_map_vectors().

◆ d2xyzdxi2_map

std::vector<RealGradient> libMesh::FEMap::d2xyzdxi2_map
protected

Vector of second partial derivatives in xi: d^2(x)/d(xi)^2, d^2(y)/d(xi)^2, d^2(z)/d(xi)^2.

Definition at line 738 of file fe_map.h.

Referenced by compute_affine_map(), compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_d2xyzdxi2(), and resize_quadrature_map_vectors().

◆ d2xyzdxideta_map

std::vector<RealGradient> libMesh::FEMap::d2xyzdxideta_map
protected

Vector of mixed second partial derivatives in xi-eta: d^2(x)/d(xi)d(eta) d^2(y)/d(xi)d(eta) d^2(z)/d(xi)d(eta)

Definition at line 744 of file fe_map.h.

Referenced by compute_affine_map(), compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_d2xyzdxideta(), and resize_quadrature_map_vectors().

◆ d2xyzdxidzeta_map

std::vector<RealGradient> libMesh::FEMap::d2xyzdxidzeta_map
protected

Vector of second partial derivatives in xi-zeta: d^2(x)/d(xi)d(zeta), d^2(y)/d(xi)d(zeta), d^2(z)/d(xi)d(zeta)

Definition at line 756 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_d2xyzdxidzeta(), and resize_quadrature_map_vectors().

◆ d2xyzdzeta2_map

std::vector<RealGradient> libMesh::FEMap::d2xyzdzeta2_map
protected

Vector of second partial derivatives in zeta: d^2(x)/d(zeta)^2.

Definition at line 768 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_d2xyzdzeta2(), and resize_quadrature_map_vectors().

◆ d2zetadxyz2_map

std::vector<std::vector<Real> > libMesh::FEMap::d2zetadxyz2_map
protected

Second derivatives of "zeta" reference coordinate wrt physical coordinates.

At each qp: (zeta_{xx}, zeta_{xy}, zeta_{xz}, zeta_{yy}, zeta_{yz}, zeta_{zz})

Definition at line 845 of file fe_map.h.

Referenced by compute_inverse_map_second_derivs(), compute_single_point_map(), get_d2zetadxyz2(), and resize_quadrature_map_vectors().

◆ detadx_map

std::vector<Real> libMesh::FEMap::detadx_map
protected

Map for partial derivatives: d(eta)/d(x).

Needed for the Jacobian.

Definition at line 795 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_detadx(), and resize_quadrature_map_vectors().

◆ detady_map

std::vector<Real> libMesh::FEMap::detady_map
protected

Map for partial derivatives: d(eta)/d(y).

Needed for the Jacobian.

Definition at line 801 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_detady(), and resize_quadrature_map_vectors().

◆ detadz_map

std::vector<Real> libMesh::FEMap::detadz_map
protected

Map for partial derivatives: d(eta)/d(z).

Needed for the Jacobian.

Definition at line 807 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_detadz(), and resize_quadrature_map_vectors().

◆ dphideta_map

std::vector<std::vector<Real> > libMesh::FEMap::dphideta_map
protected

Map for the derivative, d(phi)/d(eta).

Definition at line 861 of file fe_map.h.

Referenced by compute_single_point_map(), get_dphideta_map(), and init_reference_to_physical_map().

◆ dphidxi_map

std::vector<std::vector<Real> > libMesh::FEMap::dphidxi_map
protected

Map for the derivative, d(phi)/d(xi).

Definition at line 856 of file fe_map.h.

Referenced by compute_single_point_map(), get_dphidxi_map(), and init_reference_to_physical_map().

◆ dphidzeta_map

std::vector<std::vector<Real> > libMesh::FEMap::dphidzeta_map
protected

Map for the derivative, d(phi)/d(zeta).

Definition at line 866 of file fe_map.h.

Referenced by compute_single_point_map(), get_dphidzeta_map(), and init_reference_to_physical_map().

◆ dpsideta_map

std::vector<std::vector<Real> > libMesh::FEMap::dpsideta_map
protected

Map for the derivative of the side function, d(psi)/d(eta).

Definition at line 917 of file fe_map.h.

Referenced by libMesh::FEXYZMap::compute_face_map(), compute_face_map(), get_dpsideta(), and init_face_shape_functions().

◆ dpsidxi_map

std::vector<std::vector<Real> > libMesh::FEMap::dpsidxi_map
protected

Map for the derivative of the side functions, d(psi)/d(xi).

Definition at line 911 of file fe_map.h.

Referenced by compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), get_dpsidxi(), init_edge_shape_functions(), and init_face_shape_functions().

◆ dxidx_map

std::vector<Real> libMesh::FEMap::dxidx_map
protected

Map for partial derivatives: d(xi)/d(x).

Needed for the Jacobian.

Definition at line 776 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_dxidx(), and resize_quadrature_map_vectors().

◆ dxidy_map

std::vector<Real> libMesh::FEMap::dxidy_map
protected

Map for partial derivatives: d(xi)/d(y).

Needed for the Jacobian.

Definition at line 782 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_dxidy(), and resize_quadrature_map_vectors().

◆ dxidz_map

std::vector<Real> libMesh::FEMap::dxidz_map
protected

Map for partial derivatives: d(xi)/d(z).

Needed for the Jacobian.

Definition at line 788 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_dxidz(), and resize_quadrature_map_vectors().

◆ dxyzdeta_map

std::vector<RealGradient> libMesh::FEMap::dxyzdeta_map
protected

◆ dxyzdxi_map

std::vector<RealGradient> libMesh::FEMap::dxyzdxi_map
protected

◆ dxyzdzeta_map

std::vector<RealGradient> libMesh::FEMap::dxyzdzeta_map
protected

Vector of partial derivatives: d(x)/d(zeta), d(y)/d(zeta), d(z)/d(zeta)

Definition at line 730 of file fe_map.h.

Referenced by compute_affine_map(), compute_null_map(), compute_single_point_map(), dxdzeta_map(), dydzeta_map(), dzdzeta_map(), get_dxyzdzeta(), and resize_quadrature_map_vectors().

◆ dzetadx_map

std::vector<Real> libMesh::FEMap::dzetadx_map
protected

Map for partial derivatives: d(zeta)/d(x).

Needed for the Jacobian.

Definition at line 814 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_dzetadx(), and resize_quadrature_map_vectors().

◆ dzetady_map

std::vector<Real> libMesh::FEMap::dzetady_map
protected

Map for partial derivatives: d(zeta)/d(y).

Needed for the Jacobian.

Definition at line 820 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_dzetady(), and resize_quadrature_map_vectors().

◆ dzetadz_map

std::vector<Real> libMesh::FEMap::dzetadz_map
protected

Map for partial derivatives: d(zeta)/d(z).

Needed for the Jacobian.

Definition at line 826 of file fe_map.h.

Referenced by compute_affine_map(), compute_inverse_map_second_derivs(), compute_null_map(), compute_single_point_map(), get_dzetadz(), and resize_quadrature_map_vectors().

◆ jac

std::vector<Real> libMesh::FEMap::jac
protected

Jacobian values at quadrature points.

Definition at line 967 of file fe_map.h.

Referenced by compute_affine_map(), compute_null_map(), compute_single_point_map(), get_jacobian(), and resize_quadrature_map_vectors().

◆ jacobian_tolerance

Real libMesh::FEMap::jacobian_tolerance
protected

The Jacobian tolerance used for determining when the mapping fails.

The mapping is determined to fail if jac <= jacobian_tolerance. If not set by the user, this number defaults to 0

Definition at line 1011 of file fe_map.h.

Referenced by compute_single_point_map(), and set_jacobian_tolerance().

◆ JxW

std::vector<Real> libMesh::FEMap::JxW
protected

◆ normals

std::vector<Point> libMesh::FEMap::normals
protected

Normal vectors on boundary at quadrature points.

Definition at line 952 of file fe_map.h.

Referenced by compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), and get_normals().

◆ phi_map

std::vector<std::vector<Real> > libMesh::FEMap::phi_map
protected

Map for the shape function phi.

Definition at line 851 of file fe_map.h.

Referenced by compute_affine_map(), compute_single_point_map(), get_phi_map(), and init_reference_to_physical_map().

◆ psi_map

std::vector<std::vector<Real> > libMesh::FEMap::psi_map
protected

Map for the side shape functions, psi.

Definition at line 905 of file fe_map.h.

Referenced by compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), get_psi(), init_edge_shape_functions(), and init_face_shape_functions().

◆ tangents

std::vector<std::vector<Point> > libMesh::FEMap::tangents
protected

Tangent vectors on boundary at quadrature points.

Definition at line 947 of file fe_map.h.

Referenced by compute_edge_map(), libMesh::FEXYZMap::compute_face_map(), compute_face_map(), and get_tangents().

◆ xyz

std::vector<Point> libMesh::FEMap::xyz
protected

The documentation for this class was generated from the following files:
libMesh::FEMap::d2psidxi2_map
std::vector< std::vector< Real > > d2psidxi2_map
Map for the second derivatives (in xi) of the side shape functions.
Definition: fe_map.h:926
libMesh::FEInterface::shape_function
static shape_ptr shape_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:838
libMesh::FEMap::dzetadx_map
std::vector< Real > dzetadx_map
Map for partial derivatives: d(zeta)/d(x).
Definition: fe_map.h:814
libMesh::FEMap::calculate_d2xyz
bool calculate_d2xyz
Should we calculate mapping hessians?
Definition: fe_map.h:995
libMesh::FEMap::d2xyzdetadzeta_map
std::vector< RealGradient > d2xyzdetadzeta_map
Vector of mixed second partial derivatives in eta-zeta: d^2(x)/d(eta)d(zeta) d^2(y)/d(eta)d(zeta) d^2...
Definition: fe_map.h:762
libMesh::FEMap::dydeta_map
Real dydeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.
Definition: fe_map.h:675
libMesh::FEInterface::n_shape_functions
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:446
libMesh::MeshTools::Subdivision::find_one_ring
void find_one_ring(const Tri3Subdivision *elem, std::vector< const Node * > &nodes)
Determines the 1-ring of element elem, and writes it to the nodes vector.
Definition: mesh_subdivision_support.C:31
libMesh::FEMap::dpsidxi_map
std::vector< std::vector< Real > > dpsidxi_map
Map for the derivative of the side functions, d(psi)/d(xi).
Definition: fe_map.h:911
libMesh::FEMap::map_deriv
static Point map_deriv(const unsigned int dim, const Elem *elem, const unsigned int j, const Point &reference_point)
Definition: fe_map.C:2076
libMesh::FEMap::dzetadz_map
std::vector< Real > dzetadz_map
Map for partial derivatives: d(zeta)/d(z).
Definition: fe_map.h:826
libMesh::FEMap::dxidy_map
std::vector< Real > dxidy_map
Map for partial derivatives: d(xi)/d(y).
Definition: fe_map.h:782
libMesh::FEMap::calculate_xyz
bool calculate_xyz
Should we calculate physical point locations?
Definition: fe_map.h:983
libMesh::FEMap::d2xyzdxi2_map
std::vector< RealGradient > d2xyzdxi2_map
Vector of second partial derivatives in xi: d^2(x)/d(xi)^2, d^2(y)/d(xi)^2, d^2(z)/d(xi)^2.
Definition: fe_map.h:738
libMesh::index_range
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:106
B
Definition: assembly.h:38
libMesh::XYZ
Definition: enum_fe_family.h:46
libMesh::FEMap::jac
std::vector< Real > jac
Jacobian values at quadrature points.
Definition: fe_map.h:967
libMesh::RealTensor
RealTensorValue RealTensor
Definition: hp_coarsentest.h:50
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::FEMap::dxdzeta_map
Real dxdzeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.
Definition: fe_map.h:691
libMesh::FEMap::d2phidxideta_map
std::vector< std::vector< Real > > d2phidxideta_map
Map for the second derivative, d^2(phi)/d(xi)d(eta).
Definition: fe_map.h:878
libMesh::FEMap::JxW
std::vector< Real > JxW
Jacobian*Weight values at quadrature points.
Definition: fe_map.h:972
libMesh::FEMap::dxdxi_map
Real dxdxi_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.
Definition: fe_map.h:643
libMesh::FEInterface::shape_deriv_function
static shape_deriv_ptr shape_deriv_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:916
std::sqrt
MetaPhysicL::DualNumber< T, D > sqrt(const MetaPhysicL::DualNumber< T, D > &in)
libMesh::FEMap::map
static Point map(const unsigned int dim, const Elem *elem, const Point &reference_point)
Definition: fe_map.C:2043
libMesh::FEMap::_elem_nodes
std::vector< const Node * > _elem_nodes
Work vector for compute_affine_map()
Definition: fe_map.h:1023
libMesh::FEMap::dydzeta_map
Real dydzeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.
Definition: fe_map.h:699
libMesh::FEMap::dxyzdeta_map
std::vector< RealGradient > dxyzdeta_map
Vector of partial derivatives: d(x)/d(eta), d(y)/d(eta), d(z)/d(eta)
Definition: fe_map.h:724
libMesh::FEMap::dydxi_map
Real dydxi_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.
Definition: fe_map.h:651
libMesh::FEInterface::shape_ptr
Real(* shape_ptr)(const Elem *elem, const Order o, const unsigned int i, const Point &p, const bool add_p_level)
Definition: fe_interface.h:296
libMesh::RealVectorValue
VectorValue< Real > RealVectorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
Definition: hp_coarsentest.h:46
libMesh::FEMap::normals
std::vector< Point > normals
Normal vectors on boundary at quadrature points.
Definition: fe_map.h:952
libMesh::FEMap::xyz
std::vector< Point > xyz
The spatial locations of the quadrature points.
Definition: fe_map.h:712
libMesh::RATIONAL_BERNSTEIN_MAP
Definition: enum_elem_type.h:84
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
libMesh::FEMap::compute_affine_map
virtual void compute_affine_map(const unsigned int dim, const std::vector< Real > &qw, const Elem *elem)
Compute the jacobian and some other additional data fields.
Definition: fe_map.C:1277
libMesh::FEMap::d2phidxi2_map
std::vector< std::vector< Real > > d2phidxi2_map
Map for the second derivative, d^2(phi)/d(xi)^2.
Definition: fe_map.h:873
libMesh::FEMap::compute_null_map
virtual void compute_null_map(const unsigned int dim, const std::vector< Real > &qw)
Assign a fake jacobian and some other additional data fields.
Definition: fe_map.C:1358
libMesh::FEMap::detady_map
std::vector< Real > detady_map
Map for partial derivatives: d(eta)/d(y).
Definition: fe_map.h:801
libMesh::FEMap::compute_single_point_map
void compute_single_point_map(const unsigned int dim, const std::vector< Real > &qw, const Elem *elem, unsigned int p, const std::vector< const Node * > &elem_nodes, bool compute_second_derivatives)
Compute the jacobian and some other additional data fields at the single point with index p.
Definition: fe_map.C:446
libMesh::FEMap::compute_face_map
virtual void compute_face_map(int dim, const std::vector< Real > &qw, const Elem *side)
Same as compute_map, but for a side.
Definition: fe_boundary.C:587
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::FEMap::determine_calculations
void determine_calculations()
Determine which values are to be calculated.
Definition: fe_map.h:621
libMesh::InfFEMap::inverse_map
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: inf_fe_map.C:88
libMesh::FEMap::d2xyzdxidzeta_map
std::vector< RealGradient > d2xyzdxidzeta_map
Vector of second partial derivatives in xi-zeta: d^2(x)/d(xi)d(zeta), d^2(y)/d(xi)d(zeta),...
Definition: fe_map.h:756
libMesh::FEMap::inverse_map
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_map.C:1622
libMesh::FEMap::resize_quadrature_map_vectors
void resize_quadrature_map_vectors(const unsigned int dim, unsigned int n_qp)
A utility function for use by compute_*_map.
Definition: fe_map.C:1196
libMesh::FEMap::dpsideta_map
std::vector< std::vector< Real > > dpsideta_map
Map for the derivative of the side function, d(psi)/d(eta).
Definition: fe_map.h:917
libMesh::FEMap::detadz_map
std::vector< Real > detadz_map
Map for partial derivatives: d(eta)/d(z).
Definition: fe_map.h:807
libMesh::FEMap::tangents
std::vector< std::vector< Point > > tangents
Tangent vectors on boundary at quadrature points.
Definition: fe_map.h:947
libMesh::RealTensorValue
TensorValue< Real > RealTensorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
Definition: hp_coarsentest.h:48
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:526
libMesh::FEMap::psi_map
std::vector< std::vector< Real > > psi_map
Map for the side shape functions, psi.
Definition: fe_map.h:905
libMesh::FEMap::d2phidzeta2_map
std::vector< std::vector< Real > > d2phidzeta2_map
Map for the second derivative, d^2(phi)/d(zeta)^2.
Definition: fe_map.h:898
libMesh::FEMap::calculate_dxyz
bool calculate_dxyz
Should we calculate mapping gradients?
Definition: fe_map.h:988
libMesh::FEMap::d2xyzdzeta2_map
std::vector< RealGradient > d2xyzdzeta2_map
Vector of second partial derivatives in zeta: d^2(x)/d(zeta)^2.
Definition: fe_map.h:768
libMesh::FEMap::dphideta_map
std::vector< std::vector< Real > > dphideta_map
Map for the derivative, d(phi)/d(eta).
Definition: fe_map.h:861
libMesh::FEMap::d2xidxyz2_map
std::vector< std::vector< Real > > d2xidxyz2_map
Second derivatives of "xi" reference coordinate wrt physical coordinates.
Definition: fe_map.h:833
libMesh::FEMap::d2psideta2_map
std::vector< std::vector< Real > > d2psideta2_map
Map for the second derivatives (in eta) of the side shape functions.
Definition: fe_map.h:940
libMesh::FEMap::dzdeta_map
Real dzdeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.
Definition: fe_map.h:683
libMesh::FEMap::map_fe_type
static FEFamily map_fe_type(const Elem &elem)
Definition: fe_map.C:47
libMesh::RATIONAL_BERNSTEIN
Definition: enum_fe_family.h:64
libMesh::FEMap::d2etadxyz2_map
std::vector< std::vector< Real > > d2etadxyz2_map
Second derivatives of "eta" reference coordinate wrt physical coordinates.
Definition: fe_map.h:839
n_nodes
const dof_id_type n_nodes
Definition: tecplot_io.C:68
A
static PetscErrorCode Mat * A
Definition: petscdmlibmeshimpl.C:1026
libMesh::FEInterface::shape_deriv_ptr
Real(* shape_deriv_ptr)(const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Definition: fe_interface.h:345
libMesh::FEMap::detadx_map
std::vector< Real > detadx_map
Map for partial derivatives: d(eta)/d(x).
Definition: fe_map.h:795
libMesh::FEMap::jacobian_tolerance
Real jacobian_tolerance
The Jacobian tolerance used for determining when the mapping fails.
Definition: fe_map.h:1011
libMesh::FEAbstract::on_reference_element
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_abstract.C:606
libMesh::LAGRANGE_MAP
Definition: enum_elem_type.h:83
libMesh::FEMap::dphidxi_map
std::vector< std::vector< Real > > dphidxi_map
Map for the derivative, d(phi)/d(xi).
Definition: fe_map.h:856
libMesh::FEMap::compute_inverse_map_second_derivs
void compute_inverse_map_second_derivs(unsigned p)
A helper function used by FEMap::compute_single_point_map() to compute second derivatives of the inve...
Definition: fe_map.C:1502
libMesh::InfFEMap::map
static Point map(const unsigned int dim, const Elem *inf_elem, const Point &reference_point)
Definition: inf_fe_map.C:40
libMesh::FEMap::d2phidetadzeta_map
std::vector< std::vector< Real > > d2phidetadzeta_map
Map for the second derivative, d^2(phi)/d(eta)d(zeta).
Definition: fe_map.h:893
libMesh::FEFamily
FEFamily
Definition: enum_fe_family.h:34
libMesh::FEMap::phi_map
std::vector< std::vector< Real > > phi_map
Map for the shape function phi.
Definition: fe_map.h:851
libMesh::FEMap::dxyzdzeta_map
std::vector< RealGradient > dxyzdzeta_map
Vector of partial derivatives: d(x)/d(zeta), d(y)/d(zeta), d(z)/d(zeta)
Definition: fe_map.h:730
libMesh::FEInterface::shape_second_deriv_ptr
Real(* shape_second_deriv_ptr)(const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Definition: fe_interface.h:409
libMesh::FEMap::dzetady_map
std::vector< Real > dzetady_map
Map for partial derivatives: d(zeta)/d(y).
Definition: fe_map.h:820
libMesh::FEMap::dphidzeta_map
std::vector< std::vector< Real > > dphidzeta_map
Map for the derivative, d(phi)/d(zeta).
Definition: fe_map.h:866
libMesh::FEMap::d2phideta2_map
std::vector< std::vector< Real > > d2phideta2_map
Map for the second derivative, d^2(phi)/d(eta)^2.
Definition: fe_map.h:888
libMesh::FEMap::d2zetadxyz2_map
std::vector< std::vector< Real > > d2zetadxyz2_map
Second derivatives of "zeta" reference coordinate wrt physical coordinates.
Definition: fe_map.h:845
libMesh::FEMap::d2phidxidzeta_map
std::vector< std::vector< Real > > d2phidxidzeta_map
Map for the second derivative, d^2(phi)/d(xi)d(zeta).
Definition: fe_map.h:883
libMesh::err
OStreamProxy err
libMesh::LAGRANGE
Definition: enum_fe_family.h:36
libMesh::FEMap::dzdxi_map
Real dzdxi_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.
Definition: fe_map.h:659
libMesh::FEMap::dzdzeta_map
Real dzdzeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.
Definition: fe_map.h:707
libMesh::TRI3SUBDIVISION
Definition: enum_elem_type.h:69
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::FEMap::dxdeta_map
Real dxdeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected.
Definition: fe_map.h:667
libMesh::FEMap::dxyzdxi_map
std::vector< RealGradient > dxyzdxi_map
Vector of partial derivatives: d(x)/d(xi), d(y)/d(xi), d(z)/d(xi)
Definition: fe_map.h:718
libMesh::FEMap::dxidx_map
std::vector< Real > dxidx_map
Map for partial derivatives: d(xi)/d(x).
Definition: fe_map.h:776
libMesh::FEInterface::shape_second_deriv_function
static shape_second_deriv_ptr shape_second_deriv_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:989
libMesh::FEMap::dxidz_map
std::vector< Real > dxidz_map
Map for partial derivatives: d(xi)/d(z).
Definition: fe_map.h:788
libMesh::FEMap::calculations_started
bool calculations_started
Have calculations with this object already been started? Then all get_* functions should already have...
Definition: fe_map.h:978
libMesh::FEMap::d2xyzdxideta_map
std::vector< RealGradient > d2xyzdxideta_map
Vector of mixed second partial derivatives in xi-eta: d^2(x)/d(xi)d(eta) d^2(y)/d(xi)d(eta) d^2(z)/d(...
Definition: fe_map.h:744
libMesh::FEMap::d2xyzdeta2_map
std::vector< RealGradient > d2xyzdeta2_map
Vector of second partial derivatives in eta: d^2(x)/d(eta)^2.
Definition: fe_map.h:750
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::FEMap::d2psidxideta_map
std::vector< std::vector< Real > > d2psidxideta_map
Map for the second (cross) derivatives in xi, eta of the side shape functions.
Definition: fe_map.h:933
libMesh::FEMap::curvatures
std::vector< Real > curvatures
The mean curvature (= one half the sum of the principal curvatures) on the boundary at the quadrature...
Definition: fe_map.h:960