dune-localfunctions  2.3.1
l2interpolation.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_L2INTERPOLATION_HH
4 #define DUNE_L2INTERPOLATION_HH
5 
6 #include <dune/geometry/topologyfactory.hh>
7 #include <dune/geometry/quadraturerules.hh>
8 
10 
11 namespace Dune
12 {
28  template< class B, class Q, bool onb >
29  struct LocalL2Interpolation;
30 
31  template< class B, class Q >
33  {
35 
36  public:
37  typedef B Basis;
38  typedef Q Quadrature;
39 
40  static const unsigned int dimension = Basis::dimension;
41 
42  template< class Function, class DofField >
43  void interpolate ( const Function &function, std::vector< DofField > &coefficients ) const
44  {
45  typedef typename Quadrature::iterator Iterator;
46  typedef FieldVector< DofField, Basis::dimRange > RangeVector;
47 
48  const unsigned int size = basis().size();
49  static std::vector< RangeVector > basisValues( size );
50 
51  coefficients.resize( size );
52  basisValues.resize( size );
53  for( unsigned int i = 0; i < size; ++i )
54  coefficients[ i ] = Zero< DofField >();
55 
56  const Iterator end = quadrature().end();
57  for( Iterator it = quadrature().begin(); it != end; ++it )
58  {
59  basis().evaluate( it->position(), basisValues );
60  typename Function::RangeType val;
61  function.evaluate( field_cast<typename Function::DomainType::field_type>(it->position()), val );
62  RangeVector factor = field_cast< DofField >( val );
63  factor *= field_cast< DofField >( it->weight() );
64  for( unsigned int i = 0; i < size; ++i )
65  coefficients[ i ] += factor * basisValues[ i ];
66  }
67  }
68 
69  const Basis &basis () const
70  {
71  return basis_;
72  }
73 
74  const Quadrature &quadrature () const
75  {
76  return quadrature_;
77  }
78 
79  protected:
81  : basis_( basis ),
82  quadrature_( quadrature )
83  {}
84 
85  const Basis &basis_;
87  };
88 
89  template< class B, class Q >
90  struct LocalL2Interpolation<B,Q,true>
91  : public LocalL2InterpolationBase<B,Q>
92  {
94  template< class BasisFactory, bool onb >
96  using typename Base::Basis;
97  using typename Base::Quadrature;
98  private:
99  LocalL2Interpolation ( const typename Base::Basis &basis, const typename Base::Quadrature &quadrature )
100  : Base(basis,quadrature)
101  {}
102  };
103  template< class B, class Q >
104  struct LocalL2Interpolation<B,Q,false>
105  : public LocalL2InterpolationBase<B,Q>
106  {
108  template< class BasisFactory, bool onb >
110  using typename Base::Basis;
111  using typename Base::Quadrature;
112  template< class Function, class DofField >
113  void interpolate ( const Function &function, std::vector< DofField > &coefficients ) const
114  {
115  const unsigned size = Base::basis().size();
116  Base::interpolate(function,val_);
117  coefficients.resize( size );
118  for (unsigned int i=0; i<size; ++i)
119  {
120  coefficients[i] = 0;
121  for (unsigned int j=0; j<size; ++j)
122  {
123  coefficients[i] += field_cast<DofField>(massMatrix_(i,j)*val_[j]);
124  }
125  }
126  }
127  private:
128  LocalL2Interpolation ( const typename Base::Basis &basis, const typename Base::Quadrature &quadrature )
129  : Base(basis,quadrature),
130  val_(basis.size()),
131  massMatrix_()
132  {
133  typedef FieldVector< Field, Base::Basis::dimRange > RangeVector;
134  typedef typename Base::Quadrature::iterator Iterator;
135  const unsigned size = basis.size();
136  std::vector< RangeVector > basisValues( size );
137 
138  massMatrix_.resize( size,size );
139  for (unsigned int i=0; i<size; ++i)
140  for (unsigned int j=0; j<size; ++j)
141  massMatrix_(i,j) = 0;
142  const Iterator end = Base::quadrature().end();
143  for( Iterator it = Base::quadrature().begin(); it != end; ++it )
144  {
145  Base::basis().evaluate( it->position(), basisValues );
146  for (unsigned int i=0; i<size; ++i)
147  for (unsigned int j=0; j<size; ++j)
148  massMatrix_(i,j) += (basisValues[i]*basisValues[j])*it->weight();
149  }
150  if ( !massMatrix_.invert() )
151  {
152  DUNE_THROW(MathError, "Mass matrix singular in LocalL2Interpolation");
153  }
154 
155  }
156  typedef typename Base::Basis::StorageField Field;
157  typedef FieldVector< Field, Base::Basis::dimRange > RangeVector;
158  typedef LFEMatrix<Field> MassMatrix;
159  mutable std::vector<Field> val_;
160  MassMatrix massMatrix_;
161  };
162 
168  template< class BasisFactory, bool onb >
169  struct LocalL2InterpolationFactory;
170  template< class BasisFactory, bool onb >
172  {
173  static const unsigned int dimension = BasisFactory::dimension;
174  // typedef typename BasisFactory::StorageField Field;
175  typedef double Field;
176  typedef QuadratureRule<Field,dimension> Quadrature;
177  typedef QuadratureRules<Field,dimension> QuadratureProvider;
178 
179  typedef typename BasisFactory::Key Key;
180  typedef typename BasisFactory::Object Basis;
181  typedef LocalL2Interpolation< Basis, Quadrature, onb > LocalInterpolation;
182  typedef const LocalInterpolation Object;
184  };
185 
186  template< class BasisFactory, bool onb >
188  public TopologyFactory< LocalL2InterpolationFactoryTraits<BasisFactory,onb> >
189  {
191  static const unsigned int dimension = Traits::dimension;
192  typedef typename Traits::Key Key;
193  typedef typename Traits::Basis Basis;
194  typedef typename Traits::Object Object;
195  typedef typename Traits::Field Field;
196  typedef typename Traits::Quadrature Quadrature;
197 
198  template< class Topology >
199  static Object *createObject ( const Key &key )
200  {
201  Dune::GeometryType gt(Topology::id, Topology::dimension);
202  const Basis *basis = BasisFactory::template create< Topology >( key );
203  const Quadrature & quadrature = Traits::QuadratureProvider::rule(gt, 2*basis->order()+1);
204  return new Object( *basis, quadrature );
205  }
206  static void release ( Object *object )
207  {
208  const Basis &basis = object->basis();
209  BasisFactory::release( &basis );
210  delete object;
211  }
212  };
213 
214 }
215 
216 #endif // #ifndef DUNE_L2INTERPOLATION_HH