DUNE PDELab (git)

nedelecbasis.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
9
10#include <array>
12
15
16#include <dune/localfunctions/common/localfiniteelementvariant.hh>
17#include <dune/localfunctions/nedelec.hh>
18
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>
23
24namespace Dune::Functions
25{
26
27namespace Impl
28{
29 template<typename GV, int dim, typename R, std::size_t order>
30 class Nedelec1stKindLocalFiniteElementMap
31 {
32 using D = typename GV::ctype;
33 constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
34
35 using CubeFiniteElement = Nedelec1stKindCubeLocalFiniteElement<D,R,dim,order>;
36 using SimplexFiniteElement = Nedelec1stKindSimplexLocalFiniteElement<D,R,dim,order>;
37
38 public:
39
40 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
41
42 constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId; // meaningless if hasFixedElementType is false
43 constexpr static GeometryType type = GeometryType(topologyId, GV::dimension);
44
45 using FiniteElement = std::conditional_t<hasFixedElementType,
46 std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
47 LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
48
49 static std::size_t numVariants(GeometryType type)
50 {
51 if (order!=1) // I am not sure whether the formula below is correct for all orders.
52 DUNE_THROW(NotImplemented, "Only Nedelec elements of order 1 are implemented!");
53
54 auto numEdges = referenceElement<D,dim>(type).size(dim-1);
55 return power(2,numEdges);
56 }
57
58 Nedelec1stKindLocalFiniteElementMap(const GV& gv)
59 : elementMapper_(gv, mcmgElementLayout()),
60 orientation_(gv.size(0))
61 {
62 // create all variants
63 if constexpr (hasFixedElementType)
64 {
65 variants_.resize(numVariants(type));
66 for (size_t i = 0; i < numVariants(type); i++)
67 variants_[i] = FiniteElement(i);
68 }
69 else
70 {
71 // for mixed grids add offset for cubes
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);
77 }
78
79
80 // compute orientation for all elements
81 const auto& indexSet = gv.indexSet();
82
83 for(const auto& element : elements(gv))
84 {
85 const auto& refElement = referenceElement(element);
86 auto elementIndex = elementMapper_.index(element);
87 orientation_[elementIndex] = 0;
88
89 for (std::size_t i=0; i<element.subEntities(dim-1); i++)
90 {
91 // Local vertex indices within the element
92 auto localV0 = refElement.subEntity(i,dim-1, 0,dim);
93 auto localV1 = refElement.subEntity(i,dim-1, 1,dim);
94
95 // Global vertex indices within the grid
96 auto globalV0 = indexSet.subIndex(element,localV0,dim);
97 auto globalV1 = indexSet.subIndex(element,localV1,dim);
98
99 if ( (localV0<localV1 && globalV0>globalV1) || (localV0>localV1 && globalV0<globalV1) )
100 orientation_[elementIndex] |= (1 << i);
101 }
102 // for mixed grids add offset for cubes
103 if constexpr (!hasFixedElementType)
104 if (element.type().isCube())
105 orientation_[elementIndex] += numVariants(GeometryTypes::simplex(dim));
106 }
107 }
108
109 template<class Element>
110 const auto& find(const Element& element) const
111 {
112 return variants_[orientation_[elementMapper_.index(element)]];
113 }
114
115 private:
116 std::vector<FiniteElement> variants_;
118 std::vector<unsigned short> orientation_;
119 };
120
121
122} // namespace Impl
123
124
125// *****************************************************************************
126// This is the reusable part of the basis. It contains
127//
128// NedelecPreBasis
129// NedelecNode
130//
131// The pre-basis allows to create the others and is the owner of possible shared
132// state. These components do _not_ depend on the global basis and local view
133// and can be used without a global basis.
134// *****************************************************************************
135
136template<typename GV, typename Range, std::size_t kind, int order>
137class NedelecNode;
138
139template<typename GV, typename Range, std::size_t kind, int order>
140class NedelecPreBasis :
141 public LeafPreBasisMixin< NedelecPreBasis<GV,Range,kind,order> >
142{
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>;
146
148public:
149
151 using GridView = GV;
152 using size_type = std::size_t;
153
154 using Node = NedelecNode<GV, Range, kind, order>;
155
157 NedelecPreBasis(const GridView& gv) :
158 gridView_(gv),
159 finiteElementMap_(gv),
160 mapper_(gridView_, mcmgLayout(Dim<1>{}))
161 {
162 if (kind!=1)
163 DUNE_THROW(NotImplemented, "Only Nedelec elements of the first kind are implemented!");
164
165 // There is no inherent reason why the basis shouldn't work for grids with more than two
166 // element types. Somebody simply has to sit down and implement the missing bits.
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");
169
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.");
173
174 if (order>1)
175 DUNE_THROW(NotImplemented, "Only first-order elements are implemented");
176
177 if (dim!=2 && dim!=3)
178 DUNE_THROW(NotImplemented, "Only 2d and 3d Nédélec elements are implemented");
179 }
180
181 void initializeIndices()
182 {}
183
186 const GridView& gridView() const
187 {
188 return gridView_;
189 }
190
191 /* \brief Update the stored grid view, to be called if the grid has changed */
192 void update (const GridView& gv)
193 {
194 gridView_ = gv;
195 mapper_.update(gridView_);
196 }
197
201 Node makeNode() const
202 {
203 return Node{&finiteElementMap_};
204 }
205
206 size_type dimension() const
207 {
208 return mapper_.size();
209 }
210
211 size_type maxNodeSize() const
212 {
213 size_type result = 0;
214 for (auto&& type : gridView_.indexSet().types(0))
215 {
216 size_type numEdges = referenceElement<typename GV::ctype,dim>(type).size(dim-1);
217 result = std::max(result, numEdges);
218 }
219
220 return result;
221 }
222
226 template<typename It>
227 It indices(const Node& node, It it) const
228 {
229 const auto& element = node.element();
230
231 // throw if Element is not of predefined type
232 if (not(element.type().isCube()) and not(element.type().isSimplex()))
233 DUNE_THROW(NotImplemented, "NedelecBasis only implemented for cube and simplex elements.");
234
235 for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
236 {
237 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
238 *it = { mapper_.subIndex(element, localKey.subEntity(), localKey.codim()) + localKey.index() };
239 }
240
241 return it;
242 }
243
244protected:
245 GridView gridView_;
246 FiniteElementMap finiteElementMap_;
247 Mapper mapper_;
248};
249
250
251
252template<typename GV, typename Range, size_t kind, int order>
253class NedelecNode :
254 public LeafBasisNode
255{
256 static const int dim = GV::dimension;
257
258public:
259
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,
266 Element>;
267
268 NedelecNode(const FiniteElementMap* finiteElementMap) :
269 element_(nullptr),
270 finiteElementMap_(finiteElementMap)
271 { }
272
274 const Element& element() const
275 {
276 return *element_;
277 }
278
283 const FiniteElement& finiteElement() const
284 {
285 return finiteElement_;
286 }
287
289 void bind(const Element& e)
290 {
291 element_ = &e;
292 finiteElement_.bind((finiteElementMap_->find(*element_)), e);
293 this->setSize(finiteElement_.size());
294 }
295
296protected:
297
298 FiniteElement finiteElement_;
299 const Element* element_;
300 const FiniteElementMap* finiteElementMap_;
301};
302
303
304
305namespace BasisFactory {
306
316template<std::size_t kind, std::size_t order, typename Range=double>
318{
319 return [](const auto& gridView) {
320 return NedelecPreBasis<std::decay_t<decltype(gridView)>, Range, kind, order>(gridView);
321 };
322}
323
324} // end namespace BasisFactory
325
326
327
328// *****************************************************************************
329// This is the actual global basis implementation based on the reusable parts.
330// *****************************************************************************
331
339template<typename GV, std::size_t kind, std::size_t order, typename Range=double>
340using NedelecBasis = DefaultGlobalBasis<NedelecPreBasis<GV, Range, kind, order > >;
341
342} // end namespace Dune::Functions
343
344
345#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
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:317
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)