Loading...
Searching...
No Matches
nedelecbasis.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_NEDELECBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
#include <array>
#include <dune/common/exceptions.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/localfunctions/common/localfiniteelementvariant.hh>
#include <dune/localfunctions/nedelec.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/globalvaluedlocalfiniteelement.hh>
#include <dune/functions/functionspacebases/leafprebasismappermixin.hh>
#include <dune/functions/functionspacebases/nodes.hh>
namespace Dune::Functions
{
namespace Impl
{
template<typename GV, int dim, typename R, std::size_t order>
class Nedelec1stKindLocalFiniteElementMap
{
using D = typename GV::ctype;
constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
using CubeFiniteElement = Nedelec1stKindCubeLocalFiniteElement<D,R,dim,order>;
using SimplexFiniteElement = Nedelec1stKindSimplexLocalFiniteElement<D,R,dim,order>;
public:
using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId; // meaningless if hasFixedElementType is false
constexpr static GeometryType type = GeometryType(topologyId, GV::dimension);
using FiniteElement = std::conditional_t<hasFixedElementType,
std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
static std::size_t numVariants(GeometryType type)
{
if (order!=1) // I am not sure whether the formula below is correct for all orders.
DUNE_THROW(NotImplemented, "Only Nedelec elements of order 1 are implemented!");
auto numEdges = referenceElement<D,dim>(type).size(dim-1);
return power(2,numEdges);
}
Nedelec1stKindLocalFiniteElementMap(const GV& gv)
: elementMapper_(gv, mcmgElementLayout()),
orientation_(gv.size(0))
{
update(gv);
}
void update(const GV& gv)
{
elementMapper_.update(gv);
orientation_.resize(gv.size(0));
// create all variants
if constexpr (hasFixedElementType)
{
variants_.resize(numVariants(type));
for (size_t i = 0; i < numVariants(type); i++)
variants_[i] = FiniteElement(i);
}
else
{
// for mixed grids add offset for cubes
variants_.resize(numVariants(GeometryTypes::simplex(dim)) + numVariants(GeometryTypes::cube(dim)));
for (size_t i = 0; i < numVariants(GeometryTypes::simplex(dim)); i++)
variants_[i] = SimplexFiniteElement(i);
for (size_t i = 0; i < numVariants(GeometryTypes::cube(dim)); i++)
variants_[i + numVariants(GeometryTypes::simplex(dim))] = CubeFiniteElement(i);
}
// compute orientation for all elements
const auto& indexSet = gv.indexSet();
for(const auto& element : elements(gv))
{
const auto& refElement = referenceElement(element);
auto elementIndex = elementMapper_.index(element);
orientation_[elementIndex] = 0;
for (std::size_t i=0; i<element.subEntities(dim-1); i++)
{
// Local vertex indices within the element
auto localV0 = refElement.subEntity(i,dim-1, 0,dim);
auto localV1 = refElement.subEntity(i,dim-1, 1,dim);
// Global vertex indices within the grid
auto globalV0 = indexSet.subIndex(element,localV0,dim);
auto globalV1 = indexSet.subIndex(element,localV1,dim);
if ( (localV0<localV1 && globalV0>globalV1) || (localV0>localV1 && globalV0<globalV1) )
orientation_[elementIndex] |= (1 << i);
}
// for mixed grids add offset for cubes
if constexpr (!hasFixedElementType)
if (element.type().isCube())
orientation_[elementIndex] += numVariants(GeometryTypes::simplex(dim));
}
}
template<class Element>
const auto& find(const Element& element) const
{
return variants_[orientation_[elementMapper_.index(element)]];
}
private:
std::vector<FiniteElement> variants_;
Dune::MultipleCodimMultipleGeomTypeMapper<GV> elementMapper_;
std::vector<unsigned short> orientation_;
};
} // namespace Impl
// *****************************************************************************
// This is the reusable part of the basis. It contains
//
// NedelecPreBasis
// NedelecNode
//
// 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, typename Range, std::size_t kind, int order>
class NedelecNode;
template<typename GV, typename Range, std::size_t kind, int order>
class NedelecPreBasis :
public LeafPreBasisMapperMixin<GV>
{
static const int dim = GV::dimension;
static_assert(kind==1, "Only the Nedelec basis of the first kind is currently implemented!");
using FiniteElementMap = typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
public:
using GridView = GV;
using size_type = std::size_t;
using Node = NedelecNode<GV, Range, kind, order>;
NedelecPreBasis(const GridView& gv)
: LeafPreBasisMapperMixin<GV>(gv, mcmgLayout(Dim<1>{}))
, finiteElementMap_(gv)
{
if (kind!=1)
DUNE_THROW(NotImplemented, "Only Nedelec elements of the first kind are implemented!");
// There is no inherent reason why the basis shouldn't work for grids with more than two
// element types. Somebody simply has to sit down and implement the missing bits.
if (gv.indexSet().types(0).size() > 2)
DUNE_THROW(NotImplemented, "Nédélec basis is only implemented for grids with simplex and cube elements");
for(auto type : gv.indexSet().types(0))
if (!type.isSimplex() && !type.isCube())
DUNE_THROW(NotImplemented, "Nédélec basis is only implemented for grids with simplex or cube elements.");
if (order>1)
DUNE_THROW(NotImplemented, "Only first-order elements are implemented");
if (dim!=2 && dim!=3)
DUNE_THROW(NotImplemented, "Only 2d and 3d Nédélec elements are implemented");
}
void update (const GridView& gv)
{
LeafPreBasisMapperMixin<GV>::update(gv);
finiteElementMap_.update(gv);
}
Node makeNode() const
{
return Node{&finiteElementMap_};
}
protected:
FiniteElementMap finiteElementMap_;
};
template<typename GV, typename Range, size_t kind, int order>
class NedelecNode :
public LeafBasisNode
{
static const int dim = GV::dimension;
public:
using size_type = std::size_t;
using Element = typename GV::template Codim<0>::Entity;
static_assert(kind==1, "Only Nedelec elements of the first kind are implemented!");
using FiniteElementMap = typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::CovariantPiolaTransformator,
typename FiniteElementMap::FiniteElement,
Element>;
NedelecNode(const FiniteElementMap* finiteElementMap) :
element_(nullptr),
finiteElementMap_(finiteElementMap)
{ }
const Element& element() const
{
return *element_;
}
const FiniteElement& finiteElement() const
{
return finiteElement_;
}
void bind(const Element& e)
{
element_ = &e;
finiteElement_.bind((finiteElementMap_->find(*element_)), e);
this->setSize(finiteElement_.size());
}
protected:
FiniteElement finiteElement_;
const Element* element_;
const FiniteElementMap* finiteElementMap_;
};
namespace BasisFactory {
template<std::size_t kind, std::size_t order, typename Range=double>
auto nedelec()
{
return [](const auto& gridView) {
return NedelecPreBasis<std::decay_t<decltype(gridView)>, Range, kind, order>(gridView);
};
}
} // end namespace BasisFactory
// *****************************************************************************
// This is the actual global basis implementation based on the reusable parts.
// *****************************************************************************
template<typename GV, std::size_t kind, std::size_t order, typename Range=double>
using NedelecBasis = DefaultGlobalBasis<NedelecPreBasis<GV, Range, kind, order > >;
} // end namespace Dune::Functions
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
nodes.hh
globalvaluedlocalfiniteelement.hh
leafprebasismappermixin.hh
defaultglobalbasis.hh
Dune::Functions::BasisFactory::power
auto power(ChildPreBasisFactory &&childPreBasisFactory, std::size_t k, const IndexMergingStrategy &)
Create a pre-basis factory that can build a PowerPreBasis.
Definition dynamicpowerbasis.hh:410
Dune::Functions::BasisFactory::nedelec
auto nedelec()
Create a pre-basis factory that can create a Nédélec pre-basis.
Definition nedelecbasis.hh:274
Dune::Functions::Experimental::elements
auto elements(const SubDomainGridView< HostGridView > &subDomainGridView)
ADL findable access to element range for a SubDomainGridView.
Definition subdomain.hh:487
Dune::Functions
Definition monomialset.hh:19
Dune::Functions::NedelecBasis
DefaultGlobalBasis< NedelecPreBasis< GV, Range, kind, order > > NedelecBasis
Basis of a k-th-order Nédélec finite element space.
Definition nedelecbasis.hh:297
Dune::Functions::BasisFactory
Definition argyrisbasis.hh:926
Dune::Functions::DefaultGlobalBasis
Global basis for given pre-basis.
Definition defaultglobalbasis.hh:53
Dune::Functions::LeafPreBasisMapperMixin::update
void update(const GridView &gv)
Update the stored GridView.
Definition leafprebasismappermixin.hh:101
Dune::Functions::LeafPreBasisMapperMixin::LeafPreBasisMapperMixin
LeafPreBasisMapperMixin(const GridView &gv, Dune::MCMGLayout layout)
Construct from GridView and local DOF layout.
Definition leafprebasismappermixin.hh:74
Dune::Functions::NedelecNode
Definition nedelecbasis.hh:212
Dune::Functions::NedelecNode::finiteElementMap_
const FiniteElementMap * finiteElementMap_
Definition nedelecbasis.hh:257
Dune::Functions::NedelecNode::finiteElement_
FiniteElement finiteElement_
Definition nedelecbasis.hh:255
Dune::Functions::NedelecNode::FiniteElement
Impl::GlobalValuedLocalFiniteElement< Impl::CovariantPiolaTransformator, typename FiniteElementMap::FiniteElement, Element > FiniteElement
Definition nedelecbasis.hh:221
Dune::Functions::NedelecNode::bind
void bind(const Element &e)
Bind to element.
Definition nedelecbasis.hh:246
Dune::Functions::NedelecNode::element
const Element & element() const
Return current element, throw if unbound.
Definition nedelecbasis.hh:231
Dune::Functions::NedelecNode::Element
typename GV::template Codim< 0 >::Entity Element
Definition nedelecbasis.hh:218
Dune::Functions::NedelecNode::finiteElement
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition nedelecbasis.hh:240
Dune::Functions::NedelecNode::NedelecNode
NedelecNode(const FiniteElementMap *finiteElementMap)
Definition nedelecbasis.hh:225
Dune::Functions::NedelecNode::FiniteElementMap
typename Impl::Nedelec1stKindLocalFiniteElementMap< GV, dim, Range, order > FiniteElementMap
Definition nedelecbasis.hh:220
Dune::Functions::NedelecNode::element_
const Element * element_
Definition nedelecbasis.hh:256
Dune::Functions::NedelecNode::size_type
std::size_t size_type
Definition nedelecbasis.hh:217
Dune::Functions::NedelecPreBasis
Definition nedelecbasis.hh:150
Dune::Functions::NedelecPreBasis::size_type
std::size_t size_type
Definition nedelecbasis.hh:160
Dune::Functions::NedelecPreBasis::NedelecPreBasis
NedelecPreBasis(const GridView &gv)
Constructor for a given grid view object.
Definition nedelecbasis.hh:165
Dune::Functions::NedelecPreBasis::GridView
GV GridView
The grid view that the FE space is defined on.
Definition nedelecbasis.hh:159
Dune::Functions::NedelecPreBasis::finiteElementMap_
FiniteElementMap finiteElementMap_
Definition nedelecbasis.hh:204
Dune::Functions::NedelecPreBasis::Node
NedelecNode< GV, Range, kind, order > Node
Definition nedelecbasis.hh:162
Dune::Functions::NedelecPreBasis::update
void update(const GridView &gv)
Update the stored grid view, to be called if the grid has changed.
Definition nedelecbasis.hh:189
Dune::Functions::NedelecPreBasis::makeNode
Node makeNode() const
Create tree node.
Definition nedelecbasis.hh:198
Dune::Functions::BasisNodeMixin::setSize
void setSize(const size_type size)
Definition nodes.hh:197
Dune::Functions::LeafBasisNode
Definition nodes.hh:218