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