1 #ifndef DUNE_L2INTERPOLATION_HH
2 #define DUNE_L2INTERPOLATION_HH
4 #include <dune/geometry/topologyfactory.hh>
5 #include <dune/geometry/quadraturerules/gaussquadrature.hh>
26 template<
class B,
class Q,
bool onb >
29 template<
class B,
class Q >
38 static const unsigned int dimension = Basis::dimension;
40 template<
class Function,
class DofField >
41 void interpolate (
const Function &
function, std::vector< DofField > &coefficients )
const
43 typedef typename Quadrature::Iterator Iterator;
44 typedef FieldVector< DofField, Basis::dimRange > RangeVector;
46 const unsigned int size =
basis().size();
47 static std::vector< RangeVector > basisValues( size );
49 coefficients.resize( size );
50 basisValues.resize( size );
51 for(
unsigned int i = 0; i < size; ++i )
55 for( Iterator it =
quadrature().begin(); it != end; ++it )
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 ];
87 template<
class B,
class Q >
92 template<
class BasisFactory,
bool onb >
98 :
Base(basis,quadrature)
101 template<
class B,
class Q >
106 template<
class BasisFactory,
bool onb >
110 template<
class Function,
class DofField >
111 void interpolate (
const Function &
function, std::vector< DofField > &coefficients )
const
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)
119 for (
unsigned int j=0;j<size;++j)
121 coefficients[i] +=
field_cast<DofField>(massMatrix_(i,j)*val_[j]);
126 LocalL2Interpolation (
const typename Base::Basis &basis,
const typename Base::Quadrature &quadrature )
127 : Base(basis,quadrature),
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 );
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 )
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();
148 if ( !massMatrix_.invert() )
150 DUNE_THROW(MathError,
"Mass matrix singular in LocalL2Interpolation");
155 typedef FieldVector< Field, Base::Basis::dimRange > RangeVector;
156 typedef LFEMatrix<Field> MassMatrix;
157 mutable std::vector<Field> val_;
158 MassMatrix massMatrix_;
166 template<
class BasisFactory,
bool onb >
168 template<
class BasisFactory,
bool onb >
171 static const unsigned int dimension = BasisFactory::dimension;
177 typedef typename BasisFactory::Key
Key;
178 typedef typename BasisFactory::Object
Basis;
184 template<
class BasisFactory,
bool onb >
186 public TopologyFactory< LocalL2InterpolationFactoryTraits<BasisFactory,onb> >
196 template<
class Topology >
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 );
206 const Basis &basis =
object->basis();
207 const Quadrature &quadrature =
object->quadrature();
208 BasisFactory::release( &basis );
209 Traits::QuadratureProvider::release( &quadrature );
216 #endif // #ifndef DUNE_L2INTERPOLATION_HH