dune-localfunctions  2.6-git
polynomialbasis.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_POLYNOMIALBASIS_HH
4 #define DUNE_POLYNOMIALBASIS_HH
5 
6 #include <fstream>
7 #include <numeric>
8 
9 #include <dune/common/fmatrix.hh>
10 
12 
17 
18 namespace Dune
19 {
20 
21  // PolynomialBasis
22  // ---------------
23 
61  template< class Eval, class CM, class D=double, class R=double >
63  {
65  typedef Eval Evaluator;
66 
67  public:
68  typedef CM CoefficientMatrix;
69 
70  typedef typename CoefficientMatrix::Field StorageField;
71 
72  static const unsigned int dimension = Evaluator::dimension;
73  static const unsigned int dimRange = Evaluator::dimRange*CoefficientMatrix::blockSize;
75  R,dimRange,FieldVector<R,dimRange>,
76  FieldMatrix<R,dimRange,dimension> > Traits;
77  typedef typename Evaluator::Basis Basis;
78  typedef typename Evaluator::DomainVector DomainVector;
79 
80  PolynomialBasis (const Basis &basis,
81  const CoefficientMatrix &coeffMatrix,
82  unsigned int size)
83  : basis_(basis),
84  coeffMatrix_(&coeffMatrix),
85  eval_(basis),
86  order_(basis.order()),
87  size_(size)
88  {
89  // assert(coeffMatrix_);
90  // assert(size_ <= coeffMatrix.size()); // !!!
91  }
92 
93  const Basis &basis () const
94  {
95  return basis_;
96  }
97 
98  const CoefficientMatrix &matrix () const
99  {
100  return *coeffMatrix_;
101  }
102 
103  unsigned int order () const
104  {
105  return order_;
106  }
107 
108  unsigned int size () const
109  {
110  return size_;
111  }
112 
114  void evaluateFunction (const typename Traits::DomainType& x,
115  std::vector<typename Traits::RangeType>& out) const
116  {
117  out.resize(size());
118  evaluate(x,out);
119  }
120 
122  void evaluateJacobian (const typename Traits::DomainType& x, // position
123  std::vector<typename Traits::JacobianType>& out) const // return value
124  {
125  out.resize(size());
126  jacobian(x,out);
127  }
128 
130  void partial (const std::array<unsigned int, dimension>& order,
131  const typename Traits::DomainType& in, // position
132  std::vector<typename Traits::RangeType>& out) const // return value
133  {
134  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
135  if (totalOrder == 0) {
136  evaluateFunction(in, out);
137  } else {
138  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
139  }
140  }
141 
142  template< unsigned int deriv, class F >
143  void evaluate ( const DomainVector &x, F *values ) const
144  {
145  coeffMatrix_->mult( eval_.template evaluate<deriv>( x ), size(), values);
146  }
147  template< unsigned int deriv, class DVector, class F >
148  void evaluate ( const DVector &x, F *values ) const
149  {
150  assert( DVector::dimension == dimension);
151  DomainVector bx;
152  for( int d = 0; d < dimension; ++d )
153  field_cast( x[ d ], bx[ d ] );
154  evaluate<deriv>( bx, values );
155  }
156 
157  template <bool dummy,class DVector>
158  struct Convert
159  {
160  static DomainVector apply( const DVector &x )
161  {
162  assert( DVector::dimension == dimension);
163  DomainVector bx;
164  for( unsigned int d = 0; d < dimension; ++d )
165  field_cast( x[ d ], bx[ d ] );
166  return bx;
167  }
168  };
169  template <bool dummy>
170  struct Convert<dummy,DomainVector>
171  {
172  static const DomainVector &apply( const DomainVector &x )
173  {
174  return x;
175  }
176  };
177  template< unsigned int deriv, class DVector, class RVector >
178  void evaluate ( const DVector &x, RVector &values ) const
179  {
180  assert(values.size()>=size());
181  const DomainVector &bx = Convert<true,DVector>::apply(x);
182  coeffMatrix_->mult( eval_.template evaluate<deriv>( bx ), values );
183  }
184 
185  template <class Fy>
186  void evaluate ( const DomainVector &x, std::vector<FieldVector<Fy,dimRange> > &values ) const
187  {
188  evaluate<0>(x,values);
189  }
190  template< class DVector, class RVector >
191  void evaluate ( const DVector &x, RVector &values ) const
192  {
193  assert( DVector::dimension == dimension);
194  DomainVector bx;
195  for( unsigned int d = 0; d < dimension; ++d )
196  field_cast( x[ d ], bx[ d ] );
197  evaluate<0>( bx, values );
198  }
199 
200  template< unsigned int deriv, class Vector >
201  void evaluateSingle ( const DomainVector &x, Vector &values ) const
202  {
203  assert(values.size()>=size());
204  coeffMatrix_->template mult<deriv>( eval_.template evaluate<deriv>( x ), values );
205  }
206  template< unsigned int deriv, class Fy >
207  void evaluateSingle ( const DomainVector &x,
208  std::vector< FieldVector<FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size>,dimRange> > &values) const
209  {
210  evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
211  }
212  template< unsigned int deriv, class Fy >
213  void evaluateSingle ( const DomainVector &x,
214  std::vector< FieldVector<LFETensor<Fy,dimension,deriv>,dimRange> > &values) const
215  {
216  evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
217  }
218 
219  template <class Fy>
220  void jacobian ( const DomainVector &x, std::vector<FieldMatrix<Fy,dimRange,dimension> > &values ) const
221  {
222  assert(values.size()>=size());
223  evaluateSingle<1>(x,reinterpret_cast<std::vector<FieldVector<Fy,dimRange*dimension> >&>(values));
224  }
225  template< class DVector, class RVector >
226  void jacobian ( const DVector &x, RVector &values ) const
227  {
228  assert( DVector::dimension == dimension);
229  DomainVector bx;
230  for( unsigned int d = 0; d < dimension; ++d )
231  field_cast( x[ d ], bx[ d ] );
232  jacobian( bx, values );
233  }
234 
235  template <class Fy>
236  void integrate ( std::vector<Fy> &values ) const
237  {
238  assert(values.size()>=size());
239  coeffMatrix_->mult( eval_.template integrate(), values );
240  }
241 
242  protected:
244  : basis_(other.basis_),
245  coeffMatrix_(other.coeffMatrix_),
246  eval_(basis_),
247  order_(basis_.order()),
248  size_(other.size_)
249  {}
251  const Basis &basis_;
252  const CoefficientMatrix* coeffMatrix_;
253  mutable Evaluator eval_;
254  unsigned int order_,size_;
255  };
256 
263  template< class Eval, class CM = SparseCoeffMatrix<typename Eval::Field,Eval::dimRange>,
264  class D=double, class R=double>
266  : public PolynomialBasis< Eval, CM, D, R >
267  {
268  public:
269  typedef CM CoefficientMatrix;
270 
271  private:
272  typedef Eval Evaluator;
273 
276 
277  public:
278  typedef typename Base::Basis Basis;
279 
281  : Base(basis,coeffMatrix_,0)
282  {}
283 
284  template <class Matrix>
285  void fill(const Matrix& matrix)
286  {
287  coeffMatrix_.fill(matrix);
288  this->size_ = coeffMatrix_.size();
289  }
290  template <class Matrix>
291  void fill(const Matrix& matrix,int size)
292  {
293  coeffMatrix_.fill(matrix);
294  assert(size<=coeffMatrix_.size());
295  this->size_ = size;
296  }
297 
298  private:
301  CoefficientMatrix coeffMatrix_;
302  };
303 }
304 #endif // DUNE_POLYNOMIALBASIS_HH
PolynomialBasis & operator=(const PolynomialBasis &)
void integrate(std::vector< Fy > &values) const
Definition: polynomialbasis.hh:236
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: polynomialbasis.hh:114
Definition: tensor.hh:30
static const DomainVector & apply(const DomainVector &x)
Definition: polynomialbasis.hh:172
Definition: polynomialbasis.hh:158
CoefficientMatrix::Field StorageField
Definition: polynomialbasis.hh:70
void fill(const Matrix &matrix, int size)
Definition: polynomialbasis.hh:291
static const unsigned int dimRange
Definition: polynomialbasis.hh:73
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:31
void evaluate(const DVector &x, RVector &values) const
Definition: polynomialbasis.hh:178
D DomainType
domain type
Definition: localbasis.hh:43
void fill(const Matrix &matrix)
Definition: polynomialbasis.hh:285
void evaluate(const DomainVector &x, F *values) const
Definition: polynomialbasis.hh:143
unsigned int order() const
Definition: polynomialbasis.hh:103
Evaluator::DomainVector DomainVector
Definition: polynomialbasis.hh:78
void evaluate(const DVector &x, RVector &values) const
Definition: polynomialbasis.hh:191
void evaluateSingle(const DomainVector &x, std::vector< FieldVector< FieldVector< Fy, LFETensor< Fy, dimension, deriv >::size >, dimRange > > &values) const
Definition: polynomialbasis.hh:207
void jacobian(const DomainVector &x, std::vector< FieldMatrix< Fy, dimRange, dimension > > &values) const
Definition: polynomialbasis.hh:220
PolynomialBasis(const Basis &basis, const CoefficientMatrix &coeffMatrix, unsigned int size)
Definition: polynomialbasis.hh:80
PolynomialBasisWithMatrix(const Basis &basis)
Definition: polynomialbasis.hh:280
LocalBasisTraits< D, dimension, FieldVector< D, dimension >, R, dimRange, FieldVector< R, dimRange >, FieldMatrix< R, dimRange, dimension > > Traits
Definition: polynomialbasis.hh:76
static DomainVector apply(const DVector &x)
Definition: polynomialbasis.hh:160
unsigned int size() const
Definition: polynomialbasis.hh:108
Evaluator::Basis Basis
Definition: polynomialbasis.hh:77
unsigned int order_
Definition: polynomialbasis.hh:254
CM CoefficientMatrix
Definition: polynomialbasis.hh:269
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
const CoefficientMatrix & matrix() const
Definition: polynomialbasis.hh:98
void evaluate(const DVector &x, F *values) const
Definition: polynomialbasis.hh:148
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
const Basis & basis_
Definition: polynomialbasis.hh:251
void evaluate(const DomainVector &x, std::vector< FieldVector< Fy, dimRange > > &values) const
Definition: polynomialbasis.hh:186
void jacobian(const DVector &x, RVector &values) const
Definition: polynomialbasis.hh:226
Base::Basis Basis
Definition: polynomialbasis.hh:278
Definition: polynomialbasis.hh:265
Definition: polynomialbasis.hh:62
Evaluator eval_
Definition: polynomialbasis.hh:253
unsigned int size_
Definition: polynomialbasis.hh:254
void evaluateSingle(const DomainVector &x, Vector &values) const
Definition: polynomialbasis.hh:201
const CoefficientMatrix * coeffMatrix_
Definition: polynomialbasis.hh:252
void partial(const std::array< unsigned int, dimension > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of all shape functions.
Definition: polynomialbasis.hh:130
static const unsigned int dimension
Definition: polynomialbasis.hh:72
PolynomialBasis(const PolynomialBasis &other)
Definition: polynomialbasis.hh:243
const Basis & basis() const
Definition: polynomialbasis.hh:93
CM CoefficientMatrix
Definition: polynomialbasis.hh:68
void evaluateSingle(const DomainVector &x, std::vector< FieldVector< LFETensor< Fy, dimension, deriv >, dimRange > > &values) const
Definition: polynomialbasis.hh:213
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: polynomialbasis.hh:122