DUNE-FEM (unstable)

tuple.hh
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
20namespace 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 >
57 static constexpr int size () { return std::tuple_element< j, std::tuple< std::integral_constant< int, i > ... > >::type::value; }
58
59 template< int ... j >
60 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 >
66 static constexpr int offset ()
67 {
68 return mysum( sizes( std::make_integer_sequence< int, j >() ) );
69 }
70
71 private:
72 template< int ... j >
73 static constexpr int mysum ( std::integer_sequence< int, j ... > )
74 {
75 return Std::sum( j ... );
76 }
77
78 static constexpr int mysum ( std::integer_sequence< int > )
79 {
80 return 0;
81 }
82 };
83
85 typedef std::tuple< BasisFunctionSets ... > BasisFunctionSetTupleType;
86
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;
91
93 static const int dimDomain = ContainedFunctionSpaceType::dimDomain;
94
96 static const int setSize = sizeof ... ( BasisFunctionSets );
97 static const int setIterationSize = sizeof ... ( BasisFunctionSets )-1;
98
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
109 Storage () {}
110
111 template< int i >
112 typename std::tuple_element< i, RangeStorageTupleType >::type
113 & rangeStorage() const { return std::get< i >( rangeStorage_ ); }
114
115 template< int i >
116 typename std::tuple_element< i, JacobianStorageTupleType >::type
117 & jacobianStorage() const { return std::get< i >( jacobianStorage_ ); }
118
119 template< int i >
120 typename std::tuple_element< i, HessianStorageTupleType >::type
121 & 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:
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 {
142 static constexpr const int dimRange = Std::sum( static_cast< int >( BasisFunctionSets::FunctionSpaceType::dimRange ) ... );
143
145 typedef FunctionSpace< typename ContainedFunctionSpaceType::DomainFieldType,
146 typename ContainedFunctionSpaceType::RangeFieldType,
147 dimDomain, dimRange > FunctionSpaceType;
148
149 template <class ThisRange, class Range>
150 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 >
157 static constexpr int rangeOffset ()
158 {
159 return RangeIndices::template offset< i >();
160 }
161 };
162
163 template <int dummy>
164 struct CombinationTraits< dummy, TupleSpaceSummation >
165 {
167 typedef ContainedFunctionSpaceType FunctionSpaceType;
168
169 template <class ThisRange, class Range>
170 static void apply(const int rangeOffset, const ThisRange& thisRange, Range& values)
171 {
172 // sum results
173 values += thisRange;
174 }
175
176 template< int j >
177 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
195 typedef typename CombinationType :: FunctionSpaceType FunctionSpaceType;
196
198 static const int dimRange = FunctionSpaceType::dimRange;
199
201 typedef typename FunctionSpaceType::DomainType DomainType;
202
204 typedef typename FunctionSpaceType::RangeType RangeType;
205
207 typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
208
210 typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
211
213 typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
214
215 static_assert( Std::are_all_same< typename BasisFunctionSets::EntityType ... >::value,
216 "TupleBasisFunctionSet needs common EntityType" );
218 typedef typename std::tuple_element< 0, BasisFunctionSetTupleType >::type::EntityType EntityType;
219
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;
224
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;
229
230 // default constructor, needed by localmatrix
231 TupleBasisFunctionSet () : offset_( {{ 0 }} ) {}
232
233 // constructor taking a pack of basisFunctionSets
234 TupleBasisFunctionSet ( const BasisFunctionSets & ... basisFunctionSets )
235 : basisFunctionSetTuple_( basisFunctionSets ... ),
236 offset_()
237 {
238 offset_[ 0 ] = 0;
239 Fem::ForLoop< ComputeOffset, 0, setIterationSize >::apply( offset_, basisFunctionSetTuple_ );
240 }
241
242 // constructor taking a tuple of basisfunction sets
243 TupleBasisFunctionSet ( const BasisFunctionSetTupleType &basisFunctionSetTuple )
244 : basisFunctionSetTuple_( basisFunctionSetTuple ),
245 offset_()
246 {
247 offset_[ 0 ] = 0;
248 Fem::ForLoop< ComputeOffset, 0, setIterationSize >::apply( offset_, basisFunctionSetTuple_ );
249 }
250
251 // Basis Function Set Interface Methods
252 // ------------------------------------
253
255 int order () const
256 {
257 return order( std::index_sequence_for< BasisFunctionSets ... >() );
258 }
259
261 std::size_t size () const
262 {
263 return size( std::index_sequence_for< BasisFunctionSets ... >() );
264 }
265
267 Dune::GeometryType type () const
268 {
269 return std::get< 0 >( basisFunctionSetTuple_ ).type();
270 }
271
273 bool valid () const
274 {
275 return std::get< 0 >( basisFunctionSetTuple_ ).valid();
276 }
277
279 const EntityType &entity () const
280 {
281 return std::get< 0 >( basisFunctionSetTuple_ ).entity();
282 }
283
285 const Geometry &geometry () const
286 {
287 return std::get< 0 >( basisFunctionSetTuple_ ).geometry();
288 }
289
291 const ReferenceElementType &referenceElement () const
292 {
293 return std::get< 0 >( basisFunctionSetTuple_ ).referenceElement();
294 }
295
297 template< class Point, class DofVector >
298 void evaluateAll ( const Point &x, const DofVector &dofs, RangeType &value ) const
299 {
300 Fem::ForLoop< EvaluateAll, 0, setIterationSize >::apply( x, dofs, value, offset_, basisFunctionSetTuple_ );
301 }
302
304 template< class Point, class RangeArray >
305 void evaluateAll ( const Point &x, RangeArray &values ) const
306 {
307 assert( values.size() >= size() );
308 Fem::ForLoop< EvaluateAllRanges, 0, setIterationSize >::apply( x, values, storage_, offset_, basisFunctionSetTuple_ );
309 }
310
312 template< class QuadratureType, class DofVector, class RangeArray >
313 void evaluateAll ( const QuadratureType &quad, const DofVector &dofs, RangeArray &ranges ) const
314 {
315 const int nop = quad.nop();
316 for( int qp = 0; qp < nop; ++qp )
317 evaluateAll( quad[ qp ], dofs, ranges[ qp ] );
318 }
319
321 template< class Point, class DofVector >
322 void jacobianAll ( const Point &x, const DofVector &dofs, JacobianRangeType &jacobian ) const
323 {
324 Fem::ForLoop< JacobianAll, 0, setIterationSize >::apply( x, dofs, jacobian, offset_, basisFunctionSetTuple_ );
325 }
326
328 template< class Point, class JacobianRangeArray >
329 void jacobianAll ( const Point &x, JacobianRangeArray &jacobians ) const
330 {
331 assert( jacobians.size() >= size() );
332 Fem::ForLoop< JacobianAllRanges, 0, setIterationSize >::apply( x, jacobians, storage_, offset_, basisFunctionSetTuple_ );
333 }
334
336 template< class QuadratureType, class DofVector, class JacobianArray >
337 void jacobianAll ( const QuadratureType &quad, const DofVector &dofs, JacobianArray &jacobians ) const
338 {
339 const int nop = quad.nop();
340 for( int qp = 0; qp < nop; ++qp )
341 jacobianAll( quad[ qp ], dofs, jacobians[ qp ] );
342 }
343
345 template< class Point, class DofVector >
346 void hessianAll ( const Point &x, const DofVector &dofs, HessianRangeType &hessian ) const
347 {
348 Fem::ForLoop< HessianAll, 0, setIterationSize >::apply( x, dofs, hessian, offset_, basisFunctionSetTuple_ );
349 }
350
352 template< class QuadratureType, class DofVector, class HessianArray >
353 void hessianAll ( const QuadratureType &quad, const DofVector &dofs, HessianArray &hessians ) const
354 {
355 const int nop = quad.nop();
356 for( int qp = 0; qp < nop; ++qp )
357 hessianAll( quad[ qp ], dofs, hessians[ qp ] );
358 }
359
361 template< class Point, class HessianRangeArray >
362 void hessianAll ( const Point &x, HessianRangeArray &hessians ) const
363 {
364 assert( hessians.size() >= size() );
365 Fem::ForLoop< HessianAllRanges, 0, setIterationSize >::apply( x, hessians, storage_, offset_, basisFunctionSetTuple_ );
366 }
367
369 template< class QuadratureType, class Vector, class DofVector >
370 void axpy ( const QuadratureType &quad, const Vector &values, DofVector &dofs ) const
371 {
372 // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
373 const unsigned int nop = quad.nop();
374 for( unsigned int qp = 0; qp < nop; ++qp )
375 axpy( quad[ qp ], values[ qp ], dofs );
376 }
377
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
381 {
382 // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
383 const unsigned int nop = quad.nop();
384 for( unsigned int qp = 0; qp < nop; ++qp )
385 {
386 axpy( quad[ qp ], valuesA[ qp ], dofs );
387 axpy( quad[ qp ], valuesB[ qp ], dofs );
388 }
389 }
390
392 template< class Point, class DofVector >
393 void axpy ( const Point &x, const RangeType &valueFactor, DofVector &dofs ) const
394 {
395 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, valueFactor, dofs, offset_, basisFunctionSetTuple_ );
396 }
397
399 template< class Point, class DofVector >
400 void axpy ( const Point &x, const JacobianRangeType &jacobianFactor, DofVector &dofs ) const
401 {
402 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, jacobianFactor, dofs, offset_, basisFunctionSetTuple_ );
403 }
404
406 template< class Point, class DofVector >
407 void axpy ( const Point &x, const HessianRangeType &hessianFactor, DofVector &dofs ) const
408 {
409 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, hessianFactor, dofs, offset_, basisFunctionSetTuple_ );
410 }
411
413 template< class Point, class DofVector >
414 void axpy ( const Point &x, const RangeType &valueFactor, const JacobianRangeType &jacobianFactor, DofVector &dofs ) const
415 {
416 Fem::ForLoop< Axpy, 0, setIterationSize >::apply( x, valueFactor, jacobianFactor, dofs, offset_, basisFunctionSetTuple_ );
417 }
418
419 /***** NON Interface methods ****/
420
422 template< int i >
423 const typename SubBasisFunctionSet< i >::type& subBasisFunctionSet() const
424 {
425 return std::get< i >( basisFunctionSetTuple_ );
426 }
427
429 std::size_t offset ( int i ) const
430 {
431 return offset_[ i ];
432 }
433
435 static int numSubBasisFunctionSets ()
436 {
437 return setSize;
438 }
439
440 protected:
441 // unroll index sequence and take maximal order
442 template< std::size_t ... i >
443 int order ( std::index_sequence< i ... > ) const
444 {
445 return Std::max( std::get< i >( basisFunctionSetTuple_ ).order() ... );
446 }
447
448 // unroll index sequence and sum up sizes
449 template< std::size_t ... i >
450 std::size_t size ( std::index_sequence< i ... > ) const
451 {
452 return Std::sum( std::get< i >( basisFunctionSetTuple_ ).size() ... );
453 }
454
455 private:
456 BasisFunctionSetTupleType basisFunctionSetTuple_;
457 OffsetType offset_;
458
459 Storage storage_;
460 };
461
462
463
464 // ComputeOffset
465 // -------------
466
467 template< class CombineOp, class ... BasisFunctionSets >
468 template< int i >
469 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
470 ComputeOffset
471 {
472 template< class Tuple >
473 static void apply ( OffsetType &offset, const Tuple &tuple )
474 {
475 offset[ i + 1 ] = offset[ i ] + std::get< i >( tuple ).size();
476 }
477 };
478
479
480 // EvaluateAll
481 // -----------
482
483 template< class CombineOp, class ... BasisFunctionSets >
484 template< int i >
485 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
486 EvaluateAll
487 {
488 // only needed for product spaces
489 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
490
491 template< class Point, class DofVector, class Tuple >
492 static void apply ( const Point &x, const DofVector &dofVector, RangeType &values, const OffsetType &offset, const Tuple &tuple )
493 {
494 // evaluateAll for this BasisFunctionSet, with View on DofVector
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 );
499
500 // distribute result to values
501 CombinationType::apply( rangeOffset, thisRange, values );
502 }
503 };
504
505
506 // JacobianAll
507 // -----------
508
509 template< class CombineOp, class ... BasisFunctionSets >
510 template< int i >
511 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
512 JacobianAll
513 {
514 // only needed for product spaces
515 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
516
517 template< class Point, class DofVector, class Tuple >
518 static void apply ( const Point &x, const DofVector &dofVector, JacobianRangeType &values, const OffsetType &offset, const Tuple &tuple )
519 {
520 // jacobianAll for this BasisFunctionSet, with View on DofVector
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 );
525
526 // distribute result to values
527 CombinationType::apply( rangeOffset, thisJacobian, values );
528 }
529 };
530
531
532 // HessianAll
533 // ----------
534
535 template< class CombineOp, class ... BasisFunctionSets >
536 template< int i >
537 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
538 HessianAll
539 {
540 // only needed for product spaces
541 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
542
543 template< class Point, class DofVector, class Tuple >
544 static void apply ( const Point &x, const DofVector &dofVector, HessianRangeType &values, const OffsetType &offset, const Tuple &tuple )
545 {
546 // hessianAll for this BasisFunctionSet, with View on DofVector
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 );
551
552 // distribute result to values
553 CombinationType::apply( rangeOffset, thisHessian, values );
554 }
555 };
556
557
558 // EvaluateAllRanges
559 // -----------------
560
561 template< class CombineOp, class ... BasisFunctionSets >
562 template< int i >
563 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
564 EvaluateAllRanges
565 {
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;
568
569 // only needed for product spaces
570 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
571
572 template< class Point, class RangeArray, class Tuple >
573 static void apply ( const Point &x, RangeArray &values, const Storage &s, const OffsetType &offset, const Tuple &tuple )
574 {
575 std::size_t size = std::get< i >( tuple ).size();
576 ThisStorage &thisStorage = s.template rangeStorage< i >();
577 thisStorage.resize( size );
578
579 std::get< i >( tuple ).evaluateAll( x, thisStorage );
580
581 for( std::size_t j = 0; j < size; ++j )
582 {
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 ];
586 }
587 }
588 };
589
590
591 // JacobianAllRanges
592 // -----------------
593
594 template< class CombineOp, class ... BasisFunctionSets >
595 template< int i >
596 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
597 JacobianAllRanges
598 {
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;
601
602 // only needed for product spaces
603 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
604
605 template< class Point, class RangeArray, class Tuple >
606 static void apply ( const Point &x, RangeArray &values, const Storage &s, const OffsetType &offset, const Tuple &tuple )
607 {
608 std::size_t size = std::get< i >( tuple ).size();
609 ThisStorage &thisStorage = s.template jacobianStorage< i >();
610 thisStorage.resize( size );
611
612 std::get< i >( tuple ).jacobianAll( x, thisStorage );
613
614 for( std::size_t j = 0; j < size; ++j )
615 {
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 ];
619 }
620 }
621 };
622
623
624 // HessianAllRanges
625 // ----------------
626
627 template< class CombineOp, class ... BasisFunctionSets >
628 template< int i >
629 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
630 HessianAllRanges
631 {
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;
634
635 // only needed for product spaces
636 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
637
638 template< class Point, class RangeArray, class Tuple >
639 static void apply ( const Point &x, RangeArray &values, const Storage &s, const OffsetType &offset, const Tuple &tuple )
640 {
641 std::size_t size = std::get< i >( tuple ).size();
642 ThisStorage &thisStorage = s.template hessianStorage< i >();
643 thisStorage.resize( size );
644
645 std::get< i >( tuple ).hessianAll( x, thisStorage );
646
647 for( std::size_t j = 0; j < size; ++j )
648 {
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 ];
652 }
653 }
654 };
655
656
657 // Axpy
658 // ----
659
660 template< class CombineOp, class ... BasisFunctionSets >
661 template< int i >
662 struct TupleBasisFunctionSet< CombineOp, BasisFunctionSets ... >::
663 Axpy
664 {
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;
668
669 // only needed for product spaces
670 static const int rangeOffset = CombinationType :: template rangeOffset< i >();
671
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;
675
676 // axpy with range type
677 template< class Point, class DofVector, class Tuple >
678 static void apply ( const Point &x, const RangeType &factor, DofVector &dofVector, const OffsetType &offset, const Tuple &tuple )
679 {
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 );
684 }
685
686 // axpy with jacobian range type
687 template< class Point, class DofVector, class Tuple >
688 static void apply ( const Point &x, const JacobianRangeType &factor, DofVector &dofVector, const OffsetType &offset, const Tuple &tuple )
689 {
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 );
694 }
695
696 // axpy with hessian range type
697 template< class Point, class DofVector, class Tuple >
698 static void apply ( const Point &x, const HessianRangeType &factor, DofVector &dofVector, const OffsetType &offset, const Tuple &tuple )
699 {
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 );
704 }
705
706 // axpy with range and jacobian range type
707 template< class Point, 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 )
709 {
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 );
715 }
716 };
717
718 } // namespace Fem
719
720} // namespace Dune
721
722#endif // #ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_TUPLE_HH
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.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)