1 #ifndef DUNE_FEM_SCHEMES_GALERKIN_HH
2 #define DUNE_FEM_SCHEMES_GALERKIN_HH
11 #include <dune/common/hybridutilities.hh>
13 #include <dune/grid/common/rangegenerators.hh>
36 #include <dune/fempy/quadrature/fempyquadratures.hh>
51 template <
class T>
static std::true_type testSignature(
int (T::*)()
const);
52 template <
class T>
static std::true_type testSignature(
int (T::*)());
53 template <
class T>
static std::true_type testSignature(
unsigned int (T::*)()
const);
54 template <
class T>
static std::true_type testSignature(
unsigned int (T::*)());
55 template <
class T>
static std::true_type testSignature(std::size_t (T::*)()
const);
56 template <
class T>
static std::true_type testSignature(std::size_t (T::*)());
59 static decltype(testSignature(&T::order)) test(std::nullptr_t);
62 static std::false_type test(...);
64 using type = decltype(test<M>(
nullptr));
67 static int callOrder(
const F& f, std::false_type)
70 std::cerr <<
"WARNING: not order method available on " <<
typeid(F).name() <<
", defaulting to 1!" << std::endl;
76 static int callOrder(
const F& f, std::true_type)
83 static int order (
const F& f ) {
return callOrder(f, type() ); }
89 template<
class Integrands >
90 struct GalerkinOperator
92 typedef GalerkinOperator<Integrands> ThisType;
93 typedef std::conditional_t< Fem::IntegrandsTraits< Integrands >::isFull, Integrands, FullIntegrands< Integrands > > IntegrandsType;
95 typedef typename IntegrandsType::GridPartType GridPartType;
97 typedef typename GridPartType::ctype ctype;
98 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
100 template <
class Space>
101 struct QuadratureSelector
103 typedef CachingQuadrature< GridPartType, 0, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > InteriorQuadratureType;
104 typedef CachingQuadrature< GridPartType, 1, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > SurfaceQuadratureType;
110 typedef typename IntegrandsType::DomainValueType DomainValueType;
111 typedef typename IntegrandsType::RangeValueType RangeValueType;
112 typedef std::make_index_sequence< std::tuple_size< DomainValueType >::value > DomainValueIndices;
113 typedef std::make_index_sequence< std::tuple_size< RangeValueType >::value > RangeValueIndices;
116 template< std::size_t... i >
117 static auto makeDomainValueVector ( std::size_t maxNumLocalDofs, std::index_sequence< i... > )
119 return std::make_tuple( std::vector< std::tuple_element_t< i, DomainValueType > >( maxNumLocalDofs )... );
122 static auto makeDomainValueVector ( std::size_t maxNumLocalDofs )
124 return makeDomainValueVector( maxNumLocalDofs, DomainValueIndices() );
127 template< std::size_t... i >
128 static auto makeRangeValueVector ( std::size_t maxNumLocalDofs, std::index_sequence< i... > )
130 return std::make_tuple( std::vector< std::tuple_element_t< i, RangeValueType > >( maxNumLocalDofs )... );
133 static auto makeRangeValueVector ( std::size_t maxNumLocalDofs )
135 return makeRangeValueVector( maxNumLocalDofs, RangeValueIndices() );
138 typedef decltype( makeDomainValueVector( 0u ) ) DomainValueVectorType;
139 typedef decltype( makeRangeValueVector( 0u ) ) RangeValueVectorType;
141 static void resizeDomainValueVector ( DomainValueVectorType& vec,
const std::size_t size )
144 std::get< i >( vec ).resize( size );
148 static void resizeRangeValueVector ( RangeValueVectorType& vec,
const std::size_t size )
151 std::get< i >( vec ).resize( size );
155 template<
class LocalFunction,
class Quadrature >
156 static void evaluateQuadrature (
const LocalFunction &u,
const Quadrature &quad, std::vector< typename LocalFunction::RangeType > &phi )
158 u.evaluateQuadrature( quad, phi );
161 template<
class LocalFunction,
class Quadrature>
162 static void evaluateQuadrature (
const LocalFunction &u,
const Quadrature &quad, std::vector< typename LocalFunction::JacobianRangeType > &phi )
164 u.jacobianQuadrature( quad, phi );
167 template<
class LocalFunction,
class Quadrature >
168 static void evaluateQuadrature (
const LocalFunction &u,
const Quadrature &quad, std::vector< typename LocalFunction::HessianRangeType > &phi )
170 u.hessianQuadrature( quad, phi );
173 template<
class LocalFunction,
class Po
int >
176 u.evaluate( x, phi );
179 template<
class LocalFunction,
class Po
int >
182 u.jacobian( x, phi );
185 template<
class LocalFunction,
class Po
int >
191 template<
class LocalFunction,
class Point,
class... T >
192 static void value (
const LocalFunction &u,
const Point &x, std::tuple< T... > &phi )
194 Hybrid::forEach( std::index_sequence_for< T... >(), [ &u, &x, &phi ] (
auto i ) { GalerkinOperator::value( u, x, std::get< i >( phi ) ); } );
197 template<
class Basis,
class Po
int >
198 static void values (
const Basis &basis,
const Point &x, std::vector< typename Basis::RangeType > &phi )
200 basis.evaluateAll( x, phi );
203 template<
class Basis,
class Po
int >
204 static void values (
const Basis &basis,
const Point &x, std::vector< typename Basis::JacobianRangeType > &phi )
206 basis.jacobianAll( x, phi );
209 template<
class Basis,
class Po
int >
210 static void values (
const Basis &basis,
const Point &x, std::vector< typename Basis::HessianRangeType > &phi )
212 basis.hessianAll( x, phi );
215 template<
class Basis,
class Point,
class... T >
216 static void values (
const Basis &basis,
const Point &x, std::tuple< std::vector< T >... > &phi )
218 Hybrid::forEach( std::index_sequence_for< T... >(), [ &basis, &x, &phi ] (
auto i ) { GalerkinOperator::values( basis, x, std::get< i >( phi ) ); } );
221 template<
class LocalFunction,
class Po
int >
222 static DomainValueType domainValue (
const LocalFunction &u,
const Point &x )
229 static DomainValueType domainValue (
const unsigned int qpIdx, DomainValueVectorType& vec)
232 Hybrid::forEach( DomainValueIndices(), [ &qpIdx, &vec, &phi ] (
auto i ) {
233 std::get< i > ( phi ) = std::get< i >( vec )[ qpIdx ];
238 template<
class LocalFunction,
class Quadrature >
239 static void domainValue (
const LocalFunction &u,
const Quadrature& quadrature, DomainValueVectorType &result )
241 Hybrid::forEach( DomainValueIndices(), [ &u, &quadrature, &result ] (
auto i ) {
242 auto& vec = std::get< i >( result );
243 vec.resize( quadrature.nop() );
244 ThisType::evaluateQuadrature( u, quadrature, vec );
248 template<
class Phi, std::size_t... i >
249 static auto value (
const Phi &phi, std::size_t col, std::index_sequence< i... > )
251 return std::make_tuple( std::get< i >( phi )[ col ]... );
254 template<
class... T >
255 static auto value (
const std::tuple< std::vector< T >... > &phi, std::size_t col )
257 return value( phi, col, std::index_sequence_for< T... >() );
260 static void assignRange( RangeValueVectorType& ranges,
const std::size_t idx,
const RangeValueType& range )
262 Hybrid::forEach( RangeValueIndices(), [ &ranges, &idx, &range ] (
auto i ) {
263 std::get< i >( ranges )[ idx ] = std::get< i >( range );
267 static void assignRange( RangeValueVectorType& ranges,
const std::size_t idx,
const RangeValueType& range,
const W &weight )
269 Hybrid::forEach( RangeValueIndices(), [ &ranges, &idx, &range, &weight ] (
auto i ) {
270 std::get< i >( ranges )[ idx ] = std::get< i >( range );
271 std::get< i >( ranges )[ idx ] *= weight;
275 static void assignDomain( DomainValueVectorType& domains,
const std::size_t idx,
const DomainValueType& domain )
277 Hybrid::forEach( DomainValueIndices(), [ &domains, &idx, &domain ] (
auto i ) {
278 std::get< i >( domains )[ idx ] = std::get< i >( domain );
282 template <
class W,
class Quadrature>
283 static void axpyQuadrature( W& w,
const Quadrature& quadrature, RangeValueVectorType& ranges )
285 Hybrid::forEach( RangeValueIndices(), [ &w, &quadrature, &ranges ] (
auto i ) {
286 w.axpyQuadrature( quadrature, std::get< i >( ranges ) );
293 template<
class U,
class W >
294 void addInteriorIntegral (
const U &u, W &w )
const
296 if( !integrands_.init( u.entity() ) )
299 const auto geometry = u.entity().geometry();
301 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: InteriorQuadratureType InteriorQuadratureType;
302 InteriorQuadratureType quadrature( u.entity(), interiorQuadratureOrder(maxOrder(u, w)) );
305 DomainValueVectorType& domains = domainValues_;
306 domainValue( u, quadrature, domains );
308 auto& ranges = values_;
309 resizeRangeValueVector( ranges, quadrature.nop() );
312 for(
const auto qp : quadrature )
314 const ctype weight = qp.weight() * geometry.integrationElement( qp.position() );
315 assignRange( ranges, qp.index(), integrands_.interior( qp, domainValue( qp.index(), domains ) ), weight );
319 axpyQuadrature( w, quadrature, ranges );
322 template<
class U,
class J >
323 void addLinearizedInteriorIntegral (
const U &u, DomainValueVectorType &phi, J &j )
const
325 if( !integrands_.init( u.entity() ) )
328 const auto geometry = u.entity().geometry();
329 const auto &domainBasis = j.domainBasisFunctionSet();
330 const auto &rangeBasis = j.rangeBasisFunctionSet();
332 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: InteriorQuadratureType InteriorQuadratureType;
333 InteriorQuadratureType quadrature( u.entity(), interiorQuadratureOrder( maxOrder( u, domainBasis, rangeBasis )) );
334 const size_t domainSize = domainBasis.size();
335 const size_t quadNop = quadrature.nop();
337 auto& basisValues = basisValues_;
338 resizeDomainValueVector( basisValues_, domainSize );
341 DomainValueVectorType& domains = domainValues_;
342 domainValue( u, quadrature, domains );
344 rangeValues_.resize( domainSize );
345 for( std::size_t col = 0; col < domainSize; ++col )
347 resizeRangeValueVector( rangeValues_[ col ], quadNop );
351 for(
const auto qp : quadrature )
353 values( domainBasis, qp, basisValues );
354 const auto weight = qp.weight() * geometry.integrationElement( qp.position() );
355 auto integrand = integrands_.linearizedInterior( qp, domainValue( qp.index(), domains ) );
356 for( std::size_t col = 0; col < domainSize; ++col )
358 assignRange( rangeValues_[ col ], qp.index(), integrand( value( basisValues, col ) ), weight );
363 for( std::size_t col = 0; col < domainSize; ++col )
365 LocalMatrixColumn< J > jCol( j, col );
366 axpyQuadrature( jCol, quadrature, rangeValues_[ col ] );
372 template<
class Intersection,
class U,
class W >
373 void addBoundaryIntegral (
const Intersection &intersection,
const U &u, W &w )
const
375 if( !integrands_.init( intersection ) )
378 const auto geometry = intersection.geometry();
379 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
380 SurfaceQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( u, w )), SurfaceQuadratureType::INSIDE );
381 for(
const auto qp : quadrature )
383 const ctype weight = qp.weight() * geometry.integrationElement( qp.localPosition() );
385 RangeValueType integrand = integrands_.boundary( qp, domainValue( u, qp ) );
387 Hybrid::forEach( RangeValueIndices(), [ &qp, &w, &integrand, weight ] (
auto i ) {
388 std::get< i >( integrand ) *= weight;
389 w.axpy( qp, std::get< i >( integrand ) );
394 template<
class Intersection,
class U,
class J >
395 void addLinearizedBoundaryIntegral (
const Intersection &intersection,
const U &u, DomainValueVectorType &phi, J &j )
const
397 if( !integrands_.init( intersection ) )
400 const auto geometry = intersection.geometry();
401 const auto &domainBasis = j.domainBasisFunctionSet();
402 const auto &rangeBasis = j.rangeBasisFunctionSet();
404 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
405 SurfaceQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder(u, domainBasis, rangeBasis )), SurfaceQuadratureType::INSIDE );
406 for(
const auto qp : quadrature )
408 const ctype weight = qp.weight() * geometry.integrationElement( qp.localPosition() );
410 values( domainBasis, qp, phi );
411 auto integrand = integrands_.linearizedBoundary( qp, domainValue( u, qp ) );
413 for( std::size_t col = 0, cols = domainBasis.size(); col < cols; ++col )
415 LocalMatrixColumn< J > jCol( j, col );
416 RangeValueType intPhi = integrand( value( phi, col ) );
418 Hybrid::forEach( RangeValueIndices(), [ &qp, &jCol, &intPhi, weight ] (
auto i ) {
419 std::get< i >( intPhi ) *= weight;
420 jCol.axpy( qp, std::get< i >( intPhi ) );
429 template<
bool conforming,
class Intersection,
class U,
class W >
430 void addSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut, W &wIn )
const
432 const auto geometry = intersection.geometry();
434 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
435 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
436 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( uIn, uOut, wIn)),
false );
437 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
439 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
441 const auto qpIn = quadrature.inside()[ qp ];
442 const auto qpOut = quadrature.outside()[ qp ];
443 std::pair< RangeValueType, RangeValueType > integrand = integrands_.skeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
445 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &wIn, &integrand, weight ] (
auto i ) {
446 std::get< i >( integrand.first ) *= weight;
447 wIn.axpy( qpIn, std::get< i >( integrand.first ) );
452 template<
bool conforming,
class Intersection,
class U,
class W >
453 void addSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut, W &wIn, W &wOut )
const
455 const auto geometry = intersection.geometry();
456 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
457 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
458 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( uIn, uOut, wIn, wOut)),
false );
459 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
461 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
463 const auto qpIn = quadrature.inside()[ qp ];
464 const auto qpOut = quadrature.outside()[ qp ];
465 std::pair< RangeValueType, RangeValueType > integrand = integrands_.skeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
467 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &wIn, &qpOut, &wOut, &integrand, weight ] (
auto i ) {
468 std::get< i >( integrand.first ) *= weight;
469 wIn.axpy( qpIn, std::get< i >( integrand.first ) );
471 std::get< i >( integrand.second ) *= weight;
472 wOut.axpy( qpOut, std::get< i >( integrand.second ) );
477 template<
bool conforming,
class Intersection,
class U,
class J >
478 void addLinearizedSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut,
479 DomainValueVectorType &phiIn, DomainValueVectorType &phiOut, J &jInIn, J &jOutIn )
const
481 const auto &domainBasisIn = jInIn.domainBasisFunctionSet();
482 const auto &domainBasisOut = jOutIn.domainBasisFunctionSet();
484 const auto &rangeBasisIn = jInIn.rangeBasisFunctionSet();
486 const int order =
std::max( maxOrder(uIn, uOut), maxOrder( domainBasisIn, domainBasisOut, rangeBasisIn ));
488 const auto geometry = intersection.geometry();
489 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
490 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
491 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(order),
false );
492 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
494 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
496 const auto qpIn = quadrature.inside()[ qp ];
497 const auto qpOut = quadrature.outside()[ qp ];
499 values( domainBasisIn, qpIn, phiIn );
500 values( domainBasisOut, qpOut, phiOut );
502 auto integrand = integrands_.linearizedSkeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
503 for( std::size_t col = 0, cols = domainBasisIn.size(); col < cols; ++col )
505 LocalMatrixColumn< J > jInInCol( jInIn, col );
506 std::pair< RangeValueType, RangeValueType > intPhi = integrand.first( value( phiIn, col ) );
508 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jInInCol, &intPhi, weight ] (
auto i ) {
509 std::get< i >( intPhi.first ) *= weight;
510 jInInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
513 for( std::size_t col = 0, cols = domainBasisOut.size(); col < cols; ++col )
515 LocalMatrixColumn< J > jOutInCol( jOutIn, col );
516 std::pair< RangeValueType, RangeValueType > intPhi = integrand.second( value( phiOut, col ) );
518 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jOutInCol, &intPhi, weight ] (
auto i ) {
519 std::get< i >( intPhi.first ) *= weight;
520 jOutInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
526 template<
bool conforming,
class Intersection,
class U,
class J >
527 void addLinearizedSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut,
528 DomainValueVectorType &phiIn, DomainValueVectorType &phiOut, J &jInIn, J &jOutIn, J &jInOut, J &jOutOut )
const
530 const auto &domainBasisIn = jInIn.domainBasisFunctionSet();
531 const auto &domainBasisOut = jOutIn.domainBasisFunctionSet();
533 const auto &rangeBasisIn = jInIn.rangeBasisFunctionSet();
534 const auto &rangeBasisOut = jInOut.rangeBasisFunctionSet();
536 const int order =
std::max( maxOrder(uIn, uOut), maxOrder( domainBasisIn, domainBasisOut, rangeBasisIn, rangeBasisOut ));
538 const auto geometry = intersection.geometry();
539 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
540 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
541 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(order),
false );
542 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
544 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
546 const auto qpIn = quadrature.inside()[ qp ];
547 const auto qpOut = quadrature.outside()[ qp ];
549 values( domainBasisIn, qpIn, phiIn );
550 values( domainBasisOut, qpOut, phiOut );
552 auto integrand = integrands_.linearizedSkeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
553 for( std::size_t col = 0, cols = domainBasisIn.size(); col < cols; ++col )
555 LocalMatrixColumn< J > jInInCol( jInIn, col );
556 LocalMatrixColumn< J > jInOutCol( jInOut, col );
557 std::pair< RangeValueType, RangeValueType > intPhi = integrand.first( value( phiIn, col ) );
559 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jInInCol, &qpOut, &jInOutCol, &intPhi, weight ] (
auto i ) {
560 std::get< i >( intPhi.first ) *= weight;
561 jInInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
563 std::get< i >( intPhi.second ) *= weight;
564 jInOutCol.axpy( qpOut, std::get< i >( intPhi.second ) );
567 for( std::size_t col = 0, cols = domainBasisOut.size(); col < cols; ++col )
569 LocalMatrixColumn< J > jOutInCol( jOutIn, col );
570 LocalMatrixColumn< J > jOutOutCol( jOutOut, col );
571 std::pair< RangeValueType, RangeValueType > intPhi = integrand.second( value( phiOut, col ) );
573 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jOutInCol, &qpOut, &jOutOutCol, &intPhi, weight ] (
auto i ) {
574 std::get< i >( intPhi.first ) *= weight;
575 jOutInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
577 std::get< i >( intPhi.second ) *= weight;
578 jOutOutCol.axpy( qpOut, std::get< i >( intPhi.second ) );
585 template<
class Intersection,
class U,
class... W >
586 void addSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut, W &... w )
const
588 if( !integrands_.init( intersection ) )
591 if( intersection.conforming() )
592 addSkeletonIntegral< true >( intersection, uIn, uOut, w... );
594 addSkeletonIntegral< false >( intersection, uIn, uOut, w... );
597 template<
class Intersection,
class U,
class... J >
598 void addLinearizedSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut, DomainValueVectorType &phiIn, DomainValueVectorType &phiOut, J &... j )
const
600 if( !integrands_.init( intersection ) )
603 if( intersection.conforming() )
604 addLinearizedSkeletonIntegral< true >( intersection, uIn, uOut, phiIn, phiOut, j... );
606 addLinearizedSkeletonIntegral< false >( intersection, uIn, uOut, phiIn, phiOut, j... );
611 template<
class... Args >
612 explicit GalerkinOperator (
const GridPartType &gridPart, Args &&... args )
613 : gridPart_( gridPart ),
614 integrands_( std::forward< Args >( args )... ),
615 defaultInteriorOrder_( [] (const int order) {
return 2 * order; } ),
616 defaultSurfaceOrder_ ( [] (
const int order) {
return 2 * order + 1; } ),
617 interiorQuadOrder_(0), surfaceQuadOrder_(0),
622 void setCommunicate(
const bool communicate )
624 communicate_ = communicate;
627 std::cout <<
"Impl::GalerkinOperator::setCommunicate: communicate was disabled!" << std::endl;
631 void setQuadratureOrders(
unsigned int interior,
unsigned int surface)
633 interiorQuadOrder_ = interior;
634 surfaceQuadOrder_ = surface;
637 IntegrandsType &model()
const
643 template<
class Gr
idFunction,
class DiscreteFunction >
644 void evaluate (
const GridFunction &u, DiscreteFunction &w, std::false_type )
const
648 typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
656 for(
const EntityType &entity : elements( gridPart(), Partitions::interiorBorder ) )
658 uLocal.bind( entity );
659 wLocal.bind( entity );
662 if( integrands_.hasInterior() )
663 addInteriorIntegral( uLocal, wLocal );
665 if( integrands_.hasBoundary() && entity.hasBoundaryIntersections() )
667 for(
const auto &intersection : intersections( gridPart(), entity ) )
669 if( intersection.boundary() )
670 addBoundaryIntegral( intersection, uLocal, wLocal );
674 w.addLocalDofs( entity, wLocal.localDofVector() );
678 template<
class Gr
idFunction,
class DiscreteFunction >
679 void evaluate (
const GridFunction &u, DiscreteFunction &w, std::true_type )
const
683 typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
692 const auto &indexSet = gridPart().indexSet();
693 for(
const EntityType &inside : elements( gridPart(), Partitions::interiorBorder ) )
695 uInside.bind( inside );
696 wInside.bind( inside );
699 if( integrands_.hasInterior() )
700 addInteriorIntegral( uInside, wInside );
702 for(
const auto &intersection : intersections( gridPart(), inside ) )
706 if( intersection.neighbor() )
708 const EntityType &outside = intersection.outside();
711 uOutside.bind( outside );
712 if( outside.partitionType() != InteriorEntity )
713 addSkeletonIntegral( intersection, uInside, uOutside, wInside );
714 else if( indexSet.index( inside ) < indexSet.index( outside ) )
716 wOutside.bind( outside );
718 addSkeletonIntegral( intersection, uInside, uOutside, wInside, wOutside );
719 w.addLocalDofs( outside, wOutside.localDofVector() );
722 else if( intersection.boundary() )
724 if( integrands_.hasBoundary() )
725 addBoundaryIntegral( intersection, uInside, wInside );
729 w.addLocalDofs( inside, wInside.localDofVector() );
734 template<
class Gr
idFunction,
class DiscreteFunction >
735 void evaluate (
const GridFunction &u, DiscreteFunction &w )
const
737 static_assert( std::is_same< typename GridFunction::GridPartType, GridPartType >::value,
"Argument 'u' and Integrands must be defined on the same grid part." );
738 static_assert( std::is_same< typename DiscreteFunction::GridPartType, GridPartType >::value,
"Argument 'w' and Integrands must be defined on the same grid part." );
740 if( integrands_.hasSkeleton() )
741 evaluate( u, w, std::true_type() );
743 evaluate( u, w, std::false_type() );
753 template<
class Gr
idFunction,
class JacobianOperator >
754 void assemble (
const GridFunction &u, JacobianOperator &jOp, std::false_type )
const
756 typedef typename JacobianOperator::DomainSpaceType DomainSpaceType;
757 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
759 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
765 DiagonalStencil< DomainSpaceType, RangeSpaceType > stencil( jOp.domainSpace(), jOp.rangeSpace() );
766 jOp.reserve( stencil );
769 const std::size_t maxNumLocalDofs = jOp.domainSpace().blockMapper().maxNumDofs() * jOp.domainSpace().localBlockSize;
770 DomainValueVectorType phi = makeDomainValueVector( maxNumLocalDofs );
772 TemporaryLocalMatrixType jOpLocal( jOp.domainSpace(), jOp.rangeSpace() );
775 for(
const EntityType &entity : elements( gridPart(), Partitions::interiorBorder ) )
777 uLocal.bind( entity );
779 jOpLocal.init( entity, entity );
782 if( integrands_.hasInterior() )
783 addLinearizedInteriorIntegral( uLocal, phi, jOpLocal );
785 if( integrands_.hasBoundary() && entity.hasBoundaryIntersections() )
787 for(
const auto &intersection : intersections( gridPart(), entity ) )
789 if( intersection.boundary() )
790 addLinearizedBoundaryIntegral( intersection, uLocal, phi, jOpLocal );
794 jOp.addLocalMatrix( entity, entity, jOpLocal );
798 template<
class Gr
idFunction,
class JacobianOperator >
799 void assemble (
const GridFunction &u, JacobianOperator &jOp, std::true_type )
const
801 typedef typename JacobianOperator::DomainSpaceType DomainSpaceType;
802 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
804 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
810 DiagonalAndNeighborStencil< DomainSpaceType, RangeSpaceType > stencil( jOp.domainSpace(), jOp.rangeSpace() );
811 jOp.reserve( stencil );
814 const std::size_t maxNumLocalDofs = jOp.domainSpace().blockMapper().maxNumDofs() * jOp.domainSpace().localBlockSize;
815 DomainValueVectorType phiIn = makeDomainValueVector( maxNumLocalDofs );
816 DomainValueVectorType phiOut = makeDomainValueVector( maxNumLocalDofs );
818 TemporaryLocalMatrixType jOpInIn( jOp.domainSpace(), jOp.rangeSpace() ), jOpOutIn( jOp.domainSpace(), jOp.rangeSpace() );
819 TemporaryLocalMatrixType jOpInOut( jOp.domainSpace(), jOp.rangeSpace() ), jOpOutOut( jOp.domainSpace(), jOp.rangeSpace() );
822 const auto &indexSet = gridPart().indexSet();
823 for(
const EntityType &inside : elements( gridPart(), Partitions::interiorBorder ) )
827 jOpInIn.init( inside, inside );
830 if( integrands_.hasInterior() )
831 addLinearizedInteriorIntegral( uIn, phiIn, jOpInIn );
833 for(
const auto &intersection : intersections( gridPart(), inside ) )
837 if( intersection.neighbor() )
839 const EntityType &outside = intersection.outside();
841 jOpOutIn.init( outside, inside );
845 uOut.bind( outside );
847 if( outside.partitionType() != InteriorEntity )
848 addLinearizedSkeletonIntegral( intersection, uIn, uOut, phiIn, phiOut, jOpInIn, jOpOutIn );
849 else if( indexSet.index( inside ) < indexSet.index( outside ) )
851 jOpInOut.init( inside, outside );
853 jOpOutOut.init( outside, outside );
856 addLinearizedSkeletonIntegral( intersection, uIn, uOut, phiIn, phiOut, jOpInIn, jOpOutIn, jOpInOut, jOpOutOut );
858 jOp.addLocalMatrix( inside, outside, jOpInOut );
859 jOp.addLocalMatrix( outside, outside, jOpOutOut );
862 jOp.addLocalMatrix( outside, inside, jOpOutIn );
864 else if( intersection.boundary() )
866 if( integrands_.hasBoundary() )
867 addLinearizedBoundaryIntegral( intersection, uIn, phiIn, jOpInIn );
871 jOp.addLocalMatrix( inside, inside, jOpInIn );
876 template<
class Gr
idFunction,
class JacobianOperator >
877 void assemble (
const GridFunction &u, JacobianOperator &jOp )
const
879 static_assert( std::is_same< typename GridFunction::GridPartType, GridPartType >::value,
"Argument 'u' and Integrands must be defined on the same grid part." );
880 static_assert( std::is_same< typename JacobianOperator::DomainSpaceType::GridPartType, GridPartType >::value,
"Argument 'jOp' and Integrands must be defined on the same grid part." );
881 static_assert( std::is_same< typename JacobianOperator::RangeSpaceType::GridPartType, GridPartType >::value,
"Argument 'jOp' and Integrands must be defined on the same grid part." );
883 if( integrands_.hasSkeleton() )
884 assemble( u, jOp, std::true_type() );
886 assemble( u, jOp, std::false_type() );
895 const GridPartType &gridPart ()
const {
return gridPart_; }
897 unsigned int interiorQuadratureOrder(
unsigned int order)
const {
return interiorQuadOrder_ == 0 ? defaultInteriorOrder_(order) : interiorQuadOrder_; }
898 unsigned int surfaceQuadratureOrder(
unsigned int order)
const {
return surfaceQuadOrder_ == 0 ? defaultSurfaceOrder_ (order) : surfaceQuadOrder_; }
902 int maxOrder(
const U& u )
const
904 return CallOrder< U > :: order( u );
907 template<
class U,
class W >
908 int maxOrder(
const U& u,
const W& w )
const
910 return std::max( maxOrder( u ), maxOrder( w ) );
913 template<
class U,
class V,
class W >
914 int maxOrder(
const U& u,
const V& v,
const W& w )
const
916 return std::max( maxOrder( u, v ), maxOrder( w ) );
919 template<
class U,
class V,
class W,
class X >
920 int maxOrder(
const U& u,
const V& v,
const W& w,
const X& x )
const
922 return std::max( maxOrder( u, v ), maxOrder( w, x) );
926 const GridPartType &gridPart_;
927 mutable IntegrandsType integrands_;
929 mutable std::function<int(
const int)> defaultInteriorOrder_;
930 mutable std::function<int(
const int)> defaultSurfaceOrder_;
932 unsigned int interiorQuadOrder_;
933 unsigned int surfaceQuadOrder_;
935 mutable std::vector< RangeValueVectorType > rangeValues_;
936 mutable RangeValueVectorType values_;
937 mutable DomainValueVectorType basisValues_;
938 mutable DomainValueVectorType domainValues_;
951 template<
class Integrands,
class DomainFunction,
class RangeFunction = DomainFunction >
953 :
public virtual Operator< DomainFunction, RangeFunction >
958 static_assert( std::is_same< typename DomainFunctionType::GridPartType, typename RangeFunctionType::GridPartType >::value,
"DomainFunction and RangeFunction must be defined on the same grid part." );
962 template<
class... Args >
973 impl_.evaluate( u, w );
976 template<
class Gr
idFunction >
979 return impl_.evaluate( u, w );
989 Impl::GalerkinOperator< Integrands >
impl_;
997 template<
class Integrands,
class JacobianOperator >
999 :
public GalerkinOperator< Integrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType >,
1014 template<
class... Args >
1024 impl_.assemble( u, jOp );
1027 template<
class Gr
idFunction >
1030 impl_.assemble( u, jOp );
1053 template<
class Integrands,
class DomainFunction,
class RangeFunction >
1064 template<
class... Args >
1075 template <
class LinearOperator,
class ModelIntegrands >
1090 template<
class Gr
idFunction >
1096 template<
class Gr
idFunction >
1099 (*this).jacobian( u, jOp );
1108 template <
class O,
bool addDirichletBC>
1109 struct DirichletBlockSelector {
using type = void; };
1111 struct DirichletBlockSelector<O,true> {
using type =
typename O::DirichletBlockVector; };
1112 template<
class Integrands,
class LinearOperator,
class InverseOperator,
bool addDirichletBC,
1113 template <
class,
class>
class DifferentiableGalerkinOperatorImpl = DifferentiableGalerkinOperator >
1114 struct GalerkinSchemeImpl
1116 typedef InverseOperator InverseOperatorType;
1117 typedef Integrands ModelType;
1118 using DifferentiableOperatorType = std::conditional_t< addDirichletBC,
1120 DifferentiableGalerkinOperatorImpl< Integrands, LinearOperator > >;
1121 using DirichletBlockVector =
typename DirichletBlockSelector<
1123 DifferentiableGalerkinOperatorImpl< Integrands, LinearOperator >>,
1124 addDirichletBC>::type;
1126 typedef typename DifferentiableOperatorType::DomainFunctionType DomainFunctionType;
1127 typedef typename DifferentiableOperatorType::RangeFunctionType RangeFunctionType;
1128 typedef typename DifferentiableOperatorType::JacobianOperatorType LinearOperatorType;
1129 typedef typename DifferentiableOperatorType::JacobianOperatorType JacobianOperatorType;
1131 typedef RangeFunctionType DiscreteFunctionType;
1132 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeFunctionSpaceType;
1133 typedef typename RangeFunctionType::DiscreteFunctionSpaceType DomainFunctionSpaceType;
1134 typedef typename RangeFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
1136 typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
1137 typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
1140 typedef InverseOperator LinearInverseOperatorType;
1145 SolverInfo (
bool converged,
int linearIterations,
int nonlinearIterations )
1146 : converged( converged ), linearIterations( linearIterations ), nonlinearIterations( nonlinearIterations )
1150 int linearIterations, nonlinearIterations;
1153 GalerkinSchemeImpl (
const DiscreteFunctionSpaceType &dfSpace,
1154 const Integrands &integrands,
1156 : dfSpace_( dfSpace ),
1157 fullOperator_( dfSpace, dfSpace, std::move( integrands ) ),
1161 void setQuadratureOrders(
unsigned int interior,
unsigned int surface) { fullOperator().setQuadratureOrders(interior,surface); }
1163 const DifferentiableOperatorType &fullOperator()
const {
return fullOperator_; }
1164 DifferentiableOperatorType &fullOperator() {
return fullOperator_; }
1166 void constraint ( DiscreteFunctionType &u )
const {}
1168 template<
class Gr
idFunction >
1169 void operator() (
const GridFunction &u, DiscreteFunctionType &w )
const
1171 fullOperator()( u, w );
1174 void setErrorMeasure(ErrorMeasureType &errorMeasure)
const
1176 invOp_.setErrorMeasure(errorMeasure);
1179 SolverInfo solve (
const DiscreteFunctionType &rhs, DiscreteFunctionType &solution )
const
1181 DiscreteFunctionType rhs0 = rhs;
1182 setZeroConstraints( rhs0 );
1183 setModelConstraints( solution );
1185 invOp_.bind(fullOperator());
1186 invOp_( rhs0, solution );
1188 return SolverInfo( invOp_.converged(), invOp_.linearIterations(), invOp_.iterations() );
1191 SolverInfo solve ( DiscreteFunctionType &solution )
const
1193 DiscreteFunctionType bnd( solution );
1195 setModelConstraints( solution );
1196 invOp_.bind(fullOperator());
1197 invOp_( bnd, solution );
1199 return SolverInfo( invOp_.converged(), invOp_.linearIterations(), invOp_.iterations() );
1202 template<
class Gr
idFunction >
1203 void jacobian(
const GridFunction &ubar, LinearOperatorType &linearOp)
const
1205 fullOperator().jacobian( ubar, linearOp );
1208 const DiscreteFunctionSpaceType &space ()
const {
return dfSpace_; }
1209 const GridPartType &gridPart ()
const {
return space().gridPart(); }
1210 ModelType &model()
const {
return fullOperator().model(); }
1212 void setConstraints( DomainFunctionType &u )
const
1214 if constexpr (addDirichletBC)
1215 fullOperator().setConstraints( u );
1217 void setConstraints(
const typename DiscreteFunctionType::RangeType &value, DiscreteFunctionType &u )
const
1219 if constexpr (addDirichletBC)
1220 fullOperator().setConstraints( value, u );
1222 void setConstraints(
const DiscreteFunctionType &u, DiscreteFunctionType &v )
const
1224 if constexpr (addDirichletBC)
1225 fullOperator().setConstraints( u, v );
1227 void subConstraints(
const DiscreteFunctionType &u, DiscreteFunctionType &v )
const
1229 if constexpr (addDirichletBC)
1230 fullOperator().subConstraints( u, v );
1232 const auto& dirichletBlocks()
const
1234 if constexpr (addDirichletBC)
1235 return fullOperator().dirichletBlocks();
1239 void setZeroConstraints( DiscreteFunctionType &u )
const
1241 if constexpr (addDirichletBC)
1242 fullOperator().setConstraints(
typename DiscreteFunctionType::RangeType(0), u );
1244 void setModelConstraints( DiscreteFunctionType &u )
const
1246 if constexpr (addDirichletBC)
1247 fullOperator().setConstraints( u );
1249 const DiscreteFunctionSpaceType &dfSpace_;
1250 DifferentiableOperatorType fullOperator_;
1251 mutable NewtonOperatorType invOp_;
1259 template<
class Integrands,
class LinearOperator,
class InverseOperator,
bool addDirichletBC >
Definition: bindguard.hh:11
Impl::GalerkinSchemeImpl< Integrands, LinearOperator, InverseOperator, addDirichletBC, DifferentiableGalerkinOperator > GalerkinScheme
Definition: galerkin.hh:1261
typename Impl::ConstLocalFunction< GridFunction >::Type ConstLocalFunction
Definition: const.hh:517
BasicParameterReader< std::function< const std::string *(const std::string &, const std::string *) > > ParameterReader
Definition: reader.hh:315
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
static constexpr T max(T a)
Definition: utility.hh:77
FunctionSpaceType::RangeType RangeType
type of range vectors, i.e., type of function values
Definition: localfunction.hh:107
FunctionSpaceType::HessianRangeType HessianRangeType
type of the Hessian
Definition: localfunction.hh:111
FunctionSpaceType::JacobianRangeType JacobianRangeType
type of the Jacobian, i.e., type of evaluated Jacobian matrix
Definition: localfunction.hh:109
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
static ParameterContainer & container()
Definition: io/parameter.hh:193
operator providing a Jacobian through automatic differentiation
Definition: automaticdifferenceoperator.hh:86
abstract differentiable operator
Definition: differentiableoperator.hh:29
BaseType::RangeFunctionType RangeFunctionType
type of discrete function in the operator's range
Definition: differentiableoperator.hh:40
JacobianOperator JacobianOperatorType
type of linear operator modelling the operator's Jacobian
Definition: differentiableoperator.hh:35
BaseType::DomainFunctionType DomainFunctionType
type of discrete function in the operator's domain
Definition: differentiableoperator.hh:38
abstract operator
Definition: operator.hh:34
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:36
abstract affine-linear operator
Definition: operator.hh:87
Definition: dirichletwrapper.hh:27
Definition: galerkin.hh:954
GalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: galerkin.hh:963
ModelType & model() const
Definition: galerkin.hh:986
Impl::GalerkinOperator< Integrands > impl_
Definition: galerkin.hh:989
Integrands DirichletModelType
Definition: galerkin.hh:985
DomainFunction DomainFunctionType
Definition: galerkin.hh:955
RangeFunction RangeFunctionType
Definition: galerkin.hh:956
virtual void operator()(const DomainFunctionType &u, RangeFunctionType &w) const final override
Definition: galerkin.hh:971
const GridPartType & gridPart() const
Definition: galerkin.hh:982
Integrands ModelType
Definition: galerkin.hh:984
RangeFunctionType::GridPartType GridPartType
Definition: galerkin.hh:958
void setQuadratureOrders(unsigned int interior, unsigned int surface)
Definition: galerkin.hh:969
void setCommunicate(const bool communicate)
Definition: galerkin.hh:967
Definition: galerkin.hh:1001
BaseType::DomainFunctionType DomainFunctionType
Definition: galerkin.hh:1007
BaseType::GridPartType GridPartType
Definition: galerkin.hh:1012
const RangeDiscreteFunctionSpaceType & rangeSpace() const
Definition: galerkin.hh:1037
void jacobian(const GridFunction &u, JacobianOperatorType &jOp) const
Definition: galerkin.hh:1028
const DomainDiscreteFunctionSpaceType & domainSpace() const
Definition: galerkin.hh:1033
RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType
Definition: galerkin.hh:1010
virtual void jacobian(const DomainFunctionType &u, JacobianOperatorType &jOp) const final override
obtain linearization
Definition: galerkin.hh:1022
BaseType::RangeFunctionType RangeFunctionType
Definition: galerkin.hh:1008
DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: galerkin.hh:1009
const RangeDiscreteFunctionSpaceType & rSpace_
Definition: galerkin.hh:1045
JacobianOperator JacobianOperatorType
Definition: galerkin.hh:1005
const DomainDiscreteFunctionSpaceType & dSpace_
Definition: galerkin.hh:1044
DifferentiableGalerkinOperator(const DomainDiscreteFunctionSpaceType &dSpace, const RangeDiscreteFunctionSpaceType &rSpace, Args &&... args)
Definition: galerkin.hh:1015
Definition: galerkin.hh:1057
BaseType::GridPartType GridPartType
Definition: galerkin.hh:1062
AutomaticDifferenceGalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: galerkin.hh:1065
Definition: galerkin.hh:1078
ModelDifferentiableGalerkinOperator(ModelType &model, const DiscreteFunctionSpaceType &dfSpace)
Definition: galerkin.hh:1086
ModelIntegrands::ModelType ModelType
Definition: galerkin.hh:1081
LinearOperator::DomainFunctionType RangeFunctionType
Definition: galerkin.hh:1083
void apply(const GridFunction &u, RangeFunctionType &w) const
Definition: galerkin.hh:1091
DifferentiableGalerkinOperator< ModelIntegrands, LinearOperator > BaseType
Definition: galerkin.hh:1079
LinearOperator::RangeSpaceType DiscreteFunctionSpaceType
Definition: galerkin.hh:1084
void apply(const GridFunction &u, LinearOperator &jOp) const
Definition: galerkin.hh:1097
inverse operator based on a newton scheme
Definition: newtoninverseoperator.hh:209
std::function< bool(const RangeFunctionType &w, const RangeFunctionType &dw, double residualNorm) > ErrorMeasureType
Definition: newtoninverseoperator.hh:230
static int volumeOrder(const int k)
return quadrature order for volume quadratures for given polynomial order k
Definition: space/common/capabilities.hh:138
static int surfaceOrder(const int k)
return quadrature order for surface quadratures (i.e. over intersections) for given polynomial order ...
Definition: space/common/capabilities.hh:140