dune-localfunctions  2.8.0
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/common/concept.hh>
7 
8 #include <dune/geometry/quadraturerules.hh>
9 
12 
13 namespace Dune
14 {
30  template< class B, class Q, bool onb >
32 
33  template< class B, class Q >
35  {
37 
38  public:
39  typedef B Basis;
40  typedef Q Quadrature;
41 
42  static const unsigned int dimension = Basis::dimension;
43 
45  template< class Function, class DofField, std::enable_if_t<models<Impl::FunctionWithEvaluate<typename Function::DomainType, typename Function::RangeType>, Function>(), int> = 0 >
46  void interpolate ( const Function &function, std::vector< DofField > &coefficients ) const
47  {
48  typedef typename Quadrature::iterator Iterator;
49  typedef FieldVector< DofField, Basis::dimRange > RangeVector;
50 
51  const unsigned int size = basis().size();
52  static std::vector< RangeVector > basisValues( size );
53 
54  coefficients.resize( size );
55  basisValues.resize( size );
56  for( unsigned int i = 0; i < size; ++i )
57  coefficients[ i ] = Zero< DofField >();
58 
59  const Iterator end = quadrature().end();
60  for( Iterator it = quadrature().begin(); it != end; ++it )
61  {
62  basis().evaluate( it->position(), basisValues );
63  typename Function::RangeType val;
64  function.evaluate( field_cast<typename Function::DomainType::field_type>(it->position()), val );
65  RangeVector factor = field_cast< DofField >( val );
66  factor *= field_cast< DofField >( it->weight() );
67  for( unsigned int i = 0; i < size; ++i )
68  coefficients[ i ] += factor * basisValues[ i ];
69  }
70  }
71 
73  template< class Function, class DofField, std::enable_if_t<models<Impl::FunctionWithCallOperator<typename Quadrature::value_type::Vector>, Function>(), int> = 0 >
74  void interpolate ( const Function &function, std::vector< DofField > &coefficients ) const
75  {
76  typedef FieldVector< DofField, Basis::dimRange > RangeVector;
77 
78  const unsigned int size = basis().size();
79  static std::vector< RangeVector > basisValues( size );
80 
81  coefficients.resize( size );
82  basisValues.resize( size );
83  for( unsigned int i = 0; i < size; ++i )
84  coefficients[ i ] = Zero< DofField >();
85 
86  for (auto&& qp : quadrature())
87  {
88  basis().evaluate( qp.position(), basisValues );
89  auto val = function( qp.position() );
90  RangeVector factor = field_cast< DofField >( val );
91  factor *= field_cast< DofField >( qp.weight() );
92  for( unsigned int i = 0; i < size; ++i )
93  coefficients[ i ] += factor * basisValues[ i ];
94  }
95  }
96 
97  const Basis &basis () const
98  {
99  return basis_;
100  }
101 
102  const Quadrature &quadrature () const
103  {
104  return quadrature_;
105  }
106 
107  protected:
109  : basis_( basis ),
111  {}
112 
113  const Basis &basis_;
115  };
116 
117  template< class B, class Q >
118  struct LocalL2Interpolation<B,Q,true>
119  : public LocalL2InterpolationBase<B,Q>
120  {
122  template< class BasisFactory, bool onb >
124  using typename Base::Basis;
125  using typename Base::Quadrature;
126  private:
127  LocalL2Interpolation ( const typename Base::Basis &basis, const typename Base::Quadrature &quadrature )
128  : Base(basis,quadrature)
129  {}
130  };
131  template< class B, class Q >
132  struct LocalL2Interpolation<B,Q,false>
133  : public LocalL2InterpolationBase<B,Q>
134  {
136  template< class BasisFactory, bool onb >
138  using typename Base::Basis;
139  using typename Base::Quadrature;
140  template< class Function, class DofField >
141  void interpolate ( const Function &function, std::vector< DofField > &coefficients ) const
142  {
143  const unsigned size = Base::basis().size();
144  Base::interpolate(function,val_);
145  coefficients.resize( size );
146  for (unsigned int i=0; i<size; ++i)
147  {
148  coefficients[i] = 0;
149  for (unsigned int j=0; j<size; ++j)
150  {
151  coefficients[i] += field_cast<DofField>(massMatrix_(i,j)*val_[j]);
152  }
153  }
154  }
155  private:
156  LocalL2Interpolation ( const typename Base::Basis &basis, const typename Base::Quadrature &quadrature )
157  : Base(basis,quadrature),
158  val_(basis.size()),
159  massMatrix_()
160  {
161  typedef FieldVector< Field, Base::Basis::dimRange > RangeVector;
162  typedef typename Base::Quadrature::iterator Iterator;
163  const unsigned size = basis.size();
164  std::vector< RangeVector > basisValues( size );
165 
166  massMatrix_.resize( size,size );
167  for (unsigned int i=0; i<size; ++i)
168  for (unsigned int j=0; j<size; ++j)
169  massMatrix_(i,j) = 0;
170  const Iterator end = Base::quadrature().end();
171  for( Iterator it = Base::quadrature().begin(); it != end; ++it )
172  {
173  Base::basis().evaluate( it->position(), basisValues );
174  for (unsigned int i=0; i<size; ++i)
175  for (unsigned int j=0; j<size; ++j)
176  massMatrix_(i,j) += (basisValues[i]*basisValues[j])*it->weight();
177  }
178  if ( !massMatrix_.invert() )
179  {
180  DUNE_THROW(MathError, "Mass matrix singular in LocalL2Interpolation");
181  }
182 
183  }
184  typedef typename Base::Basis::StorageField Field;
185  typedef FieldVector< Field, Base::Basis::dimRange > RangeVector;
186  typedef LFEMatrix<Field> MassMatrix;
187  mutable std::vector<Field> val_;
188  MassMatrix massMatrix_;
189  };
190 
195  template< class BasisFactory, bool onb >
197  {
198  static const unsigned int dimension = BasisFactory::dimension;
199  typedef typename BasisFactory::Key Key;
200  typedef typename BasisFactory::Object Basis;
201  typedef double Field;
202  typedef QuadratureRule<Field,dimension> Quadrature;
203  typedef QuadratureRules<Field,dimension> QuadratureProvider;
205  typedef const LocalInterpolation Object;
206 
207  template< GeometryType::Id geometryId >
208  static Object *create ( const Key &key )
209  {
210  constexpr Dune::GeometryType geometry = geometryId;
211  const Basis *basis = BasisFactory::template create< geometry >( key );
212  const Quadrature & quadrature = QuadratureProvider::rule(geometry, 2*basis->order()+1);
213  return new Object( *basis, quadrature );
214  }
215  static void release ( Object *object )
216  {
217  const Basis &basis = object->basis();
218  BasisFactory::release( &basis );
219  delete object;
220  }
221  };
222 
223 }
224 
225 #endif // #ifndef DUNE_L2INTERPOLATION_HH
Definition: bdfmcube.hh:16
A class representing the zero of a given Field.
Definition: field.hh:77
A local L2 interpolation taking a test basis and a quadrature rule.
Definition: l2interpolation.hh:31
Definition: l2interpolation.hh:35
const Basis & basis() const
Definition: l2interpolation.hh:97
LocalL2InterpolationBase(const Basis &basis, const Quadrature &quadrature)
Definition: l2interpolation.hh:108
void interpolate(const Function &function, std::vector< DofField > &coefficients) const
Interpolate a function that implements void evaluate(Domain, Range&)
Definition: l2interpolation.hh:46
const Quadrature & quadrature_
Definition: l2interpolation.hh:114
const Basis & basis_
Definition: l2interpolation.hh:113
static const unsigned int dimension
Definition: l2interpolation.hh:42
Q Quadrature
Definition: l2interpolation.hh:40
const Quadrature & quadrature() const
Definition: l2interpolation.hh:102
B Basis
Definition: l2interpolation.hh:39
LocalL2InterpolationBase< B, Q > Base
Definition: l2interpolation.hh:121
void interpolate(const Function &function, std::vector< DofField > &coefficients) const
Definition: l2interpolation.hh:141
LocalL2InterpolationBase< B, Q > Base
Definition: l2interpolation.hh:135
A factory class for the local l2 interpolations taking a basis factory.
Definition: l2interpolation.hh:197
static const unsigned int dimension
Definition: l2interpolation.hh:198
static void release(Object *object)
Definition: l2interpolation.hh:215
BasisFactory::Object Basis
Definition: l2interpolation.hh:200
double Field
Definition: l2interpolation.hh:201
QuadratureRules< Field, dimension > QuadratureProvider
Definition: l2interpolation.hh:203
QuadratureRule< Field, dimension > Quadrature
Definition: l2interpolation.hh:202
static Object * create(const Key &key)
Definition: l2interpolation.hh:208
BasisFactory::Key Key
Definition: l2interpolation.hh:199
const LocalInterpolation Object
Definition: l2interpolation.hh:205
LocalL2Interpolation< Basis, Quadrature, onb > LocalInterpolation
Definition: l2interpolation.hh:204