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