3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_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>
25#include <dune/functions/functionspacebases/nodes.hh>
26#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
27#include <dune/functions/functionspacebases/flatmultiindex.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;
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,
171 const FiniteElement& find(
const EntityType& e)
const
173 return *cubeVariant_[orient_[is_->index(e)]];
181 template<
class EntityType,
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_;
211template<
typename GV,
int k>
212class RaviartThomasNode;
214template<
typename GV,
int k,
class MI>
215class RaviartThomasNodeIndexSet;
217template<
typename GV,
int k,
class MI>
218class RaviartThomasPreBasis
220 static const int dim = GV::dimension;
221 using FiniteElementMap =
typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
223 template<
typename,
int,
class>
224 friend class RaviartThomasNodeIndexSet;
230 using size_type = std::size_t;
232 using Node = RaviartThomasNode<GV, k>;
234 using IndexSet = RaviartThomasNodeIndexSet<GV, k, MI>;
237 using MultiIndex = MI;
242 RaviartThomasPreBasis(
const GridView& gv) :
244 finiteElementMap_(gv)
248 if (gv.indexSet().types(0).size() > 1)
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;
255 dofsPerCodim_ = {{dofsPerElement, dofsPerFace}};
258 void initializeIndices()
261 codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
266 const GridView& gridView()
const
272 void update (
const GridView& gv)
280 Node makeNode()
const
282 return Node{&finiteElementMap_};
291 IndexSet makeIndexSet()
const
293 return IndexSet{*
this};
296 size_type size()
const
298 return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1);
302 size_type size(
const SizePrefix prefix)
const
304 assert(prefix.size() == 0 || prefix.size() == 1);
305 return (prefix.size() == 0) ? size() : 0;
309 size_type dimension()
const
314 size_type maxNodeSize()
const
316 size_type result = 0;
317 for (
auto&& type : gridView_.indexSet().types(0))
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);
329 const GridView gridView_;
330 std::array<size_t,dim+1> codimOffset_;
331 FiniteElementMap finiteElementMap_;
333 std::array<int,dim+1> dofsPerCodim_;
338template<
typename GV,
int k>
339class RaviartThomasNode :
342 static const int dim = GV::dimension;
346 using size_type = std::size_t;
347 using Element =
typename GV::template Codim<0>::Entity;
348 using FiniteElementMap =
typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
349 using FiniteElement =
typename FiniteElementMap::FiniteElement;
351 RaviartThomasNode(
const FiniteElementMap* finiteElementMap) :
352 finiteElement_(nullptr),
354 finiteElementMap_(finiteElementMap)
358 const Element& element()
const
367 const FiniteElement& finiteElement()
const
369 return *finiteElement_;
373 void bind(
const Element& e)
376 finiteElement_ = &(finiteElementMap_->find(*element_));
377 this->setSize(finiteElement_->size());
382 const FiniteElement* finiteElement_;
383 const Element* element_;
384 const FiniteElementMap* finiteElementMap_;
387template<
typename GV,
int k,
class MI>
388class RaviartThomasNodeIndexSet
390 enum {dim = GV::dimension};
394 using size_type = std::size_t;
397 using MultiIndex = MI;
399 using PreBasis = RaviartThomasPreBasis<GV, k, MI>;
401 using Node = RaviartThomasNode<GV,k>;
403 RaviartThomasNodeIndexSet(
const PreBasis& preBasis) :
412 void bind(
const Node& node)
426 size_type size()
const
428 return node_->finiteElement().size();
436 template<
typename It>
437 It indices(It it)
const
439 const auto& gridIndexSet = preBasis_->gridView().indexSet();
440 const auto& element = node_->element();
443 if (not(element.type().isCube()) and not(element.type().isSimplex()))
446 for(std::size_t i=0, end=size(); i<end; ++i, ++it)
448 Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
452 size_t codim = localKey.
codim();
454 if (not(codim==0 or codim==1))
457 *it = { preBasis_->codimOffset_[codim] +
458 preBasis_->dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.
index() };
465 const PreBasis* preBasis_;
471namespace BasisFactory {
475template<std::
size_t k>
476class RaviartThomasPreBasisFactory
479 static const std::size_t requiredMultiIndexSize=1;
481 template<
class MultiIndex,
class Gr
idView>
482 auto makePreBasis(
const GridView& gridView)
const
484 return RaviartThomasPreBasis<GridView, k, MultiIndex>(gridView);
498template<std::
size_t k>
501 return Imp::RaviartThomasPreBasisFactory<k>();
519template<
typename GV,
int k>
520using RaviartThomasBasis = DefaultGlobalBasis<RaviartThomasPreBasis<GV, k, FlatMultiIndex<std::size_t>> >;
Describe position of one degree of freedom.
Definition: localkey.hh:21
unsigned int index() const
Return offset within subentity.
Definition: localkey.hh:66
unsigned int codim() const
Return codim of associated entity.
Definition: localkey.hh:60
unsigned int subEntity() const
Return number of associated subentity.
Definition: localkey.hh:54
Default exception for dummy implementations.
Definition: exceptions.hh:261
A Vector class with statically reserved memory.
Definition: reservedvector.hh:43
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:180
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
auto raviartThomas()
Create a pre-basis factory that can create a Raviart-Thomas pre-basis.
Definition: raviartthomasbasis.hh:499
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune namespace.
Definition: alignedallocator.hh:14
static const bool value
always a false value
Definition: typetraits.hh:128
static const bool value
always a true value
Definition: typetraits.hh:141
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196