3 #ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
7 #include <dune/common/exceptions.hh>
9 #include <dune/grid/common/capabilities.hh>
10 #include <dune/grid/common/mcmgmapper.hh>
12 #include <dune/localfunctions/common/localfiniteelementvariant.hh>
13 #include <dune/localfunctions/nedelec.hh>
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;
39 constexpr
static GeometryType type = GeometryType(topologyId, GV::dimension);
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)
55 : elementMapper_(gv, mcmgElementLayout()),
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);
68 variants_.resize(numVariants(GeometryTypes::simplex(dim)) + numVariants(GeometryTypes::cube(dim)));
69 for (
size_t i = 0; i < numVariants(GeometryTypes::simplex(dim)); i++)
70 variants_[i] = SimplexFiniteElement(i);
71 for (
size_t i = 0; i < numVariants(GeometryTypes::cube(dim)); i++)
72 variants_[i + numVariants(GeometryTypes::simplex(dim))] = CubeFiniteElement(i);
77 const auto& indexSet = gv.indexSet();
79 for(
const auto& element : elements(gv))
81 const auto& refElement = referenceElement(element);
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())
101 orientation_[elementIndex] += numVariants(GeometryTypes::simplex(dim));
105 template<
class Element>
106 const auto& find(
const Element& element)
const
108 return variants_[orientation_[elementMapper_.index(element)]];
112 std::vector<FiniteElement> variants_;
113 const Dune::MultipleCodimMultipleGeomTypeMapper<GV> elementMapper_;
114 std::vector<unsigned short> orientation_;
132 template<
typename GV,
typename Range, std::
size_t kind,
int order>
135 template<
typename GV,
typename Range, std::
size_t kind,
int order,
class MI>
138 static const int dim = GV::dimension;
139 static_assert(kind==1,
"Only the Nedelec basis of the first kind is currently implemented!");
140 using FiniteElementMap =
typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
142 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
152 using IndexSet = Impl::DefaultNodeIndexSet<NedelecPreBasis>;
166 DUNE_THROW(NotImplemented,
"Only Nedelec elements of the first kind are implemented!");
170 if (gv.indexSet().types(0).size() > 2)
171 DUNE_THROW(NotImplemented,
"Nédélec basis is only implemented for grids with simplex and cube elements");
173 for(
auto type : gv.indexSet().types(0))
174 if (!type.isSimplex() && !type.isCube())
175 DUNE_THROW(NotImplemented,
"Nédélec basis is only implemented for grids with simplex or cube elements.");
178 DUNE_THROW(NotImplemented,
"Only first-order elements are implemented");
180 if (dim!=2 && dim!=3)
181 DUNE_THROW(NotImplemented,
"Only 2d and 3d Nédélec elements are implemented");
216 [[deprecated(
"Warning: The IndexSet typedef and the makeIndexSet method are deprecated. "\
217 "As a replacement use the indices() method of the PreBasis directly.")]]
231 assert(prefix.size() == 0 || prefix.size() == 1);
232 return (prefix.size() == 0) ?
size() : 0;
243 for (
auto&& type :
gridView_.indexSet().types(0))
245 size_type numEdges = referenceElement<typename GV::ctype,dim>(type).size(dim-1);
246 result = std::max(result, numEdges);
255 template<
typename It>
258 const auto& element = node.
element();
261 if (not(element.type().isCube()) and not(element.type().isSimplex()))
262 DUNE_THROW(NotImplemented,
"NedelecBasis only implemented for cube and simplex elements.");
264 for(std::size_t i=0, end=node.
size(); i<end; ++i, ++it)
266 Dune::LocalKey localKey = node.
finiteElement().localCoefficients().localKey(i);
267 *it = {
mapper_.subIndex(element, localKey.subEntity(), localKey.codim()) + localKey.index() };
281 template<
typename GV,
typename Range,
size_t kind,
int order>
285 static const int dim = GV::dimension;
290 using Element =
typename GV::template Codim<0>::Entity;
291 static_assert(kind==1,
"Only Nedelec elements of the first kind are implemented!");
292 using FiniteElementMap =
typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
293 using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::CovariantPiolaTransformator,
294 typename FiniteElementMap::FiniteElement,
334 namespace BasisFactory {
338 template<std::
size_t kind, std::
size_t order,
typename Range>
339 class NedelecPreBasisFactory
342 static const std::size_t requiredMultiIndexSize=1;
344 template<
class MultiIndex,
class Gr
idView>
345 auto makePreBasis(
const GridView& gridView)
const
363 template<std::
size_t kind, std::
size_t order,
typename Range=
double>
366 return Impl::NedelecPreBasisFactory<kind, order, Range>();
384 template<
typename GV, std::
size_t kind, std::
size_t order,
typename Range=
double>
auto power(ChildPreBasisFactory &&childPreBasisFactory, const IndexMergingStrategy &ims)
Create a pre-basis factory that can build a PowerPreBasis.
Definition: powerbasis.hh:425
auto nedelec()
Create a pre-basis factory that can create a Nédélec pre-basis.
Definition: nedelecbasis.hh:364
Definition: polynomial.hh:11
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:47
Definition: nedelecbasis.hh:284
const FiniteElementMap * finiteElementMap_
Definition: nedelecbasis.hh:329
FiniteElement finiteElement_
Definition: nedelecbasis.hh:327
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: nedelecbasis.hh:312
Impl::GlobalValuedLocalFiniteElement< Impl::CovariantPiolaTransformator, typename FiniteElementMap::FiniteElement, Element > FiniteElement
Definition: nedelecbasis.hh:295
void bind(const Element &e)
Bind to element.
Definition: nedelecbasis.hh:318
typename GV::template Codim< 0 >::Entity Element
Definition: nedelecbasis.hh:290
NedelecNode(const FiniteElementMap *finiteElementMap)
Definition: nedelecbasis.hh:297
typename Impl::Nedelec1stKindLocalFiniteElementMap< GV, dim, Range, order > FiniteElementMap
Definition: nedelecbasis.hh:292
const Element * element_
Definition: nedelecbasis.hh:328
const Element & element() const
Return current element, throw if unbound.
Definition: nedelecbasis.hh:303
std::size_t size_type
Definition: nedelecbasis.hh:289
Definition: nedelecbasis.hh:137
size_type dimension() const
Definition: nedelecbasis.hh:235
FiniteElementMap finiteElementMap_
Definition: nedelecbasis.hh:275
Mapper mapper_
Definition: nedelecbasis.hh:276
It indices(const Node &node, It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis.
Definition: nedelecbasis.hh:256
std::size_t size_type
Definition: nedelecbasis.hh:147
size_type size() const
Definition: nedelecbasis.hh:223
Impl::DefaultNodeIndexSet< NedelecPreBasis > IndexSet
Type of created tree node index set.
Definition: nedelecbasis.hh:152
void update(const GridView &gv)
Definition: nedelecbasis.hh:195
GridView gridView_
Definition: nedelecbasis.hh:274
Dune::ReservedVector< size_type, 1 > SizePrefix
Definition: nedelecbasis.hh:157
NedelecPreBasis(const GridView &gv)
Constructor for a given grid view object.
Definition: nedelecbasis.hh:160
size_type maxNodeSize() const
Definition: nedelecbasis.hh:240
void initializeIndices()
Definition: nedelecbasis.hh:184
IndexSet makeIndexSet() const
Create tree node index set.
Definition: nedelecbasis.hh:218
GV GridView
The grid view that the FE space is defined on.
Definition: nedelecbasis.hh:146
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: nedelecbasis.hh:155
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: nedelecbasis.hh:189
Node makeNode() const
Create tree node.
Definition: nedelecbasis.hh:204
size_type size(const SizePrefix prefix) const
Return number possible values for next position in multi index.
Definition: nedelecbasis.hh:229
size_type size() const
Definition: nodes.hh:140
void setSize(const size_type size)
Definition: nodes.hh:162