dune-localfunctions 2.8.0
Loading...
Searching...
No Matches
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
18namespace 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:
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 template <class Fy>
80 using HessianFyType = FieldVector<FieldMatrix<Fy,dimension,dimension>,dimRange>;
82
84 const CoefficientMatrix &coeffMatrix,
85 unsigned int size)
86 : basis_(basis),
87 coeffMatrix_(&coeffMatrix),
88 eval_(basis),
90 size_(size)
91 {
92 // assert(coeffMatrix_);
93 // assert(size_ <= coeffMatrix.size()); // !!!
94 }
95
96 const Basis &basis () const
97 {
98 return basis_;
99 }
100
101 const CoefficientMatrix &matrix () const
102 {
103 return *coeffMatrix_;
104 }
105
106 unsigned int order () const
107 {
108 return order_;
109 }
110
111 unsigned int size () const
112 {
113 return size_;
114 }
115
117 void evaluateFunction (const typename Traits::DomainType& x,
118 std::vector<typename Traits::RangeType>& out) const
119 {
120 out.resize(size());
121 evaluate(x,out);
122 }
123
125 void evaluateJacobian (const typename Traits::DomainType& x, // position
126 std::vector<typename Traits::JacobianType>& out) const // return value
127 {
128 out.resize(size());
129 jacobian(x,out);
130 }
131
133 void evaluateHessian (const typename Traits::DomainType& x, // position
134 std::vector<HessianType>& out) const // return value
135 {
136 out.resize(size());
137 hessian(x,out);
138 }
139
141 void partial (const std::array<unsigned int, dimension>& order,
142 const typename Traits::DomainType& in, // position
143 std::vector<typename Traits::RangeType>& out) const // return value
144 {
145 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
146 if (totalOrder == 0) {
147 evaluateFunction(in, out);
148 }
149 else if (totalOrder == 1) {
150 std::vector<typename Traits::JacobianType> jacs(out.size());
151 unsigned int k;
152 for (unsigned int i=0;i<order.size();++i)
153 if (order[i]==1) k=i;
154 evaluateJacobian(in, jacs);
155 for (unsigned int i=0;i<out.size();++i)
156 for (unsigned int r=0;r<Traits::RangeType::dimension;++r)
157 out[i][r] = jacs[i][r][k];
158 }
159 else if (totalOrder == 2) {
160 std::vector<HessianType> hesss(out.size());
161 int k=-1,l=-1;
162 for (unsigned int i=0;i<order.size();++i) {
163 if (order[i]==1 && k==-1) k=i;
164 else if (order[i]==1) l=i;
165 }
166 if (l==-1) l=k;
167 evaluateHessian(in, hesss);
168 for (unsigned int i=0;i<out.size();++i)
169 for (unsigned int r=0;r<Traits::RangeType::dimension;++r)
170 out[i][r] = hesss[i][r][k][l];
171 }
172 else {
173 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
174 }
175 }
176
177 template< unsigned int deriv, class F >
178 void evaluate ( const DomainVector &x, F *values ) const
179 {
180 coeffMatrix_->mult( eval_.template evaluate<deriv>( x ), size(), values);
181 }
182 template< unsigned int deriv, class DVector, class F >
183 void evaluate ( const DVector &x, F *values ) const
184 {
185 assert( DVector::dimension == dimension);
186 DomainVector bx;
187 for( int d = 0; d < dimension; ++d )
188 field_cast( x[ d ], bx[ d ] );
189 evaluate<deriv>( bx, values );
190 }
191
192 template <bool dummy,class DVector>
193 struct Convert
194 {
195 static DomainVector apply( const DVector &x )
196 {
197 assert( DVector::dimension == dimension);
198 DomainVector bx;
199 for( unsigned int d = 0; d < dimension; ++d )
200 field_cast( x[ d ], bx[ d ] );
201 return bx;
202 }
203 };
204 template <bool dummy>
205 struct Convert<dummy,DomainVector>
206 {
207 static const DomainVector &apply( const DomainVector &x )
208 {
209 return x;
210 }
211 };
212 template< unsigned int deriv, class DVector, class RVector >
213 void evaluate ( const DVector &x, RVector &values ) const
214 {
215 assert(values.size()>=size());
217 coeffMatrix_->mult( eval_.template evaluate<deriv>( bx ), values );
218 }
219
220 template <class Fy>
221 void evaluate ( const DomainVector &x, std::vector<FieldVector<Fy,dimRange> > &values ) const
222 {
223 evaluate<0>(x,values);
224 }
225 template< class DVector, class RVector >
226 void evaluate ( 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 evaluate<0>( bx, values );
233 }
234
235 template< unsigned int deriv, class Vector >
236 void evaluateSingle ( const DomainVector &x, Vector &values ) const
237 {
238 assert(values.size()>=size());
239 coeffMatrix_->template mult<deriv>( eval_.template evaluate<deriv>( x ), values );
240 }
241 template< unsigned int deriv, class Fy >
243 std::vector< FieldVector<FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size>,dimRange> > &values) const
244 {
245 evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
246 }
247 template< unsigned int deriv, class Fy >
249 std::vector< FieldVector<LFETensor<Fy,dimension,deriv>,dimRange> > &values) const
250 {
251 evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
252 }
253
254 template <class Fy>
255 void jacobian ( const DomainVector &x,
256 std::vector<FieldMatrix<Fy,dimRange,dimension> > &values ) const
257 {
258 assert(values.size()>=size());
259 evaluateSingle<1>(x,reinterpret_cast<std::vector<FieldVector<Fy,dimRange*dimension> >&>(values));
260 }
261 template< class DVector, class RVector >
262 void jacobian ( const DVector &x, RVector &values ) const
263 {
264 assert( DVector::dimension == dimension);
265 DomainVector bx;
266 for( unsigned int d = 0; d < dimension; ++d )
267 field_cast( x[ d ], bx[ d ] );
268 jacobian( bx, values );
269 }
270 template <class Fy>
271 void hessian ( const DomainVector &x,
272 std::vector<HessianFyType<Fy>> &values ) const
273 {
274 assert(values.size()>=size());
275 // only upper part of hessians matrix is computed - so we have
276 // y[0] = FV< FV<Fy,d*(d+1)/2>, dimRange>
277 const unsigned int hsize = LFETensor<Fy,dimension,2>::size;
278 std::vector< FieldVector< FieldVector<Fy,hsize>, dimRange> > y( size() );
279 evaluateSingle<2>( x, y );
280 unsigned int q=0;
281 for (unsigned int i=0;i<size();++i)
282 for (unsigned int r=0;r<dimRange;++r)
283 for (unsigned int k=0;k<dimension;++k)
284 for (unsigned int l=k;l<dimension;++l)
285 {
286 values[i][r][k][l] = y[i][r][q];
287 values[i][r][l][k] = y[i][r][q];
288 ++q;
289 assert(q < hsize);
290 }
291 // evaluateSingle<2>(x,reinterpret_cast<std::vector<FieldVector<Fy,dimRange*dimension*dimension> >&>(values));
292 }
293 template< class DVector, class HVector >
294 void hessian ( const DVector &x, HVector &values ) const
295 {
296 assert( DVector::dimension == dimension);
297 DomainVector bx;
298 for( unsigned int d = 0; d < dimension; ++d )
299 field_cast( x[ d ], bx[ d ] );
300 hessian( bx, values );
301 }
302
303 template <class Fy>
304 void integrate ( std::vector<Fy> &values ) const
305 {
306 assert(values.size()>=size());
307 coeffMatrix_->mult( eval_.template integrate(), values );
308 }
309
310 protected:
312 : basis_(other.basis_),
314 eval_(basis_),
316 size_(other.size_)
317 {}
319 const Basis &basis_;
321 mutable Evaluator eval_;
322 unsigned int order_,size_;
323 };
324
331 template< class Eval, class CM = SparseCoeffMatrix<typename Eval::Field,Eval::dimRange>,
332 class D=double, class R=double>
334 : public PolynomialBasis< Eval, CM, D, R >
335 {
336 public:
338
339 private:
340 typedef Eval Evaluator;
341
344
345 public:
346 typedef typename Base::Basis Basis;
347
349 : Base(basis,coeffMatrix_,0)
350 {}
351
352 template <class Matrix>
353 void fill(const Matrix& matrix)
354 {
355 coeffMatrix_.fill(matrix);
356 this->size_ = coeffMatrix_.size();
357 }
358 template <class Matrix>
359 void fill(const Matrix& matrix,int size)
360 {
361 coeffMatrix_.fill(matrix);
362 assert(size<=coeffMatrix_.size());
363 this->size_ = size;
364 }
365
366 private:
369 CoefficientMatrix coeffMatrix_;
370 };
371}
372#endif // DUNE_POLYNOMIALBASIS_HH
Definition: bdfmcube.hh:16
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
Type traits for LocalBasisVirtualInterface.
Definition: common/localbasis.hh:32
D DomainType
domain type
Definition: common/localbasis.hh:43
Definition: polynomialbasis.hh:63
void evaluate(const DVector &x, RVector &values) const
Definition: polynomialbasis.hh:226
void evaluate(const DomainVector &x, std::vector< FieldVector< Fy, dimRange > > &values) const
Definition: polynomialbasis.hh:221
PolynomialBasis(const PolynomialBasis &other)
Definition: polynomialbasis.hh:311
void evaluate(const DVector &x, F *values) const
Definition: polynomialbasis.hh:183
void evaluateHessian(const typename Traits::DomainType &x, std::vector< HessianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: polynomialbasis.hh:133
CoefficientMatrix::Field StorageField
Definition: polynomialbasis.hh:70
static const unsigned int dimRange
Definition: polynomialbasis.hh:73
void jacobian(const DVector &x, RVector &values) const
Definition: polynomialbasis.hh:262
Evaluator::DomainVector DomainVector
Definition: polynomialbasis.hh:78
Evaluator::Basis Basis
Definition: polynomialbasis.hh:77
void evaluateSingle(const DomainVector &x, Vector &values) const
Definition: polynomialbasis.hh:236
void evaluateSingle(const DomainVector &x, std::vector< FieldVector< LFETensor< Fy, dimension, deriv >, dimRange > > &values) const
Definition: polynomialbasis.hh:248
void jacobian(const DomainVector &x, std::vector< FieldMatrix< Fy, dimRange, dimension > > &values) const
Definition: polynomialbasis.hh:255
const CoefficientMatrix & matrix() const
Definition: polynomialbasis.hh:101
const Basis & basis_
Definition: polynomialbasis.hh:319
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: polynomialbasis.hh:117
static const unsigned int dimension
Definition: polynomialbasis.hh:72
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: polynomialbasis.hh:125
PolynomialBasis & operator=(const PolynomialBasis &)
const CoefficientMatrix * coeffMatrix_
Definition: polynomialbasis.hh:320
void integrate(std::vector< Fy > &values) const
Definition: polynomialbasis.hh:304
unsigned int size() const
Definition: polynomialbasis.hh:111
void evaluate(const DomainVector &x, F *values) const
Definition: polynomialbasis.hh:178
void hessian(const DVector &x, HVector &values) const
Definition: polynomialbasis.hh:294
CM CoefficientMatrix
Definition: polynomialbasis.hh:68
HessianFyType< R > HessianType
Definition: polynomialbasis.hh:81
LocalBasisTraits< D, dimension, FieldVector< D, dimension >, R, dimRange, FieldVector< R, dimRange >, FieldMatrix< R, dimRange, dimension > > Traits
Definition: polynomialbasis.hh:76
unsigned int order_
Definition: polynomialbasis.hh:322
void hessian(const DomainVector &x, std::vector< HessianFyType< Fy > > &values) const
Definition: polynomialbasis.hh:271
void evaluate(const DVector &x, RVector &values) const
Definition: polynomialbasis.hh:213
PolynomialBasis(const Basis &basis, const CoefficientMatrix &coeffMatrix, unsigned int size)
Definition: polynomialbasis.hh:83
FieldVector< FieldMatrix< Fy, dimension, dimension >, dimRange > HessianFyType
Definition: polynomialbasis.hh:80
unsigned int order() const
Definition: polynomialbasis.hh:106
void evaluateSingle(const DomainVector &x, std::vector< FieldVector< FieldVector< Fy, LFETensor< Fy, dimension, deriv >::size >, dimRange > > &values) const
Definition: polynomialbasis.hh:242
const Basis & basis() const
Definition: polynomialbasis.hh:96
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:141
unsigned int size_
Definition: polynomialbasis.hh:322
Evaluator eval_
Definition: polynomialbasis.hh:321
Definition: polynomialbasis.hh:194
static DomainVector apply(const DVector &x)
Definition: polynomialbasis.hh:195
static const DomainVector & apply(const DomainVector &x)
Definition: polynomialbasis.hh:207
Definition: polynomialbasis.hh:335
PolynomialBasisWithMatrix(const Basis &basis)
Definition: polynomialbasis.hh:348
CM CoefficientMatrix
Definition: polynomialbasis.hh:337
void fill(const Matrix &matrix, int size)
Definition: polynomialbasis.hh:359
Base::Basis Basis
Definition: polynomialbasis.hh:346
void fill(const Matrix &matrix)
Definition: polynomialbasis.hh:353
Definition: tensor.hh:31