dune-localfunctions  2.5.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/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;
125  typedef OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory;
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  typedef Dune::ReferenceElements< Field, dimension > RefElements;
291  typedef Dune::ReferenceElement< Field, dimension > RefElement;
292  typedef typename RefElement::template Codim< 1 >::Geometry Geometry;
293 
294  const RefElement &refElement = RefElements::general( geoType );
295 
296  for (unsigned int f=0; f<builder_.faceSize(); ++f)
297  {
298  if (!builder_.testFaceBasis(f))
299  continue;
300  testBasisVal.resize(builder_.testFaceBasis(f)->size());
301 
302  const Geometry &geometry = refElement.template geometry< 1 >( f );
303  const Dune::GeometryType subGeoType( geometry.type().id(), dimension-1 );
304  const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*order_+2 );
305 
306  const unsigned int quadratureSize = faceQuad.size();
307  for( unsigned int qi = 0; qi < quadratureSize; ++qi )
308  {
309  if (dimension>1)
310  builder_.testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal);
311  else
312  testBasisVal[0] = 1.;
313  fillBnd( row, testBasisVal,
314  func.evaluate( geometry.global( faceQuad[qi].position() ) ),
315  builder_.normal(f), faceQuad[qi].weight(),
316  func);
317  }
318 
319  row += builder_.testFaceBasis(f)->size();
320  }
321  // element dofs
322  if (builder_.testBasis())
323  {
324  testBasisVal.resize(builder_.testBasis()->size());
325 
326  typedef Dune::QuadratureRule<Field, dimension> Quadrature;
327  typedef Dune::QuadratureRules<Field, dimension> QuadratureRules;
328  const Quadrature &elemQuad = QuadratureRules::rule( geoType, 2*order_+1 );
329 
330  const unsigned int quadratureSize = elemQuad.size();
331  for( unsigned int qi = 0; qi < quadratureSize; ++qi )
332  {
333  builder_.testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
334  fillInterior( row, testBasisVal,
335  func.evaluate(elemQuad[qi].position()),
336  elemQuad[qi].weight(),
337  func );
338  }
339 
340  row += builder_.testBasis()->size()*dimension;
341  }
342  assert(row==size());
343  }
344 
345  private:
347  template <class MVal, class RTVal,class Matrix>
348  void fillBnd (unsigned int startRow,
349  const MVal &mVal,
350  const RTVal &rtVal,
351  const FieldVector<Field,dimension> &normal,
352  const Field &weight,
353  Matrix &matrix) const
354  {
355  const unsigned int endRow = startRow+mVal.size();
356  typename RTVal::const_iterator rtiter = rtVal.begin();
357  for ( unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
358  {
359  Field cFactor = (*rtiter)*normal;
360  typename MVal::const_iterator miter = mVal.begin();
361  for (unsigned int row = startRow;
362  row!=endRow; ++miter, ++row )
363  {
364  matrix.add(row,col, (weight*cFactor)*(*miter) );
365  }
366  assert( miter == mVal.end() );
367  }
368  }
369  template <class MVal, class RTVal,class Matrix>
370  void fillInterior (unsigned int startRow,
371  const MVal &mVal,
372  const RTVal &rtVal,
373  Field weight,
374  Matrix &matrix) const
375  {
376  const unsigned int endRow = startRow+mVal.size()*dimension;
377  typename RTVal::const_iterator rtiter = rtVal.begin();
378  for ( unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
379  {
380  typename MVal::const_iterator miter = mVal.begin();
381  for (unsigned int row = startRow;
382  row!=endRow; ++miter,row+=dimension )
383  {
384  for (unsigned int i=0; i<dimension; ++i)
385  {
386  matrix.add(row+i,col, (weight*(*miter))*(*rtiter)[i] );
387  }
388  }
389  assert( miter == mVal.end() );
390  }
391  }
392 
393  Builder builder_;
394  unsigned int order_;
395  unsigned int size_;
396  };
397 
398  template < unsigned int dim, class F >
400  {
401  static const unsigned int dimension = dim;
402  typedef unsigned int Key;
405  };
406  template < unsigned int dim, class Field >
408  public TopologyFactory< RaviartThomasL2InterpolationFactoryTraits<dim,Field> >
409  {
412  typedef typename Traits::Object Object;
413  typedef typename std::remove_const<Object>::type NonConstObject;
414  template <class Topology>
415  static typename Traits::Object *createObject( const typename Traits::Key &key )
416  {
417  if ( !supports<Topology>(key) )
418  return 0;
419  NonConstObject *interpol = new NonConstObject();
420  interpol->template build<Topology>(key);
421  return interpol;
422  }
423  template< class Topology >
424  static bool supports ( const typename Traits::Key &key )
425  {
427  }
428  };
429 
430 } // namespace Dune
431 
432 #endif // #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
Definition: tensor.hh:165
An L2-based interpolation for Raviart Thomas.
Definition: raviartthomassimplexinterpolation.hh:206
Definition: raviartthomassimplexinterpolation.hh:120
Definition: polynomialbasis.hh:62
const LocalKey & localKey(const unsigned int i) const
Definition: raviartthomassimplexinterpolation.hh:49
static bool supports(const typename Traits::Key &key)
Definition: raviartthomassimplexinterpolation.hh:102
unsigned int Key
Definition: raviartthomassimplexinterpolation.hh:402
LocalCoefficientsContainer(const Setter &setter)
Definition: raviartthomassimplexinterpolation.hh:44
TestBasisFactory::Object TestBasis
Definition: raviartthomassimplexinterpolation.hh:124
RaviartThomasL2InterpolationFactory< dim, F > Factory
Definition: raviartthomassimplexinterpolation.hh:404
~RTL2InterpolationBuilder()
Definition: raviartthomassimplexinterpolation.hh:134
void interpolate(const Function &function, std::vector< Fy > &coefficients) const
Definition: raviartthomassimplexinterpolation.hh:221
TestFaceBasis * testFaceBasis(unsigned int f) const
Definition: raviartthomassimplexinterpolation.hh:150
void interpolate(const Basis &basis, Matrix &matrix) const
Definition: raviartthomassimplexinterpolation.hh:228
const Normal & normal(unsigned int f) const
Definition: raviartthomassimplexinterpolation.hh:152
RaviartThomasL2Interpolation()
Definition: raviartthomassimplexinterpolation.hh:215
void build(unsigned int order)
Definition: raviartthomassimplexinterpolation.hh:244
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory
Definition: raviartthomassimplexinterpolation.hh:125
Definition: interpolationhelper.hh:17
unsigned int faceSize() const
Definition: raviartthomassimplexinterpolation.hh:147
void setLocalKeys(std::vector< LocalKey > &keys) const
Definition: raviartthomassimplexinterpolation.hh:256
unsigned int order() const
Definition: raviartthomassimplexinterpolation.hh:235
unsigned int Key
Definition: raviartthomassimplexinterpolation.hh:74
Traits::Object Object
Definition: raviartthomassimplexinterpolation.hh:412
static Traits::Object * createObject(const typename Traits::Key &key)
Definition: raviartthomassimplexinterpolation.hh:415
std::remove_const< Object >::type NonConstObject
Definition: raviartthomassimplexinterpolation.hh:413
FieldVector< Field, dimension > Normal
Definition: raviartthomassimplexinterpolation.hh:127
RTL2InterpolationBuilder< dim, Field > Builder
Definition: raviartthomassimplexinterpolation.hh:411
Definition: raviartthomassimplexinterpolation.hh:31
TestBasis * testBasis() const
Definition: raviartthomassimplexinterpolation.hh:149
unsigned int order() const
Definition: raviartthomassimplexinterpolation.hh:145
const LocalCoefficientsContainer Object
Definition: raviartthomassimplexinterpolation.hh:73
static bool supports(const typename Traits::Key &key)
Definition: raviartthomassimplexinterpolation.hh:424
GeometryType type() const
Definition: raviartthomassimplexinterpolation.hh:143
Definition: interpolationhelper.hh:20
OrthonormalBasisFactory< dimension, Field > TestBasisFactory
Definition: raviartthomassimplexinterpolation.hh:123
RaviartThomasCoefficientsFactoryTraits< dim > Traits
Definition: raviartthomassimplexinterpolation.hh:87
unsigned int size() const
Definition: raviartthomassimplexinterpolation.hh:239
unsigned int size() const
Definition: raviartthomassimplexinterpolation.hh:55
Definition: raviartthomassimplexinterpolation.hh:70
RTL2InterpolationBuilder< dimension, Field > Builder
Definition: raviartthomassimplexinterpolation.hh:214
Definition: raviartthomassimplexinterpolation.hh:28
void build(unsigned int order)
Definition: raviartthomassimplexinterpolation.hh:155
F Field
Definition: raviartthomassimplexinterpolation.hh:213
Definition: raviartthomassimplexinterpolation.hh:399
Definition: raviartthomassimplexinterpolation.hh:38
TestFaceBasisFactory::Object TestFaceBasis
Definition: raviartthomassimplexinterpolation.hh:126
static Traits::Object * createObject(const typename Traits::Key &key)
Definition: raviartthomassimplexinterpolation.hh:90
void interpolate(typename Base::template Helper< Func, Container, type > &func) const
Definition: raviartthomassimplexinterpolation.hh:274
Definition: orthonormalbasis.hh:19
RaviartThomasCoefficientsFactory< dim > Factory
Definition: raviartthomassimplexinterpolation.hh:75
Describe position of one degree of freedom.
Definition: localkey.hh:20
RaviartThomasL2InterpolationFactoryTraits< dim, Field > Traits
Definition: raviartthomassimplexinterpolation.hh:410
const RaviartThomasL2Interpolation< dim, F > Object
Definition: raviartthomassimplexinterpolation.hh:403
unsigned int topologyId() const
Definition: raviartthomassimplexinterpolation.hh:141