7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINIBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINIBASIS_HH
11#include <dune/common/exceptions.hh>
12#include <dune/geometry/referenceelements.hh>
14#include <dune/localfunctions/common/virtualinterface.hh>
15#include <dune/localfunctions/common/virtualwrappers.hh>
17#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1cube2d.hh>
18#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1cube3d.hh>
19#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1simplex2d.hh>
20#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini2cube2d.hh>
21#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini2simplex2d.hh>
33 template<
int dim,
typename D,
typename R, std::
size_t k>
34 struct BDMSimplexLocalInfo
36 static_assert((AlwaysFalse<D>::value),
"The requested type of BDM element is not implemented, sorry!");
39 template<
typename D,
typename R>
40 struct BDMSimplexLocalInfo<2,D,R,1>
42 using FiniteElement = BDM1Simplex2DLocalFiniteElement<D,R>;
43 static const std::size_t Variants = 8;
46 template<
typename D,
typename R>
47 struct BDMSimplexLocalInfo<2,D,R,2>
49 using FiniteElement = BDM2Simplex2DLocalFiniteElement<D,R>;
50 static const std::size_t Variants = 8;
53 template<
int dim,
typename D,
typename R, std::
size_t k>
54 struct BDMCubeLocalInfo
56 static_assert((AlwaysFalse<D>::value),
"The requested type of BDM element is not implemented, sorry!");
59 template<
typename D,
typename R>
60 struct BDMCubeLocalInfo<2,D,R,1>
62 using FiniteElement = BDM1Cube2DLocalFiniteElement<D,R>;
63 static const std::size_t Variants = 16;
66 template<
typename D,
typename R>
67 struct BDMCubeLocalInfo<2,D,R,2>
69 using FiniteElement = BDM2Cube2DLocalFiniteElement<D,R>;
70 static const std::size_t Variants = 16;
73 template<
typename D,
typename R>
74 struct BDMCubeLocalInfo<3,D,R,1>
76 using FiniteElement = BDM1Cube3DLocalFiniteElement<D,R>;
77 static const std::size_t Variants = 64;
80 template<
typename GV,
int dim,
typename R, std::
size_t k>
81 class BDMLocalFiniteElementMap
83 using D =
typename GV::ctype;
84 using CubeFiniteElement =
typename BDMCubeLocalInfo<dim, D, R, k>::FiniteElement;
85 using SimplexFiniteElement =
typename BDMSimplexLocalInfo<dim, D, R, k>::FiniteElement;
89 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
90 using FiniteElement = LocalFiniteElementVirtualInterface<T>;
92 BDMLocalFiniteElementMap(
const GV& gv)
93 : is_(&(gv.indexSet())), orient_(gv.size(0))
95 cubeVariant_.resize(BDMCubeLocalInfo<dim, D, R, k>::Variants);
96 simplexVariant_.resize(BDMSimplexLocalInfo<dim, D, R, k>::Variants);
99 for (
size_t i = 0; i < cubeVariant_.size(); i++)
100 cubeVariant_[i] = std::make_shared<LocalFiniteElementVirtualImp<CubeFiniteElement> >(CubeFiniteElement(i));
102 for (
size_t i = 0; i < simplexVariant_.size(); i++)
103 simplexVariant_[i] = std::make_shared<LocalFiniteElementVirtualImp<SimplexFiniteElement> >(SimplexFiniteElement(i));
107 for(
const auto& cell : elements(gv))
109 unsigned int myId = is_->index(cell);
112 for (
const auto& intersection : intersections(gv,cell))
114 if (intersection.neighbor() && (is_->index(intersection.outside()) > myId))
115 orient_[myId] |= (1 << intersection.indexInInside());
121 template<
class EntityType>
122 const FiniteElement& find(
const EntityType& e)
const
124 if (e.type().isCube())
125 return *cubeVariant_[orient_[is_->index(e)]];
127 return *simplexVariant_[orient_[is_->index(e)]];
131 std::vector<std::shared_ptr<LocalFiniteElementVirtualImp<CubeFiniteElement> > > cubeVariant_;
132 std::vector<std::shared_ptr<LocalFiniteElementVirtualImp<SimplexFiniteElement> > > simplexVariant_;
133 const typename GV::IndexSet* is_;
134 std::vector<unsigned char> orient_;
152template<
typename GV,
int k>
153class BrezziDouglasMariniNode;
155template<
typename GV,
int k>
159 static const int dim = GV::dimension;
160 using FiniteElementMap =
typename Impl::BDMLocalFiniteElementMap<GV, dim, double, k>;
177 if (gv.indexSet().types(0).size() > 1)
178 DUNE_THROW(Dune::NotImplemented,
"Brezzi-Douglas-Marini basis is only implemented for grids with a single element type");
218 GeometryType elementType = *(
gridView_.indexSet().types(0).begin());
219 size_t numFaces = ReferenceElements<double,dim>::general(elementType).size(1);
228 template<
typename It>
231 const auto& gridIndexSet =
gridView().indexSet();
232 const auto& element = node.
element();
235 if (not(element.type().isCube()) and not(element.type().isSimplex()))
236 DUNE_THROW(Dune::NotImplemented,
"BrezziDouglasMariniBasis only implemented for cube and simplex elements.");
238 for(std::size_t i=0, end=node.
size(); i<end; ++i, ++it)
240 Dune::LocalKey localKey = node.
finiteElement().localCoefficients().localKey(i);
243 size_t subentity = localKey.subEntity();
244 size_t codim = localKey.codim();
247 dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
263template<
typename GV,
int k>
267 static const int dim = GV::dimension;
272 using Element =
typename GV::template Codim<0>::Entity;
274 using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::ContravariantPiolaTransformator,
275 typename FiniteElementMap::FiniteElement,
315namespace BasisFactory {
324template<std::
size_t k>
327 return [](
const auto& gridView) {
347template<
typename GV,
int k>
auto brezziDouglasMarini()
Create a pre-basis factory that can create a Brezzi-Douglas-Marini pre-basis.
Definition brezzidouglasmarinibasis.hh:325
Definition polynomial.hh:17
Definition brezzidouglasmarinibasis.hh:266
const FiniteElement & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition brezzidouglasmarinibasis.hh:293
typename GV::template Codim< 0 >::Entity Element
Definition brezzidouglasmarinibasis.hh:272
const FiniteElementMap * finiteElementMap_
Definition brezzidouglasmarinibasis.hh:310
std::size_t size_type
Definition brezzidouglasmarinibasis.hh:271
FiniteElement finiteElement_
Definition brezzidouglasmarinibasis.hh:308
typename Impl::BDMLocalFiniteElementMap< GV, dim, double, k > FiniteElementMap
Definition brezzidouglasmarinibasis.hh:273
const Element * element_
Definition brezzidouglasmarinibasis.hh:309
Impl::GlobalValuedLocalFiniteElement< Impl::ContravariantPiolaTransformator, typename FiniteElementMap::FiniteElement, Element > FiniteElement
Definition brezzidouglasmarinibasis.hh:276
void bind(const Element &e)
Bind to element.
Definition brezzidouglasmarinibasis.hh:299
BrezziDouglasMariniNode(const FiniteElementMap *finiteElementMap)
Definition brezzidouglasmarinibasis.hh:278
const Element & element() const
Return current element, throw if unbound.
Definition brezzidouglasmarinibasis.hh:284
Definition brezzidouglasmarinibasis.hh:158
std::size_t size_type
Definition brezzidouglasmarinibasis.hh:166
std::array< int, 2 > dofsPerCodim_
Definition brezzidouglasmarinibasis.hh:258
size_type dimension() const
Definition brezzidouglasmarinibasis.hh:209
std::array< size_t, dim+1 > codimOffset_
Definition brezzidouglasmarinibasis.hh:255
Node makeNode() const
Create tree node.
Definition brezzidouglasmarinibasis.hh:204
GV GridView
The grid view that the FE space is defined on.
Definition brezzidouglasmarinibasis.hh:165
BrezziDouglasMariniPreBasis(const GridView &gv)
Constructor for a given grid view object.
Definition brezzidouglasmarinibasis.hh:171
FiniteElementMap finiteElementMap_
Definition brezzidouglasmarinibasis.hh:256
void update(const GridView &gv)
Definition brezzidouglasmarinibasis.hh:196
void initializeIndices()
Definition brezzidouglasmarinibasis.hh:181
const GridView & gridView() const
Obtain the grid view that the basis is defined on.
Definition brezzidouglasmarinibasis.hh:190
size_type maxNodeSize() const
Definition brezzidouglasmarinibasis.hh:214
GridView gridView_
Definition brezzidouglasmarinibasis.hh:254
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 brezzidouglasmarinibasis.hh:229
Global basis for given pre-basis.
Definition defaultglobalbasis.hh:50
A generic MixIn class for PreBasis.
Definition leafprebasismixin.hh:36
size_type size() const
Definition nodes.hh:147
void setSize(const size_type size)
Definition nodes.hh:169