dune-functions  2.6-dev
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/localfunctions/raviartthomas.hh>
10 #include <dune/localfunctions/raviartthomas/raviartthomas0cube2d.hh>
11 #include <dune/localfunctions/raviartthomas/raviartthomas0cube3d.hh>
12 #include <dune/localfunctions/raviartthomas/raviartthomas02d.hh>
13 #include <dune/localfunctions/raviartthomas/raviartthomas1cube2d.hh>
14 #include <dune/localfunctions/raviartthomas/raviartthomas1cube3d.hh>
15 #include <dune/localfunctions/raviartthomas/raviartthomas12d.hh>
16 #include <dune/localfunctions/raviartthomas/raviartthomas2cube2d.hh>
17 
18 #include <dune/typetree/leafnode.hh>
19 
23 
24 namespace Dune {
25 namespace Functions {
26 
27 namespace Impl {
28 
29  template<int dim, GeometryType::BasicType basic_type, typename D, typename R, std::size_t k>
30  struct RaviartThomasLocalInfo
31  {
32  static_assert((AlwaysFalse<D>::value),"The requested type of Raviart-Thomas element is not implemented, sorry!");
33  };
34 
35  template<typename D, typename R>
36  struct RaviartThomasLocalInfo<2,GeometryType::simplex,D,R,0>
37  {
38  using FiniteElement = RT02DLocalFiniteElement<D,R>;
39  static const std::size_t Variants = 8;
40  };
41 
42  template<typename D, typename R>
43  struct RaviartThomasLocalInfo<2,GeometryType::simplex,D,R,1>
44  {
45  using FiniteElement = RT12DLocalFiniteElement<D,R>;
46  static const std::size_t Variants = 8;
47  };
48 
49  template<typename D, typename R>
50  struct RaviartThomasLocalInfo<2,GeometryType::cube,D,R,0>
51  {
52  using FiniteElement = RT0Cube2DLocalFiniteElement<D,R>;
53  static const std::size_t Variants = 16;
54  };
55 
56  template<typename D, typename R>
57  struct RaviartThomasLocalInfo<2,GeometryType::cube,D,R,1>
58  {
59  using FiniteElement = RT1Cube2DLocalFiniteElement<D,R>;
60  static const std::size_t Variants = 16;
61  };
62 
63  template<typename D, typename R>
64  struct RaviartThomasLocalInfo<2,GeometryType::cube,D,R,2>
65  {
66  using FiniteElement = RT2Cube2DLocalFiniteElement<D,R>;
67  static const std::size_t Variants = 16;
68  };
69 
70  template<typename D, typename R>
71  struct RaviartThomasLocalInfo<3,GeometryType::cube,D,R,0>
72  {
73  using FiniteElement = RT0Cube3DLocalFiniteElement<D,R>;
74  static const std::size_t Variants = 64;
75  };
76 
77  template<typename D, typename R>
78  struct RaviartThomasLocalInfo<3,GeometryType::cube,D,R,1>
79  {
80  using FiniteElement = RT1Cube3DLocalFiniteElement<D,R>;
81  static const std::size_t Variants = 64;
82  };
83 
84  template<typename GV, int dim, GeometryType::BasicType basic_type, typename D, typename R, std::size_t k>
85  class RaviartThomasLocalFiniteElementMap
86  {
87  static const std::size_t Variants = RaviartThomasLocalInfo<dim, basic_type, D, R, k>::Variants;
88 
89  public:
90 
91  using FiniteElement = typename RaviartThomasLocalInfo<dim, basic_type, D, R, k>::FiniteElement;
92 
93  RaviartThomasLocalFiniteElementMap(const GV& gv)
94  : gv_(gv), is_(&(gv_.indexSet())), orient_(gv.size(0))
95  {
96  // create all variants
97  for (size_t i = 0; i < Variants; i++)
98  variant_[i] = FiniteElement(i);
99 
100  // compute orientation for all elements
101  // loop once over the grid
102  for(const auto& cell : elements(gv))
103  {
104  unsigned int myId = is_->index(cell);
105  orient_[myId] = 0;
106 
107  for (const auto& intersection : intersections(gv,cell))
108  {
109  if (intersection.neighbor() && (is_->index(intersection.outside()) > myId))
110  orient_[myId] |= (1 << intersection.indexInInside());
111  }
112  }
113  }
114 
116  template<class EntityType>
117  const FiniteElement& find(const EntityType& e) const
118  {
119  return variant_[orient_[is_->index(e)]];
120  }
121 
122  private:
123  GV gv_;
124  std::array<FiniteElement,Variants> variant_;
125  const typename GV::IndexSet* is_;
126  std::vector<unsigned char> orient_;
127  };
128 
129 
130 } // namespace Impl
131 
132 
133 // *****************************************************************************
134 // This is the reusable part of the basis. It contains
135 //
136 // RaviartThomasPreBasis
137 // RaviartThomasNodeIndexSet
138 // RaviartThomasNode
139 //
140 // The pre-basis allows to create the others and is the owner of possible shared
141 // state. These three components do _not_ depend on the global basis or index
142 // set and can be used without a global basis.
143 // *****************************************************************************
144 
145 template<typename GV, int k, typename ST, typename TP, GeometryType::BasicType basic_type>
147 
148 template<typename GV, int k, class MI, class TP, class ST, GeometryType::BasicType basic_type>
150 
151 template<typename GV, int k, class MI, class ST, GeometryType::BasicType basic_type>
153 
154 template<typename GV, int k, class MI, class ST, GeometryType::BasicType basic_type>
156 {
157  static const int dim = GV::dimension;
158  using FiniteElementMap = typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, basic_type, typename GV::ctype, double, k>;
159 
160 private:
161 
162  template<typename, int, class, class, class, GeometryType::BasicType>
164 
165 public:
166 
168  using GridView = GV;
169  using size_type = ST;
170 
171  // Precompute the number of dofs per entity type depending on the entity's codimension and the grid's type (only valid for cube and simplex grids)
172 
173  // for 3D only for cubes k=0,1
174  // Note: dofsPerElement = dofsPerFace * dim for cubes, k=0,1
175 
176  const static int dofsPerFace = dim == 2 ? k+1 : 3*k+1;
177  const static int dofsPerElement = dim == 2 ? (basic_type == GeometryType::cube ? k*(k+1)*dim : k*dim) : k*(k+1)*(k+1)*dim;
178 
179  const std::vector<int> dofsPerCodim {dofsPerElement, dofsPerFace};
180 
181  template<class TP>
183 
184  template<class TP>
186 
188  using MultiIndex = MI;
189 
190  using SizePrefix = Dune::ReservedVector<size_type, 1>;
191 
194  gridView_(gv),
196  { }
197 
199  {
200  codimOffset_.resize(2);
201  codimOffset_[0] = 0;
202  codimOffset_[1] = codimOffset_[0] + dofsPerCodim[0] * gridView_.size(0);
203  }
204 
207  const GridView& gridView() const
208  {
209  return gridView_;
210  }
211 
212  /* \brief Update the stored grid view, to be called if the grid has changed */
213  void update (const GridView& gv)
214  {
215  gridView_ = gv;
216  }
217 
218  template<class TP>
219  Node<TP> node(const TP& tp) const
220  {
221  return Node<TP>{tp, &finiteElementMap_};
222  }
223 
224  template<class TP>
226  {
227  return IndexSet<TP>{*this};
228  }
229 
230  size_type size() const
231  {
232  return dofsPerCodim[0] * gridView_.size(0) + dofsPerCodim[1] * gridView_.size(1);
233  }
234 
236  size_type size(const SizePrefix prefix) const
237  {
238  assert(prefix.size() == 0 || prefix.size() == 1);
239  return (prefix.size() == 0) ? size() : 0;
240  }
241 
244  {
245  return size();
246  }
247 
249  {
250  return StaticPower<(k+1),GV::dimension>::power;
251  }
252 
253 protected:
255  std::vector<size_t> codimOffset_;
256  FiniteElementMap finiteElementMap_;
257 };
258 
259 
260 
261 template<typename GV, int k, typename ST, typename TP, GeometryType::BasicType basic_type>
262 class RaviartThomasNode :
263  public LeafBasisNode<ST, TP>
264 {
265  static const int dim = GV::dimension;
266  static const int maxSize = StaticPower<(k+1),GV::dimension>::power;
267 
268  using Base = LeafBasisNode<ST,TP>;
269 
270 public:
271 
272  using size_type = ST;
273  using TreePath = TP;
274  using Element = typename GV::template Codim<0>::Entity;
275  using FiniteElementMap = typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, basic_type, typename GV::ctype, double, k>;
276  using FiniteElement = typename FiniteElementMap::FiniteElement;
277 
278  RaviartThomasNode(const TreePath& treePath, const FiniteElementMap* finiteElementMap) :
279  Base(treePath),
280  finiteElement_(nullptr),
281  element_(nullptr),
282  finiteElementMap_(finiteElementMap)
283  { }
284 
286  const Element& element() const
287  {
288  return *element_;
289  }
290 
296  {
297  return *finiteElement_;
298  }
299 
301  void bind(const Element& e)
302  {
303  element_ = &e;
305  this->setSize(finiteElement_->size());
306  }
307 
308 protected:
309 
313 };
314 
315 template<typename GV, int k, class MI, class TP, class ST, GeometryType::BasicType basic_type>
317 {
318  enum {dim = GV::dimension};
319 
320 public:
321 
322  using size_type = ST;
323 
325  using MultiIndex = MI;
326 
328 
329  using Node = typename PreBasis::template Node<TP>;
330 
332  preBasis_(&preBasis)
333  {}
334 
340  void bind(const Node& node)
341  {
342  node_ = &node;
343  }
344 
347  void unbind()
348  {
349  node_ = nullptr;
350  }
351 
354  size_type size() const
355  {
356  return node_->finiteElement().size();
357  }
358 
364  template<typename It>
365  It indices(It it) const
366  {
367  const auto& gridIndexSet = preBasis_->gridView().indexSet();
368  const auto& element = node_->element();
369 
370  // throw if Element is not of predefined type
371  if (not(basic_type==GeometryType::BasicType::cube and element.type().isCube()) and
372  not(basic_type==GeometryType::BasicType::simplex and element.type().isSimplex())) DUNE_THROW(Dune::NotImplemented, "RaviartThomasNodalBasis only implemented for cube and simplex elements.");
373 
374  for(std::size_t i=0, end=size(); i<end; ++i, ++it)
375  {
376  Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
377 
378  // The dimension of the entity that the current dof is related to
379  size_t subentity = localKey.subEntity();
380  size_t codim = localKey.codim();
381 
382  if (not(codim==0 or codim==1)) DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the RaviartThomasBasis");
383 
384  *it = { preBasis_->codimOffset_[codim] +
385  preBasis_->dofsPerCodim[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
386  }
387 
388  return it;
389  }
390 
391 protected:
393  const Node* node_;
394 };
395 
396 
397 
398 namespace BasisBuilder {
399 
400 namespace Imp {
401 
402 template<std::size_t k, GeometryType::BasicType basic_type, class size_type=std::size_t>
403 class RaviartThomasPreBasisFactory
404 {
405 public:
406  static const std::size_t requiredMultiIndexSize=1;
407 
408  template<class MultiIndex, class GridView>
409  auto makePreBasis(const GridView& gridView) const
410  {
412  }
413 
414 };
415 
416 } // end namespace BasisBuilder::Imp
417 
426 template<std::size_t k, GeometryType::BasicType basic_type, class size_type=std::size_t>
427 auto rt()
428 {
429  return Imp::RaviartThomasPreBasisFactory<k, basic_type, size_type>();
430 }
431 
432 } // end namespace BasisBuilder
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, GeometryType::BasicType basic_type, class ST = std::size_t>
449 
450 } // end namespace Functions
451 } // end namespace Dune
452 
453 
454 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
Dune::Functions::RaviartThomasPreBasis::finiteElementMap_
FiniteElementMap finiteElementMap_
Definition: raviartthomasbasis.hh:256
Dune::Functions::RaviartThomasNode::finiteElementMap_
const FiniteElementMap * finiteElementMap_
Definition: raviartthomasbasis.hh:312
Dune
Definition: polynomial.hh:7
Dune::Functions::RaviartThomasNode::Element
typename GV::template Codim< 0 >::Entity Element
Definition: raviartthomasbasis.hh:274
Dune::Functions::RaviartThomasNodeIndexSet::RaviartThomasNodeIndexSet
RaviartThomasNodeIndexSet(const PreBasis &preBasis)
Definition: raviartthomasbasis.hh:331
Dune::Functions::RaviartThomasPreBasis::update
void update(const GridView &gv)
Definition: raviartthomasbasis.hh:213
Dune::Functions::RaviartThomasPreBasis::size
size_type size(const SizePrefix prefix) const
Return number possible values for next position in multi index.
Definition: raviartthomasbasis.hh:236
flatmultiindex.hh
Dune::Functions::RaviartThomasPreBasis::dofsPerFace
const static int dofsPerFace
Definition: raviartthomasbasis.hh:176
Dune::Functions::RaviartThomasNode::element_
const Element * element_
Definition: raviartthomasbasis.hh:311
Dune::Functions::RaviartThomasNode::FiniteElementMap
typename Impl::RaviartThomasLocalFiniteElementMap< GV, dim, basic_type, typename GV::ctype, double, k > FiniteElementMap
Definition: raviartthomasbasis.hh:275
Dune::Functions::RaviartThomasPreBasis::dofsPerCodim
const std::vector< int > dofsPerCodim
Definition: raviartthomasbasis.hh:179
Dune::Functions::RaviartThomasNode::bind
void bind(const Element &e)
Bind to element.
Definition: raviartthomasbasis.hh:301
Dune::Functions::RaviartThomasNode::size_type
ST size_type
Definition: raviartthomasbasis.hh:272
Dune::Functions::RaviartThomasNode::FiniteElement
typename FiniteElementMap::FiniteElement FiniteElement
Definition: raviartthomasbasis.hh:276
nodes.hh
defaultglobalbasis.hh
Dune::Functions::RaviartThomasPreBasis::size
size_type size() const
Definition: raviartthomasbasis.hh:230
Dune::Functions::RaviartThomasPreBasis::maxNodeSize
size_type maxNodeSize() const
Definition: raviartthomasbasis.hh:248
Dune::Functions::RaviartThomasNodeIndexSet::MultiIndex
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: raviartthomasbasis.hh:325
Dune::Functions::LeafBasisNode
Definition: nodes.hh:189
Dune::Functions::RaviartThomasNode::finiteElement
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: raviartthomasbasis.hh:295
Dune::Functions::RaviartThomasNode::RaviartThomasNode
RaviartThomasNode(const TreePath &treePath, const FiniteElementMap *finiteElementMap)
Definition: raviartthomasbasis.hh:278
Dune::Functions::BasisBuilder::rt
auto rt()
Create a pre-basis factory that can create a Raviart-Thomas pre-basis.
Definition: raviartthomasbasis.hh:427
Dune::Functions::RaviartThomasNodeIndexSet::node_
const Node * node_
Definition: raviartthomasbasis.hh:393
Dune::Functions::RaviartThomasPreBasis::dofsPerElement
const static int dofsPerElement
Definition: raviartthomasbasis.hh:177
Dune::Functions::RaviartThomasNodeIndexSet::indices
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:365
Dune::Functions::RaviartThomasPreBasis::codimOffset_
std::vector< size_t > codimOffset_
Definition: raviartthomasbasis.hh:255
Dune::Functions::RaviartThomasNodeIndexSet::size_type
ST size_type
Definition: raviartthomasbasis.hh:322
Dune::Functions::RaviartThomasNodeIndexSet::unbind
void unbind()
Unbind the view.
Definition: raviartthomasbasis.hh:347
Dune::Functions::RaviartThomasPreBasis::MultiIndex
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: raviartthomasbasis.hh:188
Dune::Functions::RaviartThomasPreBasis::dimension
size_type dimension() const
Definition: raviartthomasbasis.hh:243
Dune::Functions::RaviartThomasNode::TreePath
TP TreePath
Definition: raviartthomasbasis.hh:273
Dune::Functions::RaviartThomasPreBasis::RaviartThomasPreBasis
RaviartThomasPreBasis(const GridView &gv)
Constructor for a given grid view object.
Definition: raviartthomasbasis.hh:193
Dune::Functions::RaviartThomasNode::element
const Element & element() const
Return current element, throw if unbound.
Definition: raviartthomasbasis.hh:286
Dune::Functions::RaviartThomasPreBasis::GridView
GV GridView
The grid view that the FE space is defined on.
Definition: raviartthomasbasis.hh:168
Dune::Functions::RaviartThomasNodeIndexSet::size
size_type size() const
Size of subtree rooted in this node (element-local)
Definition: raviartthomasbasis.hh:354
Dune::Functions::RaviartThomasPreBasis::gridView
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: raviartthomasbasis.hh:207
Dune::Functions::RaviartThomasNodeIndexSet
Definition: raviartthomasbasis.hh:149
Dune::Functions::RaviartThomasPreBasis::size_type
ST size_type
Definition: raviartthomasbasis.hh:169
Dune::Functions::RaviartThomasNodeIndexSet::bind
void bind(const Node &node)
Bind the view to a grid element.
Definition: raviartthomasbasis.hh:340
Dune::Functions::RaviartThomasNodeIndexSet::preBasis_
const PreBasis * preBasis_
Definition: raviartthomasbasis.hh:392
Dune::Functions::RaviartThomasPreBasis::initializeIndices
void initializeIndices()
Definition: raviartthomasbasis.hh:198
Dune::Functions::RaviartThomasNode
Definition: raviartthomasbasis.hh:146
Dune::Functions::RaviartThomasPreBasis::SizePrefix
Dune::ReservedVector< size_type, 1 > SizePrefix
Definition: raviartthomasbasis.hh:190
Dune::Functions::RaviartThomasNode::finiteElement_
const FiniteElement * finiteElement_
Definition: raviartthomasbasis.hh:310
Dune::Functions::BasisBuilder::power
auto power(ChildPreBasisFactory &&childPreBasisFactory, const IndexMergingStrategy &ims)
Create a pre-basis factory that can build a PowerPreBasis.
Definition: powerbasis.hh:493
Dune::Functions::RaviartThomasPreBasis::node
Node< TP > node(const TP &tp) const
Definition: raviartthomasbasis.hh:219
Dune::Functions::RaviartThomasPreBasis::indexSet
IndexSet< TP > indexSet() const
Definition: raviartthomasbasis.hh:225
Dune::Functions::RaviartThomasPreBasis::gridView_
const GridView gridView_
Definition: raviartthomasbasis.hh:254
Dune::Functions::DefaultGlobalBasis
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:42
Dune::Functions::RaviartThomasPreBasis
Definition: raviartthomasbasis.hh:152
Dune::Functions::RaviartThomasNodeIndexSet::Node
typename PreBasis::template Node< TP > Node
Definition: raviartthomasbasis.hh:329