DUNE PDELab (git)

brezzidouglasmarinibasis.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_BREZZIDOUGLASMARINIBASIS_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINIBASIS_HH
9
10#include <array>
12#include <dune/geometry/referenceelements.hh>
13
14#include <dune/localfunctions/common/virtualinterface.hh>
15#include <dune/localfunctions/common/virtualwrappers.hh>
16
17#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1cube2d.hh>
18#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1cube3d.hh>
19#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1simplex2d.hh>
20#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini2cube2d.hh>
21#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini2simplex2d.hh>
22
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>
27
28namespace Dune {
29namespace Functions {
30
31namespace Impl {
32
33 template<int dim, typename D, typename R, std::size_t k>
34 struct BDMSimplexLocalInfo
35 {
36 static_assert((AlwaysFalse<D>::value),"The requested type of BDM element is not implemented, sorry!");
37 };
38
39 template<typename D, typename R>
40 struct BDMSimplexLocalInfo<2,D,R,1>
41 {
42 using FiniteElement = BDM1Simplex2DLocalFiniteElement<D,R>;
43 static const std::size_t Variants = 8;
44 };
45
46 template<typename D, typename R>
47 struct BDMSimplexLocalInfo<2,D,R,2>
48 {
49 using FiniteElement = BDM2Simplex2DLocalFiniteElement<D,R>;
50 static const std::size_t Variants = 8;
51 };
52
53 template<int dim, typename D, typename R, std::size_t k>
54 struct BDMCubeLocalInfo
55 {
56 static_assert((AlwaysFalse<D>::value),"The requested type of BDM element is not implemented, sorry!");
57 };
58
59 template<typename D, typename R>
60 struct BDMCubeLocalInfo<2,D,R,1>
61 {
62 using FiniteElement = BDM1Cube2DLocalFiniteElement<D,R>;
63 static const std::size_t Variants = 16;
64 };
65
66 template<typename D, typename R>
67 struct BDMCubeLocalInfo<2,D,R,2>
68 {
69 using FiniteElement = BDM2Cube2DLocalFiniteElement<D,R>;
70 static const std::size_t Variants = 16;
71 };
72
73 template<typename D, typename R>
74 struct BDMCubeLocalInfo<3,D,R,1>
75 {
76 using FiniteElement = BDM1Cube3DLocalFiniteElement<D,R>;
77 static const std::size_t Variants = 64;
78 };
79
80 template<typename GV, int dim, typename R, std::size_t k>
81 class BDMLocalFiniteElementMap
82 {
83 using D = typename GV::ctype;
84 using CubeFiniteElement = typename BDMCubeLocalInfo<dim, D, R, k>::FiniteElement;
85 using SimplexFiniteElement = typename BDMSimplexLocalInfo<dim, D, R, k>::FiniteElement;
86
87 public:
88
89 using T = LocalBasisTraits<D, dim, FieldVector<D,dim>, R, dim, FieldVector<R,dim>, FieldMatrix<D,dim,dim> >;
90 using FiniteElement = LocalFiniteElementVirtualInterface<T>;
91
92 BDMLocalFiniteElementMap(const GV& gv)
93 : is_(&(gv.indexSet())), orient_(gv.size(0))
94 {
95 cubeVariant_.resize(BDMCubeLocalInfo<dim, D, R, k>::Variants);
96 simplexVariant_.resize(BDMSimplexLocalInfo<dim, D, R, k>::Variants);
97
98 // create all variants
99 for (size_t i = 0; i < cubeVariant_.size(); i++)
100 cubeVariant_[i] = std::make_shared<LocalFiniteElementVirtualImp<CubeFiniteElement> >(CubeFiniteElement(i));
101
102 for (size_t i = 0; i < simplexVariant_.size(); i++)
103 simplexVariant_[i] = std::make_shared<LocalFiniteElementVirtualImp<SimplexFiniteElement> >(SimplexFiniteElement(i));
104
105 // compute orientation for all elements
106 // loop once over the grid
107 for(const auto& cell : elements(gv))
108 {
109 unsigned int myId = is_->index(cell);
110 orient_[myId] = 0;
111
112 for (const auto& intersection : intersections(gv,cell))
113 {
114 if (intersection.neighbor() && (is_->index(intersection.outside()) > myId))
115 orient_[myId] |= (1 << intersection.indexInInside());
116 }
117 }
118 }
119
121 template<class EntityType>
122 const FiniteElement& find(const EntityType& e) const
123 {
124 if (e.type().isCube())
125 return *cubeVariant_[orient_[is_->index(e)]];
126 else
127 return *simplexVariant_[orient_[is_->index(e)]];
128 }
129
130 private:
131 std::vector<std::shared_ptr<LocalFiniteElementVirtualImp<CubeFiniteElement> > > cubeVariant_;
132 std::vector<std::shared_ptr<LocalFiniteElementVirtualImp<SimplexFiniteElement> > > simplexVariant_;
133 const typename GV::IndexSet* is_;
134 std::vector<unsigned char> orient_;
135 };
136
137
138} // namespace Impl
139
140
141// *****************************************************************************
142// This is the reusable part of the basis. It contains
143//
144// BrezziDouglasMariniPreBasis
145// BrezziDouglasMariniNode
146//
147// The pre-basis allows to create the others and is the owner of possible shared
148// state. These components do _not_ depend on the global basis and local view
149// and can be used without a global basis.
150// *****************************************************************************
151
152template<typename GV, int k>
153class BrezziDouglasMariniNode;
154
155template<typename GV, int k>
156class BrezziDouglasMariniPreBasis :
157 public LeafPreBasisMixin< BrezziDouglasMariniPreBasis<GV,k> >
158{
159 static const int dim = GV::dimension;
160 using FiniteElementMap = typename Impl::BDMLocalFiniteElementMap<GV, dim, double, k>;
161
162public:
163
165 using GridView = GV;
166 using size_type = std::size_t;
167
168 using Node = BrezziDouglasMariniNode<GV, k>;
169
171 BrezziDouglasMariniPreBasis(const GridView& gv) :
172 gridView_(gv),
173 finiteElementMap_(gv)
174 {
175 // There is no inherent reason why the basis shouldn't work for grids with more than one
176 // element types. Somebody simply has to sit down and implement the missing bits.
177 if (gv.indexSet().types(0).size() > 1)
178 DUNE_THROW(Dune::NotImplemented, "Brezzi-Douglas-Marini basis is only implemented for grids with a single element type");
179 }
180
181 void initializeIndices()
182 {
183 codimOffset_[0] = 0;
184 codimOffset_[1] = codimOffset_[0] + dofsPerCodim_[0] * gridView_.size(0);
185 //if (dim==3) codimOffset_[2] = codimOffset_[1] + dofsPerCodim[1] * gridView_.size(1);
186 }
187
190 const GridView& gridView() const
191 {
192 return gridView_;
193 }
194
195 /* \brief Update the stored grid view, to be called if the grid has changed */
196 void update (const GridView& gv)
197 {
198 gridView_ = gv;
199 }
200
204 Node makeNode() const
205 {
206 return Node{&finiteElementMap_};
207 }
208
209 size_type dimension() const
210 {
211 return dofsPerCodim_[0] * gridView_.size(0) + dofsPerCodim_[1] * gridView_.size(1); // only 2d
212 }
213
214 size_type maxNodeSize() const
215 {
216 // The implementation currently only supports grids with a single element type.
217 // We can therefore return the actual number of dofs here.
218 GeometryType elementType = *(gridView_.indexSet().types(0).begin());
219 size_t numFaces = ReferenceElements<double,dim>::general(elementType).size(1);
220 return dofsPerCodim_[0] + dofsPerCodim_[1] * numFaces;
221 }
222
228 template<typename It>
229 It indices(const Node& node, It it) const
230 {
231 const auto& gridIndexSet = gridView().indexSet();
232 const auto& element = node.element();
233
234 // throw if element is not of predefined type
235 if (not(element.type().isCube()) and not(element.type().isSimplex()))
236 DUNE_THROW(Dune::NotImplemented, "BrezziDouglasMariniBasis only implemented for cube and simplex elements.");
237
238 for(std::size_t i=0, end=node.size(); i<end; ++i, ++it)
239 {
240 Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
241
242 // The dimension of the entity that the current dof is related to
243 size_t subentity = localKey.subEntity();
244 size_t codim = localKey.codim();
245
246 *it = { codimOffset_[codim] +
247 dofsPerCodim_[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
248 }
249
250 return it;
251 }
252
253protected:
254 GridView gridView_;
255 std::array<size_t,dim+1> codimOffset_;
256 FiniteElementMap finiteElementMap_;
257 // Number of dofs per entity type depending on the entity's codimension and type
258 std::array<int,2> dofsPerCodim_ {{dim*(k-1)*3, dim+(k-1)}};
259};
260
261
262
263template<typename GV, int k>
264class BrezziDouglasMariniNode :
265 public LeafBasisNode
266{
267 static const int dim = GV::dimension;
268
269public:
270
271 using size_type = std::size_t;
272 using Element = typename GV::template Codim<0>::Entity;
273 using FiniteElementMap = typename Impl::BDMLocalFiniteElementMap<GV, dim, double, k>;
274 using FiniteElement = Impl::GlobalValuedLocalFiniteElement<Impl::ContravariantPiolaTransformator,
275 typename FiniteElementMap::FiniteElement,
276 Element>;
277
278 BrezziDouglasMariniNode(const FiniteElementMap* finiteElementMap) :
279 element_(nullptr),
280 finiteElementMap_(finiteElementMap)
281 {}
282
284 const Element& element() const
285 {
286 return *element_;
287 }
288
293 const FiniteElement& finiteElement() const
294 {
295 return finiteElement_;
296 }
297
299 void bind(const Element& e)
300 {
301 element_ = &e;
302 finiteElement_.bind((finiteElementMap_->find(*element_)), e);
303 this->setSize(finiteElement_.size());
304 }
305
306protected:
307
308 FiniteElement finiteElement_;
309 const Element* element_;
310 const FiniteElementMap* finiteElementMap_;
311};
312
313
314
315namespace BasisFactory {
316
324template<std::size_t k>
326{
327 return [](const auto& gridView) {
328 return BrezziDouglasMariniPreBasis<std::decay_t<decltype(gridView)>, k>(gridView);
329 };
330}
331
332} // end namespace BasisFactory
333
334
335
336// *****************************************************************************
337// This is the actual global basis implementation based on the reusable parts.
338// *****************************************************************************
339
347template<typename GV, int k>
348using BrezziDouglasMariniBasis = DefaultGlobalBasis<BrezziDouglasMariniPreBasis<GV, k> >;
349
350} // end namespace Functions
351} // end namespace Dune
352
353
354#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINIBASIS_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
Default exception for dummy implementations.
Definition: exceptions.hh:263
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 brezziDouglasMarini()
Create a pre-basis factory that can create a Brezzi-Douglas-Marini pre-basis.
Definition: brezzidouglasmarinibasis.hh:325
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
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)