dune-localfunctions  2.8.0
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 
21 namespace 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_ )
128  TestFaceBasisFactory::release( f.basis_ );
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
@ value
Definition: tensor.hh:166
Describe position of one degree of freedom.
Definition: localkey.hh:21
Definition: nedelecsimplexinterpolation.hh:36
const LocalKey & localKey(const unsigned int i) const
Definition: raviartthomassimplexinterpolation.hh:46
LocalCoefficientsContainer(const Setter &setter)
Definition: nedelecsimplexinterpolation.hh:41
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 Object * create(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:429
static void release(Object *object)
Definition: raviartthomassimplexinterpolation.hh:442
static bool supports(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:438
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 Object * create(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:73
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
Definition: raviartthomassimplexinterpolation.hh:105
FieldVector< Field, dimension > Normal
Definition: raviartthomassimplexinterpolation.hh:117
TestBasis * testBasis() const
Definition: raviartthomassimplexinterpolation.hh:142
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
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
const Normal & normal(unsigned int f) const
Definition: raviartthomassimplexinterpolation.hh:148
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