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