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 15 #include <dune/fem/space/shapefunctionset/caching.hh> 16 17 // dune-fem includes 18 #include <dune/fem/space/basisfunctionset/functor.hh> 19 #include <dune/fem/space/basisfunctionset/transformation.hh> 20 #include <dune/fem/space/shapefunctionset/caching.hh> 21 #include <dune/fem/quadrature/cachingpointlist.hh> 22 #include <dune/fem/space/common/functionspace.hh> 23 #include <dune/fem/version.hh> 24 25 #include <dune/fem/space/basisfunctionset/codegen.hh> 26 #include <dune/fem/space/basisfunctionset/evaluatecaller.hh> 27 28 namespace Dune 29 { 30 31 namespace Fem 32 { 33 34 // DefaultBasisFunctionSet 35 // ----------------------- 36 37 /** 38 * \brief implementation of a basis function set for given entity 39 * 40 * \tparam Entity entity type 41 * \tparam ShapeFunctionSet shape function set 42 * 43 * \note ShapeFunctionSet must be a copyable object. For most 44 * non-trivial implementations, you may want to use a 45 * proxy, see file 46 \code 47 <dune/fem/space/shapefunctionset/proxy.hh> 48 \endcode 49 */ 50 template< class Entity, class ShapeFunctionSet > 51 class DefaultBasisFunctionSet 52 { 53 typedef DefaultBasisFunctionSet< Entity, ShapeFunctionSet > ThisType; 54 55 public: 56 //! \brief entity type 57 typedef Entity EntityType; 58 //! \brief shape function set type 59 typedef ShapeFunctionSet ShapeFunctionSetType; 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: 66 typedef typename ShapeFunctionSetType::FunctionSpaceType LocalFunctionSpaceType; 67 typedef typename LocalFunctionSpaceType::JacobianRangeType LocalJacobianRangeType; 68 typedef typename LocalFunctionSpaceType::HessianRangeType LocalHessianRangeType; 69 70 typedef typename LocalFunctionSpaceType::RangeFieldType RangeFieldType; 71 72 typedef typename EntityType::Geometry Geometry ; 73 74 typedef typename Geometry::ctype ctype; 75 public: 76 // slight misuse of struct ToLocalFunctionSpace!!! 77 //! \brief type of function space 78 typedef typename ToNewDimDomainFunctionSpace< LocalFunctionSpaceType, Geometry::coorddimension > :: Type FunctionSpaceType; 79 80 //! \brief range type 81 typedef typename FunctionSpaceType::RangeType RangeType; 82 //! \brief domain type 83 typedef typename FunctionSpaceType::DomainType DomainType; 84 //! \brief jacobian range type 85 typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType; 86 //! \brief hessian range type 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 94 //! \brief type of reference element 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 99 //! \brief constructor DefaultBasisFunctionSet()100 DefaultBasisFunctionSet () {} 101 102 //! \brief constructor DefaultBasisFunctionSet(const EntityType & entity,const ShapeFunctionSet & shapeFunctionSet=ShapeFunctionSet ())103 explicit DefaultBasisFunctionSet ( const EntityType &entity, const ShapeFunctionSet &shapeFunctionSet = ShapeFunctionSet() ) 104 : entity_( &entity ), 105 shapeFunctionSet_( shapeFunctionSet ) 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 DefaultBasisFunctionSet(const DefaultBasisFunctionSet & other)113 DefaultBasisFunctionSet ( const DefaultBasisFunctionSet &other ) 114 : entity_( other.entity_ ), 115 shapeFunctionSet_( other.shapeFunctionSet_ ) 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 operator =(const DefaultBasisFunctionSet & other)124 DefaultBasisFunctionSet &operator= ( const DefaultBasisFunctionSet &other ) 125 { 126 entity_ = other.entity_; 127 shapeFunctionSet_ = other.shapeFunctionSet_; 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 140 //! \brief return order of basis function set order() const141 int order () const { return shapeFunctionSet().order(); } 142 143 //! \brief return size of basis function set size() const144 std::size_t size () const { return shapeFunctionSet().size(); } 145 146 //! \brief return reference element referenceElement() const147 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 153 /** \brief evaluate all basis function and multiply with given 154 * values and add to dofs 155 */ 156 template< class QuadratureType, class Vector, class DofVector > axpy(const QuadratureType & quad,const Vector & values,DofVector & dofs) const157 void axpy ( const QuadratureType &quad, const Vector &values, DofVector &dofs ) const 158 { 159 axpyImpl( quad, values, dofs, values[ 0 ] ); 160 } 161 162 /** \brief evaluate all basis function and multiply with given 163 * values and add to dofs 164 * 165 * \note valuesA and valuesB can be vectors of RangeType or 166 * JacobianRangeType 167 */ 168 template< class QuadratureType, class VectorA, class VectorB, class DofVector > axpy(const QuadratureType & quad,const VectorA & valuesA,const VectorB & valuesB,DofVector & dofs) const169 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 178 /** \brief evaluate all basis function and multiply with given 179 * values and add to dofs 180 */ 181 template< class Point, class DofVector > axpy(const Point & x,const RangeType & valueFactor,DofVector & dofs) const182 void axpy ( const Point &x, const RangeType &valueFactor, DofVector &dofs ) const 183 { 184 FunctionalAxpyFunctor< RangeType, DofVector > f( valueFactor, dofs ); 185 shapeFunctionSet().evaluateEach( x, f ); 186 } 187 188 /** \brief evaluate all derivatives of all basis function and multiply with given 189 * values and add to dofs 190 */ 191 template< class Point, class DofVector > axpy(const Point & x,const JacobianRangeType & jacobianFactor,DofVector & dofs) const192 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 201 FunctionalAxpyFunctor< LocalJacobianRangeType, DofVector > f( tmpJacobianFactor, dofs ); 202 shapeFunctionSet().jacobianEach( x, f ); 203 } 204 205 /** \brief evaluate all basis function and derivatives and multiply with given 206 * values and add to dofs 207 */ 208 template< class Point, class DofVector > axpy(const Point & x,const RangeType & valueFactor,const JacobianRangeType & jacobianFactor,DofVector & dofs) const209 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 216 /** \brief Add H:D^2phi to each dof 217 */ 218 template< class Point, class DofVector > axpy(const Point & x,const HessianRangeType & hessianFactor,DofVector & dofs) const219 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 } 238 FunctionalAxpyFunctor< LocalHessianRangeType, DofVector > f( tmpHessianFactor, dofs ); 239 shapeFunctionSet().hessianEach( x, f ); 240 } 241 242 243 /** \copydoc BasisFunctionSet::evaluateAll( quad, dofs, ranges ) */ 244 template< class QuadratureType, class DofVector, class RangeArray > evaluateAll(const QuadratureType & quad,const DofVector & dofs,RangeArray & ranges) const245 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 276 //! \todo please doc me 277 template< class Point, class DofVector > evaluateAll(const Point & x,const DofVector & dofs,RangeType & value) const278 void evaluateAll ( const Point &x, const DofVector &dofs, RangeType &value ) const 279 { 280 value = RangeType( 0 ); 281 AxpyFunctor< DofVector, RangeType > f( dofs, value ); 282 shapeFunctionSet().evaluateEach( x, f ); 283 } 284 285 //! \todo please doc me 286 template< class Point, class RangeArray > evaluateAll(const Point & x,RangeArray & values) const287 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 ); 292 shapeFunctionSet().evaluateEach( x, f ); 293 } 294 295 /** \copydoc BasisFunctionSet::jacobianAll( quad, dofs, jacobians ) */ 296 template< class QuadratureType, class DofVector, class JacobianArray > jacobianAll(const QuadratureType & quad,const DofVector & dofs,JacobianArray & jacobians) const297 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 330 //! \todo please doc me 331 template< class Point, class DofVector > jacobianAll(const Point & x,const DofVector & dofs,JacobianRangeType & jacobian) const332 void jacobianAll ( const Point &x, const DofVector &dofs, JacobianRangeType &jacobian ) const 333 { 334 LocalJacobianRangeType localJacobian( RangeFieldType( 0 ) ); 335 AxpyFunctor< DofVector, LocalJacobianRangeType > f( dofs, localJacobian ); 336 shapeFunctionSet().jacobianEach( x, f ); 337 338 typedef JacobianTransformation< Geometry > Transformation; 339 const Geometry &geo = geometry(); 340 Transformation transformation( geo, coordinate( x ) ); 341 transformation( localJacobian, jacobian ); 342 } 343 344 //! \todo please doc me 345 template< class Point, class JacobianRangeArray > jacobianAll(const Point & x,JacobianRangeArray & jacobians) const346 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 ); 354 shapeFunctionSet().jacobianEach( x, f ); 355 } 356 357 //! \todo please doc me 358 template< class QuadratureType, class DofVector, class HessianArray > hessianAll(const QuadratureType & quad,const DofVector & dofs,HessianArray & hessians) const359 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 370 //! \todo please doc me 371 template< class Point, class DofVector > hessianAll(const Point & x,const DofVector & dofs,HessianRangeType & hessian) const372 void hessianAll ( const Point &x, const DofVector &dofs, HessianRangeType &hessian ) const 373 { 374 LocalHessianRangeType localHessian( typename LocalHessianRangeType::value_type( RangeFieldType( 0 ) ) ); 375 AxpyFunctor< DofVector, LocalHessianRangeType > f( dofs, localHessian ); 376 shapeFunctionSet().hessianEach( x, f ); 377 378 typedef HessianTransformation< Geometry > Transformation; 379 const Geometry &geo = geometry(); 380 Transformation transformation( geo, coordinate( x ) ); 381 transformation( localHessian, hessian ); 382 } 383 384 //! \todo please doc me 385 template< class Point, class HessianRangeArray > hessianAll(const Point & x,HessianRangeArray & hessians) const386 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 ); 393 shapeFunctionSet().hessianEach( x, f ); 394 } 395 396 //! \brief return entity entity() const397 const Entity &entity () const 398 { 399 assert( valid() ); 400 return *entity_; 401 } 402 403 //! \brief return true if entity pointer is set valid() const404 bool valid () const { return bool(entity_); } 405 406 //! \brief return geometry type type() const407 Dune::GeometryType type () const { return entity().type(); } 408 409 // Non-interface methods 410 // --------------------- 411 412 //! \brief return shape function set shapeFunctionSet() const413 const ShapeFunctionSetType &shapeFunctionSet () const { return shapeFunctionSet_; } 414 415 protected: 416 //! \brief evaluate all basis function and multiply with given values and add to dofs 417 template< class QuadratureType, class RangeArray, class DofVector > axpyImpl(const QuadratureType & quad,const RangeArray & rangeFactors,DofVector & dofs,const RangeType &) const418 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 449 //! \brief evaluate all basis function and multiply with given values and add to dofs 450 template< class QuadratureType, class JacobianArray, class DofVector > axpyImpl(const QuadratureType & quad,const JacobianArray & jacobianFactors,DofVector & dofs,const JacobianRangeType &) const451 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> rangeCache(const QuadratureType & quad) const484 const auto& rangeCache( const QuadratureType& quad ) const 485 { 486 return shapeFunctionSet().scalarShapeFunctionSet().impl().rangeCache( quad ); 487 } 488 489 template <class QuadratureType> jacobianCache(const QuadratureType & quad) const490 const auto& jacobianCache( const QuadratureType& quad ) const 491 { 492 return shapeFunctionSet().scalarShapeFunctionSet().impl().jacobianCache( quad ); 493 } 494 495 protected: geometry() const496 Geometry geometry () const { return geometry_.value(); } 497 498 const EntityType *entity_ = nullptr; 499 ShapeFunctionSetType shapeFunctionSet_; 500 std::optional< Geometry > geometry_; 501 }; 502 503 } // namespace Fem 504 505 } // namespace Dune 506 507 #endif // #ifndef DUNE_FEM_BASISFUNCTIONSET_DEFAULT_HH 508