Forked from
agnumpde / dune-fufem
367 commits behind the upstream repository.
-
Jonathan Youett authoredJonathan Youett authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
p3nodalbasis.hh 8.65 KiB
#ifndef P3_NODALBASIS_HH
#define P3_NODALBASIS_HH
/*
#warning This file is deprecated. All implementations of function space bases in dune-fufem \
are in the process of being replaced by counterparts in the new dune-functions module. \
Those are syntactically different, but semantically very close to the dune-fufem implementations. \
To get rid of this warning, replace all occurrences of the Q1NodalBasis<...> class in your code \
by DuneFunctionsBasis<Dune::Functions::PQkNodalBasis<...,3> >.
*/
/**
@file
@brief The nodal basis for the space of 3rd-order Lagrangian finite elements
@author Oliver Sander
*/
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/common/typetraits.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh>
#include <dune/fufem/functionspacebases/functionspacebasis.hh>
/** \brief The nodal basis for the space of 3rd-order Lagrangian finite elements
*
* The main task of this class is to assign global indices to the local degrees of
* freedom of a given element. A P3 Lagrangian FE space has
* - 1 dof per vertex
* - 2 dofs per edge
* - 1 dof per triangle
* - 4 dofs per quadrilateral
* - 0 dofs per tetrahedron
* - 1 dof per pyramid
* - 2 dofs per prism
* - 8 dofs per hexahedron
*
* \tparam GV Grid view that the FE space is defined on
* \tparam RT Value type of the shape functions
*/
template <class GV, class RT=double>
class P3NodalBasis :
public FunctionSpaceBasis<
GV,
RT,
typename Dune::PQkLocalFiniteElementCache<typename GV::Grid::ctype, RT, GV::dimension, 3>::FiniteElementType >
{
protected:
typedef typename GV::Grid::ctype ctype;
typedef typename Dune::PQkLocalFiniteElementCache<typename GV::Grid::ctype, RT, GV::dimension, 3> FiniteElementCache;
typedef typename FiniteElementCache::FiniteElementType LFE;
typedef FunctionSpaceBasis<GV, RT, LFE> Base;
typedef typename Base::Element Element;
/** \brief Dimension of the grid */
using Base::dim;
using Base::gridview_;
public:
typedef typename Base::GridView GridView;
typedef typename Base::ReturnType ReturnType;
typedef typename Base::LocalFiniteElement LocalFiniteElement;
typedef typename Base::LinearCombination LinearCombination;
/** \brief Constructor for a given grid view */ P3NodalBasis(const GridView& gridview) :
Base(gridview)
{
updateMapper();
}
/** \brief The number of degrees of freedom of the function space */
size_t size() const
{
switch (dim) {
case 1:
// One for each vertex, and two for each element
return gridview_.size(1) + 2*gridview_.size(0);
case 2: {
Dune::GeometryType triangle = Dune::GeometryTypes::triangle;
Dune::GeometryType quad = Dune::GeometryTypes::quadrilateral;
// One for each vertex, two for each edge,
// one for each triangle element, four for each quad element
return gridview_.size(dim) + 2*gridview_.size(1)
+ gridview_.size(triangle) + 4*gridview_.size(quad);
}
case 3: {
Dune::GeometryType triangle = Dune::GeometryTypes::triangle;
Dune::GeometryType quad = Dune::GeometryTypes::quadrilateral;
Dune::GeometryType pyramid = Dune::GeometryTypes::pyramid;
Dune::GeometryType prism = Dune::GeometryTypes::prism;
Dune::GeometryType hexahedron = Dune::GeometryTypes::cube(3);
// One for each vertex, two for each edge,
// one for each triangle element, four for each quad element
return gridview_.size(dim) + 2*gridview_.size(2)
+ gridview_.size(triangle) + 4*gridview_.size(quad)
+ gridview_.size(pyramid) + 2*gridview_.size(prism) + 8*gridview_.size(hexahedron);
}
default:
DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
}
return 0;
}
void update(const GridView& gridview)
{
Base::update(gridview);
updateMapper();
}
/** \brief Get the local finite element for a given grid element */
const LocalFiniteElement& getLocalFiniteElement(const Element& e) const
{
return cache_.get(e.type());
}
/** \brief Get the global index for the i-th local index on the given element */
int index(const Element& e, const int i) const
{
Dune::LocalKey localKey = cache_.get(e.type()).localCoefficients().localKey(i);
const typename GridView::IndexSet& indexSet = gridview_.indexSet();
// The dimension of the entity that the current dof is related to
size_t dofDim = dim - localKey.codim();
if (dofDim==0) { // vertex dof
return indexSet.subIndex(e,localKey.subEntity(),dim);
}
if (dofDim==1) { // edge dof
if (dim==1) // element dof -- any local numbering is fine
return edgeOffset_ + 2*indexSet.subIndex(e,0,0) + localKey.index();
else {
const auto& refElement = Dune::ReferenceElements<double,dim>::general(e.type());
// we have to reverse the numbering if the local triangle edge is
// not aligned with the global edge
//
// This use of ReferenceElement::subEntity is OK, because
// the local finite element is defined on the reference element
// and hence the local keys correspond to the subentity embedding
// of the reference element.
size_t v0 = indexSet.subIndex(e,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim);
size_t v1 = indexSet.subIndex(e,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim);
bool flip = (v0 > v1);
return (flip)
? edgeOffset_ + 2*indexSet.subIndex(e,localKey.subEntity(),localKey.codim()) + 1-localKey.index()
: edgeOffset_ + 2*indexSet.subIndex(e,localKey.subEntity(),localKey.codim()) + localKey.index();
}
}
if (dofDim==2) {
if (dim==2) { // element dof -- any local numbering is fine
if (e.type().isTriangle())
return triangleOffset_ + 1*indexSet.subIndex(e,0,0) + localKey.index();
else if (e.type().isQuadrilateral())
return quadrilateralOffset_ + 4*indexSet.subIndex(e,0,0) + localKey.index();
else
DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
} else {
const auto& refElement = Dune::ReferenceElements<double,dim>::general(e.type());
assert(refElement.type(localKey.subEntity(), localKey.codim()).isTriangle());
return triangleOffset_ + indexSet.subIndex(e,localKey.subEntity(),localKey.codim());
}
}
return 0;
}
protected:
void updateMapper()
{
vertexOffset_ = 0;
edgeOffset_ = vertexOffset_ + gridview_.size(dim);
triangleOffset_ = edgeOffset_ + 2*gridview_.size(dim-1);
quadrilateralOffset_ = triangleOffset_ + 1*gridview_.size(Dune::GeometryTypes::triangle);
if (dim==3) {
tetrahedronOffset_ = quadrilateralOffset_ + 4*gridview_.size(Dune::GeometryTypes::quadrilateral);
pyramidOffset_ = tetrahedronOffset_ + 0*gridview_.size(Dune::GeometryTypes::tetrahedron);
prismOffset_ = tetrahedronOffset_ + 1*gridview_.size(Dune::GeometryTypes::pyramid);
hexahedronOffset_ = tetrahedronOffset_ + 2*gridview_.size(Dune::GeometryTypes::prism);
}
}
FiniteElementCache cache_;
size_t vertexOffset_;
size_t edgeOffset_;
size_t triangleOffset_;
size_t quadrilateralOffset_;
size_t tetrahedronOffset_;
size_t pyramidOffset_;
size_t prismOffset_;
size_t hexahedronOffset_;
};
#endif