dune-localfunctions  2.2.1
polynomialbasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_POLYNOMIALBASIS_HH
2 #define DUNE_POLYNOMIALBASIS_HH
3 
4 #include <fstream>
5 
6 #include <dune/common/fmatrix.hh>
7 
9 
14 
15 namespace Dune
16 {
17 
18  // PolynomialBasis
19  // ---------------
20 
58  template< class Eval, class CM, class D=double, class R=double >
60  {
62  typedef Eval Evaluator;
63 
64  public:
65  typedef CM CoefficientMatrix;
66 
67  typedef typename CoefficientMatrix::Field StorageField;
68 
69  static const unsigned int dimension = Evaluator::dimension;
70  static const unsigned int dimRange = Evaluator::dimRange*CoefficientMatrix::blockSize;
72  R,dimRange,FieldVector<R,dimRange>,
73  FieldMatrix<R,dimRange,dimension> > Traits;
74  typedef typename Evaluator::Basis Basis;
75  typedef typename Evaluator::DomainVector DomainVector;
76 
78  const CoefficientMatrix &coeffMatrix,
79  unsigned int size)
80  : basis_(basis),
81  coeffMatrix_(&coeffMatrix),
82  eval_(basis),
83  order_(basis.order()),
84  size_(size)
85  {
86  // assert(coeffMatrix_);
87  // assert(size_ <= coeffMatrix.size()); // !!!
88  }
89 
90  const Basis &basis () const
91  {
92  return basis_;
93  }
94 
95  const CoefficientMatrix &matrix () const
96  {
97  return *coeffMatrix_;
98  }
99 
100  const unsigned int order () const
101  {
102  return order_;
103  }
104 
105  const unsigned int size () const
106  {
107  return size_;
108  }
109 
111  void evaluateFunction (const typename Traits::DomainType& x,
112  std::vector<typename Traits::RangeType>& out) const
113  {
114  out.resize(size());
115  evaluate(x,out);
116  }
117 
119  void evaluateJacobian (const typename Traits::DomainType& x, // position
120  std::vector<typename Traits::JacobianType>& out) const // return value
121  {
122  out.resize(size());
123  jacobian(x,out);
124  }
125 
126  template< unsigned int deriv, class F >
127  void evaluate ( const DomainVector &x, F *values ) const
128  {
129  coeffMatrix_->mult( eval_.template evaluate<deriv>( x ), size(), values);
130  }
131  template< unsigned int deriv, class DVector, class F >
132  void evaluate ( const DVector &x, F *values ) const
133  {
134  assert( DVector::dimension == dimension);
135  DomainVector bx;
136  for( int d = 0; d < dimension; ++d )
137  field_cast( x[ d ], bx[ d ] );
138  evaluate<deriv>( bx, values );
139  }
140 
141  template <bool dummy,class DVector>
142  struct Convert
143  {
144  static DomainVector apply( const DVector &x )
145  {
146  assert( DVector::dimension == dimension);
147  DomainVector bx;
148  for( unsigned int d = 0; d < dimension; ++d )
149  field_cast( x[ d ], bx[ d ] );
150  return bx;
151  }
152  };
153  template <bool dummy>
154  struct Convert<dummy,DomainVector>
155  {
156  static const DomainVector &apply( const DomainVector &x )
157  {
158  return x;
159  }
160  };
161  template< unsigned int deriv, class DVector, class RVector >
162  void evaluate ( const DVector &x, RVector &values ) const
163  {
164  assert(values.size()>=size());
166  coeffMatrix_->mult( eval_.template evaluate<deriv>( bx ), values );
167  }
168 
169  template <class Fy>
170  void evaluate ( const DomainVector &x, std::vector<FieldVector<Fy,dimRange> > &values ) const
171  {
172  evaluate<0>(x,values);
173  }
174  template< class DVector, class RVector >
175  void evaluate ( const DVector &x, RVector &values ) const
176  {
177  assert( DVector::dimension == dimension);
178  DomainVector bx;
179  for( unsigned int d = 0; d < dimension; ++d )
180  field_cast( x[ d ], bx[ d ] );
181  evaluate<0>( bx, values );
182  }
183 
184  template< unsigned int deriv, class Vector >
185  void evaluateSingle ( const DomainVector &x, Vector &values ) const
186  {
187  assert(values.size()>=size());
188  coeffMatrix_->template mult<deriv>( eval_.template evaluate<deriv>( x ), values );
189  }
190  template< unsigned int deriv, class Fy >
191  void evaluateSingle ( const DomainVector &x,
192  std::vector< FieldVector<FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size>,dimRange> > &values) const
193  {
194  evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
195  }
196  template< unsigned int deriv, class Fy >
197  void evaluateSingle ( const DomainVector &x,
198  std::vector< FieldVector<LFETensor<Fy,dimension,deriv>,dimRange> > &values) const
199  {
200  evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
201  }
202 
203  template <class Fy>
204  void jacobian ( const DomainVector &x, std::vector<FieldMatrix<Fy,dimRange,dimension> > &values ) const
205  {
206  assert(values.size()>=size());
207  evaluateSingle<1>(x,reinterpret_cast<std::vector<FieldVector<Fy,dimRange*dimension> >&>(values));
208  }
209  template< class DVector, class RVector >
210  void jacobian ( const DVector &x, RVector &values ) const
211  {
212  assert( DVector::dimension == dimension);
213  DomainVector bx;
214  for( unsigned int d = 0; d < dimension; ++d )
215  field_cast( x[ d ], bx[ d ] );
216  jacobian( bx, values );
217  }
218 
219  template <class Fy>
220  void integrate ( std::vector<Fy> &values ) const
221  {
222  assert(values.size()>=size());
223  coeffMatrix_->mult( eval_.template integrate(), values );
224  }
225 
226  protected:
228  : basis_(other.basis_),
229  coeffMatrix_(other.coeffMatrix_),
230  eval_(basis_),
231  order_(basis_.order()),
232  size_(other.size_)
233  {
234  }
236  const Basis &basis_;
238  mutable Evaluator eval_;
239  unsigned int order_,size_;
240  };
241 
248  template< class Eval, class CM = SparseCoeffMatrix<typename Eval::Field,Eval::dimRange>,
249  class D=double, class R=double>
251  : public PolynomialBasis< Eval, CM, D, R >
252  {
253  public:
254  typedef CM CoefficientMatrix;
255 
256  private:
257  typedef Eval Evaluator;
258 
261 
262  public:
263  typedef typename Base::Basis Basis;
264 
266  : Base(basis,coeffMatrix_,0)
267  {}
268 
269  template <class Matrix>
270  void fill(const Matrix& matrix)
271  {
272  coeffMatrix_.fill(matrix);
273  this->size_ = coeffMatrix_.size();
274  }
275  template <class Matrix>
276  void fill(const Matrix& matrix,int size)
277  {
278  coeffMatrix_.fill(matrix);
279  assert(size<=coeffMatrix_.size());
280  this->size_ = size;
281  }
282 
283  private:
286  CoefficientMatrix coeffMatrix_;
287  };
288 }
289 #endif // DUNE_POLYNOMIALBASIS_HH
290