7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
11#include <dune/common/exceptions.hh>
13#include <dune/grid/common/capabilities.hh>
14#include <dune/grid/common/mcmgmapper.hh>
16#include <dune/localfunctions/common/localfiniteelementvariant.hh>
17#include <dune/localfunctions/nedelec.hh>
19#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
20#include <dune/functions/functionspacebases/globalvaluedlocalfiniteelement.hh>
21#include <dune/functions/functionspacebases/leafprebasismixin.hh>
22#include <dune/functions/functionspacebases/nodes.hh>
24namespace Dune::Functions
29 template<
typename GV,
int dim,
typename R, std::
size_t order>
30 class Nedelec1stKindLocalFiniteElementMap
32 using D =
typename GV::ctype;
33 constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
35 using CubeFiniteElement = Nedelec1stKindCubeLocalFiniteElement<D,R,dim,order>;
36 using SimplexFiniteElement = Nedelec1stKindSimplexLocalFiniteElement<D,R,dim,order>;
40 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
42 constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId;
43 constexpr static GeometryType type = GeometryType(topologyId, GV::dimension);
45 using FiniteElement = std::conditional_t<hasFixedElementType,
46 std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
47 LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
49 static std::size_t numVariants(GeometryType type)
52 DUNE_THROW(NotImplemented,
"Only Nedelec elements of order 1 are implemented!");
54 auto numEdges = referenceElement<D,dim>(type).size(dim-1);
55 return power(2,numEdges);
58 Nedelec1stKindLocalFiniteElementMap(
const GV& gv)
59 : elementMapper_(gv, mcmgElementLayout()),
60 orientation_(gv.size(0))
63 if constexpr (hasFixedElementType)
65 variants_.resize(numVariants(type));
66 for (
size_t i = 0; i < numVariants(type); i++)
67 variants_[i] = FiniteElement(i);
72 variants_.resize(numVariants(GeometryTypes::simplex(dim)) + numVariants(GeometryTypes::cube(dim)));
73 for (
size_t i = 0; i < numVariants(GeometryTypes::simplex(dim)); i++)
74 variants_[i] = SimplexFiniteElement(i);
75 for (
size_t i = 0; i < numVariants(GeometryTypes::cube(dim)); i++)
76 variants_[i + numVariants(GeometryTypes::simplex(dim))] = CubeFiniteElement(i);
81 const auto& indexSet = gv.indexSet();
83 for(
const auto& element : elements(gv))
85 const auto& refElement = referenceElement(element);
86 auto elementIndex = elementMapper_.index(element);
87 orientation_[elementIndex] = 0;
89 for (std::size_t i=0; i<element.subEntities(dim-1); i++)
92 auto localV0 = refElement.subEntity(i,dim-1, 0,dim);
93 auto localV1 = refElement.subEntity(i,dim-1, 1,dim);
96 auto globalV0 = indexSet.subIndex(element,localV0,dim);
97 auto globalV1 = indexSet.subIndex(element,localV1,dim);
99 if ( (localV0<localV1 && globalV0>globalV1) || (localV0>localV1 && globalV0<globalV1) )
100 orientation_[elementIndex] |= (1 << i);
103 if constexpr (!hasFixedElementType)
104 if (element.type().isCube())
105 orientation_[elementIndex] += numVariants(GeometryTypes::simplex(dim));
109 template<
class Element>
110 const auto& find(
const Element& element)
const
112 return variants_[orientation_[elementMapper_.index(element)]];
116 std::vector<FiniteElement> variants_;
117 Dune::MultipleCodimMultipleGeomTypeMapper<GV> elementMapper_;
118 std::vector<unsigned short> orientation_;
136template<
typename GV,
typename Range, std::
size_t kind,
int order>
139template<
typename GV,
typename Range, std::
size_t kind,
int order>
140class NedelecPreBasis :
141 public LeafPreBasisMixin< NedelecPreBasis<GV,Range,kind,order> >
143 static const int dim = GV::dimension;
144 static_assert(kind==1,
"Only the Nedelec basis of the first kind is currently implemented!");
145 using FiniteElementMap =
typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
147 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
152 using size_type = std::size_t;
154 using Node = NedelecNode<GV, Range, kind, order>;
157 NedelecPreBasis(
const GridView& gv) :
159 finiteElementMap_(gv),
160 mapper_(gridView_, mcmgLayout(Dim<1>{}))
163 DUNE_THROW(NotImplemented,
"Only Nedelec elements of the first kind are implemented!");
167 if (gv.indexSet().types(0).size() > 2)
168 DUNE_THROW(NotImplemented,
"Nédélec basis is only implemented for grids with simplex and cube elements");
170 for(
auto type : gv.indexSet().types(0))
171 if (!type.isSimplex() && !type.isCube())
172 DUNE_THROW(NotImplemented,
"Nédélec basis is only implemented for grids with simplex or cube elements.");
175 DUNE_THROW(NotImplemented,
"Only first-order elements are implemented");
177 if (dim!=2 && dim!=3)
178 DUNE_THROW(NotImplemented,
"Only 2d and 3d Nédélec elements are implemented");
181 void initializeIndices()
186 const GridView& gridView()
const
192 void update (
const GridView& gv)
195 mapper_.update(gridView_);
201 Node makeNode()
const
203 return Node{&finiteElementMap_};
206 size_type dimension()
const
208 return mapper_.size();
211 size_type maxNodeSize()
const
213 size_type result = 0;
214 for (
auto&& type : gridView_.indexSet().types(0))
216 size_type numEdges = referenceElement<typename GV::ctype,dim>(type).size(dim-1);
217 result = std::max(result, numEdges);
226 template<
typename It>
227 It indices(
const Node& node, It it)
const
229 const auto& element = node.element();
232 if (not(element.type().isCube()) and not(element.type().isSimplex()))
233 DUNE_THROW(NotImplemented,
"NedelecBasis only implemented for cube and simplex elements.");
235 for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
237 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
238 *it = { mapper_.subIndex(element, localKey.subEntity(), localKey.codim()) + localKey.index() };
246 FiniteElementMap finiteElementMap_;
252template<
typename GV,
typename Range,
size_t kind,
int order>
256 static const int dim = GV::dimension;
260 using size_type = std::size_t;
261 using Element =
typename GV::template Codim<0>::Entity;
262 static_assert(kind==1,
"Only Nedelec elements of the first kind are implemented!");
263 using FiniteElementMap =
typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
264 using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::CovariantPiolaTransformator,
265 typename FiniteElementMap::FiniteElement,
268 NedelecNode(
const FiniteElementMap* finiteElementMap) :
270 finiteElementMap_(finiteElementMap)
274 const Element& element()
const
283 const FiniteElement& finiteElement()
const
285 return finiteElement_;
289 void bind(
const Element& e)
292 finiteElement_.bind((finiteElementMap_->find(*element_)), e);
293 this->setSize(finiteElement_.size());
298 FiniteElement finiteElement_;
299 const Element* element_;
300 const FiniteElementMap* finiteElementMap_;
305namespace BasisFactory {
316template<std::
size_t kind, std::
size_t order,
typename Range=
double>
319 return [](
const auto& gridView) {
320 return NedelecPreBasis<std::decay_t<
decltype(gridView)>, Range, kind, order>(gridView);
339template<
typename GV, std::
size_t kind, std::
size_t order,
typename Range=
double>
340using NedelecBasis = DefaultGlobalBasis<NedelecPreBasis<GV, Range, kind, order > >;
auto power(ChildPreBasisFactory &&childPreBasisFactory, std::size_t k, const IndexMergingStrategy &)
Create a pre-basis factory that can build a PowerPreBasis.
Definition: dynamicpowerbasis.hh:409
auto nedelec()
Create a pre-basis factory that can create a Nédélec pre-basis.
Definition: nedelecbasis.hh:317