DUNE PDELab (git)

raviartthomasbasis.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
9
10#include <array>
12
15
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>
26
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>
31
32namespace Dune {
33namespace Functions {
34
35namespace Impl {
36
37 template<int dim, typename D, typename R, std::size_t k>
38 struct RaviartThomasSimplexLocalInfo
39 {
40 // Dummy type, must be something that we can have a std::unique_ptr to
41 using FiniteElement = void*;
42 };
43
44 template<typename D, typename R>
45 struct RaviartThomasSimplexLocalInfo<2,D,R,0>
46 {
47 using FiniteElement = RT02DLocalFiniteElement<D,R>;
48 };
49
50 template<typename D, typename R>
51 struct RaviartThomasSimplexLocalInfo<2,D,R,1>
52 {
53 using FiniteElement = RT12DLocalFiniteElement<D,R>;
54 };
55
56 template<typename D, typename R>
57 struct RaviartThomasSimplexLocalInfo<3,D,R,0>
58 {
59 using FiniteElement = RT03DLocalFiniteElement<D,R>;
60 };
61
62 template<int dim, typename D, typename R, std::size_t k>
63 struct RaviartThomasCubeLocalInfo
64 {
65 // Dummy type, must be something that we can have a std::unique_ptr to
66 using FiniteElement = void*;
67 };
68
69 template<typename D, typename R>
70 struct RaviartThomasCubeLocalInfo<2,D,R,0>
71 {
72 using FiniteElement = RT0Cube2DLocalFiniteElement<D,R>;
73 };
74
75 template<typename D, typename R>
76 struct RaviartThomasCubeLocalInfo<2,D,R,1>
77 {
78 using FiniteElement = RT1Cube2DLocalFiniteElement<D,R>;
79 };
80
81 template<typename D, typename R>
82 struct RaviartThomasCubeLocalInfo<2,D,R,2>
83 {
84 using FiniteElement = RT2Cube2DLocalFiniteElement<D,R>;
85 };
86
87 template<typename D, typename R>
88 struct RaviartThomasCubeLocalInfo<3,D,R,0>
89 {
90 using FiniteElement = RT0Cube3DLocalFiniteElement<D,R>;
91 };
92
93 template<typename D, typename R>
94 struct RaviartThomasCubeLocalInfo<3,D,R,1>
95 {
96 using FiniteElement = RT1Cube3DLocalFiniteElement<D,R>;
97 };
98
99 template<typename GV, int dim, typename R, std::size_t k>
100 class RaviartThomasLocalFiniteElementMap
101 {
102 using D = typename GV::ctype;
103 constexpr static bool hasFixedElementType = Capabilities::hasSingleGeometryType<typename GV::Grid>::v;
104
105 using CubeFiniteElement = typename RaviartThomasCubeLocalInfo<dim, D, R, k>::FiniteElement;
106 using SimplexFiniteElement = typename RaviartThomasSimplexLocalInfo<dim, D, R, k>::FiniteElement;
107
108 public:
109
110 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
111
112 constexpr static unsigned int topologyId = Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId; // meaningless if hasFixedElementType is false
113 constexpr static GeometryType type = GeometryType(topologyId, GV::dimension);
114
115 using FiniteElement = std::conditional_t<hasFixedElementType,
116 std::conditional_t<type.isCube(),CubeFiniteElement,SimplexFiniteElement>,
117 LocalFiniteElementVariant<CubeFiniteElement, SimplexFiniteElement> >;
118
119 // Each element facet can have its orientation reversed, hence there are
120 // 2^#facets different variants.
121 static std::size_t numVariants(GeometryType type)
122 {
123 auto numFacets = referenceElement<D,dim>(type).size(1);
124 return power(2,numFacets);
125 }
126
127 RaviartThomasLocalFiniteElementMap(const GV& gv)
128 : elementMapper_(gv, mcmgElementLayout()),
129 orient_(gv.size(0))
130 {
131 if constexpr (hasFixedElementType)
132 {
133 variants_.resize(numVariants(type));
134 for (size_t i = 0; i < numVariants(type); i++)
135 variants_[i] = FiniteElement(i);
136 }
137 else
138 {
139 // for mixed grids add offset for cubes
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);
145 }
146
147 for(const auto& cell : elements(gv))
148 {
149 unsigned int myId = elementMapper_.index(cell);
150 orient_[myId] = 0;
151
152 for (const auto& intersection : intersections(gv,cell))
153 {
154 if (intersection.neighbor() && (elementMapper_.index(intersection.outside()) > myId))
155 orient_[myId] |= (1 << intersection.indexInInside());
156 }
157
158 // for mixed grids add offset for cubes
159 if constexpr (!hasFixedElementType)
160 if (cell.type().isCube())
161 orient_[myId] += numVariants(GeometryTypes::simplex(dim));
162 }
163 }
164
165 template<class EntityType>
166 const FiniteElement& find(const EntityType& e) const
167 {
168 return variants_[orient_[elementMapper_.index(e)]];
169 }
170
171 private:
172 std::vector<FiniteElement> variants_;
174 std::vector<unsigned char> orient_;
175 };
176
177
178} // namespace Impl
179
180
181// *****************************************************************************
182// This is the reusable part of the basis. It contains
183//
184// RaviartThomasPreBasis
185// RaviartThomasNode
186//
187// The pre-basis allows to create the others and is the owner of possible shared
188// state. These components do _not_ depend on the global basis and local view
189// and can be used without a global basis.
190// *****************************************************************************
191
192template<typename GV, int k>
193class RaviartThomasNode;
194
195template<typename GV, int k>
196class RaviartThomasPreBasis :
197 public LeafPreBasisMixin< RaviartThomasPreBasis<GV,k> >
198{
199 static const int dim = GV::dimension;
200 using FiniteElementMap = typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, double, k>;
201
202public:
203
205 using GridView = GV;
206 using size_type = std::size_t;
207
208 using Node = RaviartThomasNode<GV, k>;
209
211 RaviartThomasPreBasis(const GridView& gv) :
212 gridView_(gv),
213 finiteElementMap_(gv)
214 {
215 // Currently there are some unresolved bugs with hybrid grids and higher order Raviart-Thomas elements
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");
218
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.");
222
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 ;
226
227 dofsPerCodim_ = {{dofsPerElement, dofsPerFace}};
228 }
229
230 void initializeIndices()
231 {
232 codimOffset_[0] = 0;
233 codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
234 }
235
238 const GridView& gridView() const
239 {
240 return gridView_;
241 }
242
243 /* \brief Update the stored grid view, to be called if the grid has changed */
244 void update (const GridView& gv)
245 {
246 gridView_ = gv;
247 }
248
252 Node makeNode() const
253 {
254 return Node{&finiteElementMap_};
255 }
256
257 size_type dimension() const
258 {
259 return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1);
260 }
261
262 size_type maxNodeSize() const
263 {
264 size_type result = 0;
265 for (auto&& type : gridView_.indexSet().types(0))
266 {
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);
271 }
272
273 return result;
274 }
275
281 template<typename It>
282 It indices(const Node& node, It it) const
283 {
284 const auto& gridIndexSet = gridView().indexSet();
285 const auto& element = node.element();
286
287 // throw if Element is not of predefined type
288 if (not(element.type().isCube()) and not(element.type().isSimplex()))
289 DUNE_THROW(Dune::NotImplemented, "RaviartThomasBasis only implemented for cube and simplex elements.");
290
291 for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
292 {
293 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
294
295 // The dimension of the entity that the current dof is related to
296 size_t subentity = localKey.subEntity();
297 size_t codim = localKey.codim();
298
299 if (not(codim==0 or codim==1))
300 DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the RaviartThomasBasis");
301
302 *it = { codimOffset_[codim] +
303 dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
304 }
305
306 return it;
307 }
308
309protected:
310 GridView gridView_;
311 std::array<size_t,dim+1> codimOffset_;
312 FiniteElementMap finiteElementMap_;
313 // Number of dofs per entity type depending on the entity's codimension and type
314 std::array<int,dim+1> dofsPerCodim_;
315};
316
317
318
319template<typename GV, int k>
320class RaviartThomasNode :
321 public LeafBasisNode
322{
323 static const int dim = GV::dimension;
324
325public:
326
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,
332 Element>;
333
334 RaviartThomasNode(const FiniteElementMap* finiteElementMap) :
335 element_(nullptr),
336 finiteElementMap_(finiteElementMap)
337 { }
338
340 const Element& element() const
341 {
342 return *element_;
343 }
344
349 const FiniteElement& finiteElement() const
350 {
351 return finiteElement_;
352 }
353
355 void bind(const Element& e)
356 {
357 element_ = &e;
358 finiteElement_.bind((finiteElementMap_->find(*element_)), e);
359 this->setSize(finiteElement_.size());
360 }
361
362protected:
363
364 FiniteElement finiteElement_;
365 const Element* element_;
366 const FiniteElementMap* finiteElementMap_;
367};
368
369namespace BasisFactory {
370
378template<std::size_t k>
380{
381 return [](const auto& gridView) {
382 return RaviartThomasPreBasis<std::decay_t<decltype(gridView)>, k>(gridView);
383 };
384}
385
386} // end namespace BasisFactory
387
388
389
390// *****************************************************************************
391// This is the actual global basis implementation based on the reusable parts.
392// *****************************************************************************
393
401template<typename GV, int k>
402using RaviartThomasBasis = DefaultGlobalBasis<RaviartThomasPreBasis<GV, k> >;
403
404} // end namespace Functions
405} // end namespace Dune
406
407
408#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
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:379
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)