1#ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_TUPLE_HH
2#define DUNE_FEM_SPACE_BASISFUNCTIONSET_TUPLE_HH
10#include <dune/fem/common/forloop.hh>
14#include <dune/fem/common/utility.hh>
16#include <dune/fem/space/basisfunctionset/basisfunctionset.hh>
17#include <dune/fem/space/combinedspace/subobjects.hh>
18#include <dune/fem/storage/subvector.hh>
34 struct TupleSpaceProduct {};
35 struct TupleSpaceSummation {};
37 template<
class CombineOp,
class ... BasisFunctionSets >
38 class TupleBasisFunctionSet
40 typedef TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... > ThisType;
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;
57 static constexpr int size () {
return std::tuple_element< j, std::tuple< std::integral_constant< int, i > ... > >::type::value; }
60 static constexpr std::integer_sequence< int, size< j >() ... > sizes ( std::integer_sequence< int, j ... > )
62 return std::integer_sequence< int, size< j >() ... >();
66 static constexpr int offset ()
68 return mysum( sizes( std::make_integer_sequence< int, j >() ) );
73 static constexpr int mysum ( std::integer_sequence< int, j ... > )
75 return Std::sum( j ... );
78 static constexpr int mysum ( std::integer_sequence< int > )
85 typedef std::tuple< BasisFunctionSets ... > BasisFunctionSetTupleType;
87 static_assert( Std::are_all_same<
typename BasisFunctionSets::DomainType ... >::value,
88 "TupleBasisFunctionSet needs common DomainType" );
90 typedef typename std::tuple_element< 0, BasisFunctionSetTupleType >::type::FunctionSpaceType ContainedFunctionSpaceType;
93 static const int dimDomain = ContainedFunctionSpaceType::dimDomain;
96 static const int setSize =
sizeof ... ( BasisFunctionSets );
97 static const int setIterationSize =
sizeof ... ( BasisFunctionSets )-1;
100 typedef std::array< std::size_t, setSize + 1 > OffsetType;
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;
112 typename std::tuple_element< i, RangeStorageTupleType >::type
113 & rangeStorage()
const {
return std::get< i >( rangeStorage_ ); }
116 typename std::tuple_element< i, JacobianStorageTupleType >::type
117 & jacobianStorage()
const {
return std::get< i >( jacobianStorage_ ); }
120 typename std::tuple_element< i, HessianStorageTupleType >::type
121 & hessianStorage()
const {
return std::get< i >( hessianStorage_ ); }
124 mutable RangeStorageTupleType rangeStorage_;
125 mutable JacobianStorageTupleType jacobianStorage_;
126 mutable HessianStorageTupleType hessianStorage_;
132 typedef Indices< BasisFunctionSets::FunctionSpaceType::dimRange ... > RangeIndices;
135 template <
int dummy,
class CombOp>
136 struct CombinationTraits;
139 struct CombinationTraits< dummy, TupleSpaceProduct >
142 static constexpr const int dimRange = Std::sum(
static_cast< int >( BasisFunctionSets::FunctionSpaceType::dimRange ) ... );
145 typedef FunctionSpace<
typename ContainedFunctionSpaceType::DomainFieldType,
146 typename ContainedFunctionSpaceType::RangeFieldType,
149 template <
class ThisRange,
class Range>
150 static void apply(
const int rangeOffset,
const ThisRange& thisRange, Range& values)
153 std::copy( thisRange.begin(), thisRange.end(), values.begin() + rangeOffset );
157 static constexpr int rangeOffset ()
159 return RangeIndices::template offset< i >();
164 struct CombinationTraits< dummy, TupleSpaceSummation >
169 template <
class ThisRange,
class Range>
170 static void apply(
const int rangeOffset,
const ThisRange& thisRange, Range& values)
177 static constexpr int rangeOffset ()
184 typedef CombinationTraits< 0, CombineOp > CombinationType;
189 struct SubBasisFunctionSet
191 typedef typename std::tuple_element< i, BasisFunctionSetTupleType >::type type;
216 "TupleBasisFunctionSet needs common EntityType" );
218 typedef typename std::tuple_element< 0, BasisFunctionSetTupleType >::type::EntityType EntityType;
220 static_assert( Std::are_all_same<
typename BasisFunctionSets::Geometry ... >::value,
221 "TupleBasisFunctionSet needs common Geometry" );
223 typedef typename std::tuple_element< 0, BasisFunctionSetTupleType >::type::Geometry Geometry;
225 static_assert( Std::are_all_same<
typename BasisFunctionSets::ReferenceElementType ... >::value,
226 "TupleBasisFunctionSet needs ReferenceElementType" );
228 typedef typename std::tuple_element< 0, BasisFunctionSetTupleType >::type::ReferenceElementType ReferenceElementType;
231 TupleBasisFunctionSet () : offset_( {{ 0 }} ) {}
234 TupleBasisFunctionSet (
const BasisFunctionSets & ... basisFunctionSets )
235 : basisFunctionSetTuple_( basisFunctionSets ... ),
239 Fem::ForLoop< ComputeOffset, 0, setIterationSize >::apply( offset_, basisFunctionSetTuple_ );
243 TupleBasisFunctionSet (
const BasisFunctionSetTupleType &basisFunctionSetTuple )
244 : basisFunctionSetTuple_( basisFunctionSetTuple ),
248 Fem::ForLoop< ComputeOffset, 0, setIterationSize >::apply( offset_, basisFunctionSetTuple_ );
257 return order( std::index_sequence_for< BasisFunctionSets ... >() );
261 std::size_t
size ()
const
263 return size( std::index_sequence_for< BasisFunctionSets ... >() );
269 return std::get< 0 >( basisFunctionSetTuple_ ).type();
275 return std::get< 0 >( basisFunctionSetTuple_ ).valid();
279 const EntityType &entity ()
const
281 return std::get< 0 >( basisFunctionSetTuple_ ).entity();
285 const Geometry &geometry ()
const
287 return std::get< 0 >( basisFunctionSetTuple_ ).geometry();
293 return std::get< 0 >( basisFunctionSetTuple_ ).referenceElement();
297 template<
class Po
int,
class DofVector >
298 void evaluateAll (
const Point &x,
const DofVector &dofs, RangeType &value )
const
300 Fem::ForLoop< EvaluateAll, 0, setIterationSize >::apply( x, dofs, value, offset_, basisFunctionSetTuple_ );
304 template<
class Po
int,
class RangeArray >
305 void evaluateAll (
const Point &x, RangeArray &values )
const
307 assert( values.size() >=
size() );
308 Fem::ForLoop< EvaluateAllRanges, 0, setIterationSize >::apply( x, values, storage_, offset_, basisFunctionSetTuple_ );
312 template<
class QuadratureType,
class DofVector,
class RangeArray >
313 void evaluateAll (
const QuadratureType &quad,
const DofVector &dofs, RangeArray &ranges )
const
315 const int nop = quad.nop();
316 for(
int qp = 0; qp < nop; ++qp )
317 evaluateAll( quad[ qp ], dofs, ranges[ qp ] );
321 template<
class Po
int,
class DofVector >
322 void jacobianAll (
const Point &x,
const DofVector &dofs, JacobianRangeType &jacobian )
const
324 Fem::ForLoop< JacobianAll, 0, setIterationSize >::apply( x, dofs, jacobian, offset_, basisFunctionSetTuple_ );
328 template<
class Po
int,
class JacobianRangeArray >
329 void jacobianAll (
const Point &x, JacobianRangeArray &jacobians )
const
331 assert( jacobians.size() >=
size() );
332 Fem::ForLoop< JacobianAllRanges, 0, setIterationSize >::apply( x, jacobians, storage_, offset_, basisFunctionSetTuple_ );
336 template<
class QuadratureType,
class DofVector,
class JacobianArray >
337 void jacobianAll (
const QuadratureType &quad,
const DofVector &dofs, JacobianArray &jacobians )
const
339 const int nop = quad.nop();
340 for(
int qp = 0; qp < nop; ++qp )
341 jacobianAll( quad[ qp ], dofs, jacobians[ qp ] );
345 template<
class Po
int,
class DofVector >
346 void hessianAll (
const Point &x,
const DofVector &dofs, HessianRangeType &hessian )
const
348 Fem::ForLoop< HessianAll, 0, setIterationSize >::apply( x, dofs, hessian, offset_, basisFunctionSetTuple_ );
352 template<
class QuadratureType,
class DofVector,
class HessianArray >
353 void hessianAll (
const QuadratureType &quad,
const DofVector &dofs, HessianArray &hessians )
const
355 const int nop = quad.nop();
356 for(
int qp = 0; qp < nop; ++qp )
357 hessianAll( quad[ qp ], dofs, hessians[ qp ] );
361 template<
class Po
int,
class HessianRangeArray >
362 void hessianAll (
const Point &x, HessianRangeArray &hessians )
const
364 assert( hessians.size() >=
size() );
365 Fem::ForLoop< HessianAllRanges, 0, setIterationSize >::apply( x, hessians, storage_, offset_, basisFunctionSetTuple_ );
369 template<
class QuadratureType,
class Vector,
class DofVector >
370 void axpy (
const QuadratureType &quad,
const Vector &values, DofVector &dofs )
const
373 const unsigned int nop = quad.nop();
374 for(
unsigned int qp = 0; qp < nop; ++qp )
375 axpy( quad[ qp ], values[ qp ], dofs );
379 template<
class QuadratureType,
class VectorA,
class VectorB,
class DofVector >
380 void axpy (
const QuadratureType &quad,
const VectorA &valuesA,
const VectorB &valuesB, DofVector &dofs )
const
383 const unsigned int nop = quad.nop();
384 for(
unsigned int qp = 0; qp < nop; ++qp )
386 axpy( quad[ qp ], valuesA[ qp ], dofs );
387 axpy( quad[ qp ], valuesB[ qp ], dofs );
392 template<
class Po
int,
class DofVector >
393 void axpy (
const Point &x,
const RangeType &valueFactor, DofVector &dofs )
const
395 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, valueFactor, dofs, offset_, basisFunctionSetTuple_ );
399 template<
class Po
int,
class DofVector >
400 void axpy (
const Point &x,
const JacobianRangeType &jacobianFactor, DofVector &dofs )
const
402 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, jacobianFactor, dofs, offset_, basisFunctionSetTuple_ );
406 template<
class Po
int,
class DofVector >
407 void axpy (
const Point &x,
const HessianRangeType &hessianFactor, DofVector &dofs )
const
409 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, hessianFactor, dofs, offset_, basisFunctionSetTuple_ );
413 template<
class Po
int,
class DofVector >
414 void axpy (
const Point &x,
const RangeType &valueFactor,
const JacobianRangeType &jacobianFactor, DofVector &dofs )
const
416 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, valueFactor, jacobianFactor, dofs, offset_, basisFunctionSetTuple_ );
423 const typename SubBasisFunctionSet< i >::type& subBasisFunctionSet()
const
425 return std::get< i >( basisFunctionSetTuple_ );
429 std::size_t offset (
int i )
const
435 static int numSubBasisFunctionSets ()
442 template< std::size_t ... i >
443 int order ( std::index_sequence< i ... > )
const
445 return Std::max( std::get< i >( basisFunctionSetTuple_ ).order() ... );
449 template< std::size_t ... i >
450 std::size_t
size ( std::index_sequence< i ... > )
const
452 return Std::sum( std::get< i >( basisFunctionSetTuple_ ).
size() ... );
456 BasisFunctionSetTupleType basisFunctionSetTuple_;
467 template<
class CombineOp,
class ... BasisFunctionSets >
469 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
472 template<
class Tuple >
473 static void apply ( OffsetType &offset,
const Tuple &tuple )
475 offset[ i + 1 ] = offset[ i ] + std::get< i >( tuple ).size();
483 template<
class CombineOp,
class ... BasisFunctionSets >
485 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
489 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
491 template<
class Po
int,
class DofVector,
class Tuple >
492 static void apply (
const Point &x,
const DofVector &dofVector, RangeType &values,
const OffsetType &offset,
const Tuple &tuple )
495 std::size_t
size = std::get< i >( tuple ).size();
496 SubVector< const DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper(
size, offset[ i ] ) );
497 typename std::tuple_element< i, BasisFunctionSetTupleType >::type::RangeType thisRange;
498 std::get< i >( tuple ).evaluateAll( x, subDofVector, thisRange );
501 CombinationType::apply( rangeOffset, thisRange, values );
509 template<
class CombineOp,
class ... BasisFunctionSets >
511 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
515 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
517 template<
class Po
int,
class DofVector,
class Tuple >
518 static void apply (
const Point &x,
const DofVector &dofVector, JacobianRangeType &values,
const OffsetType &offset,
const Tuple &tuple )
521 std::size_t
size = std::get< i >( tuple ).size();
522 SubVector< const DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper(
size, offset[ i ] ) );
523 typename std::tuple_element< i, BasisFunctionSetTupleType >::type::JacobianRangeType thisJacobian;
524 std::get< i >( tuple ).jacobianAll( x, subDofVector, thisJacobian );
527 CombinationType::apply( rangeOffset, thisJacobian, values );
535 template<
class CombineOp,
class ... BasisFunctionSets >
537 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
541 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
543 template<
class Po
int,
class DofVector,
class Tuple >
544 static void apply (
const Point &x,
const DofVector &dofVector, HessianRangeType &values,
const OffsetType &offset,
const Tuple &tuple )
547 std::size_t
size = std::get< i >( tuple ).size();
548 SubVector< const DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper(
size, offset[ i ] ) );
549 typename std::tuple_element< i, BasisFunctionSetTupleType >::type::HessianRangeType thisHessian;
550 std::get< i >( tuple ).hessianAll( x, subDofVector, thisHessian );
553 CombinationType::apply( rangeOffset, thisHessian, values );
561 template<
class CombineOp,
class ... BasisFunctionSets >
563 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
566 typedef typename std::tuple_element< i, typename Storage::RangeStorageTupleType >::type ThisStorage;
567 static const int thisDimRange = std::tuple_element< i, BasisFunctionSetTupleType >::type::FunctionSpaceType::dimRange;
570 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
572 template<
class Po
int,
class RangeArray,
class Tuple >
573 static void apply (
const Point &x, RangeArray &values,
const Storage &s,
const OffsetType &offset,
const Tuple &tuple )
575 std::size_t
size = std::get< i >( tuple ).size();
576 ThisStorage &thisStorage = s.template rangeStorage< i >();
577 thisStorage.resize(
size );
579 std::get< i >( tuple ).evaluateAll( x, thisStorage );
581 for( std::size_t j = 0; j <
size; ++j )
583 values[ j + offset[ i ] ] = RangeType( 0.0 );
584 for(
int r = 0; r < thisDimRange; ++r )
585 values[ j + offset[ i ] ][ r + rangeOffset ] = thisStorage[ j ][ r ];
594 template<
class CombineOp,
class ... BasisFunctionSets >
596 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
599 typedef typename std::tuple_element< i, typename Storage::JacobianStorageTupleType >::type ThisStorage;
600 static const int thisDimRange = std::tuple_element< i, BasisFunctionSetTupleType >::type::FunctionSpaceType::dimRange;
603 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
605 template<
class Po
int,
class RangeArray,
class Tuple >
606 static void apply (
const Point &x, RangeArray &values,
const Storage &s,
const OffsetType &offset,
const Tuple &tuple )
608 std::size_t
size = std::get< i >( tuple ).size();
609 ThisStorage &thisStorage = s.template jacobianStorage< i >();
610 thisStorage.resize(
size );
612 std::get< i >( tuple ).jacobianAll( x, thisStorage );
614 for( std::size_t j = 0; j <
size; ++j )
616 values[ j + offset[ i ] ] = JacobianRangeType( RangeFieldType( 0.0 ) );
617 for(
int r = 0; r < thisDimRange; ++r )
618 values[ j + offset[ i ] ][ r + rangeOffset ] = thisStorage[ j ][ r ];
627 template<
class CombineOp,
class ... BasisFunctionSets >
629 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
632 typedef typename std::tuple_element< i, typename Storage::HessianStorageTupleType >::type ThisStorage;
633 static const int thisDimRange = std::tuple_element< i, BasisFunctionSetTupleType >::type::FunctionSpaceType::dimRange;
636 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
638 template<
class Po
int,
class RangeArray,
class Tuple >
639 static void apply (
const Point &x, RangeArray &values,
const Storage &s,
const OffsetType &offset,
const Tuple &tuple )
641 std::size_t
size = std::get< i >( tuple ).size();
642 ThisStorage &thisStorage = s.template hessianStorage< i >();
643 thisStorage.resize(
size );
645 std::get< i >( tuple ).hessianAll( x, thisStorage );
647 for( std::size_t j = 0; j <
size; ++j )
649 values[ j + offset[ i ] ] =
typename HessianRangeType::value_type( 0.0 );
650 for(
int r = 0; r < thisDimRange; ++r )
651 values[ j + offset[ i ] ][ r + rangeOffset ] = thisStorage[ j ][ r ];
660 template<
class CombineOp,
class ... BasisFunctionSets >
662 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
665 typedef typename std::tuple_element< i, BasisFunctionSetTupleType >::type::RangeType ThisRangeType;
666 typedef typename std::tuple_element< i, BasisFunctionSetTupleType >::type::JacobianRangeType ThisJacobianRangeType;
667 typedef typename std::tuple_element< i, BasisFunctionSetTupleType >::type::HessianRangeType ThisHessianRangeType;
670 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
672 typedef SubObject< const RangeType, const ThisRangeType, rangeOffset > SubRangeType;
673 typedef SubObject< const JacobianRangeType, const ThisJacobianRangeType, rangeOffset > SubJacobianRangeType;
674 typedef SubObject< const HessianRangeType, const ThisHessianRangeType, rangeOffset > SubHessianRangeType;
677 template<
class Po
int,
class DofVector,
class Tuple >
678 static void apply (
const Point &x,
const RangeType &factor, DofVector &dofVector,
const OffsetType &offset,
const Tuple &tuple )
680 std::size_t
size = std::get< i >( tuple ).size();
681 SubVector< DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper(
size, offset[ i ] ) );
682 SubRangeType subFactor( factor );
683 std::get< i >( tuple ).axpy( x, (ThisRangeType) subFactor, subDofVector );
687 template<
class Po
int,
class DofVector,
class Tuple >
688 static void apply (
const Point &x,
const JacobianRangeType &factor, DofVector &dofVector,
const OffsetType &offset,
const Tuple &tuple )
690 std::size_t
size = std::get< i >( tuple ).size();
691 SubVector< DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper(
size, offset[ i ] ) );
692 SubJacobianRangeType subFactor( factor );
693 std::get< i >( tuple ).axpy( x, (ThisJacobianRangeType) subFactor, subDofVector );
697 template<
class Po
int,
class DofVector,
class Tuple >
698 static void apply (
const Point &x,
const HessianRangeType &factor, DofVector &dofVector,
const OffsetType &offset,
const Tuple &tuple )
700 std::size_t
size = std::get< i >( tuple ).size();
701 SubVector< DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper(
size, offset[ i ] ) );
702 SubHessianRangeType subFactor( factor );
703 std::get< i >( tuple ).axpy( x, (ThisHessianRangeType) subFactor, subDofVector );
707 template<
class Po
int,
class DofVector,
class Tuple >
708 static void apply (
const Point &x,
const RangeType &rangeFactor,
const JacobianRangeType &jacobianFactor, DofVector &dofVector,
const OffsetType &offset,
const Tuple &tuple )
710 std::size_t
size = std::get< i >( tuple ).size();
711 SubVector< DofVector, OffsetSubMapper > subDofVector( dofVector, OffsetSubMapper(
size, offset[ i ] ) );
712 SubRangeType subRangeFactor( rangeFactor );
713 SubJacobianRangeType subJacobianFactor( jacobianFactor );
714 std::get< i >( tuple ).axpy( x, (ThisRangeType) subRangeFactor, (ThisJacobianRangeType) subJacobianFactor, subDofVector );
ImplementationDefined EntityType
entity type
Definition: basisfunctionsets.hh:29
Definition: explicitfieldvector.hh:75
FunctionSpaceTraits::RangeType RangeType
Type of range vector (using type of range field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:71
@ dimRange
dimension of range vector space
Definition: functionspaceinterface.hh:48
FunctionSpaceTraits::LinearMappingType JacobianRangeType
Intrinsic type used for the jacobian values has a Dune::FieldMatrix type interface.
Definition: functionspaceinterface.hh:75
FunctionSpaceTraits::DomainType DomainType
Type of domain vector (using type of domain field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:67
FunctionSpaceTraits::RangeFieldType RangeFieldType
Intrinsic type used for values in the range field (usually a double)
Definition: functionspaceinterface.hh:63
A vector valued function space.
Definition: functionspace.hh:60
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
A unique label for each type of element that can occur in a grid.