dune-functions  2.7.1
brezzidouglasmarinibasis.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_BREZZIDOUGLASMARINIBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINIBASIS_HH
5 
6 #include <array>
7 #include <dune/common/exceptions.hh>
8 #include <dune/geometry/referenceelements.hh>
9 
10 #include <dune/localfunctions/common/virtualinterface.hh>
11 #include <dune/localfunctions/common/virtualwrappers.hh>
12 
13 #include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1cube2d.hh>
14 #include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1cube3d.hh>
15 #include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1simplex2d.hh>
16 #include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini2cube2d.hh>
17 #include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini2simplex2d.hh>
18 
19 #include <dune/typetree/leafnode.hh>
20 
24 
25 namespace Dune {
26 namespace Functions {
27 
28 namespace Impl {
29 
30  template<int dim, typename D, typename R, std::size_t k>
31  struct BDMSimplexLocalInfo
32  {
33  static_assert((AlwaysFalse<D>::value),"The requested type of BDM element is not implemented, sorry!");
34  };
35 
36  template<typename D, typename R>
37  struct BDMSimplexLocalInfo<2,D,R,1>
38  {
39  using FiniteElement = BDM1Simplex2DLocalFiniteElement<D,R>;
40  static const std::size_t Variants = 8;
41  };
42 
43  template<typename D, typename R>
44  struct BDMSimplexLocalInfo<2,D,R,2>
45  {
46  using FiniteElement = BDM2Simplex2DLocalFiniteElement<D,R>;
47  static const std::size_t Variants = 8;
48  };
49 
50  template<int dim, typename D, typename R, std::size_t k>
51  struct BDMCubeLocalInfo
52  {
53  static_assert((AlwaysFalse<D>::value),"The requested type of BDM element is not implemented, sorry!");
54  };
55 
56  template<typename D, typename R>
57  struct BDMCubeLocalInfo<2,D,R,1>
58  {
59  using FiniteElement = BDM1Cube2DLocalFiniteElement<D,R>;
60  static const std::size_t Variants = 16;
61  };
62 
63  template<typename D, typename R>
64  struct BDMCubeLocalInfo<2,D,R,2>
65  {
66  using FiniteElement = BDM2Cube2DLocalFiniteElement<D,R>;
67  static const std::size_t Variants = 16;
68  };
69 
70  template<typename D, typename R>
71  struct BDMCubeLocalInfo<3,D,R,1>
72  {
73  using FiniteElement = BDM1Cube3DLocalFiniteElement<D,R>;
74  static const std::size_t Variants = 64;
75  };
76 
77  template<typename GV, int dim, typename R, std::size_t k>
78  class BDMLocalFiniteElementMap
79  {
80  using D = typename GV::ctype;
81  using CubeFiniteElement = typename BDMCubeLocalInfo<dim, D, R, k>::FiniteElement;
82  using SimplexFiniteElement = typename BDMSimplexLocalInfo<dim, D, R, k>::FiniteElement;
83 
84  public:
85 
86  using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
87  using FiniteElement = LocalFiniteElementVirtualInterface<T>;
88 
89  BDMLocalFiniteElementMap(const GV& gv)
90  : is_(&(gv.indexSet())), orient_(gv.size(0))
91  {
92  cubeVariant_.resize(BDMCubeLocalInfo<dim, D, R, k>::Variants);
93  simplexVariant_.resize(BDMSimplexLocalInfo<dim, D, R, k>::Variants);
94 
95  // create all variants
96  for (size_t i = 0; i < cubeVariant_.size(); i++)
97  cubeVariant_[i] = std::make_unique<LocalFiniteElementVirtualImp<CubeFiniteElement> >(CubeFiniteElement(i));
98 
99  for (size_t i = 0; i < simplexVariant_.size(); i++)
100  simplexVariant_[i] = std::make_unique<LocalFiniteElementVirtualImp<SimplexFiniteElement> >(SimplexFiniteElement(i));
101 
102  // compute orientation for all elements
103  // loop once over the grid
104  for(const auto& cell : elements(gv))
105  {
106  unsigned int myId = is_->index(cell);
107  orient_[myId] = 0;
108 
109  for (const auto& intersection : intersections(gv,cell))
110  {
111  if (intersection.neighbor() && (is_->index(intersection.outside()) > myId))
112  orient_[myId] |= (1 << intersection.indexInInside());
113  }
114  }
115  }
116 
118  template<class EntityType>
119  const FiniteElement& find(const EntityType& e) const
120  {
121  if (e.type().isCube())
122  return *cubeVariant_[orient_[is_->index(e)]];
123  else
124  return *simplexVariant_[orient_[is_->index(e)]];
125  }
126 
127  private:
128  std::vector<std::unique_ptr<LocalFiniteElementVirtualImp<CubeFiniteElement> > > cubeVariant_;
129  std::vector<std::unique_ptr<LocalFiniteElementVirtualImp<SimplexFiniteElement> > > simplexVariant_;
130  const typename GV::IndexSet* is_;
131  std::vector<unsigned char> orient_;
132  };
133 
134 
135 } // namespace Impl
136 
137 
138 // *****************************************************************************
139 // This is the reusable part of the basis. It contains
140 //
141 // BrezziDouglasMariniPreBasis
142 // BrezziDouglasMariniNodeIndexSet
143 // BrezziDouglasMariniNode
144 //
145 // The pre-basis allows to create the others and is the owner of possible shared
146 // state. These three components do _not_ depend on the global basis or index
147 // set and can be used without a global basis.
148 // *****************************************************************************
149 
150 template<typename GV, int k>
151 class BrezziDouglasMariniNode;
152 
153 template<typename GV, int k, class MI>
154 class BrezziDouglasMariniNodeIndexSet;
155 
156 template<typename GV, int k, class MI>
158 {
159  static const int dim = GV::dimension;
160  using FiniteElementMap = typename Impl::BDMLocalFiniteElementMap<GV, dim, double, k>;
161 
162  template<typename, int, class>
164 
165 public:
166 
168  using GridView = GV;
169  using size_type = std::size_t;
170 
172 
174 
176  using MultiIndex = MI;
177 
178  using SizePrefix = Dune::ReservedVector<size_type, 1>;
179 
182  gridView_(gv),
184  {
185  // There is no inherent reason why the basis shouldn't work for grids with more than one
186  // element types. Somebody simply has to sit down and implement the missing bits.
187  if (gv.indexSet().types(0).size() > 1)
188  DUNE_THROW(Dune::NotImplemented, "Brezzi-Douglas-Marini basis is only implemented for grids with a single element type");
189  }
190 
192  {
193  codimOffset_[0] = 0;
194  codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
195  //if (dim==3) codimOffset_[2] = codimOffset_[1] + dofsPerCodim[1] * gridView_.size(1);
196  }
197 
200  const GridView& gridView() const
201  {
202  return gridView_;
203  }
204 
205  /* \brief Update the stored grid view, to be called if the grid has changed */
206  void update (const GridView& gv)
207  {
208  gridView_ = gv;
209  }
210 
214  Node makeNode() const
215  {
216  return Node{&finiteElementMap_};
217  }
218 
226  {
227  return IndexSet{*this};
228  }
229 
230  size_type size() const
231  {
232  return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1); // only 2d
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  // The implementation currently only supports grids with a single element type.
251  // We can therefore return the actual number of dofs here.
252  GeometryType elementType = *(gridView_.indexSet().types(0).begin());
253  size_t numFaces = ReferenceElements<double,dim>::general(elementType).size(1);
254  return dofsPerCodim_[0] + dofsPerCodim_[1] * numFaces;
255  }
256 
257 protected:
259  std::array<size_t,dim+1> codimOffset_;
260  FiniteElementMap finiteElementMap_;
261  // Number of dofs per entity type depending on the entity's codimension and type
262  std::array<int,2> dofsPerCodim_ {{dim*(k-1)*3, dim+(k-1)}};
263 };
264 
265 
266 
267 template<typename GV, int k>
269  public LeafBasisNode
270 {
271  static const int dim = GV::dimension;
272 
273 public:
274 
275  using size_type = std::size_t;
276  using Element = typename GV::template Codim<0>::Entity;
277  using FiniteElementMap = typename Impl::BDMLocalFiniteElementMap<GV, dim, double, k>;
278  using FiniteElement = typename FiniteElementMap::FiniteElement;
279 
280  BrezziDouglasMariniNode(const FiniteElementMap* finiteElementMap) :
281  finiteElement_(nullptr),
282  element_(nullptr),
283  finiteElementMap_(finiteElementMap)
284  {}
285 
287  const Element& element() const
288  {
289  return *element_;
290  }
291 
297  {
298  return *finiteElement_;
299  }
300 
302  void bind(const Element& e)
303  {
304  element_ = &e;
306  this->setSize(finiteElement_->size());
307  }
308 
309 protected:
310 
314 };
315 
316 
317 
318 template<typename GV, int k, class MI>
320 {
321  enum {dim = GV::dimension};
322 
323 public:
324 
325  using size_type = std::size_t;
326 
328  using MultiIndex = MI;
329 
331 
333 
335  preBasis_(&preBasis)
336  {}
337 
343  void bind(const Node& node)
344  {
345  node_ = &node;
346  }
347 
350  void unbind()
351  {
352  node_ = nullptr;
353  }
354 
357  size_type size() const
358  {
359  return node_->finiteElement().size();
360  }
361 
367  template<typename It>
368  It indices(It it) const
369  {
370  const auto& gridIndexSet = preBasis_->gridView().indexSet();
371  const auto& element = node_->element();
372 
373  // throw if element is not of predefined type
374  if (not(element.type().isCube()) and not(element.type().isSimplex()))
375  DUNE_THROW(Dune::NotImplemented, "BrezziDouglasMariniBasis only implemented for cube and simplex elements.");
376 
377  for(std::size_t i=0, end=size(); i<end; ++i, ++it)
378  {
379  Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
380 
381  // The dimension of the entity that the current dof is related to
382  size_t subentity = localKey.subEntity();
383  size_t codim = localKey.codim();
384 
385  *it = { preBasis_->codimOffset_[codim] +
386  preBasis_->dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
387  }
388 
389  return it;
390  }
391 
392 protected:
394  const Node* node_;
395 };
396 
397 
398 
399 namespace BasisFactory {
400 
401 namespace Imp {
402 
403 template<std::size_t k>
404 class BrezziDouglasMariniPreBasisFactory
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
412  {
413  return {gridView};
414  }
415 };
416 
417 } // end namespace BasisFactory::Imp
418 
426 template<std::size_t k>
428 {
429  return Imp::BrezziDouglasMariniPreBasisFactory<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_BREZZIDOUGLASMARINIBASIS_HH
auto brezziDouglasMarini()
Create a pre-basis factory that can create a Brezzi-Douglas-Marini pre-basis.
Definition: brezzidouglasmarinibasis.hh:427
Definition: polynomial.hh:10
Definition: brezzidouglasmarinibasis.hh:270
typename FiniteElementMap::FiniteElement FiniteElement
Definition: brezzidouglasmarinibasis.hh:278
typename GV::template Codim< 0 >::Entity Element
Definition: brezzidouglasmarinibasis.hh:276
const FiniteElementMap * finiteElementMap_
Definition: brezzidouglasmarinibasis.hh:313
const Element & element() const
Return current element, throw if unbound.
Definition: brezzidouglasmarinibasis.hh:287
std::size_t size_type
Definition: brezzidouglasmarinibasis.hh:275
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition: brezzidouglasmarinibasis.hh:296
typename Impl::BDMLocalFiniteElementMap< GV, dim, double, k > FiniteElementMap
Definition: brezzidouglasmarinibasis.hh:277
const Element * element_
Definition: brezzidouglasmarinibasis.hh:312
void bind(const Element &e)
Bind to element.
Definition: brezzidouglasmarinibasis.hh:302
BrezziDouglasMariniNode(const FiniteElementMap *finiteElementMap)
Definition: brezzidouglasmarinibasis.hh:280
const FiniteElement * finiteElement_
Definition: brezzidouglasmarinibasis.hh:311
Definition: brezzidouglasmarinibasis.hh:320
const Node * node_
Definition: brezzidouglasmarinibasis.hh:394
size_type size() const
Size of subtree rooted in this node (element-local)
Definition: brezzidouglasmarinibasis.hh:357
void bind(const Node &node)
Bind the view to a grid element.
Definition: brezzidouglasmarinibasis.hh:343
const PreBasis * preBasis_
Definition: brezzidouglasmarinibasis.hh:393
void unbind()
Unbind the view.
Definition: brezzidouglasmarinibasis.hh:350
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: brezzidouglasmarinibasis.hh:328
It indices(It it) const
Maps from subtree index set [0..size-1] to a globally unique multi index in global basis.
Definition: brezzidouglasmarinibasis.hh:368
std::size_t size_type
Definition: brezzidouglasmarinibasis.hh:325
BrezziDouglasMariniNodeIndexSet(const PreBasis &preBasis)
Definition: brezzidouglasmarinibasis.hh:334
Definition: brezzidouglasmarinibasis.hh:158
BrezziDouglasMariniPreBasis(const GridView &gv)
Constructor for a given grid view object.
Definition: brezzidouglasmarinibasis.hh:181
size_type maxNodeSize() const
Definition: brezzidouglasmarinibasis.hh:248
Node makeNode() const
Create tree node.
Definition: brezzidouglasmarinibasis.hh:214
std::array< int, 2 > dofsPerCodim_
Definition: brezzidouglasmarinibasis.hh:262
void update(const GridView &gv)
Definition: brezzidouglasmarinibasis.hh:206
Dune::ReservedVector< size_type, 1 > SizePrefix
Definition: brezzidouglasmarinibasis.hh:178
FiniteElementMap finiteElementMap_
Definition: brezzidouglasmarinibasis.hh:260
size_type size() const
Definition: brezzidouglasmarinibasis.hh:230
IndexSet makeIndexSet() const
Create tree node index set.
Definition: brezzidouglasmarinibasis.hh:225
const GridView gridView_
Definition: brezzidouglasmarinibasis.hh:258
GV GridView
The grid view that the FE space is defined on.
Definition: brezzidouglasmarinibasis.hh:168
void initializeIndices()
Definition: brezzidouglasmarinibasis.hh:191
std::array< size_t, dim+1 > codimOffset_
Definition: brezzidouglasmarinibasis.hh:259
size_type size(const SizePrefix prefix) const
Return number possible values for next position in multi index.
Definition: brezzidouglasmarinibasis.hh:236
std::size_t size_type
Definition: brezzidouglasmarinibasis.hh:169
size_type dimension() const
Definition: brezzidouglasmarinibasis.hh:243
MI MultiIndex
Type used for global numbering of the basis vectors.
Definition: brezzidouglasmarinibasis.hh:176
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition: brezzidouglasmarinibasis.hh:200
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:47
void setSize(const size_type size)
Definition: nodes.hh:160
Definition: nodes.hh:182