dune-functions  2.9.0
nedelecbasis.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
5 
6 #include <array>
7 #include <dune/common/exceptions.hh>
8 
9 #include <dune/grid/common/capabilities.hh>
10 #include <dune/grid/common/mcmgmapper.hh>
11 
12 #include <dune/localfunctions/common/localfiniteelementvariant.hh>
13 #include <dune/localfunctions/nedelec.hh>
14 
18 
19 namespace Dune::Functions
20 {
21 
22 namespace Impl
23 {
24  template<typename GV, int dim, typename R, std::size_t order>
25  class Nedelec1stKindLocalFiniteElementMap
26  {
27  using D = typename GV::ctype;
28  constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
29 
30  using CubeFiniteElement = Nedelec1stKindCubeLocalFiniteElement<D,R,dim,order>;
31  using SimplexFiniteElement = Nedelec1stKindSimplexLocalFiniteElement<D,R,dim,order>;
32 
33  public:
34 
35  using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
36 
37  constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId; // meaningless if hasFixedElementType is false
38  constexpr static GeometryType type = GeometryType(topologyId, GV::dimension);
39 
40  using FiniteElement = std::conditional_t<hasFixedElementType,
41  std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
42  LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
43 
44  static std::size_t numVariants(GeometryType type)
45  {
46  if (order!=1) // I am not sure whether the formula below is correct for all orders.
47  DUNE_THROW(NotImplemented, "Only Nedelec elements of order 1 are implemented!");
48 
49  auto numEdges = referenceElement<D,dim>(type).size(dim-1);
50  return power(2,numEdges);
51  }
52 
53  Nedelec1stKindLocalFiniteElementMap(const GV& gv)
54  : elementMapper_(gv, mcmgElementLayout()),
55  orientation_(gv.size(0))
56  {
57  // create all variants
58  if constexpr (hasFixedElementType)
59  {
60  variants_.resize(numVariants(type));
61  for (size_t i = 0; i < numVariants(type); i++)
62  variants_[i] = FiniteElement(i);
63  }
64  else
65  {
66  // for mixed grids add offset for cubes
67  variants_.resize(numVariants(GeometryTypes::simplex(dim)) + numVariants(GeometryTypes::cube(dim)));
68  for (size_t i = 0; i < numVariants(GeometryTypes::simplex(dim)); i++)
69  variants_[i] = SimplexFiniteElement(i);
70  for (size_t i = 0; i < numVariants(GeometryTypes::cube(dim)); i++)
71  variants_[i + numVariants(GeometryTypes::simplex(dim))] = CubeFiniteElement(i);
72  }
73 
74 
75  // compute orientation for all elements
76  const auto& indexSet = gv.indexSet();
77 
78  for(const auto& element : elements(gv))
79  {
80  const auto& refElement = referenceElement(element);
81  auto elementIndex = elementMapper_.index(element);
82  orientation_[elementIndex] = 0;
83 
84  for (std::size_t i=0; i<element.subEntities(dim-1); i++)
85  {
86  // Local vertex indices within the element
87  auto localV0 = refElement.subEntity(i,dim-1, 0,dim);
88  auto localV1 = refElement.subEntity(i,dim-1, 1,dim);
89 
90  // Global vertex indices within the grid
91  auto globalV0 = indexSet.subIndex(element,localV0,dim);
92  auto globalV1 = indexSet.subIndex(element,localV1,dim);
93 
94  if ( (localV0<localV1 && globalV0>globalV1) || (localV0>localV1 && globalV0<globalV1) )
95  orientation_[elementIndex] |= (1 << i);
96  }
97  // for mixed grids add offset for cubes
98  if constexpr (!hasFixedElementType)
99  if (element.type().isCube())
100  orientation_[elementIndex] += numVariants(GeometryTypes::simplex(dim));
101  }
102  }
103 
104  template<class Element>
105  const auto& find(const Element& element) const
106  {
107  return variants_[orientation_[elementMapper_.index(element)]];
108  }
109 
110  private:
111  std::vector<FiniteElement> variants_;
112  const Dune::MultipleCodimMultipleGeomTypeMapper<GV> elementMapper_;
113  std::vector<unsigned short> orientation_;
114  };
115 
116 
117 } // namespace Impl
118 
119 
120 // *****************************************************************************
121 // This is the reusable part of the basis. It contains
122 //
123 // NedelecPreBasis
124 // NedelecNode
125 //
126 // The pre-basis allows to create the others and is the owner of possible shared
127 // state. These components do _not_ depend on the global basis and local view
128 // and can be used without a global basis.
129 // *****************************************************************************
130 
131 template<typename GV, typename Range, std::size_t kind, int order>
132 class NedelecNode;
133 
134 template<typename GV, typename Range, std::size_t kind, int order>
136 {
137  static const int dim = GV::dimension;
138  static_assert(kind==1, "Only the Nedelec basis of the first kind is currently implemented!");
139  using FiniteElementMap = typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
140 
141  using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
142 public:
143 
145  using GridView = GV;
146  using size_type = std::size_t;
147 
149 
150  static constexpr size_type maxMultiIndexSize = 1;
151  static constexpr size_type minMultiIndexSize = 1;
152  static constexpr size_type multiIndexBufferSize = 1;
153 
156  gridView_(gv),
157  finiteElementMap_(gv),
158  mapper_(gridView_, mcmgLayout(Dim<1>{}))
159  {
160  if (kind!=1)
161  DUNE_THROW(NotImplemented, "Only Nedelec elements of the first kind are implemented!");
162 
163  // There is no inherent reason why the basis shouldn't work for grids with more than two
164  // element types. Somebody simply has to sit down and implement the missing bits.
165  if (gv.indexSet().types(0).size() > 2)
166  DUNE_THROW(NotImplemented, "Nédélec basis is only implemented for grids with simplex and cube elements");
167 
168  for(auto type : gv.indexSet().types(0))
169  if (!type.isSimplex() && !type.isCube())
170  DUNE_THROW(NotImplemented, "Nédélec basis is only implemented for grids with simplex or cube elements.");
171 
172  if (order>1)
173  DUNE_THROW(NotImplemented, "Only first-order elements are implemented");
174 
175  if (dim!=2 && dim!=3)
176  DUNE_THROW(NotImplemented, "Only 2d and 3d Nédélec elements are implemented");
177  }
178 
180  {}
181 
184  const GridView& gridView() const
185  {
186  return gridView_;
187  }
188 
189  /* \brief Update the stored grid view, to be called if the grid has changed */
190  void update (const GridView& gv)
191  {
192  gridView_ = gv;
193  mapper_.update(gridView_);
194  }
195 
199  Node makeNode() const
200  {
201  return Node{&finiteElementMap_};
202  }
203 
204  size_type size() const
205  {
206  return mapper_.size();
207  }
208 
210  template<class SizePrefix>
211  size_type size(const SizePrefix& prefix) const
212  {
213  assert(prefix.size() == 0 || prefix.size() == 1);
214  return (prefix.size() == 0) ? size() : 0;
215  }
216 
218  {
219  return size();
220  }
221 
223  {
224  size_type result = 0;
225  for (auto&& type : gridView_.indexSet().types(0))
226  {
227  size_type numEdges = referenceElement<typename GV::ctype,dim>(type).size(dim-1);
228  result = std::max(result, numEdges);
229  }
230 
231  return result;
232  }
233 
237  template<typename It>
238  It indices(const Node& node, It it) const
239  {
240  const auto& element = node.element();
241 
242  // throw if Element is not of predefined type
243  if (not(element.type().isCube()) and not(element.type().isSimplex()))
244  DUNE_THROW(NotImplemented, "NedelecBasis only implemented for cube and simplex elements.");
245 
246  for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
247  {
248  Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
249  *it = { mapper_.subIndex(element, localKey.subEntity(), localKey.codim()) + localKey.index() };
250  }
251 
252  return it;
253  }
254 
255 protected:
257  FiniteElementMap finiteElementMap_;
258  Mapper mapper_;
259 };
260 
261 
262 
263 template<typename GV, typename Range, size_t kind, int order>
264 class NedelecNode :
265  public LeafBasisNode
266 {
267  static const int dim = GV::dimension;
268 
269 public:
270 
271  using size_type = std::size_t;
272  using Element = typename GV::template Codim<0>::Entity;
273  static_assert(kind==1, "Only Nedelec elements of the first kind are implemented!");
274  using FiniteElementMap = typename Impl::Nedelec1stKindLocalFiniteElementMap<GV, dim, Range, order>;
275  using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::CovariantPiolaTransformator,
276  typename FiniteElementMap::FiniteElement,
277  Element>;
278 
279  NedelecNode(const FiniteElementMap* finiteElementMap) :
280  element_(nullptr),
281  finiteElementMap_(finiteElementMap)
282  { }
283 
285  const Element& element() const
286  {
287  return *element_;
288  }
289 
295  {
296  return finiteElement_;
297  }
298 
300  void bind(const Element& e)
301  {
302  element_ = &e;
303  finiteElement_.bind((finiteElementMap_->find(*element_)), e);
304  this->setSize(finiteElement_.size());
305  }
306 
307 protected:
308 
312 };
313 
314 
315 
316 namespace BasisFactory {
317 
327 template<std::size_t kind, std::size_t order, typename Range=double>
328 auto nedelec()
329 {
330  return [](const auto& gridView) {
331  return NedelecPreBasis<std::decay_t<decltype(gridView)>, Range, kind, order>(gridView);
332  };
333 }
334 
335 } // end namespace BasisFactory
336 
337 
338 
339 // *****************************************************************************
340 // This is the actual global basis implementation based on the reusable parts.
341 // *****************************************************************************
342 
350 template<typename GV, std::size_t kind, std::size_t order, typename Range=double>
352 
353 } // end namespace Dune::Functions
354 
355 
356 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
auto power(ChildPreBasisFactory &&childPreBasisFactory, const IndexMergingStrategy &)
Create a pre-basis factory that can build a PowerPreBasis.
Definition: powerbasis.hh:369
auto nedelec()
Create a pre-basis factory that can create a Nédélec pre-basis.
Definition: nedelecbasis.hh:328
Definition: polynomial.hh:11
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:46
Definition: nedelecbasis.hh:266
const FiniteElementMap * finiteElementMap_
Definition: nedelecbasis.hh:311
FiniteElement finiteElement_
Definition: nedelecbasis.hh:309
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: nedelecbasis.hh:294
Impl::GlobalValuedLocalFiniteElement< Impl::CovariantPiolaTransformator, typename FiniteElementMap::FiniteElement, Element > FiniteElement
Definition: nedelecbasis.hh:277
void bind(const Element &e)
Bind to element.
Definition: nedelecbasis.hh:300
typename GV::template Codim< 0 >::Entity Element
Definition: nedelecbasis.hh:272
NedelecNode(const FiniteElementMap *finiteElementMap)
Definition: nedelecbasis.hh:279
typename Impl::Nedelec1stKindLocalFiniteElementMap< GV, dim, Range, order > FiniteElementMap
Definition: nedelecbasis.hh:274
const Element * element_
Definition: nedelecbasis.hh:310
const Element & element() const
Return current element, throw if unbound.
Definition: nedelecbasis.hh:285
std::size_t size_type
Definition: nedelecbasis.hh:271
Definition: nedelecbasis.hh:136
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:238
size_type size(const SizePrefix &prefix) const
Return number possible values for next position in multi index.
Definition: nedelecbasis.hh:211
std::size_t size_type
Definition: nedelecbasis.hh:146
NedelecPreBasis(const GridView &gv)
Constructor for a given grid view object.
Definition: nedelecbasis.hh:155
GV GridView
The grid view that the FE space is defined on.
Definition: nedelecbasis.hh:145
size_type dimension() const
Definition: nedelecbasis.hh:217
static constexpr size_type maxMultiIndexSize
Definition: nedelecbasis.hh:150
static constexpr size_type multiIndexBufferSize
Definition: nedelecbasis.hh:152
GridView gridView_
Definition: nedelecbasis.hh:256
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: nedelecbasis.hh:184
FiniteElementMap finiteElementMap_
Definition: nedelecbasis.hh:257
static constexpr size_type minMultiIndexSize
Definition: nedelecbasis.hh:151
void initializeIndices()
Definition: nedelecbasis.hh:179
void update(const GridView &gv)
Definition: nedelecbasis.hh:190
size_type maxNodeSize() const
Definition: nedelecbasis.hh:222
Node makeNode() const
Create tree node.
Definition: nedelecbasis.hh:199
size_type size() const
Definition: nedelecbasis.hh:204
Mapper mapper_
Definition: nedelecbasis.hh:258
size_type size() const
Definition: nodes.hh:142
void setSize(const size_type size)
Definition: nodes.hh:164
Definition: nodes.hh:186