3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
12#include <dune/localfunctions/common/localfiniteelementvariant.hh>
13#include <dune/localfunctions/nedelec.hh>
15#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
16#include <dune/functions/functionspacebases/globalvaluedlocalfiniteelement.hh>
17#include <dune/functions/functionspacebases/leafprebasismixin.hh>
18#include <dune/functions/functionspacebases/nodes.hh>
20namespace Dune::Functions
25 template<
typename GV,
int dim,
typename R, std::
size_t order>
26 class Nedelec1stKindLocalFiniteElementMap
28 using D =
typename GV::ctype;
29 constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
31 using CubeFiniteElement = Nedelec1stKindCubeLocalFiniteElement<D,R,dim,order>;
32 using SimplexFiniteElement = Nedelec1stKindSimplexLocalFiniteElement<D,R,dim,order>;
36 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
38 constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId;
41 using FiniteElement = std::conditional_t<hasFixedElementType,
42 std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
43 LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
45 static std::size_t numVariants(GeometryType type)
48 DUNE_THROW(NotImplemented,
"Only Nedelec elements of order 1 are implemented!");
50 auto numEdges = referenceElement<D,dim>(type).size(dim-1);
51 return power(2,numEdges);
54 Nedelec1stKindLocalFiniteElementMap(
const GV& gv)
56 orientation_(gv.
size(0))
59 if constexpr (hasFixedElementType)
61 variants_.resize(numVariants(type));
62 for (
size_t i = 0; i < numVariants(type); i++)
63 variants_[i] = FiniteElement(i);
70 variants_[i] = SimplexFiniteElement(i);
77 const auto& indexSet = gv.indexSet();
79 for(
const auto& element : elements(gv))
82 auto elementIndex = elementMapper_.index(element);
83 orientation_[elementIndex] = 0;
85 for (std::size_t i=0; i<element.subEntities(dim-1); i++)
88 auto localV0 = refElement.subEntity(i,dim-1, 0,dim);
89 auto localV1 = refElement.subEntity(i,dim-1, 1,dim);
92 auto globalV0 = indexSet.subIndex(element,localV0,dim);
93 auto globalV1 = indexSet.subIndex(element,localV1,dim);
95 if ( (localV0<localV1 && globalV0>globalV1) || (localV0>localV1 && globalV0<globalV1) )
96 orientation_[elementIndex] |= (1 << i);
99 if constexpr (!hasFixedElementType)
100 if (element.type().isCube())
105 template<
class Element>
106 const auto& find(
const Element& element)
const
108 return variants_[orientation_[elementMapper_.index(element)]];
112 std::vector<FiniteElement> variants_;
114 std::vector<unsigned short> orientation_;
132template<
typename GV,
typename Range, std::
size_t kind,
int order>
135template<
typename GV,
typename Range, std::
size_t kind,
int order>
136class NedelecPreBasis :
137 public LeafPreBasisMixin< NedelecPreBasis<GV,Range,kind,order> >
139 static const int dim = GV::dimension;
140 static_assert(kind==1,
"Only the Nedelec basis of the first kind is currently implemented!");
141 using FiniteElementMap =
typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
148 using size_type = std::size_t;
150 using Node = NedelecNode<GV, Range, kind, order>;
153 NedelecPreBasis(
const GridView& gv) :
155 finiteElementMap_(gv),
159 DUNE_THROW(NotImplemented,
"Only Nedelec elements of the first kind are implemented!");
163 if (gv.indexSet().types(0).size() > 2)
164 DUNE_THROW(NotImplemented,
"Nédélec basis is only implemented for grids with simplex and cube elements");
166 for(
auto type : gv.indexSet().types(0))
167 if (!type.isSimplex() && !type.isCube())
168 DUNE_THROW(NotImplemented,
"Nédélec basis is only implemented for grids with simplex or cube elements.");
171 DUNE_THROW(NotImplemented,
"Only first-order elements are implemented");
173 if (dim!=2 && dim!=3)
174 DUNE_THROW(NotImplemented,
"Only 2d and 3d Nédélec elements are implemented");
177 void initializeIndices()
182 const GridView& gridView()
const
188 void update (
const GridView& gv)
191 mapper_.update(gridView_);
197 Node makeNode()
const
199 return Node{&finiteElementMap_};
202 size_type dimension()
const
204 return mapper_.size();
207 size_type maxNodeSize()
const
209 size_type result = 0;
210 for (
auto&& type : gridView_.indexSet().types(0))
212 size_type numEdges = referenceElement<typename GV::ctype,dim>(type).size(dim-1);
213 result =
std::max(result, numEdges);
222 template<
typename It>
223 It indices(
const Node& node, It it)
const
225 const auto& element = node.element();
228 if (not(element.type().isCube()) and not(element.type().isSimplex()))
229 DUNE_THROW(NotImplemented,
"NedelecBasis only implemented for cube and simplex elements.");
231 for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
233 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
234 *it = { mapper_.subIndex(element, localKey.
subEntity(), localKey.
codim()) + localKey.
index() };
242 FiniteElementMap finiteElementMap_;
248template<
typename GV,
typename Range,
size_t kind,
int order>
252 static const int dim = GV::dimension;
256 using size_type = std::size_t;
257 using Element =
typename GV::template Codim<0>::Entity;
258 static_assert(kind==1,
"Only Nedelec elements of the first kind are implemented!");
259 using FiniteElementMap =
typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
260 using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::CovariantPiolaTransformator,
261 typename FiniteElementMap::FiniteElement,
264 NedelecNode(
const FiniteElementMap* finiteElementMap) :
266 finiteElementMap_(finiteElementMap)
270 const Element& element()
const
279 const FiniteElement& finiteElement()
const
281 return finiteElement_;
285 void bind(
const Element& e)
288 finiteElement_.bind((finiteElementMap_->find(*element_)), e);
289 this->setSize(finiteElement_.size());
294 FiniteElement finiteElement_;
295 const Element* element_;
296 const FiniteElementMap* finiteElementMap_;
301namespace BasisFactory {
312template<std::
size_t kind, std::
size_t order,
typename Range=
double>
315 return [](
const auto& gridView) {
316 return NedelecPreBasis<std::decay_t<
decltype(gridView)>, Range, kind, order>(gridView);
335template<
typename GV, std::
size_t kind, std::
size_t order,
typename Range=
double>
336using NedelecBasis = DefaultGlobalBasis<NedelecPreBasis<GV, Range, kind, order > >;
Describe position of one degree of freedom.
Definition: localkey.hh:24
constexpr unsigned int index() const noexcept
Return offset within subentity.
Definition: localkey.hh:70
constexpr unsigned int codim() const noexcept
Return codim of associated entity.
Definition: localkey.hh:63
constexpr unsigned int subEntity() const noexcept
Return number of associated subentity.
Definition: localkey.hh:56
Implementation class for a multiple codim and multiple geometry type mapper.
Definition: mcmgmapper.hh:129
A set of traits classes to store static information about grid implementation.
A few common exception classes.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
auto nedelec()
Create a pre-basis factory that can create a Nédélec pre-basis.
Definition: nedelecbasis.hh:313
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
MCMGLayout mcmgLayout(Codim< codim >)
layout for entities of codimension codim
Definition: mcmgmapper.hh:72
MCMGLayout mcmgElementLayout()
layout for elements (codim-0 entities)
Definition: mcmgmapper.hh:97
Mapper for multiple codim and multiple geometry types.
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75