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<>
UniquePtr< FEGenericBase< Real > > build (const unsigned int dim, const FEType &fet)
 
template<>
UniquePtr< FEGenericBase< RealGradient > > build (const unsigned int dim, const FEType &fet)
 
template<>
UniquePtr< FEGenericBase< Real > > build_InfFE (const unsigned int dim, const FEType &fet)
 
template<>
UniquePtr< 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
 
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 UniquePtr< FEGenericBasebuild (const unsigned int dim, const FEType &type)
 
static UniquePtr< 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
 
UniquePtr< QBasebase_qrule
 
UniquePtr< QBaseradial_qrule
 
UniquePtr< Elembase_elem
 
UniquePtr< FEBasebase_fe
 
FEType current_fe_type
 
UniquePtr< 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
 
UniquePtr< 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.reset(FEBase::build(Dim-1, fet).release());
62 }
unsigned int _n_total_qp
Definition: inf_fe.h:739
InfMapType inf_map
Definition: fe_type.h:257
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:249
UniquePtr< FEBase > base_fe
Definition: inf_fe.h:771
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
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::libmesh_assert(), 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.get());
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.reset(QBase::build(q->type(), qrule_dim-1, base_int_order).release());
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:77
UniquePtr< QBase > radial_qrule
Definition: inf_fe.h:757
OrderWrapper radial_order
Definition: fe_type.h:236
UniquePtr< QBase > base_qrule
Definition: inf_fe.h:751
libmesh_assert(j)
UniquePtr< FEBase > base_fe
Definition: inf_fe.h:771
static UniquePtr< QBase > build(const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
template<>
UniquePtr< FEGenericBase< Real > > libMesh::FEGenericBase< Real >::build ( const unsigned int  dim,
const FEType fet 
)
inherited

Definition at line 186 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.

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

Definition at line 387 of file fe_base.C.

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

389 {
390  switch (dim)
391  {
392  // 0D
393  case 0:
394  {
395  switch (fet.family)
396  {
397  case LAGRANGE_VEC:
398  return UniquePtr<FEVectorBase>(new FELagrangeVec<0>(fet));
399 
400  default:
401  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
402  }
403  }
404  case 1:
405  {
406  switch (fet.family)
407  {
408  case LAGRANGE_VEC:
409  return UniquePtr<FEVectorBase>(new FELagrangeVec<1>(fet));
410 
411  default:
412  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
413  }
414  }
415  case 2:
416  {
417  switch (fet.family)
418  {
419  case LAGRANGE_VEC:
420  return UniquePtr<FEVectorBase>(new FELagrangeVec<2>(fet));
421 
422  case NEDELEC_ONE:
423  return UniquePtr<FEVectorBase>(new FENedelecOne<2>(fet));
424 
425  default:
426  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
427  }
428  }
429  case 3:
430  {
431  switch (fet.family)
432  {
433  case LAGRANGE_VEC:
434  return UniquePtr<FEVectorBase>(new FELagrangeVec<3>(fet));
435 
436  case NEDELEC_ONE:
437  return UniquePtr<FEVectorBase>(new FENedelecOne<3>(fet));
438 
439  default:
440  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
441  }
442  }
443 
444  default:
445  libmesh_error_msg("Invalid dimension dim = " << dim);
446  } // switch(dim)
447 
448  libmesh_error_msg("We'll never get here!");
449  return UniquePtr<FEVectorBase>();
450 }
FEFamily family
Definition: fe_type.h:203
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
const unsigned int dim
Definition: fe_abstract.h:517
template<typename OutputType>
static UniquePtr<FEGenericBase> libMesh::FEGenericBase< OutputType >::build_InfFE ( const unsigned int  dim,
const FEType type 
)
staticinherited

Builds a specific infinite element type. A UniquePtr<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

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

Definition at line 463 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.

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

Definition at line 670 of file fe_base.C.

672 {
673  // No vector types defined... YET.
674  libmesh_not_implemented();
675  return UniquePtr<FEVectorBase>();
676 }
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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 804 of file fe_base.C.

References std::abs(), libMesh::C_ONE, libMesh::Elem::child_ptr(), 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::TensorTools::inner_product(), libMesh::FEInterface::inverse_map(), libMesh::Elem::is_child_on_edge(), libMesh::Elem::is_child_on_side(), libMesh::Elem::is_vertex(), libMesh::libmesh_assert(), libmesh_nullptr, libMesh::Elem::max_descendant_p_level(), libMesh::Elem::n_children(), libMesh::FEInterface::n_dofs(), libMesh::FEInterface::n_dofs_at_node(), libMesh::Elem::n_edges(), n_nodes, libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::DofMap::old_dof_indices(), libMesh::FEType::order, libMesh::Elem::p_level(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), 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().

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

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

1356 {
1357  Ue.resize(0);
1358 
1359  for (unsigned int v=0; v != dof_map.n_variables(); ++v)
1360  {
1361  DenseVector<Number> Usub;
1362 
1363  coarsened_dof_values(old_vector, dof_map, elem, Usub,
1364  use_old_dof_indices);
1365 
1366  Ue.append (Usub);
1367  }
1368 }
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:804
unsigned int n_variables() const
Definition: dof_map.h:477
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::libmesh_assert(), 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
UniquePtr< QBase > radial_qrule
Definition: inf_fe.h:757
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::vector< unsigned int > _base_node_index
Definition: inf_fe.h:702
OrderWrapper radial_order
Definition: fe_type.h:236
static ElemType get_elem_type(const ElemType type)
std::vector< std::vector< Real > > mode
Definition: inf_fe.h:639
UniquePtr< QBase > base_qrule
Definition: inf_fe.h:751
UniquePtr< Elem > base_elem
Definition: inf_fe.h:763
libmesh_assert(j)
std::vector< std::vector< OutputShape > > phi
Definition: fe_base.h:499
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
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
std::vector< unsigned int > _radial_shape_index
Definition: inf_fe.h:712
UniquePtr< FEBase > base_fe
Definition: inf_fe.h:771
std::vector< unsigned int > _base_shape_index
Definition: inf_fe.h:722
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 239 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::libmesh_assert(), 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().

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

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

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

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

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

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

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

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

742 {
743  this->calculations_started = true;
744 
745  // If the user forgot to request anything, we'll be safe and
746  // calculate everything:
747 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
748  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_d2phi
749  && !this->calculate_curl_phi && !this->calculate_div_phi)
750  {
751  this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = this->calculate_dphiref = true;
753  {
754  this->calculate_curl_phi = true;
755  this->calculate_div_phi = true;
756  }
757  }
758 #else
759  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_curl_phi && !this->calculate_div_phi)
760  {
761  this->calculate_phi = this->calculate_dphi = this->calculate_dphiref = true;
763  {
764  this->calculate_curl_phi = true;
765  this->calculate_div_phi = true;
766  }
767  }
768 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
769 
770  // Request whichever terms are necessary from the FEMap
771  if (this->calculate_phi)
772  this->_fe_trans->init_map_phi(*this);
773 
774  if (this->calculate_dphiref)
775  this->_fe_trans->init_map_dphi(*this);
776 
777 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
778  if (this->calculate_d2phi)
779  this->_fe_trans->init_map_d2phi(*this);
780 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
781 }
FEFamily family
Definition: fe_type.h:203
static FEFieldType field_type(const FEType &fe_type)
UniquePtr< FETransformationBase< OutputType > > _fe_trans
Definition: fe_base.h:494
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 75 of file inf_fe_map_eval.C.

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

Definition at line 76 of file inf_fe_map_eval.C.

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

Definition at line 77 of file inf_fe_map_eval.C.

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

Definition at line 316 of file inf_fe_legendre_eval.C.

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

Definition at line 317 of file inf_fe_legendre_eval.C.

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

Definition at line 318 of file inf_fe_legendre_eval.C.

318 { return legendre_eval(v, i); }
template<>
Real libMesh::InfFE< 1, JACOBI_30_00, CARTESIAN >::eval ( 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(v, i); }
template<>
Real libMesh::InfFE< 2, JACOBI_30_00, CARTESIAN >::eval ( 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(v, i); }
template<>
Real libMesh::InfFE< 3, JACOBI_30_00, CARTESIAN >::eval ( 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(v, i); }
template<>
Real libMesh::InfFE< 1, JACOBI_20_00, CARTESIAN >::eval ( 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(v, i); }
template<>
Real libMesh::InfFE< 2, JACOBI_20_00, CARTESIAN >::eval ( 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(v, i); }
template<>
Real libMesh::InfFE< 3, JACOBI_20_00, CARTESIAN >::eval ( 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(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 2613 of file inf_fe_lagrange_eval.C.

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

Definition at line 2614 of file inf_fe_lagrange_eval.C.

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

Definition at line 2615 of file inf_fe_lagrange_eval.C.

2615 { 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 81 of file inf_fe_map_eval.C.

81 { 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 82 of file inf_fe_map_eval.C.

82 { 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 83 of file inf_fe_map_eval.C.

83 { 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 322 of file inf_fe_legendre_eval.C.

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

Definition at line 323 of file inf_fe_legendre_eval.C.

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

Definition at line 324 of file inf_fe_legendre_eval.C.

324 { 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 469 of file inf_fe_jacobi_30_00_eval.C.

469 { 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 470 of file inf_fe_jacobi_30_00_eval.C.

470 { 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 471 of file inf_fe_jacobi_30_00_eval.C.

471 { 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 474 of file inf_fe_jacobi_20_00_eval.C.

474 { 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 475 of file inf_fe_jacobi_20_00_eval.C.

475 { 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 476 of file inf_fe_jacobi_20_00_eval.C.

476 { 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 2619 of file inf_fe_lagrange_eval.C.

2619 { 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 2620 of file inf_fe_lagrange_eval.C.

2620 { 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 2621 of file inf_fe_lagrange_eval.C.

2621 { 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().

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

Definition at line 377 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().

378  { return this->_fe_map->get_curvatures();}
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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().

372  calculate_d2phi = calculate_dphiref = true; return d2phideta2; }
libmesh_assert(j)
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().

libmesh_assert(j)
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.

300  calculate_d2phi = calculate_dphiref = true; return d2phidx2; }
std::vector< std::vector< OutputShape > > d2phidx2
Definition: fe_base.h:587
libmesh_assert(j)
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.

308  calculate_d2phi = calculate_dphiref = true; return d2phidxdy; }
libmesh_assert(j)
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.

316  calculate_d2phi = calculate_dphiref = true; return d2phidxdz; }
std::vector< std::vector< OutputShape > > d2phidxdz
Definition: fe_base.h:597
libmesh_assert(j)
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().

348  calculate_d2phi = calculate_dphiref = true; return d2phidxi2; }
libmesh_assert(j)
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().

std::vector< std::vector< OutputShape > > d2phidxideta
Definition: fe_base.h:562
libmesh_assert(j)
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().

libmesh_assert(j)
std::vector< std::vector< OutputShape > > d2phidxidzeta
Definition: fe_base.h:567