3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINIBASIS_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINIBASIS_HH
8#include <dune/geometry/referenceelements.hh>
10#include <dune/localfunctions/common/virtualinterface.hh>
11#include <dune/localfunctions/common/virtualwrappers.hh>
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>
19#include <dune/typetree/leafnode.hh>
21#include <dune/functions/functionspacebases/nodes.hh>
22#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
23#include <dune/functions/functionspacebases/flatmultiindex.hh>
30 template<
int dim,
typename D,
typename R, std::
size_t k>
31 struct BDMSimplexLocalInfo
33 static_assert((
AlwaysFalse<D>::value),
"The requested type of BDM element is not implemented, sorry!");
36 template<
typename D,
typename R>
37 struct BDMSimplexLocalInfo<2,D,R,1>
39 using FiniteElement = BDM1Simplex2DLocalFiniteElement<D,R>;
40 static const std::size_t Variants = 8;
43 template<
typename D,
typename R>
44 struct BDMSimplexLocalInfo<2,D,R,2>
46 using FiniteElement = BDM2Simplex2DLocalFiniteElement<D,R>;
47 static const std::size_t Variants = 8;
50 template<
int dim,
typename D,
typename R, std::
size_t k>
51 struct BDMCubeLocalInfo
53 static_assert((
AlwaysFalse<D>::value),
"The requested type of BDM element is not implemented, sorry!");
56 template<
typename D,
typename R>
57 struct BDMCubeLocalInfo<2,D,R,1>
59 using FiniteElement = BDM1Cube2DLocalFiniteElement<D,R>;
60 static const std::size_t Variants = 16;
63 template<
typename D,
typename R>
64 struct BDMCubeLocalInfo<2,D,R,2>
66 using FiniteElement = BDM2Cube2DLocalFiniteElement<D,R>;
67 static const std::size_t Variants = 16;
70 template<
typename D,
typename R>
71 struct BDMCubeLocalInfo<3,D,R,1>
73 using FiniteElement = BDM1Cube3DLocalFiniteElement<D,R>;
74 static const std::size_t Variants = 64;
77 template<
typename GV,
int dim,
typename R, std::
size_t k>
78 class BDMLocalFiniteElementMap
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;
86 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
87 using FiniteElement = LocalFiniteElementVirtualInterface<T>;
89 BDMLocalFiniteElementMap(
const GV& gv)
90 : is_(&(gv.indexSet())), orient_(gv.size(0))
92 cubeVariant_.resize(BDMCubeLocalInfo<dim, D, R, k>::Variants);
93 simplexVariant_.resize(BDMSimplexLocalInfo<dim, D, R, k>::Variants);
96 for (
size_t i = 0; i < cubeVariant_.size(); i++)
97 cubeVariant_[i] = std::make_unique<LocalFiniteElementVirtualImp<CubeFiniteElement> >(CubeFiniteElement(i));
99 for (
size_t i = 0; i < simplexVariant_.size(); i++)
100 simplexVariant_[i] = std::make_unique<LocalFiniteElementVirtualImp<SimplexFiniteElement> >(SimplexFiniteElement(i));
104 for(
const auto& cell : elements(gv))
106 unsigned int myId = is_->index(cell);
109 for (
const auto& intersection : intersections(gv,cell))
111 if (intersection.neighbor() && (is_->index(intersection.outside()) > myId))
112 orient_[myId] |= (1 << intersection.indexInInside());
118 template<
class EntityType>
119 const FiniteElement& find(
const EntityType& e)
const
121 if (e.type().isCube())
122 return *cubeVariant_[orient_[is_->index(e)]];
124 return *simplexVariant_[orient_[is_->index(e)]];
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_;
150template<
typename GV,
int k>
151class BrezziDouglasMariniNode;
153template<
typename GV,
int k,
class MI>
154class BrezziDouglasMariniNodeIndexSet;
156template<
typename GV,
int k,
class MI>
157class BrezziDouglasMariniPreBasis
159 static const int dim = GV::dimension;
160 using FiniteElementMap =
typename Impl::BDMLocalFiniteElementMap<GV, dim, double, k>;
162 template<
typename,
int,
class>
163 friend class BrezziDouglasMariniNodeIndexSet;
169 using size_type = std::size_t;
171 using Node = BrezziDouglasMariniNode<GV, k>;
173 using IndexSet = BrezziDouglasMariniNodeIndexSet<GV, k, MI>;
176 using MultiIndex = MI;
181 BrezziDouglasMariniPreBasis(
const GridView& gv) :
183 finiteElementMap_(gv)
187 if (gv.indexSet().types(0).size() > 1)
191 void initializeIndices()
194 codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
200 const GridView& gridView()
const
206 void update (
const GridView& gv)
214 Node makeNode()
const
216 return Node{&finiteElementMap_};
225 IndexSet makeIndexSet()
const
227 return IndexSet{*
this};
230 size_type size()
const
232 return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1);
236 size_type size(
const SizePrefix prefix)
const
238 assert(prefix.size() == 0 || prefix.size() == 1);
239 return (prefix.size() == 0) ? size() : 0;
243 size_type dimension()
const
248 size_type maxNodeSize()
const
252 GeometryType elementType = *(gridView_.indexSet().types(0).begin());
254 return dofsPerCodim_[0] + dofsPerCodim_[1] * numFaces;
258 const GridView gridView_;
259 std::array<size_t,dim+1> codimOffset_;
260 FiniteElementMap finiteElementMap_;
262 std::array<int,2> dofsPerCodim_ {{dim*(k-1)*3, dim+(k-1)}};
267template<
typename GV,
int k>
268class BrezziDouglasMariniNode :
271 static const int dim = GV::dimension;
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;
280 BrezziDouglasMariniNode(
const FiniteElementMap* finiteElementMap) :
281 finiteElement_(nullptr),
283 finiteElementMap_(finiteElementMap)
287 const Element& element()
const
296 const FiniteElement& finiteElement()
const
298 return *finiteElement_;
302 void bind(
const Element& e)
305 finiteElement_ = &(finiteElementMap_->find(*element_));
306 this->setSize(finiteElement_->size());
311 const FiniteElement* finiteElement_;
312 const Element* element_;
313 const FiniteElementMap* finiteElementMap_;
318template<
typename GV,
int k,
class MI>
319class BrezziDouglasMariniNodeIndexSet
321 enum {dim = GV::dimension};
325 using size_type = std::size_t;
328 using MultiIndex = MI;
330 using PreBasis = BrezziDouglasMariniPreBasis<GV, k, MI>;
332 using Node = BrezziDouglasMariniNode<GV, k>;
334 BrezziDouglasMariniNodeIndexSet(
const PreBasis& preBasis) :
343 void bind(
const Node& node)
357 size_type size()
const
359 return node_->finiteElement().size();
367 template<
typename It>
368 It indices(It it)
const
370 const auto& gridIndexSet = preBasis_->gridView().indexSet();
371 const auto& element = node_->element();
374 if (not(element.type().isCube()) and not(element.type().isSimplex()))
377 for(std::size_t i=0, end=size(); i<end; ++i, ++it)
379 Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
383 size_t codim = localKey.
codim();
385 *it = { preBasis_->codimOffset_[codim] +
386 preBasis_->dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.
index() };
393 const PreBasis* preBasis_;
399namespace BasisFactory {
403template<std::
size_t k>
404class BrezziDouglasMariniPreBasisFactory
407 static const std::size_t requiredMultiIndexSize=1;
409 template<
class MultiIndex,
class Gr
idView>
410 auto makePreBasis(
const GridView& gridView)
const
411 -> BrezziDouglasMariniPreBasis<GridView, k, MultiIndex>
426template<std::
size_t k>
429 return Imp::BrezziDouglasMariniPreBasisFactory<k>();
447template<
typename GV,
int k>
448using BrezziDouglasMariniBasis = DefaultGlobalBasis<BrezziDouglasMariniPreBasis<GV, k, FlatMultiIndex<std::size_t>> >;
Describe position of one degree of freedom.
Definition: localkey.hh:21
unsigned int index() const
Return offset within subentity.
Definition: localkey.hh:66
unsigned int codim() const
Return codim of associated entity.
Definition: localkey.hh:60
unsigned int subEntity() const
Return number of associated subentity.
Definition: localkey.hh:54
Default exception for dummy implementations.
Definition: exceptions.hh:261
A Vector class with statically reserved memory.
Definition: reservedvector.hh:43
A few common exception classes.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:180
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
auto brezziDouglasMarini()
Create a pre-basis factory that can create a Brezzi-Douglas-Marini pre-basis.
Definition: brezzidouglasmarinibasis.hh:427
Dune namespace.
Definition: alignedallocator.hh:14
static const bool value
always a false value
Definition: typetraits.hh:128
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196