libMesh::FEGenericBase< OutputType > Class Template Referenceabstract

#include <exact_error_estimator.h>

Inheritance diagram for libMesh::FEGenericBase< OutputType >:

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

virtual ~FEGenericBase ()
 
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
 
template<>
UniquePtr< FEGenericBase< Real > > build (const unsigned int dim, const FEType &fet)
 
template<>
UniquePtr< FEGenericBase< RealGradient > > build (const unsigned int dim, const FEType &fet)
 
template<>
UniquePtr< FEGenericBase< Real > > build_InfFE (const unsigned int dim, const FEType &fet)
 
template<>
UniquePtr< FEGenericBase< RealGradient > > build_InfFE (const unsigned int, const FEType &)
 
virtual void reinit (const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr)=0
 
virtual void reinit (const Elem *elem, const unsigned int side, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr)=0
 
virtual void edge_reinit (const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *pts=libmesh_nullptr, const std::vector< Real > *weights=libmesh_nullptr)=0
 
virtual void side_map (const Elem *elem, const Elem *side, const unsigned int s, const std::vector< Point > &reference_side_points, std::vector< Point > &reference_points)=0
 
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
 
virtual void attach_quadrature_rule (QBase *q)=0
 
virtual unsigned int n_shape_functions () const =0
 
virtual unsigned int n_quadrature_points () const =0
 
ElemType get_type () const
 
unsigned int get_p_level () const
 
FEType get_fe_type () const
 
Order get_order () const
 
void set_fe_order (int new_order)
 
virtual FEContinuity get_continuity () const =0
 
virtual bool is_hierarchic () const =0
 
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 UniquePtr< FEGenericBasebuild (const unsigned int dim, const FEType &type)
 
static UniquePtr< FEGenericBasebuild_InfFE (const unsigned int dim, const FEType &type)
 
static void compute_proj_constraints (DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
 
static void coarsened_dof_values (const NumericVector< Number > &global_vector, const DofMap &dof_map, const Elem *coarse_elem, DenseVector< Number > &coarse_dofs, const unsigned int var, const bool use_old_dof_indices=false)
 
static void coarsened_dof_values (const NumericVector< Number > &global_vector, const DofMap &dof_map, const Elem *coarse_elem, DenseVector< Number > &coarse_dofs, const bool use_old_dof_indices=false)
 
static void compute_periodic_constraints (DofConstraints &constraints, DofMap &dof_map, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const unsigned int variable_number, const Elem *elem)
 
static bool on_reference_element (const Point &p, const ElemType t, const Real eps=TOLERANCE)
 
static void get_refspace_nodes (const ElemType t, std::vector< Point > &nodes)
 
static void compute_node_constraints (NodeConstraints &constraints, const Elem *elem)
 
static void compute_periodic_node_constraints (NodeConstraints &constraints, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const Elem *elem)
 
static void print_info (std::ostream &out=libMesh::out)
 
static std::string get_info ()
 
static unsigned int n_objects ()
 
static void enable_print_counter_info ()
 
static void disable_print_counter_info ()
 

Protected Types

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

Protected Member Functions

 FEGenericBase (const unsigned int dim, const FEType &fet)
 
virtual void init_base_shape_functions (const std::vector< Point > &qp, const Elem *e)=0
 
void determine_calculations ()
 
virtual void compute_shape_functions (const Elem *elem, const std::vector< Point > &qp)
 
virtual bool shapes_need_reinit () const =0
 
void increment_constructor_count (const std::string &name)
 
void increment_destructor_count (const std::string &name)
 

Protected Attributes

UniquePtr< FETransformationBase< OutputType > > _fe_trans
 
std::vector< std::vector< OutputShape > > phi
 
std::vector< std::vector< OutputGradient > > dphi
 
std::vector< std::vector< OutputShape > > curl_phi
 
std::vector< std::vector< OutputDivergence > > div_phi
 
std::vector< std::vector< OutputShape > > dphidxi
 
std::vector< std::vector< OutputShape > > dphideta
 
std::vector< std::vector< OutputShape > > dphidzeta
 
std::vector< std::vector< OutputShape > > dphidx
 
std::vector< std::vector< OutputShape > > dphidy
 
std::vector< std::vector< OutputShape > > dphidz
 
std::vector< std::vector< OutputTensor > > d2phi
 
std::vector< std::vector< OutputShape > > d2phidxi2
 
std::vector< std::vector< OutputShape > > d2phidxideta
 
std::vector< std::vector< OutputShape > > d2phidxidzeta
 
std::vector< std::vector< OutputShape > > d2phideta2
 
std::vector< std::vector< OutputShape > > d2phidetadzeta
 
std::vector< std::vector< OutputShape > > d2phidzeta2
 
std::vector< std::vector< OutputShape > > d2phidx2
 
std::vector< std::vector< OutputShape > > d2phidxdy
 
std::vector< std::vector< OutputShape > > d2phidxdz
 
std::vector< std::vector< OutputShape > > d2phidy2
 
std::vector< std::vector< OutputShape > > d2phidydz
 
std::vector< std::vector< OutputShape > > d2phidz2
 
std::vector< OutputGradientdphase
 
std::vector< RealGradientdweight
 
std::vector< Realweight
 
UniquePtr< FEMap_fe_map
 
const unsigned int dim
 
bool calculations_started
 
bool calculate_phi
 
bool calculate_dphi
 
bool calculate_d2phi
 
bool calculate_curl_phi
 
bool calculate_div_phi
 
bool calculate_dphiref
 
FEType fe_type
 
ElemType elem_type
 
unsigned int _p_level
 
QBaseqrule
 
bool shapes_on_quadrature
 

Static Protected Attributes

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

Friends

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

Detailed Description

template<typename OutputType>
class libMesh::FEGenericBase< OutputType >

This class forms the foundation from which generic finite elements may be derived. In the current implementation the templated derived class FE offers a wide variety of commonly used finite element concepts. Check there for details.

Use the FEGenericBase<OutputType>::build() method to create an object of any of the derived classes which is compatible with OutputType.

Author
Benjamin S. Kirk
Date
2002

Definition at line 37 of file exact_error_estimator.h.

Member Typedef Documentation

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

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

Definition at line 119 of file reference_counter.h.

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

Definition at line 123 of file fe_base.h.

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

Definition at line 121 of file fe_base.h.

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

Definition at line 124 of file fe_base.h.

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

Definition at line 127 of file fe_base.h.

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

Definition at line 125 of file fe_base.h.

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

Definition at line 126 of file fe_base.h.

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

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

Definition at line 120 of file fe_base.h.

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

Definition at line 122 of file fe_base.h.

Constructor & Destructor Documentation

template<typename OutputType >
libMesh::FEGenericBase< OutputType >::FEGenericBase ( const unsigned int  dim,
const FEType fet 
)
inlineprotected

Constructor. Optionally initializes required data structures. Protected so that this base class cannot be explicitly instantiated.

Definition at line 677 of file fe_base.h.

References libMesh::FEGenericBase< OutputType >::~FEGenericBase().

678  :
679  FEAbstract(d,fet),
681  phi(),
682  dphi(),
683  curl_phi(),
684  div_phi(),
685  dphidxi(),
686  dphideta(),
687  dphidzeta(),
688  dphidx(),
689  dphidy(),
690  dphidz()
691 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
692  ,d2phi(),
693  d2phidxi2(),
694  d2phidxideta(),
695  d2phidxidzeta(),
696  d2phideta2(),
697  d2phidetadzeta(),
698  d2phidzeta2(),
699  d2phidx2(),
700  d2phidxdy(),
701  d2phidxdz(),
702  d2phidy2(),
703  d2phidydz(),
704  d2phidz2()
705 #endif
706 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
707  ,dphase(),
708  dweight(),
709  weight()
710 #endif
711 {
712 }
713 
714 
715 
716 template <typename OutputType>
717 inline
719 {
720 }
std::vector< std::vector< OutputTensor > > d2phi
Definition: fe_base.h:552
std::vector< std::vector< OutputShape > > dphidxi
Definition: fe_base.h:519
std::vector< std::vector< OutputShape > > d2phidxdz
Definition: fe_base.h:597
std::vector< std::vector< OutputShape > > dphidzeta
Definition: fe_base.h:529
std::vector< std::vector< OutputShape > > d2phidydz
Definition: fe_base.h:607
static UniquePtr< FETransformationBase< OutputShape > > build(const FEType &type)
std::vector< std::vector< OutputShape > > d2phidxideta
Definition: fe_base.h:562
FEAbstract(const unsigned int dim, const FEType &fet)
Definition: fe_abstract.h:602
std::vector< Real > weight
Definition: fe_base.h:644
std::vector< std::vector< OutputShape > > d2phidx2
Definition: fe_base.h:587
std::vector< std::vector< OutputShape > > curl_phi
Definition: fe_base.h:509
std::vector< std::vector< OutputShape > > dphidy
Definition: fe_base.h:539
std::vector< std::vector< OutputShape > > d2phidy2
Definition: fe_base.h:602
std::vector< std::vector< OutputShape > > d2phidetadzeta
Definition: fe_base.h:577
std::vector< std::vector< OutputShape > > d2phidxidzeta
Definition: fe_base.h:567
std::vector< std::vector< OutputShape > > d2phidxdy
Definition: fe_base.h:592
std::vector< std::vector< OutputShape > > dphidx
Definition: fe_base.h:534
std::vector< std::vector< OutputShape > > phi
Definition: fe_base.h:499
std::vector< std::vector< OutputShape > > d2phideta2
Definition: fe_base.h:572
std::vector< OutputGradient > dphase
Definition: fe_base.h:630
std::vector< std::vector< OutputDivergence > > div_phi
Definition: fe_base.h:514
std::vector< std::vector< OutputGradient > > dphi
Definition: fe_base.h:504
UniquePtr< FETransformationBase< OutputType > > _fe_trans
Definition: fe_base.h:494
std::vector< std::vector< OutputShape > > d2phidz2
Definition: fe_base.h:612
std::vector< std::vector< OutputShape > > d2phidxi2
Definition: fe_base.h:557
std::vector< std::vector< OutputShape > > d2phidzeta2
Definition: fe_base.h:582
std::vector< std::vector< OutputShape > > dphidz
Definition: fe_base.h:544
std::vector< std::vector< OutputShape > > dphideta
Definition: fe_base.h:524
std::vector< RealGradient > dweight
Definition: fe_base.h:637
template<typename OutputType>
virtual libMesh::FEGenericBase< OutputType >::~FEGenericBase ( )
virtual

Member Function Documentation

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

Definition at line 186 of file fe_base.C.

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

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

Definition at line 387 of file fe_base.C.

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

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

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

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

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

Definition at line 463 of file fe_base.C.

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

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

Definition at line 670 of file fe_base.C.

672 {
673  // No vector types defined... YET.
674  libmesh_not_implemented();
675  return UniquePtr<FEVectorBase>();
676 }
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::coarsened_dof_values ( const NumericVector< Number > &  global_vector,
const DofMap dof_map,
const Elem coarse_elem,
DenseVector< Number > &  coarse_dofs,
const unsigned int  var,
const bool  use_old_dof_indices = false 
)
static

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

Definition at line 804 of file fe_base.C.

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

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

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

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

Definition at line 1351 of file fe_base.C.

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

1356 {
1357  Ue.resize(0);
1358 
1359  for (unsigned int v=0; v != dof_map.n_variables(); ++v)
1360  {
1361  DenseVector<Number> Usub;
1362 
1363  coarsened_dof_values(old_vector, dof_map, elem, Usub,
1364  use_old_dof_indices);
1365 
1366  Ue.append (Usub);
1367  }
1368 }
static void coarsened_dof_values(const NumericVector< Number > &global_vector, const DofMap &dof_map, const Elem *coarse_elem, DenseVector< Number > &coarse_dofs, const unsigned int var, const bool use_old_dof_indices=false)
Definition: fe_base.C:804
unsigned int n_variables() const
Definition: dof_map.h:477
void libMesh::FEAbstract::compute_node_constraints ( NodeConstraints constraints,
const Elem elem 
)
staticinherited

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

Definition at line 791 of file fe_abstract.C.

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

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

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

Definition at line 1658 of file fe_base.C.

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

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

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

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

Definition at line 934 of file fe_abstract.C.

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

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

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

Definition at line 1374 of file fe_base.C.

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

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

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

Implements libMesh::FEAbstract.

Reimplemented in libMesh::FEXYZ< Dim >, and libMesh::InfFE< Dim, T_radial, T_map >.

Definition at line 682 of file fe_base.C.

Referenced by libMesh::FEGenericBase< FEOutputType< T >::type >::get_Sobolev_dweight().

684 {
685  //-------------------------------------------------------------------------
686  // Compute the shape function values (and derivatives)
687  // at the Quadrature points. Note that the actual values
688  // have already been computed via init_shape_functions
689 
690  // Start logging the shape function computation
691  LOG_SCOPE("compute_shape_functions()", "FE");
692 
693  this->determine_calculations();
694 
695  if (calculate_phi)
696  this->_fe_trans->map_phi(this->dim, elem, qp, (*this), this->phi);
697 
698  if (calculate_dphi)
699  this->_fe_trans->map_dphi(this->dim, elem, qp, (*this), this->dphi,
700  this->dphidx, this->dphidy, this->dphidz);
701 
702 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
703  if (calculate_d2phi)
704  this->_fe_trans->map_d2phi(this->dim, qp, (*this), this->d2phi,
705  this->d2phidx2, this->d2phidxdy, this->d2phidxdz,
706  this->d2phidy2, this->d2phidydz, this->d2phidz2);
707 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
708 
709  // Only compute curl for vector-valued elements
711  this->_fe_trans->map_curl(this->dim, elem, qp, (*this), this->curl_phi);
712 
713  // Only compute div for vector-valued elements
715  this->_fe_trans->map_div(this->dim, elem, qp, (*this), this->div_phi);
716 }
std::vector< std::vector< OutputTensor > > d2phi
Definition: fe_base.h:552
std::vector< std::vector< OutputShape > > d2phidxdz
Definition: fe_base.h:597
std::vector< std::vector< OutputShape > > d2phidydz
Definition: fe_base.h:607
std::vector< std::vector< OutputShape > > d2phidx2
Definition: fe_base.h:587
std::vector< std::vector< OutputShape > > curl_phi
Definition: fe_base.h:509
std::vector< std::vector< OutputShape > > dphidy
Definition: fe_base.h:539
std::vector< std::vector< OutputShape > > d2phidy2
Definition: fe_base.h:602
std::vector< std::vector< OutputShape > > d2phidxdy
Definition: fe_base.h:592
std::vector< std::vector< OutputShape > > dphidx
Definition: fe_base.h:534
std::vector< std::vector< OutputShape > > phi
Definition: fe_base.h:499
const unsigned int dim
Definition: fe_abstract.h:517
std::vector< std::vector< OutputDivergence > > div_phi
Definition: fe_base.h:514
void determine_calculations()
Definition: fe_base.C:741
std::vector< std::vector< OutputGradient > > dphi
Definition: fe_base.h:504
UniquePtr< FETransformationBase< OutputType > > _fe_trans
Definition: fe_base.h:494
std::vector< std::vector< OutputShape > > d2phidz2
Definition: fe_base.h:612
std::vector< std::vector< OutputShape > > dphidz
Definition: fe_base.h:544
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::determine_calculations ( )
protected

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

Definition at line 741 of file fe_base.C.

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

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

742 {
743  this->calculations_started = true;
744 
745  // If the user forgot to request anything, we'll be safe and
746  // calculate everything:
747 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
748  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_d2phi
749  && !this->calculate_curl_phi && !this->calculate_div_phi)
750  {
751  this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = this->calculate_dphiref = true;
753  {
754  this->calculate_curl_phi = true;
755  this->calculate_div_phi = true;
756  }
757  }
758 #else
759  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_curl_phi && !this->calculate_div_phi)
760  {
761  this->calculate_phi = this->calculate_dphi = this->calculate_dphiref = true;
763  {
764  this->calculate_curl_phi = true;
765  this->calculate_div_phi = true;
766  }
767  }
768 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
769 
770  // Request whichever terms are necessary from the FEMap
771  if (this->calculate_phi)
772  this->_fe_trans->init_map_phi(*this);
773 
774  if (this->calculate_dphiref)
775  this->_fe_trans->init_map_dphi(*this);
776 
777 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
778  if (this->calculate_d2phi)
779  this->_fe_trans->init_map_d2phi(*this);
780 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
781 }
FEFamily family
Definition: fe_type.h:203
static FEFieldType field_type(const FEType &fe_type)
UniquePtr< FETransformationBase< OutputType > > _fe_trans
Definition: fe_base.h:494
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited
virtual void libMesh::FEAbstract::edge_reinit ( const Elem elem,
const unsigned int  edge,
const Real  tolerance = TOLERANCE,
const std::vector< Point > *  pts = libmesh_nullptr,
const std::vector< Real > *  weights = libmesh_nullptr 
)
pure virtualinherited

Reinitializes all the physical element-dependent data based on the edge of the element elem. The tolerance paremeter is passed to the involved call to inverse_map(). By default the element data are computed at the quadrature points specified by the quadrature rule qrule, but any set of points on the reference edge element may be specified in the optional argument pts.

Implemented in libMesh::InfFE< Dim, T_radial, T_map >, libMesh::FE< Dim, T >, libMesh::FE< 2, SUBDIVISION >, libMesh::FE< Dim, HIERARCHIC >, libMesh::FE< Dim, SCALAR >, libMesh::FE< Dim, L2_LAGRANGE >, libMesh::FE< Dim, NEDELEC_ONE >, libMesh::FE< Dim, HERMITE >, libMesh::FE< Dim, CLOUGH >, libMesh::FE< Dim, MONOMIAL >, libMesh::FE< Dim, XYZ >, libMesh::FE< Dim, LAGRANGE >, libMesh::FE< Dim, L2_HIERARCHIC >, and libMesh::FE< Dim, LAGRANGE_VEC >.

void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

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

Definition at line 101 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

102 {
103  _enable_print_counter = true;
104  return;
105 }
virtual FEContinuity libMesh::FEAbstract::get_continuity ( ) const
pure virtualinherited
Returns
The continuity level of the finite element.

Implemented in libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::InfFE< Dim, T_radial, T_map >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< 2, SUBDIVISION >, libMesh::FE< Dim, HIERARCHIC >, libMesh::FE< Dim, SCALAR >, libMesh::FE< Dim, L2_LAGRANGE >, libMesh::FE< Dim, NEDELEC_ONE >, libMesh::FE< Dim, HERMITE >, libMesh::FE< Dim, CLOUGH >, libMesh::FE< Dim, MONOMIAL >, libMesh::FE< Dim, XYZ >, libMesh::FE< Dim, LAGRANGE >, libMesh::FE< Dim, L2_HIERARCHIC >, libMesh::FE< Dim, LAGRANGE_VEC >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, and libMesh::FE< Dim, T >.

Referenced by libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator()(), and libMesh::FEAbstract::set_fe_order().

template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_curl_phi ( ) const
inline
Returns
The curl of the shape function at the quadrature points.

Definition at line 224 of file fe_base.h.

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

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

Definition at line 377 of file fe_abstract.h.

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

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

Definition at line 370 of file fe_base.h.

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

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

Definition at line 378 of file fe_base.h.

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

libmesh_assert(j)
std::vector< std::vector< OutputShape > > d2phidetadzeta
Definition: fe_base.h:577
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidx2 ( ) const
inline
Returns
The shape function second derivatives at the quadrature points.

Definition at line 298 of file fe_base.h.

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

Definition at line 306 of file fe_base.h.

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

Definition at line 314 of file fe_base.h.

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

Definition at line 346 of file fe_base.h.

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

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

Definition at line 354 of file fe_base.h.

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

std::vector< std::vector< OutputShape > > d2phidxideta
Definition: fe_base.h:562
libmesh_assert(j)
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidxidzeta ( ) const
inline
Returns
The shape function second derivatives at the quadrature points, in reference coordinates

Definition at line 362 of file fe_base.h.

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

libmesh_assert(j)
std::vector< std::vector< OutputShape > > d2phidxidzeta
Definition: fe_base.h:567
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_d2phidy2 ( ) const
inline
Returns
The shape function second derivatives at the quadrature points.

Definition at line 322 of file fe_base.h.

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

Definition at line 330 of file fe_base.h.

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

Definition at line 338 of file fe_base.h.

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

Definition at line 386 of file fe_base.h.

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

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

Definition at line 264 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

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

Definition at line 294 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

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

Definition at line 258 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

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

Definition at line 280 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

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

Definition at line 288 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

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

Definition at line 272 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

273  { return this->_fe_map->get_d2xyzdzeta2(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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 324 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

325  { return this->_fe_map->get_detadx(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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 331 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

332  { return this->_fe_map->get_detady(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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 338 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

339  { return this->_fe_map->get_detadz(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
template<typename OutputType>
const std::vector<std::vector<OutputDivergence> >& libMesh::FEGenericBase< OutputType >::get_div_phi ( ) const
inline
Returns
The divergence of the shape function at the quadrature points.

Definition at line 232 of file fe_base.h.

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

234  calculate_div_phi = calculate_dphiref = true; return div_phi; }
libmesh_assert(j)
std::vector< std::vector< OutputDivergence > > div_phi
Definition: fe_base.h:514
template<typename OutputType>
const std::vector<OutputGradient>& libMesh::FEGenericBase< OutputType >::get_dphase ( ) const
inline
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 produces correct element matrices for a mesh using both finite and infinite elements.

Definition at line 404 of file fe_base.h.

405  { return dphase; }
std::vector< OutputGradient > dphase
Definition: fe_base.h:630
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_dphideta ( ) const
inline
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_dphidx ( ) const
inline
Returns
The shape function x-derivative at the quadrature points.

Definition at line 240 of file fe_base.h.

242  calculate_dphi = calculate_dphiref = true; return dphidx; }
libmesh_assert(j)
std::vector< std::vector< OutputShape > > dphidx
Definition: fe_base.h:534
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_dphidxi ( ) const
inline
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_dphidy ( ) const
inline
Returns
The shape function y-derivative at the quadrature points.

Definition at line 248 of file fe_base.h.

250  calculate_dphi = calculate_dphiref = true; return dphidy; }
std::vector< std::vector< OutputShape > > dphidy
Definition: fe_base.h:539
libmesh_assert(j)
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_dphidz ( ) const
inline
Returns
The shape function z-derivative at the quadrature points.

Definition at line 256 of file fe_base.h.

258  calculate_dphi = calculate_dphiref = true; return dphidz; }
libmesh_assert(j)
std::vector< std::vector< OutputShape > > dphidz
Definition: fe_base.h:544
template<typename OutputType>
const std::vector<std::vector<OutputShape> >& libMesh::FEGenericBase< OutputType >::get_dphidzeta ( ) const
inline
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 303 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

304  { return this->_fe_map->get_dxidx(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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 310 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

311  { return this->_fe_map->get_dxidy(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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 317 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

318  { return this->_fe_map->get_dxidz(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
const std::vector<RealGradient>& libMesh::FEAbstract::get_dxyzdeta ( ) const
inlineinherited
Returns
The element tangents in eta-direction at the quadrature points.

Definition at line 245 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

246  { return this->_fe_map->get_dxyzdeta(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
const std::vector<RealGradient>& libMesh::FEAbstract::get_dxyzdxi ( ) const
inlineinherited
Returns
The element tangents in xi-direction at the quadrature points.

Definition at line 238 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

239  { return this->_fe_map->get_dxyzdxi(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
const std::vector<RealGradient>& libMesh::FEAbstract::get_dxyzdzeta ( ) const
inlineinherited
Returns
The element tangents in zeta-direction at the quadrature points.

Definition at line 252 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

253  { return _fe_map->get_dxyzdzeta(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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 345 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

346  { return this->_fe_map->get_dzetadx(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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 352 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

353  { return this->_fe_map->get_dzetady(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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 359 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

360  { return this->_fe_map->get_dzetadz(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
FEFamily libMesh::FEAbstract::get_family ( ) const
inlineinherited
Returns
The finite element family of this element.

Definition at line 441 of file fe_abstract.h.

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

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

441 { return fe_type.family; }
FEFamily family
Definition: fe_type.h:203
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 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
const std::vector<Point>& libMesh::FEAbstract::get_normals ( ) const
inlineinherited
Order libMesh::FEAbstract::get_order ( ) const
inlineinherited
Returns
The approximation order of the finite element.

Definition at line 420 of file fe_abstract.h.

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

420 { return static_cast<Order>(fe_type.order + _p_level); }
unsigned int _p_level
Definition: fe_abstract.h:573
OrderWrapper order
Definition: fe_type.h:197
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 410 of file fe_abstract.h.

References libMesh::FEAbstract::_p_level.

410 { return _p_level; }
unsigned int _p_level
Definition: fe_abstract.h:573
void libMesh::FEAbstract::get_refspace_nodes ( const ElemType  t,
std::vector< Point > &  nodes 
)
staticinherited
Returns
The reference space coordinates of nodes based on the element type.

Definition at line 259 of file fe_abstract.C.

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

260 {
261  switch(itemType)
262  {
263  case EDGE2:
264  {
265  nodes.resize(2);
266  nodes[0] = Point (-1.,0.,0.);
267  nodes[1] = Point (1.,0.,0.);
268  return;
269  }
270  case EDGE3:
271  {
272  nodes.resize(3);
273  nodes[0] = Point (-1.,0.,0.);
274  nodes[1] = Point (1.,0.,0.);
275  nodes[2] = Point (0.,0.,0.);
276  return;
277  }
278  case TRI3:
279  case TRISHELL3:
280  {
281  nodes.resize(3);
282  nodes[0] = Point (0.,0.,0.);
283  nodes[1] = Point (1.,0.,0.);
284  nodes[2] = Point (0.,1.,0.);
285  return;
286  }
287  case TRI6:
288  {
289  nodes.resize(6);
290  nodes[0] = Point (0.,0.,0.);
291  nodes[1] = Point (1.,0.,0.);
292  nodes[2] = Point (0.,1.,0.);
293  nodes[3] = Point (.5,0.,0.);
294  nodes[4] = Point (.5,.5,0.);
295  nodes[5] = Point (0.,.5,0.);
296  return;
297  }
298  case QUAD4:
299  case QUADSHELL4:
300  {
301  nodes.resize(4);
302  nodes[0] = Point (-1.,-1.,0.);
303  nodes[1] = Point (1.,-1.,0.);
304  nodes[2] = Point (1.,1.,0.);
305  nodes[3] = Point (-1.,1.,0.);
306  return;
307  }
308  case QUAD8:
309  {
310  nodes.resize(8);
311  nodes[0] = Point (-1.,-1.,0.);
312  nodes[1] = Point (1.,-1.,0.);
313  nodes[2] = Point (1.,1.,0.);
314  nodes[3] = Point (-1.,1.,0.);
315  nodes[4] = Point (0.,-1.,0.);
316  nodes[5] = Point (1.,0.,0.);
317  nodes[6] = Point (0.,1.,0.);
318  nodes[7] = Point (-1.,0.,0.);
319  return;
320  }
321  case QUAD9:
322  {
323  nodes.resize(9);
324  nodes[0] = Point (-1.,-1.,0.);
325  nodes[1] = Point (1.,-1.,0.);
326  nodes[2] = Point (1.,1.,0.);
327  nodes[3] = Point (-1.,1.,0.);
328  nodes[4] = Point (0.,-1.,0.);
329  nodes[5] = Point (1.,0.,0.);
330  nodes[6] = Point (0.,1.,0.);
331  nodes[7] = Point (-1.,0.,0.);
332  nodes[8] = Point (0.,0.,0.);
333  return;
334  }
335  case TET4:
336  {
337  nodes.resize(4);
338  nodes[0] = Point (0.,0.,0.);
339  nodes[1] = Point (1.,0.,0.);
340  nodes[2] = Point (0.,1.,0.);
341  nodes[3] = Point (0.,0.,1.);
342  return;
343  }
344  case TET10:
345  {
346  nodes.resize(10);
347  nodes[0] = Point (0.,0.,0.);
348  nodes[1] = Point (1.,0.,0.);
349  nodes[2] = Point (0.,1.,0.);
350  nodes[3] = Point (0.,0.,1.);
351  nodes[4] = Point (.5,0.,0.);
352  nodes[5] = Point (.5,.5,0.);
353  nodes[6] = Point (0.,.5,0.);
354  nodes[7] = Point (0.,0.,.5);
355  nodes[8] = Point (.5,0.,.5);
356  nodes[9] = Point (0.,.5,.5);
357  return;
358  }
359  case HEX8:
360  {
361  nodes.resize(8);
362  nodes[0] = Point (-1.,-1.,-1.);
363  nodes[1] = Point (1.,-1.,-1.);
364  nodes[2] = Point (1.,1.,-1.);
365  nodes[3] = Point (-1.,1.,-1.);
366  nodes[4] = Point (-1.,-1.,1.);
367  nodes[5] = Point (1.,-1.,1.);
368  nodes[6] = Point (1.,1.,1.);
369  nodes[7] = Point (-1.,1.,1.);
370  return;
371  }
372  case HEX20:
373  {
374  nodes.resize(20);
375  nodes[0] = Point (-1.,-1.,-1.);
376  nodes[1] = Point (1.,-1.,-1.);
377  nodes[2] = Point (1.,1.,-1.);
378  nodes[3] = Point (-1.,1.,-1.);
379  nodes[4] = Point (-1.,-1.,1.);
380  nodes[5] = Point (1.,-1.,1.);
381  nodes[6] = Point (1.,1.,1.);
382  nodes[7] = Point (-1.,1.,1.);
383  nodes[8] = Point (0.,-1.,-1.);
384  nodes[9] = Point (1.,0.,-1.);
385  nodes[10] = Point (0.,1.,-1.);
386  nodes[11] = Point (-1.,0.,-1.);
387  nodes[12] = Point (-1.,-1.,0.);
388  nodes[13] = Point (1.,-1.,0.);
389  nodes[14] = Point (1.,1.,0.);
390  nodes[15] = Point (-1.,1.,0.);
391  nodes[16] = Point (0.,-1.,1.);
392  nodes[17] = Point (1.,0.,1.);
393  nodes[18] = Point (0.,1.,1.);
394  nodes[19] = Point (-1.,0.,1.);
395  return;
396  }
397  case HEX27:
398  {
399  nodes.resize(27);
400  nodes[0] = Point (-1.,-1.,-1.);
401  nodes[1] = Point (1.,-1.,-1.);
402  nodes[2] = Point (1.,1.,-1.);
403  nodes[3] = Point (-1.,1.,-1.);
404  nodes[4] = Point (-1.,-1.,1.);
405  nodes[5] = Point (1.,-1.,1.);
406  nodes[6] = Point (1.,1.,1.);
407  nodes[7] = Point (-1.,1.,1.);
408  nodes[8] = Point (0.,-1.,-1.);
409  nodes[9] = Point (1.,0.,-1.);
410  nodes[10] = Point (0.,1.,-1.);
411  nodes[11] = Point (-1.,0.,-1.);
412  nodes[12] = Point (-1.,-1.,0.);
413  nodes[13] = Point (1.,-1.,0.);
414  nodes[14] = Point (1.,1.,0.);
415  nodes[15] = Point (-1.,1.,0.);
416  nodes[16] = Point (0.,-1.,1.);
417  nodes[17] = Point (1.,0.,1.);
418  nodes[18] = Point (0.,1.,1.);
419  nodes[19] = Point (-1.,0.,1.);
420  nodes[20] = Point (0.,0.,-1.);
421  nodes[21] = Point (0.,-1.,0.);
422  nodes[22] = Point (1.,0.,0.);
423  nodes[23] = Point (0.,1.,0.);
424  nodes[24] = Point (-1.,0.,0.);
425  nodes[25] = Point (0.,0.,1.);
426  nodes[26] = Point (0.,0.,0.);
427  return;
428  }
429  case PRISM6:
430  {
431  nodes.resize(6);
432  nodes[0] = Point (0.,0.,-1.);
433  nodes[1] = Point (1.,0.,-1.);
434  nodes[2] = Point (0.,1.,-1.);
435  nodes[3] = Point (0.,0.,1.);
436  nodes[4] = Point (1.,0.,1.);
437  nodes[5] = Point (0.,1.,1.);
438  return;
439  }
440  case PRISM15:
441  {
442  nodes.resize(15);
443  nodes[0] = Point (0.,0.,-1.);
444  nodes[1] = Point (1.,0.,-1.);
445  nodes[2] = Point (0.,1.,-1.);
446  nodes[3] = Point (0.,0.,1.);
447  nodes[4] = Point (1.,0.,1.);
448  nodes[5] = Point (0.,1.,1.);
449  nodes[6] = Point (.5,0.,-1.);
450  nodes[7] = Point (.5,.5,-1.);
451  nodes[8] = Point (0.,.5,-1.);
452  nodes[9] = Point (0.,0.,0.);
453  nodes[10] = Point (1.,0.,0.);
454  nodes[11] = Point (0.,1.,0.);
455  nodes[12] = Point (.5,0.,1.);
456  nodes[13] = Point (.5,.5,1.);
457  nodes[14] = Point (0.,.5,1.);
458  return;
459  }
460  case PRISM18:
461  {
462  nodes.resize(18);
463  nodes[0] = Point (0.,0.,-1.);
464  nodes[1] = Point (1.,0.,-1.);
465  nodes[2] = Point (0.,1.,-1.);
466  nodes[3] = Point (0.,0.,1.);
467  nodes[4] = Point (1.,0.,1.);
468  nodes[5] = Point (0.,1.,1.);
469  nodes[6] = Point (.5,0.,-1.);
470  nodes[7] = Point (.5,.5,-1.);
471  nodes[8] = Point (0.,.5,-1.);
472  nodes[9] = Point (0.,0.,0.);
473  nodes[10] = Point (1.,0.,0.);
474  nodes[11] = Point (0.,1.,0.);
475  nodes[12] = Point (.5,0.,1.);
476  nodes[13] = Point (.5,.5,1.);
477  nodes[14] = Point (0.,.5,1.);
478  nodes[15] = Point (.5,0.,0.);
479  nodes[16] = Point (.5,.5,0.);
480  nodes[17] = Point (0.,.5,0.);
481  return;
482  }
483  case PYRAMID5:
484  {
485  nodes.resize(5);
486  nodes[0] = Point (-1.,-1.,0.);
487  nodes[1] = Point (1.,-1.,0.);
488  nodes[2] = Point (1.,1.,0.);
489  nodes[3] = Point (-1.,1.,0.);
490  nodes[4] = Point (0.,0.,1.);
491  return;
492  }
493  case PYRAMID13:
494  {
495  nodes.resize(13);
496 
497  // base corners
498  nodes[0] = Point (-1.,-1.,0.);
499  nodes[1] = Point (1.,-1.,0.);
500  nodes[2] = Point (1.,1.,0.);
501  nodes[3] = Point (-1.,1.,0.);
502 
503  // apex
504  nodes[4] = Point (0.,0.,1.);
505 
506  // base midedge
507  nodes[5] = Point (0.,-1.,0.);
508  nodes[6] = Point (1.,0.,0.);
509  nodes[7] = Point (0.,1.,0.);
510  nodes[8] = Point (-1,0.,0.);
511 
512  // lateral midedge
513  nodes[9] = Point (-.5,-.5,.5);
514  nodes[10] = Point (.5,-.5,.5);
515  nodes[11] = Point (.5,.5,.5);
516  nodes[12] = Point (-.5,.5,.5);
517 
518  return;
519  }
520  case PYRAMID14:
521  {
522  nodes.resize(14);
523 
524  // base corners
525  nodes[0] = Point (-1.,-1.,0.);
526  nodes[1] = Point (1.,-1.,0.);
527  nodes[2] = Point (1.,1.,0.);
528  nodes[3] = Point (-1.,1.,0.);
529 
530  // apex
531  nodes[4] = Point (0.,0.,1.);
532 
533  // base midedge
534  nodes[5] = Point (0.,-1.,0.);
535  nodes[6] = Point (1.,0.,0.);
536  nodes[7] = Point (0.,1.,0.);
537  nodes[8] = Point (-1,0.,0.);
538 
539  // lateral midedge
540  nodes[9] = Point (-.5,-.5,.5);
541  nodes[10] = Point (.5,-.5,.5);
542  nodes[11] = Point (.5,.5,.5);
543  nodes[12] = Point (-.5,.5,.5);
544 
545  // base center
546  nodes[13] = Point (0.,0.,0.);
547 
548  return;
549  }
550 
551  default:
552  libmesh_error_msg("ERROR: Unknown element type " << itemType);
553  }
554 }
template<typename OutputType>
const std::vector<RealGradient>& libMesh::FEGenericBase< OutputType >::get_Sobolev_dweight ( ) const
inline
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 428 of file fe_base.h.

429  { return dweight; }
std::vector< RealGradient > dweight
Definition: fe_base.h:637
template<typename OutputType>
const std::vector<Real>& libMesh::FEGenericBase< OutputType >::get_Sobolev_weight ( ) const
inline
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 produces correct element matrices for a mesh using both finite and infinite elements.

Definition at line 420 of file fe_base.h.

421  { return weight; }
std::vector< Real > weight
Definition: fe_base.h:644
const std::vector<std::vector<Point> >& libMesh::FEAbstract::get_tangents ( ) const
inlineinherited
Returns
The tangent vectors for face integration.

Definition at line 365 of file fe_abstract.h.

References libMesh::FEAbstract::_fe_map.

366  { return this->_fe_map->get_tangents(); }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
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 404 of file fe_abstract.h.

References libMesh::FEAbstract::elem_type.

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

404 { return elem_type; }
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 185 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCounter::n_objects(), and libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

186 {
187  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
188  std::pair<unsigned int, unsigned int> & p = _counts[name];
189 
190  p.first++;
191 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
spin_mutex spin_mtx
Definition: threads.C:29
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 198 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCounter::n_objects(), and libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

199 {
200  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
201  std::pair<unsigned int, unsigned int> & p = _counts[name];
202 
203  p.second++;
204 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
spin_mutex spin_mtx
Definition: threads.C:29
virtual bool libMesh::FEAbstract::is_hierarchic ( ) const
pure virtualinherited
Returns
true if the finite element's higher order shape functions are hierarchic

Implemented in libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::InfFE< Dim, T_radial, T_map >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< 2, SUBDIVISION >, libMesh::FE< Dim, HIERARCHIC >, libMesh::FE< Dim, SCALAR >, libMesh::FE< Dim, L2_LAGRANGE >, libMesh::FE< Dim, NEDELEC_ONE >, libMesh::FE< Dim, HERMITE >, libMesh::FE< Dim, CLOUGH >, libMesh::FE< Dim, MONOMIAL >, libMesh::FE< Dim, XYZ >, libMesh::FE< Dim, LAGRANGE >, libMesh::FE< Dim, L2_HIERARCHIC >, libMesh::FE< Dim, LAGRANGE_VEC >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, and libMesh::FE< Dim, T >.

Referenced by libMesh::FEAbstract::set_fe_order().

static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited
virtual unsigned int libMesh::FEAbstract::n_shape_functions ( ) const
pure virtualinherited
bool libMesh::FEAbstract::on_reference_element ( const Point p,
const ElemType  t,
const Real  eps = TOLERANCE 
)
staticinherited
Returns
true if the point p is located on the reference element for element type t, false otherwise. Since we are doing floating point comparisons here the parameter eps can be specified to indicate a tolerance. For example, $ x \le 1 $ becomes $ x \le 1 + \epsilon $.

Definition at line 556 of file fe_abstract.C.

References libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::INFHEX8, libMesh::INFPRISM6, libMesh::NODEELEM, libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID5, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::QUADSHELL4, libMesh::Real, libMesh::TET10, libMesh::TET4, libMesh::TRI3, libMesh::TRI6, and libMesh::TRISHELL3.

Referenced by libMesh::FEInterface::ifem_on_reference_element(), libMesh::FE< Dim, T >::inverse_map(), and libMesh::FEInterface::on_reference_element().

557 {
558  libmesh_assert_greater_equal (eps, 0.);
559 
560  const Real xi = p(0);
561 #if LIBMESH_DIM > 1
562  const Real eta = p(1);
563 #else
564  const Real eta = 0.;
565 #endif
566 #if LIBMESH_DIM > 2
567  const Real zeta = p(2);
568 #else
569  const Real zeta = 0.;
570 #endif
571 
572  switch (t)
573  {
574  case NODEELEM:
575  {
576  return (!xi && !eta && !zeta);
577  }
578  case EDGE2:
579  case EDGE3:
580  case EDGE4:
581  {
582  // The reference 1D element is [-1,1].
583  if ((xi >= -1.-eps) &&
584  (xi <= 1.+eps))
585  return true;
586 
587  return false;
588  }
589 
590 
591  case TRI3:
592  case TRISHELL3:
593  case TRI6:
594  {
595  // The reference triangle is isocoles
596  // and is bound by xi=0, eta=0, and xi+eta=1.
597  if ((xi >= 0.-eps) &&
598  (eta >= 0.-eps) &&
599  ((xi + eta) <= 1.+eps))
600  return true;
601 
602  return false;
603  }
604 
605 
606  case QUAD4:
607  case QUADSHELL4:
608  case QUAD8:
609  case QUAD9:
610  {
611  // The reference quadrilateral element is [-1,1]^2.
612  if ((xi >= -1.-eps) &&
613  (xi <= 1.+eps) &&
614  (eta >= -1.-eps) &&
615  (eta <= 1.+eps))
616  return true;
617 
618  return false;
619  }
620 
621 
622  case TET4:
623  case TET10:
624  {
625  // The reference tetrahedral is isocoles
626  // and is bound by xi=0, eta=0, zeta=0,
627  // and xi+eta+zeta=1.
628  if ((xi >= 0.-eps) &&
629  (eta >= 0.-eps) &&
630  (zeta >= 0.-eps) &&
631  ((xi + eta + zeta) <= 1.+eps))
632  return true;
633 
634  return false;
635  }
636 
637 
638  case HEX8:
639  case HEX20:
640  case HEX27:
641  {
642  /*
643  if ((xi >= -1.) &&
644  (xi <= 1.) &&
645  (eta >= -1.) &&
646  (eta <= 1.) &&
647  (zeta >= -1.) &&
648  (zeta <= 1.))
649  return true;
650  */
651 
652  // The reference hexahedral element is [-1,1]^3.
653  if ((xi >= -1.-eps) &&
654  (xi <= 1.+eps) &&
655  (eta >= -1.-eps) &&
656  (eta <= 1.+eps) &&
657  (zeta >= -1.-eps) &&
658  (zeta <= 1.+eps))
659  {
660  // libMesh::out << "Strange Point:\n";
661  // p.print();
662  return true;
663  }
664 
665  return false;
666  }
667 
668  case PRISM6:
669  case PRISM15:
670  case PRISM18:
671  {
672  // Figure this one out...
673  // inside the reference triange with zeta in [-1,1]
674  if ((xi >= 0.-eps) &&
675  (eta >= 0.-eps) &&
676  (zeta >= -1.-eps) &&
677  (zeta <= 1.+eps) &&
678  ((xi + eta) <= 1.+eps))
679  return true;
680 
681  return false;
682  }
683 
684 
685  case PYRAMID5:
686  case PYRAMID13:
687  case PYRAMID14:
688  {
689  // Check that the point is on the same side of all the faces
690  // by testing whether:
691  //
692  // n_i.(x - x_i) <= 0
693  //
694  // for each i, where:
695  // n_i is the outward normal of face i,
696  // x_i is a point on face i.
697  if ((-eta - 1. + zeta <= 0.+eps) &&
698  ( xi - 1. + zeta <= 0.+eps) &&
699  ( eta - 1. + zeta <= 0.+eps) &&
700  ( -xi - 1. + zeta <= 0.+eps) &&
701  ( zeta >= 0.-eps))
702  return true;
703 
704  return false;
705  }
706 
707 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
708  case INFHEX8:
709  {
710  // The reference infhex8 is a [-1,1]^3.
711  if ((xi >= -1.-eps) &&
712  (xi <= 1.+eps) &&
713  (eta >= -1.-eps) &&
714  (eta <= 1.+eps) &&
715  (zeta >= -1.-eps) &&
716  (zeta <= 1.+eps))
717  {
718  return true;
719  }
720  return false;
721  }
722 
723  case INFPRISM6:
724  {
725  // inside the reference triange with zeta in [-1,1]
726  if ((xi >= 0.-eps) &&
727  (eta >= 0.-eps) &&
728  (zeta >= -1.-eps) &&
729  (zeta <= 1.+eps) &&
730  ((xi + eta) <= 1.+eps))
731  {
732  return true;
733  }
734 
735  return false;
736  }
737 #endif
738 
739  default:
740  libmesh_error_msg("ERROR: Unknown element type " << t);
741  }
742 
743  // If we get here then the point is _not_ in the
744  // reference element. Better return false.
745 
746  return false;
747 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::print_d2phi ( std::ostream &  os) const
virtual

Prints the value of each shape function's second derivatives at each quadrature point.

Implements libMesh::FEAbstract.

Definition at line 789 of file fe_base.C.

Referenced by libMesh::FEGenericBase< FEOutputType< T >::type >::get_Sobolev_dweight().

790 {
791  for (std::size_t i=0; i<dphi.size(); ++i)
792  for (std::size_t j=0; j<dphi[i].size(); ++j)
793  os << " d2phi[" << i << "][" << j << "]=" << d2phi[i][j];
794 }
std::vector< std::vector< OutputTensor > > d2phi
Definition: fe_base.h:552
std::vector< std::vector< OutputGradient > > dphi
Definition: fe_base.h:504
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::print_dphi ( std::ostream &  os) const
virtual

Prints the value of each shape function's derivative at each quadrature point.

Implements libMesh::FEAbstract.

Definition at line 731 of file fe_base.C.

Referenced by libMesh::FEGenericBase< FEOutputType< T >::type >::get_Sobolev_dweight().

732 {
733  for (std::size_t i=0; i<dphi.size(); ++i)
734  for (std::size_t j=0; j<dphi[i].size(); ++j)
735  os << " dphi[" << i << "][" << j << "]=" << dphi[i][j];
736 }
std::vector< std::vector< OutputGradient > > dphi
Definition: fe_base.h:504
void libMesh::ReferenceCounter::print_info ( std::ostream &  out = libMesh::out)
staticinherited

Prints the reference information, by default to libMesh::out.

Definition at line 88 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

Referenced by libMesh::LibMeshInit::LibMeshInit().

89 {
91  out_stream << ReferenceCounter::get_info();
92 }
static std::string get_info()
void libMesh::FEAbstract::print_info ( std::ostream &  os) const
inherited

Prints all the relevant information about the current element.

Definition at line 764 of file fe_abstract.C.

References libMesh::FEAbstract::print_dphi(), libMesh::FEAbstract::print_JxW(), libMesh::FEAbstract::print_phi(), and libMesh::FEAbstract::print_xyz().

Referenced by libMesh::FEAbstract::get_fe_map(), and libMesh::operator<<().

765 {
766  os << "phi[i][j]: Shape function i at quadrature pt. j" << std::endl;
767  this->print_phi(os);
768 
769  os << "dphi[i][j]: Shape function i's gradient at quadrature pt. j" << std::endl;
770  this->print_dphi(os);
771 
772  os << "XYZ locations of the quadrature pts." << std::endl;
773  this->print_xyz(os);
774 
775  os << "Values of JxW at the quadrature pts." << std::endl;
776  this->print_JxW(os);
777 }
void print_JxW(std::ostream &os) const
Definition: fe_abstract.C:751
void print_xyz(std::ostream &os) const
Definition: fe_abstract.C:758
virtual void print_phi(std::ostream &os) const =0
virtual void print_dphi(std::ostream &os) const =0
void libMesh::FEAbstract::print_JxW ( std::ostream &  os) const
inherited

Prints the Jacobian times the weight for each quadrature point.

Definition at line 751 of file fe_abstract.C.

References libMesh::FEAbstract::_fe_map.

Referenced by libMesh::FEAbstract::get_fe_map(), and libMesh::FEAbstract::print_info().

752 {
753  this->_fe_map->print_JxW(os);
754 }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
template<typename OutputType >
void libMesh::FEGenericBase< OutputType >::print_phi ( std::ostream &  os) const
virtual

Prints the value of each shape function at each quadrature point.

Implements libMesh::FEAbstract.

Definition at line 720 of file fe_base.C.

Referenced by libMesh::FEGenericBase< FEOutputType< T >::type >::get_Sobolev_dweight().

721 {
722  for (std::size_t i=0; i<phi.size(); ++i)
723  for (std::size_t j=0; j<phi[i].size(); ++j)
724  os << " phi[" << i << "][" << j << "]=" << phi[i][j] << std::endl;
725 }
std::vector< std::vector< OutputShape > > phi
Definition: fe_base.h:499
void libMesh::FEAbstract::print_xyz ( std::ostream &  os) const
inherited

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

Definition at line 758 of file fe_abstract.C.

References libMesh::FEAbstract::_fe_map.

Referenced by libMesh::FEAbstract::get_fe_map(), and libMesh::FEAbstract::print_info().

759 {
760  this->_fe_map->print_xyz(os);
761 }
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:511
virtual void libMesh::FEAbstract::reinit ( const Elem elem,
const std::vector< Point > *const  pts = libmesh_nullptr,
const std::vector< Real > *const  weights = libmesh_nullptr 
)
pure virtualinherited

This is at the core of this class. Use this for each new element in the mesh. Reinitializes the requested physical element-dependent data based on the current element elem. By default the element data are computed at the quadrature points specified by the quadrature rule qrule, but any set of points on the reference element may be specified in the optional argument pts.

Note
The FE classes decide which data to initialize based on which accessor functions such as get_phi() or get_d2phi() have been called, so all such accessors should be called before the first reinit().

Implemented in libMesh::FEXYZ< Dim >, libMesh::FESubdivision, libMesh::InfFE< Dim, T_radial, T_map >, libMesh::FE< Dim, T >, libMesh::FE< 2, SUBDIVISION >, libMesh::FE< Dim, HIERARCHIC >, libMesh::FE< Dim, SCALAR >, libMesh::FE< Dim, L2_LAGRANGE >, libMesh::FE< Dim, NEDELEC_ONE >, libMesh::FE< Dim, HERMITE >, libMesh::FE< Dim, CLOUGH >, libMesh::FE< Dim, MONOMIAL >, libMesh::FE< Dim, XYZ >, libMesh::FE< Dim, LAGRANGE >, libMesh::FE< Dim, L2_HIERARCHIC >, and libMesh::FE< Dim, LAGRANGE_VEC >.

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::FEMContext::build_new_fe(), libMesh::System::calculate_norm(), libMesh::ExactErrorEstimator::find_squared_element_error(), and libMesh::JumpErrorEstimator::reinit_sides().

virtual void libMesh::FEAbstract::reinit ( const Elem elem,
const unsigned int  side,
const Real  tolerance = TOLERANCE,
const std::vector< Point > *const  pts = libmesh_nullptr,
const std::vector< Real > *const  weights = libmesh_nullptr 
)
pure virtualinherited

Reinitializes all the physical element-dependent data based on the side of the element elem. The tolerance paremeter is passed to the involved call to inverse_map(). By default the element data are computed at the quadrature points specified by the quadrature rule qrule, but any set of points on the reference side element may be specified in the optional argument pts.

Implemented in libMesh::FEXYZ< Dim >, libMesh::FESubdivision, libMesh::InfFE< Dim, T_radial, T_map >, libMesh::FE< Dim, T >, libMesh::FE< 2, SUBDIVISION >, libMesh::FE< Dim, HIERARCHIC >, libMesh::FE< Dim, SCALAR >, libMesh::FE< Dim, L2_LAGRANGE >, libMesh::FE< Dim, NEDELEC_ONE >, libMesh::FE< Dim, HERMITE >, libMesh::FE< Dim, CLOUGH >, libMesh::FE< Dim, MONOMIAL >, libMesh::FE< Dim, XYZ >, libMesh::FE< Dim, LAGRANGE >, libMesh::FE< Dim, L2_HIERARCHIC >, libMesh::FE< Dim, LAGRANGE_VEC >, libMesh::FEXYZ< Dim >, and libMesh::FEXYZ< Dim >.

void libMesh::FEAbstract::set_fe_order ( int  new_order)
inlineinherited

Sets the base FE order of the finite element.

Definition at line 425 of file fe_abstract.h.

References libMesh::FEAbstract::fe_type, libMesh::FEAbstract::get_continuity(), libMesh::FEAbstract::is_hierarchic(), and libMesh::FEType::order.

425 { fe_type.order = new_order; }
OrderWrapper order
Definition: fe_type.h:197
virtual bool libMesh::FEAbstract::shapes_need_reinit ( ) const
protectedpure virtualinherited
Returns
true when the shape functions (for this FEFamily) depend on the particular element, and therefore needs to be re-initialized for each new element. false otherwise.

Implemented in libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::InfFE< Dim, T_radial, T_map >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< 2, SUBDIVISION >, libMesh::FE< Dim, HIERARCHIC >, libMesh::FE< Dim, SCALAR >, libMesh::FE< Dim, L2_LAGRANGE >, libMesh::FE< Dim, NEDELEC_ONE >, libMesh::FE< Dim, HERMITE >, libMesh::FE< Dim, CLOUGH >, libMesh::FE< Dim, MONOMIAL >, libMesh::FE< Dim, XYZ >, libMesh::FE< Dim, LAGRANGE >, libMesh::FE< Dim, L2_HIERARCHIC >, libMesh::FE< Dim, LAGRANGE_VEC >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, libMesh::FE< Dim, T >, and libMesh::FE< Dim, T >.

virtual void libMesh::FEAbstract::side_map ( const Elem elem,
const Elem side,
const unsigned int  s,
const std::vector< Point > &  reference_side_points,
std::vector< Point > &  reference_points 
)
pure virtualinherited

Friends And Related Function Documentation

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

Make all InfFE<Dim,T_radial,T_map> classes friends so that they can safely used FE<Dim-1,T_base> through a FEGenericBase * as base approximation.

Definition at line 658 of file fe_base.h.

Member Data Documentation

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts
staticprotectedinherited
bool libMesh::ReferenceCounter::_enable_print_counter = true
staticprotectedinherited

Flag to control whether reference count information is printed when print_info is called.

Definition at line 143 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().

template<typename OutputType>