dune-localfunctions 3.0-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
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
81 const CoefficientMatrix &coeffMatrix,
82 unsigned int size)
83 : basis_(basis),
84 coeffMatrix_(&coeffMatrix),
85 eval_(basis),
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());
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 >
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 >
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_),
246 eval_(basis_),
248 size_(other.size_)
249 {}
251 const Basis &basis_;
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:
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
Definition brezzidouglasmarini1cube2d.hh:14
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 localbasis.hh:38
D DomainType
domain type
Definition localbasis.hh:49
Definition polynomialbasis.hh:63
void evaluate(const DVector &x, RVector &values) const
Definition polynomialbasis.hh:191
void evaluate(const DomainVector &x, std::vector< FieldVector< Fy, dimRange > > &values) const
Definition polynomialbasis.hh:186
PolynomialBasis(const PolynomialBasis &other)
Definition polynomialbasis.hh:243
void evaluate(const DVector &x, F *values) const
Definition polynomialbasis.hh:148
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:226
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:201
void evaluateSingle(const DomainVector &x, std::vector< FieldVector< LFETensor< Fy, dimension, deriv >, dimRange > > &values) const
Definition polynomialbasis.hh:213
void jacobian(const DomainVector &x, std::vector< FieldMatrix< Fy, dimRange, dimension > > &values) const
Definition polynomialbasis.hh:220
const CoefficientMatrix & matrix() const
Definition polynomialbasis.hh:98
const Basis & basis_
Definition polynomialbasis.hh:251
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition polynomialbasis.hh:114
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:122
PolynomialBasis & operator=(const PolynomialBasis &)
const CoefficientMatrix * coeffMatrix_
Definition polynomialbasis.hh:252
void integrate(std::vector< Fy > &values) const
Definition polynomialbasis.hh:236
unsigned int size() const
Definition polynomialbasis.hh:108
void evaluate(const DomainVector &x, F *values) const
Definition polynomialbasis.hh:143
CM CoefficientMatrix
Definition polynomialbasis.hh:68
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:254
void evaluate(const DVector &x, RVector &values) const
Definition polynomialbasis.hh:178
PolynomialBasis(const Basis &basis, const CoefficientMatrix &coeffMatrix, unsigned int size)
Definition polynomialbasis.hh:80
unsigned int order() const
Definition polynomialbasis.hh:103
void evaluateSingle(const DomainVector &x, std::vector< FieldVector< FieldVector< Fy, LFETensor< Fy, dimension, deriv >::size >, dimRange > > &values) const
Definition polynomialbasis.hh:207
const Basis & basis() const
Definition polynomialbasis.hh:93
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
unsigned int size_
Definition polynomialbasis.hh:254
Evaluator eval_
Definition polynomialbasis.hh:253
Definition polynomialbasis.hh:159
static DomainVector apply(const DVector &x)
Definition polynomialbasis.hh:160
static const DomainVector & apply(const DomainVector &x)
Definition polynomialbasis.hh:172
Definition polynomialbasis.hh:267
PolynomialBasisWithMatrix(const Basis &basis)
Definition polynomialbasis.hh:280
CM CoefficientMatrix
Definition polynomialbasis.hh:269
void fill(const Matrix &matrix, int size)
Definition polynomialbasis.hh:291
Base::Basis Basis
Definition polynomialbasis.hh:278
void fill(const Matrix &matrix)
Definition polynomialbasis.hh:285
Definition tensor.hh:31