Loading...
Searching...
No Matches
lagrangebasis.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
#include <type_traits>
#include <dune/common/exceptions.hh>
#include <dune/localfunctions/lagrange/cache.hh>
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/leafprebasismixin.hh>
#include <dune/grid/common/capabilities.hh>
namespace Dune {
namespace Functions {
// *****************************************************************************
// This is the reusable part of the LagrangeBasis. It contains
//
// LagrangePreBasis
// LagrangeNode
//
// The pre-basis allows to create the others and is the owner of possible shared
// state. These components do _not_ depend on the global basis and local view
// and can be used without a global basis.
// *****************************************************************************
template<typename GV, int k, typename R=double>
class LagrangeNode;
template<typename GV, int k, typename R=double>
class LagrangePreBasis;
template<typename GV, int k, typename R>
class LagrangePreBasis :
public LeafPreBasisMixin< LagrangePreBasis<GV,k,R> >
{
static const int dim = GV::dimension;
static const bool useDynamicOrder = (k<0);
public:
using GridView = GV;
using size_type = std::size_t;
using Node = LagrangeNode<GV, k, R>;
LagrangePreBasis(const GridView& gv)
: LagrangePreBasis(gv, std::numeric_limits<unsigned int>::max())
{}
LagrangePreBasis(const GridView& gv, unsigned int order) :
gridView_(gv), order_(order)
{
if (!useDynamicOrder && order!=std::numeric_limits<unsigned int>::max())
DUNE_THROW(RangeError, "Template argument k has to be -1 when supplying a run-time order!");
for (int i=0; i<=dim; i++)
{
dofsPerCube_[i] = computeDofsPerCube(i);
dofsPerSimplex_[i] = computeDofsPerSimplex(i);
}
dofsPerPrism_ = computeDofsPerPrism();
dofsPerPyramid_ = computeDofsPerPyramid();
}
void initializeIndices()
{
vertexOffset_ = 0;
edgeOffset_ = vertexOffset_ + dofsPerCube(0) * ((size_type)gridView_.size(dim));
if (dim>=2)
{
triangleOffset_ = edgeOffset_ + dofsPerCube(1) * ((size_type) gridView_.size(dim-1));
quadrilateralOffset_ = triangleOffset_ + dofsPerSimplex(2) * ((size_type)gridView_.size(Dune::GeometryTypes::triangle));
}
if (dim==3) {
tetrahedronOffset_ = quadrilateralOffset_ + dofsPerCube(2) * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral));
prismOffset_ = tetrahedronOffset_ + dofsPerSimplex(3) * ((size_type)gridView_.size(Dune::GeometryTypes::tetrahedron));
hexahedronOffset_ = prismOffset_ + dofsPerPrism() * ((size_type)gridView_.size(Dune::GeometryTypes::prism));
pyramidOffset_ = hexahedronOffset_ + dofsPerCube(3) * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
}
}
const GridView& gridView() const
{
return gridView_;
}
void update (const GridView& gv)
{
gridView_ = gv;
}
Node makeNode() const
{
return Node{order_};
}
size_type dimension() const
{
switch (dim)
{
case 1:
return dofsPerCube(0) * ((size_type)gridView_.size(dim))
+ dofsPerCube(1) * ((size_type)gridView_.size(dim-1));
case 2:
{
return dofsPerCube(0) * ((size_type)gridView_.size(dim))
+ dofsPerCube(1) * ((size_type)gridView_.size(dim-1))
+ dofsPerSimplex(2) * ((size_type)gridView_.size(Dune::GeometryTypes::triangle))
+ dofsPerCube(2) * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral));
}
case 3:
{
return dofsPerCube(0) * ((size_type)gridView_.size(dim))
+ dofsPerCube(1) * ((size_type)gridView_.size(dim-1))
+ dofsPerSimplex(2) * ((size_type)gridView_.size(Dune::GeometryTypes::triangle))
+ dofsPerCube(2) * ((size_type)gridView_.size(Dune::GeometryTypes::quadrilateral))
+ dofsPerSimplex(3) * ((size_type)gridView_.size(Dune::GeometryTypes::tetrahedron))
+ dofsPerPyramid() * ((size_type)gridView_.size(Dune::GeometryTypes::pyramid))
+ dofsPerPrism() * ((size_type)gridView_.size(Dune::GeometryTypes::prism))
+ dofsPerCube(3) * ((size_type)gridView_.size(Dune::GeometryTypes::hexahedron));
}
}
DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
}
size_type maxNodeSize() const
{
// That cast to unsigned int is necessary because GV::dimension is an enum,
// which is not recognized by the power method as an integer type...
return power(order()+1, (unsigned int)GV::dimension);
}
template<typename It>
It indices(const Node& node, It it) const
{
for (size_type i = 0, end = node.finiteElement().size() ; i < end ; ++it, ++i)
{
Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
const auto& gridIndexSet = gridView().indexSet();
const auto& element = node.element();
// The dimension of the entity that the current dof is related to
auto dofDim = dim - localKey.codim();
// Test for a vertex dof
// The test for k==1 is redundant, but having it here allows the compiler to conclude
// at compile-time that the dofDim==0 case is the only one that will ever happen.
// This leads to measurable speed-up: see
// https://gitlab.dune-project.org/staging/dune-functions/issues/30
if (k==1 || dofDim==0) {
*it = {{ (size_type)(gridIndexSet.subIndex(element,localKey.subEntity(),dim)) }};
continue;
}
if (dofDim==1)
{ // edge dof
if (dim==1) // element dof -- any local numbering is fine
{
*it = {{ edgeOffset_
+ dofsPerCube(1) * ((size_type)gridIndexSet.subIndex(element,0,0))
+ localKey.index() }};
continue;
}
else
{
const auto refElement
= Dune::referenceElement<double,dim>(element.type());
// We have to reverse the numbering if the local element edge is
// not aligned with the global edge.
auto v0 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim);
auto v1 = (size_type)gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim);
bool flip = (v0 > v1);
*it = {{ (flip)
? edgeOffset_
+ dofsPerCube(1)*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
+ (dofsPerCube(1)-1)-localKey.index()
: edgeOffset_
+ dofsPerCube(1)*((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim()))
+ localKey.index() }};
continue;
}
}
if (dofDim==2)
{
if (dim==2) // element dof -- any local numbering is fine
{
if (element.type().isTriangle())
{
*it = {{ triangleOffset_ + dofsPerSimplex(2)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
continue;
}
else if (element.type().isQuadrilateral())
{
*it = {{ quadrilateralOffset_ + dofsPerCube(2)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
continue;
}
else
DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals");
} else
{
const auto refElement
= Dune::referenceElement<double,dim>(element.type());
if (order()>3)
DUNE_THROW(Dune::NotImplemented, "LagrangeBasis for 3D grids is only implemented if k<=3");
if (order()==3 and !refElement.type(localKey.subEntity(), localKey.codim()).isTriangle())
DUNE_THROW(Dune::NotImplemented, "LagrangeBasis for 3D grids with k==3 is only implemented if the grid is a simplex grid");
*it = {{ triangleOffset_ + ((size_type)gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())) }};
continue;
}
}
if (dofDim==3)
{
if (dim==3) // element dof -- any local numbering is fine
{
if (element.type().isTetrahedron())
{
*it = {{ tetrahedronOffset_ + dofsPerSimplex(3)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
continue;
}
else if (element.type().isHexahedron())
{
*it = {{ hexahedronOffset_ + dofsPerCube(3)*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
continue;
}
else if (element.type().isPrism())
{
*it = {{ prismOffset_ + dofsPerPrism()*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
continue;
}
else if (element.type().isPyramid())
{
*it = {{ pyramidOffset_ + dofsPerPyramid()*((size_type)gridIndexSet.subIndex(element,0,0)) + localKey.index() }};
continue;
}
else
DUNE_THROW(Dune::NotImplemented, "3d elements have to be tetrahedra, hexahedra, prisms, or pyramids");
} else
DUNE_THROW(Dune::NotImplemented, "Grids of dimension larger than 3 are no supported");
}
DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the LagrangeBasis");
}
return it;
}
unsigned int order() const
{
return (useDynamicOrder) ? order_ : k;
}
protected:
GridView gridView_;
// Run-time order, only valid if k<0
unsigned int order_;
size_type dofsPerSimplex(std::size_t simplexDim) const
{
return useDynamicOrder ? dofsPerSimplex_[simplexDim] : computeDofsPerSimplex(simplexDim);
}
size_type dofsPerCube(std::size_t cubeDim) const
{
return useDynamicOrder ? dofsPerCube_[cubeDim] : computeDofsPerCube(cubeDim);
}
size_type dofsPerPrism() const
{
return useDynamicOrder ? dofsPerPrism_ : computeDofsPerPrism();
}
size_type dofsPerPyramid() const
{
return useDynamicOrder ? dofsPerPyramid_ : computeDofsPerPyramid();
}
size_type computeDofsPerSimplex(std::size_t simplexDim) const
{
return order() == 0 ? (dim == simplexDim ? 1 : 0) : Dune::binomial(std::size_t(order()-1),simplexDim);
}
size_type computeDofsPerCube(std::size_t cubeDim) const
{
return order() == 0 ? (dim == cubeDim ? 1 : 0) : Dune::power(order()-1, cubeDim);
}
size_type computeDofsPerPrism() const
{
return order() == 0 ? (dim == 3 ? 1 : 0) : (order()-1)*(order()-1)*(order()-2)/2;
}
size_type computeDofsPerPyramid() const
{
return order() == 0 ? (dim == 3 ? 1 : 0) : (order()-2)*(order()-1)*(2*order()-3)/6;
}
// When the order is given at run-time, the following numbers are pre-computed:
std::array<size_type,dim+1> dofsPerSimplex_;
std::array<size_type,dim+1> dofsPerCube_;
size_type dofsPerPrism_;
size_type dofsPerPyramid_;
size_type vertexOffset_;
size_type edgeOffset_;
size_type triangleOffset_;
size_type quadrilateralOffset_;
size_type tetrahedronOffset_;
size_type pyramidOffset_;
size_type prismOffset_;
size_type hexahedronOffset_;
};
template<typename GV, int k, typename R>
class LagrangeNode :
public LeafBasisNode
{
static constexpr int dim = GV::dimension;
static constexpr bool useDynamicOrder = (k<0);
// Compute the GeometryType id in case the grid has only a single GeometryType
static constexpr GeometryType::Id geometryTypeId()
{
if constexpr(Dune::Capabilities::hasSingleGeometryType<typename GV::Grid>::v)
return GeometryType(Dune::Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId, GV::dimension);
else
return GeometryType::Id(~0u);
}
// Select the static LFECache if k >= 0, else the dynamic LFECache
using FiniteElementCache = std::conditional_t<(useDynamicOrder),
DynamicLagrangeLocalFiniteElementCache<typename GV::ctype,R,dim>,
StaticLagrangeLocalFiniteElementCache<geometryTypeId(),typename GV::ctype,R,dim,std::max(k,0)>
>;
// Construct the FiniteElementCache depending on whether the order is dynamic or static
static auto makeFiniteElementCache(unsigned int order)
{
if constexpr (useDynamicOrder)
return FiniteElementCache{order};
else
return FiniteElementCache{};
}
public:
using size_type = std::size_t;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElement = typename FiniteElementCache::FiniteElementType;
LagrangeNode() :
LagrangeNode(k)
{}
LagrangeNode(unsigned int order) :
order_(order),
cache_(makeFiniteElementCache(order)),
finiteElement_(nullptr),
element_(nullptr)
{}
const Element& element() const
{
return *element_;
}
const FiniteElement& finiteElement() const
{
return *finiteElement_;
}
void bind(const Element& e)
{
element_ = &e;
finiteElement_ = &(cache_.get(element_->type()));
this->setSize(finiteElement_->size());
}
protected:
unsigned int order() const
{
return order_;
}
// Run-time order, only valid if k<0
unsigned int order_;
FiniteElementCache cache_;
const FiniteElement* finiteElement_;
const Element* element_;
};
namespace BasisFactory {
template<std::size_t k, typename R=double>
auto lagrange()
{
return [](const auto& gridView) {
return LagrangePreBasis<std::decay_t<decltype(gridView)>, k, R>(gridView);
};
}
template<typename R=double>
auto lagrange(int order)
{
return [=](const auto& gridView) {
return LagrangePreBasis<std::decay_t<decltype(gridView)>, -1, R>(gridView, order);
};
}
} // end namespace BasisFactory
template<typename GV, int k=-1, typename R=double>
using LagrangeBasis = DefaultGlobalBasis<LagrangePreBasis<GV, k, R> >;
} // end namespace Functions
} // end namespace Dune
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
nodes.hh
defaultglobalbasis.hh
leafprebasismixin.hh
Dune::Functions::BasisFactory::lagrange
auto lagrange()
Create a pre-basis factory that can create a Lagrange pre-basis.
Definition lagrangebasis.hh:477
Dune::Functions::LagrangeBasis
DefaultGlobalBasis< LagrangePreBasis< GV, k, R > > LagrangeBasis
Nodal basis of a scalar k-th-order Lagrangean finite element space.
Definition lagrangebasis.hh:532
Dune
Definition monomialset.hh:19
Dune::Functions
Definition monomialset.hh:19
Dune::Functions::BasisFactory
Definition argyrisbasis.hh:926
Dune::Functions::DefaultGlobalBasis
Global basis for given pre-basis.
Definition defaultglobalbasis.hh:53
Dune::Functions::LagrangeNode
Definition lagrangebasis.hh:379
Dune::Functions::LagrangeNode::LagrangeNode
LagrangeNode(unsigned int order)
Constructor with a run-time order.
Definition lagrangebasis.hh:419
Dune::Functions::LagrangeNode::order
unsigned int order() const
Definition lagrangebasis.hh:451
Dune::Functions::LagrangeNode::finiteElement
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition lagrangebasis.hh:436
Dune::Functions::LagrangeNode< GV, 2 >::element_
const Element * element_
Definition lagrangebasis.hh:461
Dune::Functions::LagrangeNode::element
const Element & element() const
Return current element, throw if unbound.
Definition lagrangebasis.hh:427
Dune::Functions::LagrangeNode< GV, 2 >::cache_
FiniteElementCache cache_
Definition lagrangebasis.hh:459
Dune::Functions::LagrangeNode::FiniteElement
typename FiniteElementCache::FiniteElementType FiniteElement
Definition lagrangebasis.hh:411
Dune::Functions::LagrangeNode::bind
void bind(const Element &e)
Bind to element.
Definition lagrangebasis.hh:442
Dune::Functions::LagrangeNode::Element
typename GV::template Codim< 0 >::Entity Element
Definition lagrangebasis.hh:410
Dune::Functions::LagrangeNode< GV, 2 >::finiteElement_
const FiniteElement * finiteElement_
Definition lagrangebasis.hh:460
Dune::Functions::LagrangeNode< GV, 2 >::order_
unsigned int order_
Definition lagrangebasis.hh:457
Dune::Functions::LagrangeNode::size_type
std::size_t size_type
Definition lagrangebasis.hh:409
Dune::Functions::LagrangeNode::LagrangeNode
LagrangeNode()
Constructor without order (uses the compile-time value).
Definition lagrangebasis.hh:414
Dune::Functions::LagrangePreBasis
A pre-basis for a PQ-lagrange bases with given order.
Definition lagrangebasis.hh:65
Dune::Functions::LagrangePreBasis::dofsPerPrism
size_type dofsPerPrism() const
Definition lagrangebasis.hh:325
Dune::Functions::LagrangePreBasis::computeDofsPerCube
size_type computeDofsPerCube(std::size_t cubeDim) const
Number of degrees of freedom assigned to a cube (without the ones assigned to its faces!...
Definition lagrangebasis.hh:342
Dune::Functions::LagrangePreBasis::computeDofsPerSimplex
size_type computeDofsPerSimplex(std::size_t simplexDim) const
Number of degrees of freedom assigned to a simplex (without the ones assigned to its faces!...
Definition lagrangebasis.hh:336
Dune::Functions::LagrangePreBasis::computeDofsPerPrism
size_type computeDofsPerPrism() const
Definition lagrangebasis.hh:347
Dune::Functions::LagrangePreBasis< GV, 1 >::order_
unsigned int order_
Definition lagrangebasis.hh:311
Dune::Functions::LagrangePreBasis< GV, 1 >::edgeOffset_
size_type edgeOffset_
Definition lagrangebasis.hh:364
Dune::Functions::LagrangePreBasis::Node
LagrangeNode< GV, k, R > Node
Template mapping root tree path to type of created tree node.
Definition lagrangebasis.hh:78
Dune::Functions::LagrangePreBasis< GV, 1 >::dofsPerSimplex_
std::array< size_type, dim+1 > dofsPerSimplex_
Definition lagrangebasis.hh:358
Dune::Functions::LagrangePreBasis::indices
It indices(const Node &node, It it) const
Definition lagrangebasis.hh:184
Dune::Functions::LagrangePreBasis< GV, 1 >::vertexOffset_
size_type vertexOffset_
Definition lagrangebasis.hh:363
Dune::Functions::LagrangePreBasis::dofsPerSimplex
size_type dofsPerSimplex(std::size_t simplexDim) const
Number of degrees of freedom assigned to a simplex (without the ones assigned to its faces!...
Definition lagrangebasis.hh:314
Dune::Functions::LagrangePreBasis::size_type
std::size_t size_type
Type used for indices and size information.
Definition lagrangebasis.hh:75
Dune::Functions::LagrangePreBasis< GV, 1 >::pyramidOffset_
size_type pyramidOffset_
Definition lagrangebasis.hh:368
Dune::Functions::LagrangePreBasis< GV, 1 >::prismOffset_
size_type prismOffset_
Definition lagrangebasis.hh:369
Dune::Functions::LagrangePreBasis< GV, 1 >::tetrahedronOffset_
size_type tetrahedronOffset_
Definition lagrangebasis.hh:367
Dune::Functions::LagrangePreBasis::initializeIndices
void initializeIndices()
Initialize the global indices.
Definition lagrangebasis.hh:102
Dune::Functions::LagrangePreBasis::computeDofsPerPyramid
size_type computeDofsPerPyramid() const
Definition lagrangebasis.hh:352
Dune::Functions::LagrangePreBasis::LagrangePreBasis
LagrangePreBasis(const GridView &gv, unsigned int order)
Constructor for a given grid view object and run-time order.
Definition lagrangebasis.hh:86
Dune::Functions::LagrangePreBasis< GV, 1 >::dofsPerPyramid_
size_type dofsPerPyramid_
Definition lagrangebasis.hh:361
Dune::Functions::LagrangePreBasis::dimension
size_type dimension() const
Get the total dimension of the space spanned by this basis.
Definition lagrangebasis.hh:146
Dune::Functions::LagrangePreBasis::LagrangePreBasis
LagrangePreBasis(const GridView &gv)
Constructor for a given grid view object with compile-time order.
Definition lagrangebasis.hh:81
Dune::Functions::LagrangePreBasis::update
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition lagrangebasis.hh:132
Dune::Functions::LagrangePreBasis::GridView
GV GridView
The grid view that the FE basis is defined on.
Definition lagrangebasis.hh:72
Dune::Functions::LagrangePreBasis< GV, 1 >::quadrilateralOffset_
size_type quadrilateralOffset_
Definition lagrangebasis.hh:366
Dune::Functions::LagrangePreBasis::gridView
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition lagrangebasis.hh:126
Dune::Functions::LagrangePreBasis< GV, 1 >::gridView_
GridView gridView_
Definition lagrangebasis.hh:308
Dune::Functions::LagrangePreBasis::makeNode
Node makeNode() const
Create tree node.
Definition lagrangebasis.hh:140
Dune::Functions::LagrangePreBasis< GV, 1 >::order
unsigned int order() const
Definition lagrangebasis.hh:302
Dune::Functions::LagrangePreBasis< GV, 1 >::dofsPerCube_
std::array< size_type, dim+1 > dofsPerCube_
Definition lagrangebasis.hh:359
Dune::Functions::LagrangePreBasis::maxNodeSize
size_type maxNodeSize() const
Get the maximal number of DOFs associated to node for any element.
Definition lagrangebasis.hh:176
Dune::Functions::LagrangePreBasis< GV, 1 >::dofsPerPrism_
size_type dofsPerPrism_
Definition lagrangebasis.hh:360
Dune::Functions::LagrangePreBasis::dofsPerCube
size_type dofsPerCube(std::size_t cubeDim) const
Number of degrees of freedom assigned to a cube (without the ones assigned to its faces!...
Definition lagrangebasis.hh:320
Dune::Functions::LagrangePreBasis< GV, 1 >::triangleOffset_
size_type triangleOffset_
Definition lagrangebasis.hh:365
Dune::Functions::LagrangePreBasis< GV, 1 >::hexahedronOffset_
size_type hexahedronOffset_
Definition lagrangebasis.hh:370
Dune::Functions::LagrangePreBasis::dofsPerPyramid
size_type dofsPerPyramid() const
Definition lagrangebasis.hh:330
Dune::Functions::LeafPreBasisMixin
A generic MixIn class for PreBasis.
Definition leafprebasismixin.hh:36
Dune::Functions::BasisNodeMixin::setSize
void setSize(const size_type size)
Definition nodes.hh:197
Dune::Functions::LeafBasisNode
Definition nodes.hh:218