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>
23#include <dune/functions/functionspacebases/globalvaluedlocalfiniteelement.hh>
24#include <dune/functions/functionspacebases/nodes.hh>
25#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
26#include <dune/functions/functionspacebases/leafprebasismixin.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>
156class BrezziDouglasMariniPreBasis :
157 public LeafPreBasisMixin< BrezziDouglasMariniPreBasis<GV,k> >
159 static const int dim = GV::dimension;
160 using FiniteElementMap =
typename Impl::BDMLocalFiniteElementMap<GV, dim, double, k>;
166 using size_type = std::size_t;
168 using Node = BrezziDouglasMariniNode<GV, k>;
171 BrezziDouglasMariniPreBasis(
const GridView& gv) :
173 finiteElementMap_(gv)
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");
181 void initializeIndices()
184 codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
190 const GridView& gridView()
const
196 void update (
const GridView& gv)
204 Node makeNode()
const
206 return Node{&finiteElementMap_};
209 size_type dimension()
const
211 return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1);
214 size_type maxNodeSize()
const
218 GeometryType elementType = *(gridView_.indexSet().types(0).begin());
219 size_t numFaces = ReferenceElements<double,dim>::general(elementType).size(1);
220 return dofsPerCodim_[0] + dofsPerCodim_[1] * numFaces;
228 template<
typename It>
229 It indices(
const Node& node, It it)
const
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();
246 *it = { codimOffset_[codim] +
247 dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
255 std::array<size_t,dim+1> codimOffset_;
256 FiniteElementMap finiteElementMap_;
258 std::array<int,2> dofsPerCodim_ {{dim*(k-1)*3, dim+(k-1)}};
263template<
typename GV,
int k>
264class BrezziDouglasMariniNode :
267 static const int dim = GV::dimension;
271 using size_type = std::size_t;
272 using Element =
typename GV::template Codim<0>::Entity;
273 using FiniteElementMap =
typename Impl::BDMLocalFiniteElementMap<GV, dim, double, k>;
274 using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::ContravariantPiolaTransformator,
275 typename FiniteElementMap::FiniteElement,
278 BrezziDouglasMariniNode(
const FiniteElementMap* finiteElementMap) :
280 finiteElementMap_(finiteElementMap)
284 const Element& element()
const
293 const FiniteElement& finiteElement()
const
295 return finiteElement_;
299 void bind(
const Element& e)
302 finiteElement_.bind((finiteElementMap_->find(*element_)), e);
303 this->setSize(finiteElement_.size());
308 FiniteElement finiteElement_;
309 const Element* element_;
310 const FiniteElementMap* finiteElementMap_;
315namespace BasisFactory {
324template<std::
size_t k>
327 return [](
const auto& gridView) {
328 return BrezziDouglasMariniPreBasis<std::decay_t<
decltype(gridView)>, k>(gridView);
347template<
typename GV,
int k>
348using BrezziDouglasMariniBasis = DefaultGlobalBasis<BrezziDouglasMariniPreBasis<GV, 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