dune-fem 2.8.0
Loading...
Searching...
No Matches
space/basisfunctionset/default.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_BASISFUNCTIONSET_DEFAULT_HH
2#define DUNE_FEM_BASISFUNCTIONSET_DEFAULT_HH
3
4// C++ includes
5#include <cassert>
6#include <cstddef>
7#include <memory>
8#include <type_traits>
9#include <utility>
10
11// dune-geometry includes
12#include <dune/geometry/referenceelements.hh>
13#include <dune/geometry/type.hh>
14
16
17// dune-fem includes
23#include <dune/fem/version.hh>
24
27
28namespace Dune
29{
30
31 namespace Fem
32 {
33
34 // DefaultBasisFunctionSet
35 // -----------------------
36
50 template< class Entity, class ShapeFunctionSet >
52 {
54
55 public:
57 typedef Entity EntityType;
60
61 // if underlying shape function set was created with storage CodegenStorage
62 // then this value should be true (see selectcaching.hh)
63 static constexpr bool codegenShapeFunctionSet = detail::IsCodegenShapeFunctionSet< ShapeFunctionSetType >::value;
64
65 protected:
69
71
72 typedef typename EntityType::Geometry Geometry ;
73
74 typedef typename Geometry::ctype ctype;
75 public:
76 // slight misuse of struct ToLocalFunctionSpace!!!
79
81 typedef typename FunctionSpaceType::RangeType RangeType;
83 typedef typename FunctionSpaceType::DomainType DomainType;
85 typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
87 typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
88
89 typedef typename FunctionSpaceType::ScalarFunctionSpaceType ScalarFunctionSpaceType;
90
91 typedef typename ScalarFunctionSpaceType::RangeType ScalarRangeType;
92 typedef typename ScalarFunctionSpaceType::JacobianRangeType ScalarJacobianRangeType;
93
95 typedef std::decay_t< decltype( Dune::ReferenceElements< ctype, Geometry::coorddimension >::general( std::declval< const Dune::GeometryType & >() ) ) > ReferenceElementType;
96
97 static const int pointSetId = detail::SelectPointSetId< ShapeFunctionSetType >::value;
98
101
104 : entity_( &entity ),
106 {
107 // Note that this should be geometry_ = entity.geometry()
108 // But Dune::Geometries are not assignable ...
109 geometry_.reset();
110 geometry_.emplace( entity.geometry() );
111 }
112
114 : entity_( other.entity_ ),
116 {
117 // Note that this should be geometry_ = entity.geometry()
118 // But Dune::Geometries are not assignable ...
119 geometry_.reset();
120 if( other.geometry_ )
121 geometry_.emplace( other.geometry_.value() );
122 }
123
125 {
126 entity_ = other.entity_;
128
129 // Note that this should be geometry_ = entity.geometry()
130 // But Dune::Geometries are not assignable ...
131 geometry_.reset();
132 if( other.geometry_ )
133 geometry_.emplace( other.geometry_.value() );
134 return *this;
135 }
136
137 // Basis Function Set Interface Methods
138 // ------------------------------------
139
141 int order () const { return shapeFunctionSet().order(); }
142
144 std::size_t size () const { return shapeFunctionSet().size(); }
145
147 auto referenceElement () const
148 -> decltype( Dune::ReferenceElements< ctype, Geometry::coorddimension >::general( std::declval< const Dune::GeometryType & >() ) )
149 {
150 return Dune::ReferenceElements< ctype, Geometry::coorddimension >::general( type() );
151 }
152
156 template< class QuadratureType, class Vector, class DofVector >
157 void axpy ( const QuadratureType &quad, const Vector &values, DofVector &dofs ) const
158 {
159 axpyImpl( quad, values, dofs, values[ 0 ] );
160 }
161
168 template< class QuadratureType, class VectorA, class VectorB, class DofVector >
169 void axpy ( const QuadratureType &quad, const VectorA &valuesA, const VectorB &valuesB, DofVector &dofs ) const
170 {
171 assert( valuesA.size() > 0 );
172 assert( valuesB.size() > 0 );
173
174 axpyImpl( quad, valuesA, dofs, valuesA[ 0 ] );
175 axpyImpl( quad, valuesB, dofs, valuesB[ 0 ] );
176 }
177
181 template< class Point, class DofVector >
182 void axpy ( const Point &x, const RangeType &valueFactor, DofVector &dofs ) const
183 {
186 }
187
191 template< class Point, class DofVector >
192 void axpy ( const Point &x, const JacobianRangeType &jacobianFactor, DofVector &dofs ) const
193 {
194 typedef typename Geometry::JacobianInverseTransposed GeometryJacobianInverseTransposedType;
195 const Geometry &geo = geometry();
196 const GeometryJacobianInverseTransposedType &gjit = geo.jacobianInverseTransposed( coordinate( x ) );
197 LocalJacobianRangeType tmpJacobianFactor;
198 for( int r = 0; r < FunctionSpaceType::dimRange; ++r )
199 gjit.mtv( jacobianFactor[ r ], tmpJacobianFactor[ r ] );
200
203 }
204
208 template< class Point, class DofVector >
209 void axpy ( const Point &x, const RangeType &valueFactor, const JacobianRangeType &jacobianFactor,
210 DofVector &dofs ) const
211 {
212 axpy( x, valueFactor, dofs );
213 axpy( x, jacobianFactor, dofs );
214 }
215
218 template< class Point, class DofVector >
219 void axpy ( const Point &x, const HessianRangeType &hessianFactor, DofVector &dofs ) const
220 {
221 typedef typename Geometry::JacobianInverseTransposed GeometryJacobianInverseTransposedType;
222 const Geometry &geo = geometry();
223 const GeometryJacobianInverseTransposedType &gjit = geo.jacobianInverseTransposed( coordinate( x ) );
224 LocalHessianRangeType tmpHessianFactor( RangeFieldType(0) );
225 // don't know how to work directly with the DiagonalMatrix
226 // returned from YaspGrid since gjit[k][l] is not correctly
227 // implemented for k!=l so convert first into a dense matrix...
228 Dune::FieldMatrix<double,DomainType::dimension,DomainType::dimension>
229 G = gjit;
230 DomainType Hg;
231 for( int r = 0; r < FunctionSpaceType::dimRange; ++r )
232 for( int j = 0; j < gjit.cols; ++j )
233 for( int k = 0; k < gjit.cols; ++k )
234 {
235 hessianFactor[r].mv(G[k],Hg);
236 tmpHessianFactor[r][j][k] += Hg * G[j];
237 }
240 }
241
242
244 template< class QuadratureType, class DofVector, class RangeArray >
245 void evaluateAll ( const QuadratureType &quad, const DofVector &dofs, RangeArray &ranges ) const
246 {
247 assert( ranges.size() >= quad.nop() );
248
249 // if shape function set supports codegen and quadrature supports caching
250 if constexpr ( codegenShapeFunctionSet && std::is_base_of< CachingInterface, QuadratureType > :: value)
251 {
252 typedef Codegen :: EvaluateCallerInterfaceTraits< QuadratureType, RangeArray, DofVector > Traits;
253 typedef Codegen :: EvaluateCallerInterface< Traits > BaseEvaluationType;
254
255 // get base function evaluate caller
256 const auto& baseEval = BaseEvaluationType::storage( *this, rangeCache( quad ), quad );
257
258 // true if implementation exists, otherwise this is a nullptr
259 if( baseEval )
260 {
261 baseEval->evaluateRanges( quad, dofs, ranges );
262 return ;
263 }
264 }
265
266 {
267 // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
268 const unsigned int nop = quad.nop();
269 for( unsigned int qp = 0; qp < nop; ++qp )
270 {
271 evaluateAll( quad[ qp ], dofs, ranges[ qp ] );
272 }
273 }
274 }
275
277 template< class Point, class DofVector >
278 void evaluateAll ( const Point &x, const DofVector &dofs, RangeType &value ) const
279 {
280 value = RangeType( 0 );
283 }
284
286 template< class Point, class RangeArray >
287 void evaluateAll ( const Point &x, RangeArray &values ) const
288 {
289 assert( values.size() >= size() );
290 std::fill( values.begin(), values.end(), RangeType(0) );
291 AssignFunctor< RangeArray > f( values );
293 }
294
296 template< class QuadratureType, class DofVector, class JacobianArray >
297 void jacobianAll ( const QuadratureType &quad, const DofVector &dofs, JacobianArray &jacobians ) const
298 {
299 assert( jacobians.size() >= quad.nop() );
300
301 // if shape function set supports codegen and quadrature supports caching
302 if constexpr ( codegenShapeFunctionSet && std::is_base_of< CachingInterface, QuadratureType > :: value)
303 {
304 typedef Codegen :: EvaluateCallerInterfaceTraits< QuadratureType, JacobianArray, DofVector, Geometry > Traits;
305 typedef Codegen :: EvaluateCallerInterface< Traits > BaseEvaluationType;
306
307 // get base function evaluate caller (calls axpyRanges)
308 const auto& baseEval = BaseEvaluationType::storage( *this, jacobianCache( quad ), quad );
309
310 // true if implementation exists
311 if( baseEval )
312 {
313 // call appropriate axpyJacobian method
314 const Geometry &geo = geometry();
315 baseEval->evaluateJacobians( quad, geo, dofs, jacobians );
316 return ;
317 }
318 }
319
320 {
321 // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
322 const unsigned int nop = quad.nop();
323 for( unsigned int qp = 0; qp < nop; ++qp )
324 {
325 jacobianAll( quad[ qp ], dofs, jacobians[ qp ] );
326 }
327 }
328 }
329
331 template< class Point, class DofVector >
332 void jacobianAll ( const Point &x, const DofVector &dofs, JacobianRangeType &jacobian ) const
333 {
334 LocalJacobianRangeType localJacobian( RangeFieldType( 0 ) );
337
338 typedef JacobianTransformation< Geometry > Transformation;
339 const Geometry &geo = geometry();
340 Transformation transformation( geo, coordinate( x ) );
341 transformation( localJacobian, jacobian );
342 }
343
345 template< class Point, class JacobianRangeArray >
346 void jacobianAll ( const Point &x, JacobianRangeArray &jacobians ) const
347 {
348 assert( jacobians.size() >= size() );
349 typedef JacobianTransformation< Geometry > Transformation;
350 const Geometry &geo = geometry();
351
352 Transformation transformation( geo, coordinate( x ) );
353 AssignFunctor< JacobianRangeArray, Transformation > f( jacobians, transformation );
355 }
356
358 template< class QuadratureType, class DofVector, class HessianArray >
359 void hessianAll ( const QuadratureType &quad, const DofVector &dofs, HessianArray &hessians ) const
360 {
361 assert( hessians.size() >= quad.nop() );
362 // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
363 const unsigned int nop = quad.nop();
364 for( unsigned int qp = 0; qp < nop; ++qp )
365 {
366 hessianAll( quad[ qp ], dofs, hessians[ qp ] );
367 }
368 }
369
371 template< class Point, class DofVector >
372 void hessianAll ( const Point &x, const DofVector &dofs, HessianRangeType &hessian ) const
373 {
374 LocalHessianRangeType localHessian( typename LocalHessianRangeType::value_type( RangeFieldType( 0 ) ) );
377
378 typedef HessianTransformation< Geometry > Transformation;
379 const Geometry &geo = geometry();
380 Transformation transformation( geo, coordinate( x ) );
381 transformation( localHessian, hessian );
382 }
383
385 template< class Point, class HessianRangeArray >
386 void hessianAll ( const Point &x, HessianRangeArray &hessians ) const
387 {
388 assert( hessians.size() >= size() );
389 typedef HessianTransformation< Geometry > Transformation;
390 const Geometry &geo = geometry();
391 Transformation transformation( geo, coordinate( x ) );
392 AssignFunctor< HessianRangeArray, Transformation > f( hessians, transformation );
394 }
395
397 const Entity &entity () const
398 {
399 assert( valid() );
400 return *entity_;
401 }
402
404 bool valid () const { return bool(entity_); }
405
407 Dune::GeometryType type () const { return entity().type(); }
408
409 // Non-interface methods
410 // ---------------------
411
414
415 protected:
417 template< class QuadratureType, class RangeArray, class DofVector >
418 void axpyImpl ( const QuadratureType &quad, const RangeArray &rangeFactors, DofVector &dofs, const RangeType& ) const
419 {
420 assert( rangeFactors.size() >= quad.nop() );
421
422 // if shape function set supports codegen and quadrature supports caching
423 if constexpr ( codegenShapeFunctionSet && std::is_base_of< CachingInterface, QuadratureType > :: value)
424 {
425 typedef Codegen :: EvaluateCallerInterfaceTraits< QuadratureType, RangeArray, DofVector > Traits;
426 typedef Codegen :: EvaluateCallerInterface< Traits > BaseEvaluationType;
427
428 // get base function evaluate caller (calls axpyRanges)
429 const auto& baseEval = BaseEvaluationType::storage( *this, rangeCache( quad ), quad );
430
431 // true if implementation exists
432 if( baseEval )
433 {
434 baseEval->axpyRanges( quad, rangeFactors, dofs );
435 return ;
436 }
437 }
438
439 {
440 // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
441 const unsigned int nop = quad.nop();
442 for( unsigned int qp = 0; qp < nop; ++qp )
443 {
444 axpy( quad[ qp ], rangeFactors[ qp ], dofs );
445 }
446 }
447 }
448
450 template< class QuadratureType, class JacobianArray, class DofVector >
451 void axpyImpl ( const QuadratureType &quad, const JacobianArray &jacobianFactors, DofVector &dofs, const JacobianRangeType& ) const
452 {
453 assert( jacobianFactors.size() >= quad.nop() );
454 // if shape function set supports codegen and quadrature supports caching
455 if constexpr ( codegenShapeFunctionSet && std::is_base_of< CachingInterface, QuadratureType > :: value)
456 {
457 typedef Codegen :: EvaluateCallerInterfaceTraits< QuadratureType, JacobianArray, DofVector, Geometry > Traits;
458 typedef Codegen :: EvaluateCallerInterface< Traits > BaseEvaluationType;
459
460 // get base function evaluate caller (calls axpyRanges)
461 const auto& baseEval = BaseEvaluationType::storage( *this, jacobianCache( quad ), quad );
462
463 // true if implementation exists
464 if( baseEval )
465 {
466 // call appropriate axpyRanges method
467 const Geometry &geo = geometry();
468 baseEval->axpyJacobians( quad, geo, jacobianFactors, dofs );
469 return ;
470 }
471 }
472
473 {
474 // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
475 const unsigned int nop = quad.nop();
476 for( unsigned int qp = 0; qp < nop; ++qp )
477 {
478 axpy( quad[ qp ], jacobianFactors[ qp ], dofs );
479 }
480 }
481 }
482
483 template <class QuadratureType>
484 const auto& rangeCache( const QuadratureType& quad ) const
485 {
486 return shapeFunctionSet().scalarShapeFunctionSet().impl().rangeCache( quad );
487 }
488
489 template <class QuadratureType>
490 const auto& jacobianCache( const QuadratureType& quad ) const
491 {
492 return shapeFunctionSet().scalarShapeFunctionSet().impl().jacobianCache( quad );
493 }
494
495 protected:
496 Geometry geometry () const { return geometry_.value(); }
497
498 const EntityType *entity_ = nullptr;
500 std::optional< Geometry > geometry_;
501 };
502
503 } // namespace Fem
504
505} // namespace Dune
506
507#endif // #ifndef DUNE_FEM_BASISFUNCTIONSET_DEFAULT_HH
STL namespace.
Definition: bindguard.hh:11
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
static const Point & coordinate(const Point &x)
Definition: coordinate.hh:14
Definition: explicitfieldvector.hh:75
Definition: misc/functor.hh:31
Definition: space/basisfunctionset/default.hh:52
EntityType::Geometry Geometry
Definition: space/basisfunctionset/default.hh:72
FunctionSpaceType::DomainType DomainType
domain type
Definition: space/basisfunctionset/default.hh:83
void evaluateAll(const Point &x, const DofVector &dofs, RangeType &value) const
Definition: space/basisfunctionset/default.hh:278
LocalFunctionSpaceType::HessianRangeType LocalHessianRangeType
Definition: space/basisfunctionset/default.hh:68
DefaultBasisFunctionSet()
constructor
Definition: space/basisfunctionset/default.hh:100
FunctionSpaceType::HessianRangeType HessianRangeType
hessian range type
Definition: space/basisfunctionset/default.hh:87
const EntityType * entity_
Definition: space/basisfunctionset/default.hh:498
ToNewDimDomainFunctionSpace< LocalFunctionSpaceType, Geometry::coorddimension >::Type FunctionSpaceType
type of function space
Definition: space/basisfunctionset/default.hh:78
void axpy(const QuadratureType &quad, const VectorA &valuesA, const VectorB &valuesB, DofVector &dofs) const
evaluate all basis function and multiply with given values and add to dofs
Definition: space/basisfunctionset/default.hh:169
const ShapeFunctionSetType & shapeFunctionSet() const
return shape function set
Definition: space/basisfunctionset/default.hh:413
auto referenceElement() const -> decltype(Dune::ReferenceElements< ctype, Geometry::coorddimension >::general(std::declval< const Dune::GeometryType & >()))
return reference element
Definition: space/basisfunctionset/default.hh:147
void axpy(const Point &x, const RangeType &valueFactor, const JacobianRangeType &jacobianFactor, DofVector &dofs) const
evaluate all basis function and derivatives and multiply with given values and add to dofs
Definition: space/basisfunctionset/default.hh:209
static constexpr bool codegenShapeFunctionSet
Definition: space/basisfunctionset/default.hh:63
std::optional< Geometry > geometry_
Definition: space/basisfunctionset/default.hh:500
FunctionSpaceType::RangeType RangeType
range type
Definition: space/basisfunctionset/default.hh:81
ShapeFunctionSetType::FunctionSpaceType LocalFunctionSpaceType
Definition: space/basisfunctionset/default.hh:66
LocalFunctionSpaceType::RangeFieldType RangeFieldType
Definition: space/basisfunctionset/default.hh:70
LocalFunctionSpaceType::JacobianRangeType LocalJacobianRangeType
Definition: space/basisfunctionset/default.hh:67
void axpy(const Point &x, const JacobianRangeType &jacobianFactor, DofVector &dofs) const
evaluate all derivatives of all basis function and multiply with given values and add to dofs
Definition: space/basisfunctionset/default.hh:192
void jacobianAll(const Point &x, JacobianRangeArray &jacobians) const
Definition: space/basisfunctionset/default.hh:346
bool valid() const
return true if entity pointer is set
Definition: space/basisfunctionset/default.hh:404
DefaultBasisFunctionSet(const EntityType &entity, const ShapeFunctionSet &shapeFunctionSet=ShapeFunctionSet())
constructor
Definition: space/basisfunctionset/default.hh:103
const auto & jacobianCache(const QuadratureType &quad) const
Definition: space/basisfunctionset/default.hh:490
void axpy(const Point &x, const HessianRangeType &hessianFactor, DofVector &dofs) const
Add H:D^2phi to each dof.
Definition: space/basisfunctionset/default.hh:219
FunctionSpaceType::JacobianRangeType JacobianRangeType
jacobian range type
Definition: space/basisfunctionset/default.hh:85
void hessianAll(const Point &x, HessianRangeArray &hessians) const
Definition: space/basisfunctionset/default.hh:386
void evaluateAll(const QuadratureType &quad, const DofVector &dofs, RangeArray &ranges) const
Definition: space/basisfunctionset/default.hh:245
ShapeFunctionSetType shapeFunctionSet_
Definition: space/basisfunctionset/default.hh:499
ScalarFunctionSpaceType::RangeType ScalarRangeType
Definition: space/basisfunctionset/default.hh:91
std::decay_t< decltype(Dune::ReferenceElements< ctype, Geometry::coorddimension >::general(std::declval< const Dune::GeometryType & >())) > ReferenceElementType
type of reference element
Definition: space/basisfunctionset/default.hh:95
DefaultBasisFunctionSet & operator=(const DefaultBasisFunctionSet &other)
Definition: space/basisfunctionset/default.hh:124
void hessianAll(const Point &x, const DofVector &dofs, HessianRangeType &hessian) const
Definition: space/basisfunctionset/default.hh:372
void hessianAll(const QuadratureType &quad, const DofVector &dofs, HessianArray &hessians) const
Definition: space/basisfunctionset/default.hh:359
DefaultBasisFunctionSet(const DefaultBasisFunctionSet &other)
Definition: space/basisfunctionset/default.hh:113
void evaluateAll(const Point &x, RangeArray &values) const
Definition: space/basisfunctionset/default.hh:287
ShapeFunctionSet ShapeFunctionSetType
shape function set type
Definition: space/basisfunctionset/default.hh:59
ScalarFunctionSpaceType::JacobianRangeType ScalarJacobianRangeType
Definition: space/basisfunctionset/default.hh:92
const Entity & entity() const
return entity
Definition: space/basisfunctionset/default.hh:397
static const int pointSetId
Definition: space/basisfunctionset/default.hh:97
std::size_t size() const
return size of basis function set
Definition: space/basisfunctionset/default.hh:144
void jacobianAll(const Point &x, const DofVector &dofs, JacobianRangeType &jacobian) const
Definition: space/basisfunctionset/default.hh:332
Geometry geometry() const
Definition: space/basisfunctionset/default.hh:496
const auto & rangeCache(const QuadratureType &quad) const
Definition: space/basisfunctionset/default.hh:484
void axpyImpl(const QuadratureType &quad, const RangeArray &rangeFactors, DofVector &dofs, const RangeType &) const
evaluate all basis function and multiply with given values and add to dofs
Definition: space/basisfunctionset/default.hh:418
Entity EntityType
entity type
Definition: space/basisfunctionset/default.hh:57
void axpyImpl(const QuadratureType &quad, const JacobianArray &jacobianFactors, DofVector &dofs, const JacobianRangeType &) const
evaluate all basis function and multiply with given values and add to dofs
Definition: space/basisfunctionset/default.hh:451
Dune::GeometryType type() const
return geometry type
Definition: space/basisfunctionset/default.hh:407
void axpy(const Point &x, const RangeType &valueFactor, DofVector &dofs) const
evaluate all basis function and multiply with given values and add to dofs
Definition: space/basisfunctionset/default.hh:182
Geometry::ctype ctype
Definition: space/basisfunctionset/default.hh:74
FunctionSpaceType::ScalarFunctionSpaceType ScalarFunctionSpaceType
Definition: space/basisfunctionset/default.hh:89
void jacobianAll(const QuadratureType &quad, const DofVector &dofs, JacobianArray &jacobians) const
Definition: space/basisfunctionset/default.hh:297
void axpy(const QuadratureType &quad, const Vector &values, DofVector &dofs) const
evaluate all basis function and multiply with given values and add to dofs
Definition: space/basisfunctionset/default.hh:157
int order() const
return order of basis function set
Definition: space/basisfunctionset/default.hh:141
Definition: space/basisfunctionset/functor.hh:108
Definition: space/basisfunctionset/functor.hh:132
Definition: transformation.hh:36
Definition: transformation.hh:92
A vector valued function space.
Definition: functionspace.hh:60
convert functions space to space with new dim domain
Definition: functionspace.hh:246
FunctionSpaceTraits::LinearMappingType JacobianRangeType
Intrinsic type used for the jacobian values has a Dune::FieldMatrix type interface.
Definition: functionspaceinterface.hh:75
FunctionSpaceTraits::RangeFieldType RangeFieldType
Intrinsic type used for values in the range field (usually a double)
Definition: functionspaceinterface.hh:63
Interface class for shape function sets.
Definition: shapefunctionset/shapefunctionset.hh:33
void hessianEach(const Point &x, Functor functor) const
evalute hessian of each shape function
void evaluateEach(const Point &x, Functor functor) const
evalute each shape function
std::size_t size() const
return number of shape functions
int order() const
return order of shape functions
void jacobianEach(const Point &x, Functor functor) const
evalute jacobian of each shape function