3 #ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
4 #define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
7 #include <dune/common/exceptions.hh>
9 #include <dune/grid/common/capabilities.hh>
11 #include <dune/localfunctions/common/virtualinterface.hh>
12 #include <dune/localfunctions/common/virtualwrappers.hh>
14 #include <dune/localfunctions/raviartthomas.hh>
15 #include <dune/localfunctions/raviartthomas/raviartthomas0cube2d.hh>
16 #include <dune/localfunctions/raviartthomas/raviartthomas0cube3d.hh>
17 #include <dune/localfunctions/raviartthomas/raviartthomas02d.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/typetree/leafnode.hh>
34 template<
int dim,
typename D,
typename R, std::
size_t k>
35 struct RaviartThomasSimplexLocalInfo
37 static_assert((AlwaysFalse<D>::value),
"The requested type of Raviart-Thomas element is not implemented, sorry!");
40 template<
typename D,
typename R>
41 struct RaviartThomasSimplexLocalInfo<2,D,R,0>
43 using FiniteElement = RT02DLocalFiniteElement<D,R>;
44 static const std::size_t Variants = 8;
47 template<
typename D,
typename R>
48 struct RaviartThomasSimplexLocalInfo<2,D,R,1>
50 using FiniteElement = RT12DLocalFiniteElement<D,R>;
51 static const std::size_t Variants = 8;
54 template<
int dim,
typename D,
typename R, std::
size_t k>
55 struct RaviartThomasCubeLocalInfo
57 static_assert((AlwaysFalse<D>::value),
"The requested type of Raviart-Thomas element is not implemented, sorry!");
60 template<
typename D,
typename R>
61 struct RaviartThomasCubeLocalInfo<2,D,R,0>
63 using FiniteElement = RT0Cube2DLocalFiniteElement<D,R>;
64 static const std::size_t Variants = 16;
67 template<
typename D,
typename R>
68 struct RaviartThomasCubeLocalInfo<2,D,R,1>
70 using FiniteElement = RT1Cube2DLocalFiniteElement<D,R>;
71 static const std::size_t Variants = 16;
74 template<
typename D,
typename R>
75 struct RaviartThomasCubeLocalInfo<2,D,R,2>
77 using FiniteElement = RT2Cube2DLocalFiniteElement<D,R>;
78 static const std::size_t Variants = 16;
81 template<
typename D,
typename R>
82 struct RaviartThomasCubeLocalInfo<3,D,R,0>
84 using FiniteElement = RT0Cube3DLocalFiniteElement<D,R>;
85 static const std::size_t Variants = 64;
88 template<
typename D,
typename R>
89 struct RaviartThomasCubeLocalInfo<3,D,R,1>
91 using FiniteElement = RT1Cube3DLocalFiniteElement<D,R>;
92 static const std::size_t Variants = 64;
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;
103 using CubeFiniteElementImp =
typename std::conditional<hasFixedElementType,
105 LocalFiniteElementVirtualImp<CubeFiniteElement> >::type;
106 using SimplexFiniteElementImp =
typename std::conditional<hasFixedElementType,
107 SimplexFiniteElement,
108 LocalFiniteElementVirtualImp<SimplexFiniteElement> >::type;
112 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
114 constexpr
static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId;
115 constexpr
static GeometryType type = GeometryType(topologyId, GV::dimension);
117 using FiniteElement =
typename std::conditional<hasFixedElementType,
118 typename std::conditional<type.isCube(),CubeFiniteElement,SimplexFiniteElement>::type,
119 LocalFiniteElementVirtualInterface<T> >::type;
121 RaviartThomasLocalFiniteElementMap(
const GV& gv)
122 : is_(&(gv.indexSet())), orient_(gv.size(0))
124 cubeVariant_.resize(RaviartThomasCubeLocalInfo<dim, D, R, k>::Variants);
125 simplexVariant_.resize(RaviartThomasSimplexLocalInfo<dim, D, R, k>::Variants);
128 for (
size_t i = 0; i < cubeVariant_.size(); i++)
129 cubeVariant_[i] = std::make_unique<CubeFiniteElementImp>(CubeFiniteElement(i));
131 for (
size_t i = 0; i < simplexVariant_.size(); i++)
132 simplexVariant_[i] = std::make_unique<SimplexFiniteElementImp>(SimplexFiniteElement(i));
136 for(
const auto& cell : elements(gv))
138 unsigned int myId = is_->index(cell);
141 for (
const auto& intersection : intersections(gv,cell))
143 if (intersection.neighbor() && (is_->index(intersection.outside()) > myId))
144 orient_[myId] |= (1 << intersection.indexInInside());
154 template<
class EntityType,
155 std::enable_if_t<!hasFixedElementType and AlwaysTrue<EntityType>::value,
int> = 0>
156 const FiniteElement& find(
const EntityType& e)
const
158 if (e.type().isCube())
159 return *cubeVariant_[orient_[is_->index(e)]];
161 return *simplexVariant_[orient_[is_->index(e)]];
169 template<
class EntityType,
170 std::enable_if_t<hasFixedElementType and type.isCube() and AlwaysTrue<EntityType>::value,
int> = 0>
171 const FiniteElement& find(
const EntityType& e)
const
173 return *cubeVariant_[orient_[is_->index(e)]];
181 template<
class EntityType,
182 std::enable_if_t<hasFixedElementType and type.isSimplex() and AlwaysTrue<EntityType>::value,
int> = 0>
183 const FiniteElement& find(
const EntityType& e)
const
185 return *simplexVariant_[orient_[is_->index(e)]];
189 std::vector<std::unique_ptr<CubeFiniteElementImp> > cubeVariant_;
190 std::vector<std::unique_ptr<SimplexFiniteElementImp> > simplexVariant_;
191 const typename GV::IndexSet* is_;
192 std::vector<unsigned char> orient_;
211 template<
typename GV,
int k>
214 template<
typename GV,
int k,
class MI>
217 template<
typename GV,
int k,
class MI>
220 static const int dim = GV::dimension;
221 using FiniteElementMap =
typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
223 template<
typename,
int,
class>
248 if (gv.indexSet().types(0).size() > 1)
249 DUNE_THROW(Dune::NotImplemented,
"Raviart-Thomas basis is only implemented for grids with a single element type");
251 GeometryType type = gv.template begin<0>()->type();
252 const static int dofsPerElement = (dim == 2) ? (type.isCube() ? k*(k+1)*dim : k*dim) : k*(k+1)*(k+1)*dim;
253 constexpr
int dofsPerFace = (dim == 2) ? k+1 : 3*k+1;
304 assert(prefix.size() == 0 || prefix.size() == 1);
305 return (prefix.size() == 0) ?
size() : 0;
317 for (
auto&& type :
gridView_.indexSet().types(0))
319 size_t numFaces = ReferenceElements<double,dim>::general(type).size(1);
320 const static int dofsPerCodim0 = (dim == 2) ? (type.isCube() ? k*(k+1)*dim : k*dim) : k*(k+1)*(k+1)*dim;
321 constexpr
int dofsPerCodim1 = (dim == 2) ? k+1 : 3*k+1;
322 result = std::max(result, dofsPerCodim0 + dofsPerCodim1 * numFaces);
338 template<
typename GV,
int k>
342 static const int dim = GV::dimension;
347 using Element =
typename GV::template Codim<0>::Entity;
348 using FiniteElementMap =
typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
387 template<
typename GV,
int k,
class MI>
390 enum {dim = GV::dimension};
436 template<
typename It>
443 if (not(element.type().isCube()) and not(element.type().isSimplex()))
444 DUNE_THROW(Dune::NotImplemented,
"RaviartThomasBasis only implemented for cube and simplex elements.");
446 for(std::size_t i=0, end=
size(); i<end; ++i, ++it)
451 size_t subentity = localKey.subEntity();
452 size_t codim = localKey.codim();
454 if (not(codim==0 or codim==1))
455 DUNE_THROW(Dune::NotImplemented,
"Grid contains elements not supported for the RaviartThomasBasis");
471 namespace BasisFactory {
475 template<std::
size_t k>
476 class RaviartThomasPreBasisFactory
479 static const std::size_t requiredMultiIndexSize=1;
481 template<
class MultiIndex,
class Gr
idView>
482 auto makePreBasis(
const GridView& gridView)
const
498 template<std::
size_t k>
501 return Imp::RaviartThomasPreBasisFactory<k>();
519 template<
typename GV,
int k>
526 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH