DUNE-FEM (unstable)

polynomialbasis.hh
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
12
13#include <dune/localfunctions/common/localbasis.hh>
14
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>
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:
70 typedef CM CoefficientMatrix;
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>,
79 typedef typename Evaluator::Basis Basis;
80 typedef typename Evaluator::DomainVector DomainVector;
81 template <class Fy>
84
85 PolynomialBasis (const Basis &basis,
86 const CoefficientMatrix &coeffMatrix,
87 unsigned int size)
88 : basis_(basis),
89 coeffMatrix_(&coeffMatrix),
90 eval_(basis),
91 order_(basis.order()),
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());
220 const DomainVector &bx = Convert<true,DVector>::apply(x);
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 >
246 void evaluateSingle ( const DomainVector &x,
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 >
252 void evaluateSingle ( const DomainVector &x,
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:
324 PolynomialBasis(const PolynomialBasis &other)
325 : basis_(other.basis_),
326 coeffMatrix_(other.coeffMatrix_),
327 eval_(basis_),
328 order_(basis_.order()),
329 size_(other.size_)
330 {}
331 PolynomialBasis &operator=(const PolynomialBasis&);
332 const Basis &basis_;
333 const CoefficientMatrix* coeffMatrix_;
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:
350 typedef CM CoefficientMatrix;
351
352 private:
353 typedef Eval Evaluator;
354
357
358 public:
359 typedef typename Base::Basis Basis;
360
361 PolynomialBasisWithMatrix (const Basis &basis)
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
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)