7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
11#include <dune/common/exceptions.hh>
13#include <dune/grid/common/capabilities.hh>
14#include <dune/grid/common/mcmgmapper.hh>
16#include <dune/localfunctions/common/localfiniteelementvariant.hh>
17#include <dune/localfunctions/raviartthomas.hh>
18#include <dune/localfunctions/raviartthomas/raviartthomas0cube2d.hh>
19#include <dune/localfunctions/raviartthomas/raviartthomas0cube3d.hh>
20#include <dune/localfunctions/raviartthomas/raviartthomas02d.hh>
21#include <dune/localfunctions/raviartthomas/raviartthomas03d.hh>
22#include <dune/localfunctions/raviartthomas/raviartthomas1cube2d.hh>
23#include <dune/localfunctions/raviartthomas/raviartthomas1cube3d.hh>
24#include <dune/localfunctions/raviartthomas/raviartthomas12d.hh>
25#include <dune/localfunctions/raviartthomas/raviartthomas2cube2d.hh>
27#include <dune/functions/functionspacebases/globalvaluedlocalfiniteelement.hh>
28#include <dune/functions/functionspacebases/nodes.hh>
29#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
30#include <dune/functions/functionspacebases/leafprebasismixin.hh>
37 template<
int dim,
typename D,
typename R, std::
size_t k>
38 struct RaviartThomasSimplexLocalInfo
41 using FiniteElement =
void*;
44 template<
typename D,
typename R>
45 struct RaviartThomasSimplexLocalInfo<2,D,R,0>
47 using FiniteElement = RT02DLocalFiniteElement<D,R>;
50 template<
typename D,
typename R>
51 struct RaviartThomasSimplexLocalInfo<2,D,R,1>
53 using FiniteElement = RT12DLocalFiniteElement<D,R>;
56 template<
typename D,
typename R>
57 struct RaviartThomasSimplexLocalInfo<3,D,R,0>
59 using FiniteElement = RT03DLocalFiniteElement<D,R>;
62 template<
int dim,
typename D,
typename R, std::
size_t k>
63 struct RaviartThomasCubeLocalInfo
66 using FiniteElement =
void*;
69 template<
typename D,
typename R>
70 struct RaviartThomasCubeLocalInfo<2,D,R,0>
72 using FiniteElement = RT0Cube2DLocalFiniteElement<D,R>;
75 template<
typename D,
typename R>
76 struct RaviartThomasCubeLocalInfo<2,D,R,1>
78 using FiniteElement = RT1Cube2DLocalFiniteElement<D,R>;
81 template<
typename D,
typename R>
82 struct RaviartThomasCubeLocalInfo<2,D,R,2>
84 using FiniteElement = RT2Cube2DLocalFiniteElement<D,R>;
87 template<
typename D,
typename R>
88 struct RaviartThomasCubeLocalInfo<3,D,R,0>
90 using FiniteElement = RT0Cube3DLocalFiniteElement<D,R>;
93 template<
typename D,
typename R>
94 struct RaviartThomasCubeLocalInfo<3,D,R,1>
96 using FiniteElement = RT1Cube3DLocalFiniteElement<D,R>;
99 template<
typename GV,
int dim,
typename R, std::
size_t k>
100 class RaviartThomasLocalFiniteElementMap
102 using D =
typename GV::ctype;
103 constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
105 using CubeFiniteElement =
typename RaviartThomasCubeLocalInfo<dim, D, R, k>::FiniteElement;
106 using SimplexFiniteElement =
typename RaviartThomasSimplexLocalInfo<dim, D, R, k>::FiniteElement;
110 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
112 constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId;
113 constexpr static GeometryType type = GeometryType(topologyId, GV::dimension);
115 using FiniteElement = std::conditional_t<hasFixedElementType,
116 std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
117 LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
121 static std::size_t numVariants(GeometryType type)
123 auto numFacets = referenceElement<D,dim>(type).size(1);
124 return power(2,numFacets);
127 RaviartThomasLocalFiniteElementMap(
const GV& gv)
128 : elementMapper_(gv, mcmgElementLayout()),
131 if constexpr (hasFixedElementType)
133 variants_.resize(numVariants(type));
134 for (
size_t i = 0; i < numVariants(type); i++)
135 variants_[i] = FiniteElement(i);
140 variants_.resize(numVariants(GeometryTypes::simplex(dim)) + numVariants(GeometryTypes::cube(dim)));
141 for (
size_t i = 0; i < numVariants(GeometryTypes::simplex(dim)); i++)
142 variants_[i] = SimplexFiniteElement(i);
143 for (
size_t i = 0; i < numVariants(GeometryTypes::cube(dim)); i++)
144 variants_[i + numVariants(GeometryTypes::simplex(dim))] = CubeFiniteElement(i);
147 for(
const auto& cell : elements(gv))
149 unsigned int myId = elementMapper_.index(cell);
152 for (
const auto& intersection : intersections(gv,cell))
154 if (intersection.neighbor() && (elementMapper_.index(intersection.outside()) > myId))
155 orient_[myId] |= (1 << intersection.indexInInside());
159 if constexpr (!hasFixedElementType)
160 if (cell.type().isCube())
161 orient_[myId] += numVariants(GeometryTypes::simplex(dim));
165 template<
class EntityType>
166 const FiniteElement& find(
const EntityType& e)
const
168 return variants_[orient_[elementMapper_.index(e)]];
172 std::vector<FiniteElement> variants_;
173 Dune::MultipleCodimMultipleGeomTypeMapper<GV> elementMapper_;
174 std::vector<unsigned char> orient_;
192template<
typename GV,
int k>
193class RaviartThomasNode;
195template<
typename GV,
int k>
196class RaviartThomasPreBasis :
197 public LeafPreBasisMixin< RaviartThomasPreBasis<GV,k> >
199 static const int dim = GV::dimension;
200 using FiniteElementMap =
typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
206 using size_type = std::size_t;
208 using Node = RaviartThomasNode<GV, k>;
211 RaviartThomasPreBasis(
const GridView& gv) :
213 finiteElementMap_(gv)
216 if (gv.indexSet().types(0).size() > 1 and k>0)
217 DUNE_THROW(Dune::NotImplemented,
"Raviart-Thomas basis with index k>0 is only implemented for grids with a single element type");
219 for(
auto type : gv.indexSet().types(0))
220 if (!type.isSimplex() && !type.isCube())
221 DUNE_THROW(Dune::NotImplemented,
"Raviart-Thomas elements are only implemented for grids with simplex or cube elements.");
223 GeometryType type = gv.template begin<0>()->type();
224 const static int dofsPerElement = type.isCube() ? ((dim == 2) ? k*(k+1)*dim : k*(k+1)*(k+1)*dim) : k*dim;
225 const static int dofsPerFace = type.isCube() ? (dim-2)*2*k+k+1 : (dim-1)*k+1 ;
227 dofsPerCodim_ = {{dofsPerElement, dofsPerFace}};
230 void initializeIndices()
233 codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
238 const GridView& gridView()
const
244 void update (
const GridView& gv)
252 Node makeNode()
const
254 return Node{&finiteElementMap_};
257 size_type dimension()
const
259 return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1);
262 size_type maxNodeSize()
const
264 size_type result = 0;
265 for (
auto&& type : gridView_.indexSet().types(0))
267 size_t numFaces = ReferenceElements<double,dim>::general(type).size(1);
268 const static int dofsPerElement = type.isCube() ? ((dim == 2) ? k*(k+1)*dim : k*(k+1)*(k+1)*dim) : k*dim;
269 const static int dofsPerFace = type.isCube() ? (dim-2)*2*k+k+1 : (dim-1)*k+1 ;
270 result = std::max(result, dofsPerElement + dofsPerFace * numFaces);
281 template<
typename It>
282 It indices(
const Node& node, It it)
const
284 const auto& gridIndexSet = gridView().indexSet();
285 const auto& element = node.element();
288 if (not(element.type().isCube()) and not(element.type().isSimplex()))
289 DUNE_THROW(Dune::NotImplemented,
"RaviartThomasBasis only implemented for cube and simplex elements.");
291 for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
293 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
296 size_t subentity = localKey.subEntity();
297 size_t codim = localKey.codim();
299 if (not(codim==0 or codim==1))
300 DUNE_THROW(Dune::NotImplemented,
"Grid contains elements not supported for the RaviartThomasBasis");
302 *it = { codimOffset_[codim] +
303 dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
311 std::array<size_t,dim+1> codimOffset_;
312 FiniteElementMap finiteElementMap_;
314 std::array<int,dim+1> dofsPerCodim_;
319template<
typename GV,
int k>
320class RaviartThomasNode :
323 static const int dim = GV::dimension;
327 using size_type = std::size_t;
328 using Element =
typename GV::template Codim<0>::Entity;
329 using FiniteElementMap =
typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
330 using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::ContravariantPiolaTransformator,
331 typename FiniteElementMap::FiniteElement,
334 RaviartThomasNode(
const FiniteElementMap* finiteElementMap) :
336 finiteElementMap_(finiteElementMap)
340 const Element& element()
const
349 const FiniteElement& finiteElement()
const
351 return finiteElement_;
355 void bind(
const Element& e)
358 finiteElement_.bind((finiteElementMap_->find(*element_)), e);
359 this->setSize(finiteElement_.size());
364 FiniteElement finiteElement_;
365 const Element* element_;
366 const FiniteElementMap* finiteElementMap_;
369namespace BasisFactory {
378template<std::
size_t k>
381 return [](
const auto& gridView) {
382 return RaviartThomasPreBasis<std::decay_t<
decltype(gridView)>, k>(gridView);
401template<
typename GV,
int k>
402using RaviartThomasBasis = DefaultGlobalBasis<RaviartThomasPreBasis<GV, k> >;
auto raviartThomas()
Create a pre-basis factory that can create a Raviart-Thomas pre-basis.
Definition: raviartthomasbasis.hh:379
auto power(ChildPreBasisFactory &&childPreBasisFactory, std::size_t k, const IndexMergingStrategy &)
Create a pre-basis factory that can build a PowerPreBasis.
Definition: dynamicpowerbasis.hh:409
Definition: polynomial.hh:17