dune-localfunctions 2.10
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// SPDX-FileCopyrightInfo: Copyright © 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_POLYNOMIALBASIS_HH
6#define DUNE_POLYNOMIALBASIS_HH
7
8#include <fstream>
9#include <numeric>
10
11#include <dune/common/fmatrix.hh>
12
14
19
20namespace Dune
21{
22
23 // PolynomialBasis
24 // ---------------
25
63 template< class Eval, class CM, class D=double, class R=double >
65 {
67 typedef Eval Evaluator;
68
69 public:
71
72 typedef typename CoefficientMatrix::Field StorageField;
73
74 static const unsigned int dimension = Evaluator::dimension;
75 static const unsigned int dimRange = Evaluator::dimRange*CoefficientMatrix::blockSize;
77 R,dimRange,FieldVector<R,dimRange>,
78 FieldMatrix<R,dimRange,dimension> > Traits;
79 typedef typename Evaluator::Basis Basis;
80 typedef typename Evaluator::DomainVector DomainVector;
81 template <class Fy>
82 using HessianFyType = FieldVector<FieldMatrix<Fy,dimension,dimension>,dimRange>;
84
86 const CoefficientMatrix &coeffMatrix,
87 unsigned int size)
88 : basis_(basis),
89 coeffMatrix_(&coeffMatrix),
90 eval_(basis),
92 size_(size)
93 {
94 // assert(coeffMatrix_);
95 // assert(size_ <= coeffMatrix.size()); // !!!
96 }
97
98 const Basis &basis () const
99 {
100 return basis_;
101 }
102
103 const CoefficientMatrix &matrix () const
104 {
105 return *coeffMatrix_;
106 }
107
108 unsigned int order () const
109 {
110 return order_;
111 }
112
113 unsigned int size () const
114 {
115 return size_;
116 }
117
119 void evaluateFunction (const typename Traits::DomainType& x,
120 std::vector<typename Traits::RangeType>& out) const
121 {
122 out.resize(size());
123 evaluate(x,out);
124 }
125
127 void evaluateJacobian (const typename Traits::DomainType& x, // position
128 std::vector<typename Traits::JacobianType>& out) const // return value
129 {
130 out.resize(size());
131 jacobian(x,out);
132 }
133
135 void evaluateHessian (const typename Traits::DomainType& x, // position
136 std::vector<HessianType>& out) const // return value
137 {
138 out.resize(size());
139 hessian(x,out);
140 }
141
143 void partial (const std::array<unsigned int, dimension>& order,
144 const typename Traits::DomainType& in, // position
145 std::vector<typename Traits::RangeType>& out) const // return value
146 {
147 out.resize(size());
148 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
149 if (totalOrder == 0) {
150 evaluateFunction(in, out);
151 }
152 else if (totalOrder == 1) {
153 std::vector<typename Traits::JacobianType> jacs(out.size());
154 unsigned int k;
155 for (unsigned int i=0;i<order.size();++i)
156 if (order[i]==1) k=i;
157 evaluateJacobian(in, jacs);
158 for (unsigned int i=0;i<out.size();++i)
159 for (unsigned int r=0;r<Traits::RangeType::dimension;++r)
160 out[i][r] = jacs[i][r][k];
161 }
162 else if (totalOrder == 2) {
163 std::vector<HessianType> hesss(out.size());
164 int k=-1,l=-1;
165 for (unsigned int i=0;i<order.size();++i) {
166 if (order[i] >= 1 && k == -1)
167 k = i;
168 else if (order[i]==1) l=i;
169 }
170 if (l==-1) l=k;
171 evaluateHessian(in, hesss);
172 for (unsigned int i=0;i<out.size();++i)
173 for (unsigned int r=0;r<Traits::RangeType::dimension;++r)
174 out[i][r] = hesss[i][r][k][l];
175 }
176 else {
177 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
178 }
179 }
180
181 template< unsigned int deriv, class F >
182 void evaluate ( const DomainVector &x, F *values ) const
183 {
184 coeffMatrix_->mult( eval_.template evaluate<deriv>( x ), size(), values);
185 }
186 template< unsigned int deriv, class DVector, class F >
187 void evaluate ( const DVector &x, F *values ) const
188 {
189 assert( DVector::dimension == dimension);
190 DomainVector bx;
191 for( int d = 0; d < dimension; ++d )
192 field_cast( x[ d ], bx[ d ] );
193 evaluate<deriv>( bx, values );
194 }
195
196 template <bool dummy,class DVector>
197 struct Convert
198 {
199 static DomainVector apply( const DVector &x )
200 {
201 assert( DVector::dimension == dimension);
202 DomainVector bx;
203 for( unsigned int d = 0; d < dimension; ++d )
204 field_cast( x[ d ], bx[ d ] );
205 return bx;
206 }
207 };
208 template <bool dummy>
209 struct Convert<dummy,DomainVector>
210 {
211 static const DomainVector &apply( const DomainVector &x )
212 {
213 return x;
214 }
215 };
216 template< unsigned int deriv, class DVector, class RVector >
217 void evaluate ( const DVector &x, RVector &values ) const
218 {
219 assert(values.size()>=size());
221 coeffMatrix_->mult( eval_.template evaluate<deriv>( bx ), values );
222 }
223
224 template <class Fy>
225 void evaluate ( const DomainVector &x, std::vector<FieldVector<Fy,dimRange> > &values ) const
226 {
227 evaluate<0>(x,values);
228 }
229 template< class DVector, class RVector >
230 void evaluate ( const DVector &x, RVector &values ) const
231 {
232 assert( DVector::dimension == dimension);
233 DomainVector bx;
234 for( unsigned int d = 0; d < dimension; ++d )
235 field_cast( x[ d ], bx[ d ] );
236 evaluate<0>( bx, values );
237 }
238
239 template< unsigned int deriv, class Vector >
240 void evaluateSingle ( const DomainVector &x, Vector &values ) const
241 {
242 assert(values.size()>=size());
243 coeffMatrix_->template mult<deriv>( eval_.template evaluate<deriv>( x ), values );
244 }
245 template< unsigned int deriv, class Fy >
247 std::vector< FieldVector<FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size>,dimRange> > &values) const
248 {
249 evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
250 }
251 template< unsigned int deriv, class Fy >
253 std::vector< FieldVector<LFETensor<Fy,dimension,deriv>,dimRange> > &values) const
254 {
255 evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
256 }
257
258 template <class Fy>
259 void jacobian ( const DomainVector &x,
260 std::vector<FieldMatrix<Fy,dimRange,dimension> > &values ) const
261 {
262 assert(values.size()>=size());
263 evaluateSingle<1>(x,reinterpret_cast<std::vector<FieldVector<Fy,dimRange*dimension> >&>(values));
264 }
265 template< class DVector, class RVector >
266 void jacobian ( const DVector &x, RVector &values ) const
267 {
268 assert( DVector::dimension == dimension);
269 DomainVector bx;
270 for( unsigned int d = 0; d < dimension; ++d )
271 field_cast( x[ d ], bx[ d ] );
272 jacobian( bx, values );
273 }
274 template <class Fy>
275 void hessian ( const DomainVector &x,
276 std::vector<HessianFyType<Fy>> &values ) const
277 {
278 assert(values.size()>=size());
279 // only upper part of hessians matrix is computed - so we have
280 // y[0] = FV< FV<Fy,d*(d+1)/2>, dimRange>
281 const unsigned int hsize = LFETensor<Fy,dimension,2>::size;
282 std::vector< FieldVector< FieldVector<Fy,hsize>, dimRange> > y( size() );
283 evaluateSingle<2>(x, y);
284 unsigned int q = 0;
285 for (unsigned int i = 0; i < size(); ++i)
286 for (unsigned int r = 0; r < dimRange; ++r)
287 {
288 q = 0;
289 // tensor-based things follow unintuitive index scheme
290 // e.g. for dim = 3, the k-l index of y is 00,01,11,02,12,22, i.e. partial derivatives
291 // are ordered: xx,xy,yy,xz,yz,zz
292
293 // Fill values 'directionwise'
294 for (unsigned int k = 0; k < dimension; ++k)
295 for (unsigned int l = 0; l <= k; ++l)
296 {
297
298 values[i][r][k][l] = y[i][r][q];
299 values[i][r][l][k] = y[i][r][q];
300 assert(q < hsize);
301 ++q;
302 }
303 }
304 // evaluateSingle<2>(x,reinterpret_cast<std::vector<FieldVector<Fy,dimRange*dimension*dimension> >&>(values));
305 }
306 template< class DVector, class HVector >
307 void hessian ( const DVector &x, HVector &values ) const
308 {
309 assert( DVector::dimension == dimension);
310 DomainVector bx;
311 for( unsigned int d = 0; d < dimension; ++d )
312 field_cast( x[ d ], bx[ d ] );
313 hessian( bx, values );
314 }
315
316 template <class Fy>
317 void integrate ( std::vector<Fy> &values ) const
318 {
319 assert(values.size()>=size());
320 coeffMatrix_->mult( eval_.template integrate(), values );
321 }
322
323 protected:
325 : basis_(other.basis_),
327 eval_(basis_),
329 size_(other.size_)
330 {}
332 const Basis &basis_;
334 mutable Evaluator eval_;
335 unsigned int order_,size_;
336 };
337
344 template< class Eval, class CM = SparseCoeffMatrix<typename Eval::Field,Eval::dimRange>,
345 class D=double, class R=double>
347 : public PolynomialBasis< Eval, CM, D, R >
348 {
349 public:
351
352 private:
353 typedef Eval Evaluator;
354
357
358 public:
359 typedef typename Base::Basis Basis;
360
362 : Base(basis,coeffMatrix_,0)
363 {}
364
365 template <class Matrix>
366 void fill(const Matrix& matrix)
367 {
368 coeffMatrix_.fill(matrix);
369 this->size_ = coeffMatrix_.size();
370 }
371 template <class Matrix>
372 void fill(const Matrix& matrix,int size)
373 {
374 coeffMatrix_.fill(matrix);
375 assert(size<=coeffMatrix_.size());
376 this->size_ = size;
377 }
378
379 private:
382 CoefficientMatrix coeffMatrix_;
383 };
384}
385#endif // DUNE_POLYNOMIALBASIS_HH
Definition bdfmcube.hh:18
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition field.hh:159
Type traits for LocalBasisVirtualInterface.
Definition common/localbasis.hh:35
D DomainType
domain type
Definition common/localbasis.hh:43
Definition polynomialbasis.hh:65
void evaluate(const DVector &x, RVector &values) const
Definition polynomialbasis.hh:230
void evaluate(const DomainVector &x, std::vector< FieldVector< Fy, dimRange > > &values) const
Definition polynomialbasis.hh:225
PolynomialBasis(const PolynomialBasis &other)
Definition polynomialbasis.hh:324
void evaluate(const DVector &x, F *values) const
Definition polynomialbasis.hh:187
void evaluateHessian(const typename Traits::DomainType &x, std::vector< HessianType > &out) const
Evaluate Jacobian of all shape functions.
Definition polynomialbasis.hh:135
CoefficientMatrix::Field StorageField
Definition polynomialbasis.hh:72
static const unsigned int dimRange
Definition polynomialbasis.hh:75
void jacobian(const DVector &x, RVector &values) const
Definition polynomialbasis.hh:266
Evaluator::DomainVector DomainVector
Definition polynomialbasis.hh:80
Evaluator::Basis Basis
Definition polynomialbasis.hh:79
void evaluateSingle(const DomainVector &x, Vector &values) const
Definition polynomialbasis.hh:240
void evaluateSingle(const DomainVector &x, std::vector< FieldVector< LFETensor< Fy, dimension, deriv >, dimRange > > &values) const
Definition polynomialbasis.hh:252
void jacobian(const DomainVector &x, std::vector< FieldMatrix< Fy, dimRange, dimension > > &values) const
Definition polynomialbasis.hh:259
const CoefficientMatrix & matrix() const
Definition polynomialbasis.hh:103
const Basis & basis_
Definition polynomialbasis.hh:332
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition polynomialbasis.hh:119
static const unsigned int dimension
Definition polynomialbasis.hh:74
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition polynomialbasis.hh:127
PolynomialBasis & operator=(const PolynomialBasis &)
const CoefficientMatrix * coeffMatrix_
Definition polynomialbasis.hh:333
void integrate(std::vector< Fy > &values) const
Definition polynomialbasis.hh:317
unsigned int size() const
Definition polynomialbasis.hh:113
void evaluate(const DomainVector &x, F *values) const
Definition polynomialbasis.hh:182
void hessian(const DVector &x, HVector &values) const
Definition polynomialbasis.hh:307
CM CoefficientMatrix
Definition polynomialbasis.hh:70
HessianFyType< R > HessianType
Definition polynomialbasis.hh:83
LocalBasisTraits< D, dimension, FieldVector< D, dimension >, R, dimRange, FieldVector< R, dimRange >, FieldMatrix< R, dimRange, dimension > > Traits
Definition polynomialbasis.hh:78
unsigned int order_
Definition polynomialbasis.hh:335
void hessian(const DomainVector &x, std::vector< HessianFyType< Fy > > &values) const
Definition polynomialbasis.hh:275
void evaluate(const DVector &x, RVector &values) const
Definition polynomialbasis.hh:217
PolynomialBasis(const Basis &basis, const CoefficientMatrix &coeffMatrix, unsigned int size)
Definition polynomialbasis.hh:85
FieldVector< FieldMatrix< Fy, dimension, dimension >, dimRange > HessianFyType
Definition polynomialbasis.hh:82
unsigned int order() const
Definition polynomialbasis.hh:108
void evaluateSingle(const DomainVector &x, std::vector< FieldVector< FieldVector< Fy, LFETensor< Fy, dimension, deriv >::size >, dimRange > > &values) const
Definition polynomialbasis.hh:246
const Basis & basis() const
Definition polynomialbasis.hh:98
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:143
unsigned int size_
Definition polynomialbasis.hh:335
Evaluator eval_
Definition polynomialbasis.hh:334
Definition polynomialbasis.hh:198
static DomainVector apply(const DVector &x)
Definition polynomialbasis.hh:199
static const DomainVector & apply(const DomainVector &x)
Definition polynomialbasis.hh:211
Definition polynomialbasis.hh:348
PolynomialBasisWithMatrix(const Basis &basis)
Definition polynomialbasis.hh:361
CM CoefficientMatrix
Definition polynomialbasis.hh:350
void fill(const Matrix &matrix, int size)
Definition polynomialbasis.hh:372
Base::Basis Basis
Definition polynomialbasis.hh:359
void fill(const Matrix &matrix)
Definition polynomialbasis.hh:366
Definition tensor.hh:33