dune-localfunctions  2.6-git
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/topologyfactory.hh>
12 #include <dune/geometry/quadraturerules.hh>
13 #include <dune/geometry/referenceelements.hh>
14 #include <dune/geometry/type.hh>
15 
20 
21 namespace Dune
22 {
23 
24  // Internal Forward Declarations
25  // -----------------------------
26 
27  template < unsigned int dim >
29 
30  template < unsigned int dim, class Field >
32 
33 
34 
35  // LocalCoefficientsContainer
36  // --------------------------
37 
39  {
41 
42  public:
43  template <class Setter>
44  LocalCoefficientsContainer ( const Setter &setter )
45  {
46  setter.setLocalKeys(localKey_);
47  }
48 
49  const LocalKey &localKey ( const unsigned int i ) const
50  {
51  assert( i < size() );
52  return localKey_[ i ];
53  }
54 
55  unsigned int size () const
56  {
57  return localKey_.size();
58  }
59 
60  private:
61  std::vector< LocalKey > localKey_;
62  };
63 
64 
65 
66  // RaviartThomasCoefficientsFactoryTraits
67  // --------------------------------------
68 
69  template < unsigned int dim >
71  {
72  static const unsigned int dimension = dim;
74  typedef unsigned int Key;
76  };
77 
78 
79 
80  // RaviartThomasCoefficientsFactory
81  // --------------------------------
82 
83  template < unsigned int dim >
85  : public TopologyFactory< RaviartThomasCoefficientsFactoryTraits< dim > >
86  {
88 
89  template< class Topology >
90  static typename Traits::Object *createObject( const typename Traits::Key &key )
91  {
92  typedef RaviartThomasL2InterpolationFactory< dim, double > InterpolationFactory;
93  if( !supports< Topology >( key ) )
94  return nullptr;
95  typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< Topology >( key );
96  typename Traits::Object *localKeys = new typename Traits::Object( *interpolation );
97  InterpolationFactory::release( interpolation );
98  return localKeys;
99  }
100 
101  template< class Topology >
102  static bool supports ( const typename Traits::Key &key )
103  {
105  }
106  };
107 
108 
109 
110  // RTL2InterpolationBuilder
111  // ------------------------
112 
113  // L2 Interpolation requires:
114  // - for element
115  // - test basis
116  // - for each face (dynamic)
117  // - test basis
118  // - normal
119  template< unsigned int dim, class Field >
121  {
122  static const unsigned int dimension = dim;
127  typedef FieldVector< Field, dimension > Normal;
128 
129  RTL2InterpolationBuilder () = default;
130 
133 
135  {
136  TestBasisFactory::release( testBasis_ );
137  for( FaceStructure &f : faceStructure_ )
138  TestFaceBasisFactory::release( f.basis_ );
139  }
140 
141  unsigned int topologyId () const { return topologyId_; }
142 
143  GeometryType type () const { return GeometryType( topologyId(), dimension ); }
144 
145  unsigned int order () const { return order_; }
146 
147  unsigned int faceSize () const { return faceSize_; }
148 
149  TestBasis *testBasis () const { return testBasis_; }
150  TestFaceBasis *testFaceBasis ( unsigned int f ) const { assert( f < faceSize() ); return faceStructure_[ f ].basis_; }
151 
152  const Normal &normal ( unsigned int f ) const { assert( f < faceSize() ); return *(faceStructure_[ f ].normal_); }
153 
154  template< class Topology >
155  void build ( unsigned int order )
156  {
157  order_ = order;
158  topologyId_ = Topology::id;
159 
160  testBasis_ = (order > 0 ? TestBasisFactory::template create< Topology >( order-1 ) : nullptr);
161 
162  const auto &refElement = ReferenceElements< Field, dimension >::general( type() );
163  faceSize_ = refElement.size( 1 );
164  faceStructure_.reserve( faceSize_ );
165  for( unsigned int face = 0; face < faceSize_; ++face )
166  {
167  TestFaceBasis *faceBasis = Impl::IfTopology< CreateFaceBasis, dimension-1 >::apply( refElement.type( face, 1 ).id(), order );
168  faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal( face ) );
169  }
170  assert( faceStructure_.size() == faceSize_ );
171  }
172 
173  private:
174  struct FaceStructure
175  {
176  FaceStructure( TestFaceBasis *tfb, const Normal &n )
177  : basis_( tfb ), normal_( &n )
178  {}
179 
180  TestFaceBasis *basis_;
181  const Dune::FieldVector< Field, dimension > *normal_;
182  };
183 
184  template< class FaceTopology >
185  struct CreateFaceBasis
186  {
187  static TestFaceBasis *apply ( unsigned int order ) { return TestFaceBasisFactory::template create< FaceTopology >( order ); }
188  };
189 
190  std::vector< FaceStructure > faceStructure_;
191  TestBasis *testBasis_ = nullptr;
192  unsigned int topologyId_, order_, faceSize_;
193  };
194 
195 
196 
197  // RaviartThomasL2Interpolation
198  // ----------------------------
199 
205  template< unsigned int dimension, class F>
207  : public InterpolationHelper< F ,dimension >
208  {
211 
212  public:
213  typedef F Field;
216  : order_(0),
217  size_(0)
218  {}
219 
220  template< class Function, class Fy >
221  void interpolate ( const Function &function, std::vector< Fy > &coefficients ) const
222  {
223  coefficients.resize(size());
224  typename Base::template Helper<Function,std::vector<Fy>,true> func( function,coefficients );
225  interpolate(func);
226  }
227  template< class Basis, class Matrix >
228  void interpolate ( const Basis &basis, Matrix &matrix ) const
229  {
230  matrix.resize( size(), basis.size() );
231  typename Base::template Helper<Basis,Matrix,false> func( basis,matrix );
232  interpolate(func);
233  }
234 
235  unsigned int order() const
236  {
237  return order_;
238  }
239  unsigned int size() const
240  {
241  return size_;
242  }
243  template <class Topology>
244  void build( unsigned int order )
245  {
246  size_ = 0;
247  order_ = order;
248  builder_.template build<Topology>(order_);
249  if (builder_.testBasis())
250  size_ += dimension*builder_.testBasis()->size();
251  for ( unsigned int f=0; f<builder_.faceSize(); ++f )
252  if (builder_.testFaceBasis(f))
253  size_ += builder_.testFaceBasis(f)->size();
254  }
255 
256  void setLocalKeys(std::vector< LocalKey > &keys) const
257  {
258  keys.resize(size());
259  unsigned int row = 0;
260  for (unsigned int f=0; f<builder_.faceSize(); ++f)
261  {
262  if (builder_.faceSize())
263  for (unsigned int i=0; i<builder_.testFaceBasis(f)->size(); ++i,++row)
264  keys[row] = LocalKey(f,1,i);
265  }
266  if (builder_.testBasis())
267  for (unsigned int i=0; i<builder_.testBasis()->size()*dimension; ++i,++row)
268  keys[row] = LocalKey(0,0,i);
269  assert( row == size() );
270  }
271 
272  protected:
273  template< class Func, class Container, bool type >
274  void interpolate ( typename Base::template Helper<Func,Container,type> &func ) const
275  {
276  const Dune::GeometryType geoType( builder_.topologyId(), dimension );
277 
278  std::vector< Field > testBasisVal;
279 
280  for (unsigned int i=0; i<size(); ++i)
281  for (unsigned int j=0; j<func.size(); ++j)
282  func.set(i,j,0);
283 
284  unsigned int row = 0;
285 
286  // boundary dofs:
287  typedef Dune::QuadratureRule<Field, dimension-1> FaceQuadrature;
288  typedef Dune::QuadratureRules<Field, dimension-1> FaceQuadratureRules;
289 
290  const auto &refElement = Dune::ReferenceElements< Field, dimension >::general( geoType );
291 
292  for (unsigned int f=0; f<builder_.faceSize(); ++f)
293  {
294  if (!builder_.testFaceBasis(f))
295  continue;
296  testBasisVal.resize(builder_.testFaceBasis(f)->size());
297 
298  const auto &geometry = refElement.template geometry< 1 >( f );
299  const Dune::GeometryType subGeoType( geometry.type().id(), dimension-1 );
300  const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*order_+2 );
301 
302  const unsigned int quadratureSize = faceQuad.size();
303  for( unsigned int qi = 0; qi < quadratureSize; ++qi )
304  {
305  if (dimension>1)
306  builder_.testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal);
307  else
308  testBasisVal[0] = 1.;
309  fillBnd( row, testBasisVal,
310  func.evaluate( geometry.global( faceQuad[qi].position() ) ),
311  builder_.normal(f), faceQuad[qi].weight(),
312  func);
313  }
314 
315  row += builder_.testFaceBasis(f)->size();
316  }
317  // element dofs
318  if (builder_.testBasis())
319  {
320  testBasisVal.resize(builder_.testBasis()->size());
321 
322  typedef Dune::QuadratureRule<Field, dimension> Quadrature;
323  typedef Dune::QuadratureRules<Field, dimension> QuadratureRules;
324  const Quadrature &elemQuad = QuadratureRules::rule( geoType, 2*order_+1 );
325 
326  const unsigned int quadratureSize = elemQuad.size();
327  for( unsigned int qi = 0; qi < quadratureSize; ++qi )
328  {
329  builder_.testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
330  fillInterior( row, testBasisVal,
331  func.evaluate(elemQuad[qi].position()),
332  elemQuad[qi].weight(),
333  func );
334  }
335 
336  row += builder_.testBasis()->size()*dimension;
337  }
338  assert(row==size());
339  }
340 
341  private:
343  template <class MVal, class RTVal,class Matrix>
344  void fillBnd (unsigned int startRow,
345  const MVal &mVal,
346  const RTVal &rtVal,
347  const FieldVector<Field,dimension> &normal,
348  const Field &weight,
349  Matrix &matrix) const
350  {
351  const unsigned int endRow = startRow+mVal.size();
352  typename RTVal::const_iterator rtiter = rtVal.begin();
353  for ( unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
354  {
355  Field cFactor = (*rtiter)*normal;
356  typename MVal::const_iterator miter = mVal.begin();
357  for (unsigned int row = startRow;
358  row!=endRow; ++miter, ++row )
359  {
360  matrix.add(row,col, (weight*cFactor)*(*miter) );
361  }
362  assert( miter == mVal.end() );
363  }
364  }
365  template <class MVal, class RTVal,class Matrix>
366  void fillInterior (unsigned int startRow,
367  const MVal &mVal,
368  const RTVal &rtVal,
369  Field weight,
370  Matrix &matrix) const
371  {
372  const unsigned int endRow = startRow+mVal.size()*dimension;
373  typename RTVal::const_iterator rtiter = rtVal.begin();
374  for ( unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
375  {
376  typename MVal::const_iterator miter = mVal.begin();
377  for (unsigned int row = startRow;
378  row!=endRow; ++miter,row+=dimension )
379  {
380  for (unsigned int i=0; i<dimension; ++i)
381  {
382  matrix.add(row+i,col, (weight*(*miter))*(*rtiter)[i] );
383  }
384  }
385  assert( miter == mVal.end() );
386  }
387  }
388 
389  Builder builder_;
390  unsigned int order_;
391  unsigned int size_;
392  };
393 
394  template < unsigned int dim, class F >
396  {
397  static const unsigned int dimension = dim;
398  typedef unsigned int Key;
401  };
402  template < unsigned int dim, class Field >
404  public TopologyFactory< RaviartThomasL2InterpolationFactoryTraits<dim,Field> >
405  {
408  typedef typename Traits::Object Object;
409  typedef typename std::remove_const<Object>::type NonConstObject;
410  template <class Topology>
411  static typename Traits::Object *createObject( const typename Traits::Key &key )
412  {
413  if ( !supports<Topology>(key) )
414  return 0;
415  NonConstObject *interpol = new NonConstObject();
416  interpol->template build<Topology>(key);
417  return interpol;
418  }
419  template< class Topology >
420  static bool supports ( const typename Traits::Key &key )
421  {
423  }
424  };
425 
426 } // namespace Dune
427 
428 #endif // #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
Dune::RaviartThomasL2InterpolationFactory::createObject
static Traits::Object * createObject(const typename Traits::Key &key)
Definition: raviartthomassimplexinterpolation.hh:411
Dune::RaviartThomasL2Interpolation::build
void build(unsigned int order)
Definition: raviartthomassimplexinterpolation.hh:244
Dune::RaviartThomasL2Interpolation::Builder
RTL2InterpolationBuilder< dimension, Field > Builder
Definition: raviartthomassimplexinterpolation.hh:214
Dune::InterpolationHelper::Helper< Basis, Matrix, false >
Definition: interpolationhelper.hh:73
Dune::RTL2InterpolationBuilder::TestBasisFactory
OrthonormalBasisFactory< dimension, Field > TestBasisFactory
Definition: raviartthomassimplexinterpolation.hh:123
Dune::RTL2InterpolationBuilder::~RTL2InterpolationBuilder
~RTL2InterpolationBuilder()
Definition: raviartthomassimplexinterpolation.hh:134
Dune::RaviartThomasL2InterpolationFactory::Traits
RaviartThomasL2InterpolationFactoryTraits< dim, Field > Traits
Definition: raviartthomassimplexinterpolation.hh:406
Dune::InterpolationHelper
Definition: interpolationhelper.hh:17
Dune::value
@ value
Definition: tensor.hh:165
polynomialbasis.hh
Dune::RTL2InterpolationBuilder::TestFaceBasis
TestFaceBasisFactory::Object TestFaceBasis
Definition: raviartthomassimplexinterpolation.hh:126
Dune::RTL2InterpolationBuilder::dimension
static const unsigned int dimension
Definition: raviartthomassimplexinterpolation.hh:122
Dune::RaviartThomasL2InterpolationFactoryTraits::Key
unsigned int Key
Definition: raviartthomassimplexinterpolation.hh:398
localkey.hh
Dune::RaviartThomasCoefficientsFactory::createObject
static Traits::Object * createObject(const typename Traits::Key &key)
Definition: raviartthomassimplexinterpolation.hh:90
Dune::RTL2InterpolationBuilder::order
unsigned int order() const
Definition: raviartthomassimplexinterpolation.hh:145
Dune::RTL2InterpolationBuilder::build
void build(unsigned int order)
Definition: raviartthomassimplexinterpolation.hh:155
Dune::RaviartThomasCoefficientsFactory
Definition: raviartthomassimplexinterpolation.hh:28
Dune::RaviartThomasL2Interpolation::interpolate
void interpolate(typename Base::template Helper< Func, Container, type > &func) const
Definition: raviartthomassimplexinterpolation.hh:274
Dune::OrthonormalBasisFactory
Definition: orthonormalbasis.hh:19
Dune::RTL2InterpolationBuilder::testFaceBasis
TestFaceBasis * testFaceBasis(unsigned int f) const
Definition: raviartthomassimplexinterpolation.hh:150
Dune::LocalCoefficientsContainer
Definition: raviartthomassimplexinterpolation.hh:38
Dune::RTL2InterpolationBuilder::faceSize
unsigned int faceSize() const
Definition: raviartthomassimplexinterpolation.hh:147
Dune::RaviartThomasL2Interpolation::size
unsigned int size() const
Definition: raviartthomassimplexinterpolation.hh:239
Dune::RaviartThomasL2Interpolation::RaviartThomasL2Interpolation
RaviartThomasL2Interpolation()
Definition: raviartthomassimplexinterpolation.hh:215
Dune::RaviartThomasL2InterpolationFactoryTraits
Definition: raviartthomassimplexinterpolation.hh:395
Dune::LocalCoefficientsContainer::size
unsigned int size() const
Definition: raviartthomassimplexinterpolation.hh:55
Dune::RaviartThomasL2Interpolation::interpolate
void interpolate(const Function &function, std::vector< Fy > &coefficients) const
Definition: raviartthomassimplexinterpolation.hh:221
Dune::RTL2InterpolationBuilder::TestBasis
TestBasisFactory::Object TestBasis
Definition: raviartthomassimplexinterpolation.hh:124
Dune::RTL2InterpolationBuilder::Normal
FieldVector< Field, dimension > Normal
Definition: raviartthomassimplexinterpolation.hh:127
Dune::RTL2InterpolationBuilder::type
GeometryType type() const
Definition: raviartthomassimplexinterpolation.hh:143
Dune::RaviartThomasL2InterpolationFactory::NonConstObject
std::remove_const< Object >::type NonConstObject
Definition: raviartthomassimplexinterpolation.hh:409
Dune::RaviartThomasCoefficientsFactoryTraits::Factory
RaviartThomasCoefficientsFactory< dim > Factory
Definition: raviartthomassimplexinterpolation.hh:75
Dune::RaviartThomasL2Interpolation::interpolate
void interpolate(const Basis &basis, Matrix &matrix) const
Definition: raviartthomassimplexinterpolation.hh:228
Dune::RaviartThomasCoefficientsFactory::supports
static bool supports(const typename Traits::Key &key)
Definition: raviartthomassimplexinterpolation.hh:102
Dune::RaviartThomasCoefficientsFactory::Traits
RaviartThomasCoefficientsFactoryTraits< dim > Traits
Definition: raviartthomassimplexinterpolation.hh:87
Dune::RaviartThomasL2InterpolationFactory::Object
Traits::Object Object
Definition: raviartthomassimplexinterpolation.hh:408
Dune::RaviartThomasL2Interpolation::setLocalKeys
void setLocalKeys(std::vector< LocalKey > &keys) const
Definition: raviartthomassimplexinterpolation.hh:256
Dune::RTL2InterpolationBuilder
Definition: raviartthomassimplexinterpolation.hh:120
Dune::RaviartThomasL2InterpolationFactory
Definition: raviartthomassimplexinterpolation.hh:31
Dune::RTL2InterpolationBuilder::testBasis
TestBasis * testBasis() const
Definition: raviartthomassimplexinterpolation.hh:149
Dune::RTL2InterpolationBuilder::RTL2InterpolationBuilder
RTL2InterpolationBuilder()=default
Dune::RaviartThomasL2InterpolationFactory::Builder
RTL2InterpolationBuilder< dim, Field > Builder
Definition: raviartthomassimplexinterpolation.hh:407
Dune::RaviartThomasL2InterpolationFactory::supports
static bool supports(const typename Traits::Key &key)
Definition: raviartthomassimplexinterpolation.hh:420
orthonormalbasis.hh
Dune::RaviartThomasL2InterpolationFactoryTraits::Object
const typedef RaviartThomasL2Interpolation< dim, F > Object
Definition: raviartthomassimplexinterpolation.hh:399
Dune::RaviartThomasCoefficientsFactoryTraits
Definition: raviartthomassimplexinterpolation.hh:70
Dune::RaviartThomasCoefficientsFactoryTraits::dimension
static const unsigned int dimension
Definition: raviartthomassimplexinterpolation.hh:72
Dune::InterpolationHelper::Helper
Definition: interpolationhelper.hh:20
Dune::RaviartThomasL2InterpolationFactoryTraits::Factory
RaviartThomasL2InterpolationFactory< dim, F > Factory
Definition: raviartthomassimplexinterpolation.hh:400
Dune::LocalCoefficientsContainer::localKey
const LocalKey & localKey(const unsigned int i) const
Definition: raviartthomassimplexinterpolation.hh:49
Dune::RaviartThomasL2Interpolation::order
unsigned int order() const
Definition: raviartthomassimplexinterpolation.hh:235
Dune::RTL2InterpolationBuilder::topologyId
unsigned int topologyId() const
Definition: raviartthomassimplexinterpolation.hh:141
Dune::RaviartThomasCoefficientsFactoryTraits::Object
const typedef LocalCoefficientsContainer Object
Definition: raviartthomassimplexinterpolation.hh:73
Dune::OrthonormalBasisFactory::Object
Traits::Object Object
Definition: orthonormalbasis.hh:45
Dune::LocalKey
Describe position of one degree of freedom.
Definition: localkey.hh:20
Dune::RTL2InterpolationBuilder::normal
const Normal & normal(unsigned int f) const
Definition: raviartthomassimplexinterpolation.hh:152
Dune::LocalCoefficientsContainer::LocalCoefficientsContainer
LocalCoefficientsContainer(const Setter &setter)
Definition: raviartthomassimplexinterpolation.hh:44
Dune::RaviartThomasCoefficientsFactoryTraits::Key
unsigned int Key
Definition: raviartthomassimplexinterpolation.hh:74
interpolationhelper.hh
Dune
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
Dune::RaviartThomasL2InterpolationFactoryTraits::dimension
static const unsigned int dimension
Definition: raviartthomassimplexinterpolation.hh:397
Dune::RTL2InterpolationBuilder::TestFaceBasisFactory
OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory
Definition: raviartthomassimplexinterpolation.hh:125
Dune::RaviartThomasL2Interpolation
An L2-based interpolation for Raviart Thomas.
Definition: raviartthomassimplexinterpolation.hh:206
Dune::RaviartThomasL2Interpolation::Field
F Field
Definition: raviartthomassimplexinterpolation.hh:213