libMesh::InfFE< Dim, T_radial, T_map > Class Template Reference

Base class for all the infinite geometric element types. More...

#include <fe.h>

Inheritance diagram for libMesh::InfFE< Dim, T_radial, T_map >:

Classes

class  Base
 
class  Radial
 

Public Types

typedef OutputType OutputShape
 
typedef TensorTools::IncrementRank< OutputShape >::type OutputGradient
 
typedef TensorTools::IncrementRank< OutputGradient >::type OutputTensor
 
typedef TensorTools::DecrementRank< OutputShape >::type OutputDivergence
 
typedef TensorTools::MakeNumber< OutputShape >::type OutputNumber
 
typedef TensorTools::IncrementRank< OutputNumber >::type OutputNumberGradient
 
typedef TensorTools::IncrementRank< OutputNumberGradient >::type OutputNumberTensor
 
typedef TensorTools::DecrementRank< OutputNumber >::type OutputNumberDivergence
 

Public Member Functions

 InfFE (const FEType &fet)
 
 ~InfFE ()
 
virtual FEContinuity get_continuity () const libmesh_override
 
virtual bool is_hierarchic () const libmesh_override
 
virtual void reinit (const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
 
virtual void reinit (const Elem *elem, const unsigned int side, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
 
virtual void edge_reinit (const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
 
virtual void side_map (const Elem *, const Elem *, const unsigned int, const std::vector< Point > &, std::vector< Point > &) libmesh_override
 
virtual void attach_quadrature_rule (QBase *q) libmesh_override
 
virtual unsigned int n_shape_functions () const libmesh_override
 
virtual unsigned int n_quadrature_points () const libmesh_override
 
template<>
std::unique_ptr< FEGenericBase< Real > > build (const unsigned int dim, const FEType &fet)
 
template<>
std::unique_ptr< FEGenericBase< RealGradient > > build (const unsigned int dim, const FEType &fet)
 
template<>
std::unique_ptr< FEGenericBase< Real > > build_InfFE (const unsigned int dim, const FEType &fet)
 
template<>
std::unique_ptr< FEGenericBase< RealGradient > > build_InfFE (const unsigned int, const FEType &)
 
const std::vector< std::vector< OutputShape > > & get_phi () const
 
const std::vector< std::vector< OutputGradient > > & get_dphi () const
 
const std::vector< std::vector< OutputShape > > & get_curl_phi () const
 
const std::vector< std::vector< OutputDivergence > > & get_div_phi () const
 
const std::vector< std::vector< OutputShape > > & get_dphidx () const
 
const std::vector< std::vector< OutputShape > > & get_dphidy () const
 
const std::vector< std::vector< OutputShape > > & get_dphidz () const
 
const std::vector< std::vector< OutputShape > > & get_dphidxi () const
 
const std::vector< std::vector< OutputShape > > & get_dphideta () const
 
const std::vector< std::vector< OutputShape > > & get_dphidzeta () const
 
const std::vector< std::vector< OutputTensor > > & get_d2phi () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidx2 () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidxdy () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidxdz () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidy2 () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidydz () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidz2 () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidxi2 () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidxideta () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidxidzeta () const
 
const std::vector< std::vector< OutputShape > > & get_d2phideta2 () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidetadzeta () const
 
const std::vector< std::vector< OutputShape > > & get_d2phidzeta2 () const
 
const std::vector< OutputGradient > & get_dphase () const
 
const std::vector< Real > & get_Sobolev_weight () const
 
const std::vector< RealGradient > & get_Sobolev_dweight () const
 
void print_phi (std::ostream &os) const
 
void print_dphi (std::ostream &os) const
 
void print_d2phi (std::ostream &os) const
 
unsigned int get_dim () const
 
const std::vector< Point > & get_xyz () 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< Point > > & get_tangents () const
 
const std::vector< Point > & get_normals () const
 
const std::vector< Real > & get_curvatures () const
 
ElemType get_type () const
 
unsigned int get_p_level () const
 
FEType get_fe_type () const
 
Order get_order () const
 
void set_fe_order (int new_order)
 
FEFamily get_family () const
 
const FEMapget_fe_map () const
 
void print_JxW (std::ostream &os) const
 
void print_xyz (std::ostream &os) const
 
void print_info (std::ostream &os) const
 

Static Public Member Functions

static Real shape (const FEType &fet, const ElemType t, const unsigned int i, const Point &p)
 
static Real shape (const FEType &fet, const Elem *elem, const unsigned int i, const Point &p)
 
static void compute_data (const FEType &fe_t, const Elem *inf_elem, FEComputeData &data)
 
static unsigned int n_shape_functions (const FEType &fet, const ElemType t)
 
static unsigned int n_dofs (const FEType &fet, const ElemType inf_elem_type)
 
static unsigned int n_dofs_at_node (const FEType &fet, const ElemType inf_elem_type, const unsigned int n)
 
static unsigned int n_dofs_per_elem (const FEType &fet, const ElemType inf_elem_type)
 
static void nodal_soln (const FEType &fet, const Elem *elem, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
 
static Point map (const Elem *inf_elem, const Point &reference_point)
 
static Point inverse_map (const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
 
static void inverse_map (const Elem *elem, const std::vector< Point > &physical_points, std::vector< Point > &reference_points, const Real tolerance=TOLERANCE, const bool secure=true)
 
static std::unique_ptr< FEGenericBasebuild (const unsigned int dim, const FEType &type)
 
static std::unique_ptr< FEGenericBasebuild_InfFE (const unsigned int dim, const FEType &type)
 
static void compute_proj_constraints (DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
 
static void coarsened_dof_values (const NumericVector< Number > &global_vector, const DofMap &dof_map, const Elem *coarse_elem, DenseVector< Number > &coarse_dofs, const unsigned int var, const bool use_old_dof_indices=false)
 
static void coarsened_dof_values (const NumericVector< Number > &global_vector, const DofMap &dof_map, const Elem *coarse_elem, DenseVector< Number > &coarse_dofs, const bool use_old_dof_indices=false)
 
static void compute_periodic_constraints (DofConstraints &constraints, DofMap &dof_map, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const unsigned int variable_number, const Elem *elem)
 
static bool on_reference_element (const Point &p, const ElemType t, const Real eps=TOLERANCE)
 
static void get_refspace_nodes (const ElemType t, std::vector< Point > &nodes)
 
static void compute_node_constraints (NodeConstraints &constraints, const Elem *elem)
 
static void compute_periodic_node_constraints (NodeConstraints &constraints, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const Elem *elem)
 
static void print_info (std::ostream &out=libMesh::out)
 
static std::string get_info ()
 
static unsigned int n_objects ()
 
static void enable_print_counter_info ()
 
static void disable_print_counter_info ()
 

Protected Types

typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 

Protected Member Functions

void update_base_elem (const Elem *inf_elem)
 
virtual void init_base_shape_functions (const std::vector< Point > &, const Elem *) libmesh_override
 
void init_radial_shape_functions (const Elem *inf_elem, const std::vector< Point > *radial_pts=libmesh_nullptr)
 
void init_shape_functions (const std::vector< Point > &radial_qp, const std::vector< Point > &base_qp, const Elem *inf_elem)
 
void init_face_shape_functions (const std::vector< Point > &qp, const Elem *side)
 
void combine_base_radial (const Elem *inf_elem)
 
virtual void compute_shape_functions (const Elem *, const std::vector< Point > &) libmesh_override
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order o, unsigned i)
 
template<>
Real eval (Real v, Order o, unsigned i)
 
template<>
Real eval (Real v, Order o, unsigned i)
 
template<>
Real eval_deriv (Real v, Order o, unsigned i)
 
template<>
Real eval_deriv (Real v, Order o, unsigned i)
 
template<>
Real eval_deriv (Real v, Order o, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
template<>
Real eval_deriv (Real v, Order, unsigned i)
 
void determine_calculations ()
 
void increment_constructor_count (const std::string &name)
 
void increment_destructor_count (const std::string &name)
 

Static Protected Member Functions

static Real eval (Real v, Order o_radial, unsigned int i)
 
static Real eval_deriv (Real v, Order o_radial, unsigned int i)
 
static void compute_node_indices (const ElemType inf_elem_type, const unsigned int outer_node_index, unsigned int &base_node, unsigned int &radial_node)
 
static void compute_node_indices_fast (const ElemType inf_elem_type, const unsigned int outer_node_index, unsigned int &base_node, unsigned int &radial_node)
 
static void compute_shape_indices (const FEType &fet, const ElemType inf_elem_type, const unsigned int i, unsigned int &base_shape, unsigned int &radial_shape)
 

Protected Attributes

std::vector< Realdist
 
std::vector< Realdweightdv
 
std::vector< Realsom
 
std::vector< Realdsomdv
 
std::vector< std::vector< Real > > mode
 
std::vector< std::vector< Real > > dmodedv
 
std::vector< std::vector< Real > > radial_map
 
std::vector< std::vector< Real > > dradialdv_map
 
std::vector< Realdphasedxi
 
std::vector< Realdphasedeta
 
std::vector< Realdphasedzeta
 
std::vector< unsigned int > _radial_node_index
 
std::vector< unsigned int > _base_node_index
 
std::vector< unsigned int > _radial_shape_index
 
std::vector< unsigned int > _base_shape_index
 
unsigned int _n_total_approx_sf
 
unsigned int _n_total_qp
 
std::vector< Real_total_qrule_weights
 
std::unique_ptr< QBasebase_qrule
 
std::unique_ptr< QBaseradial_qrule
 
std::unique_ptr< Elembase_elem
 
std::unique_ptr< FEBasebase_fe
 
FEType current_fe_type
 
std::unique_ptr< FETransformationBase< OutputType > > _fe_trans
 
std::vector< std::vector< OutputShape > > phi
 
std::vector< std::vector< OutputGradient > > dphi
 
std::vector< std::vector< OutputShape > > curl_phi
 
std::vector< std::vector< OutputDivergence > > div_phi
 
std::vector< std::vector< OutputShape > > dphidxi
 
std::vector< std::vector< OutputShape > > dphideta
 
std::vector< std::vector< OutputShape > > dphidzeta
 
std::vector< std::vector< OutputShape > > dphidx
 
std::vector< std::vector< OutputShape > > dphidy
 
std::vector< std::vector< OutputShape > > dphidz
 
std::vector< std::vector< OutputTensor > > d2phi
 
std::vector< std::vector< OutputShape > > d2phidxi2
 
std::vector< std::vector< OutputShape > > d2phidxideta
 
std::vector< std::vector< OutputShape > > d2phidxidzeta
 
std::vector< std::vector< OutputShape > > d2phideta2
 
std::vector< std::vector< OutputShape > > d2phidetadzeta
 
std::vector< std::vector< OutputShape > > d2phidzeta2
 
std::vector< std::vector< OutputShape > > d2phidx2
 
std::vector< std::vector< OutputShape > > d2phidxdy
 
std::vector< std::vector< OutputShape > > d2phidxdz
 
std::vector< std::vector< OutputShape > > d2phidy2
 
std::vector< std::vector< OutputShape > > d2phidydz
 
std::vector< std::vector< OutputShape > > d2phidz2
 
std::vector< OutputGradientdphase
 
std::vector< RealGradientdweight
 
std::vector< Realweight
 
std::unique_ptr< FEMap_fe_map
 
const unsigned int dim
 
bool calculations_started
 
bool calculate_phi
 
bool calculate_dphi
 
bool calculate_d2phi
 
bool calculate_curl_phi
 
bool calculate_div_phi
 
bool calculate_dphiref
 
FEType fe_type
 
ElemType elem_type
 
unsigned int _p_level
 
QBaseqrule
 
bool shapes_on_quadrature
 

Static Protected Attributes

static Counts _counts
 
static Threads::atomic< unsigned int > _n_objects
 
static Threads::spin_mutex _mutex
 
static bool _enable_print_counter = true
 

Private Member Functions

virtual bool shapes_need_reinit () const libmesh_override
 

Static Private Attributes

static ElemType _compute_node_indices_fast_current_elem_type = INVALID_ELEM
 
static bool _warned_for_nodal_soln = false
 
static bool _warned_for_shape = false
 

Friends

template<unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
class InfFE
 

Detailed Description

template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
class libMesh::InfFE< Dim, T_radial, T_map >

Base class for all the infinite geometric element types.

A specific instantiation of the FEBase class. This class is templated, and specific template instantiations will result in different Infinite Element families, similar to the FE class. InfFE builds a FE<Dim-1,T_base>, and most of the requests related to the base are handed over to this object. All methods related to the radial part are collected in the nested class Radial. Similarly, most of the static methods concerning base approximation are contained in Base.

Having different shape approximation families in radial direction introduces the requirement for an additional Order in this class. Therefore, the FEType internals change when infinite elements are enabled. When the specific infinite element type is not known at compile time, use the FEBase::build() member to create abstract (but still optimized) infinite elements at run time.

The node numbering scheme is the one from the current infinite element. Each node in the base holds exactly the same number of dofs as an adjacent conventional FE would contain. The nodes further out hold the additional dof necessary for radial approximation. The order of the outer nodes' components is such that the radial shapes have highest priority, followed by the base shapes.

Author
Daniel Dreyer
Date
2003

Definition at line 40 of file fe.h.

Member Typedef Documentation

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts
protectedinherited

Data structure to log the information. The log is identified by the class name.

Definition at line 119 of file reference_counter.h.

template<typename OutputType>
typedef TensorTools::DecrementRank<OutputShape>::type libMesh::FEGenericBase< OutputType >::OutputDivergence
inherited

Definition at line 123 of file fe_base.h.

template<typename OutputType>
typedef TensorTools::IncrementRank<OutputShape>::type libMesh::FEGenericBase< OutputType >::OutputGradient
inherited

Definition at line 121 of file fe_base.h.

template<typename OutputType>
typedef TensorTools::MakeNumber<OutputShape>::type libMesh::FEGenericBase< OutputType >::OutputNumber
inherited

Definition at line 124 of file fe_base.h.

template<typename OutputType>
typedef TensorTools::DecrementRank<OutputNumber>::type libMesh::FEGenericBase< OutputType >::OutputNumberDivergence
inherited

Definition at line 127 of file fe_base.h.

template<typename OutputType>
typedef TensorTools::IncrementRank<OutputNumber>::type libMesh::FEGenericBase< OutputType >::OutputNumberGradient
inherited

Definition at line 125 of file fe_base.h.

template<typename OutputType>
typedef TensorTools::IncrementRank<OutputNumberGradient>::type libMesh::FEGenericBase< OutputType >::OutputNumberTensor
inherited

Definition at line 126 of file fe_base.h.

template<typename OutputType>
typedef OutputType libMesh::FEGenericBase< OutputType >::OutputShape
inherited

Convenient typedefs for gradients of output, hessians of output, and potentially-complex-valued versions of same.

Definition at line 120 of file fe_base.h.

template<typename OutputType>
typedef TensorTools::IncrementRank<OutputGradient>::type libMesh::FEGenericBase< OutputType >::OutputTensor
inherited

Definition at line 122 of file fe_base.h.

Constructor & Destructor Documentation

template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
libMesh::InfFE< Dim, T_radial, T_map >::InfFE ( const FEType fet)
explicit

Constructor and empty destructor. Initializes some data structures. Builds a FE<Dim-1,T_base> object to handle approximation in the base, so that there is no need to template InfFE<Dim,T_radial,T_map> also with respect to the base approximation T_base.

The same remarks concerning compile-time optimization for FE also hold for InfFE. Use the FEBase::build_InfFE(const unsigned int, const FEType &) method to build specific instantiations of InfFE at run time.

Definition at line 35 of file inf_fe.C.

References libMesh::InfFE< Dim, T_radial, T_map >::base_fe, libMesh::FEGenericBase< OutputType >::build(), libMesh::FEAbstract::fe_type, libMesh::FEType::inf_map, and libMesh::FEType::radial_family.

35  :
36  FEBase (Dim, fet),
37 
39  _n_total_qp (0),
40 
41  // initialize the current_fe_type to all the same
42  // values as \p fet (since the FE families and coordinate
43  // map type should not change), but use an invalid order
44  // for the radial part (since this is the only order
45  // that may change!).
46  // the data structures like \p phi etc are not initialized
47  // through the constructor, but through reinit()
48  current_fe_type (FEType(fet.order,
49  fet.family,
51  fet.radial_family,
52  fet.inf_map))
53 
54 {
55  // Sanity checks
56  libmesh_assert_equal_to (T_radial, fe_type.radial_family);
57  libmesh_assert_equal_to (T_map, fe_type.inf_map);
58 
59  // build the base_fe object
60  if (Dim != 1)
61  base_fe = FEBase::build(Dim-1, fet);
62 }
unsigned int _n_total_qp
Definition: inf_fe.h:739
std::unique_ptr< FEBase > base_fe
Definition: inf_fe.h:771
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
InfMapType inf_map
Definition: fe_type.h:258
FEType current_fe_type
Definition: inf_fe.h:781
FEGenericBase< Real > FEBase
unsigned int _n_total_approx_sf
Definition: inf_fe.h:733
FEFamily radial_family
Definition: fe_type.h:250
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
libMesh::InfFE< Dim, T_radial, T_map >::~InfFE ( )
inline

Member Function Documentation

template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void libMesh::InfFE< Dim, T_radial, T_map >::attach_quadrature_rule ( QBase q)
virtual

The use of quadrature rules with the InfFE class is somewhat different from the approach of the FE class. While the FE class requires an appropriately initialized quadrature rule object, and simply uses it, the InfFE class requires only the quadrature rule object of the current FE class. From this QBase *, it determines the necessary data, and builds two appropriate quadrature classes, one for radial, and another for base integration, using the convenient QBase::build() method.

Implements libMesh::FEAbstract.

Definition at line 67 of file inf_fe.C.

References libMesh::InfFE< Dim, T_radial, T_map >::base_fe, libMesh::InfFE< Dim, T_radial, T_map >::base_qrule, libMesh::QBase::build(), libMesh::FEAbstract::fe_type, libMesh::QBase::get_dim(), libMesh::OrderWrapper::get_order(), libMesh::QBase::get_order(), libMesh::FEAbstract::qrule, libMesh::FEType::radial_order, libMesh::InfFE< Dim, T_radial, T_map >::radial_qrule, and libMesh::QBase::type().

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::side_map().

68 {
69  libmesh_assert(q);
70  libmesh_assert(base_fe);
71 
72  const Order base_int_order = q->get_order();
73  const Order radial_int_order = static_cast<Order>(2 * (static_cast<unsigned int>(fe_type.radial_order.get_order()) + 1) +2);
74  const unsigned int qrule_dim = q->get_dim();
75 
76  if (Dim != 1)
77  {
78  // build a Dim-1 quadrature rule of the type that we received
79  base_qrule = QBase::build(q->type(), qrule_dim-1, base_int_order);
80  base_fe->attach_quadrature_rule(base_qrule.get());
81  }
82 
83  // in radial direction, always use Gauss quadrature
84  radial_qrule.reset(new QGauss(1, radial_int_order));
85 
86  // currently not used. But maybe helpful to store the QBase *
87  // with which we initialized our own quadrature rules
88  qrule = q;
89 }
int get_order() const
Definition: fe_type.h:78
std::unique_ptr< QBase > radial_qrule
Definition: inf_fe.h:757
OrderWrapper radial_order
Definition: fe_type.h:237
std::unique_ptr< FEBase > base_fe
Definition: inf_fe.h:771
static std::unique_ptr< QBase > build(const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
std::unique_ptr< QBase > base_qrule
Definition: inf_fe.h:751
template<>
std::unique_ptr< FEGenericBase< Real > > libMesh::FEGenericBase< Real >::build ( const unsigned int  dim,
const FEType fet 
)
inherited

Definition at line 185 of file fe_base.C.

References libMesh::BERNSTEIN, libMesh::CLOUGH, libMesh::FEType::family, libMesh::HERMITE, libMesh::HIERARCHIC, libMesh::L2_HIERARCHIC, libMesh::L2_LAGRANGE, libMesh::LAGRANGE, libMesh::MONOMIAL, libMesh::SCALAR, libMesh::SUBDIVISION, libMesh::SZABAB, and libMesh::XYZ.

187 {
188  switch (dim)
189  {
190  // 0D
191  case 0:
192  {
193  switch (fet.family)
194  {
195  case CLOUGH:
196  return libmesh_make_unique<FE<0,CLOUGH>>(fet);
197 
198  case HERMITE:
199  return libmesh_make_unique<FE<0,HERMITE>>(fet);
200 
201  case LAGRANGE:
202  return libmesh_make_unique<FE<0,LAGRANGE>>(fet);
203 
204  case L2_LAGRANGE:
205  return libmesh_make_unique<FE<0,L2_LAGRANGE>>(fet);
206 
207  case HIERARCHIC:
208  return libmesh_make_unique<FE<0,HIERARCHIC>>(fet);
209 
210  case L2_HIERARCHIC:
211  return libmesh_make_unique<FE<0,L2_HIERARCHIC>>(fet);
212 
213  case MONOMIAL:
214  return libmesh_make_unique<FE<0,MONOMIAL>>(fet);
215 
216 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
217  case SZABAB:
218  return libmesh_make_unique<FE<0,SZABAB>>(fet);
219 
220  case BERNSTEIN:
221  return libmesh_make_unique<FE<0,BERNSTEIN>>(fet);
222 #endif
223 
224  case XYZ:
225  return libmesh_make_unique<FEXYZ<0>>(fet);
226 
227  case SCALAR:
228  return libmesh_make_unique<FEScalar<0>>(fet);
229 
230  default:
231  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
232  }
233  }
234  // 1D
235  case 1:
236  {
237  switch (fet.family)
238  {
239  case CLOUGH:
240  return libmesh_make_unique<FE<1,CLOUGH>>(fet);
241 
242  case HERMITE:
243  return libmesh_make_unique<FE<1,HERMITE>>(fet);
244 
245  case LAGRANGE:
246  return libmesh_make_unique<FE<1,LAGRANGE>>(fet);
247 
248  case L2_LAGRANGE:
249  return libmesh_make_unique<FE<1,L2_LAGRANGE>>(fet);
250 
251  case HIERARCHIC:
252  return libmesh_make_unique<FE<1,HIERARCHIC>>(fet);
253 
254  case L2_HIERARCHIC:
255  return libmesh_make_unique<FE<1,L2_HIERARCHIC>>(fet);
256 
257  case MONOMIAL:
258  return libmesh_make_unique<FE<1,MONOMIAL>>(fet);
259 
260 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
261  case SZABAB:
262  return libmesh_make_unique<FE<1,SZABAB>>(fet);
263 
264  case BERNSTEIN:
265  return libmesh_make_unique<FE<1,BERNSTEIN>>(fet);
266 #endif
267 
268  case XYZ:
269  return libmesh_make_unique<FEXYZ<1>>(fet);
270 
271  case SCALAR:
272  return libmesh_make_unique<FEScalar<1>>(fet);
273 
274  default:
275  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
276  }
277  }
278 
279 
280  // 2D
281  case 2:
282  {
283  switch (fet.family)
284  {
285  case CLOUGH:
286  return libmesh_make_unique<FE<2,CLOUGH>>(fet);
287 
288  case HERMITE:
289  return libmesh_make_unique<FE<2,HERMITE>>(fet);
290 
291  case LAGRANGE:
292  return libmesh_make_unique<FE<2,LAGRANGE>>(fet);
293 
294  case L2_LAGRANGE:
295  return libmesh_make_unique<FE<2,L2_LAGRANGE>>(fet);
296 
297  case HIERARCHIC:
298  return libmesh_make_unique<FE<2,HIERARCHIC>>(fet);
299 
300  case L2_HIERARCHIC:
301  return libmesh_make_unique<FE<2,L2_HIERARCHIC>>(fet);
302 
303  case MONOMIAL:
304  return libmesh_make_unique<FE<2,MONOMIAL>>(fet);
305 
306 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
307  case SZABAB:
308  return libmesh_make_unique<FE<2,SZABAB>>(fet);
309 
310  case BERNSTEIN:
311  return libmesh_make_unique<FE<2,BERNSTEIN>>(fet);
312 #endif
313 
314  case XYZ:
315  return libmesh_make_unique<FEXYZ<2>>(fet);
316 
317  case SCALAR:
318  return libmesh_make_unique<FEScalar<2>>(fet);
319 
320  case SUBDIVISION:
321  return libmesh_make_unique<FESubdivision>(fet);
322 
323  default:
324  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
325  }
326  }
327 
328 
329  // 3D
330  case 3:
331  {
332  switch (fet.family)
333  {
334  case CLOUGH:
335  libmesh_error_msg("ERROR: Clough-Tocher elements currently only support 1D and 2D");
336 
337  case HERMITE:
338  return libmesh_make_unique<FE<3,HERMITE>>(fet);
339 
340  case LAGRANGE:
341  return libmesh_make_unique<FE<3,LAGRANGE>>(fet);
342 
343  case L2_LAGRANGE:
344  return libmesh_make_unique<FE<3,L2_LAGRANGE>>(fet);
345 
346  case HIERARCHIC:
347  return libmesh_make_unique<FE<3,HIERARCHIC>>(fet);
348 
349  case L2_HIERARCHIC:
350  return libmesh_make_unique<FE<3,L2_HIERARCHIC>>(fet);
351 
352  case MONOMIAL:
353  return libmesh_make_unique<FE<3,MONOMIAL>>(fet);
354 
355 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
356  case SZABAB:
357  return libmesh_make_unique<FE<3,SZABAB>>(fet);
358 
359  case BERNSTEIN:
360  return libmesh_make_unique<FE<3,BERNSTEIN>>(fet);
361 #endif
362 
363  case XYZ:
364  return libmesh_make_unique<FEXYZ<3>>(fet);
365 
366  case SCALAR:
367  return libmesh_make_unique<FEScalar<3>>(fet);
368 
369  default:
370  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
371  }
372  }
373 
374  default:
375  libmesh_error_msg("Invalid dimension dim = " << dim);
376  }
377 }
FEFamily family
Definition: fe_type.h:204
const unsigned int dim
Definition: fe_abstract.h:523
template<>
std::unique_ptr< FEGenericBase< RealGradient > > libMesh::FEGenericBase< RealGradient >::build ( const unsigned int  dim,
const FEType fet 
)
inherited

Definition at line 383 of file fe_base.C.

References libMesh::FEType::family, libMesh::LAGRANGE_VEC, and libMesh::NEDELEC_ONE.

385 {
386  switch (dim)
387  {
388  // 0D
389  case 0:
390  {
391  switch (fet.family)
392  {
393  case LAGRANGE_VEC:
394  return libmesh_make_unique<FELagrangeVec<0>>(fet);
395 
396  default:
397  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
398  }
399  }
400  case 1:
401  {
402  switch (fet.family)
403  {
404  case LAGRANGE_VEC:
405  return libmesh_make_unique<FELagrangeVec<1>>(fet);
406 
407  default:
408  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
409  }
410  }
411  case 2:
412  {
413  switch (fet.family)
414  {
415  case LAGRANGE_VEC:
416  return libmesh_make_unique<FELagrangeVec<2>>(fet);
417 
418  case NEDELEC_ONE:
419  return libmesh_make_unique<FENedelecOne<2>>(fet);
420 
421  default:
422  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
423  }
424  }
425  case 3:
426  {
427  switch (fet.family)
428  {
429  case LAGRANGE_VEC:
430  return libmesh_make_unique<FELagrangeVec<3>>(fet);
431 
432  case NEDELEC_ONE:
433  return libmesh_make_unique<FENedelecOne<3>>(fet);
434 
435  default:
436  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
437  }
438  }
439 
440  default:
441  libmesh_error_msg("Invalid dimension dim = " << dim);
442  } // switch(dim)
443 }
FEFamily family
Definition: fe_type.h:204
const unsigned int dim
Definition: fe_abstract.h:523
template<typename OutputType>
static std::unique_ptr<FEGenericBase> libMesh::FEGenericBase< OutputType >::build_InfFE ( const unsigned int  dim,
const FEType type 
)
staticinherited

Builds a specific infinite element type. A std::unique_ptr<FEGenericBase> is returned to prevent a memory leak. This way the user need not remember to delete the object.

The build call will fail if the OutputShape of this class is not compatible with the output required for the requested type

Referenced by libMesh::FEMContext::cached_fe().

template<>
std::unique_ptr< FEGenericBase< Real > > libMesh::FEGenericBase< Real >::build_InfFE ( const unsigned int  dim,
const FEType fet 
)
inherited

Definition at line 456 of file fe_base.C.

References libMesh::CARTESIAN, libMesh::FEType::inf_map, libMesh::INFINITE_MAP, libMesh::JACOBI_20_00, libMesh::JACOBI_30_00, libMesh::LAGRANGE, libMesh::LEGENDRE, and libMesh::FEType::radial_family.

458 {
459  switch (dim)
460  {
461 
462  // 1D
463  case 1:
464  {
465  switch (fet.radial_family)
466  {
467  case INFINITE_MAP:
468  libmesh_error_msg("ERROR: Can't build an infinite element with FEFamily = " << fet.radial_family);
469 
470  case JACOBI_20_00:
471  {
472  switch (fet.inf_map)
473  {
474  case CARTESIAN:
475  return libmesh_make_unique<InfFE<1,JACOBI_20_00,CARTESIAN>>(fet);
476 
477  default:
478  libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << fet.inf_map);
479  }
480  }
481 
482  case JACOBI_30_00:
483  {
484  switch (fet.inf_map)
485  {
486  case CARTESIAN:
487  return libmesh_make_unique<InfFE<1,JACOBI_30_00,CARTESIAN>>(fet);
488 
489  default:
490  libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << fet.inf_map);
491  }
492  }
493 
494  case LEGENDRE:
495  {
496  switch (fet.inf_map)
497  {
498  case CARTESIAN:
499  return libmesh_make_unique<InfFE<1,LEGENDRE,CARTESIAN>>(fet);
500 
501  default:
502  libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << fet.inf_map);
503  }
504  }
505 
506  case LAGRANGE:
507  {
508  switch (fet.inf_map)
509  {
510  case CARTESIAN:
511  return libmesh_make_unique<InfFE<1,LAGRANGE,CARTESIAN>>(fet);
512 
513  default:
514  libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << fet.inf_map);
515  }
516  }
517 
518  default:
519  libmesh_error_msg("ERROR: Bad FEType.radial_family= " << fet.radial_family);
520  }
521  }
522 
523 
524 
525 
526  // 2D
527  case 2:
528  {
529  switch (fet.radial_family)
530  {
531  case INFINITE_MAP:
532  libmesh_error_msg("ERROR: Can't build an infinite element with FEFamily = " << fet.radial_family);
533 
534  case JACOBI_20_00:
535  {
536  switch (fet.inf_map)
537  {
538  case CARTESIAN:
539  return libmesh_make_unique<InfFE<2,JACOBI_20_00,CARTESIAN>>(fet);
540 
541  default:
542  libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << fet.inf_map);
543  }
544  }
545 
546  case JACOBI_30_00:
547  {
548  switch (fet.inf_map)
549  {
550  case CARTESIAN:
551  return libmesh_make_unique<InfFE<2,JACOBI_30_00,CARTESIAN>>(fet);
552 
553  default:
554  libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << fet.inf_map);
555  }
556  }
557 
558  case LEGENDRE:
559  {
560  switch (fet.inf_map)
561  {
562  case CARTESIAN:
563  return libmesh_make_unique<InfFE<2,LEGENDRE,CARTESIAN>>(fet);
564 
565  default:
566  libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << fet.inf_map);
567  }
568  }
569 
570  case LAGRANGE:
571  {
572  switch (fet.inf_map)
573  {
574  case CARTESIAN:
575  return libmesh_make_unique<InfFE<2,LAGRANGE,CARTESIAN>>(fet);
576 
577  default:
578  libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << fet.inf_map);
579  }
580  }
581 
582  default:
583  libmesh_error_msg("ERROR: Bad FEType.radial_family= " << fet.radial_family);
584  }
585  }
586 
587 
588 
589 
590  // 3D
591  case 3:
592  {
593  switch (fet.radial_family)
594  {
595  case INFINITE_MAP:
596  libmesh_error_msg("ERROR: Don't build an infinite element with FEFamily = " << fet.radial_family);
597 
598  case JACOBI_20_00:
599  {
600  switch (fet.inf_map)
601  {
602  case CARTESIAN:
603  return libmesh_make_unique<InfFE<3,JACOBI_20_00,CARTESIAN>>(fet);
604 
605  default:
606  libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << fet.inf_map);
607  }
608  }
609 
610  case JACOBI_30_00:
611  {
612  switch (fet.inf_map)
613  {
614  case CARTESIAN:
615  return libmesh_make_unique<InfFE<3,JACOBI_30_00,CARTESIAN>>(fet);
616 
617  default:
618  libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << fet.inf_map);
619  }
620  }
621 
622  case LEGENDRE:
623  {
624  switch (fet.inf_map)
625  {
626  case CARTESIAN:
627  return libmesh_make_unique<InfFE<3,LEGENDRE,CARTESIAN>>(fet);
628 
629  default:
630  libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << fet.inf_map);
631  }
632  }
633 
634  case LAGRANGE:
635  {
636  switch (fet.inf_map)
637  {
638  case CARTESIAN:
639  return libmesh_make_unique<InfFE<3,LAGRANGE,CARTESIAN>>(fet);
640 
641  default:
642  libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << fet.inf_map);
643  }
644  }
645 
646  default:
647  libmesh_error_msg("ERROR: Bad FEType.radial_family= " << fet.radial_family);
648  }
649  }
650 
651  default:
652  libmesh_error_msg("Invalid dimension dim = " << dim);
653  }
654 }
const unsigned int dim
Definition: fe_abstract.h:523
InfMapType inf_map
Definition: fe_type.h:258
FEFamily radial_family
Definition: fe_type.h:250
template<>
std::unique_ptr< FEGenericBase< RealGradient > > libMesh::FEGenericBase< RealGradient >::build_InfFE ( const unsigned int  ,
const FEType  
)
inherited

Definition at line 660 of file fe_base.C.

662 {
663  // No vector types defined... YET.
664  libmesh_not_implemented();
665  return std::unique_ptr<FEVectorBase>();
666 }
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::coarsened_dof_values ( const NumericVector< Number > &  global_vector,
const DofMap dof_map,
const Elem coarse_elem,
DenseVector< Number > &  coarse_dofs,
const unsigned int  var,
const bool  use_old_dof_indices = false 
)
staticinherited

Creates a local projection on coarse_elem, based on the DoF values in global_vector for it's children. Computes a vector of coefficients corresponding to dof_indices for only the single given var

Definition at line 794 of file fe_base.C.

References std::abs(), libMesh::C_ONE, libMesh::Elem::child_ptr(), libMesh::Elem::child_ref_range(), libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::FEType::default_quadrature_rule(), libMesh::Elem::dim(), libMesh::DISCONTINUOUS, libMesh::DofMap::dof_indices(), libMesh::FEInterface::dofs_on_edge(), libMesh::FEInterface::dofs_on_side(), libMesh::Elem::edge_index_range(), libMesh::TensorTools::inner_product(), libMesh::FEInterface::inverse_map(), libMesh::Elem::is_child_on_edge(), libMesh::Elem::is_child_on_side(), libMesh::Elem::is_vertex(), libmesh_nullptr, libMesh::Elem::max_descendant_p_level(), libMesh::Elem::n_children(), libMesh::FEInterface::n_dofs(), libMesh::FEInterface::n_dofs_at_node(), n_nodes, libMesh::Elem::n_nodes(), libMesh::DofMap::old_dof_indices(), libMesh::FEType::order, libMesh::Elem::p_level(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::Elem::side_index_range(), libMesh::TOLERANCE, libMesh::Elem::type(), libMesh::DofMap::variable_type(), libMesh::DenseMatrix< T >::zero(), libMesh::DenseVector< T >::zero(), and libMesh::zero.

Referenced by libMesh::JumpErrorEstimator::estimate_error(), and libMesh::ExactErrorEstimator::estimate_error().

800 {
801  // Side/edge local DOF indices
802  std::vector<unsigned int> new_side_dofs, old_side_dofs;
803 
804  // FIXME: what about 2D shells in 3D space?
805  unsigned int dim = elem->dim();
806 
807  // Cache n_children(); it's a virtual call but it's const.
808  const unsigned int n_children = elem->n_children();
809 
810  // We use local FE objects for now
811  // FIXME: we should use more, external objects instead for efficiency
812  const FEType & base_fe_type = dof_map.variable_type(var);
813  std::unique_ptr<FEGenericBase<OutputShape>> fe
814  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
815  std::unique_ptr<FEGenericBase<OutputShape>> fe_coarse
816  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
817 
818  std::unique_ptr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
819  std::unique_ptr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
820  std::unique_ptr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
821  std::vector<Point> coarse_qpoints;
822 
823  // The values of the shape functions at the quadrature
824  // points
825  const std::vector<std::vector<OutputShape>> & phi_values =
826  fe->get_phi();
827  const std::vector<std::vector<OutputShape>> & phi_coarse =
828  fe_coarse->get_phi();
829 
830  // The gradients of the shape functions at the quadrature
831  // points on the child element.
832  const std::vector<std::vector<OutputGradient>> * dphi_values =
834  const std::vector<std::vector<OutputGradient>> * dphi_coarse =
836 
837  const FEContinuity cont = fe->get_continuity();
838 
839  if (cont == C_ONE)
840  {
841  const std::vector<std::vector<OutputGradient>> &
842  ref_dphi_values = fe->get_dphi();
843  dphi_values = &ref_dphi_values;
844  const std::vector<std::vector<OutputGradient>> &
845  ref_dphi_coarse = fe_coarse->get_dphi();
846  dphi_coarse = &ref_dphi_coarse;
847  }
848 
849  // The Jacobian * quadrature weight at the quadrature points
850  const std::vector<Real> & JxW =
851  fe->get_JxW();
852 
853  // The XYZ locations of the quadrature points on the
854  // child element
855  const std::vector<Point> & xyz_values =
856  fe->get_xyz();
857 
858 
859 
860  FEType fe_type = base_fe_type, temp_fe_type;
861  const ElemType elem_type = elem->type();
862  fe_type.order = static_cast<Order>(fe_type.order +
863  elem->max_descendant_p_level());
864 
865  // Number of nodes on parent element
866  const unsigned int n_nodes = elem->n_nodes();
867 
868  // Number of dofs on parent element
869  const unsigned int new_n_dofs =
870  FEInterface::n_dofs(dim, fe_type, elem_type);
871 
872  // Fixed vs. free DoFs on edge/face projections
873  std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
874  std::vector<int> free_dof(new_n_dofs, 0);
875 
878  Ue.resize(new_n_dofs); Ue.zero();
879 
880 
881  // When coarsening, in general, we need a series of
882  // projections to ensure a unique and continuous
883  // solution. We start by interpolating nodes, then
884  // hold those fixed and project edges, then
885  // hold those fixed and project faces, then
886  // hold those fixed and project interiors
887 
888  // Copy node values first
889  {
890  std::vector<dof_id_type> node_dof_indices;
891  if (use_old_dof_indices)
892  dof_map.old_dof_indices (elem, node_dof_indices, var);
893  else
894  dof_map.dof_indices (elem, node_dof_indices, var);
895 
896  unsigned int current_dof = 0;
897  for (unsigned int n=0; n!= n_nodes; ++n)
898  {
899  // FIXME: this should go through the DofMap,
900  // not duplicate dof_indices code badly!
901  const unsigned int my_nc =
902  FEInterface::n_dofs_at_node (dim, fe_type,
903  elem_type, n);
904  if (!elem->is_vertex(n))
905  {
906  current_dof += my_nc;
907  continue;
908  }
909 
910  temp_fe_type = base_fe_type;
911  // We're assuming here that child n shares vertex n,
912  // which is wrong on non-simplices right now
913  // ... but this code isn't necessary except on elements
914  // where p refinement creates more vertex dofs; we have
915  // no such elements yet.
916  /*
917  if (elem->child_ptr(n)->p_level() < elem->p_level())
918  {
919  temp_fe_type.order =
920  static_cast<Order>(temp_fe_type.order +
921  elem->child_ptr(n)->p_level());
922  }
923  */
924  const unsigned int nc =
925  FEInterface::n_dofs_at_node (dim, temp_fe_type,
926  elem_type, n);
927  for (unsigned int i=0; i!= nc; ++i)
928  {
929  Ue(current_dof) =
930  old_vector(node_dof_indices[current_dof]);
931  dof_is_fixed[current_dof] = true;
932  current_dof++;
933  }
934  }
935  }
936 
937  // In 3D, project any edge values next
938  if (dim > 2 && cont != DISCONTINUOUS)
939  for (auto e : elem->edge_index_range())
940  {
941  FEInterface::dofs_on_edge(elem, dim, fe_type,
942  e, new_side_dofs);
943 
944  // Some edge dofs are on nodes and already
945  // fixed, others are free to calculate
946  unsigned int free_dofs = 0;
947  for (std::size_t i=0; i != new_side_dofs.size(); ++i)
948  if (!dof_is_fixed[new_side_dofs[i]])
949  free_dof[free_dofs++] = i;
950  Ke.resize (free_dofs, free_dofs); Ke.zero();
951  Fe.resize (free_dofs); Fe.zero();
952  // The new edge coefficients
953  DenseVector<Number> Uedge(free_dofs);
954 
955  // Add projection terms from each child sharing
956  // this edge
957  for (unsigned int c=0; c != n_children; ++c)
958  {
959  if (!elem->is_child_on_edge(c,e))
960  continue;
961  const Elem * child = elem->child_ptr(c);
962 
963  std::vector<dof_id_type> child_dof_indices;
964  if (use_old_dof_indices)
965  dof_map.old_dof_indices (child,
966  child_dof_indices, var);
967  else
968  dof_map.dof_indices (child,
969  child_dof_indices, var);
970  const unsigned int child_n_dofs =
971  cast_int<unsigned int>
972  (child_dof_indices.size());
973 
974  temp_fe_type = base_fe_type;
975  temp_fe_type.order =
976  static_cast<Order>(temp_fe_type.order +
977  child->p_level());
978 
979  FEInterface::dofs_on_edge(child, dim,
980  temp_fe_type, e, old_side_dofs);
981 
982  // Initialize both child and parent FE data
983  // on the child's edge
984  fe->attach_quadrature_rule (qedgerule.get());
985  fe->edge_reinit (child, e);
986  const unsigned int n_qp = qedgerule->n_points();
987 
988  FEInterface::inverse_map (dim, fe_type, elem,
989  xyz_values, coarse_qpoints);
990 
991  fe_coarse->reinit(elem, &coarse_qpoints);
992 
993  // Loop over the quadrature points
994  for (unsigned int qp=0; qp<n_qp; qp++)
995  {
996  // solution value at the quadrature point
997  OutputNumber fineval = libMesh::zero;
998  // solution grad at the quadrature point
999  OutputNumberGradient finegrad;
1000 
1001  // Sum the solution values * the DOF
1002  // values at the quadrature point to
1003  // get the solution value and gradient.
1004  for (unsigned int i=0; i<child_n_dofs;
1005  i++)
1006  {
1007  fineval +=
1008  (old_vector(child_dof_indices[i])*
1009  phi_values[i][qp]);
1010  if (cont == C_ONE)
1011  finegrad += (*dphi_values)[i][qp] *
1012  old_vector(child_dof_indices[i]);
1013  }
1014 
1015  // Form edge projection matrix
1016  for (std::size_t sidei=0, freei=0; sidei != new_side_dofs.size(); ++sidei)
1017  {
1018  unsigned int i = new_side_dofs[sidei];
1019  // fixed DoFs aren't test functions
1020  if (dof_is_fixed[i])
1021  continue;
1022  for (std::size_t sidej=0, freej=0; sidej != new_side_dofs.size(); ++sidej)
1023  {
1024  unsigned int j =
1025  new_side_dofs[sidej];
1026  if (dof_is_fixed[j])
1027  Fe(freei) -=
1028  TensorTools::inner_product(phi_coarse[i][qp],
1029  phi_coarse[j][qp]) *
1030  JxW[qp] * Ue(j);
1031  else
1032  Ke(freei,freej) +=
1033  TensorTools::inner_product(phi_coarse[i][qp],
1034  phi_coarse[j][qp]) *
1035  JxW[qp];
1036  if (cont == C_ONE)
1037  {
1038  if (dof_is_fixed[j])
1039  Fe(freei) -=
1040  TensorTools::inner_product((*dphi_coarse)[i][qp],
1041  (*dphi_coarse)[j][qp]) *
1042  JxW[qp] * Ue(j);
1043  else
1044  Ke(freei,freej) +=
1045  TensorTools::inner_product((*dphi_coarse)[i][qp],
1046  (*dphi_coarse)[j][qp]) *
1047  JxW[qp];
1048  }
1049  if (!dof_is_fixed[j])
1050  freej++;
1051  }
1052  Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp],
1053  fineval) * JxW[qp];
1054  if (cont == C_ONE)
1055  Fe(freei) +=
1056  TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1057  freei++;
1058  }
1059  }
1060  }
1061  Ke.cholesky_solve(Fe, Uedge);
1062 
1063  // Transfer new edge solutions to element
1064  for (unsigned int i=0; i != free_dofs; ++i)
1065  {
1066  Number & ui = Ue(new_side_dofs[free_dof[i]]);
1067  libmesh_assert(std::abs(ui) < TOLERANCE ||
1068  std::abs(ui - Uedge(i)) < TOLERANCE);
1069  ui = Uedge(i);
1070  dof_is_fixed[new_side_dofs[free_dof[i]]] = true;
1071  }
1072  }
1073 
1074  // Project any side values (edges in 2D, faces in 3D)
1075  if (dim > 1 && cont != DISCONTINUOUS)
1076  for (auto s : elem->side_index_range())
1077  {
1078  FEInterface::dofs_on_side(elem, dim, fe_type,
1079  s, new_side_dofs);
1080 
1081  // Some side dofs are on nodes/edges and already
1082  // fixed, others are free to calculate
1083  unsigned int free_dofs = 0;
1084  for (std::size_t i=0; i != new_side_dofs.size(); ++i)
1085  if (!dof_is_fixed[new_side_dofs[i]])
1086  free_dof[free_dofs++] = i;
1087  Ke.resize (free_dofs, free_dofs); Ke.zero();
1088  Fe.resize (free_dofs); Fe.zero();
1089  // The new side coefficients
1090  DenseVector<Number> Uside(free_dofs);
1091 
1092  // Add projection terms from each child sharing
1093  // this side
1094  for (unsigned int c=0; c != n_children; ++c)
1095  {
1096  if (!elem->is_child_on_side(c,s))
1097  continue;
1098  const Elem * child = elem->child_ptr(c);
1099 
1100  std::vector<dof_id_type> child_dof_indices;
1101  if (use_old_dof_indices)
1102  dof_map.old_dof_indices (child,
1103  child_dof_indices, var);
1104  else
1105  dof_map.dof_indices (child,
1106  child_dof_indices, var);
1107  const unsigned int child_n_dofs =
1108  cast_int<unsigned int>
1109  (child_dof_indices.size());
1110 
1111  temp_fe_type = base_fe_type;
1112  temp_fe_type.order =
1113  static_cast<Order>(temp_fe_type.order +
1114  child->p_level());
1115 
1116  FEInterface::dofs_on_side(child, dim,
1117  temp_fe_type, s, old_side_dofs);
1118 
1119  // Initialize both child and parent FE data
1120  // on the child's side
1121  fe->attach_quadrature_rule (qsiderule.get());
1122  fe->reinit (child, s);
1123  const unsigned int n_qp = qsiderule->n_points();
1124 
1125  FEInterface::inverse_map (dim, fe_type, elem,
1126  xyz_values, coarse_qpoints);
1127 
1128  fe_coarse->reinit(elem, &coarse_qpoints);
1129 
1130  // Loop over the quadrature points
1131  for (unsigned int qp=0; qp<n_qp; qp++)
1132  {
1133  // solution value at the quadrature point
1134  OutputNumber fineval = libMesh::zero;
1135  // solution grad at the quadrature point
1136  OutputNumberGradient finegrad;
1137 
1138  // Sum the solution values * the DOF
1139  // values at the quadrature point to
1140  // get the solution value and gradient.
1141  for (unsigned int i=0; i<child_n_dofs;
1142  i++)
1143  {
1144  fineval +=
1145  old_vector(child_dof_indices[i]) *
1146  phi_values[i][qp];
1147  if (cont == C_ONE)
1148  finegrad += (*dphi_values)[i][qp] *
1149  old_vector(child_dof_indices[i]);
1150  }
1151 
1152  // Form side projection matrix
1153  for (std::size_t sidei=0, freei=0; sidei != new_side_dofs.size(); ++sidei)
1154  {
1155  unsigned int i = new_side_dofs[sidei];
1156  // fixed DoFs aren't test functions
1157  if (dof_is_fixed[i])
1158  continue;
1159  for (std::size_t sidej=0, freej=0; sidej != new_side_dofs.size(); ++sidej)
1160  {
1161  unsigned int j =
1162  new_side_dofs[sidej];
1163  if (dof_is_fixed[j])
1164  Fe(freei) -=
1165  TensorTools::inner_product(phi_coarse[i][qp],
1166  phi_coarse[j][qp]) *
1167  JxW[qp] * Ue(j);
1168  else
1169  Ke(freei,freej) +=
1170  TensorTools::inner_product(phi_coarse[i][qp],
1171  phi_coarse[j][qp]) *
1172  JxW[qp];
1173  if (cont == C_ONE)
1174  {
1175  if (dof_is_fixed[j])
1176  Fe(freei) -=
1177  TensorTools::inner_product((*dphi_coarse)[i][qp],
1178  (*dphi_coarse)[j][qp]) *
1179  JxW[qp] * Ue(j);
1180  else
1181  Ke(freei,freej) +=
1182  TensorTools::inner_product((*dphi_coarse)[i][qp],
1183  (*dphi_coarse)[j][qp]) *
1184  JxW[qp];
1185  }
1186  if (!dof_is_fixed[j])
1187  freej++;
1188  }
1189  Fe(freei) += TensorTools::inner_product(fineval, phi_coarse[i][qp]) * JxW[qp];
1190  if (cont == C_ONE)
1191  Fe(freei) +=
1192  TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1193  freei++;
1194  }
1195  }
1196  }
1197  Ke.cholesky_solve(Fe, Uside);
1198 
1199  // Transfer new side solutions to element
1200  for (unsigned int i=0; i != free_dofs; ++i)
1201  {
1202  Number & ui = Ue(new_side_dofs[free_dof[i]]);
1203  libmesh_assert(std::abs(ui) < TOLERANCE ||
1204  std::abs(ui - Uside(i)) < TOLERANCE);
1205  ui = Uside(i);
1206  dof_is_fixed[new_side_dofs[free_dof[i]]] = true;
1207  }
1208  }
1209 
1210  // Project the interior values, finally
1211 
1212  // Some interior dofs are on nodes/edges/sides and
1213  // already fixed, others are free to calculate
1214  unsigned int free_dofs = 0;
1215  for (unsigned int i=0; i != new_n_dofs; ++i)
1216  if (!dof_is_fixed[i])
1217  free_dof[free_dofs++] = i;
1218  Ke.resize (free_dofs, free_dofs); Ke.zero();
1219  Fe.resize (free_dofs); Fe.zero();
1220  // The new interior coefficients
1221  DenseVector<Number> Uint(free_dofs);
1222 
1223  // Add projection terms from each child
1224  for (auto & child : elem->child_ref_range())
1225  {
1226  std::vector<dof_id_type> child_dof_indices;
1227  if (use_old_dof_indices)
1228  dof_map.old_dof_indices (&child,
1229  child_dof_indices, var);
1230  else
1231  dof_map.dof_indices (&child,
1232  child_dof_indices, var);
1233  const unsigned int child_n_dofs =
1234  cast_int<unsigned int>
1235  (child_dof_indices.size());
1236 
1237  // Initialize both child and parent FE data
1238  // on the child's quadrature points
1239  fe->attach_quadrature_rule (qrule.get());
1240  fe->reinit (&child);
1241  const unsigned int n_qp = qrule->n_points();
1242 
1243  FEInterface::inverse_map (dim, fe_type, elem,
1244  xyz_values, coarse_qpoints);
1245 
1246  fe_coarse->reinit(elem, &coarse_qpoints);
1247 
1248  // Loop over the quadrature points
1249  for (unsigned int qp=0; qp<n_qp; qp++)
1250  {
1251  // solution value at the quadrature point
1252  OutputNumber fineval = libMesh::zero;
1253  // solution grad at the quadrature point
1254  OutputNumberGradient finegrad;
1255 
1256  // Sum the solution values * the DOF
1257  // values at the quadrature point to
1258  // get the solution value and gradient.
1259  for (unsigned int i=0; i<child_n_dofs; i++)
1260  {
1261  fineval +=
1262  (old_vector(child_dof_indices[i]) *
1263  phi_values[i][qp]);
1264  if (cont == C_ONE)
1265  finegrad += (*dphi_values)[i][qp] *
1266  old_vector(child_dof_indices[i]);
1267  }
1268 
1269  // Form interior projection matrix
1270  for (unsigned int i=0, freei=0;
1271  i != new_n_dofs; ++i)
1272  {
1273  // fixed DoFs aren't test functions
1274  if (dof_is_fixed[i])
1275  continue;
1276  for (unsigned int j=0, freej=0; j !=
1277  new_n_dofs; ++j)
1278  {
1279  if (dof_is_fixed[j])
1280  Fe(freei) -=
1281  TensorTools::inner_product(phi_coarse[i][qp],
1282  phi_coarse[j][qp]) *
1283  JxW[qp] * Ue(j);
1284  else
1285  Ke(freei,freej) +=
1286  TensorTools::inner_product(phi_coarse[i][qp],
1287  phi_coarse[j][qp]) *
1288  JxW[qp];
1289  if (cont == C_ONE)
1290  {
1291  if (dof_is_fixed[j])
1292  Fe(freei) -=
1293  TensorTools::inner_product((*dphi_coarse)[i][qp],
1294  (*dphi_coarse)[j][qp]) *
1295  JxW[qp] * Ue(j);
1296  else
1297  Ke(freei,freej) +=
1298  TensorTools::inner_product((*dphi_coarse)[i][qp],
1299  (*dphi_coarse)[j][qp]) *
1300  JxW[qp];
1301  }
1302  if (!dof_is_fixed[j])
1303  freej++;
1304  }
1305  Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp], fineval) *
1306  JxW[qp];
1307  if (cont == C_ONE)
1308  Fe(freei) += TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1309  freei++;
1310  }
1311  }
1312  }
1313  Ke.cholesky_solve(Fe, Uint);
1314 
1315  // Transfer new interior solutions to element
1316  for (unsigned int i=0; i != free_dofs; ++i)
1317  {
1318  Number & ui = Ue(free_dof[i]);
1319  libmesh_assert(std::abs(ui) < TOLERANCE ||
1320  std::abs(ui - Uint(i)) < TOLERANCE);
1321  ui = Uint(i);
1322  // We should be fixing all dofs by now; no need to keep track of
1323  // that unless we're debugging
1324 #ifndef NDEBUG
1325  dof_is_fixed[free_dof[i]] = true;
1326 #endif
1327  }
1328 
1329 #ifndef NDEBUG
1330  // Make sure every DoF got reached!
1331  for (unsigned int i=0; i != new_n_dofs; ++i)
1332  libmesh_assert(dof_is_fixed[i]);
1333 #endif
1334 }
static void dofs_on_edge(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int e, std::vector< unsigned int > &di)
Definition: fe_interface.C:493
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
double abs(double a)
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:411
virtual void zero() libmesh_override
Definition: dense_matrix.h:792
unsigned int p_level() const
Definition: elem.h:2423
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1719
static void dofs_on_side(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int s, std::vector< unsigned int > &di)
Definition: fe_interface.C:480
void resize(const unsigned int n)
Definition: dense_vector.h:350
virtual void zero() libmesh_override
Definition: dense_vector.h:374
The base class for all geometric element types.
Definition: elem.h:90
const class libmesh_nullptr_t libmesh_nullptr
TensorTools::IncrementRank< OutputNumber >::type OutputNumberGradient
Definition: fe_base.h:125
OrderWrapper order
Definition: fe_type.h:198
static const Real TOLERANCE
void old_dof_indices(const Elem *const elem, std::vector< dof_id_type > &di, const unsigned int vn=libMesh::invalid_uint) const
Definition: dof_map.C:2385
const Number zero
Definition: libmesh.h:178
SimpleRange< ChildRefIter > child_ref_range()
Definition: elem.h:1700
const dof_id_type n_nodes
Definition: tecplot_io.C:67
const unsigned int dim
Definition: fe_abstract.h:523
const Elem * child_ptr(unsigned int i) const
Definition: elem.h:2446
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:547
TensorTools::MakeNumber< OutputShape >::type OutputNumber
Definition: fe_base.h:124
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition: fe_interface.C:430
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:31
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
A matrix object used for finite element assembly and numerics.
Definition: dense_matrix.h:54
unsigned int n_points() const
Definition: quadrature.h:114
boostcopy::enable_if_c< ScalarTraits< T >::value &&ScalarTraits< T2 >::value, typename CompareTypes< T, T2 >::supertype >::type inner_product(const T &a, const T2 &b)
Definition: tensor_tools.h:47
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1924
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::coarsened_dof_values ( const NumericVector< Number > &  global_vector,
const DofMap dof_map,
const Elem coarse_elem,
DenseVector< Number > &  coarse_dofs,
const bool  use_old_dof_indices = false 
)
staticinherited

Creates a local projection on coarse_elem, based on the DoF values in global_vector for it's children. Computes a vector of coefficients corresponding to all dof_indices.

Definition at line 1340 of file fe_base.C.

References libMesh::DenseVector< T >::append(), libMesh::DofMap::n_variables(), and libMesh::DenseVector< T >::resize().

1345 {
1346  Ue.resize(0);
1347 
1348  for (unsigned int v=0; v != dof_map.n_variables(); ++v)
1349  {
1350  DenseVector<Number> Usub;
1351 
1352  coarsened_dof_values(old_vector, dof_map, elem, Usub,
1353  use_old_dof_indices);
1354 
1355  Ue.append (Usub);
1356  }
1357 }
static void coarsened_dof_values(const NumericVector< Number > &global_vector, const DofMap &dof_map, const Elem *coarse_elem, DenseVector< Number > &coarse_dofs, const unsigned int var, const bool use_old_dof_indices=false)
Definition: fe_base.C:794
unsigned int n_variables() const
Definition: dof_map.h:486
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void libMesh::InfFE< Dim, T_radial, T_map >::combine_base_radial ( const Elem inf_elem)
protected

Combines the shape functions, which were formed in init_shape_functions(Elem *), with geometric data. Has to be called every time the geometric configuration changes. Afterward, the fields are ready to be used to compute global derivatives, the jacobian etc, see FEAbstract::compute_map().

Definition at line 715 of file inf_fe.C.

References libMesh::InfFE< Dim, T_radial, T_map >::_base_node_index, libMesh::InfFE< Dim, T_radial, T_map >::_base_shape_index, libMesh::FEAbstract::_fe_map, libMesh::InfFE< Dim, T_radial, T_map >::_radial_node_index, libMesh::InfFE< Dim, T_radial, T_map >::_radial_shape_index, libMesh::InfFE< Dim, T_radial, T_map >::base_elem, libMesh::InfFE< Dim, T_radial, T_map >::base_fe, libMesh::InfFE< Dim, T_radial, T_map >::base_qrule, libMesh::InfFE< Dim, T_radial, T_map >::dist, libMesh::InfFE< Dim, T_radial, T_map >::dmodedv, libMesh::InfFE< Dim, T_radial, T_map >::dphasedeta, libMesh::InfFE< Dim, T_radial, T_map >::dphasedxi, libMesh::InfFE< Dim, T_radial, T_map >::dphasedzeta, libMesh::FEGenericBase< OutputType >::dphideta, libMesh::FEGenericBase< OutputType >::dphidxi, libMesh::FEGenericBase< OutputType >::dphidzeta, libMesh::InfFE< Dim, T_radial, T_map >::dradialdv_map, libMesh::InfFE< Dim, T_radial, T_map >::dsomdv, libMesh::FEAbstract::fe_type, libMesh::InfFE< Dim, T_radial, T_map >::Base::get_elem_type(), libMesh::InfFE< Dim, T_radial, T_map >::mode, libMesh::InfFE< Dim, T_radial, T_map >::Radial::n_dofs(), libMesh::Elem::origin(), libMesh::FEGenericBase< OutputType >::phi, libMesh::InfFE< Dim, T_radial, T_map >::radial_map, libMesh::FEType::radial_order, libMesh::InfFE< Dim, T_radial, T_map >::radial_qrule, libMesh::InfFE< Dim, T_radial, T_map >::som, and libMesh::Elem::type().

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::init_base_shape_functions(), and libMesh::InfFE< Dim, T_radial, T_map >::reinit().

716 {
717  libmesh_assert(inf_elem);
718  // at least check whether the base element type is correct.
719  // otherwise this version of computing dist would give problems
720  libmesh_assert_equal_to (base_elem->type(), Base::get_elem_type(inf_elem->type()));
721 
722  // Start logging the combination of radial and base parts
723  LOG_SCOPE("combine_base_radial()", "InfFE");
724 
725  // zero the phase, since it is to be summed up
726  std::fill (dphasedxi.begin(), dphasedxi.end(), 0.);
727  std::fill (dphasedeta.begin(), dphasedeta.end(), 0.);
728  std::fill (dphasedzeta.begin(), dphasedzeta.end(), 0.);
729 
730 
731  const unsigned int n_base_mapping_sf = cast_int<unsigned int>(dist.size());
732  const Point origin = inf_elem->origin();
733 
734  // for each new infinite element, compute the radial distances
735  for (unsigned int n=0; n<n_base_mapping_sf; n++)
736  dist[n] = Point(base_elem->point(n) - origin).norm();
737 
738 
739  switch (Dim)
740  {
741  // 1D
742  case 1:
743  {
744  libmesh_not_implemented();
745  break;
746  }
747 
748  // 2D
749  case 2:
750  {
751  libmesh_not_implemented();
752  break;
753  }
754 
755  // 3D
756  case 3:
757  {
758  // fast access to the approximation and mapping shapes of base_fe
759  const std::vector<std::vector<Real>> & S = base_fe->phi;
760  const std::vector<std::vector<Real>> & Ss = base_fe->dphidxi;
761  const std::vector<std::vector<Real>> & St = base_fe->dphideta;
762  const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
763  const std::vector<std::vector<Real>> & Ss_map = (base_fe->get_fe_map()).get_dphidxi_map();
764  const std::vector<std::vector<Real>> & St_map = (base_fe->get_fe_map()).get_dphideta_map();
765 
766  const unsigned int n_radial_qp = som.size();
767  if (radial_qrule)
768  libmesh_assert_equal_to(n_radial_qp, radial_qrule->n_points());
769  const unsigned int n_base_qp = S_map[0].size();
770  if (base_qrule)
771  libmesh_assert_equal_to(n_base_qp, base_qrule->n_points());
772 
773  const unsigned int n_total_mapping_sf =
774  cast_int<unsigned int>(radial_map.size()) * n_base_mapping_sf;
775 
776  const unsigned int n_total_approx_sf = Radial::n_dofs(fe_type.radial_order) * base_fe->n_shape_functions();
777 
778 
779  // compute the phase term derivatives
780  {
781  unsigned int tp=0;
782  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
783  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
784  {
785  // sum over all base shapes, to get the average distance
786  for (unsigned int i=0; i<n_base_mapping_sf; i++)
787  {
788  dphasedxi[tp] += Ss_map[i][bp] * dist[i] * radial_map [1][rp];
789  dphasedeta[tp] += St_map[i][bp] * dist[i] * radial_map [1][rp];
790  dphasedzeta[tp] += S_map [i][bp] * dist[i] * dradialdv_map[1][rp];
791  }
792 
793  tp++;
794 
795  } // loop radial and base qps
796  }
797 
798  libmesh_assert_equal_to (phi.size(), n_total_approx_sf);
799  libmesh_assert_equal_to (dphidxi.size(), n_total_approx_sf);
800  libmesh_assert_equal_to (dphideta.size(), n_total_approx_sf);
801  libmesh_assert_equal_to (dphidzeta.size(), n_total_approx_sf);
802 
803  // compute the overall approximation shape functions,
804  // pick the appropriate radial and base shapes through using
805  // _base_shape_index and _radial_shape_index
806  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
807  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
808  for (unsigned int ti=0; ti<n_total_approx_sf; ti++) // over _all_ approx_sf
809  {
810  // let the index vectors take care of selecting the appropriate base/radial shape
811  const unsigned int bi = _base_shape_index [ti];
812  const unsigned int ri = _radial_shape_index[ti];
813  phi [ti][bp+rp*n_base_qp] = S [bi][bp] * mode[ri][rp] * som[rp];
814  dphidxi [ti][bp+rp*n_base_qp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
815  dphideta [ti][bp+rp*n_base_qp] = St[bi][bp] * mode[ri][rp] * som[rp];
816  dphidzeta[ti][bp+rp*n_base_qp] = S [bi][bp]
817  * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
818  }
819 
820  std::vector<std::vector<Real>> & phi_map = this->_fe_map->get_phi_map();
821  std::vector<std::vector<Real>> & dphidxi_map = this->_fe_map->get_dphidxi_map();
822  std::vector<std::vector<Real>> & dphideta_map = this->_fe_map->get_dphideta_map();
823  std::vector<std::vector<Real>> & dphidzeta_map = this->_fe_map->get_dphidzeta_map();
824 
825  libmesh_assert_equal_to (phi_map.size(), n_total_mapping_sf);
826  libmesh_assert_equal_to (dphidxi_map.size(), n_total_mapping_sf);
827  libmesh_assert_equal_to (dphideta_map.size(), n_total_mapping_sf);
828  libmesh_assert_equal_to (dphidzeta_map.size(), n_total_mapping_sf);
829 
830  // compute the overall mapping functions,
831  // pick the appropriate radial and base entries through using
832  // _base_node_index and _radial_node_index
833  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
834  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
835  for (unsigned int ti=0; ti<n_total_mapping_sf; ti++) // over all mapping shapes
836  {
837  // let the index vectors take care of selecting the appropriate base/radial mapping shape
838  const unsigned int bi = _base_node_index [ti];
839  const unsigned int ri = _radial_node_index[ti];
840  phi_map [ti][bp+rp*n_base_qp] = S_map [bi][bp] * radial_map [ri][rp];
841  dphidxi_map [ti][bp+rp*n_base_qp] = Ss_map[bi][bp] * radial_map [ri][rp];
842  dphideta_map [ti][bp+rp*n_base_qp] = St_map[bi][bp] * radial_map [ri][rp];
843  dphidzeta_map[ti][bp+rp*n_base_qp] = S_map [bi][bp] * dradialdv_map[ri][rp];
844  }
845 
846  break;
847  }
848 
849  default:
850  libmesh_error_msg("Unsupported Dim = " << Dim);
851  }
852 }
std::vector< std::vector< OutputShape > > dphidxi
Definition: fe_base.h:519
std::vector< std::vector< OutputShape > > dphidzeta
Definition: fe_base.h:529
std::vector< Real > dsomdv
Definition: inf_fe.h:633
std::vector< std::vector< Real > > dradialdv_map
Definition: inf_fe.h:656
std::vector< std::vector< Real > > radial_map
Definition: inf_fe.h:650
std::vector< Real > dphasedxi
Definition: inf_fe.h:663
std::unique_ptr< QBase > radial_qrule
Definition: inf_fe.h:757
std::vector< unsigned int > _base_node_index
Definition: inf_fe.h:702
OrderWrapper radial_order
Definition: fe_type.h:237
static ElemType get_elem_type(const ElemType type)
std::vector< std::vector< Real > > mode
Definition: inf_fe.h:639
std::vector< std::vector< OutputShape > > phi
Definition: fe_base.h:499
std::unique_ptr< Elem > base_elem
Definition: inf_fe.h:763
std::unique_ptr< FEBase > base_fe
Definition: inf_fe.h:771
std::vector< Real > dphasedzeta
Definition: inf_fe.h:677
std::vector< Real > dphasedeta
Definition: inf_fe.h:670
static unsigned int n_dofs(const Order o_radial)
Definition: inf_fe.h:147
std::vector< unsigned int > _radial_node_index
Definition: inf_fe.h:692
std::vector< Real > som
Definition: inf_fe.h:628
std::unique_ptr< QBase > base_qrule
Definition: inf_fe.h:751
std::vector< unsigned int > _radial_shape_index
Definition: inf_fe.h:712
std::vector< unsigned int > _base_shape_index
Definition: inf_fe.h:722
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:517
std::vector< Real > dist
Definition: inf_fe.h:611
std::vector< std::vector< Real > > dmodedv
Definition: inf_fe.h:645
std::vector< std::vector< OutputShape > > dphideta
Definition: fe_base.h:524
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void libMesh::InfFE< Dim, T_radial, T_map >::compute_data ( const FEType fe_t,
const Elem inf_elem,
FEComputeData data 
)
static

Generalized version of shape(), takes an Elem *. The data contains both input and output parameters. For frequency domain simulations, the complex-valued shape is returned. In time domain both the computed shape, and the phase is returned.

Note
The phase (proportional to the distance of the Point data.p from the envelope) is actually a measure how far into the future the results are.

Definition at line 241 of file inf_fe_static.C.

References libMesh::Elem::build_side_ptr(), libMesh::InfFE< Dim, T_radial, T_map >::Radial::decay(), libMesh::InfFE< Dim, T_radial, T_map >::eval(), libMesh::FEComputeData::frequency, libMesh::INFEDGE2, libMesh::Elem::origin(), libMesh::FEComputeData::p, libMesh::FEComputeData::phase, libMesh::pi, libMesh::Elem::point(), libMesh::FEType::radial_order, libMesh::Real, libMesh::FEComputeData::shape, libMesh::FE< Dim, T >::shape(), libMesh::FEInterface::shape(), libMesh::FEComputeData::speed, and libMesh::Elem::type().

Referenced by libMesh::FEInterface::ifem_compute_data(), and libMesh::InfFE< Dim, T_radial, T_map >::~InfFE().

244 {
245  libmesh_assert(inf_elem);
246  libmesh_assert_not_equal_to (Dim, 0);
247 
248  const Order o_radial (fet.radial_order);
249  const Order radial_mapping_order (Radial::mapping_order());
250  const Point & p (data.p);
251  const Real v (p(Dim-1));
252  std::unique_ptr<const Elem> base_el (inf_elem->build_side_ptr(0));
253 
254  /*
255  * compute \p interpolated_dist containing the mapping-interpolated
256  * distance of the base point to the origin. This is the same
257  * for all shape functions. Set \p interpolated_dist to 0, it
258  * is added to.
259  */
260  Real interpolated_dist = 0.;
261  switch (Dim)
262  {
263  case 1:
264  {
265  libmesh_assert_equal_to (inf_elem->type(), INFEDGE2);
266  interpolated_dist = Point(inf_elem->point(0) - inf_elem->point(1)).norm();
267  break;
268  }
269 
270  case 2:
271  {
272  const unsigned int n_base_nodes = base_el->n_nodes();
273 
274  const Point origin = inf_elem->origin();
275  const Order base_mapping_order (base_el->default_order());
276  const ElemType base_mapping_elem_type (base_el->type());
277 
278  // interpolate the base nodes' distances
279  for (unsigned int n=0; n<n_base_nodes; n++)
280  interpolated_dist += Point(base_el->point(n) - origin).norm()
281  * FE<1,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
282  break;
283  }
284 
285  case 3:
286  {
287  const unsigned int n_base_nodes = base_el->n_nodes();
288 
289  const Point origin = inf_elem->origin();
290  const Order base_mapping_order (base_el->default_order());
291  const ElemType base_mapping_elem_type (base_el->type());
292 
293  // interpolate the base nodes' distances
294  for (unsigned int n=0; n<n_base_nodes; n++)
295  interpolated_dist += Point(base_el->point(n) - origin).norm()
296  * FE<2,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
297  break;
298  }
299 
300  default:
301  libmesh_error_msg("Unknown Dim = " << Dim);
302  }
303 
304 
305 
306 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
307 
308  // assumption on time-harmonic behavior
309  Number sign_i (0,1.);
310 
311  // the wave number
312  const Number wavenumber = 2. * libMesh::pi * data.frequency / data.speed;
313 
314  // the exponent for time-harmonic behavior
315  const Number exponent = sign_i /* imaginary unit */
316  * wavenumber /* k (can be complex) */
317  * interpolated_dist /* together with next line: */
318  * InfFE<Dim,INFINITE_MAP,T_map>::eval(v, radial_mapping_order, 1); /* phase(s,t,v) */
319 
320  const Number time_harmonic = exp(exponent); /* e^(sign*i*k*phase(s,t,v)) */
321 
322  /*
323  * compute \p shape for all dof in the element
324  */
325  if (Dim > 1)
326  {
327  const unsigned int n_dof = n_dofs (fet, inf_elem->type());
328  data.shape.resize(n_dof);
329 
330  for (unsigned int i=0; i<n_dof; i++)
331  {
332  // compute base and radial shape indices
333  unsigned int i_base, i_radial;
334  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
335 
336  data.shape[i] = (InfFE<Dim,T_radial,T_map>::Radial::decay(v) /* (1.-v)/2. in 3D */
337  * FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p) /* S_n(s,t) */
338  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)) /* L_n(v) */
339  * time_harmonic; /* e^(sign*i*k*phase(s,t,v) */
340  }
341  }
342  else
343  libmesh_error_msg("compute_data() for 1-dimensional InfFE not implemented.");
344 
345 #else
346 
347  const Real speed = data.speed;
348 
349  /*
350  * This is quite weird: the phase is actually
351  * a measure how advanced the pressure is that
352  * we compute. In other words: the further away
353  * the node \p data.p is, the further we look into
354  * the future...
355  */
356  data.phase = interpolated_dist /* phase(s,t,v)/c */
357  * InfFE<Dim,INFINITE_MAP,T_map>::eval(v, radial_mapping_order, 1) / speed;
358 
359  if (Dim > 1)
360  {
361  const unsigned int n_dof = n_dofs (fet, inf_elem->type());
362  data.shape.resize(n_dof);
363 
364  for (unsigned int i=0; i<n_dof; i++)
365  {
366  // compute base and radial shape indices
367  unsigned int i_base, i_radial;
368  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
369 
370  data.shape[i] = InfFE<Dim,T_radial,T_map>::Radial::decay(v) /* (1.-v)/2. in 3D */
371  * FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p) /* S_n(s,t) */
372  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial); /* L_n(v) */
373  }
374  }
375  else
376  libmesh_error_msg("compute_data() for 1-dimensional InfFE not implemented.");
377 
378 #endif
379 }
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
static unsigned int n_dofs(const FEType &fet, const ElemType inf_elem_type)
Definition: inf_fe_static.C:55
static Real decay(const Real v)
Definition: inf_fe.h:826
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition: fe_interface.C:614
static Real eval(Real v, Order o_radial, unsigned int i)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static void compute_shape_indices(const FEType &fet, const ElemType inf_elem_type, const unsigned int i, unsigned int &base_shape, unsigned int &radial_shape)
IterBase * data
static Order mapping_order()
Definition: inf_fe.h:131
const Real pi
Definition: libmesh.h:172
void libMesh::FEAbstract::compute_node_constraints ( NodeConstraints constraints,
const Elem elem 
)
staticinherited

Computes the nodal constraint contributions (for non-conforming adapted meshes), using Lagrange geometry

Definition at line 793 of file fe_abstract.C.

References std::abs(), libMesh::Elem::build_side_ptr(), libMesh::Elem::default_order(), libMesh::Elem::dim(), libMesh::FEAbstract::fe_type, libMesh::FEInterface::inverse_map(), libMesh::LAGRANGE, libMesh::Elem::level(), libmesh_nullptr, libMesh::FEInterface::n_dofs(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::parent(), libMesh::Real, libMesh::remote_elem, libMesh::FEInterface::shape(), libMesh::Elem::side_index_range(), libMesh::Threads::spin_mtx, and libMesh::Elem::subactive().

795 {
796  libmesh_assert(elem);
797 
798  const unsigned int Dim = elem->dim();
799 
800  // Only constrain elements in 2,3D.
801  if (Dim == 1)
802  return;
803 
804  // Only constrain active and ancestor elements
805  if (elem->subactive())
806  return;
807 
808  // We currently always use LAGRANGE mappings for geometry
809  const FEType fe_type(elem->default_order(), LAGRANGE);
810 
811  std::vector<const Node *> my_nodes, parent_nodes;
812 
813  // Look at the element faces. Check to see if we need to
814  // build constraints.
815  for (auto s : elem->side_index_range())
816  if (elem->neighbor_ptr(s) != libmesh_nullptr &&
817  elem->neighbor_ptr(s) != remote_elem)
818  if (elem->neighbor_ptr(s)->level() < elem->level()) // constrain dofs shared between
819  { // this element and ones coarser
820  // than this element.
821  // Get pointers to the elements of interest and its parent.
822  const Elem * parent = elem->parent();
823 
824  // This can't happen... Only level-0 elements have NULL
825  // parents, and no level-0 elements can be at a higher
826  // level than their neighbors!
827  libmesh_assert(parent);
828 
829  const std::unique_ptr<const Elem> my_side (elem->build_side_ptr(s));
830  const std::unique_ptr<const Elem> parent_side (parent->build_side_ptr(s));
831 
832  const unsigned int n_side_nodes = my_side->n_nodes();
833 
834  my_nodes.clear();
835  my_nodes.reserve (n_side_nodes);
836  parent_nodes.clear();
837  parent_nodes.reserve (n_side_nodes);
838 
839  for (unsigned int n=0; n != n_side_nodes; ++n)
840  my_nodes.push_back(my_side->node_ptr(n));
841 
842  for (unsigned int n=0; n != n_side_nodes; ++n)
843  parent_nodes.push_back(parent_side->node_ptr(n));
844 
845  for (unsigned int my_side_n=0;
846  my_side_n < n_side_nodes;
847  my_side_n++)
848  {
849  libmesh_assert_less (my_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
850 
851  const Node * my_node = my_nodes[my_side_n];
852 
853  // The support point of the DOF
854  const Point & support_point = *my_node;
855 
856  // Figure out where my node lies on their reference element.
857  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
858  parent_side.get(),
859  support_point);
860 
861  // Compute the parent's side shape function values.
862  for (unsigned int their_side_n=0;
863  their_side_n < n_side_nodes;
864  their_side_n++)
865  {
866  libmesh_assert_less (their_side_n, FEInterface::n_dofs(Dim-1, fe_type, parent_side->type()));
867 
868  const Node * their_node = parent_nodes[their_side_n];
869  libmesh_assert(their_node);
870 
871  const Real their_value = FEInterface::shape(Dim-1,
872  fe_type,
873  parent_side->type(),
874  their_side_n,
875  mapped_point);
876 
877  const Real their_mag = std::abs(their_value);
878 #ifdef DEBUG
879  // Protect for the case u_i ~= u_j,
880  // in which case i better equal j.
881  if (their_mag > 0.999)
882  {
883  libmesh_assert_equal_to (my_node, their_node);
884  libmesh_assert_less (std::abs(their_value - 1.), 0.001);
885  }
886  else
887 #endif
888  // To make nodal constraints useful for constructing
889  // sparsity patterns faster, we need to get EVERY
890  // POSSIBLE constraint coupling identified, even if
891  // there is no coupling in the isoparametric
892  // Lagrange case.
893  if (their_mag < 1.e-5)
894  {
895  // since we may be running this method concurrently
896  // on multiple threads we need to acquire a lock
897  // before modifying the shared constraint_row object.
898  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
899 
900  // A reference to the constraint row.
901  NodeConstraintRow & constraint_row = constraints[my_node].first;
902 
903  constraint_row.insert(std::make_pair (their_node,
904  0.));
905  }
906  // To get nodal coordinate constraints right, only
907  // add non-zero and non-identity values for Lagrange
908  // basis functions.
909  else // (1.e-5 <= their_mag <= .999)
910  {
911  // since we may be running this method concurrently
912  // on multiple threads we need to acquire a lock
913  // before modifying the shared constraint_row object.
914  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
915 
916  // A reference to the constraint row.
917  NodeConstraintRow & constraint_row = constraints[my_node].first;
918 
919  constraint_row.insert(std::make_pair (their_node,
920  their_value));
921  }
922  }
923  }
924  }
925 }
double abs(double a)
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:411
const class libmesh_nullptr_t libmesh_nullptr
spin_mutex spin_mtx
Definition: threads.C:29
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition: fe_interface.C:614
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:547
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< const Node *, Real, std::less< const Node * >, Threads::scalable_allocator< std::pair< const Node *const, Real > > > NodeConstraintRow
Definition: dof_map.h:137
const RemoteElem * remote_elem
Definition: remote_elem.C:57
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void libMesh::InfFE< Dim, T_radial, T_map >::compute_node_indices ( const ElemType  inf_elem_type,
const unsigned int  outer_node_index,
unsigned int &  base_node,
unsigned int &  radial_node 
)
staticprotected

Computes the indices in the base base_node and in radial direction radial_node (either 0 or 1) associated to the node outer_node_index of an infinite element of type inf_elem_type.

Definition at line 387 of file inf_fe_static.C.

References libMesh::INFEDGE2, libMesh::INFHEX16, libMesh::INFHEX18, libMesh::INFHEX8, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, and libMesh::INFQUAD6.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::init_base_shape_functions(), and libMesh::InfFE< Dim, T_radial, T_map >::init_shape_functions().

391 {
392  switch (inf_elem_type)
393  {
394  case INFEDGE2:
395  {
396  libmesh_assert_less (outer_node_index, 2);
397  base_node = 0;
398  radial_node = outer_node_index;
399  return;
400  }
401 
402 
403  // linear base approximation, easy to determine
404  case INFQUAD4:
405  {
406  libmesh_assert_less (outer_node_index, 4);
407  base_node = outer_node_index % 2;
408  radial_node = outer_node_index / 2;
409  return;
410  }
411 
412  case INFPRISM6:
413  {
414  libmesh_assert_less (outer_node_index, 6);
415  base_node = outer_node_index % 3;
416  radial_node = outer_node_index / 3;
417  return;
418  }
419 
420  case INFHEX8:
421  {
422  libmesh_assert_less (outer_node_index, 8);
423  base_node = outer_node_index % 4;
424  radial_node = outer_node_index / 4;
425  return;
426  }
427 
428 
429  // higher order base approximation, more work necessary
430  case INFQUAD6:
431  {
432  switch (outer_node_index)
433  {
434  case 0:
435  case 1:
436  {
437  radial_node = 0;
438  base_node = outer_node_index;
439  return;
440  }
441 
442  case 2:
443  case 3:
444  {
445  radial_node = 1;
446  base_node = outer_node_index-2;
447  return;
448  }
449 
450  case 4:
451  {
452  radial_node = 0;
453  base_node = 2;
454  return;
455  }
456 
457  case 5:
458  {
459  radial_node = 1;
460  base_node = 2;
461  return;
462  }
463 
464  default:
465  libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
466  }
467  }
468 
469 
470  case INFHEX16:
471  case INFHEX18:
472  {
473  switch (outer_node_index)
474  {
475  case 0:
476  case 1:
477  case 2:
478  case 3:
479  {
480  radial_node = 0;
481  base_node = outer_node_index;
482  return;
483  }
484 
485  case 4:
486  case 5:
487  case 6:
488  case 7:
489  {
490  radial_node = 1;
491  base_node = outer_node_index-4;
492  return;
493  }
494 
495  case 8:
496  case 9:
497  case 10:
498  case 11:
499  {
500  radial_node = 0;
501  base_node = outer_node_index-4;
502  return;
503  }
504 
505  case 12:
506  case 13:
507  case 14:
508  case 15:
509  {
510  radial_node = 1;
511  base_node = outer_node_index-8;
512  return;
513  }
514 
515  case 16:
516  {
517  libmesh_assert_equal_to (inf_elem_type, INFHEX18);
518  radial_node = 0;
519  base_node = 8;
520  return;
521  }
522 
523  case 17:
524  {
525  libmesh_assert_equal_to (inf_elem_type, INFHEX18);
526  radial_node = 1;
527  base_node = 8;
528  return;
529  }
530 
531  default:
532  libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
533  }
534  }
535 
536 
537  case INFPRISM12:
538  {
539  switch (outer_node_index)
540  {
541  case 0:
542  case 1:
543  case 2:
544  {
545  radial_node = 0;
546  base_node = outer_node_index;
547  return;
548  }
549 
550  case 3:
551  case 4:
552  case 5:
553  {
554  radial_node = 1;
555  base_node = outer_node_index-3;
556  return;
557  }
558 
559  case 6:
560  case 7:
561  case 8:
562  {
563  radial_node = 0;
564  base_node = outer_node_index-3;
565  return;
566  }
567 
568  case 9:
569  case 10:
570  case 11:
571  {
572  radial_node = 1;
573  base_node = outer_node_index-6;
574  return;
575  }
576 
577  default:
578  libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
579  }
580  }
581 
582 
583  default:
584  libmesh_error_msg("ERROR: Bad infinite element type=" << inf_elem_type << ", node=" << outer_node_index);
585  }
586 }
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void libMesh::InfFE< Dim, T_radial, T_map >::compute_node_indices_fast ( const ElemType  inf_elem_type,
const unsigned int  outer_node_index,
unsigned int &  base_node,
unsigned int &  radial_node 
)
staticprotected

Does the same as compute_node_indices(), but stores the maps for the current element type. Provided the infinite element type changes seldom, this is probably faster than using compute_node_indices () alone. This is possible since the number of nodes is not likely to change.

Definition at line 594 of file inf_fe_static.C.

References libMesh::INFEDGE2, libMesh::INFHEX16, libMesh::INFHEX18, libMesh::INFHEX8, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, libMesh::INFQUAD6, libMesh::INVALID_ELEM, libMesh::invalid_uint, and n_nodes.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::init_base_shape_functions().

598 {
599  libmesh_assert_not_equal_to (inf_elem_type, INVALID_ELEM);
600 
601  static std::vector<unsigned int> _static_base_node_index;
602  static std::vector<unsigned int> _static_radial_node_index;
603 
604  /*
605  * fast counterpart to compute_node_indices(), uses local static buffers
606  * to store the index maps. The class member
607  * \p _compute_node_indices_fast_current_elem_type remembers
608  * the current element type.
609  *
610  * Note that there exist non-static members storing the
611  * same data. However, you never know what element type
612  * is currently used by the \p InfFE object, and what
613  * request is currently directed to the static \p InfFE
614  * members (which use \p compute_node_indices_fast()).
615  * So separate these.
616  *
617  * check whether the work for this elemtype has already
618  * been done. If so, use this index. Otherwise, refresh
619  * the buffer to this element type.
620  */
622  {
623  base_node = _static_base_node_index [outer_node_index];
624  radial_node = _static_radial_node_index[outer_node_index];
625  return;
626  }
627  else
628  {
629  // store the map for _all_ nodes for this element type
631 
632  unsigned int n_nodes = libMesh::invalid_uint;
633 
634  switch (inf_elem_type)
635  {
636  case INFEDGE2:
637  {
638  n_nodes = 2;
639  break;
640  }
641  case INFQUAD4:
642  {
643  n_nodes = 4;
644  break;
645  }
646  case INFQUAD6:
647  {
648  n_nodes = 6;
649  break;
650  }
651  case INFHEX8:
652  {
653  n_nodes = 8;
654  break;
655  }
656  case INFHEX16:
657  {
658  n_nodes = 16;
659  break;
660  }
661  case INFHEX18:
662  {
663  n_nodes = 18;
664  break;
665  }
666  case INFPRISM6:
667  {
668  n_nodes = 6;
669  break;
670  }
671  case INFPRISM12:
672  {
673  n_nodes = 12;
674  break;
675  }
676  default:
677  libmesh_error_msg("ERROR: Bad infinite element type=" << inf_elem_type << ", node=" << outer_node_index);
678  }
679 
680 
681  _static_base_node_index.resize (n_nodes);
682  _static_radial_node_index.resize(n_nodes);
683 
684  for (unsigned int n=0; n<n_nodes; n++)
685  compute_node_indices (inf_elem_type,
686  n,
687  _static_base_node_index [outer_node_index],
688  _static_radial_node_index[outer_node_index]);
689 
690  // and return for the specified node
691  base_node = _static_base_node_index [outer_node_index];
692  radial_node = _static_radial_node_index[outer_node_index];
693  return;
694  }
695 }
const unsigned int invalid_uint
Definition: libmesh.h:184
static void compute_node_indices(const ElemType inf_elem_type, const unsigned int outer_node_index, unsigned int &base_node, unsigned int &radial_node)
const dof_id_type n_nodes
Definition: tecplot_io.C:67
static ElemType _compute_node_indices_fast_current_elem_type
Definition: inf_fe.h:798
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::compute_periodic_constraints ( DofConstraints constraints,
DofMap dof_map,
const PeriodicBoundaries boundaries,
const MeshBase mesh,
const PointLocatorBase point_locator,
const unsigned int  variable_number,
const Elem elem 
)
staticinherited

Computes the constraint matrix contributions (for meshes with periodic boundary conditions) corresponding to variable number var_number, using generic projections.

Definition at line 1647 of file fe_base.C.

References std::abs(), libMesh::TypeVector< T >::absolute_fuzzy_equals(), libMesh::Elem::active(), libMesh::PeriodicBoundaries::boundary(), libMesh::BoundaryInfo::boundary_ids(), libMesh::C_ONE, libMesh::C_ZERO, libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::DofMap::constrain_p_dofs(), libMesh::FEType::default_quadrature_order(), libMesh::Elem::dim(), libMesh::DISCONTINUOUS, libMesh::DofMap::dof_indices(), libMesh::DofObject::dof_number(), libMesh::FEInterface::dofs_on_side(), libMesh::MeshBase::get_boundary_info(), libMesh::PeriodicBoundaryBase::get_corresponding_pos(), libMesh::Elem::hmin(), libMesh::DofObject::id(), libMesh::TensorTools::inner_product(), libMesh::DofObject::invalid_id, libMesh::invalid_uint, libMesh::FEInterface::inverse_map(), libMesh::DofMap::is_constrained_dof(), libMesh::Elem::is_edge(), libMesh::Elem::is_face(), libMesh::PeriodicBoundaryBase::is_my_variable(), libMesh::Elem::is_node_on_edge(), libMesh::Elem::is_node_on_side(), libMesh::Elem::is_vertex(), libMesh::Elem::level(), libmesh_nullptr, std::min(), libMesh::Elem::min_p_level_by_neighbor(), libMesh::DofObject::n_comp(), libMesh::Elem::n_edges(), libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::PeriodicBoundaries::neighbor(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_ptr(), libMesh::Elem::node_ref(), libMesh::Elem::p_level(), libMesh::PeriodicBoundaryBase::pairedboundary, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::BoundaryInfo::side_with_boundary_id(), libMesh::Threads::spin_mtx, swap(), libMesh::DofMap::sys_number(), libMesh::TOLERANCE, and libMesh::DofMap::variable_type().

Referenced by libMesh::FEInterface::compute_periodic_constraints(), and libMesh::FEGenericBase< OutputType >::compute_proj_constraints().

1654 {
1655  // Only bother if we truly have periodic boundaries
1656  if (boundaries.empty())
1657  return;
1658 
1659  libmesh_assert(elem);
1660 
1661  // Only constrain active elements with this method
1662  if (!elem->active())
1663  return;
1664 
1665  const unsigned int Dim = elem->dim();
1666 
1667  // We need sys_number and variable_number for DofObject methods
1668  // later
1669  const unsigned int sys_number = dof_map.sys_number();
1670 
1671  const FEType & base_fe_type = dof_map.variable_type(variable_number);
1672 
1673  // Construct FE objects for this element and its pseudo-neighbors.
1674  std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1675  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1676  const FEContinuity cont = my_fe->get_continuity();
1677 
1678  // We don't need to constrain discontinuous elements
1679  if (cont == DISCONTINUOUS)
1680  return;
1681  libmesh_assert (cont == C_ZERO || cont == C_ONE);
1682 
1683  // We'll use element size to generate relative tolerances later
1684  const Real primary_hmin = elem->hmin();
1685 
1686  std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1687  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1688 
1689  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1690  my_fe->attach_quadrature_rule (&my_qface);
1691  std::vector<Point> neigh_qface;
1692 
1693  const std::vector<Real> & JxW = my_fe->get_JxW();
1694  const std::vector<Point> & q_point = my_fe->get_xyz();
1695  const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1696  const std::vector<std::vector<OutputShape>> & neigh_phi =
1697  neigh_fe->get_phi();
1698  const std::vector<Point> * face_normals = libmesh_nullptr;
1699  const std::vector<std::vector<OutputGradient>> * dphi = libmesh_nullptr;
1700  const std::vector<std::vector<OutputGradient>> * neigh_dphi = libmesh_nullptr;
1701  std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1702  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1703 
1704  if (cont != C_ZERO)
1705  {
1706  const std::vector<Point> & ref_face_normals =
1707  my_fe->get_normals();
1708  face_normals = &ref_face_normals;
1709  const std::vector<std::vector<OutputGradient>> & ref_dphi =
1710  my_fe->get_dphi();
1711  dphi = &ref_dphi;
1712  const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1713  neigh_fe->get_dphi();
1714  neigh_dphi = &ref_neigh_dphi;
1715  }
1716 
1717  DenseMatrix<Real> Ke;
1718  DenseVector<Real> Fe;
1719  std::vector<DenseVector<Real>> Ue;
1720 
1721  // Container to catch the boundary ids that BoundaryInfo hands us.
1722  std::vector<boundary_id_type> bc_ids;
1723 
1724  // Look at the element faces. Check to see if we need to
1725  // build constraints.
1726  const unsigned short int max_ns = elem->n_sides();
1727  for (unsigned short int s = 0; s != max_ns; ++s)
1728  {
1729  if (elem->neighbor_ptr(s))
1730  continue;
1731 
1732  mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
1733 
1734  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1735  {
1736  const boundary_id_type boundary_id = *id_it;
1737  const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
1738  if (periodic && periodic->is_my_variable(variable_number))
1739  {
1740  libmesh_assert(point_locator);
1741 
1742  // Get pointers to the element's neighbor.
1743  const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s);
1744 
1745  if (neigh == libmesh_nullptr)
1746  libmesh_error_msg("PeriodicBoundaries point locator object returned NULL!");
1747 
1748  // periodic (and possibly h refinement) constraints:
1749  // constrain dofs shared between
1750  // this element and ones as coarse
1751  // as or coarser than this element.
1752  if (neigh->level() <= elem->level())
1753  {
1754  unsigned int s_neigh =
1755  mesh.get_boundary_info().side_with_boundary_id(neigh, periodic->pairedboundary);
1756  libmesh_assert_not_equal_to (s_neigh, libMesh::invalid_uint);
1757 
1758 #ifdef LIBMESH_ENABLE_AMR
1759  // Find the minimum p level; we build the h constraint
1760  // matrix with this and then constrain away all higher p
1761  // DoFs.
1762  libmesh_assert(neigh->active());
1763  const unsigned int min_p_level =
1764  std::min(elem->p_level(), neigh->p_level());
1765 
1766  // we may need to make the FE objects reinit with the
1767  // minimum shared p_level
1768  // FIXME - I hate using const_cast<> and avoiding
1769  // accessor functions; there's got to be a
1770  // better way to do this!
1771  const unsigned int old_elem_level = elem->p_level();
1772  if (old_elem_level != min_p_level)
1773  (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
1774  const unsigned int old_neigh_level = neigh->p_level();
1775  if (old_neigh_level != min_p_level)
1776  (const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
1777 #endif // #ifdef LIBMESH_ENABLE_AMR
1778 
1779  // We can do a projection with a single integration,
1780  // due to the assumption of nested finite element
1781  // subspaces.
1782  // FIXME: it might be more efficient to do nodes,
1783  // then edges, then side, to reduce the size of the
1784  // Cholesky factorization(s)
1785  my_fe->reinit(elem, s);
1786 
1787  dof_map.dof_indices (elem, my_dof_indices,
1788  variable_number);
1789  dof_map.dof_indices (neigh, neigh_dof_indices,
1790  variable_number);
1791 
1792  const unsigned int n_qp = my_qface.n_points();
1793 
1794  // Translate the quadrature points over to the
1795  // neighbor's boundary
1796  std::vector<Point> neigh_point(q_point.size());
1797  for (std::size_t i=0; i != neigh_point.size(); ++i)
1798  neigh_point[i] = periodic->get_corresponding_pos(q_point[i]);
1799 
1800  FEInterface::inverse_map (Dim, base_fe_type, neigh,
1801  neigh_point, neigh_qface);
1802 
1803  neigh_fe->reinit(neigh, &neigh_qface);
1804 
1805  // We're only concerned with DOFs whose values (and/or first
1806  // derivatives for C1 elements) are supported on side nodes
1807  FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
1808  FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
1809 
1810  // We're done with functions that examine Elem::p_level(),
1811  // so let's unhack those levels
1812 #ifdef LIBMESH_ENABLE_AMR
1813  if (elem->p_level() != old_elem_level)
1814  (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
1815  if (neigh->p_level() != old_neigh_level)
1816  (const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
1817 #endif // #ifdef LIBMESH_ENABLE_AMR
1818 
1819  const unsigned int n_side_dofs =
1820  cast_int<unsigned int>
1821  (my_side_dofs.size());
1822  libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
1823 
1824  Ke.resize (n_side_dofs, n_side_dofs);
1825  Ue.resize(n_side_dofs);
1826 
1827  // Form the projection matrix, (inner product of fine basis
1828  // functions against fine test functions)
1829  for (unsigned int is = 0; is != n_side_dofs; ++is)
1830  {
1831  const unsigned int i = my_side_dofs[is];
1832  for (unsigned int js = 0; js != n_side_dofs; ++js)
1833  {
1834  const unsigned int j = my_side_dofs[js];
1835  for (unsigned int qp = 0; qp != n_qp; ++qp)
1836  {
1837  Ke(is,js) += JxW[qp] *
1838  TensorTools::inner_product(phi[i][qp],
1839  phi[j][qp]);
1840  if (cont != C_ZERO)
1841  Ke(is,js) += JxW[qp] *
1842  TensorTools::inner_product((*dphi)[i][qp] *
1843  (*face_normals)[qp],
1844  (*dphi)[j][qp] *
1845  (*face_normals)[qp]);
1846  }
1847  }
1848  }
1849 
1850  // Form the right hand sides, (inner product of coarse basis
1851  // functions against fine test functions)
1852  for (unsigned int is = 0; is != n_side_dofs; ++is)
1853  {
1854  const unsigned int i = neigh_side_dofs[is];
1855  Fe.resize (n_side_dofs);
1856  for (unsigned int js = 0; js != n_side_dofs; ++js)
1857  {
1858  const unsigned int j = my_side_dofs[js];
1859  for (unsigned int qp = 0; qp != n_qp; ++qp)
1860  {
1861  Fe(js) += JxW[qp] *
1862  TensorTools::inner_product(neigh_phi[i][qp],
1863  phi[j][qp]);
1864  if (cont != C_ZERO)
1865  Fe(js) += JxW[qp] *
1866  TensorTools::inner_product((*neigh_dphi)[i][qp] *
1867  (*face_normals)[qp],
1868  (*dphi)[j][qp] *
1869  (*face_normals)[qp]);
1870  }
1871  }
1872  Ke.cholesky_solve(Fe, Ue[is]);
1873  }
1874 
1875  // Make sure we're not adding recursive constraints
1876  // due to the redundancy in the way we add periodic
1877  // boundary constraints
1878  //
1879  // In order for this to work while threaded or on
1880  // distributed meshes, we need a rigorous way to
1881  // avoid recursive constraints. Here it is:
1882  //
1883  // For vertex DoFs, if there is a "prior" element
1884  // (i.e. a coarser element or an equally refined
1885  // element with a lower id) on this boundary which
1886  // contains the vertex point, then we will avoid
1887  // generating constraints; the prior element (or
1888  // something prior to it) may do so. If we are the
1889  // most prior (or "primary") element on this
1890  // boundary sharing this point, then we look at the
1891  // boundary periodic to us, we find the primary
1892  // element there, and if that primary is coarser or
1893  // equal-but-lower-id, then our vertex dofs are
1894  // constrained in terms of that element.
1895  //
1896  // For edge DoFs, if there is a coarser element
1897  // on this boundary sharing this edge, then we will
1898  // avoid generating constraints (we will be
1899  // constrained indirectly via AMR constraints
1900  // connecting us to the coarser element's DoFs). If
1901  // we are the coarsest element sharing this edge,
1902  // then we generate constraints if and only if we
1903  // are finer than the coarsest element on the
1904  // boundary periodic to us sharing the corresponding
1905  // periodic edge, or if we are at equal level but
1906  // our edge nodes have higher ids than the periodic
1907  // edge nodes (sorted from highest to lowest, then
1908  // compared lexicographically)
1909  //
1910  // For face DoFs, we generate constraints if we are
1911  // finer than our periodic neighbor, or if we are at
1912  // equal level but our element id is higher than its
1913  // element id.
1914  //
1915  // If the primary neighbor is also the current elem
1916  // (a 1-element-thick mesh) then we choose which
1917  // vertex dofs to constrain via lexicographic
1918  // ordering on point locations
1919 
1920  // FIXME: This code doesn't yet properly handle
1921  // cases where multiple different periodic BCs
1922  // intersect.
1923  std::set<dof_id_type> my_constrained_dofs;
1924 
1925  // Container to catch boundary IDs handed back by BoundaryInfo.
1926  std::vector<boundary_id_type> new_bc_ids;
1927 
1928  for (unsigned int n = 0; n != elem->n_nodes(); ++n)
1929  {
1930  if (!elem->is_node_on_side(n,s))
1931  continue;
1932 
1933  const Node & my_node = elem->node_ref(n);
1934 
1935  if (elem->is_vertex(n))
1936  {
1937  // Find all boundary ids that include this
1938  // point and have periodic boundary
1939  // conditions for this variable
1940  std::set<boundary_id_type> point_bcids;
1941 
1942  for (unsigned int new_s = 0;
1943  new_s != max_ns; ++new_s)
1944  {
1945  if (!elem->is_node_on_side(n,new_s))
1946  continue;
1947 
1948  mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
1949 
1950  for (std::vector<boundary_id_type>::const_iterator
1951  new_id_it=new_bc_ids.begin(); new_id_it!=new_bc_ids.end(); ++new_id_it)
1952  {
1953  const boundary_id_type new_boundary_id = *new_id_it;
1954  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
1955  if (new_periodic && new_periodic->is_my_variable(variable_number))
1956  {
1957  point_bcids.insert(new_boundary_id);
1958  }
1959  }
1960  }
1961 
1962  // See if this vertex has point neighbors to
1963  // defer to
1964  if (primary_boundary_point_neighbor
1965  (elem, my_node, mesh.get_boundary_info(), point_bcids)
1966  != elem)
1967  continue;
1968 
1969  // Find the complementary boundary id set
1970  std::set<boundary_id_type> point_pairedids;
1971  for (std::set<boundary_id_type>::const_iterator i =
1972  point_bcids.begin(); i != point_bcids.end(); ++i)
1973  {
1974  const boundary_id_type new_boundary_id = *i;
1975  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
1976  point_pairedids.insert(new_periodic->pairedboundary);
1977  }
1978 
1979  // What do we want to constrain against?
1980  const Elem * primary_elem = libmesh_nullptr;
1981  const Elem * main_neigh = libmesh_nullptr;
1982  Point main_pt = my_node,
1983  primary_pt = my_node;
1984 
1985  for (std::set<boundary_id_type>::const_iterator i =
1986  point_bcids.begin(); i != point_bcids.end(); ++i)
1987  {
1988  // Find the corresponding periodic point and
1989  // its primary neighbor
1990  const boundary_id_type new_boundary_id = *i;
1991  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
1992 
1993  const Point neigh_pt =
1994  new_periodic->get_corresponding_pos(my_node);
1995 
1996  // If the point is getting constrained
1997  // to itself by this PBC then we don't
1998  // generate any constraints
1999  if (neigh_pt.absolute_fuzzy_equals
2000  (my_node, primary_hmin*TOLERANCE))
2001  continue;
2002 
2003  // Otherwise we'll have a constraint in
2004  // one direction or another
2005  if (!primary_elem)
2006  primary_elem = elem;
2007 
2008  const Elem * primary_neigh =
2009  primary_boundary_point_neighbor(neigh, neigh_pt,
2010  mesh.get_boundary_info(),
2011  point_pairedids);
2012 
2013  libmesh_assert(primary_neigh);
2014 
2015  if (new_boundary_id == boundary_id)
2016  {
2017  main_neigh = primary_neigh;
2018  main_pt = neigh_pt;
2019  }
2020 
2021  // Finer elements will get constrained in
2022  // terms of coarser neighbors, not the
2023  // other way around
2024  if ((primary_neigh->level() > primary_elem->level()) ||
2025 
2026  // For equal-level elements, the one with
2027  // higher id gets constrained in terms of
2028  // the one with lower id
2029  (primary_neigh->level() == primary_elem->level() &&
2030  primary_neigh->id() > primary_elem->id()) ||
2031 
2032  // On a one-element-thick mesh, we compare
2033  // points to see what side gets constrained
2034  (primary_neigh == primary_elem &&
2035  (neigh_pt > primary_pt)))
2036  continue;
2037 
2038  primary_elem = primary_neigh;
2039  primary_pt = neigh_pt;
2040  }
2041 
2042  if (!primary_elem ||
2043  primary_elem != main_neigh ||
2044  primary_pt != main_pt)
2045  continue;
2046  }
2047  else if (elem->is_edge(n))
2048  {
2049  // Find which edge we're on
2050  unsigned int e=0;
2051  for (; e != elem->n_edges(); ++e)
2052  {
2053  if (elem->is_node_on_edge(n,e))
2054  break;
2055  }
2056  libmesh_assert_less (e, elem->n_edges());
2057 
2058  // Find the edge end nodes
2059  const Node
2060  * e1 = libmesh_nullptr,
2061  * e2 = libmesh_nullptr;
2062  for (unsigned int nn = 0; nn != elem->n_nodes(); ++nn)
2063  {
2064  if (nn == n)
2065  continue;
2066 
2067  if (elem->is_node_on_edge(nn, e))
2068  {
2069  if (e1 == libmesh_nullptr)
2070  {
2071  e1 = elem->node_ptr(nn);
2072  }
2073  else
2074  {
2075  e2 = elem->node_ptr(nn);
2076  break;
2077  }
2078  }
2079  }
2080  libmesh_assert (e1 && e2);
2081 
2082  // Find all boundary ids that include this
2083  // edge and have periodic boundary
2084  // conditions for this variable
2085  std::set<boundary_id_type> edge_bcids;
2086 
2087  for (unsigned int new_s = 0;
2088  new_s != max_ns; ++new_s)
2089  {
2090  if (!elem->is_node_on_side(n,new_s))
2091  continue;
2092 
2093  // We're reusing the new_bc_ids vector created outside the loop over nodes.
2094  mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
2095 
2096  for (std::vector<boundary_id_type>::const_iterator
2097  new_id_it=new_bc_ids.begin(); new_id_it!=new_bc_ids.end(); ++new_id_it)
2098  {
2099  const boundary_id_type new_boundary_id = *new_id_it;
2100  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2101  if (new_periodic && new_periodic->is_my_variable(variable_number))
2102  {
2103  edge_bcids.insert(new_boundary_id);
2104  }
2105  }
2106  }
2107 
2108 
2109  // See if this edge has neighbors to defer to
2110  if (primary_boundary_edge_neighbor
2111  (elem, *e1, *e2, mesh.get_boundary_info(), edge_bcids)
2112  != elem)
2113  continue;
2114 
2115  // Find the complementary boundary id set
2116  std::set<boundary_id_type> edge_pairedids;
2117  for (std::set<boundary_id_type>::const_iterator i =
2118  edge_bcids.begin(); i != edge_bcids.end(); ++i)
2119  {
2120  const boundary_id_type new_boundary_id = *i;
2121  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2122  edge_pairedids.insert(new_periodic->pairedboundary);
2123  }
2124 
2125  // What do we want to constrain against?
2126  const Elem * primary_elem = libmesh_nullptr;
2127  const Elem * main_neigh = libmesh_nullptr;
2128  Point main_pt1 = *e1,
2129  main_pt2 = *e2,
2130  primary_pt1 = *e1,
2131  primary_pt2 = *e2;
2132 
2133  for (std::set<boundary_id_type>::const_iterator i =
2134  edge_bcids.begin(); i != edge_bcids.end(); ++i)
2135  {
2136  // Find the corresponding periodic edge and
2137  // its primary neighbor
2138  const boundary_id_type new_boundary_id = *i;
2139  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2140 
2141  Point neigh_pt1 = new_periodic->get_corresponding_pos(*e1),
2142  neigh_pt2 = new_periodic->get_corresponding_pos(*e2);
2143 
2144  // If the edge is getting constrained
2145  // to itself by this PBC then we don't
2146  // generate any constraints
2147  if (neigh_pt1.absolute_fuzzy_equals
2148  (*e1, primary_hmin*TOLERANCE) &&
2149  neigh_pt2.absolute_fuzzy_equals
2150  (*e2, primary_hmin*TOLERANCE))
2151  continue;
2152 
2153  // Otherwise we'll have a constraint in
2154  // one direction or another
2155  if (!primary_elem)
2156  primary_elem = elem;
2157 
2158  const Elem * primary_neigh = primary_boundary_edge_neighbor
2159  (neigh, neigh_pt1, neigh_pt2,
2160  mesh.get_boundary_info(), edge_pairedids);
2161 
2162  libmesh_assert(primary_neigh);
2163 
2164  if (new_boundary_id == boundary_id)
2165  {
2166  main_neigh = primary_neigh;
2167  main_pt1 = neigh_pt1;
2168  main_pt2 = neigh_pt2;
2169  }
2170 
2171  // If we have a one-element thick mesh,
2172  // we'll need to sort our points to get a
2173  // consistent ordering rule
2174  //
2175  // Use >= in this test to make sure that,
2176  // for angular constraints, no node gets
2177  // constrained to itself.
2178  if (primary_neigh == primary_elem)
2179  {
2180  if (primary_pt1 > primary_pt2)
2181  std::swap(primary_pt1, primary_pt2);
2182  if (neigh_pt1 > neigh_pt2)
2183  std::swap(neigh_pt1, neigh_pt2);
2184 
2185  if (neigh_pt2 >= primary_pt2)
2186  continue;
2187  }
2188 
2189  // Otherwise:
2190  // Finer elements will get constrained in
2191  // terms of coarser ones, not the other way
2192  // around
2193  if ((primary_neigh->level() > primary_elem->level()) ||
2194 
2195  // For equal-level elements, the one with
2196  // higher id gets constrained in terms of
2197  // the one with lower id
2198  (primary_neigh->level() == primary_elem->level() &&
2199  primary_neigh->id() > primary_elem->id()))
2200  continue;
2201 
2202  primary_elem = primary_neigh;
2203  primary_pt1 = neigh_pt1;
2204  primary_pt2 = neigh_pt2;
2205  }
2206 
2207  if (!primary_elem ||
2208  primary_elem != main_neigh ||
2209  primary_pt1 != main_pt1 ||
2210  primary_pt2 != main_pt2)
2211  continue;
2212  }
2213  else if (elem->is_face(n))
2214  {
2215  // If we have a one-element thick mesh,
2216  // use the ordering of the face node and its
2217  // periodic counterpart to determine what
2218  // gets constrained
2219  if (neigh == elem)
2220  {
2221  const Point neigh_pt =
2222  periodic->get_corresponding_pos(my_node);
2223  if (neigh_pt > my_node)
2224  continue;
2225  }
2226 
2227  // Otherwise:
2228  // Finer elements will get constrained in
2229  // terms of coarser ones, not the other way
2230  // around
2231  if ((neigh->level() > elem->level()) ||
2232 
2233  // For equal-level elements, the one with
2234  // higher id gets constrained in terms of
2235  // the one with lower id
2236  (neigh->level() == elem->level() &&
2237  neigh->id() > elem->id()))
2238  continue;
2239  }
2240 
2241  // If we made it here without hitting a continue
2242  // statement, then we're at a node whose dofs
2243  // should be constrained by this element's
2244  // calculations.
2245  const unsigned int n_comp =
2246  my_node.n_comp(sys_number, variable_number);
2247 
2248  for (unsigned int i=0; i != n_comp; ++i)
2249  my_constrained_dofs.insert
2250  (my_node.dof_number
2251  (sys_number, variable_number, i));
2252  }
2253 
2254  // FIXME: old code for disambiguating periodic BCs:
2255  // this is not threadsafe nor safe to run on a
2256  // non-serialized mesh.
2257  /*
2258  std::vector<bool> recursive_constraint(n_side_dofs, false);
2259 
2260  for (unsigned int is = 0; is != n_side_dofs; ++is)
2261  {
2262  const unsigned int i = neigh_side_dofs[is];
2263  const dof_id_type their_dof_g = neigh_dof_indices[i];
2264  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2265 
2266  {
2267  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2268 
2269  if (!dof_map.is_constrained_dof(their_dof_g))
2270  continue;
2271  }
2272 
2273  DofConstraintRow & their_constraint_row =
2274  constraints[their_dof_g].first;
2275 
2276  for (unsigned int js = 0; js != n_side_dofs; ++js)
2277  {
2278  const unsigned int j = my_side_dofs[js];
2279  const dof_id_type my_dof_g = my_dof_indices[j];
2280  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2281 
2282  if (their_constraint_row.count(my_dof_g))
2283  recursive_constraint[js] = true;
2284  }
2285  }
2286  */
2287 
2288  for (unsigned int js = 0; js != n_side_dofs; ++js)
2289  {
2290  // FIXME: old code path
2291  // if (recursive_constraint[js])
2292  // continue;
2293 
2294  const unsigned int j = my_side_dofs[js];
2295  const dof_id_type my_dof_g = my_dof_indices[j];
2296  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2297 
2298  // FIXME: new code path
2299  if (!my_constrained_dofs.count(my_dof_g))
2300  continue;
2301 
2302  DofConstraintRow * constraint_row;
2303 
2304  // we may be running constraint methods concurrently
2305  // on multiple threads, so we need a lock to
2306  // ensure that this constraint is "ours"
2307  {
2308  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2309 
2310  if (dof_map.is_constrained_dof(my_dof_g))
2311  continue;
2312 
2313  constraint_row = &(constraints[my_dof_g]);
2314  libmesh_assert(constraint_row->empty());
2315  }
2316 
2317  for (unsigned int is = 0; is != n_side_dofs; ++is)
2318  {
2319  const unsigned int i = neigh_side_dofs[is];
2320  const dof_id_type their_dof_g = neigh_dof_indices[i];
2321  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2322 
2323  // Periodic constraints should never be
2324  // self-constraints
2325  // libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
2326 
2327  const Real their_dof_value = Ue[is](js);
2328 
2329  if (their_dof_g == my_dof_g)
2330  {
2331  libmesh_assert_less (std::abs(their_dof_value-1.), 1.e-5);
2332  for (unsigned int k = 0; k != n_side_dofs; ++k)
2333  libmesh_assert(k == is || std::abs(Ue[k](js)) < 1.e-5);
2334  continue;
2335  }
2336 
2337  if (std::abs(their_dof_value) < 10*TOLERANCE)
2338  continue;
2339 
2340  constraint_row->insert(std::make_pair(their_dof_g,
2341  their_dof_value));
2342  }
2343  }
2344  }
2345  // p refinement constraints:
2346  // constrain dofs shared between
2347  // active elements and neighbors with
2348  // lower polynomial degrees
2349 #ifdef LIBMESH_ENABLE_AMR
2350  const unsigned int min_p_level =
2351  neigh->min_p_level_by_neighbor(elem, elem->p_level());
2352  if (min_p_level < elem->p_level())
2353  {
2354  // Adaptive p refinement of non-hierarchic bases will
2355  // require more coding
2356  libmesh_assert(my_fe->is_hierarchic());
2357  dof_map.constrain_p_dofs(variable_number, elem,
2358  s, min_p_level);
2359  }
2360 #endif // #ifdef LIBMESH_ENABLE_AMR
2361  }
2362  }
2363  }
2364 }
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const =0
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
virtual bool is_edge(const unsigned int i) const =0
const BoundaryInfo & get_boundary_info() const
Definition: mesh_base.h:118
double abs(double a)
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
bool active() const
Definition: elem.h:2258
const unsigned int invalid_uint
Definition: libmesh.h:184
unsigned int p_level() const
Definition: elem.h:2423
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1719
virtual Point get_corresponding_pos(const Point &pt) const =0
static void dofs_on_side(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int s, std::vector< unsigned int > &di)
Definition: fe_interface.C:480
virtual unsigned int n_edges() const =0
void resize(const unsigned int n)
Definition: dense_vector.h:350
unsigned int min_p_level_by_neighbor(const Elem *neighbor, unsigned int current_min) const
Definition: elem.C:1906
The base class for all geometric element types.
Definition: elem.h:90
virtual Real hmin() const
Definition: elem.C:350
PeriodicBoundaryBase * boundary(boundary_id_type id)
const class libmesh_nullptr_t libmesh_nullptr
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
static const Real TOLERANCE
virtual unsigned int n_nodes() const =0
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1968
std::vector< boundary_id_type > boundary_ids(const Node *node) const
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:1756
unsigned int side_with_boundary_id(const Elem *const elem, const boundary_id_type boundary_id) const
spin_mutex spin_mtx
Definition: threads.C:29
std::vector< std::vector< OutputShape > > phi
Definition: fe_base.h:499
int8_t boundary_id_type
Definition: id_types.h:51
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1875
Order default_quadrature_order() const
Definition: fe_type.h:333
const Node & node_ref(const unsigned int i) const
Definition: elem.h:1897
static const dof_id_type invalid_id
Definition: dof_object.h:324
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:547
virtual unsigned int n_sides() const =0
unsigned int sys_number() const
Definition: dof_map.h:1671
std::vector< std::vector< OutputGradient > > dphi
Definition: fe_base.h:504
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool absolute_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:962
virtual bool is_vertex(const unsigned int i) const =0
void swap(Iterator &lhs, Iterator &rhs)
virtual bool is_face(const unsigned int i) const =0
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
Definition: dof_map.h:89
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
virtual unsigned int dim() const =0
unsigned int level() const
Definition: elem.h:2389
Implements 1, 2, and 3D "Gaussian" quadrature rules.
bool is_my_variable(unsigned int var_num) const
Base class for all PeriodicBoundary implementations.
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:780
dof_id_type id() const
Definition: dof_object.h:632
A matrix object used for finite element assembly and numerics.
Definition: dense_matrix.h:54
long double min(long double a, double b)
A geometric point in (x,y,z) space.
Definition: point.h:38
const Elem * neighbor(boundary_id_type boundary_id, const PointLocatorBase &point_locator, const Elem *e, unsigned int side) const
boostcopy::enable_if_c< ScalarTraits< T >::value &&ScalarTraits< T2 >::value, typename CompareTypes< T, T2 >::supertype >::type inner_product(const T &a, const T2 &b)
Definition: tensor_tools.h:47
void constrain_p_dofs(unsigned int var, const Elem *elem, unsigned int s, unsigned int p)
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1924
void libMesh::FEAbstract::compute_periodic_node_constraints ( NodeConstraints constraints,
const PeriodicBoundaries boundaries,
const MeshBase mesh,
const PointLocatorBase point_locator,
const Elem elem 
)
staticinherited

Computes the node position constraint equation contributions (for meshes with periodic boundary conditions)

Definition at line 936 of file fe_abstract.C.

References libMesh::Elem::active(), libMesh::PeriodicBoundaries::boundary(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build_side_ptr(), libMesh::Elem::default_order(), libMesh::Elem::dim(), libMesh::FEAbstract::fe_type, libMesh::MeshBase::get_boundary_info(), libMesh::PeriodicBoundaryBase::get_corresponding_pos(), libMesh::invalid_uint, libMesh::FEInterface::inverse_map(), libMesh::LAGRANGE, libMesh::Elem::level(), libMesh::FEInterface::n_dofs(), libMesh::PeriodicBoundaries::neighbor(), libMesh::Elem::neighbor_ptr(), libMesh::PeriodicBoundaryBase::pairedboundary, libMesh::Real, libMesh::FEInterface::shape(), libMesh::Elem::side_index_range(), libMesh::BoundaryInfo::side_with_boundary_id(), and libMesh::Threads::spin_mtx.

941 {
942  // Only bother if we truly have periodic boundaries
943  if (boundaries.empty())
944  return;
945 
946  libmesh_assert(elem);
947 
948  // Only constrain active elements with this method
949  if (!elem->active())
950  return;
951 
952  const unsigned int Dim = elem->dim();
953 
954  // We currently always use LAGRANGE mappings for geometry
955  const FEType fe_type(elem->default_order(), LAGRANGE);
956 
957  std::vector<const Node *> my_nodes, neigh_nodes;
958 
959  // Look at the element faces. Check to see if we need to
960  // build constraints.
961  std::vector<boundary_id_type> bc_ids;
962  for (auto s : elem->side_index_range())
963  {
964  if (elem->neighbor_ptr(s))
965  continue;
966 
967  mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
968  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
969  {
970  const boundary_id_type boundary_id = *id_it;
971  const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
972  if (periodic)
973  {
974  libmesh_assert(point_locator);
975 
976  // Get pointers to the element's neighbor.
977  const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s);
978 
979  // h refinement constraints:
980  // constrain dofs shared between
981  // this element and ones as coarse
982  // as or coarser than this element.
983  if (neigh->level() <= elem->level())
984  {
985  unsigned int s_neigh =
986  mesh.get_boundary_info().side_with_boundary_id(neigh, periodic->pairedboundary);
987  libmesh_assert_not_equal_to (s_neigh, libMesh::invalid_uint);
988 
989 #ifdef LIBMESH_ENABLE_AMR
990  libmesh_assert(neigh->active());
991 #endif // #ifdef LIBMESH_ENABLE_AMR
992 
993  const std::unique_ptr<const Elem> my_side (elem->build_side_ptr(s));
994  const std::unique_ptr<const Elem> neigh_side (neigh->build_side_ptr(s_neigh));
995 
996  const unsigned int n_side_nodes = my_side->n_nodes();
997 
998  my_nodes.clear();
999  my_nodes.reserve (n_side_nodes);
1000  neigh_nodes.clear();
1001  neigh_nodes.reserve (n_side_nodes);
1002 
1003  for (unsigned int n=0; n != n_side_nodes; ++n)
1004  my_nodes.push_back(my_side->node_ptr(n));
1005 
1006  for (unsigned int n=0; n != n_side_nodes; ++n)
1007  neigh_nodes.push_back(neigh_side->node_ptr(n));
1008 
1009  // Make sure we're not adding recursive constraints
1010  // due to the redundancy in the way we add periodic
1011  // boundary constraints, or adding constraints to
1012  // nodes that already have AMR constraints
1013  std::vector<bool> skip_constraint(n_side_nodes, false);
1014 
1015  for (unsigned int my_side_n=0;
1016  my_side_n < n_side_nodes;
1017  my_side_n++)
1018  {
1019  libmesh_assert_less (my_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
1020 
1021  const Node * my_node = my_nodes[my_side_n];
1022 
1023  // Figure out where my node lies on their reference element.
1024  const Point neigh_point = periodic->get_corresponding_pos(*my_node);
1025 
1026  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
1027  neigh_side.get(),
1028  neigh_point);
1029 
1030  // If we've already got a constraint on this
1031  // node, then the periodic constraint is
1032  // redundant
1033  {
1034  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1035 
1036  if (constraints.count(my_node))
1037  {
1038  skip_constraint[my_side_n] = true;
1039  continue;
1040  }
1041  }
1042 
1043  // Compute the neighbors's side shape function values.
1044  for (unsigned int their_side_n=0;
1045  their_side_n < n_side_nodes;
1046  their_side_n++)
1047  {
1048  libmesh_assert_less (their_side_n, FEInterface::n_dofs(Dim-1, fe_type, neigh_side->type()));
1049 
1050  const Node * their_node = neigh_nodes[their_side_n];
1051 
1052  // If there's a constraint on an opposing node,
1053  // we need to see if it's constrained by
1054  // *our side* making any periodic constraint
1055  // on us recursive
1056  {
1057  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1058 
1059  if (!constraints.count(their_node))
1060  continue;
1061 
1062  const NodeConstraintRow & their_constraint_row =
1063  constraints[their_node].first;
1064 
1065  for (unsigned int orig_side_n=0;
1066  orig_side_n < n_side_nodes;
1067  orig_side_n++)
1068  {
1069  libmesh_assert_less (orig_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
1070 
1071  const Node * orig_node = my_nodes[orig_side_n];
1072 
1073  if (their_constraint_row.count(orig_node))
1074  skip_constraint[orig_side_n] = true;
1075  }
1076  }
1077  }
1078  }
1079  for (unsigned int my_side_n=0;
1080  my_side_n < n_side_nodes;
1081  my_side_n++)
1082  {
1083  libmesh_assert_less (my_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
1084 
1085  if (skip_constraint[my_side_n])
1086  continue;
1087 
1088  const Node * my_node = my_nodes[my_side_n];
1089 
1090  // Figure out where my node lies on their reference element.
1091  const Point neigh_point = periodic->get_corresponding_pos(*my_node);
1092 
1093  // Figure out where my node lies on their reference element.
1094  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
1095  neigh_side.get(),
1096  neigh_point);
1097 
1098  for (unsigned int their_side_n=0;
1099  their_side_n < n_side_nodes;
1100  their_side_n++)
1101  {
1102  libmesh_assert_less (their_side_n, FEInterface::n_dofs(Dim-1, fe_type, neigh_side->type()));
1103 
1104  const Node * their_node = neigh_nodes[their_side_n];
1105  libmesh_assert(their_node);
1106 
1107  const Real their_value = FEInterface::shape(Dim-1,
1108  fe_type,
1109  neigh_side->type(),
1110  their_side_n,
1111  mapped_point);
1112 
1113  // since we may be running this method concurrently
1114  // on multiple threads we need to acquire a lock
1115  // before modifying the shared constraint_row object.
1116  {
1117  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1118 
1119  NodeConstraintRow & constraint_row =
1120  constraints[my_node].first;
1121 
1122  constraint_row.insert(std::make_pair(their_node,
1123  their_value));
1124  }
1125  }
1126  }
1127  }
1128  }
1129  }
1130  }
1131 }
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:411
const unsigned int invalid_uint
Definition: libmesh.h:184
MeshBase & mesh
spin_mutex spin_mtx
Definition: threads.C:29
int8_t boundary_id_type
Definition: id_types.h:51
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition: fe_interface.C:614
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:547
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< const Node *, Real, std::less< const Node * >, Threads::scalable_allocator< std::pair< const Node *const, Real > > > NodeConstraintRow
Definition: dof_map.h:137
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::compute_proj_constraints ( DofConstraints constraints,
DofMap dof_map,
const unsigned int  variable_number,
const Elem elem 
)
staticinherited

Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to variable number var_number, using generic projections.

Definition at line 1363 of file fe_base.C.

References std::abs(), libMesh::Elem::active(), libMesh::C_ONE, libMesh::C_ZERO, libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::FEGenericBase< OutputType >::compute_periodic_constraints(), libMesh::DofMap::constrain_p_dofs(), libMesh::FEType::default_quadrature_order(), libMesh::Elem::dim(), libMesh::DISCONTINUOUS, libMesh::DofMap::dof_indices(), libMesh::FEInterface::dofs_on_side(), libMesh::OrderWrapper::get_order(), libMesh::TensorTools::inner_product(), libMesh::DofObject::invalid_id, libMesh::FEInterface::inverse_map(), libMesh::DofMap::is_constrained_dof(), libMesh::Elem::level(), libmesh_nullptr, std::min(), libMesh::Elem::min_p_level_by_neighbor(), libMesh::Elem::n_neighbors(), libMesh::Elem::n_nodes(), libMesh::Elem::neighbor_ptr(), libMesh::FEType::order, libMesh::Elem::p_level(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::Elem::side_index_range(), libMesh::Threads::spin_mtx, libMesh::TOLERANCE, libMesh::DofMap::variable_type(), and libMesh::Elem::which_neighbor_am_i().

Referenced by libMesh::FE< Dim, T >::compute_constraints().

1367 {
1368  libmesh_assert(elem);
1369 
1370  const unsigned int Dim = elem->dim();
1371 
1372  // Only constrain elements in 2,3D.
1373  if (Dim == 1)
1374  return;
1375 
1376  // Only constrain active elements with this method
1377  if (!elem->active())
1378  return;
1379 
1380  const FEType & base_fe_type = dof_map.variable_type(variable_number);
1381 
1382  // Construct FE objects for this element and its neighbors.
1383  std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1384  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1385  const FEContinuity cont = my_fe->get_continuity();
1386 
1387  // We don't need to constrain discontinuous elements
1388  if (cont == DISCONTINUOUS)
1389  return;
1390  libmesh_assert (cont == C_ZERO || cont == C_ONE);
1391 
1392  std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1393  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1394 
1395  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1396  my_fe->attach_quadrature_rule (&my_qface);
1397  std::vector<Point> neigh_qface;
1398 
1399  const std::vector<Real> & JxW = my_fe->get_JxW();
1400  const std::vector<Point> & q_point = my_fe->get_xyz();
1401  const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1402  const std::vector<std::vector<OutputShape>> & neigh_phi =
1403  neigh_fe->get_phi();
1404  const std::vector<Point> * face_normals = libmesh_nullptr;
1405  const std::vector<std::vector<OutputGradient>> * dphi = libmesh_nullptr;
1406  const std::vector<std::vector<OutputGradient>> * neigh_dphi = libmesh_nullptr;
1407 
1408  std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1409  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1410 
1411  if (cont != C_ZERO)
1412  {
1413  const std::vector<Point> & ref_face_normals =
1414  my_fe->get_normals();
1415  face_normals = &ref_face_normals;
1416  const std::vector<std::vector<OutputGradient>> & ref_dphi =
1417  my_fe->get_dphi();
1418  dphi = &ref_dphi;
1419  const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1420  neigh_fe->get_dphi();
1421  neigh_dphi = &ref_neigh_dphi;
1422  }
1423 
1424  DenseMatrix<Real> Ke;
1425  DenseVector<Real> Fe;
1426  std::vector<DenseVector<Real>> Ue;
1427 
1428  // Look at the element faces. Check to see if we need to
1429  // build constraints.
1430  for (auto s : elem->side_index_range())
1431  if (elem->neighbor_ptr(s) != libmesh_nullptr)
1432  {
1433  // Get pointers to the element's neighbor.
1434  const Elem * neigh = elem->neighbor_ptr(s);
1435 
1436  // h refinement constraints:
1437  // constrain dofs shared between
1438  // this element and ones coarser
1439  // than this element.
1440  if (neigh->level() < elem->level())
1441  {
1442  unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
1443  libmesh_assert_less (s_neigh, neigh->n_neighbors());
1444 
1445  // Find the minimum p level; we build the h constraint
1446  // matrix with this and then constrain away all higher p
1447  // DoFs.
1448  libmesh_assert(neigh->active());
1449  const unsigned int min_p_level =
1450  std::min(elem->p_level(), neigh->p_level());
1451 
1452  // we may need to make the FE objects reinit with the
1453  // minimum shared p_level
1454  const unsigned int old_elem_level = elem->p_level();
1455  if (elem->p_level() != min_p_level)
1456  my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() - old_elem_level + min_p_level);
1457  const unsigned int old_neigh_level = neigh->p_level();
1458  if (old_neigh_level != min_p_level)
1459  neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() - old_neigh_level + min_p_level);
1460 
1461  my_fe->reinit(elem, s);
1462 
1463  // This function gets called element-by-element, so there
1464  // will be a lot of memory allocation going on. We can
1465  // at least minimize this for the case of the dof indices
1466  // by efficiently preallocating the requisite storage.
1467  // n_nodes is not necessarily n_dofs, but it is better
1468  // than nothing!
1469  my_dof_indices.reserve (elem->n_nodes());
1470  neigh_dof_indices.reserve (neigh->n_nodes());
1471 
1472  dof_map.dof_indices (elem, my_dof_indices,
1473  variable_number,
1474  min_p_level);
1475  dof_map.dof_indices (neigh, neigh_dof_indices,
1476  variable_number,
1477  min_p_level);
1478 
1479  const unsigned int n_qp = my_qface.n_points();
1480 
1481  FEInterface::inverse_map (Dim, base_fe_type, neigh,
1482  q_point, neigh_qface);
1483 
1484  neigh_fe->reinit(neigh, &neigh_qface);
1485 
1486  // We're only concerned with DOFs whose values (and/or first
1487  // derivatives for C1 elements) are supported on side nodes
1488  FEType elem_fe_type = base_fe_type;
1489  if (old_elem_level != min_p_level)
1490  elem_fe_type.order = base_fe_type.order.get_order() - old_elem_level + min_p_level;
1491  FEType neigh_fe_type = base_fe_type;
1492  if (old_neigh_level != min_p_level)
1493  neigh_fe_type.order = base_fe_type.order.get_order() - old_neigh_level + min_p_level;
1494  FEInterface::dofs_on_side(elem, Dim, elem_fe_type, s, my_side_dofs);
1495  FEInterface::dofs_on_side(neigh, Dim, neigh_fe_type, s_neigh, neigh_side_dofs);
1496 
1497  const unsigned int n_side_dofs =
1498  cast_int<unsigned int>(my_side_dofs.size());
1499  libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
1500 
1501  Ke.resize (n_side_dofs, n_side_dofs);
1502  Ue.resize(n_side_dofs);
1503 
1504  // Form the projection matrix, (inner product of fine basis
1505  // functions against fine test functions)
1506  for (unsigned int is = 0; is != n_side_dofs; ++is)
1507  {
1508  const unsigned int i = my_side_dofs[is];
1509  for (unsigned int js = 0; js != n_side_dofs; ++js)
1510  {
1511  const unsigned int j = my_side_dofs[js];
1512  for (unsigned int qp = 0; qp != n_qp; ++qp)
1513  {
1514  Ke(is,js) += JxW[qp] * TensorTools::inner_product(phi[i][qp], phi[j][qp]);
1515  if (cont != C_ZERO)
1516  Ke(is,js) += JxW[qp] *
1517  TensorTools::inner_product((*dphi)[i][qp] *
1518  (*face_normals)[qp],
1519  (*dphi)[j][qp] *
1520  (*face_normals)[qp]);
1521  }
1522  }
1523  }
1524 
1525  // Form the right hand sides, (inner product of coarse basis
1526  // functions against fine test functions)
1527  for (unsigned int is = 0; is != n_side_dofs; ++is)
1528  {
1529  const unsigned int i = neigh_side_dofs[is];
1530  Fe.resize (n_side_dofs);
1531  for (unsigned int js = 0; js != n_side_dofs; ++js)
1532  {
1533  const unsigned int j = my_side_dofs[js];
1534  for (unsigned int qp = 0; qp != n_qp; ++qp)
1535  {
1536  Fe(js) += JxW[qp] *
1537  TensorTools::inner_product(neigh_phi[i][qp],
1538  phi[j][qp]);
1539  if (cont != C_ZERO)
1540  Fe(js) += JxW[qp] *
1541  TensorTools::inner_product((*neigh_dphi)[i][qp] *
1542  (*face_normals)[qp],
1543  (*dphi)[j][qp] *
1544  (*face_normals)[qp]);
1545  }
1546  }
1547  Ke.cholesky_solve(Fe, Ue[is]);
1548  }
1549 
1550  for (unsigned int js = 0; js != n_side_dofs; ++js)
1551  {
1552  const unsigned int j = my_side_dofs[js];
1553  const dof_id_type my_dof_g = my_dof_indices[j];
1554  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
1555 
1556  // Hunt for "constraining against myself" cases before
1557  // we bother creating a constraint row
1558  bool self_constraint = false;
1559  for (unsigned int is = 0; is != n_side_dofs; ++is)
1560  {
1561  const unsigned int i = neigh_side_dofs[is];
1562  const dof_id_type their_dof_g = neigh_dof_indices[i];
1563  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1564 
1565  if (their_dof_g == my_dof_g)
1566  {
1567 #ifndef NDEBUG
1568  const Real their_dof_value = Ue[is](js);
1569  libmesh_assert_less (std::abs(their_dof_value-1.),
1570  10*TOLERANCE);
1571 
1572  for (unsigned int k = 0; k != n_side_dofs; ++k)
1573  libmesh_assert(k == is ||
1574  std::abs(Ue[k](js)) <
1575  10*TOLERANCE);
1576 #endif
1577 
1578  self_constraint = true;
1579  break;
1580  }
1581  }
1582 
1583  if (self_constraint)
1584  continue;
1585 
1586  DofConstraintRow * constraint_row;
1587 
1588  // we may be running constraint methods concurrently
1589  // on multiple threads, so we need a lock to
1590  // ensure that this constraint is "ours"
1591  {
1592  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1593 
1594  if (dof_map.is_constrained_dof(my_dof_g))
1595  continue;
1596 
1597  constraint_row = &(constraints[my_dof_g]);
1598  libmesh_assert(constraint_row->empty());
1599  }
1600 
1601  for (unsigned int is = 0; is != n_side_dofs; ++is)
1602  {
1603  const unsigned int i = neigh_side_dofs[is];
1604  const dof_id_type their_dof_g = neigh_dof_indices[i];
1605  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1606  libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
1607 
1608  const Real their_dof_value = Ue[is](js);
1609 
1610  if (std::abs(their_dof_value) < 10*TOLERANCE)
1611  continue;
1612 
1613  constraint_row->insert(std::make_pair(their_dof_g,
1614  their_dof_value));
1615  }
1616  }
1617 
1618  my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + old_elem_level - min_p_level);
1619  neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + old_neigh_level - min_p_level);
1620  }
1621 
1622  // p refinement constraints:
1623  // constrain dofs shared between
1624  // active elements and neighbors with
1625  // lower polynomial degrees
1626  const unsigned int min_p_level =
1627  neigh->min_p_level_by_neighbor(elem, elem->p_level());
1628  if (min_p_level < elem->p_level())
1629  {
1630  // Adaptive p refinement of non-hierarchic bases will
1631  // require more coding
1632  libmesh_assert(my_fe->is_hierarchic());
1633  dof_map.constrain_p_dofs(variable_number, elem,
1634  s, min_p_level);
1635  }
1636  }
1637 }
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
int get_order() const
Definition: fe_type.h:78
double abs(double a)
bool active() const
Definition: elem.h:2258
unsigned int p_level() const
Definition: elem.h:2423
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1719
static void dofs_on_side(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int s, std::vector< unsigned int > &di)
Definition: fe_interface.C:480
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2084
void resize(const unsigned int n)
Definition: dense_vector.h:350
unsigned int which_neighbor_am_i(const Elem *e) const
Definition: elem.h:2182
unsigned int min_p_level_by_neighbor(const Elem *neighbor, unsigned int current_min) const
Definition: elem.C:1906
The base class for all geometric element types.
Definition: elem.h:90
unsigned int n_neighbors() const
Definition: elem.h:614
const class libmesh_nullptr_t libmesh_nullptr
OrderWrapper order
Definition: fe_type.h:198
static const Real TOLERANCE
virtual unsigned int n_nodes() const =0
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1968
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:1756
spin_mutex spin_mtx
Definition: threads.C:29
std::vector< std::vector< OutputShape > > phi
Definition: fe_base.h:499
Order default_quadrature_order() const
Definition: fe_type.h:333
static const dof_id_type invalid_id
Definition: dof_object.h:324
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:547
std::vector< std::vector< OutputGradient > > dphi
Definition: fe_base.h:504
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
Definition: dof_map.h:89
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
virtual unsigned int dim() const =0
unsigned int level() const
Definition: elem.h:2389
Implements 1, 2, and 3D "Gaussian" quadrature rules.
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
A matrix object used for finite element assembly and numerics.
Definition: dense_matrix.h:54
long double min(long double a, double b)
boostcopy::enable_if_c< ScalarTraits< T >::value &&ScalarTraits< T2 >::value, typename CompareTypes< T, T2 >::supertype >::type inner_product(const T &a, const T2 &b)
Definition: tensor_tools.h:47
void constrain_p_dofs(unsigned int var, const Elem *elem, unsigned int s, unsigned int p)
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1924
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void libMesh::InfFE< Dim, T_radial, T_map >::compute_shape_functions ( const Elem ,
const std::vector< Point > &   
)
protectedvirtual

After having updated the jacobian and the transformation from local to global coordinates in FEAbstract::compute_map(), the first derivatives of the shape functions are transformed to global coordinates, giving dphi, dphidx/y/z, dphasedx/y/z, dweight. This method should barely be re-defined in derived classes, but still should be usable for children. Therefore, keep it protected.

Reimplemented from libMesh::FEGenericBase< OutputType >.

Definition at line 860 of file inf_fe.C.

References libMesh::FEAbstract::_fe_map, libMesh::InfFE< Dim, T_radial, T_map >::_n_total_qp, libMesh::FEAbstract::dim, libMesh::FEGenericBase< OutputType >::dphase, libMesh::InfFE< Dim, T_radial, T_map >::dphasedeta, libMesh::InfFE< Dim, T_radial, T_map >::dphasedxi, libMesh::InfFE< Dim, T_radial, T_map >::dphasedzeta, libMesh::FEGenericBase< OutputType >::dphi, libMesh::FEGenericBase< OutputType >::dphideta, libMesh::FEGenericBase< OutputType >::dphidx, libMesh::FEGenericBase< OutputType >::dphidxi, libMesh::FEGenericBase< OutputType >::dphidy, libMesh::FEGenericBase< OutputType >::dphidz, libMesh::FEGenericBase< OutputType >::dphidzeta, libMesh::FEGenericBase< OutputType >::dweight, libMesh::InfFE< Dim, T_radial, T_map >::dweightdv, and libMesh::FEGenericBase< OutputType >::phi.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::init_base_shape_functions(), and libMesh::InfFE< Dim, T_radial, T_map >::reinit().

862 {
863  // Start logging the overall computation of shape functions
864  LOG_SCOPE("compute_shape_functions()", "InfFE");
865 
866  const unsigned int n_total_qp = _n_total_qp;
867 
868  // Compute the shape function values (and derivatives)
869  // at the Quadrature points. Note that the actual values
870  // have already been computed via init_shape_functions
871 
872  // Compute the value of the derivative shape function i at quadrature point p
873  switch (dim)
874  {
875 
876  case 1:
877  {
878  libmesh_not_implemented();
879  break;
880  }
881 
882  case 2:
883  {
884  libmesh_not_implemented();
885  break;
886  }
887 
888  case 3:
889  {
890  const std::vector<Real> & dxidx_map = this->_fe_map->get_dxidx();
891  const std::vector<Real> & dxidy_map = this->_fe_map->get_dxidy();
892  const std::vector<Real> & dxidz_map = this->_fe_map->get_dxidz();
893 
894  const std::vector<Real> & detadx_map = this->_fe_map->get_detadx();
895  const std::vector<Real> & detady_map = this->_fe_map->get_detady();
896  const std::vector<Real> & detadz_map = this->_fe_map->get_detadz();
897 
898  const std::vector<Real> & dzetadx_map = this->_fe_map->get_dzetadx();
899  const std::vector<Real> & dzetady_map = this->_fe_map->get_dzetady();
900  const std::vector<Real> & dzetadz_map = this->_fe_map->get_dzetadz();
901 
902  // These are _all_ shape functions of this infinite element
903  for (std::size_t i=0; i<phi.size(); i++)
904  for (unsigned int p=0; p<n_total_qp; p++)
905  {
906  // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
907  dphi[i][p](0) =
908  dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
909  dphideta[i][p]*detadx_map[p] +
910  dphidzeta[i][p]*dzetadx_map[p]);
911 
912  // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
913  dphi[i][p](1) =
914  dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
915  dphideta[i][p]*detady_map[p] +
916  dphidzeta[i][p]*dzetady_map[p]);
917 
918  // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
919  dphi[i][p](2) =
920  dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
921  dphideta[i][p]*detadz_map[p] +
922  dphidzeta[i][p]*dzetadz_map[p]);
923  }
924 
925 
926  // This is the derivative of the phase term of this infinite element
927  for (unsigned int p=0; p<n_total_qp; p++)
928  {
929  // the derivative of the phase term
930  dphase[p](0) = (dphasedxi[p] * dxidx_map[p] +
931  dphasedeta[p] * detadx_map[p] +
932  dphasedzeta[p] * dzetadx_map[p]);
933 
934  dphase[p](1) = (dphasedxi[p] * dxidy_map[p] +
935  dphasedeta[p] * detady_map[p] +
936  dphasedzeta[p] * dzetady_map[p]);
937 
938  dphase[p](2) = (dphasedxi[p] * dxidz_map[p] +
939  dphasedeta[p] * detadz_map[p] +
940  dphasedzeta[p] * dzetadz_map[p]);
941 
942  // the derivative of the radial weight - varies only in radial direction,
943  // therefore dweightdxi = dweightdeta = 0.
944  dweight[p](0) = dweightdv[p] * dzetadx_map[p];
945 
946  dweight[p](1) = dweightdv[p] * dzetady_map[p];
947 
948  dweight[p](2) = dweightdv[p] * dzetadz_map[p];
949  }
950 
951  break;
952  }
953 
954  default:
955  libmesh_error_msg("Unsupported dim = " << dim);
956  }
957 }
std::vector< std::vector< OutputShape > > dphidxi
Definition: fe_base.h:519
std::vector< std::vector< OutputShape > > dphidzeta
Definition: fe_base.h:529
std::vector< Real > dphasedxi
Definition: inf_fe.h:663
std::vector< Real > dweightdv
Definition: inf_fe.h:619
unsigned int _n_total_qp
Definition: inf_fe.h:739
std::vector< std::vector< OutputShape > > dphidy
Definition: fe_base.h:539
std::vector< std::vector< OutputShape > > dphidx
Definition: fe_base.h:534
std::vector< std::vector< OutputShape > > phi
Definition: fe_base.h:499
const unsigned int dim
Definition: fe_abstract.h:523
std::vector< OutputGradient > dphase
Definition: fe_base.h:630
std::vector< Real > dphasedzeta
Definition: inf_fe.h:677
std::vector< Real > dphasedeta
Definition: inf_fe.h:670
std::vector< std::vector< OutputGradient > > dphi
Definition: fe_base.h:504
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:517
std::vector< std::vector< OutputShape > > dphidz
Definition: fe_base.h:544
std::vector< std::vector< OutputShape > > dphideta
Definition: fe_base.h:524
std::vector< RealGradient > dweight
Definition: fe_base.h:637
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void libMesh::InfFE< Dim, T_radial, T_map >::compute_shape_indices ( const FEType fet,
const ElemType  inf_elem_type,
const unsigned int  i,
unsigned int &  base_shape,
unsigned int &  radial_shape 
)
staticprotected

Computes the indices of shape functions in the base base_shape and in radial direction radial_shape (0 in the base, $ \ge 1 $ further out) associated to the shape with global index i of an infinite element of type inf_elem_type.

Definition at line 703 of file inf_fe_static.C.

References libMesh::CARTESIAN, libMesh::OrderWrapper::get_order(), libMesh::INFEDGE2, libMesh::INFHEX16, libMesh::INFHEX18, libMesh::INFHEX8, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, libMesh::INFQUAD6, libMesh::INSTANTIATE_INF_FE_MBRF(), libMesh::invalid_uint, libMesh::FEInterface::n_dofs_at_node(), libMesh::FEInterface::n_dofs_per_elem(), libMesh::FEType::radial_order, and libMesh::Real.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::init_base_shape_functions(), and libMesh::InfFE< Dim, T_radial, T_map >::init_shape_functions().

708 {
709 
710  /*
711  * An example is provided: the numbers in comments refer to
712  * a fictitious InfHex18. The numbers are chosen as exemplary
713  * values. There is currently no base approximation that
714  * requires this many dof's at nodes, sides, faces and in the element.
715  *
716  * the order of the shape functions is heavily related with the
717  * order the dofs are assigned in \p DofMap::distributed_dofs().
718  * Due to the infinite elements with higher-order base approximation,
719  * some more effort is necessary.
720  *
721  * numbering scheme:
722  * 1. all vertices in the base, assign node->n_comp() dofs to each vertex
723  * 2. all vertices further out: innermost loop: radial shapes,
724  * then the base approximation shapes
725  * 3. all side nodes in the base, assign node->n_comp() dofs to each side node
726  * 4. all side nodes further out: innermost loop: radial shapes,
727  * then the base approximation shapes
728  * 5. (all) face nodes in the base, assign node->n_comp() dofs to each face node
729  * 6. (all) face nodes further out: innermost loop: radial shapes,
730  * then the base approximation shapes
731  * 7. element-associated dof in the base
732  * 8. element-associated dof further out
733  */
734 
735  const unsigned int radial_order = static_cast<unsigned int>(fet.radial_order.get_order()); // 4
736  const unsigned int radial_order_p_one = radial_order+1; // 5
737 
738  const ElemType base_elem_type (Base::get_elem_type(inf_elem_type)); // QUAD9
739 
740  // assume that the number of dof is the same for all vertices
741  unsigned int n_base_vertices = libMesh::invalid_uint; // 4
742  const unsigned int n_base_vertex_dof = FEInterface::n_dofs_at_node (Dim-1, fet, base_elem_type, 0);// 2
743 
744  unsigned int n_base_side_nodes = libMesh::invalid_uint; // 4
745  unsigned int n_base_side_dof = libMesh::invalid_uint; // 3
746 
747  unsigned int n_base_face_nodes = libMesh::invalid_uint; // 1
748  unsigned int n_base_face_dof = libMesh::invalid_uint; // 5
749 
750  const unsigned int n_base_elem_dof = FEInterface::n_dofs_per_elem (Dim-1, fet, base_elem_type);// 9
751 
752 
753  switch (inf_elem_type)
754  {
755  case INFEDGE2:
756  {
757  n_base_vertices = 1;
758  n_base_side_nodes = 0;
759  n_base_face_nodes = 0;
760  n_base_side_dof = 0;
761  n_base_face_dof = 0;
762  break;
763  }
764 
765  case INFQUAD4:
766  {
767  n_base_vertices = 2;
768  n_base_side_nodes = 0;
769  n_base_face_nodes = 0;
770  n_base_side_dof = 0;
771  n_base_face_dof = 0;
772  break;
773  }
774 
775  case INFQUAD6:
776  {
777  n_base_vertices = 2;
778  n_base_side_nodes = 1;
779  n_base_face_nodes = 0;
780  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
781  n_base_face_dof = 0;
782  break;
783  }
784 
785  case INFHEX8:
786  {
787  n_base_vertices = 4;
788  n_base_side_nodes = 0;
789  n_base_face_nodes = 0;
790  n_base_side_dof = 0;
791  n_base_face_dof = 0;
792  break;
793  }
794 
795  case INFHEX16:
796  {
797  n_base_vertices = 4;
798  n_base_side_nodes = 4;
799  n_base_face_nodes = 0;
800  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
801  n_base_face_dof = 0;
802  break;
803  }
804 
805  case INFHEX18:
806  {
807  n_base_vertices = 4;
808  n_base_side_nodes = 4;
809  n_base_face_nodes = 1;
810  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
811  n_base_face_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, 8);
812  break;
813  }
814 
815 
816  case INFPRISM6:
817  {
818  n_base_vertices = 3;
819  n_base_side_nodes = 0;
820  n_base_face_nodes = 0;
821  n_base_side_dof = 0;
822  n_base_face_dof = 0;
823  break;
824  }
825 
826  case INFPRISM12:
827  {
828  n_base_vertices = 3;
829  n_base_side_nodes = 3;
830  n_base_face_nodes = 0;
831  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
832  n_base_face_dof = 0;
833  break;
834  }
835 
836  default:
837  libmesh_error_msg("Unrecognized inf_elem_type = " << inf_elem_type);
838  }
839 
840 
841  {
842  // these are the limits describing the intervals where the shape function lies
843  const unsigned int n_dof_at_base_vertices = n_base_vertices*n_base_vertex_dof; // 8
844  const unsigned int n_dof_at_all_vertices = n_dof_at_base_vertices*radial_order_p_one; // 40
845 
846  const unsigned int n_dof_at_base_sides = n_base_side_nodes*n_base_side_dof; // 12
847  const unsigned int n_dof_at_all_sides = n_dof_at_base_sides*radial_order_p_one; // 60
848 
849  const unsigned int n_dof_at_base_face = n_base_face_nodes*n_base_face_dof; // 5
850  const unsigned int n_dof_at_all_faces = n_dof_at_base_face*radial_order_p_one; // 25
851 
852 
853  // start locating the shape function
854  if (i < n_dof_at_base_vertices) // range of i: 0..7
855  {
856  // belongs to vertex in the base
857  radial_shape = 0;
858  base_shape = i;
859  }
860 
861  else if (i < n_dof_at_all_vertices) // range of i: 8..39
862  {
863  /* belongs to vertex in the outer shell
864  *
865  * subtract the number of dof already counted,
866  * so that i_offset contains only the offset for the base
867  */
868  const unsigned int i_offset = i - n_dof_at_base_vertices; // 0..31
869 
870  // first the radial dof are counted, then the base dof
871  radial_shape = (i_offset % radial_order) + 1;
872  base_shape = i_offset / radial_order;
873  }
874 
875  else if (i < n_dof_at_all_vertices+n_dof_at_base_sides) // range of i: 40..51
876  {
877  // belongs to base, is a side node
878  radial_shape = 0;
879  base_shape = i - radial_order * n_dof_at_base_vertices; // 8..19
880  }
881 
882  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides) // range of i: 52..99
883  {
884  // belongs to side node in the outer shell
885  const unsigned int i_offset = i - (n_dof_at_all_vertices
886  + n_dof_at_base_sides); // 0..47
887  radial_shape = (i_offset % radial_order) + 1;
888  base_shape = (i_offset / radial_order) + n_dof_at_base_vertices;
889  }
890 
891  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_base_face) // range of i: 100..104
892  {
893  // belongs to the node in the base face
894  radial_shape = 0;
895  base_shape = i - radial_order*(n_dof_at_base_vertices
896  + n_dof_at_base_sides); // 20..24
897  }
898 
899  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces) // range of i: 105..124
900  {
901  // belongs to the node in the outer face
902  const unsigned int i_offset = i - (n_dof_at_all_vertices
903  + n_dof_at_all_sides
904  + n_dof_at_base_face); // 0..19
905  radial_shape = (i_offset % radial_order) + 1;
906  base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides;
907  }
908 
909  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces+n_base_elem_dof) // range of i: 125..133
910  {
911  // belongs to the base and is an element associated shape
912  radial_shape = 0;
913  base_shape = i - (n_dof_at_all_vertices
914  + n_dof_at_all_sides
915  + n_dof_at_all_faces); // 0..8
916  }
917 
918  else // range of i: 134..169
919  {
920  libmesh_assert_less (i, n_dofs(fet, inf_elem_type));
921  // belongs to the outer shell and is an element associated shape
922  const unsigned int i_offset = i - (n_dof_at_all_vertices
923  + n_dof_at_all_sides
924  + n_dof_at_all_faces
925  + n_base_elem_dof); // 0..19
926  radial_shape = (i_offset % radial_order) + 1;
927  base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides + n_dof_at_base_face;
928  }
929  }
930 
931  return;
932 }
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:461
const unsigned int invalid_uint
Definition: libmesh.h:184
static ElemType get_elem_type(const ElemType type)
static unsigned int n_dofs(const FEType &fet, const ElemType inf_elem_type)
Definition: inf_fe_static.C:55
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition: fe_interface.C:430
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::determine_calculations ( )
protectedinherited

Determine which values are to be calculated, for both the FE itself and for the FEMap.

Definition at line 731 of file fe_base.C.

References libMesh::FEInterface::field_type(), and libMesh::TYPE_VECTOR.

Referenced by libMesh::FEGenericBase< FEOutputType< T >::type >::get_Sobolev_dweight(), and libMesh::InfFE< Dim, T_radial, T_map >::reinit().

732 {
733  this->calculations_started = true;
734 
735  // If the user forgot to request anything, we'll be safe and
736  // calculate everything:
737 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
738  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_d2phi
739  && !this->calculate_curl_phi && !this->calculate_div_phi)
740  {
741  this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = this->calculate_dphiref = true;
743  {
744  this->calculate_curl_phi = true;
745  this->calculate_div_phi = true;
746  }
747  }
748 #else
749  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_curl_phi && !this->calculate_div_phi)
750  {
751  this->calculate_phi = this->calculate_dphi = this->calculate_dphiref = true;
753  {
754  this->calculate_curl_phi = true;
755  this->calculate_div_phi = true;
756  }
757  }
758 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
759 
760  // Request whichever terms are necessary from the FEMap
761  if (this->calculate_phi)
762  this->_fe_trans->init_map_phi(*this);
763 
764  if (this->calculate_dphiref)
765  this->_fe_trans->init_map_dphi(*this);
766 
767 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
768  if (this->calculate_d2phi)
769  this->_fe_trans->init_map_d2phi(*this);
770 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
771 }
FEFamily family
Definition: fe_type.h:204
std::unique_ptr< FETransformationBase< OutputType > > _fe_trans
Definition: fe_base.h:494
static FEFieldType field_type(const FEType &fe_type)
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited
template<unsigned int Dim, FEFamily T_radial, InfMapType T_base>
void libMesh::InfFE< Dim, T_radial, T_base >::edge_reinit ( const Elem elem,
const unsigned int  edge,
const Real  tolerance = TOLERANCE,
const std::vector< Point > *const  pts = libmesh_nullptr,
const std::vector< Real > *const  weights = libmesh_nullptr 
)
virtual

Not implemented yet. Reinitializes all the physical element-dependent data based on the edge of an infinite element.

Implements libMesh::FEAbstract.

Definition at line 104 of file inf_fe_boundary.C.

References libmesh_nullptr.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::is_hierarchic().

109 {
110  // We don't do this for 1D elements!
111  //libmesh_assert_not_equal_to (Dim, 1);
112  libmesh_not_implemented_msg("ERROR: Edge conditions for infinite elements not implemented!");
113 
114  if (pts != libmesh_nullptr)
115  libmesh_not_implemented_msg("ERROR: User-specified points for infinite elements not implemented!");
116 }
const class libmesh_nullptr_t libmesh_nullptr
void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

Methods to enable/disable the reference counter output from print_info()

Definition at line 101 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::ReferenceCounter::n_objects().

102 {
103  _enable_print_counter = true;
104  return;
105 }
template<>
Real libMesh::InfFE< 1, INFINITE_MAP, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 69 of file inf_fe_map_eval.C.

69 { return infinite_map_eval(v, i); }
template<>
Real libMesh::InfFE< 2, INFINITE_MAP, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 70 of file inf_fe_map_eval.C.

70 { return infinite_map_eval(v, i); }
template<>
Real libMesh::InfFE< 3, INFINITE_MAP, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 71 of file inf_fe_map_eval.C.

71 { return infinite_map_eval(v, i); }
template<>
Real libMesh::InfFE< 1, LEGENDRE, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 310 of file inf_fe_legendre_eval.C.

310 { return legendre_eval(v, i); }
template<>
Real libMesh::InfFE< 2, LEGENDRE, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 311 of file inf_fe_legendre_eval.C.

311 { return legendre_eval(v, i); }
template<>
Real libMesh::InfFE< 3, LEGENDRE, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 312 of file inf_fe_legendre_eval.C.

312 { return legendre_eval(v, i); }
template<>
Real libMesh::InfFE< 1, JACOBI_30_00, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 457 of file inf_fe_jacobi_30_00_eval.C.

457 { return jacobi_30_00_eval(v, i); }
template<>
Real libMesh::InfFE< 2, JACOBI_30_00, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 458 of file inf_fe_jacobi_30_00_eval.C.

458 { return jacobi_30_00_eval(v, i); }
template<>
Real libMesh::InfFE< 3, JACOBI_30_00, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 459 of file inf_fe_jacobi_30_00_eval.C.

459 { return jacobi_30_00_eval(v, i); }
template<>
Real libMesh::InfFE< 1, JACOBI_20_00, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 462 of file inf_fe_jacobi_20_00_eval.C.

462 { return jacobi_20_00_eval(v, i); }
template<>
Real libMesh::InfFE< 2, JACOBI_20_00, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 463 of file inf_fe_jacobi_20_00_eval.C.

463 { return jacobi_20_00_eval(v, i); }
template<>
Real libMesh::InfFE< 3, JACOBI_20_00, CARTESIAN >::eval ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 464 of file inf_fe_jacobi_20_00_eval.C.

464 { return jacobi_20_00_eval(v, i); }
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
static Real libMesh::InfFE< Dim, T_radial, T_map >::eval ( Real  v,
Order  o_radial,
unsigned int  i 
)
staticprotected
Returns
The value of the $ i^{th} $ polynomial evaluated at v. This method provides the approximation in radial direction for the overall shape functions, which is defined in InfFE::shape(). This method is allowed to be static, since it is independent of dimension and base_family. It is templated, though, w.r.t. to radial FEFamily.
The value of the $ i^{th} $ mapping shape function in radial direction evaluated at v when T_radial == INFINITE_MAP. Currently, only one specific mapping shape is used. Namely the one by Marques JMMC, Owen DRJ: Infinite elements in quasi-static materially nonlinear problems, Computers and Structures, 1984.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::compute_data(), libMesh::InfFE< Dim, T_radial, T_map >::init_radial_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::inverse_map(), libMesh::InfFE< Dim, T_radial, T_map >::n_quadrature_points(), and libMesh::InfFE< Dim, T_radial, T_map >::shape().

template<>
Real libMesh::InfFE< 1, LAGRANGE, CARTESIAN >::eval ( Real  v,
Order  o,
unsigned  i 
)
protected

Definition at line 2607 of file inf_fe_lagrange_eval.C.

2607 { return lagrange_eval(v, o, i); }
template<>
Real libMesh::InfFE< 2, LAGRANGE, CARTESIAN >::eval ( Real  v,
Order  o,
unsigned  i 
)
protected

Definition at line 2608 of file inf_fe_lagrange_eval.C.

2608 { return lagrange_eval(v, o, i); }
template<>
Real libMesh::InfFE< 3, LAGRANGE, CARTESIAN >::eval ( Real  v,
Order  o,
unsigned  i 
)
protected

Definition at line 2609 of file inf_fe_lagrange_eval.C.

2609 { return lagrange_eval(v, o, i); }
template<>
Real libMesh::InfFE< 1, INFINITE_MAP, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 75 of file inf_fe_map_eval.C.

75 { return infinite_map_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 2, INFINITE_MAP, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 76 of file inf_fe_map_eval.C.

76 { return infinite_map_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 3, INFINITE_MAP, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 77 of file inf_fe_map_eval.C.

77 { return infinite_map_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 1, LEGENDRE, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 316 of file inf_fe_legendre_eval.C.

316 { return legendre_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 2, LEGENDRE, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 317 of file inf_fe_legendre_eval.C.

317 { return legendre_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 3, LEGENDRE, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 318 of file inf_fe_legendre_eval.C.

318 { return legendre_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 1, JACOBI_30_00, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 463 of file inf_fe_jacobi_30_00_eval.C.

463 { return jacobi_30_00_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 2, JACOBI_30_00, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 464 of file inf_fe_jacobi_30_00_eval.C.

464 { return jacobi_30_00_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 3, JACOBI_30_00, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 465 of file inf_fe_jacobi_30_00_eval.C.

465 { return jacobi_30_00_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 1, JACOBI_20_00, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 468 of file inf_fe_jacobi_20_00_eval.C.

468 { return jacobi_20_00_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 2, JACOBI_20_00, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 469 of file inf_fe_jacobi_20_00_eval.C.

469 { return jacobi_20_00_eval_deriv(v, i); }
template<>
Real libMesh::InfFE< 3, JACOBI_20_00, CARTESIAN >::eval_deriv ( Real  v,
Order  ,
unsigned  i 
)
protected

Definition at line 470 of file inf_fe_jacobi_20_00_eval.C.

470 { return jacobi_20_00_eval_deriv(v, i); }
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
static Real libMesh::InfFE< Dim, T_radial, T_map >::eval_deriv ( Real  v,
Order  o_radial,
unsigned int  i 
)
staticprotected
Returns
The value of the first derivative of the $ i^{th} $ polynomial at coordinate v. See eval for details.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::init_radial_shape_functions(), libMesh::InfFE< Dim, T_radial, T_map >::inverse_map(), and libMesh::InfFE< Dim, T_radial, T_map >::n_quadrature_points().

template<>
Real libMesh::InfFE< 1, LAGRANGE, CARTESIAN >::eval_deriv ( Real  v,
Order  o,
unsigned  i 
)
protected

Definition at line 2613 of file inf_fe_lagrange_eval.C.

2613 { return lagrange_eval_deriv(v, o, i); }
template<>
Real libMesh::InfFE< 2, LAGRANGE, CARTESIAN >::eval_deriv ( Real  v,
Order  o,
unsigned  i 
)
protected

Definition at line 2614 of file inf_fe_lagrange_eval.C.

2614 { return lagrange_eval_deriv(v, o, i); }
template<>
Real libMesh::InfFE< 3, LAGRANGE, CARTESIAN >::eval_deriv ( Real  v,
Order  o,
unsigned  i 
)
protected

Definition at line 2615 of file inf_fe_lagrange_eval.C.

2615 { return lagrange_eval_deriv(v, o, i); }
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
virtual FEContinuity libMesh::InfFE< Dim, T_radial, T_map >::get_continuity ( ) const
inlinevirtual
Returns
The continuity of the element.

Implements libMesh::FEAbstract.

Definition at line 332 of file inf_fe.h.

References libMesh::C_ZERO.

333  { return C_ZERO; } // FIXME - is this true??
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_curl_phi ( ) const
inlineinherited
Returns
The curl of the shape function at the quadrature points.

Definition at line 224 of file fe_base.h.

Referenced by libMesh::ExactSolution::_compute_error().

225  { libmesh_assert(!calculations_started || calculate_curl_phi);
226  calculate_curl_phi = calculate_dphiref = true; return curl_phi; }
std::vector< std::vector< OutputShape > > curl_phi
Definition: fe_base.h:509
const std::vector<Real>& libMesh::FEAbstract::get_curvatures ( ) const
inlineinherited
Returns
The curvatures for use in face integration.

Definition at line 383 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map, libMesh::FEAbstract::attach_quadrature_rule(), libMesh::FEAbstract::n_quadrature_points(), and libMesh::FEAbstract::n_shape_functions().

384  { return this->_fe_map->get_curvatures();}
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:517
template<typename OutputType>
const std::vector<std::vector<OutputTensor> >& libMesh::FEGenericBase< OutputType >::get_d2phi ( ) const
inlineinherited
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phideta2 ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 370 of file fe_base.h.

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

371  { libmesh_assert(!calculations_started || calculate_d2phi);
372  calculate_d2phi = calculate_dphiref = true; return d2phideta2; }
std::vector< std::vector< OutputShape > > d2phideta2
Definition: fe_base.h:572
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidetadzeta ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 378 of file fe_base.h.

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

379  { libmesh_assert(!calculations_started || calculate_d2phi);
std::vector< std::vector< OutputShape > > d2phidetadzeta
Definition: fe_base.h:577
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidx2 ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points.

Definition at line 298 of file fe_base.h.

299  { libmesh_assert(!calculations_started || calculate_d2phi);
300  calculate_d2phi = calculate_dphiref = true; return d2phidx2; }
std::vector< std::vector< OutputShape > > d2phidx2
Definition: fe_base.h:587
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidxdy ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points.

Definition at line 306 of file fe_base.h.

307  { libmesh_assert(!calculations_started || calculate_d2phi);
308  calculate_d2phi = calculate_dphiref = true; return d2phidxdy; }
std::vector< std::vector< OutputShape > > d2phidxdy
Definition: fe_base.h:592
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidxdz ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points.

Definition at line 314 of file fe_base.h.

315  { libmesh_assert(!calculations_started || calculate_d2phi);
316  calculate_d2phi = calculate_dphiref = true; return d2phidxdz; }
std::vector< std::vector< OutputShape > > d2phidxdz
Definition: fe_base.h:597
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidxi2 ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 346 of file fe_base.h.

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

347  { libmesh_assert(!calculations_started || calculate_d2phi);
348  calculate_d2phi = calculate_dphiref = true; return d2phidxi2; }
std::vector< std::vector< OutputShape > > d2phidxi2
Definition: fe_base.h:557
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidxideta ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 354 of file fe_base.h.

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

355  { libmesh_assert(!calculations_started || calculate_d2phi);
std::vector< std::vector< OutputShape > > d2phidxideta
Definition: fe_base.h:562
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidxidzeta ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 362 of file fe_base.h.

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

363  { libmesh_assert(!calculations_started || calculate_d2phi);
std::vector< std::vector< OutputShape > > d2phidxidzeta
Definition: fe_base.h:567
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidy2 ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points.

Definition at line 322 of file fe_base.h.

323  { libmesh_assert(!calculations_started || calculate_d2phi);
324  calculate_d2phi = calculate_dphiref = true; return d2phidy2; }
std::vector< std::vector< OutputShape > > d2phidy2
Definition: fe_base.h:602
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidydz ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points.

Definition at line 330 of file fe_base.h.

331  { libmesh_assert(!calculations_started || calculate_d2phi);
332  calculate_d2phi = calculate_dphiref = true; return d2phidydz; }
std::vector< std::vector< OutputShape > > d2phidydz
Definition: fe_base.h:607
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidz2 ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points.

Definition at line 338 of file fe_base.h.

339  { libmesh_assert(!calculations_started || calculate_d2phi);
340  calculate_d2phi = calculate_dphiref = true; return d2phidz2; }
std::vector< std::vector< OutputShape > > d2phidz2
Definition: fe_base.h:612
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidzeta2 ( ) const
inlineinherited
Returns
The shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 386 of file fe_base.h.

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

387  { libmesh_assert(!calculations_started || calculate_d2phi);
388  calculate_d2phi = calculate_dphiref = true; return d2phidzeta2; }
std::vector< std::vector< OutputShape > > d2phidzeta2
Definition: fe_base.h:582
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdeta2 ( ) const
inlineinherited
Returns
The second partial derivatives in eta.

Definition at line 270 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

271  { return this->_fe_map->get_d2xyzdeta2(); }
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:517
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdetadzeta ( ) const
inlineinherited
Returns
The second partial derivatives in eta-zeta.

Definition at line 300 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

301  { return this->_fe_map->get_d2xyzdetadzeta(); }
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:517
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdxi2 ( ) const
inlineinherited
Returns
The second partial derivatives in xi.

Definition at line 264 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

265  { return this->_fe_map->get_d2xyzdxi2(); }
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:517
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdxideta ( ) const
inlineinherited
Returns
The second partial derivatives in xi-eta.

Definition at line 286 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

287  { return this->_fe_map->get_d2xyzdxideta(); }
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:517
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdxidzeta ( ) const
inlineinherited
Returns
The second partial derivatives in xi-zeta.

Definition at line 294 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

295  { return this->_fe_map->get_d2xyzdxidzeta(); }
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:517
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdzeta2 ( ) const
inlineinherited
Returns
The second partial derivatives in zeta.

Definition at line 278 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

279  { return this->_fe_map->get_d2xyzdzeta2(); }
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:517
const std::vector<