1#ifndef DUNE_FEM_SCHEMES_GALERKIN_HH
2#define DUNE_FEM_SCHEMES_GALERKIN_HH
12#include <dune/common/hybridutilities.hh>
13#include <dune/common/timer.hh>
15#include <dune/grid/common/rangegenerators.hh>
41#include <dune/fempy/quadrature/fempyquadratures.hh>
56 template <
class T>
static std::true_type testSignature(
int (T::*)()
const);
57 template <
class T>
static std::true_type testSignature(
int (T::*)());
58 template <
class T>
static std::true_type testSignature(
unsigned int (T::*)()
const);
59 template <
class T>
static std::true_type testSignature(
unsigned int (T::*)());
60 template <
class T>
static std::true_type testSignature(std::size_t (T::*)()
const);
61 template <
class T>
static std::true_type testSignature(std::size_t (T::*)());
64 static decltype(testSignature(&T::order)) test(std::nullptr_t);
67 static std::false_type test(...);
69 using type =
decltype(test<M>(
nullptr));
72 static int callOrder(
const F& f, std::false_type)
75 std::cerr <<
"WARNING: not order method available on " <<
typeid(F).name() <<
", defaulting to 1!" << std::endl;
81 static int callOrder(
const F& f, std::true_type)
88 static int order (
const F& f ) {
return callOrder(f, type() ); }
94 template<
class Integrands >
95 struct GalerkinOperator
97 typedef GalerkinOperator<Integrands> ThisType;
98 typedef std::conditional_t< Fem::IntegrandsTraits< Integrands >::isFull, Integrands, FullIntegrands< Integrands > > IntegrandsType;
100 typedef typename IntegrandsType::GridPartType GridPartType;
102 typedef typename GridPartType::ctype ctype;
103 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
107 template <
class Space>
108 struct QuadratureSelector
110 typedef CachingQuadrature< GridPartType, 0, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > InteriorQuadratureType;
111 typedef CachingQuadrature< GridPartType, 1, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > SurfaceQuadratureType;
117 template<
class... Args >
118 explicit GalerkinOperator (
const GridPartType &gridPart, Args &&... args )
119 : gridPart_( gridPart ),
120 integrands_(
std::forward< Args >( args )... ),
121 defaultInteriorOrder_( [] (const int order) {
return 2 * order; } ),
122 defaultSurfaceOrder_ ( [] (
const int order) {
return 2 * order + 1; } ),
123 interiorQuadOrder_(0), surfaceQuadOrder_(0),
124 gridSizeInterior_( 0 )
129 typedef typename IntegrandsType::DomainValueType DomainValueType;
130 typedef typename IntegrandsType::RangeValueType RangeValueType;
131 typedef std::make_index_sequence< std::tuple_size< DomainValueType >::value > DomainValueIndices;
132 typedef std::make_index_sequence< std::tuple_size< RangeValueType >::value > RangeValueIndices;
135 template< std::size_t... i >
136 static auto makeDomainValueVector ( std::size_t maxNumLocalDofs, std::index_sequence< i... > )
138 return std::make_tuple( std::vector< std::tuple_element_t< i, DomainValueType > >( maxNumLocalDofs )... );
141 static auto makeDomainValueVector ( std::size_t maxNumLocalDofs )
143 return makeDomainValueVector( maxNumLocalDofs, DomainValueIndices() );
146 template< std::size_t... i >
147 static auto makeRangeValueVector ( std::size_t maxNumLocalDofs, std::index_sequence< i... > )
149 return std::make_tuple( std::vector< std::tuple_element_t< i, RangeValueType > >( maxNumLocalDofs )... );
152 static auto makeRangeValueVector ( std::size_t maxNumLocalDofs )
154 return makeRangeValueVector( maxNumLocalDofs, RangeValueIndices() );
157 typedef decltype( makeDomainValueVector( 0u ) ) DomainValueVectorType;
158 typedef decltype( makeRangeValueVector( 0u ) ) RangeValueVectorType;
160 static void resizeDomainValueVector ( DomainValueVectorType& vec,
const std::size_t size )
163 std::get< i >( vec ).resize( size );
167 static void resizeRangeValueVector ( RangeValueVectorType& vec,
const std::size_t size )
170 std::get< i >( vec ).resize( size );
174 template<
class LocalFunction,
class Quadrature >
175 static void evaluateQuadrature (
const LocalFunction &u,
const Quadrature &quad, std::vector< typename LocalFunction::RangeType > &phi )
177 u.evaluateQuadrature( quad, phi );
180 template<
class LocalFunction,
class Quadrature>
181 static void evaluateQuadrature (
const LocalFunction &u,
const Quadrature &quad, std::vector< typename LocalFunction::JacobianRangeType > &phi )
183 u.jacobianQuadrature( quad, phi );
186 template<
class LocalFunction,
class Quadrature >
187 static void evaluateQuadrature (
const LocalFunction &u,
const Quadrature &quad, std::vector< typename LocalFunction::HessianRangeType > &phi )
189 u.hessianQuadrature( quad, phi );
192 template<
class LocalFunction,
class Po
int >
195 u.evaluate( x, phi );
198 template<
class LocalFunction,
class Po
int >
201 u.jacobian( x, phi );
204 template<
class LocalFunction,
class Po
int >
210 template<
class LocalFunction,
class Point,
class... T >
211 static void value (
const LocalFunction &u,
const Point &x, std::tuple< T... > &phi )
213 Hybrid::forEach( std::index_sequence_for< T... >(), [ &u, &x, &phi ] (
auto i ) { GalerkinOperator::value( u, x, std::get< i >( phi ) ); } );
216 template<
class Basis,
class Po
int >
217 static void values (
const Basis &basis,
const Point &x, std::vector< typename Basis::RangeType > &phi )
219 basis.evaluateAll( x, phi );
222 template<
class Basis,
class Po
int >
223 static void values (
const Basis &basis,
const Point &x, std::vector< typename Basis::JacobianRangeType > &phi )
225 basis.jacobianAll( x, phi );
228 template<
class Basis,
class Po
int >
229 static void values (
const Basis &basis,
const Point &x, std::vector< typename Basis::HessianRangeType > &phi )
231 basis.hessianAll( x, phi );
234 template<
class Basis,
class Point,
class... T >
235 static void values (
const Basis &basis,
const Point &x, std::tuple< std::vector< T >... > &phi )
237 Hybrid::forEach( std::index_sequence_for< T... >(), [ &basis, &x, &phi ] (
auto i ) { GalerkinOperator::values( basis, x, std::get< i >( phi ) ); } );
240 template<
class LocalFunction,
class Po
int >
241 static DomainValueType domainValue (
const LocalFunction &u,
const Point &x )
248 static DomainValueType domainValue (
const unsigned int qpIdx, DomainValueVectorType& vec)
251 Hybrid::forEach( DomainValueIndices(), [ &qpIdx, &vec, &phi ] (
auto i ) {
252 std::get< i > ( phi ) = std::get< i >( vec )[ qpIdx ];
257 template<
class LocalFunction,
class Quadrature >
258 static void domainValue (
const LocalFunction &u,
const Quadrature& quadrature, DomainValueVectorType &result )
260 Hybrid::forEach( DomainValueIndices(), [ &u, &quadrature, &result ] (
auto i ) {
261 auto& vec = std::get< i >( result );
262 vec.resize( quadrature.nop() );
263 ThisType::evaluateQuadrature( u, quadrature, vec );
267 template<
class Phi, std::size_t... i >
268 static auto value (
const Phi &phi, std::size_t col, std::index_sequence< i... > )
270 return std::make_tuple( std::get< i >( phi )[ col ]... );
273 template<
class... T >
274 static auto value (
const std::tuple< std::vector< T >... > &phi, std::size_t col )
276 return value( phi, col, std::index_sequence_for< T... >() );
279 static void assignRange( RangeValueVectorType& ranges,
const std::size_t idx,
const RangeValueType& range )
281 Hybrid::forEach( RangeValueIndices(), [ &ranges, &idx, &range ] (
auto i ) {
282 std::get< i >( ranges )[ idx ] = std::get< i >( range );
286 static void assignRange( RangeValueVectorType& ranges,
const std::size_t idx,
const RangeValueType& range,
const W &weight )
288 Hybrid::forEach( RangeValueIndices(), [ &ranges, &idx, &range, &weight ] (
auto i ) {
289 std::get< i >( ranges )[ idx ] = std::get< i >( range );
290 std::get< i >( ranges )[ idx ] *= weight;
294 static void assignDomain( DomainValueVectorType& domains,
const std::size_t idx,
const DomainValueType& domain )
296 Hybrid::forEach( DomainValueIndices(), [ &domains, &idx, &domain ] (
auto i ) {
297 std::get< i >( domains )[ idx ] = std::get< i >( domain );
301 template <
class W,
class Quadrature>
302 static void axpyQuadrature( W& w,
const Quadrature& quadrature, RangeValueVectorType& ranges )
304 Hybrid::forEach( RangeValueIndices(), [ &w, &quadrature, &ranges ] (
auto i ) {
305 w.axpyQuadrature( quadrature, std::get< i >( ranges ) );
312 template<
class U,
class W >
313 void addInteriorIntegral (
const U &u, W &w )
const
315 if( !integrands().init( u.entity() ) )
318 const auto geometry = u.entity().geometry();
320 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: InteriorQuadratureType InteriorQuadratureType;
321 const InteriorQuadratureType quadrature( u.entity(), interiorQuadratureOrder(maxOrder(u, w)) );
324 DomainValueVectorType& domains = domainValues_;
325 domainValue( u, quadrature, domains );
327 auto& ranges = values_;
328 resizeRangeValueVector( ranges, quadrature.nop() );
331 for(
const auto qp : quadrature )
333 const ctype weight = qp.weight() * geometry.integrationElement( qp.position() );
334 assignRange( ranges, qp.index(), integrands().interior( qp, domainValue( qp.index(), domains ) ), weight );
338 axpyQuadrature( w, quadrature, ranges );
339 integrands().unbind();
342 template<
class U,
class J >
343 void addLinearizedInteriorIntegral (
const U &u, DomainValueVectorType &phi, J &j )
const
345 if( !integrands().init( u.entity() ) )
348 const auto geometry = u.entity().geometry();
349 const auto &domainBasis = j.domainBasisFunctionSet();
350 const auto &rangeBasis = j.rangeBasisFunctionSet();
352 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: InteriorQuadratureType InteriorQuadratureType;
353 const InteriorQuadratureType quadrature( u.entity(), interiorQuadratureOrder( maxOrder( u, domainBasis, rangeBasis )) );
354 const size_t domainSize = domainBasis.size();
355 const size_t quadNop = quadrature.nop();
357 auto& basisValues = basisValues_;
358 resizeDomainValueVector( basisValues, domainSize );
361 auto& rangeValues = rangeValues_;
362 DomainValueVectorType& domains = domainValues_;
363 domainValue( u, quadrature, domains );
365 rangeValues.resize( domainSize );
366 for( std::size_t col = 0; col < domainSize; ++col )
368 resizeRangeValueVector( rangeValues[ col ], quadNop );
372 for(
const auto qp : quadrature )
374 values( domainBasis, qp, basisValues );
375 const auto weight = qp.weight() * geometry.integrationElement( qp.position() );
376 auto integrand = integrands().linearizedInterior( qp, domainValue( qp.index(), domains ) );
377 for( std::size_t col = 0; col < domainSize; ++col )
379 assignRange( rangeValues[ col ], qp.index(), integrand( value( basisValues, col ) ), weight );
384 for( std::size_t col = 0; col < domainSize; ++col )
386 LocalMatrixColumn< J > jCol( j, col );
387 axpyQuadrature( jCol, quadrature, rangeValues[ col ] );
389 integrands().unbind();
394 template<
class Intersection,
class U,
class W >
395 void addBoundaryIntegral (
const Intersection &intersection,
const U &u, W &w )
const
397 if( !integrands().init( intersection ) )
400 const auto geometry = intersection.geometry();
401 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
402 const SurfaceQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( u, w )), SurfaceQuadratureType::INSIDE );
403 for(
const auto qp : quadrature )
405 const ctype weight = qp.weight() * geometry.integrationElement( qp.localPosition() );
407 RangeValueType integrand = integrands().boundary( qp, domainValue( u, qp ) );
409 Hybrid::forEach( RangeValueIndices(), [ &qp, &w, &integrand, weight ] (
auto i ) {
410 std::get< i >( integrand ) *= weight;
411 w.axpy( qp, std::get< i >( integrand ) );
414 integrands().unbind();
417 template<
class Intersection,
class U,
class J >
418 void addLinearizedBoundaryIntegral (
const Intersection &intersection,
const U &u, DomainValueVectorType &phi, J &j )
const
420 if( !integrands().init( intersection ) )
423 const auto geometry = intersection.geometry();
424 const auto &domainBasis = j.domainBasisFunctionSet();
425 const auto &rangeBasis = j.rangeBasisFunctionSet();
427 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
428 const SurfaceQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder(u, domainBasis, rangeBasis )), SurfaceQuadratureType::INSIDE );
429 for(
const auto qp : quadrature )
431 const ctype weight = qp.weight() * geometry.integrationElement( qp.localPosition() );
433 values( domainBasis, qp, phi );
434 auto integrand = integrands().linearizedBoundary( qp, domainValue( u, qp ) );
436 for( std::size_t col = 0, cols = domainBasis.size(); col < cols; ++col )
438 LocalMatrixColumn< J > jCol( j, col );
439 RangeValueType intPhi = integrand( value( phi, col ) );
441 Hybrid::forEach( RangeValueIndices(), [ &qp, &jCol, &intPhi, weight ] (
auto i ) {
442 std::get< i >( intPhi ) *= weight;
443 jCol.axpy( qp, std::get< i >( intPhi ) );
447 integrands().unbind();
453 template<
bool conforming,
class Intersection,
class U,
class W >
454 void addSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut, W &wIn )
const
456 const auto geometry = intersection.geometry();
458 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
459 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
460 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( uIn, uOut, wIn)),
false );
461 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
463 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
465 const auto qpIn = quadrature.inside()[ qp ];
466 const auto qpOut = quadrature.outside()[ qp ];
467 std::pair< RangeValueType, RangeValueType > integrand = integrands().skeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
469 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &wIn, &integrand, weight ] (
auto i ) {
470 std::get< i >( integrand.first ) *= weight;
471 wIn.axpy( qpIn, std::get< i >( integrand.first ) );
476 template<
bool conforming,
class Intersection,
class U,
class W >
477 void addSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut, W &wIn, W &wOut )
const
479 const auto geometry = intersection.geometry();
480 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
481 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
482 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( uIn, uOut, wIn, wOut)),
false );
483 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
485 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
487 const auto qpIn = quadrature.inside()[ qp ];
488 const auto qpOut = quadrature.outside()[ qp ];
489 std::pair< RangeValueType, RangeValueType > integrand = integrands().skeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
491 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &wIn, &qpOut, &wOut, &integrand, weight ] (
auto i ) {
492 std::get< i >( integrand.first ) *= weight;
493 wIn.axpy( qpIn, std::get< i >( integrand.first ) );
495 std::get< i >( integrand.second ) *= weight;
496 wOut.axpy( qpOut, std::get< i >( integrand.second ) );
501 template<
bool conforming,
class Intersection,
class U,
class J >
502 void addLinearizedSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut,
503 DomainValueVectorType &phiIn, DomainValueVectorType &phiOut, J &jInIn, J &jOutIn )
const
505 const auto &domainBasisIn = jInIn.domainBasisFunctionSet();
506 const auto &domainBasisOut = jOutIn.domainBasisFunctionSet();
508 const auto &rangeBasisIn = jInIn.rangeBasisFunctionSet();
510 const int order =
std::max( maxOrder(uIn, uOut), maxOrder( domainBasisIn, domainBasisOut, rangeBasisIn ));
512 const auto geometry = intersection.geometry();
513 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
514 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
515 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(order),
false );
516 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
518 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
520 const auto qpIn = quadrature.inside()[ qp ];
521 const auto qpOut = quadrature.outside()[ qp ];
523 values( domainBasisIn, qpIn, phiIn );
524 values( domainBasisOut, qpOut, phiOut );
526 auto integrand = integrands().linearizedSkeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
527 for( std::size_t col = 0, cols = domainBasisIn.size(); col < cols; ++col )
529 LocalMatrixColumn< J > jInInCol( jInIn, col );
530 std::pair< RangeValueType, RangeValueType > intPhi = integrand.first( value( phiIn, col ) );
532 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jInInCol, &intPhi, weight ] (
auto i ) {
533 std::get< i >( intPhi.first ) *= weight;
534 jInInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
537 for( std::size_t col = 0, cols = domainBasisOut.size(); col < cols; ++col )
539 LocalMatrixColumn< J > jOutInCol( jOutIn, col );
540 std::pair< RangeValueType, RangeValueType > intPhi = integrand.second( value( phiOut, col ) );
542 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jOutInCol, &intPhi, weight ] (
auto i ) {
543 std::get< i >( intPhi.first ) *= weight;
544 jOutInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
550 template<
bool conforming,
class Intersection,
class U,
class J >
551 void addLinearizedSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut,
552 DomainValueVectorType &phiIn, DomainValueVectorType &phiOut, J &jInIn, J &jOutIn, J &jInOut, J &jOutOut )
const
554 const auto &domainBasisIn = jInIn.domainBasisFunctionSet();
555 const auto &domainBasisOut = jOutIn.domainBasisFunctionSet();
557 const auto &rangeBasisIn = jInIn.rangeBasisFunctionSet();
558 const auto &rangeBasisOut = jInOut.rangeBasisFunctionSet();
560 const int order =
std::max( maxOrder(uIn, uOut), maxOrder( domainBasisIn, domainBasisOut, rangeBasisIn, rangeBasisOut ));
562 const auto geometry = intersection.geometry();
563 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
564 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
565 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(order),
false );
566 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
568 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
570 const auto qpIn = quadrature.inside()[ qp ];
571 const auto qpOut = quadrature.outside()[ qp ];
573 values( domainBasisIn, qpIn, phiIn );
574 values( domainBasisOut, qpOut, phiOut );
576 auto integrand = integrands().linearizedSkeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
577 for( std::size_t col = 0, cols = domainBasisIn.size(); col < cols; ++col )
579 LocalMatrixColumn< J > jInInCol( jInIn, col );
580 LocalMatrixColumn< J > jInOutCol( jInOut, col );
581 std::pair< RangeValueType, RangeValueType > intPhi = integrand.first( value( phiIn, col ) );
583 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jInInCol, &qpOut, &jInOutCol, &intPhi, weight ] (
auto i ) {
584 std::get< i >( intPhi.first ) *= weight;
585 jInInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
587 std::get< i >( intPhi.second ) *= weight;
588 jInOutCol.axpy( qpOut, std::get< i >( intPhi.second ) );
591 for( std::size_t col = 0, cols = domainBasisOut.size(); col < cols; ++col )
593 LocalMatrixColumn< J > jOutInCol( jOutIn, col );
594 LocalMatrixColumn< J > jOutOutCol( jOutOut, col );
595 std::pair< RangeValueType, RangeValueType > intPhi = integrand.second( value( phiOut, col ) );
597 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jOutInCol, &qpOut, &jOutOutCol, &intPhi, weight ] (
auto i ) {
598 std::get< i >( intPhi.first ) *= weight;
599 jOutInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
601 std::get< i >( intPhi.second ) *= weight;
602 jOutOutCol.axpy( qpOut, std::get< i >( intPhi.second ) );
609 template<
class Intersection,
class U,
class... W >
610 void addSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut, W &... w )
const
612 if( !integrands().init( intersection ) )
615 if( intersection.conforming() )
616 addSkeletonIntegral< true >( intersection, uIn, uOut, w... );
618 addSkeletonIntegral< false >( intersection, uIn, uOut, w... );
619 integrands().unbind();
622 template<
class Intersection,
class U,
class... J >
623 void addLinearizedSkeletonIntegral (
const Intersection &intersection,
const U &uIn,
const U &uOut, DomainValueVectorType &phiIn, DomainValueVectorType &phiOut, J &... j )
const
625 if( !integrands().init( intersection ) )
628 if( intersection.conforming() )
629 addLinearizedSkeletonIntegral< true >( intersection, uIn, uOut, phiIn, phiOut, j... );
631 addLinearizedSkeletonIntegral< false >( intersection, uIn, uOut, phiIn, phiOut, j... );
632 integrands().unbind();
635 void setQuadratureOrders(
unsigned int interior,
unsigned int surface)
637 interiorQuadOrder_ = interior;
638 surfaceQuadOrder_ = surface;
641 IntegrandsType& model()
const
647 IntegrandsType& integrands()
const
652 template<
class Gr
idFunction,
class DiscreteFunction,
class Iterators,
class Functor >
653 void evaluate (
const GridFunction &u, DiscreteFunction &w,
const Iterators& iterators,
654 Functor& addLocalDofs, std::false_type )
const
656 typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
660 gridSizeInterior_ = 0;
662 const auto end = iterators.end();
663 for(
auto it = iterators.begin(); it != end; ++it )
668 const EntityType entity = *it ;
670 auto uGuard =
bindGuard( uLocal, entity );
671 auto wGuard =
bindGuard( wLocal, entity );
674 if( integrands().hasInterior() )
675 addInteriorIntegral( uLocal, wLocal );
677 if( integrands().hasBoundary() && entity.hasBoundaryIntersections() )
679 for(
const auto &intersection : intersections( gridPart(), entity ) )
681 if( intersection.boundary() )
682 addBoundaryIntegral( intersection, uLocal, wLocal );
688 addLocalDofs( entity, wLocal );
692 template<
class Gr
idFunction,
class DiscreteFunction,
class Iterators,
class Functor >
693 void evaluate (
const GridFunction &u, DiscreteFunction &w,
const Iterators& iterators,
694 Functor& addLocalDofs, std::true_type )
const
699 typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
702 gridSizeInterior_ = 0;
704 const auto &indexSet = gridPart().indexSet();
706 const auto end = iterators.end();
707 for(
auto it = iterators.begin(); it != end; ++it )
710 const EntityType inside = *it ;
715 auto uGuard =
bindGuard( uInside, inside );
716 auto wGuard =
bindGuard( wInside, inside );
719 if( integrands().hasInterior() )
720 addInteriorIntegral( uInside, wInside );
722 for(
const auto &intersection : intersections( gridPart(), inside ) )
726 if( intersection.neighbor() )
728 const EntityType outside = intersection.outside();
730 if( outside.partitionType() != InteriorEntity )
732 auto uOutGuard =
bindGuard( uOutside, outside );
733 addSkeletonIntegral( intersection, uInside, uOutside, wInside );
735 else if( indexSet.index( inside ) < indexSet.index( outside ) )
737 auto uOutGuard =
bindGuard( uOutside, outside );
738 auto wOutGuard =
bindGuard( wOutside, outside );
740 addSkeletonIntegral( intersection, uInside, uOutside, wInside, wOutside );
744 addLocalDofs( outside, wOutside );
747 else if( intersection.boundary() )
749 if( integrands().hasBoundary() )
750 addBoundaryIntegral( intersection, uInside, wInside );
754 addLocalDofs( inside, wInside );
758 template <
class Space>
761 typedef typename Space::EntityType EntityType;
762 template <
class Iterators>
763 InsideEntity(
const Space &space,
const Iterators& iterators)
764 : space_(space), dofThread_(space.size(),-1)
765 , thread_(MPIManager::thread())
767 const auto& mapper = space_.blockMapper();
768 for (
const auto &entity : space_)
770 int t=iterators.threadParallel(entity);
771 mapper.mapEach(entity, [
this, t ] (
int local,
auto global )
772 { dofThread_[global] = (dofThread_[global]==t || dofThread_[global]==-1)?
776 bool operator()(
const EntityType &entity)
const
778 bool needsLocking =
false;
779 space_.blockMapper().mapEach(entity,
780 [
this, &needsLocking ] (
int local,
auto global )
781 { needsLocking = (needsLocking || dofThread_[global]!=thread_); });
782 return !needsLocking;
785 std::vector<int> dofThread_;
789 template <
class DiscreteFunction>
790 struct AddLocalEvaluate
792 AddLocalEvaluate(DiscreteFunction &w)
794 template <
class LocalDofs>
795 void operator () (
const EntityType& entity,
const LocalDofs& wLocal )
const
797 w_.addLocalDofs( entity, wLocal.localDofVector() );
799 DiscreteFunction &w_;
802 template <
class DiscreteFunction>
803 struct AddLocalEvaluateLocked :
public AddLocalEvaluate<DiscreteFunction>
805 typedef AddLocalEvaluate<DiscreteFunction> BaseType;
807 std::shared_mutex& mutex_;
808 InsideEntity<typename DiscreteFunction::DiscreteFunctionSpaceType> inside_;
810 template <
class Iterators>
811 AddLocalEvaluateLocked(DiscreteFunction &w, std::shared_mutex& mtx,
const Iterators &iterators)
812 : BaseType(w), mutex_(mtx), inside_(w.space(),iterators) {}
814 template <
class LocalDofs>
815 void operator () (
const EntityType& entity,
const LocalDofs& wLocal )
const
820 std::shared_lock<std::shared_mutex> guard ( mutex_ );
821 BaseType::operator()( entity, wLocal );
826 std::lock_guard<std::shared_mutex> guard ( mutex_ );
827 BaseType::operator()( entity, wLocal );
832 template<
class Gr
idFunction,
class DiscreteFunction,
class Iterators,
class Functor >
833 void evaluate (
const GridFunction &u, DiscreteFunction &w,
const Iterators& iterators,
834 Functor& addLocalDofs )
const
836 static_assert( std::is_same< typename GridFunction::GridPartType, GridPartType >::value,
"Argument 'u' and Integrands must be defined on the same grid part." );
837 static_assert( std::is_same< typename DiscreteFunction::GridPartType, GridPartType >::value,
"Argument 'w' and Integrands must be defined on the same grid part." );
839 typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
843 if( integrands().hasSkeleton() )
844 evaluate( u, w, iterators, addLocalDofs, std::true_type() );
846 evaluate( u, w, iterators, addLocalDofs, std::false_type() );
850 template<
class Gr
idFunction,
class DiscreteFunction,
class Iterators >
851 void evaluate (
const GridFunction &u, DiscreteFunction &w,
const Iterators& iterators, std::shared_mutex& mtx )
const
853 AddLocalEvaluateLocked<DiscreteFunction> addLocalEvaluate(w,mtx,iterators);
854 evaluate( u, w, iterators, addLocalEvaluate );
857 template<
class Gr
idFunction,
class DiscreteFunction,
class Iterators >
858 void evaluate (
const GridFunction &u, DiscreteFunction &w,
const Iterators& iterators )
const
860 AddLocalEvaluate<DiscreteFunction> addLocalEvaluate(w);
861 evaluate( u, w, iterators, addLocalEvaluate );
864 template<
class T,
int length>
869 FiniteStack () : _f(0) {}
872 bool empty ()
const {
return _f <= 0; }
875 bool full ()
const {
return (_f >= length); }
878 void clear() { _f = 0; }
881 void push (
const T& t)
883 assert ( _f < length );
900 int size ()
const {
return _f; }
911 template <
class JacobianOperator>
912 struct AddLocalAssemble
914 typedef typename JacobianOperator::DomainSpaceType DomainSpaceType;
915 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
916 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
917 JacobianOperator &jOp_;
918 std::vector< TemporaryLocalMatrixType > jOpLocal_;
920 FiniteStack< TemporaryLocalMatrixType*, 12 > jOpLocalFinalized_;
921 FiniteStack< TemporaryLocalMatrixType*, 12 > jOpLocalFree_;
923 std::size_t locked, notLocked, timesLocked;
924 AddLocalAssemble(JacobianOperator& jOp)
926 , jOpLocal_(12, TemporaryLocalMatrixType(jOp_.domainSpace(), jOp_.rangeSpace()))
927 , jOpLocalFinalized_()
929 , locked(0), notLocked(0), timesLocked(0)
931 for(
auto& jOpLocal : jOpLocal_ )
932 jOpLocalFree_.push( &jOpLocal );
935 TemporaryLocalMatrixType& bind(
const EntityType& dE,
const EntityType& rE)
937 assert( ! jOpLocalFree_.empty() );
938 TemporaryLocalMatrixType& lop = *(jOpLocalFree_.pop());
944 void unbind(TemporaryLocalMatrixType &lop)
947 jOp_.addLocalMatrix( lop.domainEntity(), lop.rangeEntity(), lop );
949 jOpLocalFree_.push( &lop );
954 locked += jOpLocalFinalized_.size();
955 while ( ! jOpLocalFinalized_.empty() )
957 TemporaryLocalMatrixType &lop = *(jOpLocalFinalized_.pop());
958 jOp_.addLocalMatrix( lop.domainEntity(), lop.rangeEntity(), lop );
960 jOpLocalFree_.push( &lop );
965 template <
class JacobianOperator>
966 struct AddLocalAssembleLocked :
public AddLocalAssemble<JacobianOperator>
968 typedef AddLocalAssemble<JacobianOperator> BaseType;
969 typedef typename BaseType::TemporaryLocalMatrixType TemporaryLocalMatrixType;
970 using BaseType::jOpLocalFinalized_;
971 using BaseType::jOpLocalFree_;
973 std::shared_mutex& mutex_;
974 InsideEntity<typename JacobianOperator::DomainSpaceType> insideDomain_;
975 InsideEntity<typename JacobianOperator::RangeSpaceType> insideRange_;
977 template <
class Iterators>
978 AddLocalAssembleLocked(JacobianOperator &jOp, std::shared_mutex &mtx,
const Iterators &iterators)
981 , insideDomain_(jOp.domainSpace(),iterators)
982 , insideRange_(jOp.rangeSpace(),iterators)
988 ++BaseType::timesLocked;
989 std::lock_guard<std::shared_mutex> guard ( mutex_ );
990 BaseType::finalize();
993 TemporaryLocalMatrixType& bind(
const EntityType& dE,
const EntityType& rE)
995 if ( jOpLocalFree_.empty() )
999 return BaseType::bind(dE,rE);
1002 void unbind(TemporaryLocalMatrixType &lop)
1011 if ( insideDomain_(lop.domainEntity()) &&
1012 insideRange_(lop.rangeEntity()) )
1014 std::shared_lock<std::shared_mutex> guard ( mutex_ );
1015 BaseType::unbind(lop);
1019 jOpLocalFinalized_.push( &lop );
1024 template<
class Gr
idFunction,
class JacobianOperator,
class Iterators,
class Functor >
1025 void assemble (
const GridFunction &u, JacobianOperator &jOp,
const Iterators& iterators,
1026 Functor& addLocalMatrix, std::false_type )
const
1028 typedef typename JacobianOperator::DomainSpaceType DomainSpaceType;
1029 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
1031 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
1033 const std::size_t maxNumLocalDofs = jOp.domainSpace().blockMapper().maxNumDofs() * jOp.domainSpace().localBlockSize;
1034 DomainValueVectorType phi = makeDomainValueVector( maxNumLocalDofs );
1038 gridSizeInterior_ = 0;
1040 const auto end = iterators.end();
1041 for(
auto it = iterators.begin(); it != end; ++it )
1044 ++gridSizeInterior_;
1046 const EntityType entity = *it;
1048 auto guard =
bindGuard( uLocal, entity );
1050 TemporaryLocalMatrixType& jOpLocal = addLocalMatrix.bind(entity, entity );
1052 if( integrands().hasInterior() )
1053 addLinearizedInteriorIntegral( uLocal, phi, jOpLocal );
1055 if( integrands().hasBoundary() && entity.hasBoundaryIntersections() )
1057 for(
const auto &intersection : intersections( gridPart(), entity ) )
1059 if( intersection.boundary() )
1060 addLinearizedBoundaryIntegral( intersection, uLocal, phi, jOpLocal );
1064 addLocalMatrix.unbind(jOpLocal);
1066 addLocalMatrix.finalize();
1069 template<
class Gr
idFunction,
class JacobianOperator,
class Iterators,
class Functor >
1070 void assemble (
const GridFunction &u, JacobianOperator &jOp,
const Iterators& iterators,
1071 Functor& addLocalMatrix, std::true_type )
const
1073 typedef typename JacobianOperator::DomainSpaceType DomainSpaceType;
1074 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
1076 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
1078 const std::size_t maxNumLocalDofs = jOp.domainSpace().blockMapper().maxNumDofs() * jOp.domainSpace().localBlockSize;
1079 DomainValueVectorType phiIn = makeDomainValueVector( maxNumLocalDofs );
1080 DomainValueVectorType phiOut = makeDomainValueVector( maxNumLocalDofs );
1085 gridSizeInterior_ = 0;
1087 const auto &indexSet = gridPart().indexSet();
1089 const auto end = iterators.end();
1090 for(
auto it = iterators.begin(); it != end; ++it )
1093 ++gridSizeInterior_;
1095 const EntityType inside = *it;
1097 auto uiGuard =
bindGuard( uIn, inside );
1099 TemporaryLocalMatrixType& jOpInIn = addLocalMatrix.bind(inside, inside );
1101 if( integrands().hasInterior() )
1102 addLinearizedInteriorIntegral( uIn, phiIn, jOpInIn );
1104 for(
const auto &intersection : intersections( gridPart(), inside ) )
1108 if( intersection.neighbor() )
1110 const EntityType &outside = intersection.outside();
1112 TemporaryLocalMatrixType &jOpOutIn = addLocalMatrix.bind( outside, inside );
1114 auto uoGuard =
bindGuard( uOut, outside );
1116 if( outside.partitionType() != InteriorEntity )
1117 addLinearizedSkeletonIntegral( intersection, uIn, uOut, phiIn, phiOut, jOpInIn, jOpOutIn );
1118 else if( indexSet.index( inside ) < indexSet.index( outside ) )
1120 TemporaryLocalMatrixType &jOpInOut = addLocalMatrix.bind( inside, outside );
1121 TemporaryLocalMatrixType &jOpOutOut = addLocalMatrix.bind( outside, outside );
1123 addLinearizedSkeletonIntegral( intersection, uIn, uOut, phiIn, phiOut, jOpInIn, jOpOutIn, jOpInOut, jOpOutOut );
1125 addLocalMatrix.unbind(jOpInOut);
1126 addLocalMatrix.unbind(jOpOutOut);
1129 addLocalMatrix.unbind(jOpOutIn);
1131 else if( intersection.boundary() )
1133 if( integrands().hasBoundary() )
1134 addLinearizedBoundaryIntegral( intersection, uIn, phiIn, jOpInIn );
1137 addLocalMatrix.unbind(jOpInIn);
1139 addLocalMatrix.finalize();
1142 template<
class Gr
idFunction,
class JacobianOperator,
class Iterators,
class Functor >
1143 void assemble (
const GridFunction &u, JacobianOperator &jOp,
const Iterators& iterators,
1144 Functor& addLocalMatrix,
int )
const
1146 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
1148 static_assert( std::is_same< typename GridFunction::GridPartType, GridPartType >::value,
"Argument 'u' and Integrands must be defined on the same grid part." );
1149 static_assert( std::is_same< typename JacobianOperator::DomainSpaceType::GridPartType, GridPartType >::value,
"Argument 'jOp' and Integrands must be defined on the same grid part." );
1150 static_assert( std::is_same< typename JacobianOperator::RangeSpaceType::GridPartType, GridPartType >::value,
"Argument 'jOp' and Integrands must be defined on the same grid part." );
1156 if( integrands().hasSkeleton() )
1157 assemble( u, jOp, iterators, addLocalMatrix, std::true_type() );
1159 assemble( u, jOp, iterators, addLocalMatrix, std::false_type() );
1163 template<
class Gr
idFunction,
class JacobianOperator,
class Iterators >
1164 void assemble (
const GridFunction &u, JacobianOperator &jOp,
const Iterators& iterators, std::shared_mutex& mtx)
const
1166 AddLocalAssembleLocked<JacobianOperator> addLocalAssemble( jOp, mtx, iterators);
1167 assemble( u, jOp, iterators, addLocalAssemble, 10 );
1169 std::lock_guard guard ( mtx );
1171 << addLocalAssemble.locked <<
" " << addLocalAssemble.notLocked <<
" "
1172 << addLocalAssemble.timesLocked << std::endl;
1176 template<
class Gr
idFunction,
class JacobianOperator,
class Iterators >
1177 void assemble (
const GridFunction &u, JacobianOperator &jOp,
const Iterators& iterators )
const
1179 AddLocalAssemble<JacobianOperator> addLocalAssemble(jOp);
1180 assemble( u, jOp, iterators, addLocalAssemble, 10 );
1184 const GridPartType &gridPart ()
const {
return gridPart_; }
1186 unsigned int interiorQuadratureOrder(
unsigned int order)
const {
return interiorQuadOrder_ == 0 ? defaultInteriorOrder_(order) : interiorQuadOrder_; }
1187 unsigned int surfaceQuadratureOrder(
unsigned int order)
const {
return surfaceQuadOrder_ == 0 ? defaultSurfaceOrder_ (order) : surfaceQuadOrder_; }
1189 std::size_t gridSizeInterior()
const {
return gridSizeInterior_; }
1193 int maxOrder(
const U& u )
const
1195 return CallOrder< U > :: order( u );
1198 template<
class U,
class W >
1199 int maxOrder(
const U& u,
const W& w )
const
1201 return std::max( maxOrder( u ), maxOrder( w ) );
1204 template<
class U,
class V,
class W >
1205 int maxOrder(
const U& u,
const V& v,
const W& w )
const
1207 return std::max( maxOrder( u, v ), maxOrder( w ) );
1210 template<
class U,
class V,
class W,
class X >
1211 int maxOrder(
const U& u,
const V& v,
const W& w,
const X& x )
const
1213 return std::max( maxOrder( u, v ), maxOrder( w, x) );
1217 const GridPartType &gridPart_;
1219 mutable IntegrandsType integrands_;
1221 mutable std::function<int(
const int)> defaultInteriorOrder_;
1222 mutable std::function<int(
const int)> defaultSurfaceOrder_;
1224 unsigned int interiorQuadOrder_;
1225 unsigned int surfaceQuadOrder_;
1227 mutable std::vector< RangeValueVectorType > rangeValues_;
1228 mutable RangeValueVectorType values_;
1229 mutable DomainValueVectorType basisValues_;
1230 mutable DomainValueVectorType domainValues_;
1232 mutable std::size_t gridSizeInterior_;
1243 template<
class Integrands,
class DomainFunction,
class RangeFunction = DomainFunction >
1245 :
public virtual Operator< DomainFunction, RangeFunction >
1252 static_assert( std::is_same< typename DomainFunctionType::GridPartType, typename RangeFunctionType::GridPartType >::value,
"DomainFunction and RangeFunction must be defined on the same grid part." );
1257 template<
class... Args >
1271 std::cout <<
"GalerkinOperator::setCommunicate: communicate was disabled!" << std::endl;
1278 for(
size_t i=0; i<size; ++i )
1287 template<
class Gr
idFunction >
1309 for(
size_t i=0; i<size; ++i )
1317 template <
class Gr
idFunction >
1323 std::shared_mutex mutex;
1325 auto doEval = [
this, &u, &w, &mutex] ()
1332 MPIManager :: run ( doEval );
1364 template<
class Integrands,
class JacobianOperator >
1366 :
public GalerkinOperator< Integrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType >,
1381 template<
class... Args >
1397 template<
class Gr
idFunction >
1435 template <
class Gr
idFunction >
1444 std::shared_mutex mutex;
1446 auto doAssemble = [
this, &u, &jOp, &mutex] ()
1453 MPIManager :: run ( doAssemble );
1470 jOp.flushAssembly();
1490 template<
class Integrands,
class DomainFunction,
class RangeFunction >
1501 template<
class... Args >
1512 template <
class LinearOperator,
class ModelIntegrands >
1527 template<
class Gr
idFunction >
1533 template<
class Gr
idFunction >
1536 (*this).jacobian( u, jOp );
1545 template <
class O,
bool addDirichletBC>
1546 struct DirichletBlockSelector {
using type = void; };
1548 struct DirichletBlockSelector<O,true> {
using type =
typename O::DirichletBlockVector; };
1549 template<
class Integrands,
class LinearOperator,
class InverseOperator,
bool addDirichletBC,
1550 template <
class,
class>
class DifferentiableGalerkinOperatorImpl = DifferentiableGalerkinOperator >
1551 struct GalerkinSchemeImpl
1553 typedef InverseOperator InverseOperatorType;
1554 typedef Integrands ModelType;
1555 using DifferentiableOperatorType = std::conditional_t< addDirichletBC,
1557 DifferentiableGalerkinOperatorImpl< Integrands, LinearOperator > >;
1558 using DirichletBlockVector =
typename DirichletBlockSelector<
1560 DifferentiableGalerkinOperatorImpl< Integrands, LinearOperator >>,
1561 addDirichletBC>::type;
1563 typedef typename DifferentiableOperatorType::DomainFunctionType DomainFunctionType;
1564 typedef typename DifferentiableOperatorType::RangeFunctionType RangeFunctionType;
1565 typedef typename DifferentiableOperatorType::JacobianOperatorType LinearOperatorType;
1566 typedef typename DifferentiableOperatorType::JacobianOperatorType JacobianOperatorType;
1568 typedef RangeFunctionType DiscreteFunctionType;
1569 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeFunctionSpaceType;
1570 typedef typename RangeFunctionType::DiscreteFunctionSpaceType DomainFunctionSpaceType;
1571 typedef typename RangeFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
1573 typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
1574 typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
1577 typedef InverseOperator LinearInverseOperatorType;
1582 SolverInfo (
bool converged,
int linearIterations,
int nonlinearIterations )
1583 : converged( converged ), linearIterations( linearIterations ), nonlinearIterations( nonlinearIterations )
1587 int linearIterations, nonlinearIterations;
1590 GalerkinSchemeImpl (
const DiscreteFunctionSpaceType &dfSpace,
1591 const Integrands &integrands,
1593 : dfSpace_( dfSpace ),
1594 fullOperator_( dfSpace, dfSpace,
std::move( integrands ) ),
1598 void setQuadratureOrders(
unsigned int interior,
unsigned int surface) { fullOperator().setQuadratureOrders(interior,surface); }
1600 const DifferentiableOperatorType &fullOperator()
const {
return fullOperator_; }
1601 DifferentiableOperatorType &fullOperator() {
return fullOperator_; }
1603 void constraint ( DiscreteFunctionType &u )
const {}
1605 template<
class Gr
idFunction >
1606 void operator() (
const GridFunction &u, DiscreteFunctionType &w )
const
1608 fullOperator()( u, w );
1611 void setErrorMeasure(ErrorMeasureType &errorMeasure)
const
1613 invOp_.setErrorMeasure(errorMeasure);
1616 SolverInfo solve (
const DiscreteFunctionType &rhs, DiscreteFunctionType &solution )
const
1618 DiscreteFunctionType rhs0 = rhs;
1619 setZeroConstraints( rhs0 );
1620 setModelConstraints( solution );
1622 invOp_.bind(fullOperator());
1623 invOp_( rhs0, solution );
1625 return SolverInfo( invOp_.converged(), invOp_.linearIterations(), invOp_.iterations() );
1628 SolverInfo solve ( DiscreteFunctionType &solution )
const
1630 DiscreteFunctionType bnd( solution );
1632 setModelConstraints( solution );
1633 invOp_.bind(fullOperator());
1634 invOp_( bnd, solution );
1636 return SolverInfo( invOp_.converged(), invOp_.linearIterations(), invOp_.iterations() );
1639 template<
class Gr
idFunction >
1640 void jacobian(
const GridFunction &ubar, LinearOperatorType &linearOp)
const
1642 fullOperator().jacobian( ubar, linearOp );
1645 const DiscreteFunctionSpaceType &space ()
const {
return dfSpace_; }
1646 const GridPartType &gridPart ()
const {
return space().gridPart(); }
1647 ModelType &model()
const {
return fullOperator().model(); }
1649 void setConstraints( DomainFunctionType &u )
const
1651 if constexpr (addDirichletBC)
1652 fullOperator().setConstraints( u );
1654 void setConstraints(
const typename DiscreteFunctionType::RangeType &value, DiscreteFunctionType &u )
const
1656 if constexpr (addDirichletBC)
1657 fullOperator().setConstraints( value, u );
1659 void setConstraints(
const DiscreteFunctionType &u, DiscreteFunctionType &v )
const
1661 if constexpr (addDirichletBC)
1662 fullOperator().setConstraints( u, v );
1664 void subConstraints(
const DiscreteFunctionType &u, DiscreteFunctionType &v )
const
1666 if constexpr (addDirichletBC)
1667 fullOperator().subConstraints( u, v );
1669 const auto& dirichletBlocks()
const
1671 if constexpr (addDirichletBC)
1672 return fullOperator().dirichletBlocks();
1676 void setZeroConstraints( DiscreteFunctionType &u )
const
1678 if constexpr (addDirichletBC)
1679 fullOperator().setConstraints(
typename DiscreteFunctionType::RangeType(0), u );
1681 void setModelConstraints( DiscreteFunctionType &u )
const
1683 if constexpr (addDirichletBC)
1684 fullOperator().setConstraints( u );
1686 const DiscreteFunctionSpaceType &dfSpace_;
1687 DifferentiableOperatorType fullOperator_;
1688 mutable NewtonOperatorType invOp_;
1696 template<
class Integrands,
class LinearOperator,
class InverseOperator,
bool addDirichletBC >
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
Definition: bindguard.hh:11
Impl::GalerkinSchemeImpl< Integrands, LinearOperator, InverseOperator, addDirichletBC, DifferentiableGalerkinOperator > GalerkinScheme
Definition: galerkin.hh:1698
typename Impl::ConstLocalFunction< GridFunction >::Type ConstLocalFunction
Definition: const.hh:604
BasicParameterReader< std::function< const std::string *(const std::string &, const std::string *) > > ParameterReader
Definition: reader.hh:315
static auto bindGuard(Object &object, Args &&... args) -> std::enable_if_t< isBindable< Object, Args... >::value, BindGuard< Object > >
Definition: bindguard.hh:67
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
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 ParameterContainer & container()
Definition: io/parameter.hh:193
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
Exception thrown when a code segment that is supposed to be only accessed in single thread mode is ac...
Definition: mpimanager.hh:43
static int thread()
return thread number
Definition: mpimanager.hh:424
static int numThreads()
return number of current threads
Definition: mpimanager.hh:421
void update()
update internal list of iterators
Definition: threaditerator.hh:88
size_t size() const
return number of threads
Definition: threadsafevalue.hh:48
operator providing a Jacobian through automatic differentiation
Definition: automaticdifferenceoperator.hh:86
abstract differentiable operator
Definition: differentiableoperator.hh:29
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
void update()
clear previously computed entries such that a re-compute happens when used again
Definition: stencil.hh:136
Definition: dirichletwrapper.hh:27
Definition: galerkin.hh:1246
GalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: galerkin.hh:1258
const GridPartType & gridPart() const
Definition: galerkin.hh:1293
void evaluate(const GridFunction &u, RangeFunctionType &w) const
Definition: galerkin.hh:1318
ThreadIterator< GridPartType > ThreadIteratorType
Definition: galerkin.hh:1255
ModelType & model() const
Definition: galerkin.hh:1297
Impl::GalerkinOperator< Integrands > GalerkinOperatorImplType
Definition: galerkin.hh:1250
std::size_t gridSizeInterior_
Definition: galerkin.hh:1355
Integrands DirichletModelType
Definition: galerkin.hh:1296
ThreadSafeValue< GalerkinOperatorImplType > impl_
Definition: galerkin.hh:1353
DomainFunction DomainFunctionType
Definition: galerkin.hh:1247
ThreadIteratorType iterators_
Definition: galerkin.hh:1352
const GalerkinOperatorImplType & impl() const
Definition: galerkin.hh:1298
RangeFunction RangeFunctionType
Definition: galerkin.hh:1248
virtual void operator()(const DomainFunctionType &u, RangeFunctionType &w) const final override
application operator
Definition: galerkin.hh:1282
std::size_t gatherGridSizeInterior() const
Definition: galerkin.hh:1304
bool communicate_
Definition: galerkin.hh:1356
Integrands ModelType
Definition: galerkin.hh:1295
std::size_t gridSizeInterior() const
Definition: galerkin.hh:1300
RangeFunctionType::GridPartType GridPartType
Definition: galerkin.hh:1254
void setQuadratureOrders(unsigned int interior, unsigned int surface)
Definition: galerkin.hh:1275
void setCommunicate(const bool communicate)
Definition: galerkin.hh:1266
Definition: galerkin.hh:1368
int domainSpaceSequence_
Definition: galerkin.hh:1482
BaseType::DomainFunctionType DomainFunctionType
Definition: galerkin.hh:1374
BaseType::GridPartType GridPartType
Definition: galerkin.hh:1379
void jacobian(const GridFunction &u, JacobianOperatorType &jOp) const
Definition: galerkin.hh:1398
RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType
Definition: galerkin.hh:1377
virtual void jacobian(const DomainFunctionType &u, JacobianOperatorType &jOp) const final override
obtain linearization
Definition: galerkin.hh:1392
void prepare(JacobianOperatorType &jOp) const
Definition: galerkin.hh:1415
BaseType::RangeFunctionType RangeFunctionType
Definition: galerkin.hh:1375
int rangeSpaceSequence_
Definition: galerkin.hh:1482
std::size_t gridSizeInterior_
Definition: galerkin.hh:1355
const RangeDiscreteFunctionSpaceType & rangeSpace() const
Definition: galerkin.hh:1407
DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: galerkin.hh:1376
const RangeDiscreteFunctionSpaceType & rSpace_
Definition: galerkin.hh:1478
ThreadIteratorType iterators_
Definition: galerkin.hh:1352
const GalerkinOperatorImplType & impl() const
Definition: galerkin.hh:1298
std::size_t gatherGridSizeInterior() const
Definition: galerkin.hh:1304
DiagonalAndNeighborStencil< DomainDiscreteFunctionSpaceType, RangeDiscreteFunctionSpaceType > stencilDAN_
Definition: galerkin.hh:1480
JacobianOperator JacobianOperatorType
Definition: galerkin.hh:1372
void assemble(const GridFunction &u, JacobianOperatorType &jOp) const
Definition: galerkin.hh:1436
const DomainDiscreteFunctionSpaceType & dSpace_
Definition: galerkin.hh:1477
const DomainDiscreteFunctionSpaceType & domainSpace() const
Definition: galerkin.hh:1403
DifferentiableGalerkinOperator(const DomainDiscreteFunctionSpaceType &dSpace, const RangeDiscreteFunctionSpaceType &rSpace, Args &&... args)
Definition: galerkin.hh:1382
DiagonalStencil< DomainDiscreteFunctionSpaceType, RangeDiscreteFunctionSpaceType > stencilD_
Definition: galerkin.hh:1481
Definition: galerkin.hh:1494
BaseType::GridPartType GridPartType
Definition: galerkin.hh:1499
AutomaticDifferenceGalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: galerkin.hh:1502
Definition: galerkin.hh:1515
ModelDifferentiableGalerkinOperator(ModelType &model, const DiscreteFunctionSpaceType &dfSpace)
Definition: galerkin.hh:1523
ModelIntegrands::ModelType ModelType
Definition: galerkin.hh:1518
LinearOperator::DomainFunctionType RangeFunctionType
Definition: galerkin.hh:1520
void apply(const GridFunction &u, RangeFunctionType &w) const
Definition: galerkin.hh:1528
DifferentiableGalerkinOperator< ModelIntegrands, LinearOperator > BaseType
Definition: galerkin.hh:1516
LinearOperator::RangeSpaceType DiscreteFunctionSpaceType
Definition: galerkin.hh:1521
void apply(const GridFunction &u, LinearOperator &jOp) const
Definition: galerkin.hh:1534
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