1 #ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_TUPLE_HH 2 #define DUNE_FEM_SPACE_BASISFUNCTIONSET_TUPLE_HH 3 4 #include <algorithm> 5 #include <array> 6 #include <tuple> 7 #include <utility> 8 #include <vector> 9 10 #include <dune/fem/common/forloop.hh> 11 12 #include <dune/geometry/type.hh> 13 14 #include <dune/fem/common/utility.hh> 15 16 #include <dune/fem/space/basisfunctionset/basisfunctionset.hh> 17 #include <dune/fem/space/combinedspace/subobjects.hh> 18 #include <dune/fem/storage/subvector.hh> 19 20 namespace Dune 21 { 22 23 namespace Fem 24 { 25 26 // TupleBasisFunctionSet 27 // --------------------- 28 // 29 // TupleDiscreteFunctionSpace combination operation 30 // describes the way in which the spaces have been combined, i.e. 31 // Product: V = V_1 x V_2 x ... 32 // Summation: V = V_1 + V_2 + ... 33 34 struct TupleSpaceProduct {}; 35 struct TupleSpaceSummation {}; 36 37 template< class CombineOp, class ... BasisFunctionSets > 38 class TupleBasisFunctionSet 39 { 40 typedef TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... > ThisType; 41 42 // Functor Helper structs 43 template< int > struct ComputeOffset; 44 template< int > struct EvaluateAll; 45 template< int > struct JacobianAll; 46 template< int > struct HessianAll; 47 template< int > struct EvaluateAllRanges; 48 template< int > struct JacobianAllRanges; 49 template< int > struct HessianAllRanges; 50 template< int > struct Axpy; 51 52 // helper class to compute static overall size and offsets 53 template< int ... i > 54 struct Indices 55 { 56 template< int j > sizeDune::Fem::TupleBasisFunctionSet::Indices57 static constexpr int size () { return std::tuple_element< j, std::tuple< std::integral_constant< int, i > ... > >::type::value; } 58 59 template< int ... j > sizesDune::Fem::TupleBasisFunctionSet::Indices60 static constexpr std::integer_sequence< int, size< j >() ... > sizes ( std::integer_sequence< int, j ... > ) 61 { 62 return std::integer_sequence< int, size< j >() ... >(); 63 } 64 65 template< int j > offsetDune::Fem::TupleBasisFunctionSet::Indices66 static constexpr int offset () 67 { 68 return mysum( sizes( std::make_integer_sequence< int, j >() ) ); 69 } 70 71 private: 72 template< int ... j > mysumDune::Fem::TupleBasisFunctionSet::Indices73 static constexpr int mysum ( std::integer_sequence< int, j ... > ) 74 { 75 return Std::sum( j ... ); 76 } 77 mysumDune::Fem::TupleBasisFunctionSet::Indices78 static constexpr int mysum ( std::integer_sequence< int > ) 79 { 80 return 0; 81 } 82 }; 83 84 //! type of tuple 85 typedef std::tuple< BasisFunctionSets ... > BasisFunctionSetTupleType; 86 87 static_assert( Std::are_all_same< typename BasisFunctionSets::DomainType ... >::value, 88 "TupleBasisFunctionSet needs common DomainType" ); 89 //! get type of first function space, to obtain dimDomain and Field Types 90 typedef typename std::tuple_element< 0, BasisFunctionSetTupleType >::type::FunctionSpaceType ContainedFunctionSpaceType; 91 92 //! size of domain space 93 static const int dimDomain = ContainedFunctionSpaceType::dimDomain; 94 95 //! number of BasisFunctionsSets in the tuple 96 static const int setSize = sizeof ... ( BasisFunctionSets ); 97 static const int setIterationSize = sizeof ... ( BasisFunctionSets )-1; 98 99 //! type of offset array 100 typedef std::array< std::size_t, setSize + 1 > OffsetType; 101 102 // storage for i-th sub basisfunction 103 struct Storage 104 { 105 typedef std::tuple< std::vector< typename BasisFunctionSets::RangeType > ... > RangeStorageTupleType; 106 typedef std::tuple< std::vector< typename BasisFunctionSets::JacobianRangeType > ... > JacobianStorageTupleType; 107 typedef std::tuple< std::vector< typename BasisFunctionSets::HessianRangeType > ... > HessianStorageTupleType; 108 StorageDune::Fem::TupleBasisFunctionSet::Storage109 Storage () {} 110 111 template< int i > 112 typename std::tuple_element< i, RangeStorageTupleType >::type rangeStorageDune::Fem::TupleBasisFunctionSet::Storage113 & rangeStorage() const { return std::get< i >( rangeStorage_ ); } 114 115 template< int i > 116 typename std::tuple_element< i, JacobianStorageTupleType >::type jacobianStorageDune::Fem::TupleBasisFunctionSet::Storage117 & jacobianStorage() const { return std::get< i >( jacobianStorage_ ); } 118 119 template< int i > 120 typename std::tuple_element< i, HessianStorageTupleType >::type hessianStorageDune::Fem::TupleBasisFunctionSet::Storage121 & hessianStorage() const { return std::get< i >( hessianStorage_ ); } 122 123 private: 124 mutable RangeStorageTupleType rangeStorage_; 125 mutable JacobianStorageTupleType jacobianStorage_; 126 mutable HessianStorageTupleType hessianStorage_; 127 }; 128 129 130 public: 131 //! helper class to compute static rangeoffsets 132 typedef Indices< BasisFunctionSets::FunctionSpaceType::dimRange ... > RangeIndices; 133 134 protected: 135 template <int dummy, class CombOp> 136 struct CombinationTraits; 137 138 template <int dummy> 139 struct CombinationTraits< dummy, TupleSpaceProduct > 140 { 141 //! size of range space for product space 142 static constexpr const int dimRange = Std::sum( static_cast< int >( BasisFunctionSets::FunctionSpaceType::dimRange ) ... ); 143 144 //! type of analytical combined product function space 145 typedef FunctionSpace< typename ContainedFunctionSpaceType::DomainFieldType, 146 typename ContainedFunctionSpaceType::RangeFieldType, 147 dimDomain, dimRange > FunctionSpaceType; 148 149 template <class ThisRange, class Range> applyDune::Fem::TupleBasisFunctionSet::CombinationTraits150 static void apply(const int rangeOffset, const ThisRange& thisRange, Range& values) 151 { 152 // scatter result to correct place in larger result vector 153 std::copy( thisRange.begin(), thisRange.end(), values.begin() + rangeOffset ); 154 } 155 156 template< int i > rangeOffsetDune::Fem::TupleBasisFunctionSet::CombinationTraits157 static constexpr int rangeOffset () 158 { 159 return RangeIndices::template offset< i >(); 160 } 161 }; 162 163 template <int dummy> 164 struct CombinationTraits< dummy, TupleSpaceSummation > 165 { 166 //! type of analytical is same as contained space since its a sum 167 typedef ContainedFunctionSpaceType FunctionSpaceType; 168 169 template <class ThisRange, class Range> applyDune::Fem::TupleBasisFunctionSet::CombinationTraits170 static void apply(const int rangeOffset, const ThisRange& thisRange, Range& values) 171 { 172 // sum results 173 values += thisRange; 174 } 175 176 template< int j > rangeOffsetDune::Fem::TupleBasisFunctionSet::CombinationTraits177 static constexpr int rangeOffset () 178 { 179 return 0; 180 } 181 }; 182 183 // type of accumulation of result after loop over basis sets 184 typedef CombinationTraits< 0, CombineOp > CombinationType; 185 186 public: 187 // export type of i-th subbasisfunction set 188 template< int i > 189 struct SubBasisFunctionSet 190 { 191 typedef typename std::tuple_element< i, BasisFunctionSetTupleType >::type type; 192 }; 193 194 //! type of analytical combined function space 195 typedef typename CombinationType :: FunctionSpaceType FunctionSpaceType; 196 197 //! size of domain space 198 static const int dimRange = FunctionSpaceType::dimRange; 199 200 //! type of Domain Vector 201 typedef typename FunctionSpaceType::DomainType DomainType; 202 203 //! type of Range Vector 204 typedef typename FunctionSpaceType::RangeType RangeType; 205 206 //! type of Range Vector field 207 typedef typename FunctionSpaceType::RangeFieldType RangeFieldType; 208 209 //! type of Jacobian Vector/Matrix 210 typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType; 211 212 //! type of Hessian Matrix 213 typedef typename FunctionSpaceType::HessianRangeType HessianRangeType; 214 215 static_assert( Std::are_all_same< typename BasisFunctionSets::EntityType ... >::value, 216 "TupleBasisFunctionSet needs common EntityType" ); 217 //! type of Entity the basis function set is initialized on 218 typedef typename std::tuple_element< 0, BasisFunctionSetTupleType >::type::EntityType EntityType; 219 220 static_assert( Std::are_all_same< typename BasisFunctionSets::ReferenceElementType ... >::value, 221 "TupleBasisFunctionSet needs ReferenceElementType" ); 222 //! type of reference element for this BasisFunctionSet 223 typedef typename std::tuple_element< 0, BasisFunctionSetTupleType >::type::ReferenceElementType ReferenceElementType; 224 225 // default constructor, needed by localmatrix TupleBasisFunctionSet()226 TupleBasisFunctionSet () : offset_( {{ 0 }} ) {} 227 228 // constructor taking a pack of basisFunctionSets TupleBasisFunctionSet(const BasisFunctionSets &...basisFunctionSets)229 TupleBasisFunctionSet ( const BasisFunctionSets & ... basisFunctionSets ) 230 : basisFunctionSetTuple_( basisFunctionSets ... ), 231 offset_() 232 { 233 offset_[ 0 ] = 0; 234 Fem::ForLoop< ComputeOffset, 0, setIterationSize >::apply( offset_, basisFunctionSetTuple_ ); 235 } 236 237 // constructor taking a tuple of basisfunction sets TupleBasisFunctionSet(const BasisFunctionSetTupleType & basisFunctionSetTuple)238 TupleBasisFunctionSet ( const BasisFunctionSetTupleType &basisFunctionSetTuple ) 239 : basisFunctionSetTuple_( basisFunctionSetTuple ), 240 offset_() 241 { 242 offset_[ 0 ] = 0; 243 Fem::ForLoop< ComputeOffset, 0, setIterationSize >::apply( offset_, basisFunctionSetTuple_ ); 244 } 245 246 // Basis Function Set Interface Methods 247 // ------------------------------------ 248 249 //! \brief return order of basis function set, maximal order in the tupleset order() const250 int order () const 251 { 252 return order( std::index_sequence_for< BasisFunctionSets ... >() ); 253 } 254 255 //! \brief return size of basis function set size() const256 std::size_t size () const 257 { 258 return size( std::index_sequence_for< BasisFunctionSets ... >() ); 259 } 260 261 //! \copydoc BasisFunctionSet::type type() const262 Dune::GeometryType type () const 263 { 264 return std::get< 0 >( basisFunctionSetTuple_ ).type(); 265 } 266 267 //! \copydoc BasisFunctionSet::valid valid() const268 bool valid () const 269 { 270 return std::get< 0 >( basisFunctionSetTuple_ ).valid(); 271 } 272 273 //! \copydoc BasisFunctionSet::entity entity() const274 const EntityType &entity () const 275 { 276 return std::get< 0 >( basisFunctionSetTuple_ ).entity(); 277 } 278 279 //! \copydoc BasisFunctionSet::entity referenceElement() const280 const ReferenceElementType &referenceElement () const 281 { 282 return std::get< 0 >( basisFunctionSetTuple_ ).referenceElement(); 283 } 284 285 //! \copydoc BasisFunctionSet::evaluateAll( x, dofs, value ) 286 template< class Point, class DofVector > evaluateAll(const Point & x,const DofVector & dofs,RangeType & value) const287 void evaluateAll ( const Point &x, const DofVector &dofs, RangeType &value ) const 288 { 289 Fem::ForLoop< EvaluateAll, 0, setIterationSize >::apply( x, dofs, value, offset_, basisFunctionSetTuple_ ); 290 } 291 292 //! \copydoc BasisFunctionSet::evaluateAll( x, values ) 293 template< class Point, class RangeArray > evaluateAll(const Point & x,RangeArray & values) const294 void evaluateAll ( const Point &x, RangeArray &values ) const 295 { 296 assert( values.size() >= size() ); 297 Fem::ForLoop< EvaluateAllRanges, 0, setIterationSize >::apply( x, values, storage_, offset_, basisFunctionSetTuple_ ); 298 } 299 300 //! \copydoc BasisFunctionSet::evaluateAll( quad, dofs, ranges ) 301 template< class QuadratureType, class DofVector, class RangeArray > evaluateAll(const QuadratureType & quad,const DofVector & dofs,RangeArray & ranges) const302 void evaluateAll ( const QuadratureType &quad, const DofVector &dofs, RangeArray &ranges ) const 303 { 304 const int nop = quad.nop(); 305 for( int qp = 0; qp < nop; ++qp ) 306 evaluateAll( quad[ qp ], dofs, ranges[ qp ] ); 307 } 308 309 //! \copydoc BasisFunctionSet::jacobianAll( x, dofs, jacobian ) 310 template< class Point, class DofVector > jacobianAll(const Point & x,const DofVector & dofs,JacobianRangeType & jacobian) const311 void jacobianAll ( const Point &x, const DofVector &dofs, JacobianRangeType &jacobian ) const 312 { 313 Fem::ForLoop< JacobianAll, 0, setIterationSize >::apply( x, dofs, jacobian, offset_, basisFunctionSetTuple_ ); 314 } 315 316 //! \copydoc BasisFunctionSet::jacobianAll( x, dofs, jacobians ) 317 template< class Point, class JacobianRangeArray > jacobianAll(const Point & x,JacobianRangeArray & jacobians) const318 void jacobianAll ( const Point &x, JacobianRangeArray &jacobians ) const 319 { 320 assert( jacobians.size() >= size() ); 321 Fem::ForLoop< JacobianAllRanges, 0, setIterationSize >::apply( x, jacobians, storage_, offset_, basisFunctionSetTuple_ ); 322 } 323 324 //! \brief evaluate the jacobian of all basis functions and store the result in the jacobians array 325 template< class QuadratureType, class DofVector, class JacobianArray > jacobianAll(const QuadratureType & quad,const DofVector & dofs,JacobianArray & jacobians) const326 void jacobianAll ( const QuadratureType &quad, const DofVector &dofs, JacobianArray &jacobians ) const 327 { 328 const int nop = quad.nop(); 329 for( int qp = 0; qp < nop; ++qp ) 330 jacobianAll( quad[ qp ], dofs, jacobians[ qp ] ); 331 } 332 333 //! \copydoc BasisFunctionSet::hessianAll( x, dofs, hessian ) 334 template< class Point, class DofVector > hessianAll(const Point & x,const DofVector & dofs,HessianRangeType & hessian) const335 void hessianAll ( const Point &x, const DofVector &dofs, HessianRangeType &hessian ) const 336 { 337 Fem::ForLoop< HessianAll, 0, setIterationSize >::apply( x, dofs, hessian, offset_, basisFunctionSetTuple_ ); 338 } 339 340 //! \brief evaluate the hessian of all basis functions and store the result in the hessians array 341 template< class QuadratureType, class DofVector, class HessianArray > hessianAll(const QuadratureType & quad,const DofVector & dofs,HessianArray & hessians) const342 void hessianAll ( const QuadratureType &quad, const DofVector &dofs, HessianArray &hessians ) const 343 { 344 const int nop = quad.nop(); 345 for( int qp = 0; qp < nop; ++qp ) 346 hessianAll( quad[ qp ], dofs, hessians[ qp ] ); 347 } 348 349 //! \copydoc BasisFunctionSet::hessianAll( x, hessians ) 350 template< class Point, class HessianRangeArray > hessianAll(const Point & x,HessianRangeArray & hessians) const351 void hessianAll ( const Point &x, HessianRangeArray &hessians ) const 352 { 353 assert( hessians.size() >= size() ); 354 Fem::ForLoop< HessianAllRanges, 0, setIterationSize >::apply( x, hessians, storage_, offset_, basisFunctionSetTuple_ ); 355 } 356 357 //! \copydoc BasisFunctionSet::axpy( quad, values, dofs ) 358 template< class QuadratureType, class Vector, class DofVector > axpy(const QuadratureType & quad,const Vector & values,DofVector & dofs) const359 void axpy ( const QuadratureType &quad, const Vector &values, DofVector &dofs ) const 360 { 361 // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector 362 const unsigned int nop = quad.nop(); 363 for( unsigned int qp = 0; qp < nop; ++qp ) 364 axpy( quad[ qp ], values[ qp ], dofs ); 365 } 366 367 //! \copydoc BasisFunctionSet::axpy( quad, valuesA, valuesB, dofs ) 368 template< class QuadratureType, class VectorA, class VectorB, class DofVector > axpy(const QuadratureType & quad,const VectorA & valuesA,const VectorB & valuesB,DofVector & dofs) const369 void axpy ( const QuadratureType &quad, const VectorA &valuesA, const VectorB &valuesB, DofVector &dofs ) const 370 { 371 // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector 372 const unsigned int nop = quad.nop(); 373 for( unsigned int qp = 0; qp < nop; ++qp ) 374 { 375 axpy( quad[ qp ], valuesA[ qp ], dofs ); 376 axpy( quad[ qp ], valuesB[ qp ], dofs ); 377 } 378 } 379 380 //! \copydoc BasisFunctionSet::axpy( x, valueFactor, dofs ) 381 template< class Point, class DofVector > axpy(const Point & x,const RangeType & valueFactor,DofVector & dofs) const382 void axpy ( const Point &x, const RangeType &valueFactor, DofVector &dofs ) const 383 { 384 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, valueFactor, dofs, offset_, basisFunctionSetTuple_ ); 385 } 386 387 //! \copydoc BasisFunctionSet::axpy( x, jacobianFactor, dofs ) 388 template< class Point, class DofVector > axpy(const Point & x,const JacobianRangeType & jacobianFactor,DofVector & dofs) const389 void axpy ( const Point &x, const JacobianRangeType &jacobianFactor, DofVector &dofs ) const 390 { 391 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, jacobianFactor, dofs, offset_, basisFunctionSetTuple_ ); 392 } 393 394 //! \copydoc BasisFunctionSet::axpy( x, hessianFactor, dofs ) 395 template< class Point, class DofVector > axpy(const Point & x,const HessianRangeType & hessianFactor,DofVector & dofs) const396 void axpy ( const Point &x, const HessianRangeType &hessianFactor, DofVector &dofs ) const 397 { 398 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, hessianFactor, dofs, offset_, basisFunctionSetTuple_ ); 399 } 400 401 //! \copydoc BasisFunctionSet::axpy( x, valueFactor, jacobianFactor, dofs ) 402 template< class Point, class DofVector > axpy(const Point & x,const RangeType & valueFactor,const JacobianRangeType & jacobianFactor,DofVector & dofs) const403 void axpy ( const Point &x, const RangeType &valueFactor, const JacobianRangeType &jacobianFactor, DofVector &dofs ) const 404 { 405 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, valueFactor, jacobianFactor, dofs, offset_, basisFunctionSetTuple_ ); 406 } 407 408 /***** NON Interface methods ****/ 409 410 //! return i-th subbasisfunctionSet 411 template< int i > subBasisFunctionSet() const412 const typename SubBasisFunctionSet< i >::type& subBasisFunctionSet() const 413 { 414 return std::get< i >( basisFunctionSetTuple_ ); 415 } 416 417 //! return offset of the i-th subbasisfunctionSet in the whole set offset(int i) const418 std::size_t offset ( int i ) const 419 { 420 return offset_[ i ]; 421 } 422 423 //! return number of subBasisFunctionSets numSubBasisFunctionSets()424 static int numSubBasisFunctionSets () 425 { 426 return setSize; 427 } 428 429 protected: 430 // unroll index sequence and take maximal order 431 template< std::size_t ... i > order(std::index_sequence<i...>) const432 int order ( std::index_sequence< i ... > ) const 433 { 434 return Std::max( std::get< i >( basisFunctionSetTuple_ ).order() ... ); 435 } 436 437 // unroll index sequence and sum up sizes 438 template< std::size_t ... i > size(std::index_sequence<i...>) const439 std::size_t size ( std::index_sequence< i ... > ) const 440 { 441 return Std::sum( std::get< i >( basisFunctionSetTuple_ ).size() ... ); 442 } 443 444 private: 445 BasisFunctionSetTupleType basisFunctionSetTuple_; 446 OffsetType offset_; 447 448 Storage storage_; 449 }; 450 451 452 453 // ComputeOffset 454 // ------------- 455 456 template< class CombineOp, class ... BasisFunctionSets > 457 template< int i > 458 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >:: 459 ComputeOffset 460 { 461 template< class Tuple > applyDune::Fem::TupleBasisFunctionSet::ComputeOffset462 static void apply ( OffsetType &offset, const Tuple &tuple ) 463 { 464 offset[ i + 1 ] = offset[ i ] + std::get< i >( tuple ).size(); 465 } 466 }; 467 468 469 // EvaluateAll 470 // ----------- 471 472 template< class CombineOp, class ... BasisFunctionSets > 473 template< int i > 474 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >:: 475 EvaluateAll 476 { 477 // only needed for product spaces 478 static const int rangeOffset = CombinationType :: template rangeOffset< i >(); 479 480 template< class Point, class DofVector, class Tuple > 481 static void apply ( const Point &x, const DofVector &dofVector, RangeType &values, const OffsetType &offset, const Tuple &tuple ) 482 { 483 // evaluateAll for this BasisFunctionSet, with View on DofVector 484 std::size_t size = std::get< i >( tuple ).size(); 485 SubVector< const DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper( size, offset[ i ] ) ); 486 typename std::tuple_element< i, BasisFunctionSetTupleType >::type::RangeType thisRange; 487 std::get< i >( tuple ).evaluateAll( x, subDofVector, thisRange ); 488 489 // distribute result to values 490 CombinationType::apply( rangeOffset, thisRange, values ); 491 } 492 }; 493 494 495 // JacobianAll 496 // ----------- 497 498 template< class CombineOp, class ... BasisFunctionSets > 499 template< int i > 500 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >:: 501 JacobianAll 502 { 503 // only needed for product spaces 504 static const int rangeOffset = CombinationType :: template rangeOffset< i >(); 505 506 template< class Point, class DofVector, class Tuple > 507 static void apply ( const Point &x, const DofVector &dofVector, JacobianRangeType &values, const OffsetType &offset, const Tuple &tuple ) 508 { 509 // jacobianAll for this BasisFunctionSet, with View on DofVector 510 std::size_t size = std::get< i >( tuple ).size(); 511 SubVector< const DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper( size, offset[ i ] ) ); 512 typename std::tuple_element< i, BasisFunctionSetTupleType >::type::JacobianRangeType thisJacobian; 513 std::get< i >( tuple ).jacobianAll( x, subDofVector, thisJacobian ); 514 515 // distribute result to values 516 CombinationType::apply( rangeOffset, thisJacobian, values ); 517 } 518 }; 519 520 521 // HessianAll 522 // ---------- 523 524 template< class CombineOp, class ... BasisFunctionSets > 525 template< int i > 526 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >:: 527 HessianAll 528 { 529 // only needed for product spaces 530 static const int rangeOffset = CombinationType :: template rangeOffset< i >(); 531 532 template< class Point, class DofVector, class Tuple > 533 static void apply ( const Point &x, const DofVector &dofVector, HessianRangeType &values, const OffsetType &offset, const Tuple &tuple ) 534 { 535 // hessianAll for this BasisFunctionSet, with View on DofVector 536 std::size_t size = std::get< i >( tuple ).size(); 537 SubVector< const DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper( size, offset[ i ] ) ); 538 typename std::tuple_element< i, BasisFunctionSetTupleType >::type::HessianRangeType thisHessian; 539 std::get< i >( tuple ).hessianAll( x, subDofVector, thisHessian ); 540 541 // distribute result to values 542 CombinationType::apply( rangeOffset, thisHessian, values ); 543 } 544 }; 545 546 547 // EvaluateAllRanges 548 // ----------------- 549 550 template< class CombineOp, class ... BasisFunctionSets > 551 template< int i > 552 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >:: 553 EvaluateAllRanges 554 { 555 typedef typename std::tuple_element< i, typename Storage::RangeStorageTupleType >::type ThisStorage; 556 static const int thisDimRange = std::tuple_element< i, BasisFunctionSetTupleType >::type::FunctionSpaceType::dimRange; 557 558 // only needed for product spaces 559 static const int rangeOffset = CombinationType :: template rangeOffset< i >(); 560 561 template< class Point, class RangeArray, class Tuple > 562 static void apply ( const Point &x, RangeArray &values, const Storage &s, const OffsetType &offset, const Tuple &tuple ) 563 { 564 std::size_t size = std::get< i >( tuple ).size(); 565 ThisStorage &thisStorage = s.template rangeStorage< i >(); 566 thisStorage.resize( size ); 567 568 std::get< i >( tuple ).evaluateAll( x, thisStorage ); 569 570 for( std::size_t j = 0; j < size; ++j ) 571 { 572 values[ j + offset[ i ] ] = RangeType( 0.0 ); 573 for( int r = 0; r < thisDimRange; ++r ) 574 values[ j + offset[ i ] ][ r + rangeOffset ] = thisStorage[ j ][ r ]; 575 } 576 } 577 }; 578 579 580 // JacobianAllRanges 581 // ----------------- 582 583 template< class CombineOp, class ... BasisFunctionSets > 584 template< int i > 585 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >:: 586 JacobianAllRanges 587 { 588 typedef typename std::tuple_element< i, typename Storage::JacobianStorageTupleType >::type ThisStorage; 589 static const int thisDimRange = std::tuple_element< i, BasisFunctionSetTupleType >::type::FunctionSpaceType::dimRange; 590 591 // only needed for product spaces 592 static const int rangeOffset = CombinationType :: template rangeOffset< i >(); 593 594 template< class Point, class RangeArray, class Tuple > 595 static void apply ( const Point &x, RangeArray &values, const Storage &s, const OffsetType &offset, const Tuple &tuple ) 596 { 597 std::size_t size = std::get< i >( tuple ).size(); 598 ThisStorage &thisStorage = s.template jacobianStorage< i >(); 599 thisStorage.resize( size ); 600 601 std::get< i >( tuple ).jacobianAll( x, thisStorage ); 602 603 for( std::size_t j = 0; j < size; ++j ) 604 { 605 values[ j + offset[ i ] ] = JacobianRangeType( RangeFieldType( 0.0 ) ); 606 for( int r = 0; r < thisDimRange; ++r ) 607 values[ j + offset[ i ] ][ r + rangeOffset ] = thisStorage[ j ][ r ]; 608 } 609 } 610 }; 611 612 613 // HessianAllRanges 614 // ---------------- 615 616 template< class CombineOp, class ... BasisFunctionSets > 617 template< int i > 618 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >:: 619 HessianAllRanges 620 { 621 typedef typename std::tuple_element< i, typename Storage::HessianStorageTupleType >::type ThisStorage; 622 static const int thisDimRange = std::tuple_element< i, BasisFunctionSetTupleType >::type::FunctionSpaceType::dimRange; 623 624 // only needed for product spaces 625 static const int rangeOffset = CombinationType :: template rangeOffset< i >(); 626 627 template< class Point, class RangeArray, class Tuple > 628 static void apply ( const Point &x, RangeArray &values, const Storage &s, const OffsetType &offset, const Tuple &tuple ) 629 { 630 std::size_t size = std::get< i >( tuple ).size(); 631 ThisStorage &thisStorage = s.template hessianStorage< i >(); 632 thisStorage.resize( size ); 633 634 std::get< i >( tuple ).hessianAll( x, thisStorage ); 635 636 for( std::size_t j = 0; j < size; ++j ) 637 { 638 values[ j + offset[ i ] ] = typename HessianRangeType::value_type( 0.0 ); 639 for( int r = 0; r < thisDimRange; ++r ) 640 values[ j + offset[ i ] ][ r + rangeOffset ] = thisStorage[ j ][ r ]; 641 } 642 } 643 }; 644 645 646 // Axpy 647 // ---- 648 649 template< class CombineOp, class ... BasisFunctionSets > 650 template< int i > 651 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >:: 652 Axpy 653 { 654 typedef typename std::tuple_element< i, BasisFunctionSetTupleType >::type::RangeType ThisRangeType; 655 typedef typename std::tuple_element< i, BasisFunctionSetTupleType >::type::JacobianRangeType ThisJacobianRangeType; 656 typedef typename std::tuple_element< i, BasisFunctionSetTupleType >::type::HessianRangeType ThisHessianRangeType; 657 658 // only needed for product spaces 659 static const int rangeOffset = CombinationType :: template rangeOffset< i >(); 660 661 typedef SubObject< const RangeType, const ThisRangeType, rangeOffset > SubRangeType; 662 typedef SubObject< const JacobianRangeType, const ThisJacobianRangeType, rangeOffset > SubJacobianRangeType; 663 typedef SubObject< const HessianRangeType, const ThisHessianRangeType, rangeOffset > SubHessianRangeType; 664 665 // axpy with range type 666 template< class Point, class DofVector, class Tuple > applyDune::Fem::TupleBasisFunctionSet::Axpy667 static void apply ( const Point &x, const RangeType &factor, DofVector &dofVector, const OffsetType &offset, const Tuple &tuple ) 668 { 669 std::size_t size = std::get< i >( tuple ).size(); 670 SubVector< DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper( size, offset[ i ] ) ); 671 SubRangeType subFactor( factor ); 672 std::get< i >( tuple ).axpy( x, (ThisRangeType) subFactor, subDofVector ); 673 } 674 675 // axpy with jacobian range type 676 template< class Point, class DofVector, class Tuple > applyDune::Fem::TupleBasisFunctionSet::Axpy677 static void apply ( const Point &x, const JacobianRangeType &factor, DofVector &dofVector, const OffsetType &offset, const Tuple &tuple ) 678 { 679 std::size_t size = std::get< i >( tuple ).size(); 680 SubVector< DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper( size, offset[ i ] ) ); 681 SubJacobianRangeType subFactor( factor ); 682 std::get< i >( tuple ).axpy( x, (ThisJacobianRangeType) subFactor, subDofVector ); 683 } 684 685 // axpy with hessian range type 686 template< class Point, class DofVector, class Tuple > applyDune::Fem::TupleBasisFunctionSet::Axpy687 static void apply ( const Point &x, const HessianRangeType &factor, DofVector &dofVector, const OffsetType &offset, const Tuple &tuple ) 688 { 689 std::size_t size = std::get< i >( tuple ).size(); 690 SubVector< DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper( size, offset[ i ] ) ); 691 SubHessianRangeType subFactor( factor ); 692 std::get< i >( tuple ).axpy( x, (ThisHessianRangeType) subFactor, subDofVector ); 693 } 694 695 // axpy with range and jacobian range type 696 template< class Point, class DofVector, class Tuple > applyDune::Fem::TupleBasisFunctionSet::Axpy697 static void apply ( const Point &x, const RangeType &rangeFactor, const JacobianRangeType &jacobianFactor, DofVector &dofVector, const OffsetType &offset, const Tuple &tuple ) 698 { 699 std::size_t size = std::get< i >( tuple ).size(); 700 SubVector< DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper( size, offset[ i ] ) ); 701 SubRangeType subRangeFactor( rangeFactor ); 702 SubJacobianRangeType subJacobianFactor( jacobianFactor ); 703 std::get< i >( tuple ).axpy( x, (ThisRangeType) subRangeFactor, (ThisJacobianRangeType) subJacobianFactor, subDofVector ); 704 } 705 }; 706 707 } // namespace Fem 708 709 } // namespace Dune 710 711 #endif // #ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_TUPLE_HH 712