dune-functions  2.7.1
raviartthomasbasis.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_RAVIARTTHOMASBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
5 
6 #include <array>
7 #include <dune/common/exceptions.hh>
8 
9 #include <dune/grid/common/capabilities.hh>
10 
11 #include <dune/localfunctions/common/virtualinterface.hh>
12 #include <dune/localfunctions/common/virtualwrappers.hh>
13 
14 #include <dune/localfunctions/raviartthomas.hh>
15 #include <dune/localfunctions/raviartthomas/raviartthomas0cube2d.hh>
16 #include <dune/localfunctions/raviartthomas/raviartthomas0cube3d.hh>
17 #include <dune/localfunctions/raviartthomas/raviartthomas02d.hh>
18 #include <dune/localfunctions/raviartthomas/raviartthomas1cube2d.hh>
19 #include <dune/localfunctions/raviartthomas/raviartthomas1cube3d.hh>
20 #include <dune/localfunctions/raviartthomas/raviartthomas12d.hh>
21 #include <dune/localfunctions/raviartthomas/raviartthomas2cube2d.hh>
22 
23 #include <dune/typetree/leafnode.hh>
24 
28 
29 namespace Dune {
30 namespace Functions {
31 
32 namespace Impl {
33 
34  template<int dim, typename D, typename R, std::size_t k>
35  struct RaviartThomasSimplexLocalInfo
36  {
37  static_assert((AlwaysFalse<D>::value),"The requested type of Raviart-Thomas element is not implemented, sorry!");
38  };
39 
40  template<typename D, typename R>
41  struct RaviartThomasSimplexLocalInfo<2,D,R,0>
42  {
43  using FiniteElement = RT02DLocalFiniteElement<D,R>;
44  static const std::size_t Variants = 8;
45  };
46 
47  template<typename D, typename R>
48  struct RaviartThomasSimplexLocalInfo<2,D,R,1>
49  {
50  using FiniteElement = RT12DLocalFiniteElement<D,R>;
51  static const std::size_t Variants = 8;
52  };
53 
54  template<int dim, typename D, typename R, std::size_t k>
55  struct RaviartThomasCubeLocalInfo
56  {
57  static_assert((AlwaysFalse<D>::value),"The requested type of Raviart-Thomas element is not implemented, sorry!");
58  };
59 
60  template<typename D, typename R>
61  struct RaviartThomasCubeLocalInfo<2,D,R,0>
62  {
63  using FiniteElement = RT0Cube2DLocalFiniteElement<D,R>;
64  static const std::size_t Variants = 16;
65  };
66 
67  template<typename D, typename R>
68  struct RaviartThomasCubeLocalInfo<2,D,R,1>
69  {
70  using FiniteElement = RT1Cube2DLocalFiniteElement<D,R>;
71  static const std::size_t Variants = 16;
72  };
73 
74  template<typename D, typename R>
75  struct RaviartThomasCubeLocalInfo<2,D,R,2>
76  {
77  using FiniteElement = RT2Cube2DLocalFiniteElement<D,R>;
78  static const std::size_t Variants = 16;
79  };
80 
81  template<typename D, typename R>
82  struct RaviartThomasCubeLocalInfo<3,D,R,0>
83  {
84  using FiniteElement = RT0Cube3DLocalFiniteElement<D,R>;
85  static const std::size_t Variants = 64;
86  };
87 
88  template<typename D, typename R>
89  struct RaviartThomasCubeLocalInfo<3,D,R,1>
90  {
91  using FiniteElement = RT1Cube3DLocalFiniteElement<D,R>;
92  static const std::size_t Variants = 64;
93  };
94 
95  template<typename GV, int dim, typename R, std::size_t k>
96  class RaviartThomasLocalFiniteElementMap
97  {
98  using D = typename GV::ctype;
99  constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
100 
101  using CubeFiniteElement = typename RaviartThomasCubeLocalInfo<dim, D, R, k>::FiniteElement;
102  using SimplexFiniteElement = typename RaviartThomasSimplexLocalInfo<dim, D, R, k>::FiniteElement;
103  using CubeFiniteElementImp = typename std::conditional<hasFixedElementType,
104  CubeFiniteElement,
105  LocalFiniteElementVirtualImp<CubeFiniteElement> >::type;
106  using SimplexFiniteElementImp = typename std::conditional<hasFixedElementType,
107  SimplexFiniteElement,
108  LocalFiniteElementVirtualImp<SimplexFiniteElement> >::type;
109 
110  public:
111 
112  using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
113 
114  constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId; // meaningless if hasFixedElementType is false
115  constexpr static GeometryType type = GeometryType(topologyId, GV::dimension);
116 
117  using FiniteElement = typename std::conditional<hasFixedElementType,
118  typename std::conditional<type.isCube(),CubeFiniteElement,SimplexFiniteElement>::type,
119  LocalFiniteElementVirtualInterface<T> >::type;
120 
121  RaviartThomasLocalFiniteElementMap(const GV& gv)
122  : is_(&(gv.indexSet())), orient_(gv.size(0))
123  {
124  cubeVariant_.resize(RaviartThomasCubeLocalInfo<dim, D, R, k>::Variants);
125  simplexVariant_.resize(RaviartThomasSimplexLocalInfo<dim, D, R, k>::Variants);
126 
127  // create all variants
128  for (size_t i = 0; i < cubeVariant_.size(); i++)
129  cubeVariant_[i] = std::make_unique<CubeFiniteElementImp>(CubeFiniteElement(i));
130 
131  for (size_t i = 0; i < simplexVariant_.size(); i++)
132  simplexVariant_[i] = std::make_unique<SimplexFiniteElementImp>(SimplexFiniteElement(i));
133 
134  // compute orientation for all elements
135  // loop once over the grid
136  for(const auto& cell : elements(gv))
137  {
138  unsigned int myId = is_->index(cell);
139  orient_[myId] = 0;
140 
141  for (const auto& intersection : intersections(gv,cell))
142  {
143  if (intersection.neighbor() && (is_->index(intersection.outside()) > myId))
144  orient_[myId] |= (1 << intersection.indexInInside());
145  }
146  }
147  }
148 
154  template<class EntityType,
155  std::enable_if_t<!hasFixedElementType and AlwaysTrue<EntityType>::value,int> = 0>
156  const FiniteElement& find(const EntityType& e) const
157  {
158  if (e.type().isCube())
159  return *cubeVariant_[orient_[is_->index(e)]];
160  else
161  return *simplexVariant_[orient_[is_->index(e)]];
162  }
163 
169  template<class EntityType,
170  std::enable_if_t<hasFixedElementType and type.isCube() and AlwaysTrue<EntityType>::value,int> = 0>
171  const FiniteElement& find(const EntityType& e) const
172  {
173  return *cubeVariant_[orient_[is_->index(e)]];
174  }
175 
181  template<class EntityType,
182  std::enable_if_t<hasFixedElementType and type.isSimplex() and AlwaysTrue<EntityType>::value,int> = 0>
183  const FiniteElement& find(const EntityType& e) const
184  {
185  return *simplexVariant_[orient_[is_->index(e)]];
186  }
187 
188  private:
189  std::vector<std::unique_ptr<CubeFiniteElementImp> > cubeVariant_;
190  std::vector<std::unique_ptr<SimplexFiniteElementImp> > simplexVariant_;
191  const typename GV::IndexSet* is_;
192  std::vector<unsigned char> orient_;
193  };
194 
195 
196 } // namespace Impl
197 
198 
199 // *****************************************************************************
200 // This is the reusable part of the basis. It contains
201 //
202 // RaviartThomasPreBasis
203 // RaviartThomasNodeIndexSet
204 // RaviartThomasNode
205 //
206 // The pre-basis allows to create the others and is the owner of possible shared
207 // state. These three components do _not_ depend on the global basis or index
208 // set and can be used without a global basis.
209 // *****************************************************************************
210 
211 template<typename GV, int k>
212 class RaviartThomasNode;
213 
214 template<typename GV, int k, class MI>
215 class RaviartThomasNodeIndexSet;
216 
217 template<typename GV, int k, class MI>
219 {
220  static const int dim = GV::dimension;
221  using FiniteElementMap = typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
222 
223  template<typename, int, class>
225 
226 public:
227 
229  using GridView = GV;
230  using size_type = std::size_t;
231 
233 
235 
237  using MultiIndex = MI;
238 
239  using SizePrefix = Dune::ReservedVector<size_type, 1>;
240 
243  gridView_(gv),
245  {
246  // There is no inherent reason why the basis shouldn't work for grids with more than one
247  // element types. Somebody simply has to sit down and implement the missing bits.
248  if (gv.indexSet().types(0).size() > 1)
249  DUNE_THROW(Dune::NotImplemented, "Raviart-Thomas basis is only implemented for grids with a single element type");
250 
251  GeometryType type = gv.template begin<0>()->type();
252  const static int dofsPerElement = (dim == 2) ? (type.isCube() ? k*(k+1)*dim : k*dim) : k*(k+1)*(k+1)*dim;
253  constexpr int dofsPerFace = (dim == 2) ? k+1 : 3*k+1;
254 
255  dofsPerCodim_ = {{dofsPerElement, dofsPerFace}};
256  }
257 
259  {
260  codimOffset_[0] = 0;
261  codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
262  }
263 
266  const GridView& gridView() const
267  {
268  return gridView_;
269  }
270 
271  /* \brief Update the stored grid view, to be called if the grid has changed */
272  void update (const GridView& gv)
273  {
274  gridView_ = gv;
275  }
276 
280  Node makeNode() const
281  {
282  return Node{&finiteElementMap_};
283  }
284 
292  {
293  return IndexSet{*this};
294  }
295 
296  size_type size() const
297  {
298  return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1);
299  }
300 
302  size_type size(const SizePrefix prefix) const
303  {
304  assert(prefix.size() == 0 || prefix.size() == 1);
305  return (prefix.size() == 0) ? size() : 0;
306  }
307 
310  {
311  return size();
312  }
313 
315  {
316  size_type result = 0;
317  for (auto&& type : gridView_.indexSet().types(0))
318  {
319  size_t numFaces = ReferenceElements<double,dim>::general(type).size(1);
320  const static int dofsPerCodim0 = (dim == 2) ? (type.isCube() ? k*(k+1)*dim : k*dim) : k*(k+1)*(k+1)*dim;
321  constexpr int dofsPerCodim1 = (dim == 2) ? k+1 : 3*k+1;
322  result = std::max(result, dofsPerCodim0 + dofsPerCodim1 * numFaces);
323  }
324 
325  return result;
326  }
327 
328 protected:
330  std::array<size_t,dim+1> codimOffset_;
331  FiniteElementMap finiteElementMap_;
332  // Number of dofs per entity type depending on the entity's codimension and type
333  std::array<int,dim+1> dofsPerCodim_;
334 };
335 
336 
337 
338 template<typename GV, int k>
340  public LeafBasisNode
341 {
342  static const int dim = GV::dimension;
343 
344 public:
345 
346  using size_type = std::size_t;
347  using Element = typename GV::template Codim<0>::Entity;
348  using FiniteElementMap = typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
349  using FiniteElement = typename FiniteElementMap::FiniteElement;
350 
351  RaviartThomasNode(const FiniteElementMap* finiteElementMap) :
352  finiteElement_(nullptr),
353  element_(nullptr),
354  finiteElementMap_(finiteElementMap)
355  { }
356 
358  const Element& element() const
359  {
360  return *element_;
361  }
362 
368  {
369  return *finiteElement_;
370  }
371 
373  void bind(const Element& e)
374  {
375  element_ = &e;
377  this->setSize(finiteElement_->size());
378  }
379 
380 protected:
381 
385 };
386 
387 template<typename GV, int k, class MI>
389 {
390  enum {dim = GV::dimension};
391 
392 public:
393 
394  using size_type = std::size_t;
395 
397  using MultiIndex = MI;
398 
400 
402 
404  preBasis_(&preBasis)
405  {}
406 
412  void bind(const Node& node)
413  {
414  node_ = &node;
415  }
416 
419  void unbind()
420  {
421  node_ = nullptr;
422  }
423 
426  size_type size() const
427  {
428  return node_->finiteElement().size();
429  }
430 
436  template<typename It>
437  It indices(It it) const
438  {
439  const auto& gridIndexSet = preBasis_->gridView().indexSet();
440  const auto& element = node_->element();
441 
442  // throw if Element is not of predefined type
443  if (not(element.type().isCube()) and not(element.type().isSimplex()))
444  DUNE_THROW(Dune::NotImplemented, "RaviartThomasBasis only implemented for cube and simplex elements.");
445 
446  for(std::size_t i=0, end=size(); i<end; ++i, ++it)
447  {
448  Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
449 
450  // The dimension of the entity that the current dof is related to
451  size_t subentity = localKey.subEntity();
452  size_t codim = localKey.codim();
453 
454  if (not(codim==0 or codim==1))
455  DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the RaviartThomasBasis");
456 
457  *it = { preBasis_->codimOffset_[codim] +
458  preBasis_->dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
459  }
460 
461  return it;
462  }
463 
464 protected:
466  const Node* node_;
467 };
468 
469 
470 
471 namespace BasisFactory {
472 
473 namespace Imp {
474 
475 template<std::size_t k>
476 class RaviartThomasPreBasisFactory
477 {
478 public:
479  static const std::size_t requiredMultiIndexSize=1;
480 
481  template<class MultiIndex, class GridView>
482  auto makePreBasis(const GridView& gridView) const
483  {
485  }
486 
487 };
488 
489 } // end namespace BasisFactory::Imp
490 
498 template<std::size_t k>
500 {
501  return Imp::RaviartThomasPreBasisFactory<k>();
502 }
503 
504 } // end namespace BasisFactory
505 
506 
507 
508 // *****************************************************************************
509 // This is the actual global basis implementation based on the reusable parts.
510 // *****************************************************************************
511 
519 template<typename GV, int k>
521 
522 } // end namespace Functions
523 } // end namespace Dune
524 
525 
526 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
auto raviartThomas()
Create a pre-basis factory that can create a Raviart-Thomas pre-basis.
Definition: raviartthomasbasis.hh:499
Definition: polynomial.hh:10
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:47
std::size_t size_type
Definition: nodes.hh:125
void setSize(const size_type size)
Definition: nodes.hh:160
Definition: nodes.hh:182
Definition: raviartthomasbasis.hh:341
typename Impl::RaviartThomasLocalFiniteElementMap< GV, dim, double, k > FiniteElementMap
Definition: raviartthomasbasis.hh:348
void bind(const Element &e)
Bind to element.
Definition: raviartthomasbasis.hh:373
const Element & element() const
Return current element, throw if unbound.
Definition: raviartthomasbasis.hh:358
typename GV::template Codim< 0 >::Entity Element
Definition: raviartthomasbasis.hh:347
const Element * element_
Definition: raviartthomasbasis.hh:383
typename FiniteElementMap::FiniteElement FiniteElement
Definition: raviartthomasbasis.hh:349
const FiniteElement * finiteElement_
Definition: raviartthomasbasis.hh:382
RaviartThomasNode(const FiniteElementMap *finiteElementMap)
Definition: raviartthomasbasis.hh:351
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: raviartthomasbasis.hh:367
const FiniteElementMap * finiteElementMap_
Definition: raviartthomasbasis.hh:384
Definition: raviartthomasbasis.hh:389
It indices(It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis.
Definition: raviartthomasbasis.hh:437
RaviartThomasNodeIndexSet(const PreBasis &preBasis)
Definition: raviartthomasbasis.hh:403
void bind(const Node &node)
Bind the view to a grid element.
Definition: raviartthomasbasis.hh:412
const PreBasis * preBasis_
Definition: raviartthomasbasis.hh:465
void unbind()
Unbind the view.
Definition: raviartthomasbasis.hh:419
std::size_t size_type
Definition: raviartthomasbasis.hh:394
size_type size() const
Size of subtree rooted in this node (element-local)
Definition: raviartthomasbasis.hh:426
const Node * node_
Definition: raviartthomasbasis.hh:466
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: raviartthomasbasis.hh:397
Definition: raviartthomasbasis.hh:219
Node makeNode() const
Create tree node.
Definition: raviartthomasbasis.hh:280
std::size_t size_type
Definition: raviartthomasbasis.hh:230
std::array< size_t, dim+1 > codimOffset_
Definition: raviartthomasbasis.hh:330
Dune::ReservedVector< size_type, 1 > SizePrefix
Definition: raviartthomasbasis.hh:239
void update(const GridView &gv)
Definition: raviartthomasbasis.hh:272
void initializeIndices()
Definition: raviartthomasbasis.hh:258
RaviartThomasPreBasis(const GridView &gv)
Constructor for a given grid view object.
Definition: raviartthomasbasis.hh:242
std::array< int, dim+1 > dofsPerCodim_
Definition: raviartthomasbasis.hh:333
GV GridView
The grid view that the FE space is defined on.
Definition: raviartthomasbasis.hh:229
FiniteElementMap finiteElementMap_
Definition: raviartthomasbasis.hh:331
size_type dimension() const
Definition: raviartthomasbasis.hh:309
size_type size() const
Definition: raviartthomasbasis.hh:296
const GridView gridView_
Definition: raviartthomasbasis.hh:329
size_type maxNodeSize() const
Definition: raviartthomasbasis.hh:314
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: raviartthomasbasis.hh:266
IndexSet makeIndexSet() const
Create tree node index set.
Definition: raviartthomasbasis.hh:291
size_type size(const SizePrefix prefix) const
Return number possible values for next position in multi index.
Definition: raviartthomasbasis.hh:302
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: raviartthomasbasis.hh:237