1#ifndef DUNE_FEM_SCHEMES_INTEGRANDS_HH
2#define DUNE_FEM_SCHEMES_INTEGRANDS_HH
12#include <dune/common/ftraits.hh>
14#include <dune/fempy/quadrature/cachingpoint.hh>
15#include <dune/fempy/quadrature/elementpoint.hh>
31 namespace IntegrandsTraits
34 template<
class Integrands >
35 using DomainValueType =
typename std::decay_t< decltype( std::ref( std::declval< const Integrands & >() ).
get() ) >::DomainValueType;
37 template<
class Integrands >
38 using RangeValueType =
typename std::decay_t< decltype( std::ref( std::declval< const Integrands & >() ).
get() ) >::RangeValueType;
40 template<
class Integrands >
41 using GridPartType =
typename std::decay_t< decltype( std::ref( std::declval< const Integrands & >() ).
get() ) >::GridPartType;
44 template<
class Integrands >
45 using EntityType =
typename GridPartType< Integrands >::template Codim< 0 >::EntityType;
47 template<
class Integrands >
48 using IntersectionType =
typename GridPartType< Integrands >::IntersectionType;
51 template<
class Integrands >
52 using InteriorQuadratureType = CachingQuadrature< GridPartType< Integrands >, 0 >;
54 template<
class Integrands >
55 using SurfaceQuadratureType = CachingQuadrature< GridPartType< Integrands >, 1 >;
58 template<
class Integrands >
59 using InteriorQuadraturePointType = QuadraturePointWrapper< InteriorQuadratureType< Integrands > >;
61 template<
class Integrands >
62 using SurfaceQuadraturePointType = QuadraturePointWrapper< SurfaceQuadratureType< Integrands > >;
65 template<
class Integrands >
66 std::true_type interior (
const Integrands &,
decltype( std::declval< const Integrands & >().interior( std::declval<
const InteriorQuadraturePointType< Integrands > & >(), std::declval<
const DomainValueType< Integrands > & >() ) ) * =
nullptr );
68 std::false_type interior ( ... );
70 template< class Integrands, std::enable_if_t< std::is_same< decltype( std::declval< const Integrands & >().hasInterior() ),
bool >::value,
int > = 0 >
71 std::true_type hasInterior (
const Integrands & );
73 std::false_type hasInterior ( ... );
76 template<
class Integrands >
77 std::true_type boundary (
const Integrands &,
decltype( std::declval< const Integrands & >().boundary( std::declval<
const SurfaceQuadraturePointType< Integrands > & >(), std::declval<
const DomainValueType< Integrands > & >() ) ) * =
nullptr );
79 std::false_type boundary ( ... );
81 template< class Integrands, std::enable_if_t< std::is_same< decltype( std::declval< const Integrands & >().hasBoundary() ),
bool >::value,
int > = 0 >
82 std::true_type hasBoundary (
const Integrands & );
84 std::false_type hasBoundary ( ... );
87 template<
class Integrands >
88 std::true_type skeleton (
const Integrands &,
decltype( std::declval< const Integrands & >().skeleton( std::declval<
const SurfaceQuadraturePointType< Integrands > & >(), std::declval<
const DomainValueType< Integrands > & >(), std::declval<
const SurfaceQuadraturePointType< Integrands > & >(), std::declval<
const DomainValueType< Integrands > & >() ) ) * =
nullptr );
90 std::false_type skeleton ( ... );
92 template< class Integrands, std::enable_if_t< std::is_same< decltype( std::declval< const Integrands & >().hasSkeleton() ),
bool >::value,
int > = 0 >
93 std::true_type hasSkeleton (
const Integrands & );
95 std::false_type hasSkeleton ( ... );
97 template<
class Integrands >
98 using Get = std::decay_t< decltype( std::ref( std::declval< const Integrands & >() ).
get() ) >;
100 template<
class Integrands >
101 using RRangeType =
typename Get<Integrands>::RRangeType;
102 template<
class Integrands >
103 using DirichletComponentType =
typename Get<Integrands>::DirichletComponentType;
108 template<
class Integrands >
113 typedef Impl::IntegrandsTraits::GridPartType< Integrands >
GridPartType;
115 typedef Impl::IntegrandsTraits::EntityType< Integrands >
EntityType;
118 static const bool interior =
decltype( Impl::IntegrandsTraits::interior( std::declval< const Integrands & >() ) )::value;
119 static const bool hasInterior =
decltype( Impl::IntegrandsTraits::hasInterior( std::declval< const Integrands & >() ) )::value;
121 static const bool boundary =
decltype( Impl::IntegrandsTraits::boundary( std::declval< const Integrands & >() ) )::value;
122 static const bool hasBoundary =
decltype( Impl::IntegrandsTraits::hasBoundary( std::declval< const Integrands & >() ) )::value;
124 static const bool skeleton =
decltype( Impl::IntegrandsTraits::skeleton( std::declval< const Integrands & >() ) )::value;
125 static const bool hasSkeleton =
decltype( Impl::IntegrandsTraits::hasSkeleton( std::declval< const Integrands & >() ) )::value;
127 static_assert( (!
hasInterior ||
interior),
"Existence of method 'hasInterior' implies existence of method interior." );
128 static_assert( (!
hasBoundary ||
boundary),
"Existence of method 'hasBoundary' implies existence of method boundary." );
129 static_assert( (!
hasSkeleton ||
skeleton),
"Existence of method 'hasSkeleton' implies existence of method skeleton." );
133 typedef Impl::IntegrandsTraits::RRangeType< Integrands >
RRangeType;
142 template<
class Integrands >
153 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasInterior,
int > = 0 >
159 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasInterior,
int > = 0 >
165 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::interior,
int > = 0 >
171 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::interior,
int > = 0 >
177 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::interior,
int > = 0 >
183 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::interior,
int > = 0 >
189 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasBoundary,
int > = 0 >
195 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasBoundary,
int > = 0 >
201 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::boundary,
int > = 0 >
207 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::boundary,
int > = 0 >
213 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::boundary,
int > = 0 >
219 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::boundary,
int > = 0 >
225 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasSkeleton,
int > = 0 >
231 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasSkeleton,
int > = 0 >
237 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::skeleton,
int > = 0 >
240 return integrands.skeleton( xIn, uIn, xOut, uOut );
243 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::skeleton,
int > = 0 >
249 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::skeleton,
int > = 0 >
252 return integrands.linearizedSkeleton( xIn, uIn, xOut, uOut );
255 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::skeleton,
int > = 0 >
259 return std::make_pair( zero, zero );
263 template<
class... Args >
276 template<
class Po
int >
282 template<
class Po
int >
285 return linearizedInterior(
integrands(), x, u );
290 template<
class Po
int >
296 template<
class Po
int >
299 return linearizedBoundary(
integrands(), x, u );
304 template<
class Po
int >
307 return skeleton(
integrands(), xIn, uIn, xOut, uOut );
310 template<
class Po
int >
313 return linearizedSkeleton(
integrands(), xIn, uIn, xOut, uOut );
337 return integrands().isDirichletIntersection(inter,dirichletComponent);
339 template <
class Po
int>
355 template <
class FT,
int r>
356 struct GetDimRange<
Dune::FieldVector<FT,r>>
358 typedef Dune::FieldVector<FT,r> type;
359 static const int value = r;
361 template <
class FT,
int r,
int c>
362 struct GetDimRange<
Dune::FieldMatrix<FT,r,c>>
364 typedef Dune::FieldVector<FT,r> type;
365 static const int value = r;
367 template <
class FT,
int r,
int c>
368 struct GetDimRange<
Dune::Fem::ExplicitFieldVector<Dune::FieldMatrix<FT,c,c>,r>>
370 typedef Dune::FieldVector<FT,r> type;
371 static const int value = r;
375 template<
class Gr
idPart,
class DomainValue,
class RangeValue = DomainValue >
385 typedef typename GridPartType::template Codim< 0 >::EntityType
EntityType;
388 using RRangeType =
typename detail::GetDimRange<std::tuple_element_t<0,RangeValueType>>::type;
390 typedef typename EntityType::Geometry::LocalCoordinate
DomainType;
393 typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
395 typedef FemPy::CachingPoint< LocalCoordinateType, 0 > InteriorCachingPointType;
396 typedef FemPy::ElementPoint< LocalCoordinateType, 0 > InteriorElementPointType;
397 typedef FemPy::CachingPoint< LocalCoordinateType, 1 > SurfaceCachingPointType;
398 typedef FemPy::ElementPoint< LocalCoordinateType, 1 > SurfaceElementPointType;
410 using LinearizationPair = std::pair< Linearization< std::pair< R, R > >, Linearization< std::pair< R, R > > >;
414 virtual ~Interface () {}
415 virtual Interface *clone ()
const = 0;
417 virtual bool init (
const EntityType &entity ) = 0;
419 virtual void unbind ( ) = 0;
421 virtual bool hasInterior ()
const = 0;
424 virtual Linearization< RangeValueType > linearizedInterior (
const InteriorCachingPointType &x,
const DomainValueType &u )
const = 0;
425 virtual Linearization< RangeValueType > linearizedInterior (
const InteriorElementPointType &x,
const DomainValueType &u )
const = 0;
427 virtual bool hasBoundary ()
const = 0;
430 virtual Linearization< RangeValueType > linearizedBoundary (
const SurfaceCachingPointType &x,
const DomainValueType &u )
const = 0;
431 virtual Linearization< RangeValueType > linearizedBoundary (
const SurfaceElementPointType &x,
const DomainValueType &u )
const = 0;
433 virtual bool hasSkeleton ()
const = 0;
434 virtual std::pair< RangeValueType, RangeValueType > skeleton (
const SurfaceCachingPointType &xIn,
const DomainValueType &uIn,
const SurfaceCachingPointType &xOut,
const DomainValueType &uOut )
const = 0;
435 virtual std::pair< RangeValueType, RangeValueType > skeleton (
const SurfaceElementPointType &xIn,
const DomainValueType &uIn,
const SurfaceElementPointType &xOut,
const DomainValueType &uOut )
const = 0;
436 virtual LinearizationPair< RangeValueType > linearizedSkeleton (
const SurfaceCachingPointType &xIn,
const DomainValueType &uIn,
const SurfaceCachingPointType &xOut,
const DomainValueType &uOut )
const = 0;
437 virtual LinearizationPair< RangeValueType > linearizedSkeleton (
const SurfaceElementPointType &xIn,
const DomainValueType &uIn,
const SurfaceElementPointType &xOut,
const DomainValueType &uOut )
const = 0;
439 virtual bool hasDirichletBoundary ()
const = 0;
444 template<
class Impl >
445 struct Implementation final
448 Implementation ( Impl impl ) : impl_(
std::move( impl ) )
451 virtual Interface *clone ()
const override {
return new Implementation( *
this ); }
453 virtual bool init (
const EntityType &entity )
override {
return impl().init( entity ); }
454 virtual bool init (
const IntersectionType &intersection )
override {
return impl().init( intersection ); }
455 virtual void unbind ( )
override { impl().unbind( ); }
457 virtual bool hasInterior ()
const override {
return impl().hasInterior(); }
458 virtual RangeValueType interior (
const InteriorCachingPointType &x,
const DomainValueType &u )
const override {
return impl().interior( asQP( x ), u ); }
459 virtual RangeValueType interior (
const InteriorElementPointType &x,
const DomainValueType &u )
const override {
return impl().interior( asQP( x ), u ); }
460 virtual Linearization< RangeValueType > linearizedInterior (
const InteriorCachingPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedInterior( asQP( x ), u ); }
461 virtual Linearization< RangeValueType > linearizedInterior (
const InteriorElementPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedInterior( asQP( x ), u ); }
463 virtual bool hasBoundary ()
const override {
return impl().hasBoundary(); }
464 virtual RangeValueType boundary (
const SurfaceCachingPointType &x,
const DomainValueType &u )
const override {
return impl().boundary( asQP( x ), u ); }
465 virtual RangeValueType boundary (
const SurfaceElementPointType &x,
const DomainValueType &u )
const override {
return impl().boundary( asQP( x ), u ); }
466 virtual Linearization< RangeValueType > linearizedBoundary (
const SurfaceCachingPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedBoundary( asQP( x ), u ); }
467 virtual Linearization< RangeValueType > linearizedBoundary (
const SurfaceElementPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedBoundary( asQP( x ), u ); }
469 virtual bool hasSkeleton ()
const override {
return impl().hasSkeleton(); }
470 virtual std::pair< RangeValueType, RangeValueType > skeleton (
const SurfaceCachingPointType &xIn,
const DomainValueType &uIn,
const SurfaceCachingPointType &xOut,
const DomainValueType &uOut )
const override {
return impl().skeleton( asQP( xIn ), uIn, asQP( xOut ), uOut ); }
471 virtual std::pair< RangeValueType, RangeValueType > skeleton (
const SurfaceElementPointType &xIn,
const DomainValueType &uIn,
const SurfaceElementPointType &xOut,
const DomainValueType &uOut )
const override {
return impl().skeleton( asQP( xIn ), uIn, asQP( xOut ), uOut ); }
472 virtual LinearizationPair< RangeValueType > linearizedSkeleton (
const SurfaceCachingPointType &xIn,
const DomainValueType &uIn,
const SurfaceCachingPointType &xOut,
const DomainValueType &uOut )
const override {
return impl().linearizedSkeleton( asQP( xIn ), uIn, asQP( xOut ), uOut ); }
473 virtual LinearizationPair< RangeValueType > linearizedSkeleton (
const SurfaceElementPointType &xIn,
const DomainValueType &uIn,
const SurfaceElementPointType &xOut,
const DomainValueType &uOut )
const override {
return impl().linearizedSkeleton( asQP( xIn ), uIn, asQP( xOut ), uOut ); }
475 virtual bool hasDirichletBoundary ()
const override {
return impl().hasDirichletBoundary(); }
476 virtual bool isDirichletIntersection(
const IntersectionType& inter,
DirichletComponentType &dirichletComponent )
const override {
return impl().isDirichletIntersection(inter,dirichletComponent); }
477 virtual void dirichlet(
int bndId,
const DomainType &x,
RRangeType &value)
const override { impl().dirichlet(bndId,x,value); }
480 const auto &impl ()
const {
return std::cref( impl_ ).get(); }
481 auto &impl () {
return std::ref( impl_ ).get(); }
486 template<
class Integrands >
487 using isVirtualized = std::is_same< std::decay_t< decltype( std::ref( std::declval< Integrands & >() ).
get() ) >, This >;
490 template<
class Integrands, std::enable_if_t< IntegrandsTraits< std::decay_t< Integrands > >::isFull && !isVirtualized< Integrands >::value,
int > = 0 >
492 : impl_( new Implementation< Integrands >(
std::move( integrands ) ) )
495 template< class Integrands, std::enable_if_t< !IntegrandsTraits< Integrands >::isFull && !isVirtualized< Integrands >::value,
int > = 0 >
507 impl_.reset( other ? other.impl().clone() :
nullptr );
512 explicit operator bool ()
const {
return static_cast< bool >( impl_ ); }
520 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
523 return impl().interior( InteriorCachingPointType( x ), u );
526 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
529 return impl().interior( InteriorElementPointType( x ), u );
532 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
535 return impl().linearizedInterior( InteriorCachingPointType( x ), u );
538 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
541 return impl().linearizedInterior( InteriorElementPointType( x ), u );
546 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
549 return impl().boundary( SurfaceCachingPointType( x ), u );
552 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
555 return impl().boundary( SurfaceElementPointType( x ), u );
558 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
561 return impl().linearizedBoundary( SurfaceCachingPointType( x ), u );
564 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
567 return impl().linearizedBoundary( SurfaceElementPointType( x ), u );
572 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
575 return impl().skeleton( SurfaceCachingPointType( xIn ), uIn, SurfaceCachingPointType( xOut ), uOut );
578 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
581 return impl().skeleton( SurfaceElementPointType( xIn ), uIn, SurfaceElementPointType( xOut ), uOut );
584 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
587 return impl().linearizedSkeleton( SurfaceCachingPointType( xIn ), uIn, SurfaceCachingPointType( xOut ), uOut );
590 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
593 return impl().linearizedSkeleton( SurfaceElementPointType( xIn ), uIn, SurfaceElementPointType( xOut ), uOut );
598 return impl().hasDirichletBoundary();
602 return impl().isDirichletIntersection(inter,dirichletComponent);
606 return impl().dirichlet(bndId,x,value);
610 const Interface &impl ()
const { assert( impl_ );
return *impl_; }
611 Interface &impl () { assert( impl_ );
return *impl_; }
613 std::unique_ptr< Interface > impl_;
621 template<
class Model >
627 typedef typename GridPartType::template Codim< 0 >::EntityType
EntityType;
642 return (intersection.boundary() &&
model().hasNeumanBoundary() &&
model().
init( intersection.inside() ));
647 template<
class Po
int >
651 model().source( x, std::get< 0 >( u ), std::get< 1 >( u ), source );
653 model().flux( x, std::get< 0 >( u ), std::get< 1 >( u ), dFlux );
654 return std::make_tuple( source, dFlux );
657 template<
class Po
int >
662 model().linSource( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), source );
664 model().linFlux( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), dFlux );
665 return std::make_tuple( source, dFlux );
669 template<
class Po
int >
673 model().alpha( x, std::get< 0 >( u ), alpha );
674 return std::make_tuple( alpha, 0 );
677 template<
class Po
int >
682 model().linAlpha( std::get< 0 >( u ), x, std::get< 0 >( phi ), alpha );
683 return std::make_tuple( alpha, 0 );
687 const Model &
model ()
const {
return *model_; }
690 const Model *model_ =
nullptr;
698 template<
class Model >
704 typedef typename GridPartType::template Codim< 0 >::EntityType
EntityType;
716 : model_( &
model ), penalty_( penalty )
721 intersection_ =
nullptr;
722 return model().init( entity );
727 intersection_ = &intersection;
728 if( intersection.boundary() )
730 const EntityType inside = intersection.inside();
731 beta_ = penalty_ * intersection.geometry().volume() / inside.geometry().volume();
734 else if( intersection.neighbor() )
736 const auto volIn = intersection.inside().geometry().volume();
737 const auto volOut = intersection.outside().geometry().volume();
738 beta_ = penalty_ * intersection.geometry().volume() /
std::min( volIn, volOut );
748 template<
class Po
int >
752 model().source( x, std::get< 0 >( u ), std::get< 1 >( u ), source );
754 model().flux( x, std::get< 0 >( u ), std::get< 1 >( u ), dFlux );
758 template<
class Po
int >
763 model().linSource( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), source );
765 model().linFlux( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), dFlux );
770 template<
class Po
int >
774 model().alpha( x, std::get< 0 >( u ), alpha );
778 template<
class Po
int >
783 model().linAlpha( std::get< 0 >( u ), x, std::get< 0 >( phi ), alpha );
788 template<
class Po
int >
791 const EntityType inside = intersection().inside();
792 const EntityType outside = intersection().outside();
795 const auto normal = intersection().unitOuterNormal( xIn.localPosition() );
798 std::get< 0 >( uJump ) = std::get< 0 >( uOut ) - std::get< 0 >( uIn );
799 for(
int i = 0; i < RangeType::dimension; ++i )
800 std::get< 1 >( uJump )[ i ].axpy( std::get< 0 >( uJump )[ i ], normal );
802 model().init( outside );
804 model().flux( xOut, std::get< 0 >( uOut ), std::get< 1 >( uJump ), dFluxPrimeOut );
805 model().flux( xOut, std::get< 0 >( uOut ), 0, dFluxOut );
806 dFluxPrimeOut -= dFluxOut;
808 model().flux( xOut, std::get< 0 >( uOut ), std::get< 1 >( uOut ), dFluxOut );
810 model().init( inside );
812 model().flux( xIn, std::get< 0 >( uIn ), std::get< 1 >( uJump ), dFluxPrimeIn );
813 model().flux( xIn, std::get< 0 >( uIn ), 0, dFluxIn );
814 dFluxPrimeIn -= dFluxIn;
816 model().flux( xIn, std::get< 0 >( uIn ), std::get< 1 >( uIn ), dFluxIn );
821 dFluxIn.usmv( -half, normal, int0 );
823 dFluxPrimeIn *= -half;
824 dFluxPrimeOut *= -half;
829 template<
class Po
int >
832 const auto normal = intersection().unitOuterNormal( xIn.localPosition() );
835 std::get< 0 >( uJump ) = std::get< 0 >( uOut ) - std::get< 0 >( uIn );
836 for(
int i = 0; i < RangeType::dimension; ++i )
837 std::get< 1 >( uJump )[ i ].axpy( std::get< 0 >( uJump )[ i ], normal );
839 auto intIn = [
this, xIn, uIn, xOut, uOut, normal, uJump ] (
const DomainValueType &phiIn ) {
840 const EntityType inside = intersection().inside();
841 const EntityType outside = intersection().outside();
846 std::get< 0 >( phiJump ) -= std::get< 0 >( phiIn );
847 for(
int i = 0; i < RangeType::dimension; ++i )
848 std::get< 1 >( phiJump )[ i ].axpy( std::get< 0 >( phiJump )[ i ], normal );
850 model().init( outside );
852 model().linFlux( std::get< 0 >( uOut ), std::get< 1 >( uJump ), xOut, 0, std::get< 1 >( phiJump ), dFluxPrimeOut );
854 model().init( inside );
856 model().linFlux( std::get< 0 >( uIn ), std::get< 1 >( uJump ), xIn, std::get< 0 >( phiIn ), std::get< 1 >( phiJump ), dFluxPrimeIn );
857 model().linFlux( std::get< 0 >( uIn ), 0, xIn, std::get< 0 >( phiIn ), 0, dFluxIn );
858 dFluxPrimeIn -= dFluxIn;
860 model().linFlux( std::get< 0 >( uIn ), std::get< 1 >( uIn ), xIn, std::get< 0 >( phiIn ), std::get< 1 >( phiIn ), dFluxIn );
862 RangeType int0 = std::get< 0 >( phiJump );
864 dFluxIn.usmv( -half, normal, int0 );
866 dFluxPrimeIn *= -half;
867 dFluxPrimeOut *= -half;
872 auto intOut = [
this, xIn, uIn, xOut, uOut, normal, uJump ] (
const DomainValueType &phiOut ) {
873 const EntityType inside = intersection().inside();
874 const EntityType outside = intersection().outside();
879 std::get< 0 >( phiJump ) = std::get< 0 >( phiOut );
880 for(
int i = 0; i < RangeType::dimension; ++i )
881 std::get< 1 >( phiJump )[ i ].axpy( std::get< 0 >( phiJump )[ i ], normal );
883 model().init( outside );
885 model().linFlux( std::get< 0 >( uOut ), std::get< 1 >( uJump ), xOut, std::get< 0 >( phiOut ), std::get< 1 >( phiJump ), dFluxPrimeOut );
886 model().linFlux( std::get< 0 >( uOut ), 0, xOut, std::get< 0 >( phiOut ), 0, dFluxOut );
887 dFluxPrimeOut -= dFluxOut;
889 model().linFlux( std::get< 0 >( uOut ), std::get< 1 >( uOut ), xOut, std::get< 0 >( phiOut ), std::get< 1 >( phiOut ), dFluxOut );
891 model().init( inside );
893 model().linFlux( std::get< 0 >( uIn ), std::get< 1 >( uJump ), xIn, 0, std::get< 1 >( phiJump ), dFluxPrimeIn );
895 RangeType int0 = std::get< 0 >( phiJump );
897 dFluxOut.usmv( -half, normal, int0 );
899 dFluxPrimeIn *= -half;
900 dFluxPrimeOut *= -half;
905 return std::make_pair( intIn, intOut );
908 const Model &
model ()
const {
return *model_; }
911 const IntersectionType &intersection ()
const { assert( intersection_ );
return *intersection_; }
913 const Model *model_ =
nullptr;
double min(const Dune::Fem::Double &v, const double p)
Definition: double.hh:953
Definition: bindguard.hh:11
std::tuple_element< i, Tuple >::type & get(Dune::TypeIndexedTuple< Tuple, Types > &tuple)
Definition: typeindexedtuple.hh:122
wrapper for a (Quadrature,int) pair
Definition: quadrature.hh:43
Definition: integrands.hh:110
static const bool boundary
Definition: integrands.hh:121
static const bool hasBoundary
Definition: integrands.hh:122
Impl::IntegrandsTraits::DomainValueType< Integrands > DomainValueType
Definition: integrands.hh:111
static const bool hasSkeleton
Definition: integrands.hh:125
Impl::IntegrandsTraits::EntityType< Integrands > EntityType
Definition: integrands.hh:115
Impl::IntegrandsTraits::RRangeType< Integrands > RRangeType
Definition: integrands.hh:133
static const bool interior
Definition: integrands.hh:118
Impl::IntegrandsTraits::GridPartType< Integrands > GridPartType
Definition: integrands.hh:113
Impl::IntegrandsTraits::IntersectionType< Integrands > IntersectionType
Definition: integrands.hh:116
static const bool hasInterior
Definition: integrands.hh:119
static const bool isFull
Definition: integrands.hh:131
Impl::IntegrandsTraits::RangeValueType< Integrands > RangeValueType
Definition: integrands.hh:112
static const bool skeleton
Definition: integrands.hh:124
Impl::IntegrandsTraits::DirichletComponentType< Integrands > DirichletComponentType
Definition: integrands.hh:134
Definition: integrands.hh:144
Integrands::type RealIntegrands
Definition: integrands.hh:317
FullIntegrands(Args &&... args)
Definition: integrands.hh:264
IntegrandsTraits< Integrands >::GridPartType GridPartType
Definition: integrands.hh:147
RealIntegrands rInt_
Definition: integrands.hh:326
void unbind()
Definition: integrands.hh:272
bool init(const IntersectionType &intersection)
Definition: integrands.hh:271
const RealIntegrands & integrands() const
Definition: integrands.hh:321
bool hasInterior() const
Definition: integrands.hh:274
bool hasSkeleton() const
Definition: integrands.hh:302
auto linearizedBoundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:297
std::pair< RangeValueType, RangeValueType > skeleton(const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:305
bool hasBoundary() const
Definition: integrands.hh:288
RealIntegrands & integrands()
Definition: integrands.hh:320
bool isDirichletIntersection(const IntersectionType &inter, DirichletComponentType &dirichletComponent) const
Definition: integrands.hh:335
IntegrandsTraits< Integrands >::IntersectionType IntersectionType
Definition: integrands.hh:150
IntegrandsTraits< Integrands >::EntityType EntityType
Definition: integrands.hh:149
bool hasDirichletBoundary() const
Definition: integrands.hh:331
bool init(const EntityType &entity)
Definition: integrands.hh:270
IntegrandsTraits< Integrands >::DirichletComponentType DirichletComponentType
Definition: integrands.hh:330
void dirichlet(int bndId, const Point &x, RRangeType &value) const
Definition: integrands.hh:340
IntegrandsTraits< Integrands >::RRangeType RRangeType
Definition: integrands.hh:329
RangeValueType interior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:277
auto linearizedSkeleton(const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:311
auto linearizedInterior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:283
RangeValueType boundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:291
IntegrandsTraits< Integrands >::DomainValueType DomainValueType
Definition: integrands.hh:145
IntegrandsTraits< Integrands >::RangeValueType RangeValueType
Definition: integrands.hh:146
Integrands integrands_
Definition: integrands.hh:325
Definition: integrands.hh:377
RangeValueType boundary(const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u) const
Definition: integrands.hh:547
DomainValue DomainValueType
Definition: integrands.hh:382
RangeValueType interior(const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u) const
Definition: integrands.hh:521
VirtualizedIntegrands & operator=(const This &other)
Definition: integrands.hh:505
auto linearizedBoundary(const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u) const
Definition: integrands.hh:559
GridPart GridPartType
Definition: integrands.hh:381
VirtualizedIntegrands(Integrands integrands)
Definition: integrands.hh:491
VirtualizedIntegrands(This &&)=default
std::array< int, RRangeType::dimension > DirichletComponentType
Definition: integrands.hh:389
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: integrands.hh:385
void unbind()
Definition: integrands.hh:516
typename detail::GetDimRange< std::tuple_element_t< 0, RangeValueType > >::type RRangeType
Definition: integrands.hh:388
VirtualizedIntegrands(const This &other)
Definition: integrands.hh:500
bool hasSkeleton() const
Definition: integrands.hh:570
auto linearizedSkeleton(const Fem::QuadraturePointWrapper< Quadrature > &xIn, const DomainValueType &uIn, const Fem::QuadraturePointWrapper< Quadrature > &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:585
EntityType::Geometry::LocalCoordinate DomainType
Definition: integrands.hh:390
bool init(const IntersectionType &intersection)
Definition: integrands.hh:515
bool isDirichletIntersection(const IntersectionType &inter, DirichletComponentType &dirichletComponent) const
Definition: integrands.hh:600
RangeValue RangeValueType
Definition: integrands.hh:383
bool hasInterior() const
Definition: integrands.hh:518
GridPartType::IntersectionType IntersectionType
Definition: integrands.hh:386
bool hasBoundary() const
Definition: integrands.hh:544
auto linearizedInterior(const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u) const
Definition: integrands.hh:533
std::pair< RangeValueType, RangeValueType > skeleton(const Fem::QuadraturePointWrapper< Quadrature > &xIn, const DomainValueType &uIn, const Fem::QuadraturePointWrapper< Quadrature > &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:573
bool hasDirichletBoundary() const
Definition: integrands.hh:596
bool init(const EntityType &entity)
Definition: integrands.hh:514
void dirichlet(int bndId, const DomainType &x, RRangeType &value) const
Definition: integrands.hh:604
Definition: integrands.hh:623
Model::RangeType RangeType
Definition: integrands.hh:630
std::tuple< RangeType, JacobianRangeType > DomainValueType
Definition: integrands.hh:633
auto linearizedInterior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:658
bool init(const EntityType &entity)
Definition: integrands.hh:638
void unbind()
Definition: integrands.hh:645
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: integrands.hh:627
RangeValueType boundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:670
bool init(const IntersectionType &intersection)
Definition: integrands.hh:640
Model ModelType
Definition: integrands.hh:624
GridPartType::IntersectionType IntersectionType
Definition: integrands.hh:628
Model::GridPartType GridPartType
Definition: integrands.hh:625
RangeValueType interior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:648
const Model & model() const
Definition: integrands.hh:687
ConservationLawModelIntegrands(const Model &model)
Definition: integrands.hh:636
std::tuple< RangeType, JacobianRangeType > RangeValueType
Definition: integrands.hh:634
Model::JacobianRangeType JacobianRangeType
Definition: integrands.hh:631
auto linearizedBoundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:678
Definition: integrands.hh:700
auto linearizedBoundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:779
Model::JacobianRangeType JacobianRangeType
Definition: integrands.hh:708
bool init(const IntersectionType &intersection)
Definition: integrands.hh:725
std::tuple< RangeType, JacobianRangeType > RangeValueType
Definition: integrands.hh:713
const Model & model() const
Definition: integrands.hh:908
GridPartType::IntersectionType IntersectionType
Definition: integrands.hh:705
FieldTraits< RangeType >::field_type RangeFieldType
Definition: integrands.hh:710
Model::RangeType RangeType
Definition: integrands.hh:707
RangeValueType boundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:771
RangeValueType interior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:749
Model::GridPartType GridPartType
Definition: integrands.hh:702
bool init(const EntityType &entity)
Definition: integrands.hh:719
auto linearizedSkeleton(const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:830
std::pair< RangeValueType, RangeValueType > skeleton(const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:789
auto linearizedInterior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:759
std::tuple< RangeType, JacobianRangeType > DomainValueType
Definition: integrands.hh:712
DGConservationLawModelIntegrands(const Model &model, RangeFieldType penalty)
Definition: integrands.hh:715
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: integrands.hh:704
Model ModelType
Definition: integrands.hh:701
void unbind()
Definition: integrands.hh:743