dune-localfunctions 2.8.0
Loading...
Searching...
No Matches
raviartthomassimplexinterpolation.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
4#define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
5
6#include <fstream>
7#include <utility>
8
9#include <dune/common/exceptions.hh>
10
11#include <dune/geometry/quadraturerules.hh>
12#include <dune/geometry/referenceelements.hh>
13#include <dune/geometry/type.hh>
14#include <dune/geometry/typeindex.hh>
15
20
21namespace Dune
22{
23
24 // Internal Forward Declarations
25 // -----------------------------
26
27 template < unsigned int dim, class Field >
28 struct RaviartThomasL2InterpolationFactory;
29
30
31
32 // LocalCoefficientsContainer
33 // --------------------------
34
35 class LocalCoefficientsContainer
36 {
37 typedef LocalCoefficientsContainer This;
38
39 public:
40 template <class Setter>
41 LocalCoefficientsContainer ( const Setter &setter )
42 {
43 setter.setLocalKeys(localKey_);
44 }
45
46 const LocalKey &localKey ( const unsigned int i ) const
47 {
48 assert( i < size() );
49 return localKey_[ i ];
50 }
51
52 std::size_t size () const
53 {
54 return localKey_.size();
55 }
56
57 private:
58 std::vector< LocalKey > localKey_;
59 };
60
61
62
63 // RaviartThomasCoefficientsFactory
64 // --------------------------------
65
66 template < unsigned int dim >
68 {
69 typedef std::size_t Key;
71
72 template< GeometryType::Id geometryId >
73 static Object *create( const Key &key )
74 {
75 typedef RaviartThomasL2InterpolationFactory< dim, double > InterpolationFactory;
76 if( !supports< geometryId >( key ) )
77 return nullptr;
78 typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< geometryId >( key );
79 Object *localKeys = new Object( *interpolation );
80 InterpolationFactory::release( interpolation );
81 return localKeys;
82 }
83
84 template< GeometryType::Id geometryId >
85 static bool supports ( const Key &key )
86 {
87 return GeometryType(geometryId).isSimplex();
88 }
89 static void release( Object *object ) { delete object; }
90 };
91
92
93
94 // RTL2InterpolationBuilder
95 // ------------------------
96
97 // L2 Interpolation requires:
98 // - for element
99 // - test basis
100 // - for each face (dynamic)
101 // - test basis
102 // - normal
103 template< unsigned int dim, class Field >
105 {
106 static const unsigned int dimension = dim;
107
108 // for the dofs associated to the element
111
112 // for the dofs associated to the faces
115
116 // the normals of the faces
117 typedef FieldVector< Field, dimension > Normal;
118
120
123
125 {
126 TestBasisFactory::release( testBasis_ );
127 for( FaceStructure &f : faceStructure_ )
129 }
130
131 [[deprecated("Use type().id() instead.")]]
132 unsigned int topologyId () const { return type().id(); }
133
134 GeometryType type () const { return geometry_; }
135
136 std::size_t order () const { return order_; }
137
138 // number of faces
139 unsigned int faceSize () const { return faceSize_; }
140
141 // basis associated to the element
142 TestBasis *testBasis () const { return testBasis_; }
143
144 // basis associated to face f
145 TestFaceBasis *testFaceBasis ( unsigned int f ) const { assert( f < faceSize() ); return faceStructure_[ f ].basis_; }
146
147 // normal of face f
148 const Normal &normal ( unsigned int f ) const { assert( f < faceSize() ); return *(faceStructure_[ f ].normal_); }
149
150 template< GeometryType::Id geometryId >
151 void build ( std::size_t order )
152 {
153 constexpr GeometryType geometry = geometryId;
154 geometry_ = geometry;
155 order_ = order;
156
157 testBasis_ = (order > 0 ? TestBasisFactory::template create< geometry >( order-1 ) : nullptr);
158
159 const auto &refElement = ReferenceElements< Field, dimension >::general( type() );
160 faceSize_ = refElement.size( 1 );
161 faceStructure_.reserve( faceSize_ );
162 for( unsigned int face = 0; face < faceSize_; ++face )
163 {
164 /* For simplices or cubes of arbitrary dimension you could just use
165 *
166 * ```
167 * GeometryType faceGeometry = Impl::getBase(geometry_);
168 * TestFaceBasis *faceBasis = TestFaceBasisFactory::template create< faceGeometry >( order );
169 * ```
170 *
171 * For i.e. Prisms and Pyramids in 3d this does not work because they contain squares and triangles as faces.
172 * And depending on the dynamic face index a different face geometry is needed.
173 *
174 */
175 TestFaceBasis *faceBasis = Impl::toGeometryTypeIdConstant<dimension-1>(refElement.type( face, 1 ), [&](auto faceGeometryTypeId) {
176 return TestFaceBasisFactory::template create< decltype(faceGeometryTypeId)::value >( order );
177 });
178 faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal( face ) );
179 }
180 assert( faceStructure_.size() == faceSize_ );
181 }
182
183 private:
184 struct FaceStructure
185 {
186 FaceStructure( TestFaceBasis *tfb, const Normal &n )
187 : basis_( tfb ), normal_( &n )
188 {}
189
190 TestFaceBasis *basis_;
191 const Dune::FieldVector< Field, dimension > *normal_;
192 };
193
194 std::vector< FaceStructure > faceStructure_;
195 TestBasis *testBasis_ = nullptr;
196 GeometryType geometry_;
197 unsigned int faceSize_;
198 std::size_t order_;
199 };
200
201
202
203 // RaviartThomasL2Interpolation
204 // ----------------------------
205
211 template< unsigned int dimension, class F>
213 : public InterpolationHelper< F ,dimension >
214 {
217
218 public:
219 typedef F Field;
222 : order_(0),
223 size_(0)
224 {}
225
226 template< class Function, class Vector >
227 auto interpolate ( const Function &function, Vector &coefficients ) const
228 -> std::enable_if_t< std::is_same< decltype(std::declval<Vector>().resize(1) ),void >::value,void>
229 {
230 coefficients.resize(size());
231 typename Base::template Helper<Function,Vector,true> func( function,coefficients );
232 interpolate(func);
233 }
234
235 template< class Basis, class Matrix >
236 auto interpolate ( const Basis &basis, Matrix &matrix ) const
237 -> std::enable_if_t< std::is_same<
238 decltype(std::declval<Matrix>().rowPtr(0)), typename Matrix::Field* >::value,void>
239 {
240 matrix.resize( size(), basis.size() );
241 typename Base::template Helper<Basis,Matrix,false> func( basis,matrix );
242 interpolate(func);
243 }
244
245 std::size_t order() const
246 {
247 return order_;
248 }
249 std::size_t size() const
250 {
251 return size_;
252 }
253 template <GeometryType::Id geometryId>
254 void build( std::size_t order )
255 {
256 size_ = 0;
257 order_ = order;
258 builder_.template build<geometryId>(order_);
259 if (builder_.testBasis())
260 size_ += dimension*builder_.testBasis()->size();
261 for ( unsigned int f=0; f<builder_.faceSize(); ++f )
262 if (builder_.testFaceBasis(f))
263 size_ += builder_.testFaceBasis(f)->size();
264 }
265
266 void setLocalKeys(std::vector< LocalKey > &keys) const
267 {
268 keys.resize(size());
269 unsigned int row = 0;
270 for (unsigned int f=0; f<builder_.faceSize(); ++f)
271 {
272 if (builder_.faceSize())
273 for (unsigned int i=0; i<builder_.testFaceBasis(f)->size(); ++i,++row)
274 keys[row] = LocalKey(f,1,i);
275 }
276 if (builder_.testBasis())
277 for (unsigned int i=0; i<builder_.testBasis()->size()*dimension; ++i,++row)
278 keys[row] = LocalKey(0,0,i);
279 assert( row == size() );
280 }
281
282 protected:
283 template< class Func, class Container, bool type >
284 void interpolate ( typename Base::template Helper<Func,Container,type> &func ) const
285 {
286 const Dune::GeometryType geoType = builder_.type();
287
288 std::vector< Field > testBasisVal;
289
290 for (unsigned int i=0; i<size(); ++i)
291 for (unsigned int j=0; j<func.size(); ++j)
292 func.set(i,j,0);
293
294 unsigned int row = 0;
295
296 // boundary dofs:
297 typedef Dune::QuadratureRule<Field, dimension-1> FaceQuadrature;
298 typedef Dune::QuadratureRules<Field, dimension-1> FaceQuadratureRules;
299
300 const auto &refElement = Dune::ReferenceElements< Field, dimension >::general( geoType );
301
302 for (unsigned int f=0; f<builder_.faceSize(); ++f)
303 {
304 if (!builder_.testFaceBasis(f))
305 continue;
306 testBasisVal.resize(builder_.testFaceBasis(f)->size());
307
308 const auto &geometry = refElement.template geometry< 1 >( f );
309 const Dune::GeometryType subGeoType( geometry.type().id(), dimension-1 );
310 const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*order_+2 );
311
312 const unsigned int quadratureSize = faceQuad.size();
313 for( unsigned int qi = 0; qi < quadratureSize; ++qi )
314 {
315 if (dimension>1)
316 builder_.testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal);
317 else
318 testBasisVal[0] = 1.;
319 fillBnd( row, testBasisVal,
320 func.evaluate( geometry.global( faceQuad[qi].position() ) ),
321 builder_.normal(f), faceQuad[qi].weight(),
322 func);
323 }
324
325 row += builder_.testFaceBasis(f)->size();
326 }
327 // element dofs
328 if (builder_.testBasis())
329 {
330 testBasisVal.resize(builder_.testBasis()->size());
331
332 typedef Dune::QuadratureRule<Field, dimension> Quadrature;
333 typedef Dune::QuadratureRules<Field, dimension> QuadratureRules;
334 const Quadrature &elemQuad = QuadratureRules::rule( geoType, 2*order_+1 );
335
336 const unsigned int quadratureSize = elemQuad.size();
337 for( unsigned int qi = 0; qi < quadratureSize; ++qi )
338 {
339 builder_.testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
340 fillInterior( row, testBasisVal,
341 func.evaluate(elemQuad[qi].position()),
342 elemQuad[qi].weight(),
343 func );
344 }
345
346 row += builder_.testBasis()->size()*dimension;
347 }
348 assert(row==size());
349 }
350
351 private:
361 template <class MVal, class RTVal,class Matrix>
362 void fillBnd (unsigned int startRow,
363 const MVal &mVal,
364 const RTVal &rtVal,
365 const FieldVector<Field,dimension> &normal,
366 const Field &weight,
367 Matrix &matrix) const
368 {
369 const unsigned int endRow = startRow+mVal.size();
370 typename RTVal::const_iterator rtiter = rtVal.begin();
371 for ( unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
372 {
373 Field cFactor = (*rtiter)*normal;
374 typename MVal::const_iterator miter = mVal.begin();
375 for (unsigned int row = startRow;
376 row!=endRow; ++miter, ++row )
377 {
378 matrix.add(row,col, (weight*cFactor)*(*miter) );
379 }
380 assert( miter == mVal.end() );
381 }
382 }
391 template <class MVal, class RTVal,class Matrix>
392 void fillInterior (unsigned int startRow,
393 const MVal &mVal,
394 const RTVal &rtVal,
395 Field weight,
396 Matrix &matrix) const
397 {
398 const unsigned int endRow = startRow+mVal.size()*dimension;
399 typename RTVal::const_iterator rtiter = rtVal.begin();
400 for ( unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
401 {
402 typename MVal::const_iterator miter = mVal.begin();
403 for (unsigned int row = startRow;
404 row!=endRow; ++miter,row+=dimension )
405 {
406 for (unsigned int i=0; i<dimension; ++i)
407 {
408 matrix.add(row+i,col, (weight*(*miter))*(*rtiter)[i] );
409 }
410 }
411 assert( miter == mVal.end() );
412 }
413 }
414
415 Builder builder_;
416 std::size_t order_;
417 std::size_t size_;
418 };
419
420 template < unsigned int dim, class Field >
422 {
425 typedef std::size_t Key;
426 typedef typename std::remove_const<Object>::type NonConstObject;
427
428 template <GeometryType::Id geometryId>
429 static Object *create( const Key &key )
430 {
431 if ( !supports<geometryId>(key) )
432 return 0;
433 NonConstObject *interpol = new NonConstObject();
434 interpol->template build<geometryId>(key);
435 return interpol;
436 }
437 template< GeometryType::Id geometryId >
438 static bool supports ( const Key &key )
439 {
440 return GeometryType(geometryId).isSimplex();
441 }
442 static void release( Object *object ) { delete object; }
443 };
444
445} // namespace Dune
446
447#endif // #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
Definition: bdfmcube.hh:16
Describe position of one degree of freedom.
Definition: localkey.hh:21
Definition: nedelecsimplexinterpolation.hh:36
LocalCoefficientsContainer(const Setter &setter)
Definition: nedelecsimplexinterpolation.hh:41
const LocalKey & localKey(const unsigned int i) const
Definition: raviartthomassimplexinterpolation.hh:46
std::size_t size() const
Definition: nedelecsimplexinterpolation.hh:52
Definition: orthonormalbasis.hh:18
static void release(Object *object)
Definition: orthonormalbasis.hh:55
Definition: raviartthomassimplexinterpolation.hh:422
std::remove_const< Object >::type NonConstObject
Definition: raviartthomassimplexinterpolation.hh:426
static void release(Object *object)
Definition: raviartthomassimplexinterpolation.hh:442
static bool supports(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:438
static Object * create(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:429
RTL2InterpolationBuilder< dim, Field > Builder
Definition: raviartthomassimplexinterpolation.hh:423
const RaviartThomasL2Interpolation< dim, Field > Object
Definition: raviartthomassimplexinterpolation.hh:424
std::size_t Key
Definition: raviartthomassimplexinterpolation.hh:425
Definition: raviartthomassimplexinterpolation.hh:68
std::size_t Key
Definition: raviartthomassimplexinterpolation.hh:69
static void release(Object *object)
Definition: raviartthomassimplexinterpolation.hh:89
static bool supports(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:85
const LocalCoefficientsContainer Object
Definition: raviartthomassimplexinterpolation.hh:70
static Object * create(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:73
Definition: raviartthomassimplexinterpolation.hh:105
TestBasis * testBasis() const
Definition: raviartthomassimplexinterpolation.hh:142
FieldVector< Field, dimension > Normal
Definition: raviartthomassimplexinterpolation.hh:117
TestBasisFactory::Object TestBasis
Definition: raviartthomassimplexinterpolation.hh:110
TestFaceBasisFactory::Object TestFaceBasis
Definition: raviartthomassimplexinterpolation.hh:114
unsigned int faceSize() const
Definition: raviartthomassimplexinterpolation.hh:139
void build(std::size_t order)
Definition: raviartthomassimplexinterpolation.hh:151
TestFaceBasis * testFaceBasis(unsigned int f) const
Definition: raviartthomassimplexinterpolation.hh:145
GeometryType type() const
Definition: raviartthomassimplexinterpolation.hh:134
const Normal & normal(unsigned int f) const
Definition: raviartthomassimplexinterpolation.hh:148
RTL2InterpolationBuilder(const RTL2InterpolationBuilder &)=delete
OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory
Definition: raviartthomassimplexinterpolation.hh:113
RTL2InterpolationBuilder(RTL2InterpolationBuilder &&)=delete
OrthonormalBasisFactory< dimension, Field > TestBasisFactory
Definition: raviartthomassimplexinterpolation.hh:109
std::size_t order() const
Definition: raviartthomassimplexinterpolation.hh:136
static const unsigned int dimension
Definition: raviartthomassimplexinterpolation.hh:106
unsigned int topologyId() const
Definition: raviartthomassimplexinterpolation.hh:132
~RTL2InterpolationBuilder()
Definition: raviartthomassimplexinterpolation.hh:124
An L2-based interpolation for Raviart Thomas.
Definition: raviartthomassimplexinterpolation.hh:214
std::size_t order() const
Definition: raviartthomassimplexinterpolation.hh:245
RaviartThomasL2Interpolation()
Definition: raviartthomassimplexinterpolation.hh:221
void interpolate(typename Base::template Helper< Func, Container, type > &func) const
Definition: raviartthomassimplexinterpolation.hh:284
auto interpolate(const Basis &basis, Matrix &matrix) const -> std::enable_if_t< std::is_same< decltype(std::declval< Matrix >().rowPtr(0)), typename Matrix::Field * >::value, void >
Definition: raviartthomassimplexinterpolation.hh:236
RTL2InterpolationBuilder< dimension, Field > Builder
Definition: raviartthomassimplexinterpolation.hh:220
F Field
Definition: raviartthomassimplexinterpolation.hh:219
void build(std::size_t order)
Definition: raviartthomassimplexinterpolation.hh:254
auto interpolate(const Function &function, Vector &coefficients) const -> std::enable_if_t< std::is_same< decltype(std::declval< Vector >().resize(1)), void >::value, void >
Definition: raviartthomassimplexinterpolation.hh:227
std::size_t size() const
Definition: raviartthomassimplexinterpolation.hh:249
void setLocalKeys(std::vector< LocalKey > &keys) const
Definition: raviartthomassimplexinterpolation.hh:266
Definition: interpolationhelper.hh:20
Definition: interpolationhelper.hh:22
Definition: polynomialbasis.hh:63
unsigned int size() const
Definition: polynomialbasis.hh:111