dune-functions  2.8.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 
19 
20 namespace Dune::Functions
21 {
22 
23 namespace Impl
24 {
25  template<typename GV, int dim, typename R, std::size_t order>
26  class Nedelec1stKindLocalFiniteElementMap
27  {
28  using D = typename GV::ctype;
29  constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
30 
31  using CubeFiniteElement = Nedelec1stKindCubeLocalFiniteElement<D,R,dim,order>;
32  using SimplexFiniteElement = Nedelec1stKindSimplexLocalFiniteElement<D,R,dim,order>;
33 
34  public:
35 
36  using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
37 
38  constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId; // meaningless if hasFixedElementType is false
39  constexpr static GeometryType type = GeometryType(topologyId, GV::dimension);
40 
41  using FiniteElement = std::conditional_t<hasFixedElementType,
42  std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
43  LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
44 
45  static std::size_t numVariants(GeometryType type)
46  {
47  if (order!=1) // I am not sure whether the formula below is correct for all orders.
48  DUNE_THROW(NotImplemented, "Only Nedelec elements of order 1 are implemented!");
49 
50  auto numEdges = referenceElement<D,dim>(type).size(dim-1);
51  return power(2,numEdges);
52  }
53 
54  Nedelec1stKindLocalFiniteElementMap(const GV& gv)
55  : elementMapper_(gv, mcmgElementLayout()),
56  orientation_(gv.size(0))
57  {
58  // create all variants
59  if constexpr (hasFixedElementType)
60  {
61  variants_.resize(numVariants(type));
62  for (size_t i = 0; i < numVariants(type); i++)
63  variants_[i] = FiniteElement(i);
64  }
65  else
66  {
67  // for mixed grids add offset for cubes
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);
73  }
74 
75 
76  // compute orientation for all elements
77  const auto& indexSet = gv.indexSet();
78 
79  for(const auto& element : elements(gv))
80  {
81  const auto& refElement = referenceElement(element);
82  auto elementIndex = elementMapper_.index(element);
83  orientation_[elementIndex] = 0;
84 
85  for (std::size_t i=0; i<element.subEntities(dim-1); i++)
86  {
87  // Local vertex indices within the element
88  auto localV0 = refElement.subEntity(i,dim-1, 0,dim);
89  auto localV1 = refElement.subEntity(i,dim-1, 1,dim);
90 
91  // Global vertex indices within the grid
92  auto globalV0 = indexSet.subIndex(element,localV0,dim);
93  auto globalV1 = indexSet.subIndex(element,localV1,dim);
94 
95  if ( (localV0<localV1 && globalV0>globalV1) || (localV0>localV1 && globalV0<globalV1) )
96  orientation_[elementIndex] |= (1 << i);
97  }
98  // for mixed grids add offset for cubes
99  if constexpr (!hasFixedElementType)
100  if (element.type().isCube())
101  orientation_[elementIndex] += numVariants(GeometryTypes::simplex(dim));
102  }
103  }
104 
105  template<class Element>
106  const auto& find(const Element& element) const
107  {
108  return variants_[orientation_[elementMapper_.index(element)]];
109  }
110 
111  private:
112  std::vector<FiniteElement> variants_;
113  const Dune::MultipleCodimMultipleGeomTypeMapper<GV> elementMapper_;
114  std::vector<unsigned short> orientation_;
115  };
116 
117 
118 } // namespace Impl
119 
120 
121 // *****************************************************************************
122 // This is the reusable part of the basis. It contains
123 //
124 // NedelecPreBasis
125 // NedelecNode
126 //
127 // The pre-basis allows to create the others and is the owner of possible shared
128 // state. These components do _not_ depend on the global basis and local view
129 // and can be used without a global basis.
130 // *****************************************************************************
131 
132 template<typename GV, typename Range, std::size_t kind, int order>
133 class NedelecNode;
134 
135 template<typename GV, typename Range, std::size_t kind, int order, class MI>
137 {
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>;
141 
142  using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
143 public:
144 
146  using GridView = GV;
147  using size_type = std::size_t;
148 
150 
152  using IndexSet = Impl::DefaultNodeIndexSet<NedelecPreBasis>;
153 
155  using MultiIndex = MI;
156 
157  using SizePrefix = Dune::ReservedVector<size_type, 1>;
158 
161  gridView_(gv),
162  finiteElementMap_(gv),
163  mapper_(gridView_, mcmgLayout(Dim<1>{}))
164  {
165  if (kind!=1)
166  DUNE_THROW(NotImplemented, "Only Nedelec elements of the first kind are implemented!");
167 
168  // There is no inherent reason why the basis shouldn't work for grids with more than two
169  // element types. Somebody simply has to sit down and implement the missing bits.
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");
172 
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.");
176 
177  if (order>1)
178  DUNE_THROW(NotImplemented, "Only first-order elements are implemented");
179 
180  if (dim!=2 && dim!=3)
181  DUNE_THROW(NotImplemented, "Only 2d and 3d Nédélec elements are implemented");
182  }
183 
185  {}
186 
189  const GridView& gridView() const
190  {
191  return gridView_;
192  }
193 
194  /* \brief Update the stored grid view, to be called if the grid has changed */
195  void update (const GridView& gv)
196  {
197  gridView_ = gv;
198  mapper_.update(gridView_);
199  }
200 
204  Node makeNode() const
205  {
206  return Node{&finiteElementMap_};
207  }
208 
216  [[deprecated("Warning: The IndexSet typedef and the makeIndexSet method are deprecated. "\
217  "As a replacement use the indices() method of the PreBasis directly.")]]
219  {
220  return IndexSet{*this};
221  }
222 
223  size_type size() const
224  {
225  return mapper_.size();
226  }
227 
229  size_type size(const SizePrefix prefix) const
230  {
231  assert(prefix.size() == 0 || prefix.size() == 1);
232  return (prefix.size() == 0) ? size() : 0;
233  }
234 
236  {
237  return size();
238  }
239 
241  {
242  size_type result = 0;
243  for (auto&& type : gridView_.indexSet().types(0))
244  {
245  size_type numEdges = referenceElement<typename GV::ctype,dim>(type).size(dim-1);
246  result = std::max(result, numEdges);
247  }
248 
249  return result;
250  }
251 
255  template<typename It>
256  It indices(const Node& node, It it) const
257  {
258  const auto& element = node.element();
259 
260  // throw if Element is not of predefined type
261  if (not(element.type().isCube()) and not(element.type().isSimplex()))
262  DUNE_THROW(NotImplemented, "NedelecBasis only implemented for cube and simplex elements.");
263 
264  for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
265  {
266  Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
267  *it = { mapper_.subIndex(element, localKey.subEntity(), localKey.codim()) + localKey.index() };
268  }
269 
270  return it;
271  }
272 
273 protected:
275  FiniteElementMap finiteElementMap_;
276  Mapper mapper_;
277 };
278 
279 
280 
281 template<typename GV, typename Range, size_t kind, int order>
282 class NedelecNode :
283  public LeafBasisNode
284 {
285  static const int dim = GV::dimension;
286 
287 public:
288 
289  using size_type = std::size_t;
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,
295  Element>;
296 
297  NedelecNode(const FiniteElementMap* finiteElementMap) :
298  element_(nullptr),
299  finiteElementMap_(finiteElementMap)
300  { }
301 
303  const Element& element() const
304  {
305  return *element_;
306  }
307 
313  {
314  return finiteElement_;
315  }
316 
318  void bind(const Element& e)
319  {
320  element_ = &e;
321  finiteElement_.bind((finiteElementMap_->find(*element_)), e);
322  this->setSize(finiteElement_.size());
323  }
324 
325 protected:
326 
330 };
331 
332 
333 
334 namespace BasisFactory {
335 
336 namespace Impl {
337 
338 template<std::size_t kind, std::size_t order, typename Range>
339 class NedelecPreBasisFactory
340 {
341 public:
342  static const std::size_t requiredMultiIndexSize=1;
343 
344  template<class MultiIndex, class GridView>
345  auto makePreBasis(const GridView& gridView) const
346  {
348  }
349 
350 };
351 
352 } // end namespace BasisFactory::Impl
353 
363 template<std::size_t kind, std::size_t order, typename Range=double>
364 auto nedelec()
365 {
366  return Impl::NedelecPreBasisFactory<kind, order, Range>();
367 }
368 
369 } // end namespace BasisFactory
370 
371 
372 
373 // *****************************************************************************
374 // This is the actual global basis implementation based on the reusable parts.
375 // *****************************************************************************
376 
384 template<typename GV, std::size_t kind, std::size_t order, typename Range=double>
386 
387 } // end namespace Dune::Functions
388 
389 
390 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_NEDELECBASIS_HH
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
Definition: nodes.hh:184