5#ifndef DUNE_POLYNOMIALBASIS_HH
6#define DUNE_POLYNOMIALBASIS_HH
13#include <dune/localfunctions/common/localbasis.hh>
15#include <dune/localfunctions/utility/coeffmatrix.hh>
16#include <dune/localfunctions/utility/monomialbasis.hh>
17#include <dune/localfunctions/utility/multiindex.hh>
18#include <dune/localfunctions/utility/basisevaluator.hh>
63 template<
class Eval,
class CM,
class D=
double,
class R=
double >
67 typedef Eval Evaluator;
70 typedef CM CoefficientMatrix;
72 typedef typename CoefficientMatrix::Field StorageField;
74 static const unsigned int dimension = Evaluator::dimension;
75 static const unsigned int dimRange = Evaluator::dimRange*CoefficientMatrix::blockSize;
79 typedef typename Evaluator::Basis Basis;
80 typedef typename Evaluator::DomainVector DomainVector;
86 const CoefficientMatrix &coeffMatrix,
89 coeffMatrix_(&coeffMatrix),
91 order_(basis.order()),
98 const Basis &basis ()
const
103 const CoefficientMatrix &matrix ()
const
105 return *coeffMatrix_;
108 unsigned int order ()
const
113 unsigned int size ()
const
120 std::vector<typename Traits::RangeType>& out)
const
128 std::vector<typename Traits::JacobianType>& out)
const
136 std::vector<HessianType>& out)
const
143 void partial (
const std::array<unsigned int, dimension>& order,
145 std::vector<typename Traits::RangeType>& out)
const
149 if (totalOrder == 0) {
152 else if (totalOrder == 1) {
153 std::vector<typename Traits::JacobianType> jacs(out.size());
155 for (
unsigned int i=0;i<order.size();++i)
156 if (order[i]==1) k=i;
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];
162 else if (totalOrder == 2) {
163 std::vector<HessianType> hesss(out.size());
165 for (
unsigned int i=0;i<order.size();++i) {
166 if (order[i] >= 1 && k == -1)
168 else if (order[i]==1) l=i;
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];
181 template<
unsigned int deriv,
class F >
182 void evaluate (
const DomainVector &x, F *values )
const
184 coeffMatrix_->mult( eval_.template evaluate<deriv>( x ), size(), values);
186 template<
unsigned int deriv,
class DVector,
class F >
187 void evaluate (
const DVector &x, F *values )
const
189 assert( DVector::dimension == dimension);
191 for(
int d = 0; d < dimension; ++d )
193 evaluate<deriv>( bx, values );
196 template <
bool dummy,
class DVector>
199 static DomainVector apply(
const DVector &x )
201 assert( DVector::dimension == dimension);
203 for(
unsigned int d = 0; d < dimension; ++d )
208 template <
bool dummy>
209 struct Convert<dummy,DomainVector>
211 static const DomainVector &apply(
const DomainVector &x )
216 template<
unsigned int deriv,
class DVector,
class RVector >
217 void evaluate (
const DVector &x, RVector &values )
const
219 assert(values.size()>=size());
220 const DomainVector &bx = Convert<true,DVector>::apply(x);
221 coeffMatrix_->mult( eval_.template evaluate<deriv>( bx ), values );
225 void evaluate (
const DomainVector &x, std::vector<FieldVector<Fy,dimRange> > &values )
const
227 evaluate<0>(x,values);
229 template<
class DVector,
class RVector >
230 void evaluate (
const DVector &x, RVector &values )
const
232 assert( DVector::dimension == dimension);
234 for(
unsigned int d = 0; d < dimension; ++d )
236 evaluate<0>( bx, values );
239 template<
unsigned int deriv,
class Vector >
240 void evaluateSingle (
const DomainVector &x, Vector &values )
const
242 assert(values.size()>=size());
243 coeffMatrix_->template mult<deriv>( eval_.template evaluate<deriv>( x ), values );
245 template<
unsigned int deriv,
class Fy >
246 void evaluateSingle (
const DomainVector &x,
247 std::vector< FieldVector<FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size>,dimRange> > &values)
const
249 evaluateSingle<deriv>(x,
reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange
> >&>(values));
251 template<
unsigned int deriv,
class Fy >
252 void evaluateSingle (
const DomainVector &x,
253 std::vector< FieldVector<LFETensor<Fy,dimension,deriv>,dimRange> > &values)
const
255 evaluateSingle<deriv>(x,
reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange
> >&>(values));
259 void jacobian (
const DomainVector &x,
260 std::vector<FieldMatrix<Fy,dimRange,dimension> > &values )
const
262 assert(values.size()>=size());
263 evaluateSingle<1>(x,
reinterpret_cast<std::vector<FieldVector<Fy,dimRange*dimension>
>&>(values));
265 template<
class DVector,
class RVector >
266 void jacobian (
const DVector &x, RVector &values )
const
268 assert( DVector::dimension == dimension);
270 for(
unsigned int d = 0; d < dimension; ++d )
272 jacobian( bx, values );
275 void hessian (
const DomainVector &x,
276 std::vector<HessianFyType<Fy>> &values )
const
278 assert(values.size()>=size());
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);
285 for (
unsigned int i = 0; i < size(); ++i)
286 for (
unsigned int r = 0; r < dimRange; ++r)
294 for (
unsigned int k = 0; k < dimension; ++k)
295 for (
unsigned int l = 0; l <= k; ++l)
298 values[i][r][k][l] = y[i][r][q];
299 values[i][r][l][k] = y[i][r][q];
306 template<
class DVector,
class HVector >
307 void hessian (
const DVector &x, HVector &values )
const
309 assert( DVector::dimension == dimension);
311 for(
unsigned int d = 0; d < dimension; ++d )
313 hessian( bx, values );
317 void integrate ( std::vector<Fy> &values )
const
319 assert(values.size()>=size());
320 coeffMatrix_->mult( eval_.template integrate(), values );
324 PolynomialBasis(
const PolynomialBasis &other)
325 : basis_(other.basis_),
326 coeffMatrix_(other.coeffMatrix_),
328 order_(basis_.order()),
331 PolynomialBasis &operator=(
const PolynomialBasis&);
333 const CoefficientMatrix* coeffMatrix_;
334 mutable Evaluator eval_;
335 unsigned int order_,size_;
344 template<
class Eval,
class CM = SparseCoeffMatrix<
typename Eval::Field,Eval::dimRange>,
345 class D=
double,
class R=
double>
350 typedef CM CoefficientMatrix;
353 typedef Eval Evaluator;
359 typedef typename Base::Basis Basis;
362 :
Base(basis,coeffMatrix_,0)
365 template <
class Matrix>
366 void fill(
const Matrix& matrix)
368 coeffMatrix_.fill(matrix);
369 this->size_ = coeffMatrix_.size();
371 template <
class Matrix>
372 void fill(
const Matrix& matrix,
int size)
374 coeffMatrix_.fill(matrix);
375 assert(size<=coeffMatrix_.size());
382 CoefficientMatrix coeffMatrix_;
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
A generic dynamic dense matrix.
Definition: matrix.hh:561
Default exception for dummy implementations.
Definition: exceptions.hh:263
Definition: polynomialbasis.hh:348
Definition: polynomialbasis.hh:65
void evaluateHessian(const typename Traits::DomainType &x, std::vector< HessianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: polynomialbasis.hh:135
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: polynomialbasis.hh:119
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: polynomialbasis.hh:127
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
Implements a matrix constructed from a given type representing a field and compile-time given number ...
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune namespace.
Definition: alignedallocator.hh:13
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: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43