dune-functions  2.9.0
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 #include <dune/grid/common/mcmgmapper.hh>
11 
12 #include <dune/localfunctions/common/localfiniteelementvariant.hh>
13 #include <dune/localfunctions/raviartthomas.hh>
14 #include <dune/localfunctions/raviartthomas/raviartthomas0cube2d.hh>
15 #include <dune/localfunctions/raviartthomas/raviartthomas0cube3d.hh>
16 #include <dune/localfunctions/raviartthomas/raviartthomas02d.hh>
17 #include <dune/localfunctions/raviartthomas/raviartthomas03d.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 
26 
27 namespace Dune {
28 namespace Functions {
29 
30 namespace Impl {
31 
32  template<int dim, typename D, typename R, std::size_t k>
33  struct RaviartThomasSimplexLocalInfo
34  {
35  // Dummy type, must be something that we can have a std::unique_ptr to
36  using FiniteElement = void*;
37  };
38 
39  template<typename D, typename R>
40  struct RaviartThomasSimplexLocalInfo<2,D,R,0>
41  {
42  using FiniteElement = RT02DLocalFiniteElement<D,R>;
43  };
44 
45  template<typename D, typename R>
46  struct RaviartThomasSimplexLocalInfo<2,D,R,1>
47  {
48  using FiniteElement = RT12DLocalFiniteElement<D,R>;
49  };
50 
51  template<typename D, typename R>
52  struct RaviartThomasSimplexLocalInfo<3,D,R,0>
53  {
54  using FiniteElement = RT03DLocalFiniteElement<D,R>;
55  };
56 
57  template<int dim, typename D, typename R, std::size_t k>
58  struct RaviartThomasCubeLocalInfo
59  {
60  // Dummy type, must be something that we can have a std::unique_ptr to
61  using FiniteElement = void*;
62  };
63 
64  template<typename D, typename R>
65  struct RaviartThomasCubeLocalInfo<2,D,R,0>
66  {
67  using FiniteElement = RT0Cube2DLocalFiniteElement<D,R>;
68  };
69 
70  template<typename D, typename R>
71  struct RaviartThomasCubeLocalInfo<2,D,R,1>
72  {
73  using FiniteElement = RT1Cube2DLocalFiniteElement<D,R>;
74  };
75 
76  template<typename D, typename R>
77  struct RaviartThomasCubeLocalInfo<2,D,R,2>
78  {
79  using FiniteElement = RT2Cube2DLocalFiniteElement<D,R>;
80  };
81 
82  template<typename D, typename R>
83  struct RaviartThomasCubeLocalInfo<3,D,R,0>
84  {
85  using FiniteElement = RT0Cube3DLocalFiniteElement<D,R>;
86  };
87 
88  template<typename D, typename R>
89  struct RaviartThomasCubeLocalInfo<3,D,R,1>
90  {
91  using FiniteElement = RT1Cube3DLocalFiniteElement<D,R>;
92  };
93 
94  template<typename GV, int dim, typename R, std::size_t k>
95  class RaviartThomasLocalFiniteElementMap
96  {
97  using D = typename GV::ctype;
98  constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
99 
100  using CubeFiniteElement = typename RaviartThomasCubeLocalInfo<dim, D, R, k>::FiniteElement;
101  using SimplexFiniteElement = typename RaviartThomasSimplexLocalInfo<dim, D, R, k>::FiniteElement;
102 
103  public:
104 
105  using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
106 
107  constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId; // meaningless if hasFixedElementType is false
108  constexpr static GeometryType type = GeometryType(topologyId, GV::dimension);
109 
110  using FiniteElement = std::conditional_t<hasFixedElementType,
111  std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
112  LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
113 
114  // Each element facet can have its orientation reversed, hence there are
115  // 2^#facets different variants.
116  static std::size_t numVariants(GeometryType type)
117  {
118  auto numFacets = referenceElement<D,dim>(type).size(1);
119  return power(2,numFacets);
120  }
121 
122  RaviartThomasLocalFiniteElementMap(const GV& gv)
123  : elementMapper_(gv, mcmgElementLayout()),
124  orient_(gv.size(0))
125  {
126  if constexpr (hasFixedElementType)
127  {
128  variants_.resize(numVariants(type));
129  for (size_t i = 0; i < numVariants(type); i++)
130  variants_[i] = FiniteElement(i);
131  }
132  else
133  {
134  // for mixed grids add offset for cubes
135  variants_.resize(numVariants(GeometryTypes::simplex(dim)) + numVariants(GeometryTypes::cube(dim)));
136  for (size_t i = 0; i < numVariants(GeometryTypes::simplex(dim)); i++)
137  variants_[i] = SimplexFiniteElement(i);
138  for (size_t i = 0; i < numVariants(GeometryTypes::cube(dim)); i++)
139  variants_[i + numVariants(GeometryTypes::simplex(dim))] = CubeFiniteElement(i);
140  }
141 
142  for(const auto& cell : elements(gv))
143  {
144  unsigned int myId = elementMapper_.index(cell);
145  orient_[myId] = 0;
146 
147  for (const auto& intersection : intersections(gv,cell))
148  {
149  if (intersection.neighbor() && (elementMapper_.index(intersection.outside()) > myId))
150  orient_[myId] |= (1 << intersection.indexInInside());
151  }
152 
153  // for mixed grids add offset for cubes
154  if constexpr (!hasFixedElementType)
155  if (cell.type().isCube())
156  orient_[myId] += numVariants(GeometryTypes::simplex(dim));
157  }
158  }
159 
160  template<class EntityType>
161  const FiniteElement& find(const EntityType& e) const
162  {
163  return variants_[orient_[elementMapper_.index(e)]];
164  }
165 
166  private:
167  std::vector<FiniteElement> variants_;
168  const Dune::MultipleCodimMultipleGeomTypeMapper<GV> elementMapper_;
169  std::vector<unsigned char> orient_;
170  };
171 
172 
173 } // namespace Impl
174 
175 
176 // *****************************************************************************
177 // This is the reusable part of the basis. It contains
178 //
179 // RaviartThomasPreBasis
180 // RaviartThomasNode
181 //
182 // The pre-basis allows to create the others and is the owner of possible shared
183 // state. These components do _not_ depend on the global basis and local view
184 // and can be used without a global basis.
185 // *****************************************************************************
186 
187 template<typename GV, int k>
188 class RaviartThomasNode;
189 
190 template<typename GV, int k>
192 {
193  static const int dim = GV::dimension;
194  using FiniteElementMap = typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
195 
196 public:
197 
199  using GridView = GV;
200  using size_type = std::size_t;
201 
203 
204  static constexpr size_type maxMultiIndexSize = 1;
205  static constexpr size_type minMultiIndexSize = 1;
206  static constexpr size_type multiIndexBufferSize = 1;
207 
210  gridView_(gv),
212  {
213  // Currently there are some unresolved bugs with hybrid grids and higher order Raviart-Thomas elements
214  if (gv.indexSet().types(0).size() > 1 and k>0)
215  DUNE_THROW(Dune::NotImplemented, "Raviart-Thomas basis with index k>0 is only implemented for grids with a single element type");
216 
217  for(auto type : gv.indexSet().types(0))
218  if (!type.isSimplex() && !type.isCube())
219  DUNE_THROW(Dune::NotImplemented, "Raviart-Thomas elements are only implemented for grids with simplex or cube elements.");
220 
221  GeometryType type = gv.template begin<0>()->type();
222  const static int dofsPerElement = type.isCube() ? ((dim == 2) ? k*(k+1)*dim : k*(k+1)*(k+1)*dim) : k*dim;
223  const static int dofsPerFace = type.isCube() ? (dim-2)*2*k+k+1 : (dim-1)*k+1 ;
224 
225  dofsPerCodim_ = {{dofsPerElement, dofsPerFace}};
226  }
227 
229  {
230  codimOffset_[0] = 0;
231  codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
232  }
233 
236  const GridView& gridView() const
237  {
238  return gridView_;
239  }
240 
241  /* \brief Update the stored grid view, to be called if the grid has changed */
242  void update (const GridView& gv)
243  {
244  gridView_ = gv;
245  }
246 
250  Node makeNode() const
251  {
252  return Node{&finiteElementMap_};
253  }
254 
255  size_type size() const
256  {
257  return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1);
258  }
259 
261  template<class SizePrefix>
262  size_type size(const SizePrefix& prefix) const
263  {
264  assert(prefix.size() == 0 || prefix.size() == 1);
265  return (prefix.size() == 0) ? size() : 0;
266  }
267 
270  {
271  return size();
272  }
273 
275  {
276  size_type result = 0;
277  for (auto&& type : gridView_.indexSet().types(0))
278  {
279  size_t numFaces = ReferenceElements<double,dim>::general(type).size(1);
280  const static int dofsPerElement = type.isCube() ? ((dim == 2) ? k*(k+1)*dim : k*(k+1)*(k+1)*dim) : k*dim;
281  const static int dofsPerFace = type.isCube() ? (dim-2)*2*k+k+1 : (dim-1)*k+1 ;
282  result = std::max(result, dofsPerElement + dofsPerFace * numFaces);
283  }
284 
285  return result;
286  }
287 
293  template<typename It>
294  It indices(const Node& node, It it) const
295  {
296  const auto& gridIndexSet = gridView().indexSet();
297  const auto& element = node.element();
298 
299  // throw if Element is not of predefined type
300  if (not(element.type().isCube()) and not(element.type().isSimplex()))
301  DUNE_THROW(Dune::NotImplemented, "RaviartThomasBasis only implemented for cube and simplex elements.");
302 
303  for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
304  {
305  Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
306 
307  // The dimension of the entity that the current dof is related to
308  size_t subentity = localKey.subEntity();
309  size_t codim = localKey.codim();
310 
311  if (not(codim==0 or codim==1))
312  DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the RaviartThomasBasis");
313 
314  *it = { codimOffset_[codim] +
315  dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
316  }
317 
318  return it;
319  }
320 
321 protected:
323  std::array<size_t,dim+1> codimOffset_;
324  FiniteElementMap finiteElementMap_;
325  // Number of dofs per entity type depending on the entity's codimension and type
326  std::array<int,dim+1> dofsPerCodim_;
327 };
328 
329 
330 
331 template<typename GV, int k>
333  public LeafBasisNode
334 {
335  static const int dim = GV::dimension;
336 
337 public:
338 
339  using size_type = std::size_t;
340  using Element = typename GV::template Codim<0>::Entity;
341  using FiniteElementMap = typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
342  using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::ContravariantPiolaTransformator,
343  typename FiniteElementMap::FiniteElement,
344  Element>;
345 
346  RaviartThomasNode(const FiniteElementMap* finiteElementMap) :
347  element_(nullptr),
348  finiteElementMap_(finiteElementMap)
349  { }
350 
352  const Element& element() const
353  {
354  return *element_;
355  }
356 
362  {
363  return finiteElement_;
364  }
365 
367  void bind(const Element& e)
368  {
369  element_ = &e;
370  finiteElement_.bind((finiteElementMap_->find(*element_)), e);
371  this->setSize(finiteElement_.size());
372  }
373 
374 protected:
375 
379 };
380 
381 namespace BasisFactory {
382 
390 template<std::size_t k>
392 {
393  return [](const auto& gridView) {
394  return RaviartThomasPreBasis<std::decay_t<decltype(gridView)>, k>(gridView);
395  };
396 }
397 
398 } // end namespace BasisFactory
399 
400 
401 
402 // *****************************************************************************
403 // This is the actual global basis implementation based on the reusable parts.
404 // *****************************************************************************
405 
413 template<typename GV, int k>
415 
416 } // end namespace Functions
417 } // end namespace Dune
418 
419 
420 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
auto power(ChildPreBasisFactory &&childPreBasisFactory, const IndexMergingStrategy &)
Create a pre-basis factory that can build a PowerPreBasis.
Definition: powerbasis.hh:369
auto raviartThomas()
Create a pre-basis factory that can create a Raviart-Thomas pre-basis.
Definition: raviartthomasbasis.hh:391
Definition: polynomial.hh:10
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:46
size_type size() const
Definition: nodes.hh:142
std::size_t size_type
Definition: nodes.hh:128
void setSize(const size_type size)
Definition: nodes.hh:164
Definition: nodes.hh:186
Definition: raviartthomasbasis.hh:334
typename Impl::RaviartThomasLocalFiniteElementMap< GV, dim, double, k > FiniteElementMap
Definition: raviartthomasbasis.hh:341
void bind(const Element &e)
Bind to element.
Definition: raviartthomasbasis.hh:367
Impl::GlobalValuedLocalFiniteElement< Impl::ContravariantPiolaTransformator, typename FiniteElementMap::FiniteElement, Element > FiniteElement
Definition: raviartthomasbasis.hh:344
const Element & element() const
Return current element, throw if unbound.
Definition: raviartthomasbasis.hh:352
typename GV::template Codim< 0 >::Entity Element
Definition: raviartthomasbasis.hh:340
const Element * element_
Definition: raviartthomasbasis.hh:377
RaviartThomasNode(const FiniteElementMap *finiteElementMap)
Definition: raviartthomasbasis.hh:346
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: raviartthomasbasis.hh:361
FiniteElement finiteElement_
Definition: raviartthomasbasis.hh:376
const FiniteElementMap * finiteElementMap_
Definition: raviartthomasbasis.hh:378
Definition: raviartthomasbasis.hh:192
static constexpr size_type minMultiIndexSize
Definition: raviartthomasbasis.hh:205
static constexpr size_type maxMultiIndexSize
Definition: raviartthomasbasis.hh:204
Node makeNode() const
Create tree node.
Definition: raviartthomasbasis.hh:250
std::array< int, dim+1 > dofsPerCodim_
Definition: raviartthomasbasis.hh:326
void update(const GridView &gv)
Definition: raviartthomasbasis.hh:242
RaviartThomasPreBasis(const GridView &gv)
Constructor for a given grid view object.
Definition: raviartthomasbasis.hh:209
std::size_t size_type
Definition: raviartthomasbasis.hh:200
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: raviartthomasbasis.hh:294
static constexpr size_type multiIndexBufferSize
Definition: raviartthomasbasis.hh:206
FiniteElementMap finiteElementMap_
Definition: raviartthomasbasis.hh:324
size_type size() const
Definition: raviartthomasbasis.hh:255
size_type dimension() const
Definition: raviartthomasbasis.hh:269
GV GridView
The grid view that the FE space is defined on.
Definition: raviartthomasbasis.hh:199
size_type maxNodeSize() const
Definition: raviartthomasbasis.hh:274
GridView gridView_
Definition: raviartthomasbasis.hh:322
size_type size(const SizePrefix &prefix) const
Return number possible values for next position in multi index.
Definition: raviartthomasbasis.hh:262
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: raviartthomasbasis.hh:236
void initializeIndices()
Definition: raviartthomasbasis.hh:228
std::array< size_t, dim+1 > codimOffset_
Definition: raviartthomasbasis.hh:323