#include <fe.h>
Inherits FE< Dim, HIERARCHIC >.
FEHierarchic (const FEType &fet)
virtual unsigned int n_shape_functions () const
virtual FEContinuity get_continuity () const
virtual bool is_hierarchic () const
virtual void reinit (const Elem *elem, const std::vector< Point > *const pts=NULL)
virtual void reinit (const Elem *elem, const unsigned int side, const Real tolerance=TOLERANCE)
virtual void edge_reinit (const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE)
virtual void attach_quadrature_rule (QBase *q)
virtual unsigned int n_quadrature_points () const
void init_base_shape_functions (const std::vector< Point > &qp, const Elem *e)
virtual bool shapes_need_reinit () const
virtual void init_shape_functions (const std::vector< Point > &qp, const Elem *e)
void init_face_shape_functions (const std::vector< Point > &qp, const Elem *side)
void init_edge_shape_functions (const std::vector< Point > &qp, const Elem *edge)
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
ElemType get_type () const
unsigned int get_p_level () const
FEType get_fe_type () const
Order get_order () const
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 Real shape (const ElemType t, const Order o, const unsigned int i, const Point &p)
static Real shape (const Elem *elem, const Order o, const unsigned int i, const Point &p)
static Real shape_deriv (const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static Real shape_deriv (const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static Real shape_second_deriv (const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static Real shape_second_deriv (const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static void nodal_soln (const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
static unsigned int n_shape_functions (const ElemType t, const Order o)
static unsigned int n_dofs (const ElemType t, const Order o)
static unsigned int n_dofs_at_node (const ElemType t, const Order o, const unsigned int n)
static unsigned int n_dofs_per_elem (const ElemType t, const Order o)
static void dofs_on_side (const Elem *const elem, const Order o, unsigned int s, std::vector< unsigned int > &di)
static void dofs_on_edge (const Elem *const elem, const Order o, unsigned int e, std::vector< unsigned int > &di)
static Point inverse_map (const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
static void inverse_map (const Elem *elem, const std::vector< Point > &physical_points, std::vector< Point > &reference_points, const Real tolerance=TOLERANCE, const bool secure=true)
static void compute_constraints (DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
static Point map (const Elem *elem, const Point &reference_point)
static Point map_xi (const Elem *elem, const Point &reference_point)
static Point map_eta (const Elem *elem, const Point &reference_point)
static Point map_zeta (const Elem *elem, const Point &reference_point)
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 void print_info ()
static std::string get_info ()
static unsigned int n_objects ()
typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
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)
std::vector< Point > cached_nodes
unsigned int last_side
unsigned int last_edge
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
class InfFE
std::ostream & operator<< (std::ostream &os, const FEBase &fe)
Author:
Date:
Version:
Revision:
Definition at line 514 of file fe.h.
Definition at line 105 of file reference_counter.h.
Definition at line 784 of file fe.h.
:
FE<Dim,HIERARCHIC> (fet)
{
}
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(), FEBase::coarsened_dof_values(), FEBase::compute_periodic_constraints(), FEBase::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(), FEBase::build(), libMeshEnums::C_ONE, Elem::child(), DenseMatrix< T >::cholesky_solve(), FEType::default_quadrature_rule(), Elem::dim(), FEBase::dim, libMeshEnums::DISCONTINUOUS, DofMap::dof_indices(), FEInterface::dofs_on_edge(), FEInterface::dofs_on_side(), FEBase::elem_type, FEBase::fe_type, FEInterface::inverse_map(), Elem::is_child_on_edge(), Elem::is_child_on_side(), Elem::is_vertex(), FEBase::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(), FEBase::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 FEBase::compute_single_point_map(), FEBase::d2xyzdeta2_map, FEBase::d2xyzdetadzeta_map, FEBase::d2xyzdxi2_map, FEBase::d2xyzdxideta_map, FEBase::d2xyzdxidzeta_map, FEBase::d2xyzdzeta2_map, FEBase::detadx_map, FEBase::detady_map, FEBase::detadz_map, FEBase::dim, FEBase::dxidx_map, FEBase::dxidy_map, FEBase::dxidz_map, FEBase::dxyzdeta_map, FEBase::dxyzdxi_map, FEBase::dxyzdzeta_map, FEBase::dzetadx_map, FEBase::dzetady_map, FEBase::dzetadz_map, FEBase::JxW, FEBase::phi_map, Elem::point(), FEBase::resize_map_vectors(), and FEBase::xyz.
Referenced by FEBase::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 FEBase::compute_face_map(), FEBase::curvatures, FEBase::d2psidxi2_map, FEBase::d2xyzdeta2_map, FEBase::d2xyzdxi2_map, FEBase::d2xyzdxideta_map, FEBase::dim, FEBase::dpsidxi_map, FEBase::dxdxi_map(), FEBase::dxyzdeta_map, FEBase::dxyzdxi_map, FEBase::dydxi_map(), FEBase::dzdxi_map(), FEBase::JxW, FEBase::normals, Elem::point(), FEBase::psi_map, FEBase::tangents, and FEBase::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(), FEBase::curvatures, FEBase::d2psideta2_map, FEBase::d2psidxi2_map, FEBase::d2psidxideta_map, FEBase::d2xyzdeta2_map, FEBase::d2xyzdxi2_map, FEBase::d2xyzdxideta_map, FEBase::dim, FEBase::dpsideta_map, FEBase::dpsidxi_map, FEBase::dxdeta_map(), FEBase::dxdxi_map(), FEBase::dxyzdeta_map, FEBase::dxyzdxi_map, FEBase::dydeta_map(), FEBase::dydxi_map(), FEBase::dzdeta_map(), FEBase::dzdxi_map(), InfFE< Dim, T_radial, T_map >::inverse_map(), FEBase::JxW, Elem::node(), FEBase::normals, Elem::parent(), Elem::point(), FEBase::psi_map, FEBase::tangents, TypeVector< T >::unit(), and FEBase::xyz.
Referenced by FEBase::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 FEBase::calculate_d2phi, FEBase::compute_affine_map(), FEBase::compute_single_point_map(), Elem::has_affine_map(), and FEBase::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, FEBase::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(), FEBase::dphi, DofObject::invalid_id, libMesh::invalid_uint, FEInterface::inverse_map(), DofMap::is_constrained_dof(), FEBase::JxW, Elem::level(), std::min(), Elem::min_p_level_by_neighbor(), Elem::n_sides(), PeriodicBoundaries::neighbor(), Elem::neighbor(), Elem::p_level(), PeriodicBoundary::pairedboundary, FEBase::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(), FEBase::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(), FEBase::dphi, DofObject::invalid_id, FEInterface::inverse_map(), FEBase::JxW, Elem::level(), std::min(), Elem::min_p_level_by_neighbor(), Elem::n_nodes(), Elem::n_sides(), Elem::neighbor(), Elem::p_level(), FEBase::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 FEBase::calculate_d2phi, FEBase::calculate_dphi, FEBase::calculate_phi, FEBase::calculations_started, FEBase::d2phi, FEBase::d2phideta2, FEBase::d2phidetadzeta, FEBase::d2phidx2, FEBase::d2phidxdy, FEBase::d2phidxdz, FEBase::d2phidxi2, FEBase::d2phidxideta, FEBase::d2phidxidzeta, FEBase::d2phidy2, FEBase::d2phidydz, FEBase::d2phidz2, FEBase::d2phidzeta2, FEBase::detadx_map, FEBase::detady_map, FEBase::detadz_map, FEBase::dim, FEBase::dphi, FEBase::dphideta, FEBase::dphidx, FEBase::dphidxi, FEBase::dphidy, FEBase::dphidz, FEBase::dphidzeta, FEBase::dxidx_map, FEBase::dxidy_map, FEBase::dxidz_map, FEBase::dzetadx_map, FEBase::dzetady_map, and FEBase::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 FEBase::d2phideta2_map, FEBase::d2phidetadzeta_map, FEBase::d2phidxi2_map, FEBase::d2phidxideta_map, FEBase::d2phidxidzeta_map, FEBase::d2phidzeta2_map, FEBase::d2xyzdeta2_map, FEBase::d2xyzdetadzeta_map, FEBase::d2xyzdxi2_map, FEBase::d2xyzdxideta_map, FEBase::d2xyzdxidzeta_map, FEBase::d2xyzdzeta2_map, FEBase::detadx_map, FEBase::detady_map, FEBase::detadz_map, FEBase::dim, FEBase::dphideta_map, FEBase::dphidxi_map, FEBase::dphidzeta_map, FEBase::dxdeta_map(), FEBase::dxdxi_map(), FEBase::dxdzeta_map(), FEBase::dxidx_map, FEBase::dxidy_map, FEBase::dxidz_map, FEBase::dxyzdeta_map, FEBase::dxyzdxi_map, FEBase::dxyzdzeta_map, FEBase::dydeta_map(), FEBase::dydxi_map(), FEBase::dydzeta_map(), FEBase::dzdeta_map(), FEBase::dzdxi_map(), FEBase::dzdzeta_map(), FEBase::dzetadx_map, FEBase::dzetady_map, FEBase::dzetadz_map, DofObject::id(), FEBase::JxW, FEBase::phi_map, Elem::point(), and FEBase::xyz.
Referenced by FEBase::compute_affine_map(), and FEBase::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();
}
}
On a p-refined element, o should be the base order of the element.
On a p-refined element, o should be the base order of the element.
Definition at line 740 of file fe_base.h.
References FEBase::dxyzdeta_map.
Referenced by FEBase::compute_face_map(), and FEBase::compute_single_point_map().
{ return dxyzdeta_map[p](0); }
Definition at line 719 of file fe_base.h.
References FEBase::dxyzdxi_map.
Referenced by FEBase::compute_edge_map(), FEBase::compute_face_map(), and FEBase::compute_single_point_map().
{ return dxyzdxi_map[p](0); }
Definition at line 761 of file fe_base.h.
References FEBase::dxyzdzeta_map.
Referenced by FEBase::compute_single_point_map().
{ return dxyzdzeta_map[p](0); }
Definition at line 747 of file fe_base.h.
References FEBase::dxyzdeta_map.
Referenced by FEBase::compute_face_map(), and FEBase::compute_single_point_map().
{ return dxyzdeta_map[p](1); }
Definition at line 726 of file fe_base.h.
References FEBase::dxyzdxi_map.
Referenced by FEBase::compute_edge_map(), FEBase::compute_face_map(), and FEBase::compute_single_point_map().
{ return dxyzdxi_map[p](1); }
Definition at line 768 of file fe_base.h.
References FEBase::dxyzdzeta_map.
Referenced by FEBase::compute_single_point_map().
{ return dxyzdzeta_map[p](1); }
Definition at line 754 of file fe_base.h.
References FEBase::dxyzdeta_map.
Referenced by FEBase::compute_face_map(), and FEBase::compute_single_point_map().
{ return dxyzdeta_map[p](2); }
Definition at line 733 of file fe_base.h.
References FEBase::dxyzdxi_map.
Referenced by FEBase::compute_edge_map(), FEBase::compute_face_map(), and FEBase::compute_single_point_map().
{ return dxyzdxi_map[p](2); }
Definition at line 775 of file fe_base.h.
References FEBase::dxyzdzeta_map.
Referenced by FEBase::compute_single_point_map().
{ return dxyzdzeta_map[p](2); }
Definition at line 539 of file fe_base.h.
References FEBase::curvatures.
{ return curvatures;}
Definition at line 300 of file fe_base.h.
References FEBase::calculate_d2phi, FEBase::calculations_started, and FEBase::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 FEBase::calculate_d2phi, FEBase::calculations_started, and FEBase::d2phidx2.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidx2; }
Definition at line 316 of file fe_base.h.
References FEBase::calculate_d2phi, FEBase::calculations_started, and FEBase::d2phidxdy.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidxdy; }
Definition at line 324 of file fe_base.h.
References FEBase::calculate_d2phi, FEBase::calculations_started, and FEBase::d2phidxdz.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidxdz; }
Definition at line 332 of file fe_base.h.
References FEBase::calculate_d2phi, FEBase::calculations_started, and FEBase::d2phidy2.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidy2; }
Definition at line 340 of file fe_base.h.
References FEBase::calculate_d2phi, FEBase::calculations_started, and FEBase::d2phidydz.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidydz; }
Definition at line 348 of file fe_base.h.
References FEBase::calculate_d2phi, FEBase::calculations_started, and FEBase::d2phidz2.
{ libmesh_assert(!calculations_started || calculate_d2phi);
calculate_d2phi = true; return d2phidz2; }
Definition at line 384 of file fe_base.h.
References FEBase::d2xyzdeta2_map.
{ return d2xyzdeta2_map; }
Definition at line 414 of file fe_base.h.
References FEBase::d2xyzdetadzeta_map.
{ return d2xyzdetadzeta_map; }
Definition at line 378 of file fe_base.h.
References FEBase::d2xyzdxi2_map.
{ return d2xyzdxi2_map; }
Definition at line 400 of file fe_base.h.
References FEBase::d2xyzdxideta_map.
{ return d2xyzdxideta_map; }
Definition at line 408 of file fe_base.h.
References FEBase::d2xyzdxidzeta_map.
{ return d2xyzdxidzeta_map; }
Definition at line 392 of file fe_base.h.
References FEBase::d2xyzdzeta2_map.
{ return d2xyzdzeta2_map; }
Definition at line 444 of file fe_base.h.
References FEBase::detadx_map.
{ return detadx_map; }
Definition at line 451 of file fe_base.h.
References FEBase::detady_map.
{ return detady_map; }
Definition at line 458 of file fe_base.h.
References FEBase::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 FEBase::dphase.
{ return dphase; }
Definition at line 242 of file fe_base.h.
References FEBase::calculate_dphi, FEBase::calculations_started, and FEBase::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 FEBase::calculate_dphi, FEBase::calculations_started, and FEBase::dphideta.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphideta; }
Definition at line 250 of file fe_base.h.
References FEBase::calculate_dphi, FEBase::calculations_started, and FEBase::dphidx.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidx; }
Definition at line 274 of file fe_base.h.
References FEBase::calculate_dphi, FEBase::calculations_started, and FEBase::dphidxi.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidxi; }
Definition at line 258 of file fe_base.h.
References FEBase::calculate_dphi, FEBase::calculations_started, and FEBase::dphidy.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidy; }
Definition at line 266 of file fe_base.h.
References FEBase::calculate_dphi, FEBase::calculations_started, and FEBase::dphidz.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidz; }
Definition at line 290 of file fe_base.h.
References FEBase::calculate_dphi, FEBase::calculations_started, and FEBase::dphidzeta.
{ libmesh_assert(!calculations_started || calculate_dphi);
calculate_dphi = true; return dphidzeta; }
Definition at line 423 of file fe_base.h.
References FEBase::dxidx_map.
{ return dxidx_map; }
Definition at line 430 of file fe_base.h.
References FEBase::dxidy_map.
{ return dxidy_map; }
Definition at line 437 of file fe_base.h.
References FEBase::dxidz_map.
{ return dxidz_map; }
Definition at line 365 of file fe_base.h.
References FEBase::dxyzdeta_map.
{ return dxyzdeta_map; }
Definition at line 358 of file fe_base.h.
References FEBase::dxyzdxi_map.
{ return dxyzdxi_map; }
Definition at line 372 of file fe_base.h.
References FEBase::dxyzdzeta_map.
{ return dxyzdzeta_map; }
Definition at line 465 of file fe_base.h.
References FEBase::dzetadx_map.
{ return dzetadx_map; }
Definition at line 472 of file fe_base.h.
References FEBase::dzetady_map.
{ return dzetady_map; }
Definition at line 479 of file fe_base.h.
References FEBase::dzetadz_map.
{ return dzetadz_map; }
Definition at line 598 of file fe_base.h.
References FEType::family, and FEBase::fe_type.
{ return fe_type.family; }
Definition at line 577 of file fe_base.h.
References FEBase::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 FEBase::JxW.
Referenced by ExactErrorEstimator::find_squared_element_error().
{ return JxW; }
Definition at line 533 of file fe_base.h.
References FEBase::normals.
{ return normals; }
Definition at line 582 of file fe_base.h.
References FEBase::_p_level, FEBase::fe_type, and FEType::order.
{ return static_cast<Order>(fe_type.order + _p_level); }
Definition at line 572 of file fe_base.h.
References FEBase::_p_level.
Referenced by REINIT_ERROR().
{ return _p_level; }
Definition at line 227 of file fe_base.h.
References FEBase::calculate_phi, FEBase::calculations_started, and FEBase::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 FEBase::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 FEBase::weight.
{ return weight; }
Definition at line 527 of file fe_base.h.
References FEBase::tangents.
{ return tangents; }
Definition at line 566 of file fe_base.h.
References FEBase::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 FEBase::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++;
}
On a p-refined element, o should be the total order of the element.
On a p-refined element, o should be the total order of the element.
On a p-refined element, o should be the total order of the element.
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; }
On a p-refined element, o should be the total order of the element.
Definition at line 195 of file fe.h.
{ return FE<Dim,T>::n_dofs (t,o); }
On a p-refined element, o should be the base order of the element.
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 FEBase::d2phi, and FEBase::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 FEBase::dphi.
Referenced by FEBase::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 FEBase::print_dphi(), FEBase::print_JxW(), FEBase::print_phi(), and FEBase::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 FEBase::JxW.
Referenced by FEBase::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 FEBase::phi.
Referenced by FEBase::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 FEBase::xyz.
Referenced by FEBase::print_info().
{
for (unsigned int i=0; i<xyz.size(); ++i)
os << xyz[i];
}
Definition at line 372 of file fe_map.C.
References FEBase::d2xyzdeta2_map, FEBase::d2xyzdetadzeta_map, FEBase::d2xyzdxi2_map, FEBase::d2xyzdxideta_map, FEBase::d2xyzdxidzeta_map, FEBase::d2xyzdzeta2_map, FEBase::detadx_map, FEBase::detady_map, FEBase::detadz_map, FEBase::dim, FEBase::dxidx_map, FEBase::dxidy_map, FEBase::dxidz_map, FEBase::dxyzdeta_map, FEBase::dxyzdxi_map, FEBase::dxyzdzeta_map, FEBase::dzetadx_map, FEBase::dzetady_map, FEBase::dzetadz_map, FEBase::JxW, and FEBase::xyz.
Referenced by FEBase::compute_affine_map(), and FEBase::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);
}
On a p-refined element, o should be the total order of the element.
On a p-refined element, o should be the base order of the element.
On a p-refined element, o should be the base order of the element.
On a p-refined element, o should be the total order of the element.
Note: Computing second derivatives is not currently supported for all element types: C1 (Clough and Hermite), Lagrange, Hierarchic, and Monomial are supported. All other element types return an error when asked for second derivatives.
On a p-refined element, o should be the base order of the element.
Note: Computing second derivatives is not currently supported for all element types: C1 (Clough and Hermite), Lagrange, Hierarchic, and Monomial are supported. All other element types return an error when asked for second derivatives.
On a p-refined element, o should be the total order of the element.
Reimplemented from FEBase.
Definition at line 424 of file fe.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 FEBase::get_order(), FEBase::get_p_level(), and REINIT_ERROR().
Definition at line 432 of file fe.h.
Definition at line 936 of file fe_base.h.
Referenced by FEBase::compute_map(), FEBase::compute_shape_functions(), FEBase::get_d2phi(), FEBase::get_d2phidx2(), FEBase::get_d2phidxdy(), FEBase::get_d2phidxdz(), FEBase::get_d2phidy2(), FEBase::get_d2phidydz(), and FEBase::get_d2phidz2().
Definition at line 931 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), FEBase::get_dphi(), FEBase::get_dphideta(), FEBase::get_dphidx(), FEBase::get_dphidxi(), FEBase::get_dphidy(), FEBase::get_dphidz(), and FEBase::get_dphidzeta().
Definition at line 926 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_phi().
Definition at line 921 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), FEBase::get_d2phi(), FEBase::get_d2phidx2(), FEBase::get_d2phidxdy(), FEBase::get_d2phidxdz(), FEBase::get_d2phidy2(), FEBase::get_d2phidydz(), FEBase::get_d2phidz2(), FEBase::get_dphi(), FEBase::get_dphideta(), FEBase::get_dphidx(), FEBase::get_dphidxi(), FEBase::get_dphidy(), FEBase::get_dphidz(), FEBase::get_dphidzeta(), and FEBase::get_phi().
Definition at line 1192 of file fe_base.h.
Referenced by FEBase::compute_edge_map(), FEBase::compute_face_map(), and FEBase::get_curvatures().
Definition at line 979 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), FEBase::get_d2phi(), and FEBase::print_d2phi().
Definition at line 999 of file fe_base.h.
Referenced by FEBase::compute_shape_functions().
Definition at line 1086 of file fe_base.h.
Referenced by FEBase::compute_single_point_map().
Definition at line 1004 of file fe_base.h.
Referenced by FEBase::compute_shape_functions().
Definition at line 1091 of file fe_base.h.
Referenced by FEBase::compute_single_point_map().
Definition at line 1014 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_d2phidx2().
Definition at line 1019 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_d2phidxdy().
Definition at line 1024 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_d2phidxdz().
Definition at line 984 of file fe_base.h.
Referenced by FEBase::compute_shape_functions().
Definition at line 1071 of file fe_base.h.
Referenced by FEBase::compute_single_point_map().
Definition at line 989 of file fe_base.h.
Referenced by FEBase::compute_shape_functions().
Definition at line 1076 of file fe_base.h.
Referenced by FEBase::compute_single_point_map().
Definition at line 994 of file fe_base.h.
Referenced by FEBase::compute_shape_functions().
Definition at line 1081 of file fe_base.h.
Referenced by FEBase::compute_single_point_map().
Definition at line 1029 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_d2phidy2().
Definition at line 1034 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_d2phidydz().
Definition at line 1039 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_d2phidz2().
Definition at line 1009 of file fe_base.h.
Referenced by FEBase::compute_shape_functions().
Definition at line 1096 of file fe_base.h.
Referenced by FEBase::compute_single_point_map().
Definition at line 1140 of file fe_base.h.
Referenced by FEBase::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 FEBase::compute_edge_map(), FEBase::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 FEBase::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 FEBase::compute_affine_map(), FEBase::compute_edge_map(), FEBase::compute_face_map(), FEBase::compute_single_point_map(), FEBase::get_d2xyzdeta2(), and FEBase::resize_map_vectors().
Definition at line 841 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_single_point_map(), FEBase::get_d2xyzdetadzeta(), and FEBase::resize_map_vectors().
Definition at line 815 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_edge_map(), FEBase::compute_face_map(), FEBase::compute_single_point_map(), FEBase::get_d2xyzdxi2(), and FEBase::resize_map_vectors().
Definition at line 821 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_edge_map(), FEBase::compute_face_map(), FEBase::compute_single_point_map(), FEBase::get_d2xyzdxideta(), and FEBase::resize_map_vectors().
Definition at line 835 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_single_point_map(), FEBase::get_d2xyzdxidzeta(), and FEBase::resize_map_vectors().
Definition at line 847 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_single_point_map(), FEBase::get_d2xyzdzeta2(), and FEBase::resize_map_vectors().
Definition at line 876 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), FEBase::get_detadx(), and FEBase::resize_map_vectors().
Definition at line 882 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), FEBase::get_detady(), and FEBase::resize_map_vectors().
Definition at line 888 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), FEBase::get_detadz(), and FEBase::resize_map_vectors().
Definition at line 784 of file fe_base.h.
Referenced by JumpErrorEstimator::coarse_n_flux_faces_increment(), FEBase::coarsened_dof_values(), FEBase::compute_affine_map(), FEBase::compute_edge_map(), FEBase::compute_face_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), JumpErrorEstimator::estimate_error(), and FEBase::resize_map_vectors().
Definition at line 1156 of file fe_base.h.
Referenced by FEBase::get_dphase().
Definition at line 941 of file fe_base.h.
Referenced by FEBase::compute_periodic_constraints(), FEBase::compute_proj_constraints(), FEBase::compute_shape_functions(), FEBase::get_dphi(), FEBase::print_d2phi(), and FEBase::print_dphi().
Definition at line 951 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_dphideta().
Definition at line 1059 of file fe_base.h.
Referenced by FEBase::compute_single_point_map().
Definition at line 961 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_dphidx().
Definition at line 946 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_dphidxi().
Definition at line 1054 of file fe_base.h.
Referenced by FEBase::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 FEBase::compute_shape_functions(), and FEBase::get_dphidy().
Definition at line 971 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_dphidz().
Definition at line 956 of file fe_base.h.
Referenced by FEBase::compute_shape_functions(), and FEBase::get_dphidzeta().
Definition at line 1064 of file fe_base.h.
Referenced by FEBase::compute_single_point_map().
Definition at line 1119 of file fe_base.h.
Referenced by FEBase::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 FEBase::compute_edge_map(), FEBase::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 FEBase::get_Sobolev_dweight().
Definition at line 855 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), FEBase::get_dxidx(), and FEBase::resize_map_vectors().
Definition at line 861 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), FEBase::get_dxidy(), and FEBase::resize_map_vectors().
Definition at line 867 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), FEBase::get_dxidz(), and FEBase::resize_map_vectors().
Definition at line 803 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_edge_map(), FEBase::compute_face_map(), FEBase::compute_single_point_map(), FEBase::dxdeta_map(), FEBase::dydeta_map(), FEBase::dzdeta_map(), FEBase::get_dxyzdeta(), and FEBase::resize_map_vectors().
Definition at line 797 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_edge_map(), FEBase::compute_face_map(), FEBase::compute_single_point_map(), FEBase::dxdxi_map(), FEBase::dydxi_map(), FEBase::dzdxi_map(), FEBase::get_dxyzdxi(), and FEBase::resize_map_vectors().
Definition at line 809 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_single_point_map(), FEBase::dxdzeta_map(), FEBase::dydzeta_map(), FEBase::dzdzeta_map(), FEBase::get_dxyzdzeta(), and FEBase::resize_map_vectors().
Definition at line 898 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), FEBase::get_dzetadx(), and FEBase::resize_map_vectors().
Definition at line 904 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), FEBase::get_dzetady(), and FEBase::resize_map_vectors().
Definition at line 910 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_shape_functions(), FEBase::compute_single_point_map(), FEBase::get_dzetadz(), and FEBase::resize_map_vectors().
Definition at line 1209 of file fe_base.h.
Referenced by FEBase::coarsened_dof_values(), FE< Dim, T >::edge_reinit(), FEBase::get_type(), InfFE< Dim, T_radial, T_map >::reinit(), and REINIT_ERROR().
Definition at line 1203 of file fe_base.h.
Referenced by FEBase::coarsened_dof_values(), JumpErrorEstimator::estimate_error(), FE< Dim, T >::FE(), FEBase::get_family(), FEBase::get_fe_type(), FEBase::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 FEBase::coarsened_dof_values(), FEBase::compute_affine_map(), FEBase::compute_edge_map(), FEBase::compute_face_map(), FEBase::compute_periodic_constraints(), FEBase::compute_proj_constraints(), FEBase::compute_single_point_map(), FE< Dim, T >::edge_reinit(), FEBase::get_JxW(), FEBase::print_JxW(), InfFE< Dim, T_radial, T_map >::reinit(), REINIT_ERROR(), and FEBase::resize_map_vectors().
Definition at line 437 of file fe.h.
Definition at line 437 of file fe.h.
Definition at line 1185 of file fe_base.h.
Referenced by FEBase::compute_edge_map(), FEBase::compute_face_map(), and FEBase::get_normals().
Definition at line 915 of file fe_base.h.
Referenced by FEBase::compute_periodic_constraints(), FEBase::compute_proj_constraints(), FEBase::get_phi(), and FEBase::print_phi().
Definition at line 1049 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::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 FEBase::compute_edge_map(), FEBase::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 FEBase::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 FEBase::compute_edge_map(), FEBase::compute_face_map(), and FEBase::get_tangents().
Definition at line 1170 of file fe_base.h.
Referenced by FEBase::get_Sobolev_weight().
Definition at line 789 of file fe_base.h.
Referenced by FEBase::compute_affine_map(), FEBase::compute_edge_map(), FEBase::compute_face_map(), FEBase::compute_single_point_map(), FE< Dim, T >::edge_reinit(), FEBase::get_xyz(), FEBase::print_xyz(), InfFE< Dim, T_radial, T_map >::reinit(), REINIT_ERROR(), and FEBase::resize_map_vectors().
Generated automatically by Doxygen for libMesh from the source code.