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