libMesh
Classes | Public Types | Public Member Functions | Static Public Member Functions | Protected Types | Protected Member Functions | Static Protected Member Functions | Protected Attributes | Static Protected Attributes | Private Member Functions | Static Private Attributes | Friends | List of all members
libMesh::InfFE< Dim, T_radial, T_map > Class Template Reference

A specific instantiation of the FEBase class. More...

#include <fe.h>

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

Classes

class  Base
 This nested class contains most of the static methods related to the base part of an infinite element. More...
 
class  Radial
 Infinite elements are in some sense directional, compared to conventional finite elements. More...
 

Public Types

typedef OutputType OutputShape
 Convenient typedefs for gradients of output, hessians of output, and potentially-complex-valued versions of same. More...
 
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)
 Constructor and empty destructor. More...
 
 ~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
 This is at the core of this class. More...
 
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
 Not implemented yet. More...
 
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
 Not implemented yet. More...
 
virtual void side_map (const Elem *, const Elem *, const unsigned int, const std::vector< Point > &, std::vector< Point > &) libmesh_override
 Computes the reference space quadrature points on the side of an element based on the side quadrature points. More...
 
virtual void attach_quadrature_rule (QBase *q) libmesh_override
 The use of quadrature rules with the InfFE class is somewhat different from the approach of the FE class. More...
 
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
 Prints the value of each shape function at each quadrature point. More...
 
void print_dphi (std::ostream &os) const
 Prints the value of each shape function's derivative at each quadrature point. More...
 
void print_d2phi (std::ostream &os) const
 Prints the value of each shape function's second derivatives at each quadrature point. More...
 
unsigned int get_dim () const
 
const std::vector< Point > & get_xyz () const
 
const std::vector< Real > & get_JxW () const
 
const std::vector< RealGradient > & get_dxyzdxi () const
 
const std::vector< RealGradient > & get_dxyzdeta () const
 
const std::vector< RealGradient > & get_dxyzdzeta () const
 
const std::vector< RealGradient > & get_d2xyzdxi2 () const
 
const std::vector< RealGradient > & get_d2xyzdeta2 () const
 
const std::vector< RealGradient > & get_d2xyzdzeta2 () const
 
const std::vector< RealGradient > & get_d2xyzdxideta () const
 
const std::vector< RealGradient > & get_d2xyzdxidzeta () const
 
const std::vector< RealGradient > & get_d2xyzdetadzeta () const
 
const std::vector< Real > & get_dxidx () const
 
const std::vector< Real > & get_dxidy () const
 
const std::vector< Real > & get_dxidz () const
 
const std::vector< Real > & get_detadx () const
 
const std::vector< Real > & get_detady () const
 
const std::vector< Real > & get_detadz () const
 
const std::vector< Real > & get_dzetadx () const
 
const std::vector< Real > & get_dzetady () const
 
const std::vector< Real > & get_dzetadz () const
 
const std::vector< std::vector< Point > > & get_tangents () const
 
const std::vector< Point > & get_normals () const
 
const std::vector< Real > & get_curvatures () const
 
ElemType get_type () const
 
unsigned int get_p_level () const
 
FEType get_fe_type () const
 
Order get_order () const
 
void set_fe_order (int new_order)
 Sets the base FE order of the finite element. More...
 
FEFamily get_family () const
 
const FEMapget_fe_map () const
 
void print_JxW (std::ostream &os) const
 Prints the Jacobian times the weight for each quadrature point. More...
 
void print_xyz (std::ostream &os) const
 Prints the spatial location of each quadrature point (on the physical element). More...
 
void print_info (std::ostream &os) const
 Prints all the relevant information about the current element. More...
 

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)
 Generalized version of shape(), takes an Elem *. More...
 
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)
 Usually, this method would build the nodal soln from the element soln. More...
 
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)
 Takes a number points in physical space (in the physical_points vector) and finds their location on the reference element for the input element elem. More...
 
static UniquePtr< FEGenericBasebuild (const unsigned int dim, const FEType &type)
 Builds a specific finite element type. More...
 
static UniquePtr< FEGenericBasebuild_InfFE (const unsigned int dim, const FEType &type)
 Builds a specific infinite element type. More...
 
static void compute_proj_constraints (DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
 Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to variable number var_number, using generic projections. More...
 
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)
 Creates a local projection on coarse_elem, based on the DoF values in global_vector for it's children. More...
 
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)
 Creates a local projection on coarse_elem, based on the DoF values in global_vector for it's children. More...
 
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)
 Computes the constraint matrix contributions (for meshes with periodic boundary conditions) corresponding to variable number var_number, using generic projections. More...
 
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)
 Computes the nodal constraint contributions (for non-conforming adapted meshes), using Lagrange geometry. More...
 
static void compute_periodic_node_constraints (NodeConstraints &constraints, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const Elem *elem)
 Computes the node position constraint equation contributions (for meshes with periodic boundary conditions) More...
 
static void print_info (std::ostream &out=libMesh::out)
 Prints the reference information, by default to libMesh::out. More...
 
static std::string get_info ()
 Gets a string containing the reference information. More...
 
static unsigned int n_objects ()
 Prints the number of outstanding (created, but not yet destroyed) objects. More...
 
static void enable_print_counter_info ()
 Methods to enable/disable the reference counter output from print_info() More...
 
static void disable_print_counter_info ()
 

Protected Types

typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 Data structure to log the information. More...
 

Protected Member Functions

void update_base_elem (const Elem *inf_elem)
 Updates the protected member base_elem to the appropriate base element for the given inf_elem. More...
 
virtual void init_base_shape_functions (const std::vector< Point > &, const Elem *) libmesh_override
 Do not use this derived member in InfFE<Dim,T_radial,T_map>. More...
 
void init_radial_shape_functions (const Elem *inf_elem, const std::vector< Point > *radial_pts=libmesh_nullptr)
 Some of the member data only depend on the radial part of the infinite element. More...
 
void init_shape_functions (const std::vector< Point > &radial_qp, const std::vector< Point > &base_qp, const Elem *inf_elem)
 Initialize all the data fields like weight, mode, phi, dphidxi, dphideta, dphidzeta, etc. More...
 
void init_face_shape_functions (const std::vector< Point > &qp, const Elem *side)
 Not implemented yet. More...
 
void combine_base_radial (const Elem *inf_elem)
 Combines the shape functions, which were formed in init_shape_functions(Elem *), with geometric data. More...
 
virtual void compute_shape_functions (const Elem *, const std::vector< Point > &) libmesh_override
 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. More...
 
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 ()
 Determine which values are to be calculated, for both the FE itself and for the FEMap. More...
 
void increment_constructor_count (const std::string &name)
 Increments the construction counter. More...
 
void increment_destructor_count (const std::string &name)
 Increments the destruction counter. More...
 

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)
 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. More...
 
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)
 Does the same as compute_node_indices(), but stores the maps for the current element type. More...
 
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)
 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. More...
 

Protected Attributes

std::vector< Realdist
 the radial distance of the base nodes from the origin More...
 
std::vector< Realdweightdv
 the additional radial weight $ 1/{r^2} $ in local coordinates, over all quadrature points. More...
 
std::vector< Realsom
 the radial decay $ 1/r $ in local coordinates. More...
 
std::vector< Realdsomdv
 the first local derivative of the radial decay $ 1/r $ in local coordinates. More...
 
std::vector< std::vector< Real > > mode
 the radial approximation shapes in local coordinates Needed when setting up the overall shape functions. More...
 
std::vector< std::vector< Real > > dmodedv
 the first local derivative of the radial approximation shapes. More...
 
std::vector< std::vector< Real > > radial_map
 the radial mapping shapes in local coordinates More...
 
std::vector< std::vector< Real > > dradialdv_map
 the first local derivative of the radial mapping shapes More...
 
std::vector< Realdphasedxi
 the first local derivative (for 3D, the first in the base) of the phase term in local coordinates. More...
 
std::vector< Realdphasedeta
 the second local derivative (for 3D, the second in the base) of the phase term in local coordinates. More...
 
std::vector< Realdphasedzeta
 the third local derivative (for 3D, the derivative in radial direction) of the phase term in local coordinates. More...
 
std::vector< unsigned int_radial_node_index
 The internal structure of the InfFE – tensor product of base element times radial nodes – has to be determined from the node numbering of the current infinite element. More...
 
std::vector< unsigned int_base_node_index
 The internal structure of the InfFE – tensor product of base element times radial nodes – has to be determined from the node numbering of the current element. More...
 
std::vector< unsigned int_radial_shape_index
 The internal structure of the InfFE – tensor product of base element shapes times radial shapes – has to be determined from the dof numbering scheme of the current infinite element. More...
 
std::vector< unsigned int_base_shape_index
 The internal structure of the InfFE – tensor product of base element shapes times radial shapes – has to be determined from the dof numbering scheme of the current infinite element. More...
 
unsigned int _n_total_approx_sf
 The number of total approximation shape functions for the current configuration. More...
 
unsigned int _n_total_qp
 The total number of quadrature points for the current configuration. More...
 
std::vector< Real_total_qrule_weights
 this vector contains the combined integration weights, so that FEAbstract::compute_map() can still be used More...
 
UniquePtr< QBasebase_qrule
 The quadrature rule for the base element associated with the current infinite element. More...
 
UniquePtr< QBaseradial_qrule
 The quadrature rule for the base element associated with the current infinite element. More...
 
UniquePtr< Elembase_elem
 The base element associated with the current infinite element. More...
 
UniquePtr< FEBasebase_fe
 Have a FE<Dim-1,T_base> handy for base approximation. More...
 
FEType current_fe_type
 This FEType stores the characteristics for which the data structures phi, phi_map etc are currently initialized. More...
 
UniquePtr< FETransformationBase< OutputType > > _fe_trans
 Object that handles computing shape function values, gradients, etc in the physical domain. More...
 
std::vector< std::vector< OutputShape > > phi
 Shape function values. More...
 
std::vector< std::vector< OutputGradient > > dphi
 Shape function derivative values. More...
 
std::vector< std::vector< OutputShape > > curl_phi
 Shape function curl values. More...
 
std::vector< std::vector< OutputDivergence > > div_phi
 Shape function divergence values. More...
 
std::vector< std::vector< OutputShape > > dphidxi
 Shape function derivatives in the xi direction. More...
 
std::vector< std::vector< OutputShape > > dphideta
 Shape function derivatives in the eta direction. More...
 
std::vector< std::vector< OutputShape > > dphidzeta
 Shape function derivatives in the zeta direction. More...
 
std::vector< std::vector< OutputShape > > dphidx
 Shape function derivatives in the x direction. More...
 
std::vector< std::vector< OutputShape > > dphidy
 Shape function derivatives in the y direction. More...
 
std::vector< std::vector< OutputShape > > dphidz
 Shape function derivatives in the z direction. More...
 
std::vector< std::vector< OutputTensor > > d2phi
 Shape function second derivative values. More...
 
std::vector< std::vector< OutputShape > > d2phidxi2
 Shape function second derivatives in the xi direction. More...
 
std::vector< std::vector< OutputShape > > d2phidxideta
 Shape function second derivatives in the xi-eta direction. More...
 
std::vector< std::vector< OutputShape > > d2phidxidzeta
 Shape function second derivatives in the xi-zeta direction. More...
 
std::vector< std::vector< OutputShape > > d2phideta2
 Shape function second derivatives in the eta direction. More...
 
std::vector< std::vector< OutputShape > > d2phidetadzeta
 Shape function second derivatives in the eta-zeta direction. More...
 
std::vector< std::vector< OutputShape > > d2phidzeta2
 Shape function second derivatives in the zeta direction. More...
 
std::vector< std::vector< OutputShape > > d2phidx2
 Shape function second derivatives in the x direction. More...
 
std::vector< std::vector< OutputShape > > d2phidxdy
 Shape function second derivatives in the x-y direction. More...
 
std::vector< std::vector< OutputShape > > d2phidxdz
 Shape function second derivatives in the x-z direction. More...
 
std::vector< std::vector< OutputShape > > d2phidy2
 Shape function second derivatives in the y direction. More...
 
std::vector< std::vector< OutputShape > > d2phidydz
 Shape function second derivatives in the y-z direction. More...
 
std::vector< std::vector< OutputShape > > d2phidz2
 Shape function second derivatives in the z direction. More...
 
std::vector< OutputGradientdphase
 Used for certain infinite element families: the first derivatives of the phase term in global coordinates, over all quadrature points. More...
 
std::vector< RealGradientdweight
 Used for certain infinite element families: the global derivative of the additional radial weight $ 1/{r^2} $, over all quadrature points. More...
 
std::vector< Realweight
 Used for certain infinite element families: the additional radial weight $ 1/{r^2} $ in local coordinates, over all quadrature points. More...
 
UniquePtr< FEMap_fe_map
 
const unsigned int dim
 The dimensionality of the object. More...
 
bool calculations_started
 Have calculations with this object already been started? Then all get_* functions should already have been called. More...
 
bool calculate_phi
 Should we calculate shape functions? More...
 
bool calculate_dphi
 Should we calculate shape function gradients? More...
 
bool calculate_d2phi
 Should we calculate shape function hessians? More...
 
bool calculate_curl_phi
 Should we calculate shape function curls? More...
 
bool calculate_div_phi
 Should we calculate shape function divergences? More...
 
bool calculate_dphiref
 Should we calculate reference shape function gradients? More...
 
FEType fe_type
 The finite element type for this object. More...
 
ElemType elem_type
 The element type the current data structures are set up for. More...
 
unsigned int _p_level
 The p refinement level the current data structures are set up for. More...
 
QBaseqrule
 A pointer to the quadrature rule employed. More...
 
bool shapes_on_quadrature
 A flag indicating if current data structures correspond to quadrature rule points. More...
 

Static Protected Attributes

static Counts _counts
 Actually holds the data. More...
 
static Threads::atomic< unsigned int_n_objects
 The number of objects. More...
 
static Threads::spin_mutex _mutex
 Mutual exclusion object to enable thread-safe reference counting. More...
 
static bool _enable_print_counter = true
 Flag to control whether reference count information is printed when print_info is called. More...
 

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
 When compute_node_indices_fast() is used, this static variable remembers the element type for which the static variables in compute_node_indices_fast() are currently set. More...
 
static bool _warned_for_nodal_soln = false
 static members that are used to issue warning messages only once. More...
 
static bool _warned_for_shape = false
 

Friends

template<unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
class InfFE
 Make all InfFE<Dim,T_radial,T_map> classes friends of each other, so that the protected eval() may be accessed. More...
 

Detailed Description

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

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 Base class for all the infinite geometric element types.

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
The total number of quadrature points for the current configuration.
Definition: inf_fe.h:739
InfMapType inf_map
The coordinate mapping type of the infinite element.
Definition: fe_type.h:257
FEType current_fe_type
This FEType stores the characteristics for which the data structures phi, phi_map etc are currently i...
Definition: inf_fe.h:781
FEGenericBase< Real > FEBase
unsigned int _n_total_approx_sf
The number of total approximation shape functions for the current configuration.
Definition: inf_fe.h:733
FEFamily radial_family
For InfFE, family contains the radial shape family, while base_family contains the approximation type...
Definition: fe_type.h:249
UniquePtr< FEBase > base_fe
Have a FE<Dim-1,T_base> handy for base approximation.
Definition: inf_fe.h:771
FEType fe_type
The finite element type for this object.
Definition: fe_abstract.h:567
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
libMesh::InfFE< Dim, T_radial, T_map >::~InfFE ( )

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
Explicitly request the order as an int.
Definition: fe_type.h:77
UniquePtr< QBase > radial_qrule
The quadrature rule for the base element associated with the current infinite element.
Definition: inf_fe.h:757
OrderWrapper radial_order
The approximation order in the base of the infinite element.
Definition: fe_type.h:236
UniquePtr< QBase > base_qrule
The quadrature rule for the base element associated with the current infinite element.
Definition: inf_fe.h:751
libmesh_assert(j)
QBase * qrule
A pointer to the quadrature rule employed.
Definition: fe_abstract.h:584
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
UniquePtr< FEBase > base_fe
Have a FE<Dim-1,T_base> handy for base approximation.
Definition: inf_fe.h:771
static UniquePtr< QBase > build(const std::string &name, const unsigned int dim, const Order order=INVALID_ORDER)
Builds a specific quadrature rule based on the name string.
FEType fe_type
The finite element type for this object.
Definition: fe_abstract.h:567
template<typename OutputType>
static UniquePtr<FEGenericBase> libMesh::FEGenericBase< OutputType >::build ( const unsigned int  dim,
const FEType type 
)
staticinherited

Builds a specific finite 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 OutputType of this class is not compatible with the output required for the requested type

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::UniformRefinementEstimator::_estimate_error(), assemble(), LinearElasticity::assemble(), assemble_1D(), assemble_biharmonic(), assemble_cd(), assemble_elasticity(), assemble_ellipticdg(), assemble_helmholtz(), assemble_laplace(), assemble_mass(), assemble_matrices(), assemble_poisson(), assemble_shell(), assemble_stokes(), assemble_wave(), libMesh::FEMContext::cached_fe(), libMesh::System::calculate_norm(), compute_stresses(), LinearElasticity::compute_stresses(), LargeDeformationElasticity::compute_stresses(), libMesh::MeshFunction::discontinuous_gradient(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::MeshFunction::gradient(), libMesh::MeshFunction::hessian(), libMesh::InfFE< Dim, T_radial, T_map >::InfFE(), libMesh::InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), libMesh::RBEIMAssembly::initialize_fe(), integrate_function(), LaplaceYoung::jacobian(), LargeDeformationElasticity::jacobian(), main(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator()(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::BoundaryProjectSolution::operator()(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::InfFE< Dim, T_radial, T_map >::reinit(), LaplaceYoung::residual(), LargeDeformationElasticity::residual(), libMesh::HPCoarsenTest::select_refinement(), FETest< order, family, elem_type >::setUp(), and libMesh::Elem::volume().

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

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

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

Definition at line 385 of file fe_base.C.

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

387 {
388  switch (dim)
389  {
390  // 0D
391  case 0:
392  {
393  switch (fet.family)
394  {
395  case LAGRANGE_VEC:
396  return UniquePtr<FEVectorBase>(new FELagrangeVec<0>(fet));
397 
398  default:
399  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
400  }
401  }
402  case 1:
403  {
404  switch (fet.family)
405  {
406  case LAGRANGE_VEC:
407  return UniquePtr<FEVectorBase>(new FELagrangeVec<1>(fet));
408 
409  default:
410  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
411  }
412  }
413  case 2:
414  {
415  switch (fet.family)
416  {
417  case LAGRANGE_VEC:
418  return UniquePtr<FEVectorBase>(new FELagrangeVec<2>(fet));
419 
420  case NEDELEC_ONE:
421  return UniquePtr<FEVectorBase>(new FENedelecOne<2>(fet));
422 
423  default:
424  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
425  }
426  }
427  case 3:
428  {
429  switch (fet.family)
430  {
431  case LAGRANGE_VEC:
432  return UniquePtr<FEVectorBase>(new FELagrangeVec<3>(fet));
433 
434  case NEDELEC_ONE:
435  return UniquePtr<FEVectorBase>(new FENedelecOne<3>(fet));
436 
437  default:
438  libmesh_error_msg("ERROR: Bad FEType.family= " << fet.family);
439  }
440  }
441 
442  default:
443  libmesh_error_msg("Invalid dimension dim = " << dim);
444  } // switch(dim)
445 
446  libmesh_error_msg("We'll never get here!");
447  return UniquePtr<FEVectorBase>();
448 }
FELagrangeVec objects are used for working with vector-valued finite elements.
Definition: fe.h:899
FEFamily family
The type of finite element.
Definition: fe_type.h:203
FENedelecOne objects are used for working with vector-valued Nedelec finite elements of the first kin...
Definition: fe.h:923
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
const unsigned int dim
The dimensionality of the object.
Definition: fe_abstract.h:523
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

Referenced by assemble_wave(), and libMesh::FEMContext::cached_fe().

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

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

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

Definition at line 668 of file fe_base.C.

670 {
671  // No vector types defined... YET.
672  libmesh_not_implemented();
673  return UniquePtr<FEVectorBase>();
674 }
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 802 of file fe_base.C.

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

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

808 {
809  // Side/edge local DOF indices
810  std::vector<unsigned int> new_side_dofs, old_side_dofs;
811 
812  // FIXME: what about 2D shells in 3D space?
813  unsigned int dim = elem->dim();
814 
815  // Cache n_children(); it's a virtual call but it's const.
816  const unsigned int n_children = elem->n_children();
817 
818  // We use local FE objects for now
819  // FIXME: we should use more, external objects instead for efficiency
820  const FEType & base_fe_type = dof_map.variable_type(var);
822  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
824  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
825 
826  UniquePtr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
827  UniquePtr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
828  UniquePtr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
829  std::vector<Point> coarse_qpoints;
830 
831  // The values of the shape functions at the quadrature
832  // points
833  const std::vector<std::vector<OutputShape>> & phi_values =
834  fe->get_phi();
835  const std::vector<std::vector<OutputShape>> & phi_coarse =
836  fe_coarse->get_phi();
837 
838  // The gradients of the shape functions at the quadrature
839  // points on the child element.
840  const std::vector<std::vector<OutputGradient>> * dphi_values =
842  const std::vector<std::vector<OutputGradient>> * dphi_coarse =
844 
845  const FEContinuity cont = fe->get_continuity();
846 
847  if (cont == C_ONE)
848  {
849  const std::vector<std::vector<OutputGradient>> &
850  ref_dphi_values = fe->get_dphi();
851  dphi_values = &ref_dphi_values;
852  const std::vector<std::vector<OutputGradient>> &
853  ref_dphi_coarse = fe_coarse->get_dphi();
854  dphi_coarse = &ref_dphi_coarse;
855  }
856 
857  // The Jacobian * quadrature weight at the quadrature points
858  const std::vector<Real> & JxW =
859  fe->get_JxW();
860 
861  // The XYZ locations of the quadrature points on the
862  // child element
863  const std::vector<Point> & xyz_values =
864  fe->get_xyz();
865 
866 
867 
868  FEType fe_type = base_fe_type, temp_fe_type;
869  const ElemType elem_type = elem->type();
870  fe_type.order = static_cast<Order>(fe_type.order +
871  elem->max_descendant_p_level());
872 
873  // Number of nodes on parent element
874  const unsigned int n_nodes = elem->n_nodes();
875 
876  // Number of dofs on parent element
877  const unsigned int new_n_dofs =
878  FEInterface::n_dofs(dim, fe_type, elem_type);
879 
880  // Fixed vs. free DoFs on edge/face projections
881  std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
882  std::vector<int> free_dof(new_n_dofs, 0);
883 
886  Ue.resize(new_n_dofs); Ue.zero();
887 
888 
889  // When coarsening, in general, we need a series of
890  // projections to ensure a unique and continuous
891  // solution. We start by interpolating nodes, then
892  // hold those fixed and project edges, then
893  // hold those fixed and project faces, then
894  // hold those fixed and project interiors
895 
896  // Copy node values first
897  {
898  std::vector<dof_id_type> node_dof_indices;
899  if (use_old_dof_indices)
900  dof_map.old_dof_indices (elem, node_dof_indices, var);
901  else
902  dof_map.dof_indices (elem, node_dof_indices, var);
903 
904  unsigned int current_dof = 0;
905  for (unsigned int n=0; n!= n_nodes; ++n)
906  {
907  // FIXME: this should go through the DofMap,
908  // not duplicate dof_indices code badly!
909  const unsigned int my_nc =
910  FEInterface::n_dofs_at_node (dim, fe_type,
911  elem_type, n);
912  if (!elem->is_vertex(n))
913  {
914  current_dof += my_nc;
915  continue;
916  }
917 
918  temp_fe_type = base_fe_type;
919  // We're assuming here that child n shares vertex n,
920  // which is wrong on non-simplices right now
921  // ... but this code isn't necessary except on elements
922  // where p refinement creates more vertex dofs; we have
923  // no such elements yet.
924  /*
925  if (elem->child_ptr(n)->p_level() < elem->p_level())
926  {
927  temp_fe_type.order =
928  static_cast<Order>(temp_fe_type.order +
929  elem->child_ptr(n)->p_level());
930  }
931  */
932  const unsigned int nc =
933  FEInterface::n_dofs_at_node (dim, temp_fe_type,
934  elem_type, n);
935  for (unsigned int i=0; i!= nc; ++i)
936  {
937  Ue(current_dof) =
938  old_vector(node_dof_indices[current_dof]);
939  dof_is_fixed[current_dof] = true;
940  current_dof++;
941  }
942  }
943  }
944 
945  // In 3D, project any edge values next
946  if (dim > 2 && cont != DISCONTINUOUS)
947  for (auto e : elem->edge_index_range())
948  {
949  FEInterface::dofs_on_edge(elem, dim, fe_type,
950  e, new_side_dofs);
951 
952  // Some edge dofs are on nodes and already
953  // fixed, others are free to calculate
954  unsigned int free_dofs = 0;
955  for (std::size_t i=0; i != new_side_dofs.size(); ++i)
956  if (!dof_is_fixed[new_side_dofs[i]])
957  free_dof[free_dofs++] = i;
958  Ke.resize (free_dofs, free_dofs); Ke.zero();
959  Fe.resize (free_dofs); Fe.zero();
960  // The new edge coefficients
961  DenseVector<Number> Uedge(free_dofs);
962 
963  // Add projection terms from each child sharing
964  // this edge
965  for (unsigned int c=0; c != n_children; ++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 (auto s : elem->side_index_range())
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 != n_children; ++c)
1103  {
1104  if (!elem->is_child_on_side(c,s))
1105  continue;
1106  const Elem * child = elem->child_ptr(c);
1107 
1108  std::vector<dof_id_type> child_dof_indices;
1109  if (use_old_dof_indices)
1110  dof_map.old_dof_indices (child,
1111  child_dof_indices, var);
1112  else
1113  dof_map.dof_indices (child,
1114  child_dof_indices, var);
1115  const unsigned int child_n_dofs =
1116  cast_int<unsigned int>
1117  (child_dof_indices.size());
1118 
1119  temp_fe_type = base_fe_type;
1120  temp_fe_type.order =
1121  static_cast<Order>(temp_fe_type.order +
1122  child->p_level());
1123 
1124  FEInterface::dofs_on_side(child, dim,
1125  temp_fe_type, s, old_side_dofs);
1126 
1127  // Initialize both child and parent FE data
1128  // on the child's side
1129  fe->attach_quadrature_rule (qsiderule.get());
1130  fe->reinit (child, s);
1131  const unsigned int n_qp = qsiderule->n_points();
1132 
1133  FEInterface::inverse_map (dim, fe_type, elem,
1134  xyz_values, coarse_qpoints);
1135 
1136  fe_coarse->reinit(elem, &coarse_qpoints);
1137 
1138  // Loop over the quadrature points
1139  for (unsigned int qp=0; qp<n_qp; qp++)
1140  {
1141  // solution value at the quadrature point
1142  OutputNumber fineval = libMesh::zero;
1143  // solution grad at the quadrature point
1144  OutputNumberGradient finegrad;
1145 
1146  // Sum the solution values * the DOF
1147  // values at the quadrature point to
1148  // get the solution value and gradient.
1149  for (unsigned int i=0; i<child_n_dofs;
1150  i++)
1151  {
1152  fineval +=
1153  old_vector(child_dof_indices[i]) *
1154  phi_values[i][qp];
1155  if (cont == C_ONE)
1156  finegrad += (*dphi_values)[i][qp] *
1157  old_vector(child_dof_indices[i]);
1158  }
1159 
1160  // Form side projection matrix
1161  for (std::size_t sidei=0, freei=0; sidei != new_side_dofs.size(); ++sidei)
1162  {
1163  unsigned int i = new_side_dofs[sidei];
1164  // fixed DoFs aren't test functions
1165  if (dof_is_fixed[i])
1166  continue;
1167  for (std::size_t sidej=0, freej=0; sidej != new_side_dofs.size(); ++sidej)
1168  {
1169  unsigned int j =
1170  new_side_dofs[sidej];
1171  if (dof_is_fixed[j])
1172  Fe(freei) -=
1173  TensorTools::inner_product(phi_coarse[i][qp],
1174  phi_coarse[j][qp]) *
1175  JxW[qp] * Ue(j);
1176  else
1177  Ke(freei,freej) +=
1178  TensorTools::inner_product(phi_coarse[i][qp],
1179  phi_coarse[j][qp]) *
1180  JxW[qp];
1181  if (cont == C_ONE)
1182  {
1183  if (dof_is_fixed[j])
1184  Fe(freei) -=
1185  TensorTools::inner_product((*dphi_coarse)[i][qp],
1186  (*dphi_coarse)[j][qp]) *
1187  JxW[qp] * Ue(j);
1188  else
1189  Ke(freei,freej) +=
1190  TensorTools::inner_product((*dphi_coarse)[i][qp],
1191  (*dphi_coarse)[j][qp]) *
1192  JxW[qp];
1193  }
1194  if (!dof_is_fixed[j])
1195  freej++;
1196  }
1197  Fe(freei) += TensorTools::inner_product(fineval, phi_coarse[i][qp]) * JxW[qp];
1198  if (cont == C_ONE)
1199  Fe(freei) +=
1200  TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1201  freei++;
1202  }
1203  }
1204  }
1205  Ke.cholesky_solve(Fe, Uside);
1206 
1207  // Transfer new side solutions to element
1208  for (unsigned int i=0; i != free_dofs; ++i)
1209  {
1210  Number & ui = Ue(new_side_dofs[free_dof[i]]);
1212  std::abs(ui - Uside(i)) < TOLERANCE);
1213  ui = Uside(i);
1214  dof_is_fixed[new_side_dofs[free_dof[i]]] = true;
1215  }
1216  }
1217 
1218  // Project the interior values, finally
1219 
1220  // Some interior dofs are on nodes/edges/sides and
1221  // already fixed, others are free to calculate
1222  unsigned int free_dofs = 0;
1223  for (unsigned int i=0; i != new_n_dofs; ++i)
1224  if (!dof_is_fixed[i])
1225  free_dof[free_dofs++] = i;
1226  Ke.resize (free_dofs, free_dofs); Ke.zero();
1227  Fe.resize (free_dofs); Fe.zero();
1228  // The new interior coefficients
1229  DenseVector<Number> Uint(free_dofs);
1230 
1231  // Add projection terms from each child
1232  for (auto & child : elem->child_ref_range())
1233  {
1234  std::vector<dof_id_type> child_dof_indices;
1235  if (use_old_dof_indices)
1236  dof_map.old_dof_indices (&child,
1237  child_dof_indices, var);
1238  else
1239  dof_map.dof_indices (&child,
1240  child_dof_indices, var);
1241  const unsigned int child_n_dofs =
1242  cast_int<unsigned int>
1243  (child_dof_indices.size());
1244 
1245  // Initialize both child and parent FE data
1246  // on the child's quadrature points
1247  fe->attach_quadrature_rule (qrule.get());
1248  fe->reinit (&child);
1249  const unsigned int n_qp = qrule->n_points();
1250 
1251  FEInterface::inverse_map (dim, fe_type, elem,
1252  xyz_values, coarse_qpoints);
1253 
1254  fe_coarse->reinit(elem, &coarse_qpoints);
1255 
1256  // Loop over the quadrature points
1257  for (unsigned int qp=0; qp<n_qp; qp++)
1258  {
1259  // solution value at the quadrature point
1260  OutputNumber fineval = libMesh::zero;
1261  // solution grad at the quadrature point
1262  OutputNumberGradient finegrad;
1263 
1264  // Sum the solution values * the DOF
1265  // values at the quadrature point to
1266  // get the solution value and gradient.
1267  for (unsigned int i=0; i<child_n_dofs; i++)
1268  {
1269  fineval +=
1270  (old_vector(child_dof_indices[i]) *
1271  phi_values[i][qp]);
1272  if (cont == C_ONE)
1273  finegrad += (*dphi_values)[i][qp] *
1274  old_vector(child_dof_indices[i]);
1275  }
1276 
1277  // Form interior projection matrix
1278  for (unsigned int i=0, freei=0;
1279  i != new_n_dofs; ++i)
1280  {
1281  // fixed DoFs aren't test functions
1282  if (dof_is_fixed[i])
1283  continue;
1284  for (unsigned int j=0, freej=0; j !=
1285  new_n_dofs; ++j)
1286  {
1287  if (dof_is_fixed[j])
1288  Fe(freei) -=
1289  TensorTools::inner_product(phi_coarse[i][qp],
1290  phi_coarse[j][qp]) *
1291  JxW[qp] * Ue(j);
1292  else
1293  Ke(freei,freej) +=
1294  TensorTools::inner_product(phi_coarse[i][qp],
1295  phi_coarse[j][qp]) *
1296  JxW[qp];
1297  if (cont == C_ONE)
1298  {
1299  if (dof_is_fixed[j])
1300  Fe(freei) -=
1301  TensorTools::inner_product((*dphi_coarse)[i][qp],
1302  (*dphi_coarse)[j][qp]) *
1303  JxW[qp] * Ue(j);
1304  else
1305  Ke(freei,freej) +=
1306  TensorTools::inner_product((*dphi_coarse)[i][qp],
1307  (*dphi_coarse)[j][qp]) *
1308  JxW[qp];
1309  }
1310  if (!dof_is_fixed[j])
1311  freej++;
1312  }
1313  Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp], fineval) *
1314  JxW[qp];
1315  if (cont == C_ONE)
1316  Fe(freei) += TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1317  freei++;
1318  }
1319  }
1320  }
1321  Ke.cholesky_solve(Fe, Uint);
1322 
1323  // Transfer new interior solutions to element
1324  for (unsigned int i=0; i != free_dofs; ++i)
1325  {
1326  Number & ui = Ue(free_dof[i]);
1328  std::abs(ui - Uint(i)) < TOLERANCE);
1329  ui = Uint(i);
1330  // We should be fixing all dofs by now; no need to keep track of
1331  // that unless we're debugging
1332 #ifndef NDEBUG
1333  dof_is_fixed[free_dof[i]] = true;
1334 #endif
1335  }
1336 
1337 #ifndef NDEBUG
1338  // Make sure every DoF got reached!
1339  for (unsigned int i=0; i != new_n_dofs; ++i)
1340  libmesh_assert(dof_is_fixed[i]);
1341 #endif
1342 }
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)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
Definition: fe_interface.C:510
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
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
Set every element in the matrix to 0.
Definition: dense_matrix.h:792
unsigned int p_level() const
Definition: elem.h:2422
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
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)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
Definition: fe_interface.C:495
ElemType
Defines an enum for geometric element types.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const class libmesh_nullptr_t libmesh_nullptr
TensorTools::IncrementRank< OutputNumber >::type OutputNumberGradient
Definition: fe_base.h:125
OrderWrapper order
The approximation order of the element.
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
After a mesh is refined and repartitioned it is possible that the _send_list will need to be augmente...
Definition: dof_map.C:2378
const Number zero
.
Definition: libmesh.h:178
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:1699
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
The dimensionality of the object.
Definition: fe_abstract.h:523
const Elem * child_ptr(unsigned int i) const
Definition: elem.h:2445
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:569
QBase * qrule
A pointer to the quadrature rule employed.
Definition: fe_abstract.h:584
TensorTools::MakeNumber< OutputShape >::type OutputNumber
Definition: fe_base.h:124
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
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)
Resize the matrix.
Definition: dense_matrix.h:776
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
FEType fe_type
The finite element type for this object.
Definition: fe_abstract.h:567
ElemType elem_type
The element type the current data structures are set up for.
Definition: fe_abstract.h:573
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
Definition: dense_matrix.C:911
unsigned int n_points() const
Definition: quadrature.h:113
This class forms the foundation from which generic finite elements may be derived.
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
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
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 1348 of file fe_base.C.

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

1353 {
1354  Ue.resize(0);
1355 
1356  for (unsigned int v=0; v != dof_map.n_variables(); ++v)
1357  {
1358  DenseVector<Number> Usub;
1359 
1360  coarsened_dof_values(old_vector, dof_map, elem, Usub,
1361  use_old_dof_indices);
1362 
1363  Ue.append (Usub);
1364  }
1365 }
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)
Creates a local projection on coarse_elem, based on the DoF values in global_vector for it&#39;s children...
Definition: fe_base.C:802
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
Shape function derivatives in the xi direction.
Definition: fe_base.h:519
std::vector< std::vector< OutputShape > > dphidzeta
Shape function derivatives in the zeta direction.
Definition: fe_base.h:529
UniquePtr< QBase > radial_qrule
The quadrature rule for the base element associated with the current infinite element.
Definition: inf_fe.h:757
std::vector< Real > dsomdv
the first local derivative of the radial decay in local coordinates.
Definition: inf_fe.h:633
std::vector< std::vector< Real > > dradialdv_map
the first local derivative of the radial mapping shapes
Definition: inf_fe.h:656
std::vector< std::vector< Real > > radial_map
the radial mapping shapes in local coordinates
Definition: inf_fe.h:650
std::vector< Real > dphasedxi
the first local derivative (for 3D, the first in the base) of the phase term in local coordinates...
Definition: inf_fe.h:663
std::vector< unsigned int > _base_node_index
The internal structure of the InfFE – tensor product of base element times radial nodes – has to be...
Definition: inf_fe.h:702
OrderWrapper radial_order
The approximation order in the base of the infinite element.
Definition: fe_type.h:236
static ElemType get_elem_type(const ElemType type)
std::vector< std::vector< Real > > mode
the radial approximation shapes in local coordinates Needed when setting up the overall shape functio...
Definition: inf_fe.h:639
UniquePtr< QBase > base_qrule
The quadrature rule for the base element associated with the current infinite element.
Definition: inf_fe.h:751
UniquePtr< Elem > base_elem
The base element associated with the current infinite element.
Definition: inf_fe.h:763
libmesh_assert(j)
std::vector< std::vector< OutputShape > > phi
Shape function values.
Definition: fe_base.h:499
std::vector< Real > dphasedzeta
the third local derivative (for 3D, the derivative in radial direction) of the phase term in local co...
Definition: inf_fe.h:677
std::vector< Real > dphasedeta
the second local derivative (for 3D, the second in the base) of the phase term in local coordinates...
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
The internal structure of the InfFE – tensor product of base element times radial nodes – has to be...
Definition: inf_fe.h:692
std::vector< Real > som
the radial decay in local coordinates.
Definition: inf_fe.h:628
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:517
std::vector< unsigned int > _radial_shape_index
The internal structure of the InfFE – tensor product of base element shapes times radial shapes – h...
Definition: inf_fe.h:712
UniquePtr< FEBase > base_fe
Have a FE<Dim-1,T_base> handy for base approximation.
Definition: inf_fe.h:771
std::vector< unsigned int > _base_shape_index
The internal structure of the InfFE – tensor product of base element shapes times radial shapes – h...
Definition: inf_fe.h:722
FEType fe_type
The finite element type for this object.
Definition: fe_abstract.h:567
std::vector< Real > dist
the radial distance of the base nodes from the origin
Definition: inf_fe.h:611
std::vector< std::vector< Real > > dmodedv
the first local derivative of the radial approximation shapes.
Definition: inf_fe.h:645
std::vector< std::vector< OutputShape > > dphideta
Shape function derivatives in the eta direction.
Definition: fe_base.h:524
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void libMesh::InfFE< Dim, T_radial, T_map >::compute_data ( const FEType fe_t,
const Elem inf_elem,
FEComputeData data 
)
static

Generalized version of shape(), takes an Elem *.

The data contains both input and output parameters. For frequency domain simulations, the complex-valued shape is returned. In time domain both the computed shape, and the phase is returned.

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

Definition at line 241 of file inf_fe_static.C.

References libMesh::Elem::build_side_ptr(), libMesh::InfFE< Dim, T_radial, T_map >::Radial::decay(), libMesh::InfFE< Dim, T_radial, T_map >::eval(), libMesh::FEComputeData::frequency, libMesh::INFEDGE2, libMesh::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().

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

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

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

Definition at line 387 of file inf_fe_static.C.

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

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

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

Does the same as compute_node_indices(), but stores the maps for the current element type.

Provided the infinite element type changes seldom, this is probably faster than using compute_node_indices () alone. This is possible since the number of nodes is not likely to change.

Definition at line 594 of file inf_fe_static.C.

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

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

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

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

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

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

1375 {
1376  libmesh_assert(elem);
1377 
1378  const unsigned int Dim = elem->dim();
1379 
1380  // Only constrain elements in 2,3D.
1381  if (Dim == 1)
1382  return;
1383 
1384  // Only constrain active elements with this method
1385  if (!elem->active())
1386  return;
1387 
1388  const FEType & base_fe_type = dof_map.variable_type(variable_number);
1389 
1390  // Construct FE objects for this element and its neighbors.
1392  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1393  const FEContinuity cont = my_fe->get_continuity();
1394 
1395  // We don't need to constrain discontinuous elements
1396  if (cont == DISCONTINUOUS)
1397  return;
1398  libmesh_assert (cont == C_ZERO || cont == C_ONE);
1399 
1401  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1402 
1403  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1404  my_fe->attach_quadrature_rule (&my_qface);
1405  std::vector<Point> neigh_qface;
1406 
1407  const std::vector<Real> & JxW = my_fe->get_JxW();
1408  const std::vector<Point> & q_point = my_fe->get_xyz();
1409  const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1410  const std::vector<std::vector<OutputShape>> & neigh_phi =
1411  neigh_fe->get_phi();
1412  const std::vector<Point> * face_normals = libmesh_nullptr;
1413  const std::vector<std::vector<OutputGradient>> * dphi = libmesh_nullptr;
1414  const std::vector<std::vector<OutputGradient>> * neigh_dphi = libmesh_nullptr;
1415 
1416  std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1417  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1418 
1419  if (cont != C_ZERO)
1420  {
1421  const std::vector<Point> & ref_face_normals =
1422  my_fe->get_normals();
1423  face_normals = &ref_face_normals;
1424  const std::vector<std::vector<OutputGradient>> & ref_dphi =
1425  my_fe->get_dphi();
1426  dphi = &ref_dphi;
1427  const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1428  neigh_fe->get_dphi();
1429  neigh_dphi = &ref_neigh_dphi;
1430  }
1431 
1432  DenseMatrix<Real> Ke;
1433  DenseVector<Real> Fe;
1434  std::vector<DenseVector<Real>> Ue;
1435 
1436  // Look at the element faces. Check to see if we need to
1437  // build constraints.
1438  for (auto s : elem->side_index_range())
1439  if (elem->neighbor_ptr(s) != libmesh_nullptr)
1440  {
1441  // Get pointers to the element's neighbor.
1442  const Elem * neigh = elem->neighbor_ptr(s);
1443 
1444  // h refinement constraints:
1445  // constrain dofs shared between
1446  // this element and ones coarser
1447  // than this element.
1448  if (neigh->level() < elem->level())
1449  {
1450  unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
1451  libmesh_assert_less (s_neigh, neigh->n_neighbors());
1452 
1453  // Find the minimum p level; we build the h constraint
1454  // matrix with this and then constrain away all higher p
1455  // DoFs.
1456  libmesh_assert(neigh->active());
1457  const unsigned int min_p_level =
1458  std::min(elem->p_level(), neigh->p_level());
1459 
1460  // we may need to make the FE objects reinit with the
1461  // minimum shared p_level
1462  const unsigned int old_elem_level = elem->p_level();
1463  if (elem->p_level() != min_p_level)
1464  my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() - old_elem_level + min_p_level);
1465  const unsigned int old_neigh_level = neigh->p_level();
1466  if (old_neigh_level != min_p_level)
1467  neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() - old_neigh_level + min_p_level);
1468 
1469  my_fe->reinit(elem, s);
1470 
1471  // This function gets called element-by-element, so there
1472  // will be a lot of memory allocation going on. We can
1473  // at least minimize this for the case of the dof indices
1474  // by efficiently preallocating the requisite storage.
1475  // n_nodes is not necessarily n_dofs, but it is better
1476  // than nothing!
1477  my_dof_indices.reserve (elem->n_nodes());
1478  neigh_dof_indices.reserve (neigh->n_nodes());
1479 
1480  dof_map.dof_indices (elem, my_dof_indices,
1481  variable_number,
1482  min_p_level);
1483  dof_map.dof_indices (neigh, neigh_dof_indices,
1484  variable_number,
1485  min_p_level);
1486 
1487  const unsigned int n_qp = my_qface.n_points();
1488 
1489  FEInterface::inverse_map (Dim, base_fe_type, neigh,
1490  q_point, neigh_qface);
1491 
1492  neigh_fe->reinit(neigh, &neigh_qface);
1493 
1494  // We're only concerned with DOFs whose values (and/or first
1495  // derivatives for C1 elements) are supported on side nodes
1496  FEType elem_fe_type = base_fe_type;
1497  if (old_elem_level != min_p_level)
1498  elem_fe_type.order = base_fe_type.order.get_order() - old_elem_level + min_p_level;
1499  FEType neigh_fe_type = base_fe_type;
1500  if (old_neigh_level != min_p_level)
1501  neigh_fe_type.order = base_fe_type.order.get_order() - old_neigh_level + min_p_level;
1502  FEInterface::dofs_on_side(elem, Dim, elem_fe_type, s, my_side_dofs);
1503  FEInterface::dofs_on_side(neigh, Dim, neigh_fe_type, s_neigh, neigh_side_dofs);
1504 
1505  const unsigned int n_side_dofs =
1506  cast_int<unsigned int>(my_side_dofs.size());
1507  libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
1508 
1509  Ke.resize (n_side_dofs, n_side_dofs);
1510  Ue.resize(n_side_dofs);
1511 
1512  // Form the projection matrix, (inner product of fine basis
1513  // functions against fine test functions)
1514  for (unsigned int is = 0; is != n_side_dofs; ++is)
1515  {
1516  const unsigned int i = my_side_dofs[is];
1517  for (unsigned int js = 0; js != n_side_dofs; ++js)
1518  {
1519  const unsigned int j = my_side_dofs[js];
1520  for (unsigned int qp = 0; qp != n_qp; ++qp)
1521  {
1522  Ke(is,js) += JxW[qp] * TensorTools::inner_product(phi[i][qp], phi[j][qp]);
1523  if (cont != C_ZERO)
1524  Ke(is,js) += JxW[qp] *
1525  TensorTools::inner_product((*dphi)[i][qp] *
1526  (*face_normals)[qp],
1527  (*dphi)[j][qp] *
1528  (*face_normals)[qp]);
1529  }
1530  }
1531  }
1532 
1533  // Form the right hand sides, (inner product of coarse basis
1534  // functions against fine test functions)
1535  for (unsigned int is = 0; is != n_side_dofs; ++is)
1536  {
1537  const unsigned int i = neigh_side_dofs[is];
1538  Fe.resize (n_side_dofs);
1539  for (unsigned int js = 0; js != n_side_dofs; ++js)
1540  {
1541  const unsigned int j = my_side_dofs[js];
1542  for (unsigned int qp = 0; qp != n_qp; ++qp)
1543  {
1544  Fe(js) += JxW[qp] *
1545  TensorTools::inner_product(neigh_phi[i][qp],
1546  phi[j][qp]);
1547  if (cont != C_ZERO)
1548  Fe(js) += JxW[qp] *
1549  TensorTools::inner_product((*neigh_dphi)[i][qp] *
1550  (*face_normals)[qp],
1551  (*dphi)[j][qp] *
1552  (*face_normals)[qp]);
1553  }
1554  }
1555  Ke.cholesky_solve(Fe, Ue[is]);
1556  }
1557 
1558  for (unsigned int js = 0; js != n_side_dofs; ++js)
1559  {
1560  const unsigned int j = my_side_dofs[js];
1561  const dof_id_type my_dof_g = my_dof_indices[j];
1562  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
1563 
1564  // Hunt for "constraining against myself" cases before
1565  // we bother creating a constraint row
1566  bool self_constraint = false;
1567  for (unsigned int is = 0; is != n_side_dofs; ++is)
1568  {
1569  const unsigned int i = neigh_side_dofs[is];
1570  const dof_id_type their_dof_g = neigh_dof_indices[i];
1571  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1572 
1573  if (their_dof_g == my_dof_g)
1574  {
1575 #ifndef NDEBUG
1576  const Real their_dof_value = Ue[is](js);
1577  libmesh_assert_less (std::abs(their_dof_value-1.),
1578  10*TOLERANCE);
1579 
1580  for (unsigned int k = 0; k != n_side_dofs; ++k)
1581  libmesh_assert(k == is ||
1582  std::abs(Ue[k](js)) <
1583  10*TOLERANCE);
1584 #endif
1585 
1586  self_constraint = true;
1587  break;
1588  }
1589  }
1590 
1591  if (self_constraint)
1592  continue;
1593 
1594  DofConstraintRow * constraint_row;
1595 
1596  // we may be running constraint methods concurrently
1597  // on multiple threads, so we need a lock to
1598  // ensure that this constraint is "ours"
1599  {
1600  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1601 
1602  if (dof_map.is_constrained_dof(my_dof_g))
1603  continue;
1604 
1605  constraint_row = &(constraints[my_dof_g]);
1606  libmesh_assert(constraint_row->empty());
1607  }
1608 
1609  for (unsigned int is = 0; is != n_side_dofs; ++is)
1610  {
1611  const unsigned int i = neigh_side_dofs[is];
1612  const dof_id_type their_dof_g = neigh_dof_indices[i];
1613  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1614  libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
1615 
1616  const Real their_dof_value = Ue[is](js);
1617 
1618  if (std::abs(their_dof_value) < 10*TOLERANCE)
1619  continue;
1620 
1621  constraint_row->insert(std::make_pair(their_dof_g,
1622  their_dof_value));
1623  }
1624  }
1625 
1626  my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + old_elem_level - min_p_level);
1627  neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + old_neigh_level - min_p_level);
1628  }
1629 
1630  // p refinement constraints:
1631  // constrain dofs shared between
1632  // active elements and neighbors with
1633  // lower polynomial degrees
1634  const unsigned int min_p_level =
1635  neigh->min_p_level_by_neighbor(elem, elem->p_level());
1636  if (min_p_level < elem->p_level())
1637  {
1638  // Adaptive p refinement of non-hierarchic bases will
1639  // require more coding
1640  libmesh_assert(my_fe->is_hierarchic());
1641  dof_map.constrain_p_dofs(variable_number, elem,
1642  s, min_p_level);
1643  }
1644  }
1645 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
int get_order() const
Explicitly request the order as an int.
Definition: fe_type.h:77
double abs(double a)
bool active() const
Definition: elem.h:2257
unsigned int p_level() const
Definition: elem.h:2422
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
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)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
Definition: fe_interface.C:495
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2083
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
unsigned int which_neighbor_am_i(const Elem *e) const
This function tells you which neighbor e is.
Definition: elem.h:2181
unsigned int min_p_level_by_neighbor(const Elem *neighbor, unsigned int current_min) const
Definition: elem.C:2017
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
unsigned int n_neighbors() const
Definition: elem.h:613
const class libmesh_nullptr_t libmesh_nullptr
OrderWrapper order
The approximation order of the element.
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:1967
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:1734
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:29
std::vector< std::vector< OutputShape > > phi
Shape function values.
Definition: fe_base.h:499
Order default_quadrature_order() const
Definition: fe_type.h:332
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:324
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:569
std::vector< std::vector< OutputGradient > > dphi
Shape function derivative values.
Definition: fe_base.h:504
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
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
A row of the Dof constraint matrix.
Definition: dof_map.h:88
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
virtual unsigned int dim() const =0
unsigned int level() const
Definition: elem.h:2388
This class implements specific orders of Gauss quadrature.
Defines a dense vector for use in Finite Element-type computations.
long double min(long double a, double b)
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
Definition: dense_matrix.C:911
This class forms the foundation from which generic finite elements may be derived.
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)
Constrains degrees of freedom on side s of element elem which correspond to variable number var and t...
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
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
Shape function derivatives in the xi direction.
Definition: fe_base.h:519
std::vector< std::vector< OutputShape > > dphidzeta
Shape function derivatives in the zeta direction.
Definition: fe_base.h:529
std::vector< Real > dphasedxi
the first local derivative (for 3D, the first in the base) of the phase term in local coordinates...
Definition: inf_fe.h:663
std::vector< Real > dweightdv
the additional radial weight in local coordinates, over all quadrature points.
Definition: inf_fe.h:619
unsigned int _n_total_qp
The total number of quadrature points for the current configuration.
Definition: inf_fe.h:739
std::vector< std::vector< OutputShape > > dphidy
Shape function derivatives in the y direction.
Definition: fe_base.h:539
std::vector< std::vector< OutputShape > > dphidx
Shape function derivatives in the x direction.
Definition: fe_base.h:534
std::vector< std::vector< OutputShape > > phi
Shape function values.
Definition: fe_base.h:499
const unsigned int dim
The dimensionality of the object.
Definition: fe_abstract.h:523
std::vector< OutputGradient > dphase
Used for certain infinite element families: the first derivatives of the phase term in global coordin...
Definition: fe_base.h:630
std::vector< Real > dphasedzeta
the third local derivative (for 3D, the derivative in radial direction) of the phase term in local co...
Definition: inf_fe.h:677
std::vector< Real > dphasedeta
the second local derivative (for 3D, the second in the base) of the phase term in local coordinates...
Definition: inf_fe.h:670
std::vector< std::vector< OutputGradient > > dphi
Shape function derivative values.
Definition: fe_base.h:504
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:517
std::vector< std::vector< OutputShape > > dphidz
Shape function derivatives in the z direction.
Definition: fe_base.h:544
std::vector< std::vector< OutputShape > > dphideta
Shape function derivatives in the eta direction.
Definition: fe_base.h:524
std::vector< RealGradient > dweight
Used for certain infinite element families: the global derivative of the additional radial weight ...
Definition: fe_base.h:637
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void libMesh::InfFE< Dim, T_radial, T_map >::compute_shape_indices ( const FEType fet,
const ElemType  inf_elem_type,
const unsigned int  i,
unsigned int base_shape,
unsigned int radial_shape 
)
staticprotected

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

Definition at line 703 of file inf_fe_static.C.

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

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

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

740 {
741  this->calculations_started = true;
742 
743  // If the user forgot to request anything, we'll be safe and
744  // calculate everything:
745 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
746  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_d2phi
747  && !this->calculate_curl_phi && !this->calculate_div_phi)
748  {
749  this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = this->calculate_dphiref = true;
751  {
752  this->calculate_curl_phi = true;
753  this->calculate_div_phi = true;
754  }
755  }
756 #else
757  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_curl_phi && !this->calculate_div_phi)
758  {
759  this->calculate_phi = this->calculate_dphi = this->calculate_dphiref = true;
761  {
762  this->calculate_curl_phi = true;
763  this->calculate_div_phi = true;
764  }
765  }
766 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
767 
768  // Request whichever terms are necessary from the FEMap
769  if (this->calculate_phi)
770  this->_fe_trans->init_map_phi(*this);
771 
772  if (this->calculate_dphiref)
773  this->_fe_trans->init_map_dphi(*this);
774 
775 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
776  if (this->calculate_d2phi)
777  this->_fe_trans->init_map_d2phi(*this);
778 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
779 }
bool calculate_d2phi
Should we calculate shape function hessians?
Definition: fe_abstract.h:544
bool calculate_curl_phi
Should we calculate shape function curls?
Definition: fe_abstract.h:549
FEFamily family
The type of finite element.
Definition: fe_type.h:203
bool calculations_started
Have calculations with this object already been started? Then all get_* functions should already have...
Definition: fe_abstract.h:529
bool calculate_phi
Should we calculate shape functions?
Definition: fe_abstract.h:534
static FEFieldType field_type(const FEType &fe_type)
bool calculate_div_phi
Should we calculate shape function divergences?
Definition: fe_abstract.h:554
bool calculate_dphiref
Should we calculate reference shape function gradients?
Definition: fe_abstract.h:559
bool calculate_dphi
Should we calculate shape function gradients?
Definition: fe_abstract.h:539
UniquePtr< FETransformationBase< OutputType > > _fe_trans
Object that handles computing shape function values, gradients, etc in the physical domain...
Definition: fe_base.h:494
FEType fe_type
The finite element type for this object.
Definition: fe_abstract.h:567
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 107 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::LibMeshInit::LibMeshInit(), and libMesh::ReferenceCounter::n_objects().

108 {
109  _enable_print_counter = false;
110  return;
111 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...
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 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...
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<>