#include <fe_base.h>
Inherits ReferenceCountedObject< FEBase >.
Inherited by FE< Dim, T >, FE< Dim, CLOUGH >, FE< Dim, HERMITE >, FE< Dim, HIERARCHIC >, FE< Dim, LAGRANGE >, FE< Dim, MONOMIAL >, FE< Dim, SCALAR >, FE< Dim, XYZ >, and InfFE< Dim, T_radial, T_map >.
virtual ~FEBase ()
virtual void reinit (const Elem *elem, const std::vector< Point > *const pts=NULL)=0
virtual void reinit (const Elem *elem, const unsigned int side, const Real tolerance=TOLERANCE)=0
virtual void edge_reinit (const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE)=0
const std::vector< Point > & get_xyz () const
const std::vector< std::vector< Real > > & get_phi () const
const std::vector< Real > & get_JxW () const
const std::vector< std::vector< RealGradient > > & get_dphi () const
const std::vector< std::vector< Real > > & get_dphidx () const
const std::vector< std::vector< Real > > & get_dphidy () const
const std::vector< std::vector< Real > > & get_dphidz () const
const std::vector< std::vector< Real > > & get_dphidxi () const
const std::vector< std::vector< Real > > & get_dphideta () const
const std::vector< std::vector< Real > > & get_dphidzeta () const
const std::vector< std::vector< RealTensor > > & get_d2phi () const
const std::vector< std::vector< Real > > & get_d2phidx2 () const
const std::vector< std::vector< Real > > & get_d2phidxdy () const
const std::vector< std::vector< Real > > & get_d2phidxdz () const
const std::vector< std::vector< Real > > & get_d2phidy2 () const
const std::vector< std::vector< Real > > & get_d2phidydz () const
const std::vector< std::vector< Real > > & get_d2phidz2 () 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< RealGradient > & get_dphase () const
const std::vector< Real > & get_Sobolev_weight () const
const std::vector< RealGradient > & get_Sobolev_dweight () 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
virtual FEContinuity get_continuity () const =0
virtual bool is_hierarchic () const =0
FEFamily get_family () const
void print_JxW (std::ostream &os) const
void print_phi (std::ostream &os) const
void print_dphi (std::ostream &os) const
void print_d2phi (std::ostream &os) const
void print_xyz (std::ostream &os) const
void print_info (std::ostream &os) const
static AutoPtr< FEBase > build (const unsigned int dim, const FEType &type)
static AutoPtr< FEBase > build_InfFE (const unsigned int dim, const FEType &type)
static bool on_reference_element (const Point &p, const ElemType t, const Real eps=TOLERANCE)
static void compute_proj_constraints (DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
static void coarsened_dof_values (const NumericVector< Number > &global_vector, const DofMap &dof_map, const Elem *coarse_elem, DenseVector< Number > &coarse_dofs, const unsigned int var, const bool use_old_dof_indices=false)
static void compute_periodic_constraints (DofConstraints &constraints, DofMap &dof_map, PeriodicBoundaries &boundaries, const MeshBase &mesh, const unsigned int variable_number, const Elem *elem)
static std::string get_info ()
static void print_info ()
static unsigned int n_objects ()
typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
FEBase (const unsigned int dim, const FEType &fet)
virtual void init_base_shape_functions (const std::vector< Point > &qp, const Elem *e)=0
virtual void compute_map (const std::vector< Real > &qw, const Elem *e)
virtual void compute_affine_map (const std::vector< Real > &qw, const Elem *e)
void compute_single_point_map (const std::vector< Real > &qw, const Elem *e, unsigned int p)
void resize_map_vectors (unsigned int n_qp)
void compute_face_map (const std::vector< Real > &qw, const Elem *side)
void compute_edge_map (const std::vector< Real > &qw, const Elem *side)
virtual void compute_shape_functions (const Elem *)
Real dxdxi_map (const unsigned int p) const
Real dydxi_map (const unsigned int p) const
Real dzdxi_map (const unsigned int p) const
Real dxdeta_map (const unsigned int p) const
Real dydeta_map (const unsigned int p) const
Real dzdeta_map (const unsigned int p) const
Real dxdzeta_map (const unsigned int p) const
Real dydzeta_map (const unsigned int p) const
Real dzdzeta_map (const unsigned int p) const
void increment_constructor_count (const std::string &name)
void increment_destructor_count (const std::string &name)
const unsigned int dim
std::vector< Point > xyz
std::vector< RealGradient > dxyzdxi_map
std::vector< RealGradient > dxyzdeta_map
std::vector< RealGradient > dxyzdzeta_map
std::vector< RealGradient > d2xyzdxi2_map
std::vector< RealGradient > d2xyzdxideta_map
std::vector< RealGradient > d2xyzdeta2_map
std::vector< RealGradient > d2xyzdxidzeta_map
std::vector< RealGradient > d2xyzdetadzeta_map
std::vector< RealGradient > d2xyzdzeta2_map
std::vector< Real > dxidx_map
std::vector< Real > dxidy_map
std::vector< Real > dxidz_map
std::vector< Real > detadx_map
std::vector< Real > detady_map
std::vector< Real > detadz_map
std::vector< Real > dzetadx_map
std::vector< Real > dzetady_map
std::vector< Real > dzetadz_map
std::vector< std::vector< Real > > phi
bool calculations_started
bool calculate_phi
bool calculate_dphi
bool calculate_d2phi
std::vector< std::vector< RealGradient > > dphi
std::vector< std::vector< Real > > dphidxi
std::vector< std::vector< Real > > dphideta
std::vector< std::vector< Real > > dphidzeta
std::vector< std::vector< Real > > dphidx
std::vector< std::vector< Real > > dphidy
std::vector< std::vector< Real > > dphidz
std::vector< std::vector< RealTensor > > d2phi
std::vector< std::vector< Real > > d2phidxi2
std::vector< std::vector< Real > > d2phidxideta
std::vector< std::vector< Real > > d2phidxidzeta
std::vector< std::vector< Real > > d2phideta2
std::vector< std::vector< Real > > d2phidetadzeta
std::vector< std::vector< Real > > d2phidzeta2
std::vector< std::vector< Real > > d2phidx2
std::vector< std::vector< Real > > d2phidxdy
std::vector< std::vector< Real > > d2phidxdz
std::vector< std::vector< Real > > d2phidy2
std::vector< std::vector< Real > > d2phidydz
std::vector< std::vector< Real > > d2phidz2
std::vector< std::vector< Real > > phi_map
std::vector< std::vector< Real > > dphidxi_map
std::vector< std::vector< Real > > dphideta_map
std::vector< std::vector< Real > > dphidzeta_map
std::vector< std::vector< Real > > d2phidxi2_map
std::vector< std::vector< Real > > d2phidxideta_map
std::vector< std::vector< Real > > d2phidxidzeta_map
std::vector< std::vector< Real > > d2phideta2_map
std::vector< std::vector< Real > > d2phidetadzeta_map
std::vector< std::vector< Real > > d2phidzeta2_map
std::vector< std::vector< Real > > psi_map
std::vector< std::vector< Real > > dpsidxi_map
std::vector< std::vector< Real > > dpsideta_map
std::vector< std::vector< Real > > d2psidxi2_map
std::vector< std::vector< Real > > d2psidxideta_map
std::vector< std::vector< Real > > d2psideta2_map
std::vector< RealGradient > dphase
std::vector< RealGradient > dweight
std::vector< Real > weight
std::vector< std::vector< Point > > tangents
std::vector< Point > normals
std::vector< Real > curvatures
std::vector< Real > JxW
const FEType fe_type
ElemType elem_type
unsigned int _p_level
QBase * qrule
bool shapes_on_quadrature
static Counts _counts
static Threads::atomic< unsigned int > _n_objects
static Threads::spin_mutex _mutex
virtual bool shapes_need_reinit () const =0
class InfFE
std::ostream & operator<< (std::ostream &os, const FEBase &fe)
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 FEBase::build() method to create an object of any of the derived classes. Note that the amount of virtual members is kept to a minimum, and the sophisticated template scheme of FE is quite likely to offer acceptably fast code.
All calls to static members of the FE classes should be requested through the FEInterface. This interface class offers sort-of runtime polymorphism for the templated finite element classes. Even internal library classes, like DofMap, request the number of dof's through this interface class. Note that this also enables the co-existence of various element-based schemes. This class is well 'at the heart' of the library, so things in here should better remain unchanged.
Author:
Definition at line 88 of file fe_base.h.
Definition at line 105 of file reference_counter.h.
Definition at line 1261 of file fe_base.h.
:
dim(d),
calculations_started(false),
calculate_phi(false),
calculate_dphi(false),
calculate_d2phi(false),
fe_type(fet),
elem_type(INVALID_ELEM),
_p_level(0),
qrule(NULL),
shapes_on_quadrature(false)
{
}
Definition at line 1279 of file fe_base.h.
{
}
Implemented in FE< Dim, T >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
Referenced by InfFE< Dim, T_radial, T_map >::init_face_shape_functions().
Definition at line 43 of file fe_base.C.
References libMeshEnums::BERNSTEIN, libMeshEnums::CLOUGH, FEType::family, libMeshEnums::HERMITE, libMeshEnums::HIERARCHIC, libMeshEnums::LAGRANGE, libMeshEnums::MONOMIAL, libMeshEnums::SCALAR, libMeshEnums::SZABAB, and libMeshEnums::XYZ.
Referenced by ExactSolution::_compute_error(), UniformRefinementEstimator::_estimate_error(), System::calculate_norm(), coarsened_dof_values(), compute_periodic_constraints(), compute_proj_constraints(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), FEMContext::FEMContext(), MeshFunction::gradient(), MeshFunction::hessian(), InfFE< Dim, T_radial, T_map >::InfFE(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), System::ProjectVector::operator()(), PatchRecoveryErrorEstimator::EstimateError::operator()(), InfFE< Dim, T_radial, T_map >::reinit(), HPCoarsenTest::select_refinement(), DofMap::use_coupled_neighbor_dofs(), and Elem::volume().
{
// The stupid AutoPtr<FEBase> ap(); return ap;
// construct is required to satisfy IBM's xlC
switch (dim)
{
// 0D
case 0:
{
switch (fet.family)
{
case CLOUGH:
{
AutoPtr<FEBase> ap(new FE<0,CLOUGH>(fet));
return ap;
}
case HERMITE:
{
AutoPtr<FEBase> ap(new FE<0,HERMITE>(fet));
return ap;
}
case LAGRANGE:
{
AutoPtr<FEBase> ap(new FE<0,LAGRANGE>(fet));
return ap;
}
case HIERARCHIC:
{
AutoPtr<FEBase> ap(new FE<0,HIERARCHIC>(fet));
return ap;
}
case MONOMIAL:
{
AutoPtr<FEBase> ap(new FE<0,MONOMIAL>(fet));
return ap;
}
#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
case SZABAB:
{
AutoPtr<FEBase> ap(new FE<0,SZABAB>(fet));
return ap;
}
case BERNSTEIN:
{
AutoPtr<FEBase> ap(new FE<0,BERNSTEIN>(fet));
return ap;
}
#endif
case XYZ:
{
AutoPtr<FEBase> ap(new FEXYZ<0>(fet));
return ap;
}
case SCALAR:
{
AutoPtr<FEBase> ap(new FEScalar<0>(fet));
return ap;
}
default:
std::cout << 'ERROR: Bad FEType.family= ' << fet.family << std::endl;
libmesh_error();
}
}
// 1D
case 1:
{
switch (fet.family)
{
case CLOUGH:
{
AutoPtr<FEBase> ap(new FE<1,CLOUGH>(fet));
return ap;
}
case HERMITE:
{
AutoPtr<FEBase> ap(new FE<1,HERMITE>(fet));
return ap;
}
case LAGRANGE:
{
AutoPtr<FEBase> ap(new FE<1,LAGRANGE>(fet));
return ap;
}
case HIERARCHIC:
{
AutoPtr<FEBase> ap(new FE<1,HIERARCHIC>(fet));
return ap;
}
case MONOMIAL:
{
AutoPtr<FEBase> ap(new FE<1,MONOMIAL>(fet));
return ap;
}
#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
case SZABAB:
{
AutoPtr<FEBase> ap(new FE<1,SZABAB>(fet));
return ap;
}
case BERNSTEIN:
{
AutoPtr<FEBase> ap(new FE<1,BERNSTEIN>(fet));
return ap;
}
#endif
case XYZ:
{
AutoPtr<FEBase> ap(new FEXYZ<1>(fet));
return ap;
}
case SCALAR:
{
AutoPtr<FEBase> ap(new FEScalar<1>(fet));
return ap;
}
default:
std::cout << 'ERROR: Bad FEType.family= ' << fet.family << std::endl;
libmesh_error();
}
}
// 2D
case 2:
{
switch (fet.family)
{
case CLOUGH:
{
AutoPtr<FEBase> ap(new FE<2,CLOUGH>(fet));
return ap;
}
case HERMITE:
{
AutoPtr<FEBase> ap(new FE<2,HERMITE>(fet));
return ap;
}
case LAGRANGE:
{
AutoPtr<FEBase> ap(new FE<2,LAGRANGE>(fet));
return ap;
}
case HIERARCHIC:
{
AutoPtr<FEBase> ap(new FE<2,HIERARCHIC>(fet));
return ap;
}
case MONOMIAL:
{
AutoPtr<FEBase> ap(new FE<2,MONOMIAL>(fet));
return ap;
}
#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
case SZABAB:
{
AutoPtr<FEBase> ap(new FE<2,SZABAB>(fet));
return ap;
}
case BERNSTEIN:
{
AutoPtr<FEBase> ap(new FE<2,BERNSTEIN>(fet));
return ap;
}
#endif
case XYZ:
{
AutoPtr<FEBase> ap(new FEXYZ<2>(fet));
return ap;
}
case SCALAR:
{
AutoPtr<FEBase> ap(new FEScalar<2>(fet));
return ap;
}
default:
std::cout << 'ERROR: Bad FEType.family= ' << fet.family << std::endl;
libmesh_error();
}
}
// 3D
case 3:
{
switch (fet.family)
{
case CLOUGH:
{
std::cout << 'ERROR: Clough-Tocher elements currently only support 1D and 2D' <<
std::endl;
libmesh_error();
}
case HERMITE:
{
AutoPtr<FEBase> ap(new FE<3,HERMITE>(fet));
return ap;
}
case LAGRANGE:
{
AutoPtr<FEBase> ap(new FE<3,LAGRANGE>(fet));
return ap;
}
case HIERARCHIC:
{
AutoPtr<FEBase> ap(new FE<3,HIERARCHIC>(fet));
return ap;
}
case MONOMIAL:
{
AutoPtr<FEBase> ap(new FE<3,MONOMIAL>(fet));
return ap;
}
#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
case SZABAB:
{
AutoPtr<FEBase> ap(new FE<3,SZABAB>(fet));
return ap;
}
case BERNSTEIN:
{
AutoPtr<FEBase> ap(new FE<3,BERNSTEIN>(fet));
return ap;
}
#endif
case XYZ:
{
AutoPtr<FEBase> ap(new FEXYZ<3>(fet));
return ap;
}
case SCALAR:
{
AutoPtr<FEBase> ap(new FEScalar<3>(fet));
return ap;
}
default:
std::cout << 'ERROR: Bad FEType.family= ' << fet.family << std::endl;
libmesh_error();
}
}
default:
libmesh_error();
}
libmesh_error();
AutoPtr<FEBase> ap(NULL);
return ap;
}
Definition at line 339 of file fe_base.C.
References libMeshEnums::CARTESIAN, FEType::inf_map, libMeshEnums::INFINITE_MAP, libMeshEnums::JACOBI_20_00, libMeshEnums::JACOBI_30_00, libMeshEnums::LAGRANGE, libMeshEnums::LEGENDRE, and FEType::radial_family.
{
// The stupid AutoPtr<FEBase> ap(); return ap;
// construct is required to satisfy IBM's xlC
switch (dim)
{
// 1D
case 1:
{
switch (fet.radial_family)
{
case INFINITE_MAP:
{
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with FEFamily = ' << fet.radial_family << std::endl;
libmesh_error();
}
case JACOBI_20_00:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<1,JACOBI_20_00,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
case JACOBI_30_00:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<1,JACOBI_30_00,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
case LEGENDRE:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<1,LEGENDRE,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
case LAGRANGE:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<1,LAGRANGE,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
default:
std::cerr << 'ERROR: Bad FEType.radial_family= ' << fet.radial_family << std::endl;
libmesh_error();
}
}
// 2D
case 2:
{
switch (fet.radial_family)
{
case INFINITE_MAP:
{
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with FEFamily = ' << fet.radial_family << std::endl;
libmesh_error();
}
case JACOBI_20_00:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<2,JACOBI_20_00,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
case JACOBI_30_00:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<2,JACOBI_30_00,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
case LEGENDRE:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<2,LEGENDRE,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
case LAGRANGE:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<2,LAGRANGE,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
default:
std::cerr << 'ERROR: Bad FEType.radial_family= ' << fet.radial_family << std::endl;
libmesh_error();
}
}
// 3D
case 3:
{
switch (fet.radial_family)
{
case INFINITE_MAP:
{
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with FEFamily = ' << fet.radial_family << std::endl;
libmesh_error();
}
case JACOBI_20_00:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<3,JACOBI_20_00,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
case JACOBI_30_00:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<3,JACOBI_30_00,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
case LEGENDRE:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<3,LEGENDRE,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
case LAGRANGE:
{
switch (fet.inf_map)
{
case CARTESIAN:
{
AutoPtr<FEBase> ap(new InfFE<3,LAGRANGE,CARTESIAN>(fet));
return ap;
}
default:
std::cerr << 'ERROR: Don't build an infinite element ' << std::endl
<< ' with InfMapType = ' << fet.inf_map << std::endl;
libmesh_error();
}
}
default:
std::cerr << 'ERROR: Bad FEType.radial_family= ' << fet.radial_family << std::endl;
libmesh_error();
}
}
default:
libmesh_error();
}
libmesh_error();
AutoPtr<FEBase> ap(NULL);
return ap;
}
Definition at line 1118 of file fe_base.C.
References TypeVector< T >::add_scaled(), build(), libMeshEnums::C_ONE, Elem::child(), DenseMatrix< T >::cholesky_solve(), FEType::default_quadrature_rule(), Elem::dim(), dim, libMeshEnums::DISCONTINUOUS, DofMap::dof_indices(), FEInterface::dofs_on_edge(), FEInterface::dofs_on_side(), elem_type, fe_type, FEInterface::inverse_map(), Elem::is_child_on_edge(), Elem::is_child_on_side(), Elem::is_vertex(), JxW, Elem::max_descendant_p_level(), Elem::n_children(), FEInterface::n_dofs(), FEInterface::n_dofs_at_node(), Elem::n_edges(), Elem::n_nodes(), MeshTools::n_nodes(), QBase::n_points(), Elem::n_sides(), DofMap::old_dof_indices(), FEType::order, Elem::p_level(), qrule, DenseMatrix< T >::resize(), DenseVector< T >::resize(), Elem::type(), DofMap::variable_type(), libMesh::zero, DenseMatrix< T >::zero(), and DenseVector< T >::zero().
Referenced by JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), and System::ProjectVector::operator()().
{
// Side/edge DOF indices
std::vector<unsigned int> new_side_dofs, old_side_dofs;
// FIXME: what about 2D shells in 3D space?
unsigned int dim = elem->dim();
// We use local FE objects for now
// FIXME: we should use more, external objects instead for efficiency
const FEType& base_fe_type = dof_map.variable_type(var);
AutoPtr<FEBase> fe (FEBase::build(dim, base_fe_type));
AutoPtr<FEBase> fe_coarse (FEBase::build(dim, base_fe_type));
AutoPtr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
AutoPtr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
AutoPtr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
std::vector<Point> coarse_qpoints;
// The values of the shape functions at the quadrature
// points
const std::vector<std::vector<Real> >& phi_values =
fe->get_phi();
const std::vector<std::vector<Real> >& phi_coarse =
fe_coarse->get_phi();
// The gradients of the shape functions at the quadrature
// points on the child element.
const std::vector<std::vector<RealGradient> > *dphi_values =
NULL;
const std::vector<std::vector<RealGradient> > *dphi_coarse =
NULL;
const FEContinuity cont = fe->get_continuity();
if (cont == C_ONE)
{
const std::vector<std::vector<RealGradient> >&
ref_dphi_values = fe->get_dphi();
dphi_values = &ref_dphi_values;
const std::vector<std::vector<RealGradient> >&
ref_dphi_coarse = fe_coarse->get_dphi();
dphi_coarse = &ref_dphi_coarse;
}
// The Jacobian * quadrature weight at the quadrature points
const std::vector<Real>& JxW =
fe->get_JxW();
// The XYZ locations of the quadrature points on the
// child element
const std::vector<Point>& xyz_values =
fe->get_xyz();
FEType fe_type = base_fe_type, temp_fe_type;
const ElemType elem_type = elem->type();
fe_type.order = static_cast<Order>(fe_type.order +
elem->max_descendant_p_level());
// Number of nodes on parent element
const unsigned int n_nodes = elem->n_nodes();
// Number of dofs on parent element
const unsigned int new_n_dofs =
FEInterface::n_dofs(dim, fe_type, elem_type);
// Fixed vs. free DoFs on edge/face projections
std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
std::vector<int> free_dof(new_n_dofs, 0);
DenseMatrix<Real> Ke;
DenseVector<Number> Fe;
Ue.resize(new_n_dofs); Ue.zero();
// When coarsening, in general, we need a series of
// projections to ensure a unique and continuous
// solution. We start by interpolating nodes, then
// hold those fixed and project edges, then
// hold those fixed and project faces, then
// hold those fixed and project interiors
// Copy node values first
{
std::vector<unsigned int> node_dof_indices;
if (use_old_dof_indices)
dof_map.old_dof_indices (elem, node_dof_indices, var);
else
dof_map.dof_indices (elem, node_dof_indices, var);
unsigned int current_dof = 0;
for (unsigned int n=0; n!= n_nodes; ++n)
{
// FIXME: this should go through the DofMap,
// not duplicate dof_indices code badly!
const unsigned int my_nc =
FEInterface::n_dofs_at_node (dim, fe_type,
elem_type, n);
if (!elem->is_vertex(n))
{
current_dof += my_nc;
continue;
}
temp_fe_type = base_fe_type;
// We're assuming here that child n shares vertex n,
// which is wrong on non-simplices right now
// ... but this code isn't necessary except on elements
// where p refinement creates more vertex dofs; we have
// no such elements yet.
/*
if (elem->child(n)->p_level() < elem->p_level())
{
temp_fe_type.order =
static_cast<Order>(temp_fe_type.order +
elem->child(n)->p_level());
}
*/
const unsigned int nc =
FEInterface::n_dofs_at_node (dim, temp_fe_type,
elem_type, n);
for (unsigned int i=0; i!= nc; ++i)
{
Ue(current_dof) =
old_vector(node_dof_indices[current_dof]);
dof_is_fixed[current_dof] = true;
current_dof++;
}
}
}
// In 3D, project any edge values next
if (dim > 2 && cont != DISCONTINUOUS)
for (unsigned int e=0; e != elem->n_edges(); ++e)
{
FEInterface::dofs_on_edge(elem, dim, fe_type,
e, new_side_dofs);
// Some edge dofs are on nodes and already
// fixed, others are free to calculate
unsigned int free_dofs = 0;
for (unsigned int i=0; i !=
new_side_dofs.size(); ++i)
if (!dof_is_fixed[new_side_dofs[i]])
free_dof[free_dofs++] = i;
Ke.resize (free_dofs, free_dofs); Ke.zero();
Fe.resize (free_dofs); Fe.zero();
// The new edge coefficients
DenseVector<Number> Uedge(free_dofs);
// Add projection terms from each child sharing
// this edge
for (unsigned int c=0; c != elem->n_children();
++c)
{
if (!elem->is_child_on_edge(c,e))
continue;
Elem *child = elem->child(c);
std::vector<unsigned int> child_dof_indices;
if (use_old_dof_indices)
dof_map.old_dof_indices (child,
child_dof_indices, var);
else
dof_map.dof_indices (child,
child_dof_indices, var);
const unsigned int child_n_dofs = child_dof_indices.size();
temp_fe_type = base_fe_type;
temp_fe_type.order =
static_cast<Order>(temp_fe_type.order +
child->p_level());
FEInterface::dofs_on_edge(child, dim,
temp_fe_type, e, old_side_dofs);
// Initialize both child and parent FE data
// on the child's edge
fe->attach_quadrature_rule (qedgerule.get());
fe->edge_reinit (child, e);
const unsigned int n_qp = qedgerule->n_points();
FEInterface::inverse_map (dim, fe_type, elem,
xyz_values, coarse_qpoints);
fe_coarse->reinit(elem, &coarse_qpoints);
// Loop over the quadrature points
for (unsigned int qp=0; qp<n_qp; qp++)
{
// solution value at the quadrature point
Number fineval = libMesh::zero;
// solution grad at the quadrature point
Gradient finegrad;
// Sum the solution values * the DOF
// values at the quadrature point to
// get the solution value and gradient.
for (unsigned int i=0; i<child_n_dofs;
i++)
{
fineval +=
(old_vector(child_dof_indices[i])*
phi_values[i][qp]);
if (cont == C_ONE)
finegrad.add_scaled((*dphi_values)[i][qp],
old_vector(child_dof_indices[i]));
}
// Form edge projection matrix
for (unsigned int sidei=0, freei=0;
sidei != new_side_dofs.size();
++sidei)
{
unsigned int i = new_side_dofs[sidei];
// fixed DoFs aren't test functions
if (dof_is_fixed[i])
continue;
for (unsigned int sidej=0, freej=0;
sidej != new_side_dofs.size();
++sidej)
{
unsigned int j =
new_side_dofs[sidej];
if (dof_is_fixed[j])
Fe(freei) -=
phi_coarse[i][qp] *
phi_coarse[j][qp] * JxW[qp] *
Ue(j);
else
Ke(freei,freej) +=
phi_coarse[i][qp] *
phi_coarse[j][qp] * JxW[qp];
if (cont == C_ONE)
{
if (dof_is_fixed[j])
Fe(freei) -=
((*dphi_coarse)[i][qp] *
(*dphi_coarse)[j][qp]) *
JxW[qp] *
Ue(j);
else
Ke(freei,freej) +=
((*dphi_coarse)[i][qp] *
(*dphi_coarse)[j][qp])
* JxW[qp];
}
if (!dof_is_fixed[j])
freej++;
}
Fe(freei) += phi_coarse[i][qp] *
fineval * JxW[qp];
if (cont == C_ONE)
Fe(freei) +=
(finegrad * (*dphi_coarse)[i][qp]) * JxW[qp];
freei++;
}
}
}
Ke.cholesky_solve(Fe, Uedge);
// Transfer new edge solutions to element
for (unsigned int i=0; i != free_dofs; ++i)
{
Number &ui = Ue(new_side_dofs[free_dof[i]]);
libmesh_assert(std::abs(ui) < TOLERANCE ||
std::abs(ui - Uedge(i)) < TOLERANCE);
ui = Uedge(i);
dof_is_fixed[new_side_dofs[free_dof[i]]] =
true;
}
}
// Project any side values (edges in 2D, faces in 3D)
if (dim > 1 && cont != DISCONTINUOUS)
for (unsigned int s=0; s != elem->n_sides(); ++s)
{
FEInterface::dofs_on_side(elem, dim, fe_type,
s, new_side_dofs);
// Some side dofs are on nodes/edges and already
// fixed, others are free to calculate
unsigned int free_dofs = 0;
for (unsigned int i=0; i !=
new_side_dofs.size(); ++i)
if (!dof_is_fixed[new_side_dofs[i]])
free_dof[free_dofs++] = i;
Ke.resize (free_dofs, free_dofs); Ke.zero();
Fe.resize (free_dofs); Fe.zero();
// The new side coefficients
DenseVector<Number> Uside(free_dofs);
// Add projection terms from each child sharing
// this side
for (unsigned int c=0; c != elem->n_children();
++c)
{
if (!elem->is_child_on_side(c,s))
continue;
Elem *child = elem->child(c);
std::vector<unsigned int> child_dof_indices;
if (use_old_dof_indices)
dof_map.old_dof_indices (child,
child_dof_indices, var);
else
dof_map.dof_indices (child,
child_dof_indices, var);
const unsigned int child_n_dofs = child_dof_indices.size();
temp_fe_type = base_fe_type;
temp_fe_type.order =
static_cast<Order>(temp_fe_type.order +
child->p_level());
FEInterface::dofs_on_side(child, dim,
temp_fe_type, s, old_side_dofs);
// Initialize both child and parent FE data
// on the child's side
fe->attach_quadrature_rule (qsiderule.get());
fe->reinit (child, s);
const unsigned int n_qp = qsiderule->n_points();
FEInterface::inverse_map (dim, fe_type, elem,
xyz_values, coarse_qpoints);
fe_coarse->reinit(elem, &coarse_qpoints);
// Loop over the quadrature points
for (unsigned int qp=0; qp<n_qp; qp++)
{
// solution value at the quadrature point
Number fineval = libMesh::zero;
// solution grad at the quadrature point
Gradient finegrad;
// Sum the solution values * the DOF
// values at the quadrature point to
// get the solution value and gradient.
for (unsigned int i=0; i<child_n_dofs;
i++)
{
fineval +=
(old_vector(child_dof_indices[i])*
phi_values[i][qp]);
if (cont == C_ONE)
finegrad.add_scaled((*dphi_values)[i][qp],
old_vector(child_dof_indices[i]));
}
// Form side projection matrix
for (unsigned int sidei=0, freei=0;
sidei != new_side_dofs.size();
++sidei)
{
unsigned int i = new_side_dofs[sidei];
// fixed DoFs aren't test functions
if (dof_is_fixed[i])
continue;
for (unsigned int sidej=0, freej=0;
sidej != new_side_dofs.size();
++sidej)
{
unsigned int j =
new_side_dofs[sidej];
if (dof_is_fixed[j])
Fe(freei) -=
phi_coarse[i][qp] *
phi_coarse[j][qp] * JxW[qp] *
Ue(j);
else
Ke(freei,freej) +=
phi_coarse[i][qp] *
phi_coarse[j][qp] * JxW[qp];
if (cont == C_ONE)
{
if (dof_is_fixed[j])
Fe(freei) -=
((*dphi_coarse)[i][qp] *
(*dphi_coarse)[j][qp]) *
JxW[qp] *
Ue(j);
else
Ke(freei,freej) +=
((*dphi_coarse)[i][qp] *
(*dphi_coarse)[j][qp])
* JxW[qp];
}
if (!dof_is_fixed[j])
freej++;
}
Fe(freei) += (fineval * phi_coarse[i][qp]) * JxW[qp];
if (cont == C_ONE)
Fe(freei) +=
(finegrad * (*dphi_coarse)[i][qp]) * JxW[qp];
freei++;
}
}
}
Ke.cholesky_solve(Fe, Uside);
// Transfer new side solutions to element
for (unsigned int i=0; i != free_dofs; ++i)
{
Number &ui = Ue(new_side_dofs[free_dof[i]]);
libmesh_assert(std::abs(ui) < TOLERANCE ||
std::abs(ui - Uside(i)) < TOLERANCE);
ui = Uside(i);
dof_is_fixed[new_side_dofs[free_dof[i]]] =
true;
}
}
// Project the interior values, finally
// Some interior dofs are on nodes/edges/sides and
// already fixed, others are free to calculate
unsigned int free_dofs = 0;
for (unsigned int i=0; i != new_n_dofs; ++i)
if (!dof_is_fixed[i])
free_dof[free_dofs++] = i;
Ke.resize (free_dofs, free_dofs); Ke.zero();
Fe.resize (free_dofs); Fe.zero();
// The new interior coefficients
DenseVector<Number> Uint(free_dofs);
// Add projection terms from each child
for (unsigned int c=0; c != elem->n_children(); ++c)
{
Elem *child = elem->child(c);
std::vector<unsigned int> child_dof_indices;
if (use_old_dof_indices)
dof_map.old_dof_indices (child,
child_dof_indices, var);
else
dof_map.dof_indices (child,
child_dof_indices, var);
const unsigned int child_n_dofs = child_dof_indices.size();
// Initialize both child and parent FE data
// on the child's quadrature points
fe->attach_quadrature_rule (qrule.get());
fe->reinit (child);
const unsigned int n_qp = qrule->n_points();
FEInterface::inverse_map (dim, fe_type, elem,
xyz_values, coarse_qpoints);
fe_coarse->reinit(elem, &coarse_qpoints);
// Loop over the quadrature points
for (unsigned int qp=0; qp<n_qp; qp++)
{
// solution value at the quadrature point
Number fineval = libMesh::zero;
// solution grad at the quadrature point
Gradient finegrad;
// Sum the solution values * the DOF
// values at the quadrature point to
// get the solution value and gradient.
for (unsigned int i=0; i<child_n_dofs; i++)
{
fineval +=
(old_vector(child_dof_indices[i])*
phi_values[i][qp]);
if (cont == C_ONE)
finegrad.add_scaled((*dphi_values)[i][qp],
old_vector(child_dof_indices[i]));
}
// Form interior projection matrix
for (unsigned int i=0, freei=0;
i != new_n_dofs; ++i)
{
// fixed DoFs aren't test functions
if (dof_is_fixed[i])
continue;
for (unsigned int j=0, freej=0; j !=
new_n_dofs; ++j)
{
if (dof_is_fixed[j])
Fe(freei) -=
phi_coarse[i][qp] *
phi_coarse[j][qp] * JxW[qp] *
Ue(j);
else
Ke(freei,freej) +=
phi_coarse[i][qp] *
phi_coarse[j][qp] * JxW[qp];
if (cont == C_ONE)
{
if (dof_is_fixed[j])
Fe(freei) -=
((*dphi_coarse)[i][qp] *
(*dphi_coarse)[j][qp]) *
JxW[qp] * Ue(j);
else
Ke(freei,freej) +=
((*dphi_coarse)[i][qp] *
(*dphi_coarse)[j][qp]) * JxW[qp];
}
if (!dof_is_fixed[j])
freej++;
}
Fe(freei) += phi_coarse[i][qp] * fineval *
JxW[qp];
if (cont == C_ONE)
Fe(freei) += (finegrad * (*dphi_coarse)[i][qp]) * JxW[qp];
freei++;
}
}
}
Ke.cholesky_solve(Fe, Uint);
// Transfer new interior solutions to element
for (unsigned int i=0; i != free_dofs; ++i)
{
Number &ui = Ue(free_dof[i]);
libmesh_assert(std::abs(ui) < TOLERANCE ||
std::abs(ui - Uint(i)) < TOLERANCE);
ui = Uint(i);
dof_is_fixed[free_dof[i]] = true;
}
// Make sure every DoF got reached!
for (unsigned int i=0; i != new_n_dofs; ++i)
libmesh_assert(dof_is_fixed[i]);
}
Definition at line 410 of file fe_map.C.
References compute_single_point_map(), d2xyzdeta2_map, d2xyzdetadzeta_map, d2xyzdxi2_map, d2xyzdxideta_map, d2xyzdxidzeta_map, d2xyzdzeta2_map, detadx_map, detady_map, detadz_map, dim, dxidx_map, dxidy_map, dxidz_map, dxyzdeta_map, dxyzdxi_map, dxyzdzeta_map, dzetadx_map, dzetady_map, dzetadz_map, JxW, phi_map, Elem::point(), resize_map_vectors(), and xyz.
Referenced by compute_map().
{
// Start logging the map computation.
START_LOG('compute_affine_map()', 'FE');
libmesh_assert (elem != NULL);
const unsigned int n_qp = qw.size();
// Resize the vectors to hold data at the quadrature points
this->resize_map_vectors(n_qp);
// Compute map at quadrature point 0
this->compute_single_point_map(qw, elem, 0);
// Compute xyz at all other quadrature points
for (unsigned int p=1; p<n_qp; p++)
{
xyz[p].zero();
for (unsigned int i=0; i<phi_map.size(); i++) // sum over the nodes
xyz[p].add_scaled (elem->point(i), phi_map[i][p] );
}
// Copy other map data from quadrature point 0
for (unsigned int p=1; p<n_qp; p++) // for each extra quadrature point
{
dxyzdxi_map[p] = dxyzdxi_map[0];
dxidx_map[p] = dxidx_map[0];
dxidy_map[p] = dxidy_map[0];
dxidz_map[p] = dxidz_map[0];
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
// The map should be affine, so second derivatives are zero
d2xyzdxi2_map[p] = 0.;
#endif
if (this->dim > 1)
{
dxyzdeta_map[p] = dxyzdeta_map[0];
detadx_map[p] = detadx_map[0];
detady_map[p] = detady_map[0];
detadz_map[p] = detadz_map[0];
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxideta_map[p] = 0.;
d2xyzdeta2_map[p] = 0.;
#endif
if (this->dim > 2)
{
dxyzdzeta_map[p] = dxyzdzeta_map[0];
dzetadx_map[p] = dzetadx_map[0];
dzetady_map[p] = dzetady_map[0];
dzetadz_map[p] = dzetadz_map[0];
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxidzeta_map[p] = 0.;
d2xyzdetadzeta_map[p] = 0.;
d2xyzdzeta2_map[p] = 0.;
#endif
}
}
JxW[p] = JxW[0] / qw[0] * qw[p];
}
STOP_LOG('compute_affine_map()', 'FE');
}
Definition at line 623 of file fe_boundary.C.
References compute_face_map(), curvatures, d2psidxi2_map, d2xyzdeta2_map, d2xyzdxi2_map, d2xyzdxideta_map, dim, dpsidxi_map, dxdxi_map(), dxyzdeta_map, dxyzdxi_map, dydxi_map(), dzdxi_map(), JxW, normals, Elem::point(), psi_map, tangents, and xyz.
Referenced by FE< Dim, T >::edge_reinit().
{
libmesh_assert (edge != NULL);
if (dim == 2)
{
// A 2D finite element living in either 2D or 3D space.
// The edges here are the sides of the element, so the
// (misnamed) compute_face_map function does what we want
FEBase::compute_face_map(qw, edge);
return;
}
libmesh_assert (dim == 3); // 1D is unnecessary and currently unsupported
START_LOG('compute_edge_map()', 'FE');
// The number of quadrature points.
const unsigned int n_qp = qw.size();
// Resize the vectors to hold data at the quadrature points
xyz.resize(n_qp);
dxyzdxi_map.resize(n_qp);
dxyzdeta_map.resize(n_qp);
d2xyzdxi2_map.resize(n_qp);
d2xyzdxideta_map.resize(n_qp);
d2xyzdeta2_map.resize(n_qp);
tangents.resize(n_qp);
normals.resize(n_qp);
curvatures.resize(n_qp);
JxW.resize(n_qp);
// Clear the entities that will be summed
for (unsigned int p=0; p<n_qp; p++)
{
tangents[p].resize(1);
xyz[p].zero();
dxyzdxi_map[p].zero();
dxyzdeta_map[p].zero();
d2xyzdxi2_map[p].zero();
d2xyzdxideta_map[p].zero();
d2xyzdeta2_map[p].zero();
}
// compute x, dxdxi at the quadrature points
for (unsigned int i=0; i<psi_map.size(); i++) // sum over the nodes
{
const Point& edge_point = edge->point(i);
for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
{
xyz[p].add_scaled (edge_point, psi_map[i][p]);
dxyzdxi_map[p].add_scaled (edge_point, dpsidxi_map[i][p]);
d2xyzdxi2_map[p].add_scaled (edge_point, d2psidxi2_map[i][p]);
}
}
// Compute the tangents at the quadrature point
// FIXME: normals (plural!) and curvatures are uncalculated
for (unsigned int p=0; p<n_qp; p++)
{
const Point n = dxyzdxi_map[p].cross(dxyzdeta_map[p]);
tangents[p][0] = dxyzdxi_map[p].unit();
// compute the jacobian at the quadrature points
const Real jac = std::sqrt(dxdxi_map(p)*dxdxi_map(p) +
dydxi_map(p)*dydxi_map(p) +
dzdxi_map(p)*dzdxi_map(p));
libmesh_assert (jac > 0.);
JxW[p] = jac*qw[p];
}
STOP_LOG('compute_edge_map()', 'FE');
}
Definition at line 345 of file fe_boundary.C.
References TypeVector< T >::cross(), curvatures, d2psideta2_map, d2psidxi2_map, d2psidxideta_map, d2xyzdeta2_map, d2xyzdxi2_map, d2xyzdxideta_map, dim, dpsideta_map, dpsidxi_map, dxdeta_map(), dxdxi_map(), dxyzdeta_map, dxyzdxi_map, dydeta_map(), dydxi_map(), dzdeta_map(), dzdxi_map(), InfFE< Dim, T_radial, T_map >::inverse_map(), JxW, Elem::node(), normals, Elem::parent(), Elem::point(), psi_map, tangents, TypeVector< T >::unit(), and xyz.
Referenced by compute_edge_map(), InfFE< Dim, T_radial, T_map >::reinit(), and REINIT_ERROR().
{
libmesh_assert (side != NULL);
START_LOG('compute_face_map()', 'FE');
// The number of quadrature points.
const unsigned int n_qp = qw.size();
switch (dim)
{
case 1:
{
// A 1D finite element, currently assumed to be in 1D space
// This means the boundary is a '0D finite element', a
// NODEELEM.
// Resize the vectors to hold data at the quadrature points
{
xyz.resize(n_qp);
normals.resize(n_qp);
JxW.resize(n_qp);
}
// If we have no quadrature points, there's nothing else to do
if (!n_qp)
break;
// We need to look back at the full edge to figure out the normal
// vector
const Elem *elem = side->parent();
libmesh_assert (elem);
if (side->node(0) == elem->node(0))
normals[0] = Point(-1.);
else
{
libmesh_assert (side->node(0) == elem->node(1));
normals[0] = Point(1.);
}
// Calculate x at the point
libmesh_assert (psi_map.size() == 1);
// In the unlikely event we have multiple quadrature
// points, they'll be in the same place
for (unsigned int p=0; p<n_qp; p++)
{
xyz[p].zero();
xyz[p].add_scaled (side->point(0), psi_map[0][p]);
normals[p] = normals[0];
JxW[p] = 1.0*qw[p];
}
// done computing the map
break;
}
case 2:
{
// A 2D finite element living in either 2D or 3D space.
// This means the boundary is a 1D finite element, i.e.
// and EDGE2 or EDGE3.
// Resize the vectors to hold data at the quadrature points
{
xyz.resize(n_qp);
dxyzdxi_map.resize(n_qp);
d2xyzdxi2_map.resize(n_qp);
tangents.resize(n_qp);
normals.resize(n_qp);
curvatures.resize(n_qp);
JxW.resize(n_qp);
}
// Clear the entities that will be summed
// Compute the tangent & normal at the quadrature point
for (unsigned int p=0; p<n_qp; p++)
{
tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
xyz[p].zero();
dxyzdxi_map[p].zero();
d2xyzdxi2_map[p].zero();
}
// compute x, dxdxi at the quadrature points
for (unsigned int i=0; i<psi_map.size(); i++) // sum over the nodes
{
const Point& side_point = side->point(i);
for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
{
xyz[p].add_scaled (side_point, psi_map[i][p]);
dxyzdxi_map[p].add_scaled (side_point, dpsidxi_map[i][p]);
d2xyzdxi2_map[p].add_scaled(side_point, d2psidxi2_map[i][p]);
}
}
// Compute the tangent & normal at the quadrature point
for (unsigned int p=0; p<n_qp; p++)
{
// The first tangent comes from just the edge's Jacobian
tangents[p][0] = dxyzdxi_map[p].unit();
#if LIBMESH_DIM == 2
// For a 2D element living in 2D, the normal is given directly
// from the entries in the edge Jacobian.
normals[p] = (Point(dxyzdxi_map[p](1), -dxyzdxi_map[p](0), 0.)).unit();
#elif LIBMESH_DIM == 3
// For a 2D element living in 3D, there is a second tangent.
// For the second tangent, we need to refer to the full
// element's (not just the edge's) Jacobian.
const Elem *elem = side->parent();
libmesh_assert (elem != NULL);
// Inverse map xyz[p] to a reference point on the parent...
Point reference_point = FE<2,LAGRANGE>::inverse_map(elem, xyz[p]);
// Get dxyz/dxi and dxyz/deta from the parent map.
Point dx_dxi = FE<2,LAGRANGE>::map_xi (elem, reference_point);
Point dx_deta = FE<2,LAGRANGE>::map_eta(elem, reference_point);
// The second tangent vector is formed by crossing these vectors.
tangents[p][1] = dx_dxi.cross(dx_deta).unit();
// Finally, the normal in this case is given by crossing these
// two tangents.
normals[p] = tangents[p][0].cross(tangents[p][1]).unit();
#endif
// The curvature is computed via the familiar Frenet formula:
// curvature = [d^2(x) / d (xi)^2] dot [normal]
// For a reference, see:
// F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill, p. 310
//
// Note: The sign convention here is different from the
// 3D case. Concave-upward curves (smiles) have a positive
// curvature. Concave-downward curves (frowns) have a
// negative curvature. Be sure to take that into account!
const Real numerator = d2xyzdxi2_map[p] * normals[p];
const Real denominator = dxyzdxi_map[p].size_sq();
libmesh_assert (denominator != 0);
curvatures[p] = numerator / denominator;
}
// compute the jacobian at the quadrature points
for (unsigned int p=0; p<n_qp; p++)
{
const Real jac = dxyzdxi_map[p].size();
libmesh_assert (jac > 0.);
JxW[p] = jac*qw[p];
}
// done computing the map
break;
}
case 3:
{
// A 3D finite element living in 3D space.
// Resize the vectors to hold data at the quadrature points
{
xyz.resize(n_qp);
dxyzdxi_map.resize(n_qp);
dxyzdeta_map.resize(n_qp);
d2xyzdxi2_map.resize(n_qp);
d2xyzdxideta_map.resize(n_qp);
d2xyzdeta2_map.resize(n_qp);
tangents.resize(n_qp);
normals.resize(n_qp);
curvatures.resize(n_qp);
JxW.resize(n_qp);
}
// Clear the entities that will be summed
for (unsigned int p=0; p<n_qp; p++)
{
tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
xyz[p].zero();
dxyzdxi_map[p].zero();
dxyzdeta_map[p].zero();
d2xyzdxi2_map[p].zero();
d2xyzdxideta_map[p].zero();
d2xyzdeta2_map[p].zero();
}
// compute x, dxdxi at the quadrature points
for (unsigned int i=0; i<psi_map.size(); i++) // sum over the nodes
{
const Point& side_point = side->point(i);
for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
{
xyz[p].add_scaled (side_point, psi_map[i][p]);
dxyzdxi_map[p].add_scaled (side_point, dpsidxi_map[i][p]);
dxyzdeta_map[p].add_scaled(side_point, dpsideta_map[i][p]);
d2xyzdxi2_map[p].add_scaled (side_point, d2psidxi2_map[i][p]);
d2xyzdxideta_map[p].add_scaled(side_point, d2psidxideta_map[i][p]);
d2xyzdeta2_map[p].add_scaled (side_point, d2psideta2_map[i][p]);
}
}
// Compute the tangents, normal, and curvature at the quadrature point
for (unsigned int p=0; p<n_qp; p++)
{
const Point n = dxyzdxi_map[p].cross(dxyzdeta_map[p]);
normals[p] = n.unit();
tangents[p][0] = dxyzdxi_map[p].unit();
tangents[p][1] = n.cross(dxyzdxi_map[p]).unit();
// Compute curvature using the typical nomenclature
// of the first and second fundamental forms.
// For reference, see:
// 1) http://mathworld.wolfram.com/MeanCurvature.html
// (note -- they are using inward normal)
// 2) F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill
const Real L = -d2xyzdxi2_map[p] * normals[p];
const Real M = -d2xyzdxideta_map[p] * normals[p];
const Real N = -d2xyzdeta2_map[p] * normals[p];
const Real E = dxyzdxi_map[p].size_sq();
const Real F = dxyzdxi_map[p] * dxyzdeta_map[p];
const Real G = dxyzdeta_map[p].size_sq();
const Real numerator = E*N -2.*F*M + G*L;
const Real denominator = E*G - F*F;
libmesh_assert (denominator != 0.);
curvatures[p] = 0.5*numerator/denominator;
}
// compute the jacobian at the quadrature points, see
// http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/thtoc.htm
for (unsigned int p=0; p<n_qp; p++)
{
const Real g11 = (dxdxi_map(p)*dxdxi_map(p) +
dydxi_map(p)*dydxi_map(p) +
dzdxi_map(p)*dzdxi_map(p));
const Real g12 = (dxdxi_map(p)*dxdeta_map(p) +
dydxi_map(p)*dydeta_map(p) +
dzdxi_map(p)*dzdeta_map(p));
const Real g21 = g12;
const Real g22 = (dxdeta_map(p)*dxdeta_map(p) +
dydeta_map(p)*dydeta_map(p) +
dzdeta_map(p)*dzdeta_map(p));
const Real jac = std::sqrt(g11*g22 - g12*g21);
libmesh_assert (jac > 0.);
JxW[p] = jac*qw[p];
}
// done computing the map
break;
}
default:
libmesh_error();
}
STOP_LOG('compute_face_map()', 'FE');
}
Definition at line 476 of file fe_map.C.
References calculate_d2phi, compute_affine_map(), compute_single_point_map(), Elem::has_affine_map(), and resize_map_vectors().
{
if (elem->has_affine_map())
{
compute_affine_map(qw, elem);
return;
}
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
static bool curvy_second_derivative_warning = false;
if (calculate_d2phi && !curvy_second_derivative_warning)
{
std::cerr << 'WARNING: Second derivatives are not currently '
<< 'correctly calculated on non-affine elements!'
<< std::endl;
curvy_second_derivative_warning = true;
}
#endif
// Start logging the map computation.
START_LOG('compute_map()', 'FE');
libmesh_assert (elem != NULL);
const unsigned int n_qp = qw.size();
// Resize the vectors to hold data at the quadrature points
this->resize_map_vectors(n_qp);
// Compute map at all quadrature points
for (unsigned int p=0; p!=n_qp; p++)
this->compute_single_point_map(qw, elem, p);
// Stop logging the map computation.
STOP_LOG('compute_map()', 'FE');
}
Definition at line 1900 of file fe_base.C.
References Elem::active(), PeriodicBoundaries::boundary(), MeshBase::boundary_info, build(), libMeshEnums::C_ONE, libMeshEnums::C_ZERO, DenseMatrix< T >::cholesky_solve(), DofMap::constrain_p_dofs(), FEType::default_quadrature_order(), Elem::dim(), libMeshEnums::DISCONTINUOUS, DofMap::dof_indices(), FEInterface::dofs_on_side(), dphi, DofObject::invalid_id, libMesh::invalid_uint, FEInterface::inverse_map(), DofMap::is_constrained_dof(), JxW, Elem::level(), std::min(), Elem::min_p_level_by_neighbor(), Elem::n_sides(), PeriodicBoundaries::neighbor(), Elem::neighbor(), Elem::p_level(), PeriodicBoundary::pairedboundary, phi, DenseVector< T >::resize(), DenseMatrix< T >::resize(), Threads::spin_mtx, PeriodicBoundary::translation_vector, and DofMap::variable_type().
{
// Only bother if we truly have periodic boundaries
if (boundaries.empty())
return;
libmesh_assert (elem != NULL);
// Only constrain active elements with this method
if (!elem->active())
return;
const unsigned int Dim = elem->dim();
const FEType& base_fe_type = dof_map.variable_type(variable_number);
// Construct FE objects for this element and its pseudo-neighbors.
AutoPtr<FEBase> my_fe (FEBase::build(Dim, base_fe_type));
const FEContinuity cont = my_fe->get_continuity();
// We don't need to constrain discontinuous elements
if (cont == DISCONTINUOUS)
return;
libmesh_assert (cont == C_ZERO || cont == C_ONE);
AutoPtr<FEBase> neigh_fe (FEBase::build(Dim, base_fe_type));
QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
my_fe->attach_quadrature_rule (&my_qface);
std::vector<Point> neigh_qface;
const std::vector<Real>& JxW = my_fe->get_JxW();
const std::vector<Point>& q_point = my_fe->get_xyz();
const std::vector<std::vector<Real> >& phi = my_fe->get_phi();
const std::vector<std::vector<Real> >& neigh_phi =
neigh_fe->get_phi();
const std::vector<Point> *face_normals = NULL;
const std::vector<std::vector<RealGradient> > *dphi = NULL;
const std::vector<std::vector<RealGradient> > *neigh_dphi = NULL;
std::vector<unsigned int> my_dof_indices, neigh_dof_indices;
std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
if (cont != C_ZERO)
{
const std::vector<Point>& ref_face_normals =
my_fe->get_normals();
face_normals = &ref_face_normals;
const std::vector<std::vector<RealGradient> >& ref_dphi =
my_fe->get_dphi();
dphi = &ref_dphi;
const std::vector<std::vector<RealGradient> >& ref_neigh_dphi =
neigh_fe->get_dphi();
neigh_dphi = &ref_neigh_dphi;
}
DenseMatrix<Real> Ke;
DenseVector<Real> Fe;
std::vector<DenseVector<Real> > Ue;
// Look at the element faces. Check to see if we need to
// build constraints.
for (unsigned int s=0; s<elem->n_sides(); s++)
{
if (elem->neighbor(s))
continue;
unsigned int boundary_id = mesh.boundary_info->boundary_id(elem, s);
PeriodicBoundary *periodic = boundaries.boundary(boundary_id);
if (periodic)
{
// Get pointers to the element's neighbor.
const Elem* neigh = boundaries.neighbor(boundary_id, mesh, elem, s);
// h refinement constraints:
// constrain dofs shared between
// this element and ones as coarse
// as or coarser than this element.
if (neigh->level() <= elem->level())
{
unsigned int s_neigh =
mesh.boundary_info->side_with_boundary_id (neigh, periodic->pairedboundary);
libmesh_assert(s_neigh != libMesh::invalid_uint);
#ifdef LIBMESH_ENABLE_AMR
// Find the minimum p level; we build the h constraint
// matrix with this and then constrain away all higher p
// DoFs.
libmesh_assert(neigh->active());
const unsigned int min_p_level =
std::min(elem->p_level(), neigh->p_level());
// we may need to make the FE objects reinit with the
// minimum shared p_level
// FIXME - I hate using const_cast<> and avoiding
// accessor functions; there's got to be a
// better way to do this!
const unsigned int old_elem_level = elem->p_level();
if (old_elem_level != min_p_level)
(const_cast<Elem *>(elem))->hack_p_level(min_p_level);
const unsigned int old_neigh_level = neigh->p_level();
if (old_neigh_level != min_p_level)
(const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
#endif // #ifdef LIBMESH_ENABLE_AMR
my_fe->reinit(elem, s);
dof_map.dof_indices (elem, my_dof_indices,
variable_number);
dof_map.dof_indices (neigh, neigh_dof_indices,
variable_number);
const unsigned int n_qp = my_qface.n_points();
// Translate the quadrature points over to the
// neighbor's boundary
std::vector<Point> neigh_point = q_point;
for (unsigned int i=0; i != neigh_point.size(); ++i)
neigh_point[i] += periodic->translation_vector;
FEInterface::inverse_map (Dim, base_fe_type, neigh,
neigh_point, neigh_qface);
neigh_fe->reinit(neigh, &neigh_qface);
// We're only concerned with DOFs whose values (and/or first
// derivatives for C1 elements) are supported on side nodes
FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
// We're done with functions that examine Elem::p_level(),
// so let's unhack those levels
#ifdef LIBMESH_ENABLE_AMR
if (elem->p_level() != old_elem_level)
(const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
if (neigh->p_level() != old_neigh_level)
(const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
#endif // #ifdef LIBMESH_ENABLE_AMR
const unsigned int n_side_dofs = my_side_dofs.size();
libmesh_assert(n_side_dofs == neigh_side_dofs.size());
Ke.resize (n_side_dofs, n_side_dofs);
Ue.resize(n_side_dofs);
// Form the projection matrix, (inner product of fine basis
// functions against fine test functions)
for (unsigned int is = 0; is != n_side_dofs; ++is)
{
const unsigned int i = my_side_dofs[is];
for (unsigned int js = 0; js != n_side_dofs; ++js)
{
const unsigned int j = my_side_dofs[js];
for (unsigned int qp = 0; qp != n_qp; ++qp)
{
Ke(is,js) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
if (cont != C_ZERO)
Ke(is,js) += JxW[qp] * (((*dphi)[i][qp] *
(*face_normals)[qp]) *
((*dphi)[j][qp] *
(*face_normals)[qp]));
}
}
}
// Form the right hand sides, (inner product of coarse basis
// functions against fine test functions)
for (unsigned int is = 0; is != n_side_dofs; ++is)
{
const unsigned int i = neigh_side_dofs[is];
Fe.resize (n_side_dofs);
for (unsigned int js = 0; js != n_side_dofs; ++js)
{
const unsigned int j = my_side_dofs[js];
for (unsigned int qp = 0; qp != n_qp; ++qp)
{
Fe(js) += JxW[qp] * (neigh_phi[i][qp] *
phi[j][qp]);
if (cont != C_ZERO)
Fe(js) += JxW[qp] * (((*neigh_dphi)[i][qp] *
(*face_normals)[qp]) *
((*dphi)[j][qp] *
(*face_normals)[qp]));
}
}
Ke.cholesky_solve(Fe, Ue[is]);
}
// Make sure we're not adding recursive constraints
// due to the redundancy in the way we add periodic
// boundary constraints
std::vector<bool> recursive_constraint(n_side_dofs, false);
for (unsigned int is = 0; is != n_side_dofs; ++is)
{
const unsigned int i = neigh_side_dofs[is];
const unsigned int their_dof_g = neigh_dof_indices[i];
libmesh_assert(their_dof_g != DofObject::invalid_id);
if (!dof_map.is_constrained_dof(their_dof_g))
continue;
DofConstraintRow& their_constraint_row =
constraints[their_dof_g];
for (unsigned int js = 0; js != n_side_dofs; ++js)
{
const unsigned int j = my_side_dofs[js];
const unsigned int my_dof_g = my_dof_indices[j];
libmesh_assert(my_dof_g != DofObject::invalid_id);
if (their_constraint_row.count(my_dof_g))
recursive_constraint[js] = true;
}
}
for (unsigned int is = 0; is != n_side_dofs; ++is)
{
const unsigned int i = neigh_side_dofs[is];
const unsigned int their_dof_g = neigh_dof_indices[i];
libmesh_assert(their_dof_g != DofObject::invalid_id);
for (unsigned int js = 0; js != n_side_dofs; ++js)
{
if (recursive_constraint[js])
continue;
const unsigned int j = my_side_dofs[js];
const unsigned int my_dof_g = my_dof_indices[j];
libmesh_assert(my_dof_g != DofObject::invalid_id);
if (dof_map.is_constrained_dof(my_dof_g))
continue;
const Real their_dof_value = Ue[is](js);
if (their_dof_g == my_dof_g)
{
libmesh_assert(std::abs(their_dof_value-1.) < 1.e-5);
for (unsigned int k = 0; k != n_side_dofs; ++k)
libmesh_assert(k == is || std::abs(Ue[k](js)) < 1.e-5);
continue;
}
if (std::abs(their_dof_value) < 1.e-5)
continue;
// since we may be running this method concurretly
// on multiple threads we need to acquire a lock
// before modifying the shared constraint_row object.
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
DofConstraintRow& constraint_row =
constraints[my_dof_g];
constraint_row.insert(std::make_pair(their_dof_g,
their_dof_value));
}
}
}
}
// p refinement constraints:
// constrain dofs shared between
// active elements and neighbors with
// lower polynomial degrees
#ifdef LIBMESH_ENABLE_AMR
const unsigned int min_p_level =
neigh->min_p_level_by_neighbor(elem, elem->p_level());
if (min_p_level < elem->p_level())
{
// Adaptive p refinement of non-hierarchic bases will
// require more coding
libmesh_assert(my_fe->is_hierarchic());
dof_map.constrain_p_dofs(variable_number, elem,
s, min_p_level);
}
#endif // #ifdef LIBMESH_ENABLE_AMR
}
}
}
Definition at line 1659 of file fe_base.C.
References Elem::active(), build(), libMeshEnums::C_ONE, libMeshEnums::C_ZERO, DenseMatrix< T >::cholesky_solve(), DofMap::constrain_p_dofs(), FEType::default_quadrature_order(), Elem::dim(), libMeshEnums::DISCONTINUOUS, DofMap::dof_indices(), FEInterface::dofs_on_side(), dphi, DofObject::invalid_id, FEInterface::inverse_map(), JxW, Elem::level(), std::min(), Elem::min_p_level_by_neighbor(), Elem::n_nodes(), Elem::n_sides(), Elem::neighbor(), Elem::p_level(), phi, DenseVector< T >::resize(), DenseMatrix< T >::resize(), Threads::spin_mtx, DofMap::variable_type(), and Elem::which_neighbor_am_i().
Referenced by FE< Dim, T >::compute_constraints().
{
libmesh_assert (elem != NULL);
const unsigned int Dim = elem->dim();
// Only constrain elements in 2,3D.
if (Dim == 1)
return;
// Only constrain active elements with this method
if (!elem->active())
return;
const FEType& base_fe_type = dof_map.variable_type(variable_number);
// Construct FE objects for this element and its neighbors.
AutoPtr<FEBase> my_fe (FEBase::build(Dim, base_fe_type));
const FEContinuity cont = my_fe->get_continuity();
// We don't need to constrain discontinuous elements
if (cont == DISCONTINUOUS)
return;
libmesh_assert (cont == C_ZERO || cont == C_ONE);
AutoPtr<FEBase> neigh_fe (FEBase::build(Dim, base_fe_type));
QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
my_fe->attach_quadrature_rule (&my_qface);
std::vector<Point> neigh_qface;
const std::vector<Real>& JxW = my_fe->get_JxW();
const std::vector<Point>& q_point = my_fe->get_xyz();
const std::vector<std::vector<Real> >& phi = my_fe->get_phi();
const std::vector<std::vector<Real> >& neigh_phi =
neigh_fe->get_phi();
const std::vector<Point> *face_normals = NULL;
const std::vector<std::vector<RealGradient> > *dphi = NULL;
const std::vector<std::vector<RealGradient> > *neigh_dphi = NULL;
std::vector<unsigned int> my_dof_indices, neigh_dof_indices;
std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
if (cont != C_ZERO)
{
const std::vector<Point>& ref_face_normals =
my_fe->get_normals();
face_normals = &ref_face_normals;
const std::vector<std::vector<RealGradient> >& ref_dphi =
my_fe->get_dphi();
dphi = &ref_dphi;
const std::vector<std::vector<RealGradient> >& ref_neigh_dphi =
neigh_fe->get_dphi();
neigh_dphi = &ref_neigh_dphi;
}
DenseMatrix<Real> Ke;
DenseVector<Real> Fe;
std::vector<DenseVector<Real> > Ue;
// Look at the element faces. Check to see if we need to
// build constraints.
for (unsigned int s=0; s<elem->n_sides(); s++)
if (elem->neighbor(s) != NULL)
{
// Get pointers to the element's neighbor.
const Elem* neigh = elem->neighbor(s);
// h refinement constraints:
// constrain dofs shared between
// this element and ones coarser
// than this element.
if (neigh->level() < elem->level())
{
unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
libmesh_assert (s_neigh < neigh->n_neighbors());
// Find the minimum p level; we build the h constraint
// matrix with this and then constrain away all higher p
// DoFs.
libmesh_assert(neigh->active());
const unsigned int min_p_level =
std::min(elem->p_level(), neigh->p_level());
// we may need to make the FE objects reinit with the
// minimum shared p_level
// FIXME - I hate using const_cast<> and avoiding
// accessor functions; there's got to be a
// better way to do this!
const unsigned int old_elem_level = elem->p_level();
if (old_elem_level != min_p_level)
(const_cast<Elem *>(elem))->hack_p_level(min_p_level);
const unsigned int old_neigh_level = neigh->p_level();
if (old_neigh_level != min_p_level)
(const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
my_fe->reinit(elem, s);
// This function gets called element-by-element, so there
// will be a lot of memory allocation going on. We can
// at least minimize this for the case of the dof indices
// by efficiently preallocating the requisite storage.
// n_nodes is not necessarily n_dofs, but it is better
// than nothing!
my_dof_indices.reserve (elem->n_nodes());
neigh_dof_indices.reserve (neigh->n_nodes());
dof_map.dof_indices (elem, my_dof_indices,
variable_number);
dof_map.dof_indices (neigh, neigh_dof_indices,
variable_number);
const unsigned int n_qp = my_qface.n_points();
FEInterface::inverse_map (Dim, base_fe_type, neigh,
q_point, neigh_qface);
neigh_fe->reinit(neigh, &neigh_qface);
// We're only concerned with DOFs whose values (and/or first
// derivatives for C1 elements) are supported on side nodes
FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
// We're done with functions that examine Elem::p_level(),
// so let's unhack those levels
if (elem->p_level() != old_elem_level)
(const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
if (neigh->p_level() != old_neigh_level)
(const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
const unsigned int n_side_dofs = my_side_dofs.size();
libmesh_assert(n_side_dofs == neigh_side_dofs.size());
Ke.resize (n_side_dofs, n_side_dofs);
Ue.resize(n_side_dofs);
// Form the projection matrix, (inner product of fine basis
// functions against fine test functions)
for (unsigned int is = 0; is != n_side_dofs; ++is)
{
const unsigned int i = my_side_dofs[is];
for (unsigned int js = 0; js != n_side_dofs; ++js)
{
const unsigned int j = my_side_dofs[js];
for (unsigned int qp = 0; qp != n_qp; ++qp)
{
Ke(is,js) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
if (cont != C_ZERO)
Ke(is,js) += JxW[qp] * (((*dphi)[i][qp] *
(*face_normals)[qp]) *
((*dphi)[j][qp] *
(*face_normals)[qp]));
}
}
}
// Form the right hand sides, (inner product of coarse basis
// functions against fine test functions)
for (unsigned int is = 0; is != n_side_dofs; ++is)
{
const unsigned int i = neigh_side_dofs[is];
Fe.resize (n_side_dofs);
for (unsigned int js = 0; js != n_side_dofs; ++js)
{
const unsigned int j = my_side_dofs[js];
for (unsigned int qp = 0; qp != n_qp; ++qp)
{
Fe(js) += JxW[qp] * (neigh_phi[i][qp] *
phi[j][qp]);
if (cont != C_ZERO)
Fe(js) += JxW[qp] * (((*neigh_dphi)[i][qp] *
(*face_normals)[qp]) *
((*dphi)[j][qp] *
(*face_normals)[qp]));
}
}
Ke.cholesky_solve(Fe, Ue[is]);
}
for (unsigned int is = 0; is != n_side_dofs; ++is)
{
const unsigned int i = neigh_side_dofs[is];
const unsigned int their_dof_g = neigh_dof_indices[i];
libmesh_assert(their_dof_g != DofObject::invalid_id);
for (unsigned int js = 0; js != n_side_dofs; ++js)
{
const unsigned int j = my_side_dofs[js];
const unsigned int my_dof_g = my_dof_indices[j];
libmesh_assert(my_dof_g != DofObject::invalid_id);
const Real their_dof_value = Ue[is](js);
if (their_dof_g == my_dof_g)
{
libmesh_assert(std::abs(their_dof_value-1.) < 1.e-5);
for (unsigned int k = 0; k != n_side_dofs; ++k)
libmesh_assert(k == is || std::abs(Ue[k](js)) < 1.e-5);
continue;
}
if (std::abs(their_dof_value) < 1.e-5)
continue;
// since we may be running this method concurretly
// on multiple threads we need to acquire a lock
// before modifying the shared constraint_row object.
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
DofConstraintRow& constraint_row =
constraints[my_dof_g];
constraint_row.insert(std::make_pair(their_dof_g,
their_dof_value));
}
}
}
}
// p refinement constraints:
// constrain dofs shared between
// active elements and neighbors with
// lower polynomial degrees
const unsigned int min_p_level =
neigh->min_p_level_by_neighbor(elem, elem->p_level());
if (min_p_level < elem->p_level())
{
// Adaptive p refinement of non-hierarchic bases will
// require more coding
libmesh_assert(my_fe->is_hierarchic());
dof_map.constrain_p_dofs(variable_number, elem,
s, min_p_level);
}
}
}
Reimplemented in FEXYZ< Dim >, and InfFE< Dim, T_radial, T_map >.
Definition at line 632 of file fe_base.C.
References calculate_d2phi, calculate_dphi, calculate_phi, calculations_started, d2phi, d2phideta2, d2phidetadzeta, d2phidx2, d2phidxdy, d2phidxdz, d2phidxi2, d2phidxideta, d2phidxidzeta, d2phidy2, d2phidydz, d2phidz2, d2phidzeta2, detadx_map, detady_map, detadz_map, dim, dphi, dphideta, dphidx, dphidxi, dphidy, dphidz, dphidzeta, dxidx_map, dxidy_map, dxidz_map, dzetadx_map, dzetady_map, and dzetadz_map.
{
//-------------------------------------------------------------------------
// Compute the shape function values (and derivatives)
// at the Quadrature points. Note that the actual values
// have already been computed via init_shape_functions
// Start logging the shape function computation
START_LOG('compute_shape_functions()', 'FE');
calculations_started = true;
// If the user forgot to request anything, we'll be safe and
// calculate everything:
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
if (!calculate_phi && !calculate_dphi && !calculate_d2phi)
calculate_phi = calculate_dphi = calculate_d2phi = true;
#else
if (!calculate_phi && !calculate_dphi)
calculate_phi = calculate_dphi = true;
#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
// Compute the value of the derivative shape function i at quadrature point p
switch (dim)
{
case 1:
{
if (calculate_dphi)
for (unsigned int i=0; i<dphi.size(); i++)
for (unsigned int p=0; p<dphi[i].size(); p++)
{
// dphi/dx = (dphi/dxi)*(dxi/dx)
dphi[i][p](0) =
dphidx[i][p] = dphidxi[i][p]*dxidx_map[p];
dphi[i][p](1) = dphidy[i][p] = 0.;
dphi[i][p](2) = dphidz[i][p] = 0.;
}
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
if (calculate_d2phi)
for (unsigned int i=0; i<d2phi.size(); i++)
for (unsigned int p=0; p<d2phi[i].size(); p++)
{
d2phi[i][p](0,0) = d2phidx2[i][p] =
d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p];
#if LIBMESH_DIM>1
d2phi[i][p](0,1) = d2phidxdy[i][p] =
d2phi[i][p](1,0) = 0.;
d2phi[i][p](1,1) = d2phidy2[i][p] = 0.;
#if LIBMESH_DIM>2
d2phi[i][p](0,2) = d2phidxdz[i][p] =
d2phi[i][p](2,0) = 0.;
d2phi[i][p](1,2) = d2phidydz[i][p] =
d2phi[i][p](2,1) = 0.;
d2phi[i][p](2,2) = d2phidz2[i][p] = 0.;
#endif
#endif
}
#endif
// All done
break;
}
case 2:
{
if (calculate_dphi)
for (unsigned int i=0; i<dphi.size(); i++)
for (unsigned int p=0; p<dphi[i].size(); p++)
{
// dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx)
dphi[i][p](0) =
dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
dphideta[i][p]*detadx_map[p]);
// dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy)
dphi[i][p](1) =
dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
dphideta[i][p]*detady_map[p]);
// dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz)
#if LIBMESH_DIM == 3
dphi[i][p](2) = // can only assign to the Z component if LIBMESH_DIM==3
#endif
dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
dphideta[i][p]*detadz_map[p]);
}
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
if (calculate_d2phi)
for (unsigned int i=0; i<d2phi.size(); i++)
for (unsigned int p=0; p<d2phi[i].size(); p++)
{
d2phi[i][p](0,0) = d2phidx2[i][p] =
d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +
d2phideta2[i][p]*detadx_map[p]*detadx_map[p];
d2phi[i][p](0,1) = d2phidxdy[i][p] =
d2phi[i][p](1,0) =
d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +
d2phidxideta[i][p]*dxidx_map[p]*detady_map[p] +
d2phideta2[i][p]*detadx_map[p]*detady_map[p] +
d2phidxideta[i][p]*detadx_map[p]*dxidy_map[p];
d2phi[i][p](1,1) = d2phidy2[i][p] =
d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +
d2phideta2[i][p]*detady_map[p]*detady_map[p];
#if LIBMESH_DIM == 3
d2phi[i][p](0,2) = d2phidxdz[i][p] =
d2phi[i][p](2,0) =
d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +
d2phidxideta[i][p]*dxidx_map[p]*detadz_map[p] +
d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +
d2phidxideta[i][p]*detadx_map[p]*dxidz_map[p];
d2phi[i][p](1,2) = d2phidydz[i][p] =
d2phi[i][p](2,1) =
d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +
d2phidxideta[i][p]*dxidy_map[p]*detadz_map[p] +
d2phideta2[i][p]*detady_map[p]*detadz_map[p] +
d2phidxideta[i][p]*detady_map[p]*dxidz_map[p];
d2phi[i][p](2,2) = d2phidz2[i][p] =
d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +
d2phideta2[i][p]*detadz_map[p]*detadz_map[p];
#endif
}
#endif
// All done
break;
}
case 3:
{
if (calculate_dphi)
for (unsigned int i=0; i<dphi.size(); i++)
for (unsigned int p=0; p<dphi[i].size(); p++)
{
// dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
dphi[i][p](0) =
dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
dphideta[i][p]*detadx_map[p] +
dphidzeta[i][p]*dzetadx_map[p]);
// dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
dphi[i][p](1) =
dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
dphideta[i][p]*detady_map[p] +
dphidzeta[i][p]*dzetady_map[p]);
// dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
dphi[i][p](2) =
dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
dphideta[i][p]*detadz_map[p] +
dphidzeta[i][p]*dzetadz_map[p]);
}
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
if (calculate_d2phi)
for (unsigned int i=0; i<d2phi.size(); i++)
for (unsigned int p=0; p<d2phi[i].size(); p++)
{
d2phi[i][p](0,0) = d2phidx2[i][p] =
d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +
2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] +
2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] +
d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +
d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p];
d2phi[i][p](0,1) = d2phidxdy[i][p] =
d2phi[i][p](1,0) =
d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +
d2phidxideta[i][p]*dxidx_map[p]*detady_map[p] +
d2phidxidzeta[i][p]*dxidx_map[p]*dzetady_map[p] +
d2phideta2[i][p]*detadx_map[p]*detady_map[p] +
d2phidxideta[i][p]*detadx_map[p]*dxidy_map[p] +
d2phidetadzeta[i][p]*detadx_map[p]*dzetady_map[p] +
d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] +
d2phidxidzeta[i][p]*dzetadx_map[p]*dxidy_map[p] +
d2phidetadzeta[i][p]*dzetadx_map[p]*detady_map[p];
d2phi[i][p](0,2) = d2phidxdz[i][p] =
d2phi[i][p](2,0) =
d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +
d2phidxideta[i][p]*dxidx_map[p]*detadz_map[p] +
d2phidxidzeta[i][p]*dxidx_map[p]*dzetadz_map[p] +
d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +
d2phidxideta[i][p]*detadx_map[p]*dxidz_map[p] +
d2phidetadzeta[i][p]*detadx_map[p]*dzetadz_map[p] +
d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] +
d2phidxidzeta[i][p]*dzetadx_map[p]*dxidz_map[p] +
d2phidetadzeta[i][p]*dzetadx_map[p]*detadz_map[p];
d2phi[i][p](1,1) = d2phidy2[i][p] =
d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +
2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] +
2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] +
d2phideta2[i][p]*detady_map[p]*detady_map[p] +
d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p];
d2phi[i][p](1,2) = d2phidydz[i][p] =
d2phi[i][p](2,1) =
d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +
d2phidxideta[i][p]*dxidy_map[p]*detadz_map[p] +
d2phidxidzeta[i][p]*dxidy_map[p]*dzetadz_map[p] +
d2phideta2[i][p]*detady_map[p]*detadz_map[p] +
d2phidxideta[i][p]*detady_map[p]*dxidz_map[p] +
d2phidetadzeta[i][p]*detady_map[p]*dzetadz_map[p] +
d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] +
d2phidxidzeta[i][p]*dzetady_map[p]*dxidz_map[p] +
d2phidetadzeta[i][p]*dzetady_map[p]*detadz_map[p];
d2phi[i][p](2,2) = d2phidz2[i][p] =
d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +
2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] +
2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] +
d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +
d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p];
}
#endif
// All done
break;
}
default:
{
libmesh_error();
}
}
// Stop logging the shape function computation
STOP_LOG('compute_shape_functions()', 'FE');
}
Definition at line 35 of file fe_map.C.
References d2phideta2_map, d2phidetadzeta_map, d2phidxi2_map, d2phidxideta_map, d2phidxidzeta_map, d2phidzeta2_map, d2xyzdeta2_map, d2xyzdetadzeta_map, d2xyzdxi2_map, d2xyzdxideta_map, d2xyzdxidzeta_map, d2xyzdzeta2_map, detadx_map, detady_map, detadz_map, dim, dphideta_map, dphidxi_map, dphidzeta_map, dxdeta_map(), dxdxi_map(), dxdzeta_map(), dxidx_map, dxidy_map, dxidz_map, dxyzdeta_map, dxyzdxi_map, dxyzdzeta_map, dydeta_map(), dydxi_map(), dydzeta_map(), dzdeta_map(), dzdxi_map(), dzdzeta_map(), dzetadx_map, dzetady_map, dzetadz_map, DofObject::id(), JxW, phi_map, Elem::point(), and xyz.
Referenced by compute_affine_map(), and compute_map().
{
libmesh_assert (elem != NULL);
switch (this->dim)
{
//--------------------------------------------------------------------
// 1D
case 1:
{
// Clear the entities that will be summed
xyz[p].zero();
dxyzdxi_map[p].zero();
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxi2_map[p].zero();
#endif
// compute x, dx, d2x at the quadrature point
for (unsigned int i=0; i<phi_map.size(); i++) // sum over the nodes
{
// Reference to the point, helps eliminate
// exessive temporaries in the inner loop
const Point& elem_point = elem->point(i);
xyz[p].add_scaled (elem_point, phi_map[i][p] );
dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p]);
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxi2_map[p].add_scaled(elem_point, d2phidxi2_map[i][p]);
#endif
}
// Compute the jacobian
//
// 1D elements can live in 2D or 3D space.
// The transformation matrix from local->global
// coordinates is
//
// T = | dx/dxi |
// | dy/dxi |
// | dz/dxi |
//
// The generalized determinant of T (from the
// so-called 'normal' eqns.) is
// jac = 'det(T)' = sqrt(det(T'T))
//
// where T'= transpose of T, so
//
// jac = sqrt( (dx/dxi)^2 + (dy/dxi)^2 + (dz/dxi)^2 )
const Real jac = dxyzdxi_map[p].size();
if (jac <= 0.)
{
std::cerr << 'ERROR: negative Jacobian: '
<< jac
<< ' in element '
<< elem->id()
<< std::endl;
libmesh_error();
}
// The inverse Jacobian entries also come from the
// generalized inverse of T (see also the 2D element
// living in 3D code).
const Real jacm2 = 1./jac/jac;
dxidx_map[p] = jacm2*dxdxi_map(p);
dxidy_map[p] = jacm2*dydxi_map(p);
dxidz_map[p] = jacm2*dzdxi_map(p);
JxW[p] = jac*qw[p];
// done computing the map
break;
}
//--------------------------------------------------------------------
// 2D
case 2:
{
//------------------------------------------------------------------
// Compute the (x,y) values at the quadrature points,
// the Jacobian at the quadrature points
xyz[p].zero();
dxyzdxi_map[p].zero();
dxyzdeta_map[p].zero();
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxi2_map[p].zero();
d2xyzdxideta_map[p].zero();
d2xyzdeta2_map[p].zero();
#endif
// compute (x,y) at the quadrature points, derivatives once
for (unsigned int i=0; i<phi_map.size(); i++) // sum over the nodes
{
// Reference to the point, helps eliminate
// exessive temporaries in the inner loop
const Point& elem_point = elem->point(i);
xyz[p].add_scaled (elem_point, phi_map[i][p] );
dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p]);
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxi2_map[p].add_scaled (elem_point, d2phidxi2_map[i][p]);
d2xyzdxideta_map[p].add_scaled (elem_point, d2phidxideta_map[i][p]);
d2xyzdeta2_map[p].add_scaled (elem_point, d2phideta2_map[i][p]);
#endif
}
// compute the jacobian once
const Real dx_dxi = dxdxi_map(p), dx_deta = dxdeta_map(p),
dy_dxi = dydxi_map(p), dy_deta = dydeta_map(p),
dz_dxi = dzdxi_map(p), dz_deta = dzdeta_map(p);
#if LIBMESH_DIM == 2
// Compute the Jacobian. This assumes the 2D face
// lives in 2D space
//
// Symbolically, the matrix determinant is
//
// | dx/dxi dx/deta |
// jac = | dy/dxi dy/deta |
//
// jac = dx/dxi*dy/deta - dx/deta*dy/dxi
const Real jac = (dx_dxi*dy_deta - dx_deta*dy_dxi);
if (jac <= 0.)
{
std::cerr << 'ERROR: negative Jacobian: '
<< jac
<< ' in element '
<< elem->id()
<< std::endl;
libmesh_error();
}
JxW[p] = jac*qw[p];
// Compute the shape function derivatives wrt x,y at the
// quadrature points
const Real inv_jac = 1./jac;
dxidx_map[p] = dy_deta*inv_jac; //dxi/dx = (1/J)*dy/deta
dxidy_map[p] = -dx_deta*inv_jac; //dxi/dy = -(1/J)*dx/deta
detadx_map[p] = -dy_dxi* inv_jac; //deta/dx = -(1/J)*dy/dxi
detady_map[p] = dx_dxi* inv_jac; //deta/dy = (1/J)*dx/dxi
dxidz_map[p] = detadz_map[p] = 0.;
#else
// Compute the Jacobian. This assumes a 2D face in
// 3D space.
//
// The transformation matrix T from local to global
// coordinates is
//
// | dx/dxi dx/deta |
// T = | dy/dxi dy/deta |
// | dz/dxi dz/deta |
// note det(T' T) = det(T')det(T) = det(T)det(T)
// so det(T) = std::sqrt(det(T' T))
//
//----------------------------------------------
// Notes:
//
// dX = R dXi -> R'dX = R'R dXi
// (R^-1)dX = dXi [(R'R)^-1 R']dX = dXi
//
// so R^-1 = (R'R)^-1 R'
//
// and R^-1 R = (R'R)^-1 R'R = I.
//
const Real g11 = (dx_dxi*dx_dxi +
dy_dxi*dy_dxi +
dz_dxi*dz_dxi);
const Real g12 = (dx_dxi*dx_deta +
dy_dxi*dy_deta +
dz_dxi*dz_deta);
const Real g21 = g12;
const Real g22 = (dx_deta*dx_deta +
dy_deta*dy_deta +
dz_deta*dz_deta);
const Real det = (g11*g22 - g12*g21);
if (det <= 0.)
{
std::cerr << 'ERROR: negative Jacobian! '
<< ' in element '
<< elem->id()
<< std::endl;
libmesh_error();
}
const Real inv_det = 1./det;
const Real jac = std::sqrt(det);
JxW[p] = jac*qw[p];
const Real g11inv = g22*inv_det;
const Real g12inv = -g12*inv_det;
const Real g21inv = -g21*inv_det;
const Real g22inv = g11*inv_det;
dxidx_map[p] = g11inv*dx_dxi + g12inv*dx_deta;
dxidy_map[p] = g11inv*dy_dxi + g12inv*dy_deta;
dxidz_map[p] = g11inv*dz_dxi + g12inv*dz_deta;
detadx_map[p] = g21inv*dx_dxi + g22inv*dx_deta;
detady_map[p] = g21inv*dy_dxi + g22inv*dy_deta;
detadz_map[p] = g21inv*dz_dxi + g22inv*dz_deta;
#endif
// done computing the map
break;
}
//--------------------------------------------------------------------
// 3D
case 3:
{
//------------------------------------------------------------------
// Compute the (x,y,z) values at the quadrature points,
// the Jacobian at the quadrature point
// Clear the entities that will be summed
xyz[p].zero ();
dxyzdxi_map[p].zero ();
dxyzdeta_map[p].zero ();
dxyzdzeta_map[p].zero ();
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxi2_map[p].zero();
d2xyzdxideta_map[p].zero();
d2xyzdxidzeta_map[p].zero();
d2xyzdeta2_map[p].zero();
d2xyzdetadzeta_map[p].zero();
d2xyzdzeta2_map[p].zero();
#endif
// compute (x,y,z) at the quadrature points,
// dxdxi, dydxi, dzdxi,
// dxdeta, dydeta, dzdeta,
// dxdzeta, dydzeta, dzdzeta all once
for (unsigned int i=0; i<phi_map.size(); i++) // sum over the nodes
{
// Reference to the point, helps eliminate
// exessive temporaries in the inner loop
const Point& elem_point = elem->point(i);
xyz[p].add_scaled (elem_point, phi_map[i][p] );
dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p] );
dxyzdzeta_map[p].add_scaled (elem_point, dphidzeta_map[i][p]);
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxi2_map[p].add_scaled (elem_point,
d2phidxi2_map[i][p]);
d2xyzdxideta_map[p].add_scaled (elem_point,
d2phidxideta_map[i][p]);
d2xyzdxidzeta_map[p].add_scaled (elem_point,
d2phidxidzeta_map[i][p]);
d2xyzdeta2_map[p].add_scaled (elem_point,
d2phideta2_map[i][p]);
d2xyzdetadzeta_map[p].add_scaled (elem_point,
d2phidetadzeta_map[i][p]);
d2xyzdzeta2_map[p].add_scaled (elem_point,
d2phidzeta2_map[i][p]);
#endif
}
// compute the jacobian
const Real
dx_dxi = dxdxi_map(p), dy_dxi = dydxi_map(p), dz_dxi = dzdxi_map(p),
dx_deta = dxdeta_map(p), dy_deta = dydeta_map(p), dz_deta = dzdeta_map(p),
dx_dzeta = dxdzeta_map(p), dy_dzeta = dydzeta_map(p), dz_dzeta = dzdzeta_map(p);
// Symbolically, the matrix determinant is
//
// | dx/dxi dy/dxi dz/dxi |
// jac = | dx/deta dy/deta dz/deta |
// | dx/dzeta dy/dzeta dz/dzeta |
//
// jac = dx/dxi*(dy/deta*dz/dzeta - dz/deta*dy/dzeta) +
// dy/dxi*(dz/deta*dx/dzeta - dx/deta*dz/dzeta) +
// dz/dxi*(dx/deta*dy/dzeta - dy/deta*dx/dzeta)
const Real jac = (dx_dxi*(dy_deta*dz_dzeta - dz_deta*dy_dzeta) +
dy_dxi*(dz_deta*dx_dzeta - dx_deta*dz_dzeta) +
dz_dxi*(dx_deta*dy_dzeta - dy_deta*dx_dzeta));
if (jac <= 0.)
{
std::cerr << 'ERROR: negative Jacobian: '
<< jac
<< ' in element '
<< elem->id()
<< std::endl;
libmesh_error();
}
JxW[p] = jac*qw[p];
// Compute the shape function derivatives wrt x,y at the
// quadrature points
const Real inv_jac = 1./jac;
dxidx_map[p] = (dy_deta*dz_dzeta - dz_deta*dy_dzeta)*inv_jac;
dxidy_map[p] = (dz_deta*dx_dzeta - dx_deta*dz_dzeta)*inv_jac;
dxidz_map[p] = (dx_deta*dy_dzeta - dy_deta*dx_dzeta)*inv_jac;
detadx_map[p] = (dz_dxi*dy_dzeta - dy_dxi*dz_dzeta )*inv_jac;
detady_map[p] = (dx_dxi*dz_dzeta - dz_dxi*dx_dzeta )*inv_jac;
detadz_map[p] = (dy_dxi*dx_dzeta - dx_dxi*dy_dzeta )*inv_jac;
dzetadx_map[p] = (dy_dxi*dz_deta - dz_dxi*dy_deta )*inv_jac;
dzetady_map[p] = (dz_dxi*dx_deta - dx_dxi*dz_deta )*inv_jac;
dzetadz_map[p] = (dx_dxi*dy_deta - dy_dxi*dx_deta )*inv_jac;
// done computing the map
break;
}
default:
libmesh_error();
}
}
Definition at line 740 of file fe_base.h.
References dxyzdeta_map.
Referenced by compute_face_map(), and compute_single_point_map().
{ return dxyzdeta_map[p](0); }
Definition at line 719 of file fe_base.h.
References dxyzdxi_map.
Referenced by compute_edge_map(), compute_face_map(), and compute_single_point_map().
{ return dxyzdxi_map[p](0); }
Definition at line 761 of file fe_base.h.
References dxyzdzeta_map.
Referenced by compute_single_point_map().
{ return dxyzdzeta_map[p](0); }
Definition at line 747 of file fe_base.h.
References dxyzdeta_map.
Referenced by compute_face_map(), and compute_single_point_map().
{ return dxyzdeta_map[p](1); }
Definition at line 726 of file fe_base.h.
References dxyzdxi_map.
Referenced by compute_edge_map(), compute_face_map(), and compute_single_point_map().
{ return dxyzdxi_map[p](1); }
Definition at line 768 of file fe_base.h.
References dxyzdzeta_map.
Referenced by compute_single_point_map().
{ return dxyzdzeta_map[p](1); }
Definition at line 754 of file fe_base.h.
References dxyzdeta_map.
Referenced by compute_face_map(), and compute_single_point_map().
{ return dxyzdeta_map[p](2); }
Definition at line 733 of file fe_base.h.
References dxyzdxi_map.
Referenced by compute_edge_map(), compute_face_map(), and compute_single_point_map().
{ return dxyzdxi_map[p](2); }
Definition at line 775 of file fe_base.h.
References dxyzdzeta_map.
Referenced by compute_single_point_map().
{ return dxyzdzeta_map[p](2); }
Implemented in FE< Dim, T >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
Implemented in FE< Dim, T >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
Definition at line 539 of file fe_base.h.
References curvatures.
{ return curvatures;}
Definition at line 300 of file fe_base.h.
References calculate_d2phi, calculations_started, and d2phi.
Referenced by ExactErrorEstimator::find_squared_element_error().
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phi; }
Definition at line 308 of file fe_base.h.
References calculate_d2phi, calculations_started, and d2phidx2.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidx2; }
Definition at line 316 of file fe_base.h.
References calculate_d2phi, calculations_started, and d2phidxdy.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidxdy; }
Definition at line 324 of file fe_base.h.
References calculate_d2phi, calculations_started, and d2phidxdz.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidxdz; }
Definition at line 332 of file fe_base.h.
References calculate_d2phi, calculations_started, and d2phidy2.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidy2; }
Definition at line 340 of file fe_base.h.
References calculate_d2phi, calculations_started, and d2phidydz.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidydz; }
Definition at line 348 of file fe_base.h.
References calculate_d2phi, calculations_started, and d2phidz2.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidz2; }
Definition at line 384 of file fe_base.h.
References d2xyzdeta2_map.
{ return d2xyzdeta2_map; }
Definition at line 414 of file fe_base.h.
References d2xyzdetadzeta_map.
{ return d2xyzdetadzeta_map; }
Definition at line 378 of file fe_base.h.
References d2xyzdxi2_map.
{ return d2xyzdxi2_map; }
Definition at line 400 of file fe_base.h.
References d2xyzdxideta_map.
{ return d2xyzdxideta_map; }
Definition at line 408 of file fe_base.h.
References d2xyzdxidzeta_map.
{ return d2xyzdxidzeta_map; }
Definition at line 392 of file fe_base.h.
References d2xyzdzeta2_map.
{ return d2xyzdzeta2_map; }
Definition at line 444 of file fe_base.h.
References detadx_map.
{ return detadx_map; }
Definition at line 451 of file fe_base.h.
References detady_map.
{ return detady_map; }
Definition at line 458 of file fe_base.h.
References detadz_map.
{ return detadz_map; }
In case of the general finite element class FE this field is initialized to all zero, so that the variational formulation for an infinite element returns correct element matrices for a mesh using both finite and infinite elements.
Definition at line 494 of file fe_base.h.
References dphase.
{ return dphase; }
Definition at line 242 of file fe_base.h.
References calculate_dphi, calculations_started, and dphi.
Referenced by ExactErrorEstimator::find_squared_element_error().
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphi; }
Definition at line 282 of file fe_base.h.
References calculate_dphi, calculations_started, and dphideta.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphideta; }
Definition at line 250 of file fe_base.h.
References calculate_dphi, calculations_started, and dphidx.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidx; }
Definition at line 274 of file fe_base.h.
References calculate_dphi, calculations_started, and dphidxi.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidxi; }
Definition at line 258 of file fe_base.h.
References calculate_dphi, calculations_started, and dphidy.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidy; }
Definition at line 266 of file fe_base.h.
References calculate_dphi, calculations_started, and dphidz.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidz; }
Definition at line 290 of file fe_base.h.
References calculate_dphi, calculations_started, and dphidzeta.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidzeta; }
Definition at line 423 of file fe_base.h.
References dxidx_map.
{ return dxidx_map; }
Definition at line 430 of file fe_base.h.
References dxidy_map.
{ return dxidy_map; }
Definition at line 437 of file fe_base.h.
References dxidz_map.
{ return dxidz_map; }
Definition at line 365 of file fe_base.h.
References dxyzdeta_map.
{ return dxyzdeta_map; }
Definition at line 358 of file fe_base.h.
References dxyzdxi_map.
{ return dxyzdxi_map; }
Definition at line 372 of file fe_base.h.
References dxyzdzeta_map.
{ return dxyzdzeta_map; }
Definition at line 465 of file fe_base.h.
References dzetadx_map.
{ return dzetadx_map; }
Definition at line 472 of file fe_base.h.
References dzetady_map.
{ return dzetady_map; }
Definition at line 479 of file fe_base.h.
References dzetadz_map.
{ return dzetadz_map; }
Definition at line 598 of file fe_base.h.
References FEType::family, and fe_type.
{ return fe_type.family; }
Definition at line 577 of file fe_base.h.
References fe_type.
{ return fe_type; }
Definition at line 45 of file reference_counter.C.
References ReferenceCounter::_counts, and Quality::name().
Referenced by ReferenceCounter::print_info().
{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
std::ostringstream out;
out << '
<< ' ----------------------------------------------------------------------------
<< '| Reference count information |
<< ' ---------------------------------------------------------------------------- ;
for (Counts::iterator it = _counts.begin();
it != _counts.end(); ++it)
{
const std::string name(it->first);
const unsigned int creations = it->second.first;
const unsigned int destructions = it->second.second;
out << '| ' << name << ' reference count information:
<< '| Creations: ' << creations << '
<< '| Destructions: ' << destructions << ';
}
out << ' ---------------------------------------------------------------------------- ;
return out.str();
#else
return '';
#endif
}
Definition at line 235 of file fe_base.h.
References JxW.
Referenced by ExactErrorEstimator::find_squared_element_error().
{ return JxW; }
Definition at line 533 of file fe_base.h.
References normals.
{ return normals; }
Definition at line 582 of file fe_base.h.
References _p_level, fe_type, and FEType::order.
{ return static_cast<Order>(fe_type.order + _p_level); }
Definition at line 572 of file fe_base.h.
References _p_level.
Referenced by REINIT_ERROR().
{ return _p_level; }
Definition at line 227 of file fe_base.h.
References calculate_phi, calculations_started, and phi.
Referenced by ExactErrorEstimator::find_squared_element_error().
{ libmesh_assert(!calculations_started || calculate_phi);
calculate_phi = true; return phi; }
Definition at line 518 of file fe_base.h.
References dweight.
{ return dweight; }
In case of the general finite element class FE this field is initialized to all ones, so that the variational formulation for an infinite element returns correct element matrices for a mesh using both finite and infinite elements.
Definition at line 510 of file fe_base.h.
References weight.
{ return weight; }
Definition at line 527 of file fe_base.h.
References tangents.
{ return tangents; }
Definition at line 566 of file fe_base.h.
References elem_type.
Referenced by FE< Dim, T >::edge_reinit(), InfFE< Dim, T_radial, T_map >::reinit(), and REINIT_ERROR().
{ return elem_type; }
Definition at line 220 of file fe_base.h.
References xyz.
Referenced by ExactErrorEstimator::find_squared_element_error().
{ return xyz; }
Definition at line 149 of file reference_counter.h.
References ReferenceCounter::_counts, Quality::name(), and Threads::spin_mtx.
Referenced by ReferenceCountedObject< Value >::ReferenceCountedObject().
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
std::pair<unsigned int, unsigned int>& p = _counts[name];
p.first++;
}
Definition at line 167 of file reference_counter.h.
References ReferenceCounter::_counts, Quality::name(), and Threads::spin_mtx.
Referenced by ReferenceCountedObject< Value >::~ReferenceCountedObject().
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
std::pair<unsigned int, unsigned int>& p = _counts[name];
p.second++;
}
Implemented in FE< Dim, T >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
Referenced by InfFE< Dim, T_radial, T_map >::init_face_shape_functions().
Implemented in FE< Dim, T >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
Definition at line 76 of file reference_counter.h.
References ReferenceCounter::_n_objects.
Referenced by System::read_serialized_blocked_dof_objects(), and System::write_serialized_blocked_dof_objects().
{ return _n_objects; }
Implemented in FE< Dim, T >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
Implemented in FE< Dim, T >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
References libMeshEnums::EDGE2, libMeshEnums::EDGE3, libMeshEnums::EDGE4, libMeshEnums::HEX20, libMeshEnums::HEX27, libMeshEnums::HEX8, libMeshEnums::INFHEX8, libMeshEnums::INFPRISM6, libMeshEnums::PRISM15, libMeshEnums::PRISM18, libMeshEnums::PRISM6, libMeshEnums::PYRAMID5, libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, libMeshEnums::TET10, libMeshEnums::TET4, libMeshEnums::TRI3, and libMeshEnums::TRI6.
{
libmesh_assert (eps >= 0.);
const Real xi = p(0);
const Real eta = p(1);
const Real zeta = p(2);
switch (t)
{
case EDGE2:
case EDGE3:
case EDGE4:
{
// The reference 1D element is [-1,1].
if ((xi >= -1.-eps) &&
(xi <= 1.+eps))
return true;
return false;
}
case TRI3:
case TRI6:
{
// The reference triangle is isocoles
// and is bound by xi=0, eta=0, and xi+eta=1.
if ((xi >= 0.-eps) &&
(eta >= 0.-eps) &&
((xi + eta) <= 1.+eps))
return true;
return false;
}
case QUAD4:
case QUAD8:
case QUAD9:
{
// The reference quadrilateral element is [-1,1]^2.
if ((xi >= -1.-eps) &&
(xi <= 1.+eps) &&
(eta >= -1.-eps) &&
(eta <= 1.+eps))
return true;
return false;
}
case TET4:
case TET10:
{
// The reference tetrahedral is isocoles
// and is bound by xi=0, eta=0, zeta=0,
// and xi+eta+zeta=1.
if ((xi >= 0.-eps) &&
(eta >= 0.-eps) &&
(zeta >= 0.-eps) &&
((xi + eta + zeta) <= 1.+eps))
return true;
return false;
}
case HEX8:
case HEX20:
case HEX27:
{
/*
if ((xi >= -1.) &&
(xi <= 1.) &&
(eta >= -1.) &&
(eta <= 1.) &&
(zeta >= -1.) &&
(zeta <= 1.))
return true;
*/
// The reference hexahedral element is [-1,1]^3.
if ((xi >= -1.-eps) &&
(xi <= 1.+eps) &&
(eta >= -1.-eps) &&
(eta <= 1.+eps) &&
(zeta >= -1.-eps) &&
(zeta <= 1.+eps))
{
// std::cout << 'Strange Point:;
// p.print();
return true;
}
return false;
}
case PRISM6:
case PRISM15:
case PRISM18:
{
// Figure this one out...
// inside the reference triange with zeta in [-1,1]
if ((xi >= 0.-eps) &&
(eta >= 0.-eps) &&
(zeta >= -1.-eps) &&
(zeta <= 1.+eps) &&
((xi + eta) <= 1.+eps))
return true;
return false;
}
case PYRAMID5:
{
std::cerr << 'BEN: Implement this you lazy bastard!'
<< std::endl;
libmesh_error();
return false;
}
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
case INFHEX8:
{
// The reference infhex8 is a [-1,1]^3.
if ((xi >= -1.-eps) &&
(xi <= 1.+eps) &&
(eta >= -1.-eps) &&
(eta <= 1.+eps) &&
(zeta >= -1.-eps) &&
(zeta <= 1.+eps))
{
return true;
}
return false;
}
case INFPRISM6:
{
// inside the reference triange with zeta in [-1,1]
if ((xi >= 0.-eps) &&
(eta >= 0.-eps) &&
(zeta >= -1.-eps) &&
(zeta <= 1.+eps) &&
((xi + eta) <= 1.+eps))
{
return true;
}
return false;
}
#endif
default:
std::cerr << 'ERROR: Unknown element type ' << t << std::endl;
libmesh_error();
}
// If we get here then the point is _not_ in the
// reference element. Better return false.
return false;
}
Definition at line 1069 of file fe_base.C.
References d2phi, and dphi.
{
for (unsigned int i=0; i<dphi.size(); ++i)
for (unsigned int j=0; j<dphi[i].size(); ++j)
os << ' d2phi[' << i << '][' << j << ']=' << d2phi[i][j];
}
Definition at line 1057 of file fe_base.C.
References dphi.
Referenced by print_info().
{
for (unsigned int i=0; i<dphi.size(); ++i)
for (unsigned int j=0; j<dphi[i].size(); ++j)
os << ' dphi[' << i << '][' << j << ']=' << dphi[i][j];
}
Definition at line 1090 of file fe_base.C.
References print_dphi(), print_JxW(), print_phi(), and print_xyz().
Referenced by operator<<().
{
os << 'Shape functions at the Gauss pts.' << std::endl;
this->print_phi(os);
os << 'Shape function gradients at the Gauss pts.' << std::endl;
this->print_dphi(os);
os << 'XYZ locations of the Gauss pts.' << std::endl;
this->print_xyz(os);
os << 'Values of JxW at the Gauss pts.' << std::endl;
this->print_JxW(os);
}
Definition at line 83 of file reference_counter.C.
References ReferenceCounter::get_info().
{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
std::cout << ReferenceCounter::get_info();
#endif
}
Definition at line 1038 of file fe_base.C.
References JxW.
Referenced by print_info().
{
for (unsigned int i=0; i<JxW.size(); ++i)
os << JxW[i] << std::endl;
}
Definition at line 1047 of file fe_base.C.
References phi.
Referenced by print_info().
{
for (unsigned int i=0; i<phi.size(); ++i)
for (unsigned int j=0; j<phi[i].size(); ++j)
os << ' phi[' << i << '][' << j << ']=' << phi[i][j] << std::endl;
}
Definition at line 1081 of file fe_base.C.
References xyz.
Referenced by print_info().
{
for (unsigned int i=0; i<xyz.size(); ++i)
os << xyz[i];
}
Note that 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 FE< Dim, T >, FEXYZ< Dim >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
Referenced by ExactErrorEstimator::find_squared_element_error().
Implemented in FE< Dim, T >, FEXYZ< Dim >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, FE< Dim, LAGRANGE >, FEXYZ< Dim >, and FEXYZ< Dim >.
Definition at line 372 of file fe_map.C.
References d2xyzdeta2_map, d2xyzdetadzeta_map, d2xyzdxi2_map, d2xyzdxideta_map, d2xyzdxidzeta_map, d2xyzdzeta2_map, detadx_map, detady_map, detadz_map, dim, dxidx_map, dxidy_map, dxidz_map, dxyzdeta_map, dxyzdxi_map, dxyzdzeta_map, dzetadx_map, dzetady_map, dzetadz_map, JxW, and xyz.
Referenced by compute_affine_map(), and compute_map().
{
// Resize the vectors to hold data at the quadrature points
xyz.resize(n_qp);
dxyzdxi_map.resize(n_qp);
dxidx_map.resize(n_qp);
dxidy_map.resize(n_qp); // 1D element may live in 2D ...
dxidz_map.resize(n_qp); // ... or 3D
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxi2_map.resize(n_qp);
#endif
if (this->dim > 1)
{
dxyzdeta_map.resize(n_qp);
detadx_map.resize(n_qp);
detady_map.resize(n_qp);
detadz_map.resize(n_qp);
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxideta_map.resize(n_qp);
d2xyzdeta2_map.resize(n_qp);
#endif
if (this->dim > 2)
{
dxyzdzeta_map.resize (n_qp);
dzetadx_map.resize (n_qp);
dzetady_map.resize (n_qp);
dzetadz_map.resize (n_qp);
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
d2xyzdxidzeta_map.resize(n_qp);
d2xyzdetadzeta_map.resize(n_qp);
d2xyzdzeta2_map.resize(n_qp);
#endif
}
}
JxW.resize(n_qp);
}
Implemented in FE< Dim, T >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
Referenced by InfFE< Dim, T_radial, T_map >::reinit().
Reimplemented in FE< Dim, T >, InfFE< Dim, T_radial, T_map >, FE< Dim, HIERARCHIC >, FE< Dim, SCALAR >, FE< Dim, HERMITE >, FE< Dim, CLOUGH >, FE< Dim, MONOMIAL >, FE< Dim, XYZ >, and FE< Dim, LAGRANGE >.
Definition at line 1248 of file fe_base.h.
Definition at line 1108 of file fe_base.C.
{
fe.print_info(os);
return os;
}
Definition at line 110 of file reference_counter.h.
Referenced by ReferenceCounter::get_info(), ReferenceCounter::increment_constructor_count(), and ReferenceCounter::increment_destructor_count().
Definition at line 123 of file reference_counter.h.
Definition at line 118 of file reference_counter.h.
Referenced by ReferenceCounter::n_objects(), ReferenceCounter::ReferenceCounter(), and ReferenceCounter::~ReferenceCounter().
Definition at line 1215 of file fe_base.h.
Referenced by get_order(), get_p_level(), and REINIT_ERROR().
Definition at line 936 of file fe_base.h.
Referenced by compute_map(), compute_shape_functions(), get_d2phi(), get_d2phidx2(), get_d2phidxdy(), get_d2phidxdz(), get_d2phidy2(), get_d2phidydz(), and get_d2phidz2().
Definition at line 931 of file fe_base.h.
Referenced by compute_shape_functions(), get_dphi(), get_dphideta(), get_dphidx(), get_dphidxi(), get_dphidy(), get_dphidz(), and get_dphidzeta().
Definition at line 926 of file fe_base.h.
Referenced by compute_shape_functions(), and get_phi().
Definition at line 921 of file fe_base.h.
Referenced by compute_shape_functions(), get_d2phi(), get_d2phidx2(), get_d2phidxdy(), get_d2phidxdz(), get_d2phidy2(), get_d2phidydz(), get_d2phidz2(), get_dphi(), get_dphideta(), get_dphidx(), get_dphidxi(), get_dphidy(), get_dphidz(), get_dphidzeta(), and get_phi().
Definition at line 1192 of file fe_base.h.
Referenced by compute_edge_map(), compute_face_map(), and get_curvatures().
Definition at line 979 of file fe_base.h.
Referenced by compute_shape_functions(), get_d2phi(), and print_d2phi().
Definition at line 999 of file fe_base.h.
Referenced by compute_shape_functions().
Definition at line 1086 of file fe_base.h.
Referenced by compute_single_point_map().
Definition at line 1004 of file fe_base.h.
Referenced by compute_shape_functions().
Definition at line 1091 of file fe_base.h.
Referenced by compute_single_point_map().
Definition at line 1014 of file fe_base.h.
Referenced by compute_shape_functions(), and get_d2phidx2().
Definition at line 1019 of file fe_base.h.
Referenced by compute_shape_functions(), and get_d2phidxdy().
Definition at line 1024 of file fe_base.h.
Referenced by compute_shape_functions(), and get_d2phidxdz().
Definition at line 984 of file fe_base.h.
Referenced by compute_shape_functions().
Definition at line 1071 of file fe_base.h.
Referenced by compute_single_point_map().
Definition at line 989 of file fe_base.h.
Referenced by compute_shape_functions().
Definition at line 1076 of file fe_base.h.
Referenced by compute_single_point_map().
Definition at line 994 of file fe_base.h.
Referenced by compute_shape_functions().
Definition at line 1081 of file fe_base.h.
Referenced by compute_single_point_map().
Definition at line 1029 of file fe_base.h.
Referenced by compute_shape_functions(), and get_d2phidy2().
Definition at line 1034 of file fe_base.h.
Referenced by compute_shape_functions(), and get_d2phidydz().
Definition at line 1039 of file fe_base.h.
Referenced by compute_shape_functions(), and get_d2phidz2().
Definition at line 1009 of file fe_base.h.
Referenced by compute_shape_functions().
Definition at line 1096 of file fe_base.h.
Referenced by compute_single_point_map().
Definition at line 1140 of file fe_base.h.
Referenced by compute_face_map(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), and FE< Dim, T >::init_face_shape_functions().
Definition at line 1126 of file fe_base.h.
Referenced by compute_edge_map(), compute_face_map(), FE< Dim, T >::init_edge_shape_functions(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), and FE< Dim, T >::init_face_shape_functions().
Definition at line 1133 of file fe_base.h.
Referenced by compute_face_map(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), and FE< Dim, T >::init_face_shape_functions().
Definition at line 827 of file fe_base.h.
Referenced by compute_affine_map(), compute_edge_map(), compute_face_map(), compute_single_point_map(), get_d2xyzdeta2(), and resize_map_vectors().
Definition at line 841 of file fe_base.h.
Referenced by compute_affine_map(), compute_single_point_map(), get_d2xyzdetadzeta(), and resize_map_vectors().
Definition at line 815 of file fe_base.h.
Referenced by compute_affine_map(), compute_edge_map(), compute_face_map(), compute_single_point_map(), get_d2xyzdxi2(), and resize_map_vectors().
Definition at line 821 of file fe_base.h.
Referenced by compute_affine_map(), compute_edge_map(), compute_face_map(), compute_single_point_map(), get_d2xyzdxideta(), and resize_map_vectors().
Definition at line 835 of file fe_base.h.
Referenced by compute_affine_map(), compute_single_point_map(), get_d2xyzdxidzeta(), and resize_map_vectors().
Definition at line 847 of file fe_base.h.
Referenced by compute_affine_map(), compute_single_point_map(), get_d2xyzdzeta2(), and resize_map_vectors().
Definition at line 876 of file fe_base.h.
Referenced by compute_affine_map(), compute_shape_functions(), compute_single_point_map(), get_detadx(), and resize_map_vectors().
Definition at line 882 of file fe_base.h.
Referenced by compute_affine_map(), compute_shape_functions(), compute_single_point_map(), get_detady(), and resize_map_vectors().
Definition at line 888 of file fe_base.h.
Referenced by compute_affine_map(), compute_shape_functions(), compute_single_point_map(), get_detadz(), and resize_map_vectors().
Definition at line 784 of file fe_base.h.
Referenced by JumpErrorEstimator::coarse_n_flux_faces_increment(), coarsened_dof_values(), compute_affine_map(), compute_edge_map(), compute_face_map(), compute_shape_functions(), compute_single_point_map(), JumpErrorEstimator::estimate_error(), and resize_map_vectors().
Definition at line 1156 of file fe_base.h.
Definition at line 941 of file fe_base.h.
Referenced by compute_periodic_constraints(), compute_proj_constraints(), compute_shape_functions(), get_dphi(), print_d2phi(), and print_dphi().
Definition at line 951 of file fe_base.h.
Referenced by compute_shape_functions(), and get_dphideta().
Definition at line 1059 of file fe_base.h.
Referenced by compute_single_point_map().
Definition at line 961 of file fe_base.h.
Referenced by compute_shape_functions(), and get_dphidx().
Definition at line 946 of file fe_base.h.
Referenced by compute_shape_functions(), and get_dphidxi().
Definition at line 1054 of file fe_base.h.
Referenced by compute_single_point_map(), and InfFE< Dim, T_radial, T_map >::init_face_shape_functions().
Definition at line 966 of file fe_base.h.
Referenced by compute_shape_functions(), and get_dphidy().
Definition at line 971 of file fe_base.h.
Referenced by compute_shape_functions(), and get_dphidz().
Definition at line 956 of file fe_base.h.
Referenced by compute_shape_functions(), and get_dphidzeta().
Definition at line 1064 of file fe_base.h.
Referenced by compute_single_point_map().
Definition at line 1119 of file fe_base.h.
Referenced by compute_face_map(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), and FE< Dim, T >::init_face_shape_functions().
Definition at line 1113 of file fe_base.h.
Referenced by compute_edge_map(), compute_face_map(), FE< Dim, T >::init_edge_shape_functions(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), and FE< Dim, T >::init_face_shape_functions().
Definition at line 1163 of file fe_base.h.
Referenced by get_Sobolev_dweight().
Definition at line 855 of file fe_base.h.
Referenced by compute_affine_map(), compute_shape_functions(), compute_single_point_map(), get_dxidx(), and resize_map_vectors().
Definition at line 861 of file fe_base.h.
Referenced by compute_affine_map(), compute_shape_functions(), compute_single_point_map(), get_dxidy(), and resize_map_vectors().
Definition at line 867 of file fe_base.h.
Referenced by compute_affine_map(), compute_shape_functions(), compute_single_point_map(), get_dxidz(), and resize_map_vectors().
Definition at line 803 of file fe_base.h.
Referenced by compute_affine_map(), compute_edge_map(), compute_face_map(), compute_single_point_map(), dxdeta_map(), dydeta_map(), dzdeta_map(), get_dxyzdeta(), and resize_map_vectors().
Definition at line 797 of file fe_base.h.
Referenced by compute_affine_map(), compute_edge_map(), compute_face_map(), compute_single_point_map(), dxdxi_map(), dydxi_map(), dzdxi_map(), get_dxyzdxi(), and resize_map_vectors().
Definition at line 809 of file fe_base.h.
Referenced by compute_affine_map(), compute_single_point_map(), dxdzeta_map(), dydzeta_map(), dzdzeta_map(), get_dxyzdzeta(), and resize_map_vectors().
Definition at line 898 of file fe_base.h.
Referenced by compute_affine_map(), compute_shape_functions(), compute_single_point_map(), get_dzetadx(), and resize_map_vectors().
Definition at line 904 of file fe_base.h.
Referenced by compute_affine_map(), compute_shape_functions(), compute_single_point_map(), get_dzetady(), and resize_map_vectors().
Definition at line 910 of file fe_base.h.
Referenced by compute_affine_map(), compute_shape_functions(), compute_single_point_map(), get_dzetadz(), and resize_map_vectors().
Definition at line 1209 of file fe_base.h.
Referenced by coarsened_dof_values(), FE< Dim, T >::edge_reinit(), get_type(), InfFE< Dim, T_radial, T_map >::reinit(), and REINIT_ERROR().
Definition at line 1203 of file fe_base.h.
Referenced by coarsened_dof_values(), JumpErrorEstimator::estimate_error(), FE< Dim, T >::FE(), get_family(), get_fe_type(), get_order(), InfFE< Dim, T_radial, T_map >::InfFE(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), and InfFE< Dim, T_radial, T_map >::reinit().
Definition at line 1197 of file fe_base.h.
Referenced by coarsened_dof_values(), compute_affine_map(), compute_edge_map(), compute_face_map(), compute_periodic_constraints(), compute_proj_constraints(), compute_single_point_map(), FE< Dim, T >::edge_reinit(), get_JxW(), print_JxW(), InfFE< Dim, T_radial, T_map >::reinit(), REINIT_ERROR(), and resize_map_vectors().
Definition at line 1185 of file fe_base.h.
Referenced by compute_edge_map(), compute_face_map(), and get_normals().
Definition at line 915 of file fe_base.h.
Referenced by compute_periodic_constraints(), compute_proj_constraints(), get_phi(), and print_phi().
Definition at line 1049 of file fe_base.h.
Referenced by compute_affine_map(), compute_single_point_map(), and InfFE< Dim, T_radial, T_map >::init_face_shape_functions().
Definition at line 1107 of file fe_base.h.
Referenced by compute_edge_map(), compute_face_map(), FE< Dim, T >::init_edge_shape_functions(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), and FE< Dim, T >::init_face_shape_functions().
Definition at line 1220 of file fe_base.h.
Referenced by coarsened_dof_values(), FE< Dim, T >::edge_reinit(), JumpErrorEstimator::estimate_error(), InfFE< Dim, T_radial, T_map >::init_face_shape_functions(), InfFE< Dim, T_radial, T_map >::reinit(), and REINIT_ERROR().
Definition at line 1226 of file fe_base.h.
Definition at line 1180 of file fe_base.h.
Referenced by compute_edge_map(), compute_face_map(), and get_tangents().
Definition at line 1170 of file fe_base.h.
Referenced by get_Sobolev_weight().
Definition at line 789 of file fe_base.h.
Referenced by compute_affine_map(), compute_edge_map(), compute_face_map(), compute_single_point_map(), FE< Dim, T >::edge_reinit(), get_xyz(), print_xyz(), InfFE< Dim, T_radial, T_map >::reinit(), REINIT_ERROR(), and resize_map_vectors().
Generated automatically by Doxygen for libMesh from the source code.