dune-localfunctions  2.2.0
raviartthomas/interpolation.hh
Go to the documentation of this file.
1 #ifndef DUNE_RAVIARTTHOMASINTERPOLATION_HH
2 #define DUNE_RAVIARTTHOMASINTERPOLATION_HH
3 
4 #include <fstream>
5 #include <utility>
6 
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/forloop.hh>
9 
10 #include <dune/geometry/topologyfactory.hh>
11 #include <dune/geometry/referenceelements.hh>
12 #include <dune/geometry/genericgeometry/referenceelements.hh>
13 #include <dune/geometry/quadraturerules/gaussquadrature.hh>
14 
20 
21 namespace Dune
22 {
23 
24  // LocalCoefficientsContainer
25  // -------------------
26  template < unsigned int dim >
27  struct RaviartThomasCoefficientsFactory;
28  template < unsigned int dim, class Field >
29  struct RaviartThomasL2InterpolationFactory;
30 
32  {
34 
35  public:
36  template <class Setter>
37  LocalCoefficientsContainer ( const Setter &setter )
38  {
39  setter.setLocalKeys(localKey_);
40  }
41 
42  const LocalKey &localKey ( const unsigned int i ) const
43  {
44  assert( i < size() );
45  return localKey_[ i ];
46  }
47 
48  unsigned int size () const
49  {
50  return localKey_.size();
51  }
52 
53  private:
54  std::vector< LocalKey > localKey_;
55  };
56 
57  template < unsigned int dim >
59  {
60  static const unsigned int dimension = dim;
62  typedef unsigned int Key;
64  };
65  template < unsigned int dim >
67  public TopologyFactory< RaviartThomasCoefficientsFactoryTraits<dim> >
68  {
70  template <class Topology>
71  static typename Traits::Object *createObject( const typename Traits::Key &key )
72  {
73  typedef RaviartThomasL2InterpolationFactory<dim,double> InterpolationFactory;
74  if (! supports<Topology>(key) )
75  return 0;
76  typename InterpolationFactory::Object *interpol
77  = InterpolationFactory::template create<Topology>(key);
78  typename Traits::Object *localKeys = new typename Traits::Object(*interpol);
79  InterpolationFactory::release(interpol);
80  return localKeys;
81  }
82  template< class Topology >
83  static bool supports ( const typename Traits::Key &key )
84  {
86  }
87  };
88 
89  // LocalInterpolation
90  // -------------------
91 
92  // -----------------------------------------
93  // RTL2InterpolationBuilder
94  // -----------------------------------------
95  // L2 Interpolation requires:
96  // - for element
97  // - test basis
98  // - for each face (dynamic)
99  // - test basis
100  // - normal
101  template <unsigned int dim, class Field>
103  {
104  static const unsigned int dimension = dim;
109  typedef FieldVector<Field,dimension> Normal;
110 
112  {}
113 
115  {
116  TestBasisFactory::release(testBasis_);
117  for (unsigned int i=0;i<faceStructure_.size();++i)
118  TestFaceBasisFactory::release(faceStructure_[i].basis_);
119  }
120 
121  unsigned int topologyId() const
122  {
123  return topologyId_;
124  }
125  unsigned int order() const
126  {
127  return order_;
128  }
129  unsigned int faceSize() const
130  {
131  return faceSize_;
132  }
134  {
135  return testBasis_;
136  }
137  TestFaceBasis *testFaceBasis( unsigned int f ) const
138  {
139  assert( f < faceSize() );
140  return faceStructure_[f].basis_;
141  }
142  const Normal &normal( unsigned int f ) const
143  {
144  return *(faceStructure_[f].normal_);
145  }
146 
147  template <class Topology>
148  void build(unsigned int order)
149  {
150  order_ = order;
151  topologyId_ = Topology::id;
152  if (order>0)
153  testBasis_ = TestBasisFactory::template create<Topology>(order-1);
154  else
155  testBasis_ = 0;
156  const unsigned int size = GenericGeometry::Size<Topology,1>::value;
157  faceSize_ = size;
158  faceStructure_.reserve( faceSize_ );
159  ForLoop< Creator<Topology>::template GetCodim,0,size-1>::apply(order, faceStructure_ );
160  assert( faceStructure_.size() == faceSize_ );
161  }
162 
163  private:
164  struct FaceStructure
165  {
166  FaceStructure( TestFaceBasis *tfb,
167  const Normal &n )
168  : basis_(tfb), normal_(&n)
169  {
170  }
171  TestFaceBasis *basis_;
172  const Dune::FieldVector<Field,dimension> *normal_;
173  };
174  template < class Topology >
175  struct Creator
176  {
177  template < int face >
178  struct GetCodim
179  {
180  typedef typename GenericGeometry::SubTopology<Topology,1,face>::type FaceTopology;
181  static void apply( const unsigned int order,
182  std::vector<FaceStructure> &faceStructure )
183  {
184  faceStructure.push_back(
185  FaceStructure(
186  TestFaceBasisFactory::template create<FaceTopology>(order),
187  GenericGeometry::ReferenceElement<Topology,Field>::integrationOuterNormal(face)
188  ) );
189  }
190  };
191  };
192 
193  std::vector<FaceStructure> faceStructure_;
194  TestBasis *testBasis_;
195  unsigned int topologyId_, order_, faceSize_;
196  };
197 
198  // A L2 based interpolation for Raviart Thomas
199  // --------------------------------------------------
200  template< unsigned int dimension, class F>
202  : public InterpolationHelper<F,dimension>
203  {
206 
207  public:
208  typedef F Field;
211  : order_(0),
212  size_(0)
213  {
214  }
215 
216  template< class Function, class Fy >
217  void interpolate ( const Function &function, std::vector< Fy > &coefficients ) const
218  {
219  coefficients.resize(size());
220  typename Base::template Helper<Function,std::vector<Fy>,true> func( function,coefficients );
221  interpolate(func);
222  }
223  template< class Basis, class Matrix >
224  void interpolate ( const Basis &basis, Matrix &matrix ) const
225  {
226  matrix.resize( size(), basis.size() );
227  typename Base::template Helper<Basis,Matrix,false> func( basis,matrix );
228  interpolate(func);
229  }
230 
231  unsigned int order() const
232  {
233  return order_;
234  }
235  unsigned int size() const
236  {
237  return size_;
238  }
239  template <class Topology>
240  void build( unsigned int order )
241  {
242  size_ = 0;
243  order_ = order;
244  builder_.template build<Topology>(order_);
245  if (builder_.testBasis())
246  size_ += dimension*builder_.testBasis()->size();
247  for ( unsigned int f=0;f<builder_.faceSize();++f )
248  if (builder_.testFaceBasis(f))
249  size_ += builder_.testFaceBasis(f)->size();
250  }
251 
252  void setLocalKeys(std::vector< LocalKey > &keys) const
253  {
254  keys.resize(size());
255  unsigned int row = 0;
256  for (unsigned int f=0;f<builder_.faceSize();++f)
257  {
258  if (builder_.faceSize())
259  for (unsigned int i=0;i<builder_.testFaceBasis(f)->size();++i,++row)
260  keys[row] = LocalKey(f,1,i);
261  }
262  if (builder_.testBasis())
263  for (unsigned int i=0;i<builder_.testBasis()->size()*dimension;++i,++row)
264  keys[row] = LocalKey(0,0,i);
265  assert( row == size() );
266  }
267 
268  protected:
269  template< class Func, class Container, bool type >
270  void interpolate ( typename Base::template Helper<Func,Container,type> &func ) const
271  {
272  const Dune::GeometryType geoType( builder_.topologyId(), dimension );
273 
274  std::vector< Field > testBasisVal;
275 
276  for (unsigned int i=0;i<size();++i)
277  for (unsigned int j=0;j<func.size();++j)
278  func.set(i,j,0);
279 
280  unsigned int row = 0;
281 
282  // boundary dofs:
283  typedef Dune::GenericGeometry::GaussQuadratureProvider< dimension-1, Field >
284  FaceQuadratureProvider;
285 
286  typedef Dune::GenericReferenceElements< Field, dimension > RefElements;
287  typedef Dune::GenericReferenceElement< Field, dimension > RefElement;
288  typedef typename RefElement::template Codim< 1 >::Mapping Mapping;
289 
290  const RefElement &refElement = RefElements::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 Mapping &mapping = refElement.template mapping< 1 >( f );
299  const Dune::GeometryType subGeoType( mapping.type().id(), dimension-1 );
300  const typename FaceQuadratureProvider::Object *faceQuad = FaceQuadratureProvider::create( subGeoType, 2*order_+2 );
301 
302  const unsigned int quadratureSize = faceQuad->size();
303  for( unsigned int qi = 0; qi < quadratureSize; ++qi )
304  {
305  builder_.testFaceBasis(f)->template evaluate<0>(faceQuad->position(qi),testBasisVal);
306  fillBnd( row, testBasisVal,
307  func.evaluate( mapping.global( faceQuad->position(qi) ) ),
308  builder_.normal(f), faceQuad->weight(qi),
309  func);
310  }
311 
312  FaceQuadratureProvider::release( faceQuad );
313 
314  row += builder_.testFaceBasis(f)->size();
315  }
316  // element dofs
317  if (builder_.testBasis())
318  {
319  testBasisVal.resize(builder_.testBasis()->size());
320 
321  typedef Dune::GenericGeometry::GaussQuadratureProvider< dimension, Field > QuadratureProvider;
322  const typename QuadratureProvider::Object *elemQuad = QuadratureProvider::create( geoType, 2*order_+1 );
323 
324  const unsigned int quadratureSize = elemQuad->size();
325  for( unsigned int qi = 0; qi < quadratureSize; ++qi )
326  {
327  builder_.testBasis()->template evaluate<0>(elemQuad->position(qi),testBasisVal);
328  fillInterior( row, testBasisVal,
329  func.evaluate(elemQuad->position(qi)),
330  elemQuad->weight(qi),
331  func );
332  }
333 
334  QuadratureProvider::release( elemQuad );
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 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 #endif // DUNE_RAVIARTTHOMASINTERPOLATION_HH
427