3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
12#include <dune/localfunctions/common/localfiniteelementvariant.hh>
13#include <dune/localfunctions/raviartthomas.hh>
14#include <dune/localfunctions/raviartthomas/raviartthomas0cube2d.hh>
15#include <dune/localfunctions/raviartthomas/raviartthomas0cube3d.hh>
16#include <dune/localfunctions/raviartthomas/raviartthomas02d.hh>
17#include <dune/localfunctions/raviartthomas/raviartthomas03d.hh>
18#include <dune/localfunctions/raviartthomas/raviartthomas1cube2d.hh>
19#include <dune/localfunctions/raviartthomas/raviartthomas1cube3d.hh>
20#include <dune/localfunctions/raviartthomas/raviartthomas12d.hh>
21#include <dune/localfunctions/raviartthomas/raviartthomas2cube2d.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 RaviartThomasSimplexLocalInfo
37 using FiniteElement =
void*;
40 template<
typename D,
typename R>
41 struct RaviartThomasSimplexLocalInfo<2,D,R,0>
43 using FiniteElement = RT02DLocalFiniteElement<D,R>;
46 template<
typename D,
typename R>
47 struct RaviartThomasSimplexLocalInfo<2,D,R,1>
49 using FiniteElement = RT12DLocalFiniteElement<D,R>;
52 template<
typename D,
typename R>
53 struct RaviartThomasSimplexLocalInfo<3,D,R,0>
55 using FiniteElement = RT03DLocalFiniteElement<D,R>;
58 template<
int dim,
typename D,
typename R, std::
size_t k>
59 struct RaviartThomasCubeLocalInfo
62 using FiniteElement =
void*;
65 template<
typename D,
typename R>
66 struct RaviartThomasCubeLocalInfo<2,D,R,0>
68 using FiniteElement = RT0Cube2DLocalFiniteElement<D,R>;
71 template<
typename D,
typename R>
72 struct RaviartThomasCubeLocalInfo<2,D,R,1>
74 using FiniteElement = RT1Cube2DLocalFiniteElement<D,R>;
77 template<
typename D,
typename R>
78 struct RaviartThomasCubeLocalInfo<2,D,R,2>
80 using FiniteElement = RT2Cube2DLocalFiniteElement<D,R>;
83 template<
typename D,
typename R>
84 struct RaviartThomasCubeLocalInfo<3,D,R,0>
86 using FiniteElement = RT0Cube3DLocalFiniteElement<D,R>;
89 template<
typename D,
typename R>
90 struct RaviartThomasCubeLocalInfo<3,D,R,1>
92 using FiniteElement = RT1Cube3DLocalFiniteElement<D,R>;
95 template<
typename GV,
int dim,
typename R, std::
size_t k>
96 class RaviartThomasLocalFiniteElementMap
98 using D =
typename GV::ctype;
99 constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
101 using CubeFiniteElement =
typename RaviartThomasCubeLocalInfo<dim, D, R, k>::FiniteElement;
102 using SimplexFiniteElement =
typename RaviartThomasSimplexLocalInfo<dim, D, R, k>::FiniteElement;
106 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
108 constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId;
111 using FiniteElement = std::conditional_t<hasFixedElementType,
112 std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
113 LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
117 static std::size_t numVariants(GeometryType type)
119 auto numFacets = referenceElement<D,dim>(type).size(1);
120 return power(2,numFacets);
123 RaviartThomasLocalFiniteElementMap(
const GV& gv)
127 if constexpr (hasFixedElementType)
129 variants_.resize(numVariants(type));
130 for (
size_t i = 0; i < numVariants(type); i++)
131 variants_[i] = FiniteElement(i);
138 variants_[i] = SimplexFiniteElement(i);
143 for(
const auto& cell : elements(gv))
145 unsigned int myId = elementMapper_.index(cell);
148 for (
const auto& intersection : intersections(gv,cell))
150 if (intersection.neighbor() && (elementMapper_.index(intersection.outside()) > myId))
151 orient_[myId] |= (1 << intersection.indexInInside());
155 if constexpr (!hasFixedElementType)
156 if (cell.type().isCube())
161 template<
class EntityType>
162 const FiniteElement& find(
const EntityType& e)
const
164 return variants_[orient_[elementMapper_.index(e)]];
168 std::vector<FiniteElement> variants_;
170 std::vector<unsigned char> orient_;
188template<
typename GV,
int k>
189class RaviartThomasNode;
191template<
typename GV,
int k>
192class RaviartThomasPreBasis :
193 public LeafPreBasisMixin< RaviartThomasPreBasis<GV,k> >
195 static const int dim = GV::dimension;
196 using FiniteElementMap =
typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
202 using size_type = std::size_t;
204 using Node = RaviartThomasNode<GV, k>;
207 RaviartThomasPreBasis(
const GridView& gv) :
209 finiteElementMap_(gv)
212 if (gv.indexSet().types(0).size() > 1 and k>0)
215 for(
auto type : gv.indexSet().types(0))
216 if (!type.isSimplex() && !type.isCube())
220 const static int dofsPerElement = type.isCube() ? ((dim == 2) ? k*(k+1)*dim : k*(k+1)*(k+1)*dim) : k*dim;
221 const static int dofsPerFace = type.isCube() ? (dim-2)*2*k+k+1 : (dim-1)*k+1 ;
223 dofsPerCodim_ = {{dofsPerElement, dofsPerFace}};
226 void initializeIndices()
229 codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
234 const GridView& gridView()
const
240 void update (
const GridView& gv)
248 Node makeNode()
const
250 return Node{&finiteElementMap_};
253 size_type dimension()
const
255 return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1);
258 size_type maxNodeSize()
const
260 size_type result = 0;
261 for (
auto&& type : gridView_.indexSet().types(0))
264 const static int dofsPerElement = type.isCube() ? ((dim == 2) ? k*(k+1)*dim : k*(k+1)*(k+1)*dim) : k*dim;
265 const static int dofsPerFace = type.isCube() ? (dim-2)*2*k+k+1 : (dim-1)*k+1 ;
266 result =
std::max(result, dofsPerElement + dofsPerFace * numFaces);
277 template<
typename It>
278 It indices(
const Node& node, It it)
const
280 const auto& gridIndexSet = gridView().indexSet();
281 const auto& element = node.element();
284 if (not(element.type().isCube()) and not(element.type().isSimplex()))
287 for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
289 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
293 size_t codim = localKey.
codim();
295 if (not(codim==0 or codim==1))
298 *it = { codimOffset_[codim] +
299 dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.
index() };
307 std::array<size_t,dim+1> codimOffset_;
308 FiniteElementMap finiteElementMap_;
310 std::array<int,dim+1> dofsPerCodim_;
315template<
typename GV,
int k>
316class RaviartThomasNode :
319 static const int dim = GV::dimension;
323 using size_type = std::size_t;
324 using Element =
typename GV::template Codim<0>::Entity;
325 using FiniteElementMap =
typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
326 using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::ContravariantPiolaTransformator,
327 typename FiniteElementMap::FiniteElement,
330 RaviartThomasNode(
const FiniteElementMap* finiteElementMap) :
332 finiteElementMap_(finiteElementMap)
336 const Element& element()
const
345 const FiniteElement& finiteElement()
const
347 return finiteElement_;
351 void bind(
const Element& e)
354 finiteElement_.bind((finiteElementMap_->find(*element_)), e);
355 this->setSize(finiteElement_.size());
360 FiniteElement finiteElement_;
361 const Element* element_;
362 const FiniteElementMap* finiteElementMap_;
365namespace BasisFactory {
374template<std::
size_t k>
377 return [](
const auto& gridView) {
378 return RaviartThomasPreBasis<std::decay_t<
decltype(gridView)>, k>(gridView);
397template<
typename GV,
int k>
398using RaviartThomasBasis = DefaultGlobalBasis<RaviartThomasPreBasis<GV, k> >;
Describe position of one degree of freedom.
Definition: localkey.hh:24
constexpr unsigned int index() const noexcept
Return offset within subentity.
Definition: localkey.hh:70
constexpr unsigned int codim() const noexcept
Return codim of associated entity.
Definition: localkey.hh:63
constexpr unsigned int subEntity() const noexcept
Return number of associated subentity.
Definition: localkey.hh:56
Implementation class for a multiple codim and multiple geometry type mapper.
Definition: mcmgmapper.hh:129
Default exception for dummy implementations.
Definition: exceptions.hh:263
A set of traits classes to store static information about grid implementation.
A few common exception classes.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
auto raviartThomas()
Create a pre-basis factory that can create a Raviart-Thomas pre-basis.
Definition: raviartthomasbasis.hh:375
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
MCMGLayout mcmgElementLayout()
layout for elements (codim-0 entities)
Definition: mcmgmapper.hh:97
Mapper for multiple codim and multiple geometry types.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:156