libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map > Class Template Reference

#include <fe.h>

Inheritance diagram for libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >:

Classes

class  Base
 
class  Radial
 

Public Types

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

Public Member Functions

 InfFE (const FEType &fet)
 
 ~InfFE ()
 
virtual FEContinuity get_continuity () const
 
virtual bool is_hierarchic () const
 
virtual void reinit (const Elem *elem, const std::vector< Point > *const pts=NULL, const std::vector< Real > *const weights=NULL)
 
virtual void reinit (const Elem *elem, const unsigned int side, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=NULL, const std::vector< Real > *const weights=NULL)
 
virtual void edge_reinit (const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=NULL, const std::vector< Real > *const weights=NULL)
 
virtual void side_map (const Elem *, const Elem *, const unsigned int, const std::vector< Point > &, std::vector< Point > &)
 
virtual void attach_quadrature_rule (QBase *q)
 
virtual unsigned int n_shape_functions () const
 
virtual unsigned int n_quadrature_points () const
 
template<>
AutoPtr< FEGenericBase< Real > > build (const unsigned int dim, const FEType &fet)
 
template<>
AutoPtr< FEGenericBase
< RealGradient > > 
build (const unsigned int dim, const FEType &fet)
 
template<>
AutoPtr< FEGenericBase< Real > > build_InfFE (const unsigned int dim, const FEType &fet)
 
template<>
AutoPtr< FEGenericBase
< RealGradient > > 
build_InfFE (const unsigned int, const FEType &)
 
const std::vector< std::vector
< OutputShape > > & 
get_phi () const
 
const std::vector< std::vector
< OutputGradient > > & 
get_dphi () const
 
const std::vector< std::vector
< OutputShape > > & 
get_curl_phi () const
 
const std::vector< std::vector
< OutputDivergence > > & 
get_div_phi () const
 
const std::vector< std::vector
< OutputShape > > & 
get_dphidx () const
 
const std::vector< std::vector
< OutputShape > > & 
get_dphidy () const
 
const std::vector< std::vector
< OutputShape > > & 
get_dphidz () const
 
const std::vector< std::vector
< OutputShape > > & 
get_dphidxi () const
 
const std::vector< std::vector
< OutputShape > > & 
get_dphideta () const
 
const std::vector< std::vector
< OutputShape > > & 
get_dphidzeta () const
 
const std::vector< std::vector
< OutputTensor > > & 
get_d2phi () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidx2 () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidxdy () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidxdz () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidy2 () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidydz () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidz2 () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidxi2 () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidxideta () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidxidzeta () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phideta2 () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidetadzeta () const
 
const std::vector< std::vector
< OutputShape > > & 
get_d2phidzeta2 () const
 
const std::vector
< OutputGradient > & 
get_dphase () const
 
const std::vector< Real > & get_Sobolev_weight () const
 
const std::vector< RealGradient > & get_Sobolev_dweight () const
 
void print_phi (std::ostream &os) const
 
void print_dphi (std::ostream &os) const
 
void print_d2phi (std::ostream &os) const
 
const std::vector< Point > & get_xyz () const
 
const std::vector< Real > & get_JxW () const
 
const std::vector< RealGradient > & get_dxyzdxi () const
 
const std::vector< RealGradient > & get_dxyzdeta () const
 
const std::vector< RealGradient > & get_dxyzdzeta () const
 
const std::vector< RealGradient > & get_d2xyzdxi2 () const
 
const std::vector< RealGradient > & get_d2xyzdeta2 () const
 
const std::vector< RealGradient > & get_d2xyzdzeta2 () const
 
const std::vector< RealGradient > & get_d2xyzdxideta () const
 
const std::vector< RealGradient > & get_d2xyzdxidzeta () const
 
const std::vector< RealGradient > & get_d2xyzdetadzeta () const
 
const std::vector< Real > & get_dxidx () const
 
const std::vector< Real > & get_dxidy () const
 
const std::vector< Real > & get_dxidz () const
 
const std::vector< Real > & get_detadx () const
 
const std::vector< Real > & get_detady () const
 
const std::vector< Real > & get_detadz () const
 
const std::vector< Real > & get_dzetadx () const
 
const std::vector< Real > & get_dzetady () const
 
const std::vector< Real > & get_dzetadz () const
 
const std::vector< std::vector
< Point > > & 
get_tangents () const
 
const std::vector< Point > & get_normals () const
 
const std::vector< Real > & get_curvatures () const
 
ElemType get_type () const
 
unsigned int get_p_level () const
 
FEType get_fe_type () const
 
Order get_order () const
 
FEFamily get_family () const
 
const FEMapget_fe_map () const
 
void print_JxW (std::ostream &os) const
 
void print_xyz (std::ostream &os) const
 
void print_info (std::ostream &os) const
 

Static Public Member Functions

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

Protected Types

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

Protected Member Functions

void update_base_elem (const Elem *inf_elem)
 
virtual void init_base_shape_functions (const std::vector< Point > &, const Elem *)
 
void init_radial_shape_functions (const Elem *inf_elem)
 
void init_shape_functions (const Elem *inf_elem)
 
void init_face_shape_functions (const std::vector< Point > &qp, const Elem *side)
 
void combine_base_radial (const Elem *inf_elem)
 
virtual void compute_shape_functions (const Elem *, const std::vector< Point > &)
 
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 increment_constructor_count (const std::string &name)
 
void increment_destructor_count (const std::string &name)
 

Static Protected Member Functions

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

Protected Attributes

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

Static Protected Attributes

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

Private Member Functions

virtual bool shapes_need_reinit () const
 

Static Private Attributes

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

Friends

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

Detailed Description

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

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

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

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

Author
Daniel Dreyer
Date
2003

Definition at line 40 of file fe.h.

Member Typedef Documentation

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

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

Definition at line 113 of file reference_counter.h.

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

Definition at line 151 of file fe_base.h.

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

Definition at line 149 of file fe_base.h.

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

Definition at line 152 of file fe_base.h.

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

Definition at line 155 of file fe_base.h.

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

Definition at line 153 of file fe_base.h.

Definition at line 154 of file fe_base.h.

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

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

Definition at line 148 of file fe_base.h.

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

Definition at line 150 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. 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 40 of file inf_fe.C.

References libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::base_fe, libMesh::FEGenericBase< T >::build(), libMesh::FEAbstract::fe_type, libMesh::FEType::inf_map, libMesh::FEType::radial_family, and libMesh::AutoPtr< Tp >::release().

40  :
41  FEBase (Dim, fet),
42 
44  _n_total_qp (0),
45 
46  base_qrule (NULL),
47  radial_qrule (NULL),
48  base_elem (NULL),
49  base_fe (NULL),
50 
51  // initialize the current_fe_type to all the same
52  // values as \p fet (since the FE families and coordinate
53  // map type should @e not change), but use an invalid order
54  // for the radial part (since this is the only order
55  // that may change!).
56  // the data structures like \p phi etc are not initialized
57  // through the constructor, but throught reinit()
58  current_fe_type ( FEType(fet.order,
59  fet.family,
61  fet.radial_family,
62  fet.inf_map) )
63 
64 {
65  // Sanity checks
66  libmesh_assert_equal_to (T_radial, fe_type.radial_family);
67  libmesh_assert_equal_to (T_map, fe_type.inf_map);
68 
69  // build the base_fe object, handle the AutoPtr
70  if (Dim != 1)
71  {
72  AutoPtr<FEBase> ap_fb(FEBase::build(Dim-1, fet));
73  base_fe = ap_fb.release();
74  }
75 }
template<unsigned int Dim, FEFamily T_radial, InfMapType T_map>
libMesh::InfFE< Dim, T_radial, T_map >::~InfFE ( )

Desctructor. Clean up.

Definition at line 82 of file inf_fe.C.

83 {
84  // delete pointers, if necessary
85  delete base_qrule;
86  base_qrule = NULL;
87 
88  delete radial_qrule;
89  radial_qrule = NULL;
90 
91  delete base_elem;
92  base_elem = NULL;
93 
94  delete base_fe;
95  base_fe = NULL;
96 }

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 103 of file inf_fe.C.

References libMesh::QBase::build(), libMesh::QBase::get_dim(), libMesh::QBase::get_order(), libMesh::libmesh_assert(), libMesh::AutoPtr< Tp >::release(), and libMesh::QBase::type().

104 {
105  libmesh_assert(q);
107 
108  const Order base_int_order = q->get_order();
109  const Order radial_int_order = static_cast<Order>(2 * (static_cast<unsigned int>(fe_type.radial_order) + 1) +2);
110  const unsigned int qrule_dim = q->get_dim();
111 
112  if (Dim != 1)
113  {
114  // build a Dim-1 quadrature rule of the type that we received
115  AutoPtr<QBase> apq( QBase::build(q->type(), qrule_dim-1, base_int_order) );
116  base_qrule = apq.release();
118  }
119 
120  // in radial direction, always use Gauss quadrature
121  radial_qrule = new QGauss(1, radial_int_order);
122 
123  // currently not used. But maybe helpful to store the QBase*
124  // with which we initialized our own quadrature rules
125  qrule = q;
126 }
template<>
AutoPtr< FEGenericBase< Real > > libMesh::FEGenericBase< Real >::build ( const unsigned int  dim,
const FEType fet 
)
inherited

Definition at line 184 of file fe_base.C.

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

186 {
187  // The stupid AutoPtr<FEBase> ap(); return ap;
188  // construct is required to satisfy IBM's xlC
189 
190  switch (dim)
191  {
192  // 0D
193  case 0:
194  {
195  switch (fet.family)
196  {
197  case CLOUGH:
198  {
199  AutoPtr<FEBase> ap(new FE<0,CLOUGH>(fet));
200  return ap;
201  }
202 
203  case HERMITE:
204  {
205  AutoPtr<FEBase> ap(new FE<0,HERMITE>(fet));
206  return ap;
207  }
208 
209  case LAGRANGE:
210  {
211  AutoPtr<FEBase> ap(new FE<0,LAGRANGE>(fet));
212  return ap;
213  }
214 
215  case L2_LAGRANGE:
216  {
217  AutoPtr<FEBase> ap(new FE<0,L2_LAGRANGE>(fet));
218  return ap;
219  }
220 
221  case HIERARCHIC:
222  {
223  AutoPtr<FEBase> ap(new FE<0,HIERARCHIC>(fet));
224  return ap;
225  }
226 
227  case L2_HIERARCHIC:
228  {
230  return ap;
231  }
232 
233  case MONOMIAL:
234  {
235  AutoPtr<FEBase> ap(new FE<0,MONOMIAL>(fet));
236  return ap;
237  }
238 
239 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
240  case SZABAB:
241  {
242  AutoPtr<FEBase> ap(new FE<0,SZABAB>(fet));
243  return ap;
244  }
245 
246  case BERNSTEIN:
247  {
248  AutoPtr<FEBase> ap(new FE<0,BERNSTEIN>(fet));
249  return ap;
250  }
251 #endif
252 
253  case XYZ:
254  {
255  AutoPtr<FEBase> ap(new FEXYZ<0>(fet));
256  return ap;
257  }
258 
259  case SCALAR:
260  {
261  AutoPtr<FEBase> ap(new FEScalar<0>(fet));
262  return ap;
263  }
264 
265  default:
266  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
267  libmesh_error();
268  }
269  }
270  // 1D
271  case 1:
272  {
273  switch (fet.family)
274  {
275  case CLOUGH:
276  {
277  AutoPtr<FEBase> ap(new FE<1,CLOUGH>(fet));
278  return ap;
279  }
280 
281  case HERMITE:
282  {
283  AutoPtr<FEBase> ap(new FE<1,HERMITE>(fet));
284  return ap;
285  }
286 
287  case LAGRANGE:
288  {
289  AutoPtr<FEBase> ap(new FE<1,LAGRANGE>(fet));
290  return ap;
291  }
292 
293  case L2_LAGRANGE:
294  {
295  AutoPtr<FEBase> ap(new FE<1,L2_LAGRANGE>(fet));
296  return ap;
297  }
298 
299  case HIERARCHIC:
300  {
301  AutoPtr<FEBase> ap(new FE<1,HIERARCHIC>(fet));
302  return ap;
303  }
304 
305  case L2_HIERARCHIC:
306  {
308  return ap;
309  }
310 
311  case MONOMIAL:
312  {
313  AutoPtr<FEBase> ap(new FE<1,MONOMIAL>(fet));
314  return ap;
315  }
316 
317 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
318  case SZABAB:
319  {
320  AutoPtr<FEBase> ap(new FE<1,SZABAB>(fet));
321  return ap;
322  }
323 
324  case BERNSTEIN:
325  {
326  AutoPtr<FEBase> ap(new FE<1,BERNSTEIN>(fet));
327  return ap;
328  }
329 #endif
330 
331  case XYZ:
332  {
333  AutoPtr<FEBase> ap(new FEXYZ<1>(fet));
334  return ap;
335  }
336 
337  case SCALAR:
338  {
339  AutoPtr<FEBase> ap(new FEScalar<1>(fet));
340  return ap;
341  }
342 
343  default:
344  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
345  libmesh_error();
346  }
347  }
348 
349 
350  // 2D
351  case 2:
352  {
353  switch (fet.family)
354  {
355  case CLOUGH:
356  {
357  AutoPtr<FEBase> ap(new FE<2,CLOUGH>(fet));
358  return ap;
359  }
360 
361  case HERMITE:
362  {
363  AutoPtr<FEBase> ap(new FE<2,HERMITE>(fet));
364  return ap;
365  }
366 
367  case LAGRANGE:
368  {
369  AutoPtr<FEBase> ap(new FE<2,LAGRANGE>(fet));
370  return ap;
371  }
372 
373  case L2_LAGRANGE:
374  {
375  AutoPtr<FEBase> ap(new FE<2,L2_LAGRANGE>(fet));
376  return ap;
377  }
378 
379  case HIERARCHIC:
380  {
381  AutoPtr<FEBase> ap(new FE<2,HIERARCHIC>(fet));
382  return ap;
383  }
384 
385  case L2_HIERARCHIC:
386  {
388  return ap;
389  }
390 
391  case MONOMIAL:
392  {
393  AutoPtr<FEBase> ap(new FE<2,MONOMIAL>(fet));
394  return ap;
395  }
396 
397 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
398  case SZABAB:
399  {
400  AutoPtr<FEBase> ap(new FE<2,SZABAB>(fet));
401  return ap;
402  }
403 
404  case BERNSTEIN:
405  {
406  AutoPtr<FEBase> ap(new FE<2,BERNSTEIN>(fet));
407  return ap;
408  }
409 #endif
410 
411  case XYZ:
412  {
413  AutoPtr<FEBase> ap(new FEXYZ<2>(fet));
414  return ap;
415  }
416 
417  case SCALAR:
418  {
419  AutoPtr<FEBase> ap(new FEScalar<2>(fet));
420  return ap;
421  }
422  default:
423  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
424  libmesh_error();
425  }
426  }
427 
428 
429  // 3D
430  case 3:
431  {
432  switch (fet.family)
433  {
434  case CLOUGH:
435  {
436  libMesh::out << "ERROR: Clough-Tocher elements currently only support 1D and 2D"
437  << std::endl;
438  libmesh_error();
439  }
440 
441  case HERMITE:
442  {
443  AutoPtr<FEBase> ap(new FE<3,HERMITE>(fet));
444  return ap;
445  }
446 
447  case LAGRANGE:
448  {
449  AutoPtr<FEBase> ap(new FE<3,LAGRANGE>(fet));
450  return ap;
451  }
452 
453  case L2_LAGRANGE:
454  {
455  AutoPtr<FEBase> ap(new FE<3,L2_LAGRANGE>(fet));
456  return ap;
457  }
458 
459  case HIERARCHIC:
460  {
461  AutoPtr<FEBase> ap(new FE<3,HIERARCHIC>(fet));
462  return ap;
463  }
464 
465  case L2_HIERARCHIC:
466  {
468  return ap;
469  }
470 
471  case MONOMIAL:
472  {
473  AutoPtr<FEBase> ap(new FE<3,MONOMIAL>(fet));
474  return ap;
475  }
476 
477 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
478  case SZABAB:
479  {
480  AutoPtr<FEBase> ap(new FE<3,SZABAB>(fet));
481  return ap;
482  }
483 
484  case BERNSTEIN:
485  {
486  AutoPtr<FEBase> ap(new FE<3,BERNSTEIN>(fet));
487  return ap;
488  }
489 #endif
490 
491  case XYZ:
492  {
493  AutoPtr<FEBase> ap(new FEXYZ<3>(fet));
494  return ap;
495  }
496 
497  case SCALAR:
498  {
499  AutoPtr<FEBase> ap(new FEScalar<3>(fet));
500  return ap;
501  }
502 
503  default:
504  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
505  libmesh_error();
506  }
507  }
508 
509  default:
510  libmesh_error();
511  }
512 
513  libmesh_error();
514  AutoPtr<FEBase> ap(NULL);
515  return ap;
516 }
template<>
AutoPtr< FEGenericBase< RealGradient > > libMesh::FEGenericBase< RealGradient >::build ( const unsigned int  dim,
const FEType fet 
)
inherited

Definition at line 522 of file fe_base.C.

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

524 {
525  // The stupid AutoPtr<FEBase> ap(); return ap;
526  // construct is required to satisfy IBM's xlC
527 
528  switch (dim)
529  {
530  // 0D
531  case 0:
532  {
533  switch (fet.family)
534  {
535  case LAGRANGE_VEC:
536  {
538  return ap;
539  }
540  default:
541  {
542  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
543  libmesh_error();
544  }
545  }
546  }
547  case 1:
548  {
549  switch (fet.family)
550  {
551  case LAGRANGE_VEC:
552  {
554  return ap;
555  }
556  default:
557  {
558  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
559  libmesh_error();
560  }
561  }
562  }
563  case 2:
564  {
565  switch (fet.family)
566  {
567  case LAGRANGE_VEC:
568  {
570  return ap;
571  }
572  case NEDELEC_ONE:
573  {
574  AutoPtr<FEVectorBase> ap( new FENedelecOne<2>(fet) );
575  return ap;
576  }
577  default:
578  {
579  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
580  libmesh_error();
581  }
582  }
583  }
584  case 3:
585  {
586  switch (fet.family)
587  {
588  case LAGRANGE_VEC:
589  {
591  return ap;
592  }
593  case NEDELEC_ONE:
594  {
595  AutoPtr<FEVectorBase> ap( new FENedelecOne<3>(fet) );
596  return ap;
597  }
598  default:
599  {
600  libMesh::out << "ERROR: Bad FEType.family= " << fet.family << std::endl;
601  libmesh_error();
602  }
603  }
604  }
605 
606  default:
607  libmesh_error();
608 
609  } // switch(dim)
610 
611  libmesh_error();
612  AutoPtr<FEVectorBase> ap(NULL);
613  return ap;
614 }
template<typename T>
static AutoPtr<FEGenericBase> libMesh::FEGenericBase< T >::build_InfFE ( const unsigned int  dim,
const FEType type 
)
staticinherited

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

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

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

Definition at line 627 of file fe_base.C.

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

629 {
630  // The stupid AutoPtr<FEBase> ap(); return ap;
631  // construct is required to satisfy IBM's xlC
632 
633  switch (dim)
634  {
635 
636  // 1D
637  case 1:
638  {
639  switch (fet.radial_family)
640  {
641  case INFINITE_MAP:
642  {
643  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
644  << " with FEFamily = " << fet.radial_family << std::endl;
645  libmesh_error();
646  }
647 
648  case JACOBI_20_00:
649  {
650  switch (fet.inf_map)
651  {
652  case CARTESIAN:
653  {
655  return ap;
656  }
657  default:
658  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
659  << " with InfMapType = " << fet.inf_map << std::endl;
660  libmesh_error();
661  }
662  }
663 
664  case JACOBI_30_00:
665  {
666  switch (fet.inf_map)
667  {
668  case CARTESIAN:
669  {
671  return ap;
672  }
673  default:
674  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
675  << " with InfMapType = " << fet.inf_map << std::endl;
676  libmesh_error();
677  }
678  }
679 
680  case LEGENDRE:
681  {
682  switch (fet.inf_map)
683  {
684  case CARTESIAN:
685  {
687  return ap;
688  }
689  default:
690  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
691  << " with InfMapType = " << fet.inf_map << std::endl;
692  libmesh_error();
693  }
694  }
695 
696  case LAGRANGE:
697  {
698  switch (fet.inf_map)
699  {
700  case CARTESIAN:
701  {
703  return ap;
704  }
705  default:
706  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
707  << " with InfMapType = " << fet.inf_map << std::endl;
708  libmesh_error();
709  }
710  }
711 
712 
713 
714  default:
715  libMesh::err << "ERROR: Bad FEType.radial_family= " << fet.radial_family << std::endl;
716  libmesh_error();
717  }
718 
719  }
720 
721 
722 
723 
724  // 2D
725  case 2:
726  {
727  switch (fet.radial_family)
728  {
729  case INFINITE_MAP:
730  {
731  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
732  << " with FEFamily = " << fet.radial_family << std::endl;
733  libmesh_error();
734  }
735 
736  case JACOBI_20_00:
737  {
738  switch (fet.inf_map)
739  {
740  case CARTESIAN:
741  {
743  return ap;
744  }
745  default:
746  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
747  << " with InfMapType = " << fet.inf_map << std::endl;
748  libmesh_error();
749  }
750  }
751 
752  case JACOBI_30_00:
753  {
754  switch (fet.inf_map)
755  {
756  case CARTESIAN:
757  {
759  return ap;
760  }
761  default:
762  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
763  << " with InfMapType = " << fet.inf_map << std::endl;
764  libmesh_error();
765  }
766  }
767 
768  case LEGENDRE:
769  {
770  switch (fet.inf_map)
771  {
772  case CARTESIAN:
773  {
775  return ap;
776  }
777  default:
778  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
779  << " with InfMapType = " << fet.inf_map << std::endl;
780  libmesh_error();
781  }
782  }
783 
784  case LAGRANGE:
785  {
786  switch (fet.inf_map)
787  {
788  case CARTESIAN:
789  {
791  return ap;
792  }
793  default:
794  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
795  << " with InfMapType = " << fet.inf_map << std::endl;
796  libmesh_error();
797  }
798  }
799 
800 
801 
802  default:
803  libMesh::err << "ERROR: Bad FEType.radial_family= " << fet.radial_family << std::endl;
804  libmesh_error();
805  }
806 
807  }
808 
809 
810 
811 
812  // 3D
813  case 3:
814  {
815  switch (fet.radial_family)
816  {
817  case INFINITE_MAP:
818  {
819  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
820  << " with FEFamily = " << fet.radial_family << std::endl;
821  libmesh_error();
822  }
823 
824  case JACOBI_20_00:
825  {
826  switch (fet.inf_map)
827  {
828  case CARTESIAN:
829  {
831  return ap;
832  }
833  default:
834  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
835  << " with InfMapType = " << fet.inf_map << std::endl;
836  libmesh_error();
837  }
838  }
839 
840  case JACOBI_30_00:
841  {
842  switch (fet.inf_map)
843  {
844  case CARTESIAN:
845  {
847  return ap;
848  }
849  default:
850  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
851  << " with InfMapType = " << fet.inf_map << std::endl;
852  libmesh_error();
853  }
854  }
855 
856  case LEGENDRE:
857  {
858  switch (fet.inf_map)
859  {
860  case CARTESIAN:
861  {
863  return ap;
864  }
865  default:
866  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
867  << " with InfMapType = " << fet.inf_map << std::endl;
868  libmesh_error();
869  }
870  }
871 
872  case LAGRANGE:
873  {
874  switch (fet.inf_map)
875  {
876  case CARTESIAN:
877  {
879  return ap;
880  }
881  default:
882  libMesh::err << "ERROR: Don't build an infinite element " << std::endl
883  << " with InfMapType = " << fet.inf_map << std::endl;
884  libmesh_error();
885  }
886  }
887 
888 
889 
890  default:
891  libMesh::err << "ERROR: Bad FEType.radial_family= " << fet.radial_family << std::endl;
892  libmesh_error();
893  }
894  }
895 
896  default:
897  libmesh_error();
898  }
899 
900  libmesh_error();
901  AutoPtr<FEBase> ap(NULL);
902  return ap;
903 }
template<>
AutoPtr< FEGenericBase< RealGradient > > libMesh::FEGenericBase< RealGradient >::build_InfFE ( const unsigned int  ,
const FEType  
)
inherited

Definition at line 909 of file fe_base.C.

911 {
912  // No vector types defined... YET.
913  libmesh_error();
914  AutoPtr<FEVectorBase> ap(NULL);
915  return ap;
916 }
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.

Definition at line 1029 of file fe_base.C.

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

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

1035 {
1036  // Side/edge local DOF indices
1037  std::vector<unsigned int> new_side_dofs, old_side_dofs;
1038 
1039  // FIXME: what about 2D shells in 3D space?
1040  unsigned int dim = elem->dim();
1041 
1042  // We use local FE objects for now
1043  // FIXME: we should use more, external objects instead for efficiency
1044  const FEType& base_fe_type = dof_map.variable_type(var);
1046  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1048  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1049 
1050  AutoPtr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
1051  AutoPtr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
1052  AutoPtr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
1053  std::vector<Point> coarse_qpoints;
1054 
1055  // The values of the shape functions at the quadrature
1056  // points
1057  const std::vector<std::vector<OutputShape> >& phi_values =
1058  fe->get_phi();
1059  const std::vector<std::vector<OutputShape> >& phi_coarse =
1060  fe_coarse->get_phi();
1061 
1062  // The gradients of the shape functions at the quadrature
1063  // points on the child element.
1064  const std::vector<std::vector<OutputGradient> > *dphi_values =
1065  NULL;
1066  const std::vector<std::vector<OutputGradient> > *dphi_coarse =
1067  NULL;
1068 
1069  const FEContinuity cont = fe->get_continuity();
1070 
1071  if (cont == C_ONE)
1072  {
1073  const std::vector<std::vector<OutputGradient> >&
1074  ref_dphi_values = fe->get_dphi();
1075  dphi_values = &ref_dphi_values;
1076  const std::vector<std::vector<OutputGradient> >&
1077  ref_dphi_coarse = fe_coarse->get_dphi();
1078  dphi_coarse = &ref_dphi_coarse;
1079  }
1080 
1081  // The Jacobian * quadrature weight at the quadrature points
1082  const std::vector<Real>& JxW =
1083  fe->get_JxW();
1084 
1085  // The XYZ locations of the quadrature points on the
1086  // child element
1087  const std::vector<Point>& xyz_values =
1088  fe->get_xyz();
1089 
1090 
1091 
1092  FEType fe_type = base_fe_type, temp_fe_type;
1093  const ElemType elem_type = elem->type();
1094  fe_type.order = static_cast<Order>(fe_type.order +
1095  elem->max_descendant_p_level());
1096 
1097  // Number of nodes on parent element
1098  const unsigned int n_nodes = elem->n_nodes();
1099 
1100  // Number of dofs on parent element
1101  const unsigned int new_n_dofs =
1102  FEInterface::n_dofs(dim, fe_type, elem_type);
1103 
1104  // Fixed vs. free DoFs on edge/face projections
1105  std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
1106  std::vector<int> free_dof(new_n_dofs, 0);
1107 
1108  DenseMatrix<Real> Ke;
1110  Ue.resize(new_n_dofs); Ue.zero();
1111 
1112 
1113  // When coarsening, in general, we need a series of
1114  // projections to ensure a unique and continuous
1115  // solution. We start by interpolating nodes, then
1116  // hold those fixed and project edges, then
1117  // hold those fixed and project faces, then
1118  // hold those fixed and project interiors
1119 
1120  // Copy node values first
1121  {
1122  std::vector<dof_id_type> node_dof_indices;
1123  if (use_old_dof_indices)
1124  dof_map.old_dof_indices (elem, node_dof_indices, var);
1125  else
1126  dof_map.dof_indices (elem, node_dof_indices, var);
1127 
1128  unsigned int current_dof = 0;
1129  for (unsigned int n=0; n!= n_nodes; ++n)
1130  {
1131  // FIXME: this should go through the DofMap,
1132  // not duplicate dof_indices code badly!
1133  const unsigned int my_nc =
1134  FEInterface::n_dofs_at_node (dim, fe_type,
1135  elem_type, n);
1136  if (!elem->is_vertex(n))
1137  {
1138  current_dof += my_nc;
1139  continue;
1140  }
1141 
1142  temp_fe_type = base_fe_type;
1143  // We're assuming here that child n shares vertex n,
1144  // which is wrong on non-simplices right now
1145  // ... but this code isn't necessary except on elements
1146  // where p refinement creates more vertex dofs; we have
1147  // no such elements yet.
1148 /*
1149  if (elem->child(n)->p_level() < elem->p_level())
1150  {
1151  temp_fe_type.order =
1152  static_cast<Order>(temp_fe_type.order +
1153  elem->child(n)->p_level());
1154  }
1155 */
1156  const unsigned int nc =
1157  FEInterface::n_dofs_at_node (dim, temp_fe_type,
1158  elem_type, n);
1159  for (unsigned int i=0; i!= nc; ++i)
1160  {
1161  Ue(current_dof) =
1162  old_vector(node_dof_indices[current_dof]);
1163  dof_is_fixed[current_dof] = true;
1164  current_dof++;
1165  }
1166  }
1167  }
1168 
1169  // In 3D, project any edge values next
1170  if (dim > 2 && cont != DISCONTINUOUS)
1171  for (unsigned int e=0; e != elem->n_edges(); ++e)
1172  {
1173  FEInterface::dofs_on_edge(elem, dim, fe_type,
1174  e, new_side_dofs);
1175 
1176  // Some edge dofs are on nodes and already
1177  // fixed, others are free to calculate
1178  unsigned int free_dofs = 0;
1179  for (unsigned int i=0; i !=
1180  new_side_dofs.size(); ++i)
1181  if (!dof_is_fixed[new_side_dofs[i]])
1182  free_dof[free_dofs++] = i;
1183  Ke.resize (free_dofs, free_dofs); Ke.zero();
1184  Fe.resize (free_dofs); Fe.zero();
1185  // The new edge coefficients
1186  DenseVector<Number> Uedge(free_dofs);
1187 
1188  // Add projection terms from each child sharing
1189  // this edge
1190  for (unsigned int c=0; c != elem->n_children();
1191  ++c)
1192  {
1193  if (!elem->is_child_on_edge(c,e))
1194  continue;
1195  Elem *child = elem->child(c);
1196 
1197  std::vector<dof_id_type> child_dof_indices;
1198  if (use_old_dof_indices)
1199  dof_map.old_dof_indices (child,
1200  child_dof_indices, var);
1201  else
1202  dof_map.dof_indices (child,
1203  child_dof_indices, var);
1204  const unsigned int child_n_dofs =
1205  libmesh_cast_int<unsigned int>
1206  (child_dof_indices.size());
1207 
1208  temp_fe_type = base_fe_type;
1209  temp_fe_type.order =
1210  static_cast<Order>(temp_fe_type.order +
1211  child->p_level());
1212 
1213  FEInterface::dofs_on_edge(child, dim,
1214  temp_fe_type, e, old_side_dofs);
1215 
1216  // Initialize both child and parent FE data
1217  // on the child's edge
1218  fe->attach_quadrature_rule (qedgerule.get());
1219  fe->edge_reinit (child, e);
1220  const unsigned int n_qp = qedgerule->n_points();
1221 
1222  FEInterface::inverse_map (dim, fe_type, elem,
1223  xyz_values, coarse_qpoints);
1224 
1225  fe_coarse->reinit(elem, &coarse_qpoints);
1226 
1227  // Loop over the quadrature points
1228  for (unsigned int qp=0; qp<n_qp; qp++)
1229  {
1230  // solution value at the quadrature point
1231  OutputNumber fineval = libMesh::zero;
1232  // solution grad at the quadrature point
1233  OutputNumberGradient finegrad;
1234 
1235  // Sum the solution values * the DOF
1236  // values at the quadrature point to
1237  // get the solution value and gradient.
1238  for (unsigned int i=0; i<child_n_dofs;
1239  i++)
1240  {
1241  fineval +=
1242  (old_vector(child_dof_indices[i])*
1243  phi_values[i][qp]);
1244  if (cont == C_ONE)
1245  finegrad += (*dphi_values)[i][qp] *
1246  old_vector(child_dof_indices[i]);
1247  }
1248 
1249  // Form edge projection matrix
1250  for (unsigned int sidei=0, freei=0;
1251  sidei != new_side_dofs.size();
1252  ++sidei)
1253  {
1254  unsigned int i = new_side_dofs[sidei];
1255  // fixed DoFs aren't test functions
1256  if (dof_is_fixed[i])
1257  continue;
1258  for (unsigned int sidej=0, freej=0;
1259  sidej != new_side_dofs.size();
1260  ++sidej)
1261  {
1262  unsigned int j =
1263  new_side_dofs[sidej];
1264  if (dof_is_fixed[j])
1265  Fe(freei) -=
1266  TensorTools::inner_product(phi_coarse[i][qp],
1267  phi_coarse[j][qp]) *
1268  JxW[qp] * Ue(j);
1269  else
1270  Ke(freei,freej) +=
1271  TensorTools::inner_product(phi_coarse[i][qp],
1272  phi_coarse[j][qp]) *
1273  JxW[qp];
1274  if (cont == C_ONE)
1275  {
1276  if (dof_is_fixed[j])
1277  Fe(freei) -=
1278  TensorTools::inner_product((*dphi_coarse)[i][qp],
1279  (*dphi_coarse)[j][qp]) *
1280  JxW[qp] * Ue(j);
1281  else
1282  Ke(freei,freej) +=
1283  TensorTools::inner_product((*dphi_coarse)[i][qp],
1284  (*dphi_coarse)[j][qp]) *
1285  JxW[qp];
1286  }
1287  if (!dof_is_fixed[j])
1288  freej++;
1289  }
1290  Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp],
1291  fineval) * JxW[qp];
1292  if (cont == C_ONE)
1293  Fe(freei) +=
1294  TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1295  freei++;
1296  }
1297  }
1298  }
1299  Ke.cholesky_solve(Fe, Uedge);
1300 
1301  // Transfer new edge solutions to element
1302  for (unsigned int i=0; i != free_dofs; ++i)
1303  {
1304  Number &ui = Ue(new_side_dofs[free_dof[i]]);
1306  std::abs(ui - Uedge(i)) < TOLERANCE);
1307  ui = Uedge(i);
1308  dof_is_fixed[new_side_dofs[free_dof[i]]] =
1309  true;
1310  }
1311  }
1312 
1313  // Project any side values (edges in 2D, faces in 3D)
1314  if (dim > 1 && cont != DISCONTINUOUS)
1315  for (unsigned int s=0; s != elem->n_sides(); ++s)
1316  {
1317  FEInterface::dofs_on_side(elem, dim, fe_type,
1318  s, new_side_dofs);
1319 
1320  // Some side dofs are on nodes/edges and already
1321  // fixed, others are free to calculate
1322  unsigned int free_dofs = 0;
1323  for (unsigned int i=0; i !=
1324  new_side_dofs.size(); ++i)
1325  if (!dof_is_fixed[new_side_dofs[i]])
1326  free_dof[free_dofs++] = i;
1327  Ke.resize (free_dofs, free_dofs); Ke.zero();
1328  Fe.resize (free_dofs); Fe.zero();
1329  // The new side coefficients
1330  DenseVector<Number> Uside(free_dofs);
1331 
1332  // Add projection terms from each child sharing
1333  // this side
1334  for (unsigned int c=0; c != elem->n_children();
1335  ++c)
1336  {
1337  if (!elem->is_child_on_side(c,s))
1338  continue;
1339  Elem *child = elem->child(c);
1340 
1341  std::vector<dof_id_type> child_dof_indices;
1342  if (use_old_dof_indices)
1343  dof_map.old_dof_indices (child,
1344  child_dof_indices, var);
1345  else
1346  dof_map.dof_indices (child,
1347  child_dof_indices, var);
1348  const unsigned int child_n_dofs =
1349  libmesh_cast_int<unsigned int>
1350  (child_dof_indices.size());
1351 
1352  temp_fe_type = base_fe_type;
1353  temp_fe_type.order =
1354  static_cast<Order>(temp_fe_type.order +
1355  child->p_level());
1356 
1357  FEInterface::dofs_on_side(child, dim,
1358  temp_fe_type, s, old_side_dofs);
1359 
1360  // Initialize both child and parent FE data
1361  // on the child's side
1362  fe->attach_quadrature_rule (qsiderule.get());
1363  fe->reinit (child, s);
1364  const unsigned int n_qp = qsiderule->n_points();
1365 
1366  FEInterface::inverse_map (dim, fe_type, elem,
1367  xyz_values, coarse_qpoints);
1368 
1369  fe_coarse->reinit(elem, &coarse_qpoints);
1370 
1371  // Loop over the quadrature points
1372  for (unsigned int qp=0; qp<n_qp; qp++)
1373  {
1374  // solution value at the quadrature point
1375  OutputNumber fineval = libMesh::zero;
1376  // solution grad at the quadrature point
1377  OutputNumberGradient finegrad;
1378 
1379  // Sum the solution values * the DOF
1380  // values at the quadrature point to
1381  // get the solution value and gradient.
1382  for (unsigned int i=0; i<child_n_dofs;
1383  i++)
1384  {
1385  fineval +=
1386  old_vector(child_dof_indices[i]) *
1387  phi_values[i][qp];
1388  if (cont == C_ONE)
1389  finegrad += (*dphi_values)[i][qp] *
1390  old_vector(child_dof_indices[i]);
1391  }
1392 
1393  // Form side projection matrix
1394  for (unsigned int sidei=0, freei=0;
1395  sidei != new_side_dofs.size();
1396  ++sidei)
1397  {
1398  unsigned int i = new_side_dofs[sidei];
1399  // fixed DoFs aren't test functions
1400  if (dof_is_fixed[i])
1401  continue;
1402  for (unsigned int sidej=0, freej=0;
1403  sidej != new_side_dofs.size();
1404  ++sidej)
1405  {
1406  unsigned int j =
1407  new_side_dofs[sidej];
1408  if (dof_is_fixed[j])
1409  Fe(freei) -=
1410  TensorTools::inner_product(phi_coarse[i][qp],
1411  phi_coarse[j][qp]) *
1412  JxW[qp] * Ue(j);
1413  else
1414  Ke(freei,freej) +=
1415  TensorTools::inner_product(phi_coarse[i][qp],
1416  phi_coarse[j][qp]) *
1417  JxW[qp];
1418  if (cont == C_ONE)
1419  {
1420  if (dof_is_fixed[j])
1421  Fe(freei) -=
1422  TensorTools::inner_product((*dphi_coarse)[i][qp],
1423  (*dphi_coarse)[j][qp]) *
1424  JxW[qp] * Ue(j);
1425  else
1426  Ke(freei,freej) +=
1427  TensorTools::inner_product((*dphi_coarse)[i][qp],
1428  (*dphi_coarse)[j][qp]) *
1429  JxW[qp];
1430  }
1431  if (!dof_is_fixed[j])
1432  freej++;
1433  }
1434  Fe(freei) += TensorTools::inner_product(fineval, phi_coarse[i][qp]) * JxW[qp];
1435  if (cont == C_ONE)
1436  Fe(freei) +=
1437  TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1438  freei++;
1439  }
1440  }
1441  }
1442  Ke.cholesky_solve(Fe, Uside);
1443 
1444  // Transfer new side solutions to element
1445  for (unsigned int i=0; i != free_dofs; ++i)
1446  {
1447  Number &ui = Ue(new_side_dofs[free_dof[i]]);
1449  std::abs(ui - Uside(i)) < TOLERANCE);
1450  ui = Uside(i);
1451  dof_is_fixed[new_side_dofs[free_dof[i]]] =
1452  true;
1453  }
1454  }
1455 
1456  // Project the interior values, finally
1457 
1458  // Some interior dofs are on nodes/edges/sides and
1459  // already fixed, others are free to calculate
1460  unsigned int free_dofs = 0;
1461  for (unsigned int i=0; i != new_n_dofs; ++i)
1462  if (!dof_is_fixed[i])
1463  free_dof[free_dofs++] = i;
1464  Ke.resize (free_dofs, free_dofs); Ke.zero();
1465  Fe.resize (free_dofs); Fe.zero();
1466  // The new interior coefficients
1467  DenseVector<Number> Uint(free_dofs);
1468 
1469  // Add projection terms from each child
1470  for (unsigned int c=0; c != elem->n_children(); ++c)
1471  {
1472  Elem *child = elem->child(c);
1473 
1474  std::vector<dof_id_type> child_dof_indices;
1475  if (use_old_dof_indices)
1476  dof_map.old_dof_indices (child,
1477  child_dof_indices, var);
1478  else
1479  dof_map.dof_indices (child,
1480  child_dof_indices, var);
1481  const unsigned int child_n_dofs =
1482  libmesh_cast_int<unsigned int>
1483  (child_dof_indices.size());
1484 
1485  // Initialize both child and parent FE data
1486  // on the child's quadrature points
1487  fe->attach_quadrature_rule (qrule.get());
1488  fe->reinit (child);
1489  const unsigned int n_qp = qrule->n_points();
1490 
1491  FEInterface::inverse_map (dim, fe_type, elem,
1492  xyz_values, coarse_qpoints);
1493 
1494  fe_coarse->reinit(elem, &coarse_qpoints);
1495 
1496  // Loop over the quadrature points
1497  for (unsigned int qp=0; qp<n_qp; qp++)
1498  {
1499  // solution value at the quadrature point
1500  OutputNumber fineval = libMesh::zero;
1501  // solution grad at the quadrature point
1502  OutputNumberGradient finegrad;
1503 
1504  // Sum the solution values * the DOF
1505  // values at the quadrature point to
1506  // get the solution value and gradient.
1507  for (unsigned int i=0; i<child_n_dofs; i++)
1508  {
1509  fineval +=
1510  (old_vector(child_dof_indices[i]) *
1511  phi_values[i][qp]);
1512  if (cont == C_ONE)
1513  finegrad += (*dphi_values)[i][qp] *
1514  old_vector(child_dof_indices[i]);
1515  }
1516 
1517  // Form interior projection matrix
1518  for (unsigned int i=0, freei=0;
1519  i != new_n_dofs; ++i)
1520  {
1521  // fixed DoFs aren't test functions
1522  if (dof_is_fixed[i])
1523  continue;
1524  for (unsigned int j=0, freej=0; j !=
1525  new_n_dofs; ++j)
1526  {
1527  if (dof_is_fixed[j])
1528  Fe(freei) -=
1529  TensorTools::inner_product(phi_coarse[i][qp],
1530  phi_coarse[j][qp]) *
1531  JxW[qp] * Ue(j);
1532  else
1533  Ke(freei,freej) +=
1534  TensorTools::inner_product(phi_coarse[i][qp],
1535  phi_coarse[j][qp]) *
1536  JxW[qp];
1537  if (cont == C_ONE)
1538  {
1539  if (dof_is_fixed[j])
1540  Fe(freei) -=
1541  TensorTools::inner_product((*dphi_coarse)[i][qp],
1542  (*dphi_coarse)[j][qp]) *
1543  JxW[qp] * Ue(j);
1544  else
1545  Ke(freei,freej) +=
1546  TensorTools::inner_product((*dphi_coarse)[i][qp],
1547  (*dphi_coarse)[j][qp]) *
1548  JxW[qp];
1549  }
1550  if (!dof_is_fixed[j])
1551  freej++;
1552  }
1553  Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp], fineval) *
1554  JxW[qp];
1555  if (cont == C_ONE)
1556  Fe(freei) += TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1557  freei++;
1558  }
1559  }
1560  }
1561  Ke.cholesky_solve(Fe, Uint);
1562 
1563  // Transfer new interior solutions to element
1564  for (unsigned int i=0; i != free_dofs; ++i)
1565  {
1566  Number &ui = Ue(free_dof[i]);
1568  std::abs(ui - Uint(i)) < TOLERANCE);
1569  ui = Uint(i);
1570  dof_is_fixed[free_dof[i]] = true;
1571  }
1572 
1573  // Make sure every DoF got reached!
1574  for (unsigned int i=0; i != new_n_dofs; ++i)
1575  libmesh_assert(dof_is_fixed[i]);
1576 }
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. Afterwards, the fields are ready to be used to compute global derivatives, the jacobian etc, see FEAbstract::compute_map().

Start logging the combination of radial and base parts

Start logging the combination of radial and base parts

Definition at line 728 of file inf_fe.C.

References libMesh::libmesh_assert(), libMesh::Elem::origin(), libMesh::START_LOG(), libMesh::STOP_LOG(), and libMesh::Elem::type().

729 {
730  libmesh_assert(inf_elem);
731  // at least check whether the base element type is correct.
732  // otherwise this version of computing dist would give problems
733  libmesh_assert_equal_to (base_elem->type(), Base::get_elem_type(inf_elem->type()));
734 
735 
739  START_LOG("combine_base_radial()", "InfFE");
740 
741 
742  // zero the phase, since it is to be summed up
743  std::fill (dphasedxi.begin(), dphasedxi.end(), 0.);
744  std::fill (dphasedeta.begin(), dphasedeta.end(), 0.);
745  std::fill (dphasedzeta.begin(), dphasedzeta.end(), 0.);
746 
747 
748  const unsigned int n_base_mapping_sf =
749  libmesh_cast_int<unsigned int>(dist.size());
750  const Point origin = inf_elem->origin();
751 
752  // for each new infinite element, compute the radial distances
753  for (unsigned int n=0; n<n_base_mapping_sf; n++)
754  dist[n] = Point(base_elem->point(n) - origin).size();
755 
756 
757  switch (Dim)
758  {
759 
760  //------------------------------------------------------------
761  // 1D
762  case 1:
763  {
764  libmesh_not_implemented();
765  break;
766  }
767 
768 
769 
770  //------------------------------------------------------------
771  // 2D
772  case 2:
773  {
774  libmesh_not_implemented();
775  break;
776  }
777 
778 
779 
780  //------------------------------------------------------------
781  // 3D
782  case 3:
783  {
784  // fast access to the approximation and mapping shapes of base_fe
785  const std::vector<std::vector<Real> >& S = base_fe->phi;
786  const std::vector<std::vector<Real> >& Ss = base_fe->dphidxi;
787  const std::vector<std::vector<Real> >& St = base_fe->dphideta;
788  const std::vector<std::vector<Real> >& S_map = (base_fe->get_fe_map()).get_phi_map();
789  const std::vector<std::vector<Real> >& Ss_map = (base_fe->get_fe_map()).get_dphidxi_map();
790  const std::vector<std::vector<Real> >& St_map = (base_fe->get_fe_map()).get_dphideta_map();
791 
792  const unsigned int n_radial_qp = radial_qrule->n_points();
793  const unsigned int n_base_qp = base_qrule-> n_points();
794 
795  const unsigned int n_total_mapping_sf =
796  libmesh_cast_int<unsigned int>(radial_map.size()) * n_base_mapping_sf;
797 
798  const unsigned int n_total_approx_sf = Radial::n_dofs(fe_type.radial_order) * base_fe->n_shape_functions();
799 
800 
801  // compute the phase term derivatives
802  {
803  unsigned int tp=0;
804  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qp's
805  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qp's
806  {
807  // sum over all base shapes, to get the average distance
808  for (unsigned int i=0; i<n_base_mapping_sf; i++)
809  {
810  dphasedxi[tp] += Ss_map[i][bp] * dist[i] * radial_map [1][rp];
811  dphasedeta[tp] += St_map[i][bp] * dist[i] * radial_map [1][rp];
812  dphasedzeta[tp] += S_map [i][bp] * dist[i] * dradialdv_map[1][rp];
813  }
814 
815  tp++;
816 
817  } // loop radial and base qp's
818 
819  }
820 
821  libmesh_assert_equal_to (phi.size(), n_total_approx_sf);
822  libmesh_assert_equal_to (dphidxi.size(), n_total_approx_sf);
823  libmesh_assert_equal_to (dphideta.size(), n_total_approx_sf);
824  libmesh_assert_equal_to (dphidzeta.size(), n_total_approx_sf);
825 
826  // compute the overall approximation shape functions,
827  // pick the appropriate radial and base shapes through using
828  // _base_shape_index and _radial_shape_index
829  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qp's
830  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qp's
831  for (unsigned int ti=0; ti<n_total_approx_sf; ti++) // over _all_ approx_sf
832  {
833  // let the index vectors take care of selecting the appropriate base/radial shape
834  const unsigned int bi = _base_shape_index [ti];
835  const unsigned int ri = _radial_shape_index[ti];
836  phi [ti][bp+rp*n_base_qp] = S [bi][bp] * mode[ri][rp] * som[rp];
837  dphidxi [ti][bp+rp*n_base_qp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
838  dphideta [ti][bp+rp*n_base_qp] = St[bi][bp] * mode[ri][rp] * som[rp];
839  dphidzeta[ti][bp+rp*n_base_qp] = S [bi][bp]
840  * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
841  }
842 
843  std::vector<std::vector<Real> >& phi_map = this->_fe_map->get_phi_map();
844  std::vector<std::vector<Real> >& dphidxi_map = this->_fe_map->get_dphidxi_map();
845  std::vector<std::vector<Real> >& dphideta_map = this->_fe_map->get_dphideta_map();
846  std::vector<std::vector<Real> >& dphidzeta_map = this->_fe_map->get_dphidzeta_map();
847 
848  libmesh_assert_equal_to (phi_map.size(), n_total_mapping_sf);
849  libmesh_assert_equal_to (dphidxi_map.size(), n_total_mapping_sf);
850  libmesh_assert_equal_to (dphideta_map.size(), n_total_mapping_sf);
851  libmesh_assert_equal_to (dphidzeta_map.size(), n_total_mapping_sf);
852 
853  // compute the overall mapping functions,
854  // pick the appropriate radial and base entries through using
855  // _base_node_index and _radial_node_index
856  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qp's
857  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qp's
858  for (unsigned int ti=0; ti<n_total_mapping_sf; ti++) // over all mapping shapes
859  {
860  // let the index vectors take care of selecting the appropriate base/radial mapping shape
861  const unsigned int bi = _base_node_index [ti];
862  const unsigned int ri = _radial_node_index[ti];
863  phi_map [ti][bp+rp*n_base_qp] = S_map [bi][bp] * radial_map [ri][rp];
864  dphidxi_map [ti][bp+rp*n_base_qp] = Ss_map[bi][bp] * radial_map [ri][rp];
865  dphideta_map [ti][bp+rp*n_base_qp] = St_map[bi][bp] * radial_map [ri][rp];
866  dphidzeta_map[ti][bp+rp*n_base_qp] = S_map [bi][bp] * dradialdv_map[ri][rp];
867  }
868 
869 
870  break;
871  }
872 
873 
874  default:
875  libmesh_error();
876  }
877 
878 
882  STOP_LOG("combine_base_radial()", "InfFE");
883 
884 }
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 that 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. Pretty weird, hm!?

Definition at line 239 of file inf_fe_static.C.

References libMesh::Elem::build_side(), libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::Radial::decay(), libMesh::err, libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::eval(), libMesh::FEComputeData::frequency, libMeshEnums::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().

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

References std::abs(), libMesh::Elem::build_side(), libMesh::Elem::default_order(), libMesh::Elem::dim(), libMesh::FEAbstract::fe_type, libMesh::FEInterface::inverse_map(), libMeshEnums::LAGRANGE, libMesh::Elem::level(), libMesh::libmesh_assert(), libMesh::FEInterface::n_dofs(), libMesh::Elem::n_sides(), libMesh::Elem::neighbor(), libMesh::Elem::parent(), libMesh::Real, libMesh::remote_elem, libMesh::FEInterface::shape(), libMesh::Threads::spin_mtx, and libMesh::Elem::subactive().

922 {
923  libmesh_assert(elem);
924 
925  const unsigned int Dim = elem->dim();
926 
927  // Only constrain elements in 2,3D.
928  if (Dim == 1)
929  return;
930 
931  // Only constrain active and ancestor elements
932  if (elem->subactive())
933  return;
934 
935  // We currently always use LAGRANGE mappings for geometry
936  const FEType fe_type(elem->default_order(), LAGRANGE);
937 
938  std::vector<const Node*> my_nodes, parent_nodes;
939 
940  // Look at the element faces. Check to see if we need to
941  // build constraints.
942  for (unsigned int s=0; s<elem->n_sides(); s++)
943  if (elem->neighbor(s) != NULL &&
944  elem->neighbor(s) != remote_elem)
945  if (elem->neighbor(s)->level() < elem->level()) // constrain dofs shared between
946  { // this element and ones coarser
947  // than this element.
948  // Get pointers to the elements of interest and its parent.
949  const Elem* parent = elem->parent();
950 
951  // This can't happen... Only level-0 elements have NULL
952  // parents, and no level-0 elements can be at a higher
953  // level than their neighbors!
954  libmesh_assert(parent);
955 
956  const AutoPtr<Elem> my_side (elem->build_side(s));
957  const AutoPtr<Elem> parent_side (parent->build_side(s));
958 
959  const unsigned int n_side_nodes = my_side->n_nodes();
960 
961  my_nodes.clear();
962  my_nodes.reserve (n_side_nodes);
963  parent_nodes.clear();
964  parent_nodes.reserve (n_side_nodes);
965 
966  for (unsigned int n=0; n != n_side_nodes; ++n)
967  my_nodes.push_back(my_side->get_node(n));
968 
969  for (unsigned int n=0; n != n_side_nodes; ++n)
970  parent_nodes.push_back(parent_side->get_node(n));
971 
972  for (unsigned int my_side_n=0;
973  my_side_n < n_side_nodes;
974  my_side_n++)
975  {
976  libmesh_assert_less (my_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
977 
978  const Node* my_node = my_nodes[my_side_n];
979 
980  // The support point of the DOF
981  const Point& support_point = *my_node;
982 
983  // Figure out where my node lies on their reference element.
984  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
985  parent_side.get(),
986  support_point);
987 
988  // Compute the parent's side shape function values.
989  for (unsigned int their_side_n=0;
990  their_side_n < n_side_nodes;
991  their_side_n++)
992  {
993  libmesh_assert_less (their_side_n, FEInterface::n_dofs(Dim-1, fe_type, parent_side->type()));
994 
995  const Node* their_node = parent_nodes[their_side_n];
996  libmesh_assert(their_node);
997 
998  const Real their_value = FEInterface::shape(Dim-1,
999  fe_type,
1000  parent_side->type(),
1001  their_side_n,
1002  mapped_point);
1003 
1004  const Real their_mag = std::abs(their_value);
1005 #ifdef DEBUG
1006  // Protect for the case u_i ~= u_j,
1007  // in which case i better equal j.
1008  if (their_mag > 0.999)
1009  {
1010  libmesh_assert_equal_to (my_node, their_node);
1011  libmesh_assert_less (std::abs(their_value - 1.), 0.001);
1012  }
1013  else
1014 #endif
1015  // To make nodal constraints useful for constructing
1016  // sparsity patterns faster, we need to get EVERY
1017  // POSSIBLE constraint coupling identified, even if
1018  // there is no coupling in the isoparametric
1019  // Lagrange case.
1020  if (their_mag < 1.e-5)
1021  {
1022  // since we may be running this method concurretly
1023  // on multiple threads we need to acquire a lock
1024  // before modifying the shared constraint_row object.
1025  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1026 
1027  // A reference to the constraint row.
1028  NodeConstraintRow& constraint_row = constraints[my_node].first;
1029 
1030  constraint_row.insert(std::make_pair (their_node,
1031  0.));
1032  }
1033  // To get nodal coordinate constraints right, only
1034  // add non-zero and non-identity values for Lagrange
1035  // basis functions.
1036  else // (1.e-5 <= their_mag <= .999)
1037  {
1038  // since we may be running this method concurretly
1039  // on multiple threads we need to acquire a lock
1040  // before modifying the shared constraint_row object.
1041  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1042 
1043  // A reference to the constraint row.
1044  NodeConstraintRow& constraint_row = constraints[my_node].first;
1045 
1046  constraint_row.insert(std::make_pair (their_node,
1047  their_value));
1048  }
1049  }
1050  }
1051  }
1052 }
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 392 of file inf_fe_static.C.

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

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

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

617 {
618  libmesh_assert_not_equal_to (inf_elem_type, INVALID_ELEM);
619 
620  static std::vector<unsigned int> _static_base_node_index;
621  static std::vector<unsigned int> _static_radial_node_index;
622 
623  /*
624  * fast counterpart to compute_node_indices(), uses local static buffers
625  * to store the index maps. The class member
626  * \p _compute_node_indices_fast_current_elem_type remembers
627  * the current element type.
628  *
629  * Note that there exist non-static members storing the
630  * same data. However, you never know what element type
631  * is currently used by the \p InfFE object, and what
632  * request is currently directed to the static \p InfFE
633  * members (which use \p compute_node_indices_fast()).
634  * So separate these.
635  *
636  * check whether the work for this elemtype has already
637  * been done. If so, use this index. Otherwise, refresh
638  * the buffer to this element type.
639  */
641  {
642  base_node = _static_base_node_index [outer_node_index];
643  radial_node = _static_radial_node_index[outer_node_index];
644  return;
645  }
646  else
647  {
648  // store the map for _all_ nodes for this element type
650 
651  unsigned int n_nodes = libMesh::invalid_uint;
652 
653  switch (inf_elem_type)
654  {
655  case INFEDGE2:
656  {
657  n_nodes = 2;
658  break;
659  }
660  case INFQUAD4:
661  {
662  n_nodes = 4;
663  break;
664  }
665  case INFQUAD6:
666  {
667  n_nodes = 6;
668  break;
669  }
670  case INFHEX8:
671  {
672  n_nodes = 8;
673  break;
674  }
675  case INFHEX16:
676  {
677  n_nodes = 16;
678  break;
679  }
680  case INFHEX18:
681  {
682  n_nodes = 18;
683  break;
684  }
685  case INFPRISM6:
686  {
687  n_nodes = 6;
688  break;
689  }
690  case INFPRISM12:
691  {
692  n_nodes = 12;
693  break;
694  }
695  default:
696  {
697  libMesh::err << "ERROR: Bad infinite element type=" << inf_elem_type
698  << ", node=" << outer_node_index << std::endl;
699  libmesh_error();
700  break;
701  }
702  }
703 
704 
705  _static_base_node_index.resize (n_nodes);
706  _static_radial_node_index.resize(n_nodes);
707 
708  for (unsigned int n=0; n<n_nodes; n++)
709  compute_node_indices (inf_elem_type,
710  n,
711  _static_base_node_index [outer_node_index],
712  _static_radial_node_index[outer_node_index]);
713 
714  // and return for the specified node
715  base_node = _static_base_node_index [outer_node_index];
716  radial_node = _static_radial_node_index[outer_node_index];
717  return;
718  }
719 }
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 1864 of file fe_base.C.

References std::abs(), libMesh::TypeVector< T >::absolute_fuzzy_equals(), libMesh::Elem::active(), libMesh::PeriodicBoundaries::boundary(), libMesh::MeshBase::boundary_info, libMeshEnums::C_ONE, libMeshEnums::C_ZERO, libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::DofMap::constrain_p_dofs(), libMesh::FEType::default_quadrature_order(), libMesh::Elem::dim(), libMeshEnums::DISCONTINUOUS, libMesh::DofMap::dof_indices(), libMesh::DofObject::dof_number(), libMesh::FEInterface::dofs_on_side(), libMesh::err, libMesh::PeriodicBoundaryBase::get_corresponding_pos(), libMesh::Elem::get_node(), 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(), 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(), libMesh::Elem::p_level(), libMesh::PeriodicBoundaryBase::pairedboundary, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::Threads::spin_mtx, swap(), libMesh::DofMap::sys_number(), libMesh::TOLERANCE, and libMesh::DofMap::variable_type().

Referenced by libMesh::FEInterface::compute_periodic_constraints().

1871 {
1872  // Only bother if we truly have periodic boundaries
1873  if (boundaries.empty())
1874  return;
1875 
1876  libmesh_assert(elem);
1877 
1878  // Only constrain active elements with this method
1879  if (!elem->active())
1880  return;
1881 
1882  const unsigned int Dim = elem->dim();
1883 
1884  // We need sys_number and variable_number for DofObject methods
1885  // later
1886  const unsigned int sys_number = dof_map.sys_number();
1887 
1888  const FEType& base_fe_type = dof_map.variable_type(variable_number);
1889 
1890  // Construct FE objects for this element and its pseudo-neighbors.
1892  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1893  const FEContinuity cont = my_fe->get_continuity();
1894 
1895  // We don't need to constrain discontinuous elements
1896  if (cont == DISCONTINUOUS)
1897  return;
1898  libmesh_assert (cont == C_ZERO || cont == C_ONE);
1899 
1900  // We'll use element size to generate relative tolerances later
1901  const Real primary_hmin = elem->hmin();
1902 
1904  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1905 
1906  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1907  my_fe->attach_quadrature_rule (&my_qface);
1908  std::vector<Point> neigh_qface;
1909 
1910  const std::vector<Real>& JxW = my_fe->get_JxW();
1911  const std::vector<Point>& q_point = my_fe->get_xyz();
1912  const std::vector<std::vector<OutputShape> >& phi = my_fe->get_phi();
1913  const std::vector<std::vector<OutputShape> >& neigh_phi =
1914  neigh_fe->get_phi();
1915  const std::vector<Point> *face_normals = NULL;
1916  const std::vector<std::vector<OutputGradient> > *dphi = NULL;
1917  const std::vector<std::vector<OutputGradient> > *neigh_dphi = NULL;
1918  std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1919  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1920 
1921  if (cont != C_ZERO)
1922  {
1923  const std::vector<Point>& ref_face_normals =
1924  my_fe->get_normals();
1925  face_normals = &ref_face_normals;
1926  const std::vector<std::vector<OutputGradient> >& ref_dphi =
1927  my_fe->get_dphi();
1928  dphi = &ref_dphi;
1929  const std::vector<std::vector<OutputGradient> >& ref_neigh_dphi =
1930  neigh_fe->get_dphi();
1931  neigh_dphi = &ref_neigh_dphi;
1932  }
1933 
1934  DenseMatrix<Real> Ke;
1935  DenseVector<Real> Fe;
1936  std::vector<DenseVector<Real> > Ue;
1937 
1938  // Look at the element faces. Check to see if we need to
1939  // build constraints.
1940  for (unsigned int s=0; s<elem->n_sides(); s++)
1941  {
1942  if (elem->neighbor(s))
1943  continue;
1944 
1945  const std::vector<boundary_id_type>& bc_ids = mesh.boundary_info->boundary_ids (elem, s);
1946  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1947  {
1948  const boundary_id_type boundary_id = *id_it;
1949  const PeriodicBoundaryBase *periodic = boundaries.boundary(boundary_id);
1950  if (periodic && periodic->is_my_variable(variable_number))
1951  {
1952  libmesh_assert(point_locator);
1953 
1954  // Get pointers to the element's neighbor.
1955  const Elem* neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s);
1956 
1957  if (neigh == NULL)
1958  {
1959  libMesh::err << "PeriodicBoundaries point locator object returned NULL!" << std::endl;
1960  libmesh_error();
1961  }
1962 
1963  // periodic (and possibly h refinement) constraints:
1964  // constrain dofs shared between
1965  // this element and ones as coarse
1966  // as or coarser than this element.
1967  if (neigh->level() <= elem->level())
1968  {
1969  unsigned int s_neigh =
1970  mesh.boundary_info->side_with_boundary_id (neigh, periodic->pairedboundary);
1971  libmesh_assert_not_equal_to (s_neigh, libMesh::invalid_uint);
1972 
1973 #ifdef LIBMESH_ENABLE_AMR
1974  // Find the minimum p level; we build the h constraint
1975  // matrix with this and then constrain away all higher p
1976  // DoFs.
1977  libmesh_assert(neigh->active());
1978  const unsigned int min_p_level =
1979  std::min(elem->p_level(), neigh->p_level());
1980 
1981  // we may need to make the FE objects reinit with the
1982  // minimum shared p_level
1983  // FIXME - I hate using const_cast<> and avoiding
1984  // accessor functions; there's got to be a
1985  // better way to do this!
1986  const unsigned int old_elem_level = elem->p_level();
1987  if (old_elem_level != min_p_level)
1988  (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
1989  const unsigned int old_neigh_level = neigh->p_level();
1990  if (old_neigh_level != min_p_level)
1991  (const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
1992 #endif // #ifdef LIBMESH_ENABLE_AMR
1993 
1994  // We can do a projection with a single integration,
1995  // due to the assumption of nested finite element
1996  // subspaces.
1997  // FIXME: it might be more efficient to do nodes,
1998  // then edges, then side, to reduce the size of the
1999  // Cholesky factorization(s)
2000  my_fe->reinit(elem, s);
2001 
2002  dof_map.dof_indices (elem, my_dof_indices,
2003  variable_number);
2004  dof_map.dof_indices (neigh, neigh_dof_indices,
2005  variable_number);
2006 
2007  const unsigned int n_qp = my_qface.n_points();
2008 
2009  // Translate the quadrature points over to the
2010  // neighbor's boundary
2011  std::vector<Point> neigh_point(q_point.size());
2012  for (unsigned int i=0; i != neigh_point.size(); ++i)
2013  neigh_point[i] = periodic->get_corresponding_pos(q_point[i]);
2014 
2015  FEInterface::inverse_map (Dim, base_fe_type, neigh,
2016  neigh_point, neigh_qface);
2017 
2018  neigh_fe->reinit(neigh, &neigh_qface);
2019 
2020  // We're only concerned with DOFs whose values (and/or first
2021  // derivatives for C1 elements) are supported on side nodes
2022  FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
2023  FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
2024 
2025  // We're done with functions that examine Elem::p_level(),
2026  // so let's unhack those levels
2027 #ifdef LIBMESH_ENABLE_AMR
2028  if (elem->p_level() != old_elem_level)
2029  (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
2030  if (neigh->p_level() != old_neigh_level)
2031  (const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
2032 #endif // #ifdef LIBMESH_ENABLE_AMR
2033 
2034  const unsigned int n_side_dofs =
2035  libmesh_cast_int<unsigned int>
2036  (my_side_dofs.size());
2037  libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
2038 
2039  Ke.resize (n_side_dofs, n_side_dofs);
2040  Ue.resize(n_side_dofs);
2041 
2042  // Form the projection matrix, (inner product of fine basis
2043  // functions against fine test functions)
2044  for (unsigned int is = 0; is != n_side_dofs; ++is)
2045  {
2046  const unsigned int i = my_side_dofs[is];
2047  for (unsigned int js = 0; js != n_side_dofs; ++js)
2048  {
2049  const unsigned int j = my_side_dofs[js];
2050  for (unsigned int qp = 0; qp != n_qp; ++qp)
2051  {
2052  Ke(is,js) += JxW[qp] *
2053  TensorTools::inner_product(phi[i][qp],
2054  phi[j][qp]);
2055  if (cont != C_ZERO)
2056  Ke(is,js) += JxW[qp] *
2057  TensorTools::inner_product((*dphi)[i][qp] *
2058  (*face_normals)[qp],
2059  (*dphi)[j][qp] *
2060  (*face_normals)[qp]);
2061  }
2062  }
2063  }
2064 
2065  // Form the right hand sides, (inner product of coarse basis
2066  // functions against fine test functions)
2067  for (unsigned int is = 0; is != n_side_dofs; ++is)
2068  {
2069  const unsigned int i = neigh_side_dofs[is];
2070  Fe.resize (n_side_dofs);
2071  for (unsigned int js = 0; js != n_side_dofs; ++js)
2072  {
2073  const unsigned int j = my_side_dofs[js];
2074  for (unsigned int qp = 0; qp != n_qp; ++qp)
2075  {
2076  Fe(js) += JxW[qp] *
2077  TensorTools::inner_product(neigh_phi[i][qp],
2078  phi[j][qp]);
2079  if (cont != C_ZERO)
2080  Fe(js) += JxW[qp] *
2081  TensorTools::inner_product((*neigh_dphi)[i][qp] *
2082  (*face_normals)[qp],
2083  (*dphi)[j][qp] *
2084  (*face_normals)[qp]);
2085  }
2086  }
2087  Ke.cholesky_solve(Fe, Ue[is]);
2088  }
2089 
2090  // Make sure we're not adding recursive constraints
2091  // due to the redundancy in the way we add periodic
2092  // boundary constraints
2093  //
2094  // In order for this to work while threaded or on
2095  // distributed meshes, we need a rigorous way to
2096  // avoid recursive constraints. Here it is:
2097  //
2098  // For vertex DoFs, if there is a "prior" element
2099  // (i.e. a coarser element or an equally refined
2100  // element with a lower id) on this boundary which
2101  // contains the vertex point, then we will avoid
2102  // generating constraints; the prior element (or
2103  // something prior to it) may do so. If we are the
2104  // most prior (or "primary") element on this
2105  // boundary sharing this point, then we look at the
2106  // boundary periodic to us, we find the primary
2107  // element there, and if that primary is coarser or
2108  // equal-but-lower-id, then our vertex dofs are
2109  // constrained in terms of that element.
2110  //
2111  // For edge DoFs, if there is a coarser element
2112  // on this boundary sharing this edge, then we will
2113  // avoid generating constraints (we will be
2114  // constrained indirectly via AMR constraints
2115  // connecting us to the coarser element's DoFs). If
2116  // we are the coarsest element sharing this edge,
2117  // then we generate constraints if and only if we
2118  // are finer than the coarsest element on the
2119  // boundary periodic to us sharing the corresponding
2120  // periodic edge, or if we are at equal level but
2121  // our edge nodes have higher ids than the periodic
2122  // edge nodes (sorted from highest to lowest, then
2123  // compared lexicographically)
2124  //
2125  // For face DoFs, we generate constraints if we are
2126  // finer than our periodic neighbor, or if we are at
2127  // equal level but our element id is higher than its
2128  // element id.
2129  //
2130  // If the primary neighbor is also the current elem
2131  // (a 1-element-thick mesh) then we choose which
2132  // vertex dofs to constrain via lexicographic
2133  // ordering on point locations
2134 
2135  // FIXME: This code doesn't yet properly handle
2136  // cases where multiple different periodic BCs
2137  // intersect.
2138  std::set<dof_id_type> my_constrained_dofs;
2139 
2140  for (unsigned int n = 0; n != elem->n_nodes(); ++n)
2141  {
2142  if (!elem->is_node_on_side(n,s))
2143  continue;
2144 
2145  const Node* my_node = elem->get_node(n);
2146 
2147  if (elem->is_vertex(n))
2148  {
2149  // Find all boundary ids that include this
2150  // point and have periodic boundary
2151  // conditions for this variable
2152  std::set<boundary_id_type> point_bcids;
2153 
2154  for (unsigned int new_s = 0; new_s !=
2155  elem->n_sides(); ++new_s)
2156  {
2157  if (!elem->is_node_on_side(n,new_s))
2158  continue;
2159 
2160  const std::vector<boundary_id_type>
2161  new_bc_ids = mesh.boundary_info->boundary_ids (elem, s);
2162  for (std::vector<boundary_id_type>::const_iterator
2163  new_id_it=new_bc_ids.begin(); new_id_it!=new_bc_ids.end(); ++new_id_it)
2164  {
2165  const boundary_id_type new_boundary_id = *new_id_it;
2166  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2167  if (new_periodic && new_periodic->is_my_variable(variable_number))
2168  {
2169  point_bcids.insert(new_boundary_id);
2170  }
2171  }
2172  }
2173 
2174  // See if this vertex has point neighbors to
2175  // defer to
2176  if (primary_boundary_point_neighbor
2177  (elem, *my_node, *mesh.boundary_info, point_bcids) != elem)
2178  continue;
2179 
2180  // Find the complementary boundary id set
2181  std::set<boundary_id_type> point_pairedids;
2182  for (std::set<boundary_id_type>::const_iterator i =
2183  point_bcids.begin(); i != point_bcids.end(); ++i)
2184  {
2185  const boundary_id_type new_boundary_id = *i;
2186  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2187  point_pairedids.insert(new_periodic->pairedboundary);
2188  }
2189 
2190  // What do we want to constrain against?
2191  const Elem* primary_elem = NULL;
2192  const Elem* main_neigh = NULL;
2193  Point main_pt = *my_node,
2194  primary_pt = *my_node;
2195 
2196  for (std::set<boundary_id_type>::const_iterator i =
2197  point_bcids.begin(); i != point_bcids.end(); ++i)
2198  {
2199  // Find the corresponding periodic point and
2200  // its primary neighbor
2201  const boundary_id_type new_boundary_id = *i;
2202  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2203 
2204  const Point neigh_pt =
2205  new_periodic->get_corresponding_pos(*my_node);
2206 
2207  // If the point is getting constrained
2208  // to itself by this PBC then we don't
2209  // generate any constraints
2210  if (neigh_pt.absolute_fuzzy_equals
2211  (*my_node, primary_hmin*TOLERANCE))
2212  continue;
2213 
2214  // Otherwise we'll have a constraint in
2215  // one direction or another
2216  if (!primary_elem)
2217  primary_elem = elem;
2218 
2219  const Elem *primary_neigh = primary_boundary_point_neighbor
2220  (neigh, neigh_pt, *mesh.boundary_info,
2221  point_pairedids);
2222 
2223  libmesh_assert(primary_neigh);
2224 
2225  if (new_boundary_id == boundary_id)
2226  {
2227  main_neigh = primary_neigh;
2228  main_pt = neigh_pt;
2229  }
2230 
2231  // Finer elements will get constrained in
2232  // terms of coarser neighbors, not the
2233  // other way around
2234  if ((primary_neigh->level() > primary_elem->level()) ||
2235 
2236  // For equal-level elements, the one with
2237  // higher id gets constrained in terms of
2238  // the one with lower id
2239  (primary_neigh->level() == primary_elem->level() &&
2240  primary_neigh->id() > primary_elem->id()) ||
2241 
2242  // On a one-element-thick mesh, we compare
2243  // points to see what side gets constrained
2244  (primary_neigh == primary_elem &&
2245  (neigh_pt > primary_pt)))
2246  continue;
2247 
2248  primary_elem = primary_neigh;
2249  primary_pt = neigh_pt;
2250  }
2251 
2252  if (!primary_elem ||
2253  primary_elem != main_neigh ||
2254  primary_pt != main_pt)
2255  continue;
2256  }
2257  else if (elem->is_edge(n))
2258  {
2259  // Find which edge we're on
2260  unsigned int e=0;
2261  for (; e != elem->n_edges(); ++e)
2262  {
2263  if (elem->is_node_on_edge(n,e))
2264  break;
2265  }
2266  libmesh_assert_less (e, elem->n_edges());
2267 
2268  // Find the edge end nodes
2269  Node *e1 = NULL,
2270  *e2 = NULL;
2271  for (unsigned int nn = 0; nn != elem->n_nodes(); ++nn)
2272  {
2273  if (nn == n)
2274  continue;
2275 
2276  if (elem->is_node_on_edge(nn, e))
2277  {
2278  if (e1 == NULL)
2279  {
2280  e1 = elem->get_node(nn);
2281  }
2282  else
2283  {
2284  e2 = elem->get_node(nn);
2285  break;
2286  }
2287  }
2288  }
2289  libmesh_assert (e1 && e2);
2290 
2291  // Find all boundary ids that include this
2292  // edge and have periodic boundary
2293  // conditions for this variable
2294  std::set<boundary_id_type> edge_bcids;
2295 
2296  for (unsigned int new_s = 0; new_s !=
2297  elem->n_sides(); ++new_s)
2298  {
2299  if (!elem->is_node_on_side(n,new_s))
2300  continue;
2301 
2302  const std::vector<boundary_id_type>&
2303  new_bc_ids = mesh.boundary_info->boundary_ids (elem, s);
2304  for (std::vector<boundary_id_type>::const_iterator
2305  new_id_it=new_bc_ids.begin(); new_id_it!=new_bc_ids.end(); ++new_id_it)
2306  {
2307  const boundary_id_type new_boundary_id = *new_id_it;
2308  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2309  if (new_periodic && new_periodic->is_my_variable(variable_number))
2310  {
2311  edge_bcids.insert(new_boundary_id);
2312  }
2313  }
2314  }
2315 
2316 
2317  // See if this edge has neighbors to defer to
2318  if (primary_boundary_edge_neighbor
2319  (elem, *e1, *e2, *mesh.boundary_info, edge_bcids) != elem)
2320  continue;
2321 
2322  // Find the complementary boundary id set
2323  std::set<boundary_id_type> edge_pairedids;
2324  for (std::set<boundary_id_type>::const_iterator i =
2325  edge_bcids.begin(); i != edge_bcids.end(); ++i)
2326  {
2327  const boundary_id_type new_boundary_id = *i;
2328  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2329  edge_pairedids.insert(new_periodic->pairedboundary);
2330  }
2331 
2332  // What do we want to constrain against?
2333  const Elem* primary_elem = NULL;
2334  const Elem* main_neigh = NULL;
2335  Point main_pt1 = *e1,
2336  main_pt2 = *e2,
2337  primary_pt1 = *e1,
2338  primary_pt2 = *e2;
2339 
2340  for (std::set<boundary_id_type>::const_iterator i =
2341  edge_bcids.begin(); i != edge_bcids.end(); ++i)
2342  {
2343  // Find the corresponding periodic edge and
2344  // its primary neighbor
2345  const boundary_id_type new_boundary_id = *i;
2346  const PeriodicBoundaryBase *new_periodic = boundaries.boundary(new_boundary_id);
2347 
2348  Point neigh_pt1 = new_periodic->get_corresponding_pos(*e1),
2349  neigh_pt2 = new_periodic->get_corresponding_pos(*e2);
2350 
2351  // If the edge is getting constrained
2352  // to itself by this PBC then we don't
2353  // generate any constraints
2354  if (neigh_pt1.absolute_fuzzy_equals
2355  (*e1, primary_hmin*TOLERANCE) &&
2356  neigh_pt2.absolute_fuzzy_equals
2357  (*e2, primary_hmin*TOLERANCE))
2358  continue;
2359 
2360  // Otherwise we'll have a constraint in
2361  // one direction or another
2362  if (!primary_elem)
2363  primary_elem = elem;
2364 
2365  const Elem *primary_neigh = primary_boundary_edge_neighbor
2366  (neigh, neigh_pt1, neigh_pt2, *mesh.boundary_info,
2367  edge_pairedids);
2368 
2369  libmesh_assert(primary_neigh);
2370 
2371  if (new_boundary_id == boundary_id)
2372  {
2373  main_neigh = primary_neigh;
2374  main_pt1 = neigh_pt1;
2375  main_pt2 = neigh_pt2;
2376  }
2377 
2378  // If we have a one-element thick mesh,
2379  // we'll need to sort our points to get a
2380  // consistent ordering rule
2381  //
2382  // Use >= in this test to make sure that,
2383  // for angular constraints, no node gets
2384  // constrained to itself.
2385  if (primary_neigh == primary_elem)
2386  {
2387  if (primary_pt1 > primary_pt2)
2388  std::swap(primary_pt1, primary_pt2);
2389  if (neigh_pt1 > neigh_pt2)
2390  std::swap(neigh_pt1, neigh_pt2);
2391 
2392  if (neigh_pt2 >= primary_pt2)
2393  continue;
2394  }
2395 
2396  // Otherwise:
2397  // Finer elements will get constrained in
2398  // terms of coarser ones, not the other way
2399  // around
2400  if ((primary_neigh->level() > primary_elem->level()) ||
2401 
2402  // For equal-level elements, the one with
2403  // higher id gets constrained in terms of
2404  // the one with lower id
2405  (primary_neigh->level() == primary_elem->level() &&
2406  primary_neigh->id() > primary_elem->id()))
2407  continue;
2408 
2409  primary_elem = primary_neigh;
2410  primary_pt1 = neigh_pt1;
2411  primary_pt2 = neigh_pt2;
2412  }
2413 
2414  if (!primary_elem ||
2415  primary_elem != main_neigh ||
2416  primary_pt1 != main_pt1 ||
2417  primary_pt2 != main_pt2)
2418  continue;
2419  }
2420  else if (elem->is_face(n))
2421  {
2422  // If we have a one-element thick mesh,
2423  // use the ordering of the face node and its
2424  // periodic counterpart to determine what
2425  // gets constrained
2426  if (neigh == elem)
2427  {
2428  const Point neigh_pt =
2429  periodic->get_corresponding_pos(*my_node);
2430  if (neigh_pt > *my_node)
2431  continue;
2432  }
2433 
2434  // Otherwise:
2435  // Finer elements will get constrained in
2436  // terms of coarser ones, not the other way
2437  // around
2438  if ((neigh->level() > elem->level()) ||
2439 
2440  // For equal-level elements, the one with
2441  // higher id gets constrained in terms of
2442  // the one with lower id
2443  (neigh->level() == elem->level() &&
2444  neigh->id() > elem->id()))
2445  continue;
2446  }
2447 
2448  // If we made it here without hitting a continue
2449  // statement, then we're at a node whose dofs
2450  // should be constrained by this element's
2451  // calculations.
2452  const unsigned int n_comp =
2453  my_node->n_comp(sys_number, variable_number);
2454 
2455  for (unsigned int i=0; i != n_comp; ++i)
2456  my_constrained_dofs.insert
2457  (my_node->dof_number
2458  (sys_number, variable_number, i));
2459  }
2460 
2461  // FIXME: old code for disambiguating periodic BCs:
2462  // this is not threadsafe nor safe to run on a
2463  // non-serialized mesh.
2464  /*
2465  std::vector<bool> recursive_constraint(n_side_dofs, false);
2466 
2467  for (unsigned int is = 0; is != n_side_dofs; ++is)
2468  {
2469  const unsigned int i = neigh_side_dofs[is];
2470  const dof_id_type their_dof_g = neigh_dof_indices[i];
2471  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2472 
2473  {
2474  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2475 
2476  if (!dof_map.is_constrained_dof(their_dof_g))
2477  continue;
2478  }
2479 
2480  DofConstraintRow& their_constraint_row =
2481  constraints[their_dof_g].first;
2482 
2483  for (unsigned int js = 0; js != n_side_dofs; ++js)
2484  {
2485  const unsigned int j = my_side_dofs[js];
2486  const dof_id_type my_dof_g = my_dof_indices[j];
2487  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2488 
2489  if (their_constraint_row.count(my_dof_g))
2490  recursive_constraint[js] = true;
2491  }
2492  }
2493  */
2494 
2495  for (unsigned int js = 0; js != n_side_dofs; ++js)
2496  {
2497  // FIXME: old code path
2498  // if (recursive_constraint[js])
2499  // continue;
2500 
2501  const unsigned int j = my_side_dofs[js];
2502  const dof_id_type my_dof_g = my_dof_indices[j];
2503  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2504 
2505  // FIXME: new code path
2506  if (!my_constrained_dofs.count(my_dof_g))
2507  continue;
2508 
2509  DofConstraintRow* constraint_row;
2510 
2511  // we may be running constraint methods concurretly
2512  // on multiple threads, so we need a lock to
2513  // ensure that this constraint is "ours"
2514  {
2515  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2516 
2517  if (dof_map.is_constrained_dof(my_dof_g))
2518  continue;
2519 
2520  constraint_row = &(constraints[my_dof_g]);
2521  libmesh_assert(constraint_row->empty());
2522  }
2523 
2524  for (unsigned int is = 0; is != n_side_dofs; ++is)
2525  {
2526  const unsigned int i = neigh_side_dofs[is];
2527  const dof_id_type their_dof_g = neigh_dof_indices[i];
2528  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2529 
2530  // Periodic constraints should never be
2531  // self-constraints
2532  // libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
2533 
2534  const Real their_dof_value = Ue[is](js);
2535 
2536  if (their_dof_g == my_dof_g)
2537  {
2538  libmesh_assert_less (std::abs(their_dof_value-1.), 1.e-5);
2539  for (unsigned int k = 0; k != n_side_dofs; ++k)
2540  libmesh_assert(k == is || std::abs(Ue[k](js)) < 1.e-5);
2541  continue;
2542  }
2543 
2544  if (std::abs(their_dof_value) < 10*TOLERANCE)
2545  continue;
2546 
2547  constraint_row->insert(std::make_pair(their_dof_g,
2548  their_dof_value));
2549  }
2550  }
2551  }
2552  // p refinement constraints:
2553  // constrain dofs shared between
2554  // active elements and neighbors with
2555  // lower polynomial degrees
2556 #ifdef LIBMESH_ENABLE_AMR
2557  const unsigned int min_p_level =
2558  neigh->min_p_level_by_neighbor(elem, elem->p_level());
2559  if (min_p_level < elem->p_level())
2560  {
2561  // Adaptive p refinement of non-hierarchic bases will
2562  // require more coding
2563  libmesh_assert(my_fe->is_hierarchic());
2564  dof_map.constrain_p_dofs(variable_number, elem,
2565  s, min_p_level);
2566  }
2567 #endif // #ifdef LIBMESH_ENABLE_AMR
2568  }
2569  }
2570  }
2571 }
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 1063 of file fe_abstract.C.

References libMesh::Elem::active(), libMesh::PeriodicBoundaries::boundary(), libMesh::MeshBase::boundary_info, libMesh::Elem::build_side(), libMesh::Elem::default_order(), libMesh::Elem::dim(), libMesh::FEAbstract::fe_type, libMesh::PeriodicBoundaryBase::get_corresponding_pos(), libMesh::invalid_uint, libMesh::FEInterface::inverse_map(), libMeshEnums::LAGRANGE, libMesh::Elem::level(), libMesh::libmesh_assert(), libMesh::FEInterface::n_dofs(), libMesh::Elem::n_sides(), libMesh::PeriodicBoundaries::neighbor(), libMesh::Elem::neighbor(), libMesh::PeriodicBoundaryBase::pairedboundary, libMesh::Real, libMesh::FEInterface::shape(), and libMesh::Threads::spin_mtx.

1068 {
1069  // Only bother if we truly have periodic boundaries
1070  if (boundaries.empty())
1071  return;
1072 
1073  libmesh_assert(elem);
1074 
1075  // Only constrain active elements with this method
1076  if (!elem->active())
1077  return;
1078 
1079  const unsigned int Dim = elem->dim();
1080 
1081  // We currently always use LAGRANGE mappings for geometry
1082  const FEType fe_type(elem->default_order(), LAGRANGE);
1083 
1084  std::vector<const Node*> my_nodes, neigh_nodes;
1085 
1086  // Look at the element faces. Check to see if we need to
1087  // build constraints.
1088  for (unsigned int s=0; s<elem->n_sides(); s++)
1089  {
1090  if (elem->neighbor(s))
1091  continue;
1092 
1093  const std::vector<boundary_id_type>& bc_ids = mesh.boundary_info->boundary_ids (elem, s);
1094  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1095  {
1096  const boundary_id_type boundary_id = *id_it;
1097  const PeriodicBoundaryBase *periodic = boundaries.boundary(boundary_id);
1098  if (periodic)
1099  {
1100  libmesh_assert(point_locator);
1101 
1102  // Get pointers to the element's neighbor.
1103  const Elem* neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s);
1104 
1105  // h refinement constraints:
1106  // constrain dofs shared between
1107  // this element and ones as coarse
1108  // as or coarser than this element.
1109  if (neigh->level() <= elem->level())
1110  {
1111  unsigned int s_neigh =
1112  mesh.boundary_info->side_with_boundary_id (neigh, periodic->pairedboundary);
1113  libmesh_assert_not_equal_to (s_neigh, libMesh::invalid_uint);
1114 
1115 #ifdef LIBMESH_ENABLE_AMR
1116  libmesh_assert(neigh->active());
1117 #endif // #ifdef LIBMESH_ENABLE_AMR
1118 
1119  const AutoPtr<Elem> my_side (elem->build_side(s));
1120  const AutoPtr<Elem> neigh_side (neigh->build_side(s_neigh));
1121 
1122  const unsigned int n_side_nodes = my_side->n_nodes();
1123 
1124  my_nodes.clear();
1125  my_nodes.reserve (n_side_nodes);
1126  neigh_nodes.clear();
1127  neigh_nodes.reserve (n_side_nodes);
1128 
1129  for (unsigned int n=0; n != n_side_nodes; ++n)
1130  my_nodes.push_back(my_side->get_node(n));
1131 
1132  for (unsigned int n=0; n != n_side_nodes; ++n)
1133  neigh_nodes.push_back(neigh_side->get_node(n));
1134 
1135  // Make sure we're not adding recursive constraints
1136  // due to the redundancy in the way we add periodic
1137  // boundary constraints, or adding constraints to
1138  // nodes that already have AMR constraints
1139  std::vector<bool> skip_constraint(n_side_nodes, false);
1140 
1141  for (unsigned int my_side_n=0;
1142  my_side_n < n_side_nodes;
1143  my_side_n++)
1144  {
1145  libmesh_assert_less (my_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
1146 
1147  const Node* my_node = my_nodes[my_side_n];
1148 
1149  // Figure out where my node lies on their reference element.
1150  const Point neigh_point = periodic->get_corresponding_pos(*my_node);
1151 
1152  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
1153  neigh_side.get(),
1154  neigh_point);
1155 
1156  // If we've already got a constraint on this
1157  // node, then the periodic constraint is
1158  // redundant
1159  {
1160  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1161 
1162  if (constraints.count(my_node))
1163  {
1164  skip_constraint[my_side_n] = true;
1165  continue;
1166  }
1167  }
1168 
1169  // Compute the neighbors's side shape function values.
1170  for (unsigned int their_side_n=0;
1171  their_side_n < n_side_nodes;
1172  their_side_n++)
1173  {
1174  libmesh_assert_less (their_side_n, FEInterface::n_dofs(Dim-1, fe_type, neigh_side->type()));
1175 
1176  const Node* their_node = neigh_nodes[their_side_n];
1177 
1178  // If there's a constraint on an opposing node,
1179  // we need to see if it's constrained by
1180  // *our side* making any periodic constraint
1181  // on us recursive
1182  {
1183  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1184 
1185  if (!constraints.count(their_node))
1186  continue;
1187 
1188  const NodeConstraintRow& their_constraint_row =
1189  constraints[their_node].first;
1190 
1191  for (unsigned int orig_side_n=0;
1192  orig_side_n < n_side_nodes;
1193  orig_side_n++)
1194  {
1195  libmesh_assert_less (orig_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
1196 
1197  const Node* orig_node = my_nodes[orig_side_n];
1198 
1199  if (their_constraint_row.count(orig_node))
1200  skip_constraint[orig_side_n] = true;
1201  }
1202  }
1203  }
1204  }
1205  for (unsigned int my_side_n=0;
1206  my_side_n < n_side_nodes;
1207  my_side_n++)
1208  {
1209  libmesh_assert_less (my_side_n, FEInterface::n_dofs(Dim-1, fe_type, my_side->type()));
1210 
1211  if (skip_constraint[my_side_n])
1212  continue;
1213 
1214  const Node* my_node = my_nodes[my_side_n];
1215 
1216  // Figure out where my node lies on their reference element.
1217  const Point neigh_point = periodic->get_corresponding_pos(*my_node);
1218 
1219  // Figure out where my node lies on their reference element.
1220  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
1221  neigh_side.get(),
1222  neigh_point);
1223 
1224  for (unsigned int their_side_n=0;
1225  their_side_n < n_side_nodes;
1226  their_side_n++)
1227  {
1228  libmesh_assert_less (their_side_n, FEInterface::n_dofs(Dim-1, fe_type, neigh_side->type()));
1229 
1230  const Node* their_node = neigh_nodes[their_side_n];
1231  libmesh_assert(their_node);
1232 
1233  const Real their_value = FEInterface::shape(Dim-1,
1234  fe_type,
1235  neigh_side->type(),
1236  their_side_n,
1237  mapped_point);
1238 
1239  // since we may be running this method concurretly
1240  // on multiple threads we need to acquire a lock
1241  // before modifying the shared constraint_row object.
1242  {
1243  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1244 
1245  NodeConstraintRow& constraint_row =
1246  constraints[my_node].first;
1247 
1248  constraint_row.insert(std::make_pair(their_node,
1249  their_value));
1250  }
1251  }
1252  }
1253  }
1254  }
1255  }
1256  }
1257 }
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 1582 of file fe_base.C.

References std::abs(), libMesh::Elem::active(), libMeshEnums::C_ONE, libMeshEnums::C_ZERO, libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::DofMap::constrain_p_dofs(), libMesh::FEType::default_quadrature_order(), libMesh::Elem::dim(), libMeshEnums::DISCONTINUOUS, libMesh::DofMap::dof_indices(), libMesh::FEInterface::dofs_on_side(), libMesh::TensorTools::inner_product(), libMesh::DofObject::invalid_id, libMesh::FEInterface::inverse_map(), libMesh::DofMap::is_constrained_dof(), libMesh::Elem::level(), libMesh::libmesh_assert(), std::min(), libMesh::Elem::min_p_level_by_neighbor(), libMesh::Elem::n_neighbors(), libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::Elem::neighbor(), libMesh::Elem::p_level(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::Threads::spin_mtx, libMesh::TOLERANCE, libMesh::DofMap::variable_type(), and libMesh::Elem::which_neighbor_am_i().

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

1586 {
1587  libmesh_assert(elem);
1588 
1589  const unsigned int Dim = elem->dim();
1590 
1591  // Only constrain elements in 2,3D.
1592  if (Dim == 1)
1593  return;
1594 
1595  // Only constrain active elements with this method
1596  if (!elem->active())
1597  return;
1598 
1599  const FEType& base_fe_type = dof_map.variable_type(variable_number);
1600 
1601  // Construct FE objects for this element and its neighbors.
1603  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1604  const FEContinuity cont = my_fe->get_continuity();
1605 
1606  // We don't need to constrain discontinuous elements
1607  if (cont == DISCONTINUOUS)
1608  return;
1609  libmesh_assert (cont == C_ZERO || cont == C_ONE);
1610 
1612  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1613 
1614  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1615  my_fe->attach_quadrature_rule (&my_qface);
1616  std::vector<Point> neigh_qface;
1617 
1618  const std::vector<Real>& JxW = my_fe->get_JxW();
1619  const std::vector<Point>& q_point = my_fe->get_xyz();
1620  const std::vector<std::vector<OutputShape> >& phi = my_fe->get_phi();
1621  const std::vector<std::vector<OutputShape> >& neigh_phi =
1622  neigh_fe->get_phi();
1623  const std::vector<Point> *face_normals = NULL;
1624  const std::vector<std::vector<OutputGradient> > *dphi = NULL;
1625  const std::vector<std::vector<OutputGradient> > *neigh_dphi = NULL;
1626 
1627  std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1628  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1629 
1630  if (cont != C_ZERO)
1631  {
1632  const std::vector<Point>& ref_face_normals =
1633  my_fe->get_normals();
1634  face_normals = &ref_face_normals;
1635  const std::vector<std::vector<OutputGradient> >& ref_dphi =
1636  my_fe->get_dphi();
1637  dphi = &ref_dphi;
1638  const std::vector<std::vector<OutputGradient> >& ref_neigh_dphi =
1639  neigh_fe->get_dphi();
1640  neigh_dphi = &ref_neigh_dphi;
1641  }
1642 
1643  DenseMatrix<Real> Ke;
1644  DenseVector<Real> Fe;
1645  std::vector<DenseVector<Real> > Ue;
1646 
1647  // Look at the element faces. Check to see if we need to
1648  // build constraints.
1649  for (unsigned int s=0; s<elem->n_sides(); s++)
1650  if (elem->neighbor(s) != NULL)
1651  {
1652  // Get pointers to the element's neighbor.
1653  const Elem* neigh = elem->neighbor(s);
1654 
1655  // h refinement constraints:
1656  // constrain dofs shared between
1657  // this element and ones coarser
1658  // than this element.
1659  if (neigh->level() < elem->level())
1660  {
1661  unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
1662  libmesh_assert_less (s_neigh, neigh->n_neighbors());
1663 
1664  // Find the minimum p level; we build the h constraint
1665  // matrix with this and then constrain away all higher p
1666  // DoFs.
1667  libmesh_assert(neigh->active());
1668  const unsigned int min_p_level =
1669  std::min(elem->p_level(), neigh->p_level());
1670 
1671  // we may need to make the FE objects reinit with the
1672  // minimum shared p_level
1673  // FIXME - I hate using const_cast<> and avoiding
1674  // accessor functions; there's got to be a
1675  // better way to do this!
1676  const unsigned int old_elem_level = elem->p_level();
1677  if (old_elem_level != min_p_level)
1678  (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
1679  const unsigned int old_neigh_level = neigh->p_level();
1680  if (old_neigh_level != min_p_level)
1681  (const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
1682 
1683  my_fe->reinit(elem, s);
1684 
1685  // This function gets called element-by-element, so there
1686  // will be a lot of memory allocation going on. We can
1687  // at least minimize this for the case of the dof indices
1688  // by efficiently preallocating the requisite storage.
1689  // n_nodes is not necessarily n_dofs, but it is better
1690  // than nothing!
1691  my_dof_indices.reserve (elem->n_nodes());
1692  neigh_dof_indices.reserve (neigh->n_nodes());
1693 
1694  dof_map.dof_indices (elem, my_dof_indices,
1695  variable_number);
1696  dof_map.dof_indices (neigh, neigh_dof_indices,
1697  variable_number);
1698 
1699  const unsigned int n_qp = my_qface.n_points();
1700 
1701  FEInterface::inverse_map (Dim, base_fe_type, neigh,
1702  q_point, neigh_qface);
1703 
1704  neigh_fe->reinit(neigh, &neigh_qface);
1705 
1706  // We're only concerned with DOFs whose values (and/or first
1707  // derivatives for C1 elements) are supported on side nodes
1708  FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
1709  FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
1710 
1711  // We're done with functions that examine Elem::p_level(),
1712  // so let's unhack those levels
1713  if (elem->p_level() != old_elem_level)
1714  (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
1715  if (neigh->p_level() != old_neigh_level)
1716  (const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
1717 
1718  const unsigned int n_side_dofs =
1719  libmesh_cast_int<unsigned int>(my_side_dofs.size());
1720  libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
1721 
1722  Ke.resize (n_side_dofs, n_side_dofs);
1723  Ue.resize(n_side_dofs);
1724 
1725  // Form the projection matrix, (inner product of fine basis
1726  // functions against fine test functions)
1727  for (unsigned int is = 0; is != n_side_dofs; ++is)
1728  {
1729  const unsigned int i = my_side_dofs[is];
1730  for (unsigned int js = 0; js != n_side_dofs; ++js)
1731  {
1732  const unsigned int j = my_side_dofs[js];
1733  for (unsigned int qp = 0; qp != n_qp; ++qp)
1734  {
1735  Ke(is,js) += JxW[qp] * TensorTools::inner_product(phi[i][qp], phi[j][qp]);
1736  if (cont != C_ZERO)
1737  Ke(is,js) += JxW[qp] *
1738  TensorTools::inner_product((*dphi)[i][qp] *
1739  (*face_normals)[qp],
1740  (*dphi)[j][qp] *
1741  (*face_normals)[qp]);
1742  }
1743  }
1744  }
1745 
1746  // Form the right hand sides, (inner product of coarse basis
1747  // functions against fine test functions)
1748  for (unsigned int is = 0; is != n_side_dofs; ++is)
1749  {
1750  const unsigned int i = neigh_side_dofs[is];
1751  Fe.resize (n_side_dofs);
1752  for (unsigned int js = 0; js != n_side_dofs; ++js)
1753  {
1754  const unsigned int j = my_side_dofs[js];
1755  for (unsigned int qp = 0; qp != n_qp; ++qp)
1756  {
1757  Fe(js) += JxW[qp] *
1758  TensorTools::inner_product(neigh_phi[i][qp],
1759  phi[j][qp]);
1760  if (cont != C_ZERO)
1761  Fe(js) += JxW[qp] *
1762  TensorTools::inner_product((*neigh_dphi)[i][qp] *
1763  (*face_normals)[qp],
1764  (*dphi)[j][qp] *
1765  (*face_normals)[qp]);
1766  }
1767  }
1768  Ke.cholesky_solve(Fe, Ue[is]);
1769  }
1770 
1771  for (unsigned int js = 0; js != n_side_dofs; ++js)
1772  {
1773  const unsigned int j = my_side_dofs[js];
1774  const dof_id_type my_dof_g = my_dof_indices[j];
1775  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
1776 
1777  // Hunt for "constraining against myself" cases before
1778  // we bother creating a constraint row
1779  bool self_constraint = false;
1780  for (unsigned int is = 0; is != n_side_dofs; ++is)
1781  {
1782  const unsigned int i = neigh_side_dofs[is];
1783  const dof_id_type their_dof_g = neigh_dof_indices[i];
1784  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1785 
1786  if (their_dof_g == my_dof_g)
1787  {
1788 #ifndef NDEBUG
1789  const Real their_dof_value = Ue[is](js);
1790  libmesh_assert_less (std::abs(their_dof_value-1.),
1791  10*TOLERANCE);
1792 
1793  for (unsigned int k = 0; k != n_side_dofs; ++k)
1794  libmesh_assert(k == is ||
1795  std::abs(Ue[k](js)) <
1796  10*TOLERANCE);
1797 #endif
1798 
1799  self_constraint = true;
1800  break;
1801  }
1802  }
1803 
1804  if (self_constraint)
1805  continue;
1806 
1807  DofConstraintRow* constraint_row;
1808 
1809  // we may be running constraint methods concurretly
1810  // on multiple threads, so we need a lock to
1811  // ensure that this constraint is "ours"
1812  {
1813  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1814 
1815  if (dof_map.is_constrained_dof(my_dof_g))
1816  continue;
1817 
1818  constraint_row = &(constraints[my_dof_g]);
1819  libmesh_assert(constraint_row->empty());
1820  }
1821 
1822  for (unsigned int is = 0; is != n_side_dofs; ++is)
1823  {
1824  const unsigned int i = neigh_side_dofs[is];
1825  const dof_id_type their_dof_g = neigh_dof_indices[i];
1826  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1827  libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
1828 
1829  const Real their_dof_value = Ue[is](js);
1830 
1831  if (std::abs(their_dof_value) < 10*TOLERANCE)
1832  continue;
1833 
1834  constraint_row->insert(std::make_pair(their_dof_g,
1835  their_dof_value));
1836  }
1837  }
1838  }
1839  // p refinement constraints:
1840  // constrain dofs shared between
1841  // active elements and neighbors with
1842  // lower polynomial degrees
1843  const unsigned int min_p_level =
1844  neigh->min_p_level_by_neighbor(elem, elem->p_level());
1845  if (min_p_level < elem->p_level())
1846  {
1847  // Adaptive p refinement of non-hierarchic bases will
1848  // require more coding
1849  libmesh_assert(my_fe->is_hierarchic());
1850  dof_map.constrain_p_dofs(variable_number, elem,
1851  s, min_p_level);
1852  }
1853  }
1854 }
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. Overloaded method from the FEBase version.

Reimplemented from libMesh::FEGenericBase< T >.

Definition at line 892 of file inf_fe.C.

References libMesh::dim, libMesh::libmesh_assert(), libMesh::START_LOG(), and libMesh::STOP_LOG().

893 {
895 
896 
897 
898  // Start logging the overall computation of shape functions
899  START_LOG("compute_shape_functions()", "InfFE");
900 
901 
902  const unsigned int n_total_qp = _n_total_qp;
903 
904 
905  //-------------------------------------------------------------------------
906  // Compute the shape function values (and derivatives)
907  // at the Quadrature points. Note that the actual values
908  // have already been computed via init_shape_functions
909 
910  // Compute the value of the derivative shape function i at quadrature point p
911  switch (dim)
912  {
913 
914  case 1:
915  {
916  libmesh_not_implemented();
917  break;
918  }
919 
920  case 2:
921  {
922  libmesh_not_implemented();
923  break;
924  }
925 
926  case 3:
927  {
928  const std::vector<Real>& dxidx_map = this->_fe_map->get_dxidx();
929  const std::vector<Real>& dxidy_map = this->_fe_map->get_dxidy();
930  const std::vector<Real>& dxidz_map = this->_fe_map->get_dxidz();
931 
932  const std::vector<Real>& detadx_map = this->_fe_map->get_detadx();
933  const std::vector<Real>& detady_map = this->_fe_map->get_detady();
934  const std::vector<Real>& detadz_map = this->_fe_map->get_detadz();
935 
936  const std::vector<Real>& dzetadx_map = this->_fe_map->get_dzetadx();
937  const std::vector<Real>& dzetady_map = this->_fe_map->get_dzetady();
938  const std::vector<Real>& dzetadz_map = this->_fe_map->get_dzetadz();
939 
940  // These are _all_ shape functions of this infinite element
941  for (unsigned int i=0; i<phi.size(); i++)
942  for (unsigned int p=0; p<n_total_qp; p++)
943  {
944  // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
945  dphi[i][p](0) =
946  dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
947  dphideta[i][p]*detadx_map[p] +
948  dphidzeta[i][p]*dzetadx_map[p]);
949 
950  // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
951  dphi[i][p](1) =
952  dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
953  dphideta[i][p]*detady_map[p] +
954  dphidzeta[i][p]*dzetady_map[p]);
955 
956  // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
957  dphi[i][p](2) =
958  dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
959  dphideta[i][p]*detadz_map[p] +
960  dphidzeta[i][p]*dzetadz_map[p]);
961  }
962 
963 
964  // This is the derivative of the phase term of this infinite element
965  for (unsigned int p=0; p<n_total_qp; p++)
966  {
967  // the derivative of the phase term
968  dphase[p](0) = (dphasedxi[p] * dxidx_map[p] +
969  dphasedeta[p] * detadx_map[p] +
970  dphasedzeta[p] * dzetadx_map[p]);
971 
972  dphase[p](1) = (dphasedxi[p] * dxidy_map[p] +
973  dphasedeta[p] * detady_map[p] +
974  dphasedzeta[p] * dzetady_map[p]);
975 
976  dphase[p](2) = (dphasedxi[p] * dxidz_map[p] +
977  dphasedeta[p] * detadz_map[p] +
978  dphasedzeta[p] * dzetadz_map[p]);
979 
980  // the derivative of the radial weight - varies only in radial direction,
981  // therefore dweightdxi = dweightdeta = 0.
982  dweight[p](0) = dweightdv[p] * dzetadx_map[p];
983 
984  dweight[p](1) = dweightdv[p] * dzetady_map[p];
985 
986  dweight[p](2) = dweightdv[p] * dzetadz_map[p];
987 
988  }
989 
990  break;
991  }
992 
993 
994 
995  default:
996  {
997  libmesh_error();
998  }
999  }
1000 
1001 
1002 
1003  // Stop logging the overall computation of shape functions
1004  STOP_LOG("compute_shape_functions()", "InfFE");
1005 
1006 }
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 727 of file inf_fe_static.C.

References libMeshEnums::INFEDGE2, libMeshEnums::INFHEX16, libMeshEnums::INFHEX18, libMeshEnums::INFHEX8, libMeshEnums::INFPRISM12, libMeshEnums::INFPRISM6, libMeshEnums::INFQUAD4, libMeshEnums::INFQUAD6, libMesh::invalid_uint, libMesh::FEInterface::n_dofs_at_node(), libMesh::FEInterface::n_dofs_per_elem(), and libMesh::FEType::radial_order.

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

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

107 {
108  _enable_print_counter = false;
109  return;
110 }
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 = NULL,
const std::vector< Real > *const  weights = NULL 
)
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 117 of file inf_fe_boundary.C.

References libMesh::err.

122 {
123  // We don't do this for 1D elements!
124  //libmesh_assert_not_equal_to (Dim, 1);
125 
126  libMesh::err << "ERROR: Edge conditions for infinite elements "
127  << "not implemented!" << std::endl;
128  libmesh_error();
129 
130  if (pts != NULL)
131  {
132  libMesh::err << "ERROR: User-specified points for infinite elements "
133  << "not implemented!" << std::endl;
134  libmesh_error();
135  }
136 }
void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

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

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

Definition at line 81 of file inf_fe_map_eval.C.

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

Definition at line 82 of file inf_fe_map_eval.C.

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

Definition at line 83 of file inf_fe_map_eval.C.

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

Definition at line 322 of file inf_fe_legendre_eval.C.

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

Definition at line 323 of file inf_fe_legendre_eval.C.

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

Definition at line 324 of file inf_fe_legendre_eval.C.

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

Definition at line 467 of file inf_fe_jacobi_30_00_eval.C.

467 { 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 468 of file inf_fe_jacobi_30_00_eval.C.

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

469 { 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 473 of file inf_fe_jacobi_20_00_eval.C.

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

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

475 { return jacobi_20_00_eval(v, i); }
template<unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
static Real libMesh::InfFE< friend_Dim, friend_T_radial, friend_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.

Specialized for T_radial=INFINITE_MAP, this function returns the value of the $ i^{th} $ mapping shape function in radial direction evaluated at v. 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< friend_Dim, friend_T_radial, friend_T_map >::compute_data(), libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::init_radial_shape_functions(), libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::inverse_map(), and libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::shape().

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

Definition at line 2649 of file inf_fe_lagrange_eval.C.

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

Definition at line 2650 of file inf_fe_lagrange_eval.C.

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

Definition at line 2651 of file inf_fe_lagrange_eval.C.

2651 { 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 87 of file inf_fe_map_eval.C.

87 { 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 88 of file inf_fe_map_eval.C.

88 { 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 89 of file inf_fe_map_eval.C.

89 { 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 328 of file inf_fe_legendre_eval.C.

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

Definition at line 329 of file inf_fe_legendre_eval.C.

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

Definition at line 330 of file inf_fe_legendre_eval.C.

330 { 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 473 of file inf_fe_jacobi_30_00_eval.C.

473 { 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 474 of file inf_fe_jacobi_30_00_eval.C.

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

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

479 { 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 480 of file inf_fe_jacobi_20_00_eval.C.

480 { 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 481 of file inf_fe_jacobi_20_00_eval.C.

481 { return jacobi_20_00_eval_deriv(v, i); }
template<unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
static Real libMesh::InfFE< friend_Dim, friend_T_radial, friend_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< friend_Dim, friend_T_radial, friend_T_map >::init_radial_shape_functions(), and libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::inverse_map().

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

Definition at line 2655 of file inf_fe_lagrange_eval.C.

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

Definition at line 2656 of file inf_fe_lagrange_eval.C.

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

Definition at line 2657 of file inf_fe_lagrange_eval.C.

2657 { return lagrange_eval_deriv(v, o, i); }
template<unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
virtual FEContinuity libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::get_continuity ( ) const
inlinevirtual
Returns
the continuity of the element.

Implements libMesh::FEAbstract.

Definition at line 340 of file inf_fe.h.

References libMeshEnums::C_ZERO.

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

Definition at line 238 of file fe_base.h.

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

const std::vector<Real>& libMesh::FEAbstract::get_curvatures ( ) const
inlineinherited
Returns
the curvatures for use in face integration.

Definition at line 380 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

381  { return this->_fe_map->get_curvatures();}
template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phideta2 ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 384 of file fe_base.h.

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

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidetadzeta ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 392 of file fe_base.h.

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

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidx2 ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points.

Definition at line 312 of file fe_base.h.

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidxdy ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points.

Definition at line 320 of file fe_base.h.

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidxdz ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points.

Definition at line 328 of file fe_base.h.

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidxi2 ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 360 of file fe_base.h.

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

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidxideta ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 368 of file fe_base.h.

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

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidxidzeta ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 376 of file fe_base.h.

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

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidy2 ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points.

Definition at line 336 of file fe_base.h.

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidydz ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points.

Definition at line 344 of file fe_base.h.

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidz2 ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points.

Definition at line 352 of file fe_base.h.

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_d2phidzeta2 ( ) const
inlineinherited
Returns
the shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 400 of file fe_base.h.

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

const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdeta2 ( ) const
inlineinherited
Returns
the second partial derivatives in eta.

Definition at line 267 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

268  { return this->_fe_map->get_d2xyzdeta2(); }
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdetadzeta ( ) const
inlineinherited
Returns
the second partial derivatives in eta-zeta.

Definition at line 297 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

298  { return this->_fe_map->get_d2xyzdetadzeta(); }
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdxi2 ( ) const
inlineinherited
Returns
the second partial derivatives in xi.

Definition at line 261 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

262  { return this->_fe_map->get_d2xyzdxi2(); }
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdxideta ( ) const
inlineinherited
Returns
the second partial derivatives in xi-eta.

Definition at line 283 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

284  { return this->_fe_map->get_d2xyzdxideta(); }
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdxidzeta ( ) const
inlineinherited
Returns
the second partial derivatives in xi-zeta.

Definition at line 291 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

292  { return this->_fe_map->get_d2xyzdxidzeta(); }
const std::vector<RealGradient>& libMesh::FEAbstract::get_d2xyzdzeta2 ( ) const
inlineinherited
Returns
the second partial derivatives in zeta.

Definition at line 275 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

276  { return this->_fe_map->get_d2xyzdzeta2(); }
const std::vector<Real>& libMesh::FEAbstract::get_detadx ( ) const
inlineinherited
Returns
the deta/dx entry in the transformation matrix from physical to local coordinates.

Definition at line 327 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

328  { return this->_fe_map->get_detadx(); }
const std::vector<Real>& libMesh::FEAbstract::get_detady ( ) const
inlineinherited
Returns
the deta/dy entry in the transformation matrix from physical to local coordinates.

Definition at line 334 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

335  { return this->_fe_map->get_detady(); }
const std::vector<Real>& libMesh::FEAbstract::get_detadz ( ) const
inlineinherited
Returns
the deta/dz entry in the transformation matrix from physical to local coordinates.

Definition at line 341 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

342  { return this->_fe_map->get_detadz(); }
template<typename T>
const std::vector<std::vector<OutputDivergence> >& libMesh::FEGenericBase< T >::get_div_phi ( ) const
inlineinherited
Returns
the divergence of the shape function at the quadrature points.

Definition at line 246 of file fe_base.h.

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

template<typename T>
const std::vector<OutputGradient>& libMesh::FEGenericBase< T >::get_dphase ( ) const
inlineinherited
Returns
the global first derivative of the phase term which is used in infinite elements, evaluated at the quadrature points.

In case of the general finite element class FE this field is initialized to all zero, so that the variational formulation for an infinite element returns correct element matrices for a mesh using both finite and infinite elements.

Definition at line 418 of file fe_base.h.

419  { return dphase; }
template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_dphideta ( ) const
inlineinherited
template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_dphidx ( ) const
inlineinherited
Returns
the shape function x-derivative at the quadrature points.

Definition at line 254 of file fe_base.h.

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_dphidxi ( ) const
inlineinherited
template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_dphidy ( ) const
inlineinherited
Returns
the shape function y-derivative at the quadrature points.

Definition at line 262 of file fe_base.h.

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_dphidz ( ) const
inlineinherited
Returns
the shape function z-derivative at the quadrature points.

Definition at line 270 of file fe_base.h.

template<typename T>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< T >::get_dphidzeta ( ) const
inlineinherited
const std::vector<Real>& libMesh::FEAbstract::get_dxidx ( ) const
inlineinherited
Returns
the dxi/dx entry in the transformation matrix from physical to local coordinates.

Definition at line 306 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

307  { return this->_fe_map->get_dxidx(); }
const std::vector<Real>& libMesh::FEAbstract::get_dxidy ( ) const
inlineinherited
Returns
the dxi/dy entry in the transformation matrix from physical to local coordinates.

Definition at line 313 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

314  { return this->_fe_map->get_dxidy(); }
const std::vector<Real>& libMesh::FEAbstract::get_dxidz ( ) const
inlineinherited
Returns
the dxi/dz entry in the transformation matrix from physical to local coordinates.

Definition at line 320 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

321  { return this->_fe_map->get_dxidz(); }
const std::vector<RealGradient>& libMesh::FEAbstract::get_dxyzdeta ( ) const
inlineinherited
Returns
the element tangents in eta-direction at the quadrature points.

Definition at line 248 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

249  { return this->_fe_map->get_dxyzdeta(); }
const std::vector<RealGradient>& libMesh::FEAbstract::get_dxyzdxi ( ) const
inlineinherited
Returns
the element tangents in xi-direction at the quadrature points.

Definition at line 241 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

242  { return this->_fe_map->get_dxyzdxi(); }
const std::vector<RealGradient>& libMesh::FEAbstract::get_dxyzdzeta ( ) const
inlineinherited
Returns
the element tangents in zeta-direction at the quadrature points.

Definition at line 255 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

256  { return _fe_map->get_dxyzdzeta(); }
const std::vector<Real>& libMesh::FEAbstract::get_dzetadx ( ) const
inlineinherited
Returns
the dzeta/dx entry in the transformation matrix from physical to local coordinates.

Definition at line 348 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

349  { return this->_fe_map->get_dzetadx(); }
const std::vector<Real>& libMesh::FEAbstract::get_dzetady ( ) const
inlineinherited
Returns
the dzeta/dy entry in the transformation matrix from physical to local coordinates.

Definition at line 355 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

356  { return this->_fe_map->get_dzetady(); }
const std::vector<Real>& libMesh::FEAbstract::get_dzetadz ( ) const
inlineinherited
Returns
the dzeta/dz entry in the transformation matrix from physical to local coordinates.

Definition at line 362 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

363  { return this->_fe_map->get_dzetadz(); }
FEFamily libMesh::FEAbstract::get_family ( ) const
inlineinherited
Returns
the finite element family of this element.

Definition at line 439 of file fe_abstract.h.

References libMesh::FEType::family, and libMesh::FEAbstract::fe_type.

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

439 { return fe_type.family; }
FEType libMesh::FEAbstract::get_fe_type ( ) const
inlineinherited
Returns
the FE Type (approximation order and family) of the finite element.

Definition at line 418 of file fe_abstract.h.

References libMesh::FEAbstract::fe_type.

Referenced by libMesh::FEMContext::build_new_fe(), libMesh::H1FETransformation< T >::map_phi(), libMesh::HCurlFETransformation< T >::map_phi(), and libMesh::ProjectFEMSolution::operator()().

418 { return fe_type; }
std::string libMesh::ReferenceCounter::get_info ( )
staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

References libMesh::ReferenceCounter::_counts, and libMesh::Quality::name().

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

48 {
49 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50 
51  std::ostringstream oss;
52 
53  oss << '\n'
54  << " ---------------------------------------------------------------------------- \n"
55  << "| Reference count information |\n"
56  << " ---------------------------------------------------------------------------- \n";
57 
58  for (Counts::iterator it = _counts.begin();
59  it != _counts.end(); ++it)
60  {
61  const std::string name(it->first);
62  const unsigned int creations = it->second.first;
63  const unsigned int destructions = it->second.second;
64 
65  oss << "| " << name << " reference count information:\n"
66  << "| Creations: " << creations << '\n'
67  << "| Destructions: " << destructions << '\n';
68  }
69 
70  oss << " ---------------------------------------------------------------------------- \n";
71 
72  return oss.str();
73 
74 #else
75 
76  return "";
77 
78 #endif
79 }
const std::vector<Real>& libMesh::FEAbstract::get_JxW ( ) const
inlineinherited
Returns
the element Jacobian times the quadrature weight for each quadrature point.

Definition at line 234 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

Referenced by libMesh::ExactErrorEstimator::find_squared_element_error(), libMesh::FEMSystem::init_context(), and libMesh::ProjectFEMSolution::operator()().

235  { return this->_fe_map->get_JxW(); }
const std::vector<Point>& libMesh::FEAbstract::get_normals ( ) const
inlineinherited
Returns
the normal vectors for face integration.

Definition at line 374 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

375  { return this->_fe_map->get_normals(); }
Order libMesh::FEAbstract::get_order ( ) const
inlineinherited
Returns
the approximation order of the finite element.

Definition at line 423 of file fe_abstract.h.

References libMesh::FEAbstract::_p_level, libMesh::FEAbstract::fe_type, and libMesh::FEType::order.

423 { return static_cast<Order>(fe_type.order + _p_level); }
unsigned int libMesh::FEAbstract::get_p_level ( ) const
inlineinherited
Returns
the p refinement level that the current shape functions have been calculated for.

Definition at line 413 of file fe_abstract.h.

References libMesh::FEAbstract::_p_level.

413 { return _p_level; }
void libMesh::FEAbstract::get_refspace_nodes ( const ElemType  t,
std::vector< Point > &  nodes 
)
staticinherited

returns the reference space nodes coordinates given the element type

Definition at line 415 of file fe_abstract.C.

References libMeshEnums::EDGE2, libMeshEnums::EDGE3, libMesh::err, libMeshEnums::HEX20, libMeshEnums::HEX27, libMeshEnums::HEX8, libMeshEnums::PRISM15, libMeshEnums::PRISM18, libMeshEnums::PRISM6, libMeshEnums::PYRAMID14, libMeshEnums::PYRAMID5, libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, libMeshEnums::TET10, libMeshEnums::TET4, libMeshEnums::TRI3, and libMeshEnums::TRI6.

416 {
417  switch(itemType)
418  {
419  case EDGE2:
420  {
421  nodes.resize(2);
422  nodes[0] = Point (-1.,0.,0.);
423  nodes[1] = Point (1.,0.,0.);
424  return;
425  }
426  case EDGE3:
427  {
428  nodes.resize(3);
429  nodes[0] = Point (-1.,0.,0.);
430  nodes[1] = Point (1.,0.,0.);
431  nodes[2] = Point (0.,0.,0.);
432  return;
433  }
434  case TRI3:
435  {
436  nodes.resize(3);
437  nodes[0] = Point (0.,0.,0.);
438  nodes[1] = Point (1.,0.,0.);
439  nodes[2] = Point (0.,1.,0.);
440  return;
441  }
442  case TRI6:
443  {
444  nodes.resize(6);
445  nodes[0] = Point (0.,0.,0.);
446  nodes[1] = Point (1.,0.,0.);
447  nodes[2] = Point (0.,1.,0.);
448  nodes[3] = Point (.5,0.,0.);
449  nodes[4] = Point (.5,.5,0.);
450  nodes[5] = Point (0.,.5,0.);
451  return;
452  }
453  case QUAD4:
454  {
455  nodes.resize(4);
456  nodes[0] = Point (-1.,-1.,0.);
457  nodes[1] = Point (1.,-1.,0.);
458  nodes[2] = Point (1.,1.,0.);
459  nodes[3] = Point (-1.,1.,0.);
460  return;
461  }
462  case QUAD8:
463  {
464  nodes.resize(8);
465  nodes[0] = Point (-1.,-1.,0.);
466  nodes[1] = Point (1.,-1.,0.);
467  nodes[2] = Point (1.,1.,0.);
468  nodes[3] = Point (-1.,1.,0.);
469  nodes[4] = Point (0.,-1.,0.);
470  nodes[5] = Point (1.,0.,0.);
471  nodes[6] = Point (0.,1.,0.);
472  nodes[7] = Point (-1.,0.,0.);
473  return;
474  }
475  case QUAD9:
476  {
477  nodes.resize(9);
478  nodes[0] = Point (-1.,-1.,0.);
479  nodes[1] = Point (1.,-1.,0.);
480  nodes[2] = Point (1.,1.,0.);
481  nodes[3] = Point (-1.,1.,0.);
482  nodes[4] = Point (0.,-1.,0.);
483  nodes[5] = Point (1.,0.,0.);
484  nodes[6] = Point (0.,1.,0.);
485  nodes[7] = Point (-1.,0.,0.);
486  nodes[8] = Point (0.,0.,0.);
487  return;
488  }
489  case TET4:
490  {
491  nodes.resize(4);
492  nodes[0] = Point (0.,0.,0.);
493  nodes[1] = Point (1.,0.,0.);
494  nodes[2] = Point (0.,1.,0.);
495  nodes[3] = Point (0.,0.,1.);
496  return;
497  }
498  case TET10:
499  {
500  nodes.resize(10);
501  nodes[0] = Point (0.,0.,0.);
502  nodes[1] = Point (1.,0.,0.);
503  nodes[2] = Point (0.,1.,0.);
504  nodes[3] = Point (0.,0.,1.);
505  nodes[4] = Point (.5,0.,0.);
506  nodes[5] = Point (.5,.5,0.);
507  nodes[6] = Point (0.,.5,0.);
508  nodes[7] = Point (0.,0.,.5);
509  nodes[8] = Point (.5,0.,.5);
510  nodes[9] = Point (0.,.5,.5);
511  return;
512  }
513  case HEX8:
514  {
515  nodes.resize(8);
516  nodes[0] = Point (-1.,-1.,-1.);
517  nodes[1] = Point (1.,-1.,-1.);
518  nodes[2] = Point (1.,1.,-1.);
519  nodes[3] = Point (-1.,1.,-1.);
520  nodes[4] = Point (-1.,-1.,1.);
521  nodes[5] = Point (1.,-1.,1.);
522  nodes[6] = Point (1.,1.,1.);
523  nodes[7] = Point (-1.,1.,1.);
524  return;
525  }
526  case HEX20:
527  {
528  nodes.resize(20);
529  nodes[0] = Point (-1.,-1.,-1.);
530  nodes[1] = Point (1.,-1.,-1.);
531  nodes[2] = Point (1.,1.,-1.);
532  nodes[3] = Point (-1.,1.,-1.);
533  nodes[4] = Point (-1.,-1.,1.);
534  nodes[5] = Point (1.,-1.,1.);
535  nodes[6] = Point (1.,1.,1.);
536  nodes[7] = Point (-1.,1.,1.);
537  nodes[8] = Point (0.,-1.,-1.);
538  nodes[9] = Point (1.,0.,-1.);
539  nodes[10] = Point (0.,1.,-1.);
540  nodes[11] = Point (-1.,0.,-1.);
541  nodes[12] = Point (-1.,-1.,0.);
542  nodes[13] = Point (1.,-1.,0.);
543  nodes[14] = Point (1.,1.,0.);
544  nodes[15] = Point (-1.,1.,0.);
545  nodes[16] = Point (0.,-1.,1.);
546  nodes[17] = Point (1.,0.,1.);
547  nodes[18] = Point (0.,1.,1.);
548  nodes[19] = Point (-1.,0.,1.);
549  return;
550  }
551  case HEX27:
552  {
553  nodes.resize(27);
554  nodes[0] = Point (-1.,-1.,-1.);
555  nodes[1] = Point (1.,-1.,-1.);
556  nodes[2] = Point (1.,1.,-1.);
557  nodes[3] = Point (-1.,1.,-1.);
558  nodes[4] = Point (-1.,-1.,1.);
559  nodes[5] = Point (1.,-1.,1.);
560  nodes[6] = Point (1.,1.,1.);
561  nodes[7] = Point (-1.,1.,1.);
562  nodes[8] = Point (0.,-1.,-1.);
563  nodes[9] = Point (1.,0.,-1.);
564  nodes[10] = Point (0.,1.,-1.);
565  nodes[11] = Point (-1.,0.,-1.);
566  nodes[12] = Point (-1.,-1.,0.);
567  nodes[13] = Point (1.,-1.,0.);
568  nodes[14] = Point (1.,1.,0.);
569  nodes[15] = Point (-1.,1.,0.);
570  nodes[16] = Point (0.,-1.,1.);
571  nodes[17] = Point (1.,0.,1.);
572  nodes[18] = Point (0.,1.,1.);
573  nodes[19] = Point (-1.,0.,1.);
574  nodes[20] = Point (0.,0.,-1.);
575  nodes[21] = Point (0.,-1.,0.);
576  nodes[22] = Point (1.,0.,0.);
577  nodes[23] = Point (0.,1.,0.);
578  nodes[24] = Point (-1.,0.,0.);
579  nodes[25] = Point (0.,0.,1.);
580  nodes[26] = Point (0.,0.,0.);
581  return;
582  }
583  case PRISM6:
584  {
585  nodes.resize(6);
586  nodes[0] = Point (0.,0.,-1.);
587  nodes[1] = Point (1.,0.,-1.);
588  nodes[2] = Point (0.,1.,-1.);
589  nodes[3] = Point (0.,0.,1.);
590  nodes[4] = Point (1.,0.,1.);
591  nodes[5] = Point (0.,1.,1.);
592  return;
593  }
594  case PRISM15:
595  {
596  nodes.resize(15);
597  nodes[0] = Point (0.,0.,-1.);
598  nodes[1] = Point (1.,0.,-1.);
599  nodes[2] = Point (0.,1.,-1.);
600  nodes[3] = Point (0.,0.,1.);
601  nodes[4] = Point (1.,0.,1.);
602  nodes[5] = Point (0.,1.,1.);
603  nodes[6] = Point (.5,0.,-1.);
604  nodes[7] = Point (.5,.5,-1.);
605  nodes[8] = Point (0.,.5,-1.);
606  nodes[9] = Point (0.,0.,0.);
607  nodes[10] = Point (1.,0.,0.);
608  nodes[11] = Point (0.,1.,0.);
609  nodes[12] = Point (.5,0.,1.);
610  nodes[13] = Point (.5,.5,1.);
611  nodes[14] = Point (0.,.5,1.);
612  return;
613  }
614  case PRISM18:
615  {
616  nodes.resize(18);
617  nodes[0] = Point (0.,0.,-1.);
618  nodes[1] = Point (1.,0.,-1.);
619  nodes[2] = Point (0.,1.,-1.);
620  nodes[3] = Point (0.,0.,1.);
621  nodes[4] = Point (1.,0.,1.);
622  nodes[5] = Point (0.,1.,1.);
623  nodes[6] = Point (.5,0.,-1.);
624  nodes[7] = Point (.5,.5,-1.);
625  nodes[8] = Point (0.,.5,-1.);
626  nodes[9] = Point (0.,0.,0.);
627  nodes[10] = Point (1.,0.,0.);
628  nodes[11] = Point (0.,1.,0.);
629  nodes[12] = Point (.5,0.,1.);
630  nodes[13] = Point (.5,.5,1.);
631  nodes[14] = Point (0.,.5,1.);
632  nodes[15] = Point (.5,0.,0.);
633  nodes[16] = Point (.5,.5,0.);
634  nodes[17] = Point (0.,.5,0.);
635  return;
636  }
637  case PYRAMID5:
638  {
639  nodes.resize(5);
640  nodes[0] = Point (-1.,-1.,0.);
641  nodes[1] = Point (1.,-1.,0.);
642  nodes[2] = Point (1.,1.,0.);
643  nodes[3] = Point (-1.,1.,0.);
644  nodes[4] = Point (0.,0.,1.);
645  return;
646  }
647  case PYRAMID14:
648  {
649  nodes.resize(14);
650 
651  // base corners
652  nodes[0] = Point (-1.,-1.,0.);
653  nodes[1] = Point (1.,-1.,0.);
654  nodes[2] = Point (1.,1.,0.);
655  nodes[3] = Point (-1.,1.,0.);
656 
657  // apex
658  nodes[4] = Point (0.,0.,1.);
659 
660  // base midedge
661  nodes[5] = Point (0.,-1.,0.);
662  nodes[6] = Point (1.,0.,0.);
663  nodes[7] = Point (0.,1.,0.);
664  nodes[8] = Point (-1,0.,0.);
665 
666  // lateral midedge
667  nodes[9] = Point (-.5,-.5,.5);
668  nodes[10] = Point (.5,-.5,.5);
669  nodes[11] = Point (.5,.5,.5);
670  nodes[12] = Point (-.5,.5,.5);
671 
672  // base center
673  nodes[13] = Point (0.,0.,0.);
674 
675  return;
676  }
677  default:
678  {
679  libMesh::err << "ERROR: Unknown element type " << itemType << std::endl;
680  libmesh_error();
681  }
682  }
683  return;
684 }
template<typename T>
const std::vector<RealGradient>& libMesh::FEGenericBase< T >::get_Sobolev_dweight ( ) const
inlineinherited
Returns
the first global derivative of the multiplicative weight at each quadrature point. See get_Sobolev_weight() for details. In case of FE initialized to all zero.

Definition at line 442 of file fe_base.h.

443  { return dweight; }
template<typename T>
const std::vector<Real>& libMesh::FEGenericBase< T >::get_Sobolev_weight ( ) const
inlineinherited
Returns
the multiplicative weight at each quadrature point. This weight is used for certain infinite element weak formulations, so that weighted Sobolev spaces are used for the trial function space. This renders the variational form easily computable.

In case of the general finite element class FE this field is initialized to all ones, so that the variational formulation for an infinite element returns correct element matrices for a mesh using both finite and infinite elements.

Definition at line 434 of file fe_base.h.

435  { return weight; }
const std::vector<std::vector<Point> >& libMesh::FEAbstract::get_tangents ( ) const
inlineinherited
Returns
the tangent vectors for face integration.

Definition at line 368 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

369  { return this->_fe_map->get_tangents(); }
ElemType libMesh::FEAbstract::get_type ( ) const
inlineinherited
Returns
the element type that the current shape functions have been calculated for. Useful in determining when shape functions must be recomputed.

Definition at line 407 of file fe_abstract.h.

References libMesh::FEAbstract::elem_type.

407 { return elem_type; }
const std::vector<Point>& libMesh::FEAbstract::get_xyz ( ) const
inlineinherited
Returns
the xyz spatial locations of the quadrature points on the element.

Definition at line 227 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

Referenced by libMesh::ExactErrorEstimator::find_squared_element_error(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), and libMesh::ProjectFEMSolution::operator()().

228  { return this->_fe_map->get_xyz(); }
void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 163 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

164 {
165  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
166  std::pair<unsigned int, unsigned int>& p = _counts[name];
167 
168  p.first++;
169 }
void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 176 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

177 {
178  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
179  std::pair<unsigned int, unsigned int>& p = _counts[name];
180 
181  p.second++;
182 }
template<unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
virtual void libMesh::InfFE< friend_Dim, friend_T_radial, friend_T_map >::init_base_shape_functions ( const std::vector< Point > &  ,
const Elem  
)
inlineprotectedvirtual

Do not use this derived member in InfFE<Dim,T_radial,T_map>.

Implements libMesh::FEGenericBase< T >.

Definition at line 520 of file inf_fe.h.

522  { libmesh_error(); }
template<unsigned int Dim, FEFamily T_radial, InfMapType T_base>
void libMesh::InfFE< Dim, T_radial, T_base >::init_face_shape_functions ( const std::vector< Point > &  qp,
const Elem side 
)
protected

Not implemented yet. Initialize all the data fields like weight, phi, etc for the side s.

Definition at line 142 of file inf_fe_boundary.C.

References libMesh::FEGenericBase< T >::build(), libMesh::libmesh_assert(), libMesh::Elem::p_level(), libMesh::AutoPtr< Tp >::release(), and libMesh::Elem::type().

144 {
145  libmesh_assert(inf_side);
146 
147  // Currently, this makes only sense in 3-D!
148  libmesh_assert_equal_to (Dim, 3);
149 
150  // Initialiize the radial shape functions
151  this->init_radial_shape_functions(inf_side);
152 
153  // Initialize the base shape functions
154  this->update_base_elem(inf_side);
155 
156  // Initialize the base quadratur rule
157  base_qrule->init(base_elem->type(), inf_side->p_level());
158 
159  // base_fe still corresponds to the (dim-1)-dimensional base of the InfFE object,
160  // so update the fe_base.
161  {
162  libmesh_assert_equal_to (Dim, 3);
163 
164  AutoPtr<FEBase> ap_fb(FEBase::build(Dim-2, this->fe_type));
165  if (base_fe != NULL)
166  delete base_fe;
167  base_fe = ap_fb.release();
169  }
170 
171  // initialize the shape functions on the base
173  base_elem);
174 
175  // the number of quadrature points
176  const unsigned int n_radial_qp =
177  libmesh_cast_int<unsigned int>(som.size());
178  const unsigned int n_base_qp = base_qrule->n_points();
179  const unsigned int n_total_qp = n_radial_qp * n_base_qp;
180 
181  // the quadratur weigths
182  _total_qrule_weights.resize(n_total_qp);
183 
184  // now inite the shapes for boundary work
185  {
186 
187  // The element type and order to use in the base map
188  const Order base_mapping_order ( base_elem->default_order() );
189  const ElemType base_mapping_elem_type ( base_elem->type() );
190 
191  // the number of mapping shape functions
192  // (Lagrange shape functions are used for mapping in the base)
193  const unsigned int n_radial_mapping_sf =
194  libmesh_cast_int<unsigned int>(radial_map.size());
195  const unsigned int n_base_mapping_shape_functions = Base::n_base_mapping_sf(base_mapping_elem_type,
196  base_mapping_order);
197 
198  const unsigned int n_total_mapping_shape_functions =
199  n_radial_mapping_sf * n_base_mapping_shape_functions;
200 
201 
202  // initialize the node and shape numbering maps
203  {
204  _radial_node_index.resize (n_total_mapping_shape_functions);
205  _base_node_index.resize (n_total_mapping_shape_functions);
206 
207  const ElemType inf_face_elem_type (inf_side->type());
208 
209  // fill the node index map
210  for (unsigned int n=0; n<n_total_mapping_shape_functions; n++)
211  {
212  compute_node_indices (inf_face_elem_type,
213  n,
214  _base_node_index[n],
215  _radial_node_index[n]);
216 
217  libmesh_assert_less (_base_node_index[n], n_base_mapping_shape_functions);
218  libmesh_assert_less (_radial_node_index[n], n_radial_mapping_sf);
219  }
220 
221  }
222 
223  // rezise map data fields
224  {
225  std::vector<std::vector<Real> >& psi_map = this->_fe_map->get_psi();
226  std::vector<std::vector<Real> >& dpsidxi_map = this->_fe_map->get_dpsidxi();
227  std::vector<std::vector<Real> >& d2psidxi2_map = this->_fe_map->get_d2psidxi2();
228  psi_map.resize (n_total_mapping_shape_functions);
229  dpsidxi_map.resize (n_total_mapping_shape_functions);
230  d2psidxi2_map.resize (n_total_mapping_shape_functions);
231 
232  // if (Dim == 3)
233  {
234  std::vector<std::vector<Real> >& dpsideta_map = this->_fe_map->get_dpsideta();
235  std::vector<std::vector<Real> >& d2psidxideta_map = this->_fe_map->get_d2psidxideta();
236  std::vector<std::vector<Real> >& d2psideta2_map = this->_fe_map->get_d2psideta2();
237  dpsideta_map.resize (n_total_mapping_shape_functions);
238  d2psidxideta_map.resize (n_total_mapping_shape_functions);
239  d2psideta2_map.resize (n_total_mapping_shape_functions);
240  }
241 
242  for (unsigned int i=0; i<n_total_mapping_shape_functions; i++)
243  {
244  psi_map[i].resize (n_total_qp);
245  dpsidxi_map[i].resize (n_total_qp);
246  d2psidxi2_map[i].resize (n_total_qp);
247 
248  // if (Dim == 3)
249  {
250  std::vector<std::vector<Real> >& dpsideta_map = this->_fe_map->get_dpsideta();
251  std::vector<std::vector<Real> >& d2psidxideta_map = this->_fe_map->get_d2psidxideta();
252  std::vector<std::vector<Real> >& d2psideta2_map = this->_fe_map->get_d2psideta2();
253  dpsideta_map[i].resize (n_total_qp);
254  d2psidxideta_map[i].resize (n_total_qp);
255  d2psideta2_map[i].resize (n_total_qp);
256  }
257  }
258  }
259 
260 
261  // compute shape maps
262  {
263  const std::vector<std::vector<Real> >& S_map = (base_fe->get_fe_map()).get_phi_map();
264  const std::vector<std::vector<Real> >& Ss_map = (base_fe->get_fe_map()).get_dphidxi_map();
265 
266  std::vector<std::vector<Real> >& psi_map = this->_fe_map->get_psi();
267  std::vector<std::vector<Real> >& dpsidxi_map = this->_fe_map->get_dpsidxi();
268  std::vector<std::vector<Real> >& dpsideta_map = this->_fe_map->get_dpsideta();
269 
270  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qp's
271  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qp's
272  for (unsigned int ti=0; ti<n_total_mapping_shape_functions; ti++) // over all mapping shapes
273  {
274  // let the index vectors take care of selecting the appropriate base/radial mapping shape
275  const unsigned int bi =