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 >
156 return integrands.hasInterior();
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 >
168 return integrands.interior( x, u );
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 >
178 static auto linearizedInterior (
const T &integrands,
const Point &x,
const DomainValueType &u )
180 return integrands.linearizedInterior( x, u );
183 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::interior,
int > = 0 >
184 static auto linearizedInterior (
const T &integrands,
const Point &x,
const DomainValueType &u )
189 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasBoundary,
int > = 0 >
192 return integrands.hasBoundary();
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 >
204 return integrands.boundary( x, u );
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 >
214 static auto linearizedBoundary (
const T &integrands,
const Point &x,
const DomainValueType &u )
216 return integrands.linearizedBoundary( x, u );
219 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::boundary,
int > = 0 >
220 static auto linearizedBoundary (
const T &integrands,
const Point &x,
const DomainValueType &u )
225 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasSkeleton,
int > = 0 >
228 return integrands.hasSkeleton();
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 >
238 static std::pair< RangeValueType, RangeValueType > skeleton (
const T &integrands,
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
240 return integrands.skeleton( xIn, uIn, xOut, uOut );
243 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::skeleton,
int > = 0 >
244 static std::pair< RangeValueType, RangeValueType > skeleton (
const T &integrands,
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
249 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::skeleton,
int > = 0 >
250 static auto linearizedSkeleton (
const T &integrands,
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
252 return integrands.linearizedSkeleton( xIn, uIn, xOut, uOut );
255 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::skeleton,
int > = 0 >
256 static auto linearizedSkeleton (
const T &integrands,
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
259 return std::make_pair( zero, zero );
263 template<
class... Args >
273 template<
class Po
int >
276 return interior( integrands(), x, u );
279 template<
class Po
int >
282 return linearizedInterior( integrands(), x, u );
287 template<
class Po
int >
290 return boundary( integrands(), x, u );
293 template<
class Po
int >
296 return linearizedBoundary( integrands(), x, u );
301 template<
class Po
int >
304 return skeleton( integrands(), xIn, uIn, xOut, uOut );
307 template<
class Po
int >
310 return linearizedSkeleton( integrands(), xIn, uIn, xOut, uOut );
314 decltype(
auto ) integrands () {
return std::ref(
integrands_ ).get(); }
315 decltype(
auto ) integrands ()
const {
return std::ref(
integrands_ ).get(); }
328 return integrands_.get().isDirichletIntersection(inter,dirichletComponent);
330 template <
class Po
int>
346 template <
class FT,
int r>
347 struct GetDimRange<
Dune::FieldVector<FT,r>>
349 typedef Dune::FieldVector<FT,r> type;
350 static const int value = r;
352 template <
class FT,
int r,
int c>
353 struct GetDimRange<
Dune::FieldMatrix<FT,r,c>>
355 typedef Dune::FieldVector<FT,r> type;
356 static const int value = r;
358 template <
class FT,
int r,
int c>
359 struct GetDimRange<
Dune::Fem::ExplicitFieldVector<Dune::FieldMatrix<FT,c,c>,r>>
361 typedef Dune::FieldVector<FT,r> type;
362 static const int value = r;
366 template<
class Gr
idPart,
class DomainValue,
class RangeValue = DomainValue >
376 typedef typename GridPartType::template Codim< 0 >::EntityType
EntityType;
379 using RRangeType =
typename detail::GetDimRange<std::tuple_element_t<0,RangeValueType>>::type;
381 typedef typename EntityType::Geometry::LocalCoordinate
DomainType;
384 typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
386 typedef FemPy::CachingPoint< LocalCoordinateType, 0 > InteriorCachingPointType;
387 typedef FemPy::ElementPoint< LocalCoordinateType, 0 > InteriorElementPointType;
388 typedef FemPy::CachingPoint< LocalCoordinateType, 1 > SurfaceCachingPointType;
389 typedef FemPy::ElementPoint< LocalCoordinateType, 1 > SurfaceElementPointType;
401 using LinearizationPair = std::pair< Linearization< std::pair< R, R > >, Linearization< std::pair< R, R > > >;
405 virtual ~Interface () {}
406 virtual Interface *clone ()
const = 0;
424 virtual std::pair< RangeValueType, RangeValueType >
skeleton (
const SurfaceCachingPointType &xIn,
const DomainValueType &uIn,
const SurfaceCachingPointType &xOut,
const DomainValueType &uOut )
const = 0;
425 virtual std::pair< RangeValueType, RangeValueType >
skeleton (
const SurfaceElementPointType &xIn,
const DomainValueType &uIn,
const SurfaceElementPointType &xOut,
const DomainValueType &uOut )
const = 0;
434 template<
class Impl >
435 struct Implementation final
438 Implementation ( Impl impl ) : impl_( std::move( impl ) ) {}
439 virtual Interface *clone ()
const override {
return new Implementation( *
this ); }
441 virtual bool init (
const EntityType &entity )
override {
return impl().init( entity ); }
442 virtual bool init (
const IntersectionType &intersection )
override {
return impl().init( intersection ); }
444 virtual bool hasInterior ()
const override {
return impl().hasInterior(); }
447 virtual Linearization< RangeValueType >
linearizedInterior (
const InteriorCachingPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedInterior( asQP( x ), u ); }
448 virtual Linearization< RangeValueType >
linearizedInterior (
const InteriorElementPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedInterior( asQP( x ), u ); }
450 virtual bool hasBoundary ()
const override {
return impl().hasBoundary(); }
453 virtual Linearization< RangeValueType >
linearizedBoundary (
const SurfaceCachingPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedBoundary( asQP( x ), u ); }
454 virtual Linearization< RangeValueType >
linearizedBoundary (
const SurfaceElementPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedBoundary( asQP( x ), u ); }
456 virtual bool hasSkeleton ()
const override {
return impl().hasSkeleton(); }
457 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 ); }
458 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 ); }
459 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 ); }
460 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 ); }
467 const auto &impl ()
const {
return std::cref( impl_ ).get(); }
468 auto &impl () {
return std::ref( impl_ ).get(); }
473 template<
class Integrands >
474 using isVirtualized = std::is_same< std::decay_t< decltype( std::ref( std::declval< Integrands & >() ).
get() ) >, This >;
477 template<
class Integrands, std::enable_if_t< IntegrandsTraits< std::decay_t< Integrands > >::isFull && !isVirtualized< Integrands >::value,
int > = 0 >
479 : impl_( new Implementation< Integrands >( std::move( integrands ) ) )
482 template< class Integrands, std::enable_if_t< !IntegrandsTraits< Integrands >::isFull && !isVirtualized< Integrands >::value,
int > = 0 >
494 impl_.reset( other ? other.impl().clone() :
nullptr );
499 explicit operator bool ()
const {
return static_cast< bool >( impl_ ); }
506 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
509 return impl().interior( InteriorCachingPointType( x ), u );
512 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
515 return impl().interior( InteriorElementPointType( x ), u );
518 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
521 return impl().linearizedInterior( InteriorCachingPointType( x ), u );
524 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
527 return impl().linearizedInterior( InteriorElementPointType( x ), u );
532 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
535 return impl().boundary( SurfaceCachingPointType( x ), u );
538 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
541 return impl().boundary( SurfaceElementPointType( x ), u );
544 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
547 return impl().linearizedBoundary( SurfaceCachingPointType( x ), u );
550 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
553 return impl().linearizedBoundary( SurfaceElementPointType( x ), u );
558 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
561 return impl().skeleton( SurfaceCachingPointType( xIn ), uIn, SurfaceCachingPointType( xOut ), uOut );
564 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
567 return impl().skeleton( SurfaceElementPointType( xIn ), uIn, SurfaceElementPointType( xOut ), uOut );
570 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
573 return impl().linearizedSkeleton( SurfaceCachingPointType( xIn ), uIn, SurfaceCachingPointType( xOut ), uOut );
576 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
579 return impl().linearizedSkeleton( SurfaceElementPointType( xIn ), uIn, SurfaceElementPointType( xOut ), uOut );
584 return impl().hasDirichletBoundary();
588 return impl().isDirichletIntersection(inter,dirichletComponent);
592 return impl().dirichlet(bndId,x,value);
596 const Interface &impl ()
const { assert( impl_ );
return *impl_; }
597 Interface &impl () { assert( impl_ );
return *impl_; }
599 std::unique_ptr< Interface > impl_;
607 template<
class Model >
613 typedef typename GridPartType::template Codim< 0 >::EntityType
EntityType;
628 return (intersection.boundary() &&
model().hasNeumanBoundary() &&
model().
init( intersection.inside() ));
631 template<
class Po
int >
635 model().source( x, std::get< 0 >( u ), std::get< 1 >( u ), source );
637 model().diffusiveFlux( x, std::get< 0 >( u ), std::get< 1 >( u ), dFlux );
638 return std::make_tuple( source, dFlux );
641 template<
class Po
int >
646 model().linSource( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), source );
648 model().linDiffusiveFlux( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), dFlux );
649 return std::make_tuple( source, dFlux );
653 template<
class Po
int >
657 model().alpha( x, std::get< 0 >( u ), alpha );
658 return std::make_tuple( alpha, 0 );
661 template<
class Po
int >
666 model().linAlpha( std::get< 0 >( u ), x, std::get< 0 >( phi ), alpha );
667 return std::make_tuple( alpha, 0 );
671 const Model &
model ()
const {
return *model_; }
682 template<
class Model >
688 typedef typename GridPartType::template Codim< 0 >::EntityType
EntityType;
700 : model_( &
model ), penalty_( penalty )
705 intersection_ =
nullptr;
706 return model().init( entity );
711 intersection_ = &intersection;
712 if( intersection.boundary() )
714 const EntityType inside = intersection.inside();
715 beta_ = penalty_ * intersection.geometry().volume() / inside.geometry().volume();
718 else if( intersection.neighbor() )
720 const auto volIn = intersection.inside().geometry().volume();
721 const auto volOut = intersection.outside().geometry().volume();
722 beta_ = penalty_ * intersection.geometry().volume() /
std::min( volIn, volOut );
727 template<
class Po
int >
731 model().source( x, std::get< 0 >( u ), std::get< 1 >( u ), source );
733 model().diffusiveFlux( x, std::get< 0 >( u ), std::get< 1 >( u ), dFlux );
737 template<
class Po
int >
742 model().linSource( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), source );
744 model().linDiffusiveFlux( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), dFlux );
749 template<
class Po
int >
753 model().alpha( x, std::get< 0 >( u ), alpha );
757 template<
class Po
int >
762 model().linAlpha( std::get< 0 >( u ), x, std::get< 0 >( phi ), alpha );
767 template<
class Po
int >
770 const EntityType inside = intersection().inside();
771 const EntityType outside = intersection().outside();
774 const auto normal = intersection().unitOuterNormal( xIn.localPosition() );
777 std::get< 0 >( uJump ) = std::get< 0 >( uOut ) - std::get< 0 >( uIn );
778 for(
int i = 0; i < RangeType::dimension; ++i )
779 std::get< 1 >( uJump )[ i ].axpy( std::get< 0 >( uJump )[ i ], normal );
781 model().init( outside );
783 model().diffusiveFlux( xOut, std::get< 0 >( uOut ), std::get< 1 >( uJump ), dFluxPrimeOut );
784 model().diffusiveFlux( xOut, std::get< 0 >( uOut ), 0, dFluxOut );
785 dFluxPrimeOut -= dFluxOut;
787 model().diffusiveFlux( xOut, std::get< 0 >( uOut ), std::get< 1 >( uOut ), dFluxOut );
789 model().init( inside );
791 model().diffusiveFlux( xIn, std::get< 0 >( uIn ), std::get< 1 >( uJump ), dFluxPrimeIn );
792 model().diffusiveFlux( xIn, std::get< 0 >( uIn ), 0, dFluxIn );
793 dFluxPrimeIn -= dFluxIn;
795 model().diffusiveFlux( xIn, std::get< 0 >( uIn ), std::get< 1 >( uIn ), dFluxIn );
800 dFluxIn.usmv( -half, normal, int0 );
802 dFluxPrimeIn *= -half;
803 dFluxPrimeOut *= -half;
808 template<
class Po
int >
811 const auto normal = intersection().unitOuterNormal( xIn.localPosition() );
814 std::get< 0 >( uJump ) = std::get< 0 >( uOut ) - std::get< 0 >( uIn );
815 for(
int i = 0; i < RangeType::dimension; ++i )
816 std::get< 1 >( uJump )[ i ].axpy( std::get< 0 >( uJump )[ i ], normal );
818 auto intIn = [
this, xIn, uIn, xOut, uOut, normal, uJump ] (
const DomainValueType &phiIn ) {
819 const EntityType inside = intersection().inside();
820 const EntityType outside = intersection().outside();
825 std::get< 0 >( phiJump ) -= std::get< 0 >( phiIn );
826 for(
int i = 0; i < RangeType::dimension; ++i )
827 std::get< 1 >( phiJump )[ i ].axpy( std::get< 0 >( phiJump )[ i ], normal );
829 model().init( outside );
831 model().linDiffusiveFlux( std::get< 0 >( uOut ), std::get< 1 >( uJump ), xOut, 0, std::get< 1 >( phiJump ), dFluxPrimeOut );
833 model().init( inside );
835 model().linDiffusiveFlux( std::get< 0 >( uIn ), std::get< 1 >( uJump ), xIn, std::get< 0 >( phiIn ), std::get< 1 >( phiJump ), dFluxPrimeIn );
836 model().linDiffusiveFlux( std::get< 0 >( uIn ), 0, xIn, std::get< 0 >( phiIn ), 0, dFluxIn );
837 dFluxPrimeIn -= dFluxIn;
839 model().linDiffusiveFlux( std::get< 0 >( uIn ), std::get< 1 >( uIn ), xIn, std::get< 0 >( phiIn ), std::get< 1 >( phiIn ), dFluxIn );
841 RangeType int0 = std::get< 0 >( phiJump );
843 dFluxIn.usmv( -half, normal, int0 );
845 dFluxPrimeIn *= -half;
846 dFluxPrimeOut *= -half;
851 auto intOut = [
this, xIn, uIn, xOut, uOut, normal, uJump ] (
const DomainValueType &phiOut ) {
852 const EntityType inside = intersection().inside();
853 const EntityType outside = intersection().outside();
858 std::get< 0 >( phiJump ) = std::get< 0 >( phiOut );
859 for(
int i = 0; i < RangeType::dimension; ++i )
860 std::get< 1 >( phiJump )[ i ].axpy( std::get< 0 >( phiJump )[ i ], normal );
862 model().init( outside );
864 model().linDiffusiveFlux( std::get< 0 >( uOut ), std::get< 1 >( uJump ), xOut, std::get< 0 >( phiOut ), std::get< 1 >( phiJump ), dFluxPrimeOut );
865 model().linDiffusiveFlux( std::get< 0 >( uOut ), 0, xOut, std::get< 0 >( phiOut ), 0, dFluxOut );
866 dFluxPrimeOut -= dFluxOut;
868 model().linDiffusiveFlux( std::get< 0 >( uOut ), std::get< 1 >( uOut ), xOut, std::get< 0 >( phiOut ), std::get< 1 >( phiOut ), dFluxOut );
870 model().init( inside );
872 model().linDiffusiveFlux( std::get< 0 >( uIn ), std::get< 1 >( uJump ), xIn, 0, std::get< 1 >( phiJump ), dFluxPrimeIn );
874 RangeType int0 = std::get< 0 >( phiJump );
876 dFluxOut.usmv( -half, normal, int0 );
878 dFluxPrimeIn *= -half;
879 dFluxPrimeOut *= -half;
884 return std::make_pair( intIn, intOut );
887 const Model &
model ()
const {
return *model_; }
890 const IntersectionType &intersection ()
const { assert( intersection_ );
return *intersection_; }
Definition: bindguard.hh:11
std::tuple_element< i, Tuple >::type & get(Dune::TypeIndexedTuple< Tuple, Types > &tuple)
Definition: typeindexedtuple.hh:122
static constexpr T min(T a)
Definition: utility.hh:93
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
FullIntegrands(Args &&... args)
Definition: integrands.hh:264
bool init(const IntersectionType &intersection)
Definition: integrands.hh:269
bool hasInterior() const
Definition: integrands.hh:271
bool hasSkeleton() const
Definition: integrands.hh:299
auto linearizedBoundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:294
std::pair< RangeValueType, RangeValueType > skeleton(const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:302
bool hasBoundary() const
Definition: integrands.hh:285
IntegrandsTraits< Integrands >::DomainValueType DomainValueType
Definition: integrands.hh:145
IntegrandsTraits< Integrands >::RRangeType RRangeType
Definition: integrands.hh:320
bool isDirichletIntersection(const IntersectionType &inter, DirichletComponentType &dirichletComponent) const
Definition: integrands.hh:326
bool hasDirichletBoundary() const
Definition: integrands.hh:322
bool init(const EntityType &entity)
Definition: integrands.hh:268
void dirichlet(int bndId, const Point &x, RRangeType &value) const
Definition: integrands.hh:331
IntegrandsTraits< Integrands >::DirichletComponentType DirichletComponentType
Definition: integrands.hh:321
RangeValueType interior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:274
auto linearizedSkeleton(const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:308
IntegrandsTraits< Integrands >::EntityType EntityType
Definition: integrands.hh:149
auto linearizedInterior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:280
IntegrandsTraits< Integrands >::IntersectionType IntersectionType
Definition: integrands.hh:150
RangeValueType boundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:288
IntegrandsTraits< Integrands >::GridPartType GridPartType
Definition: integrands.hh:147
Integrands integrands_
Definition: integrands.hh:317
IntegrandsTraits< Integrands >::RangeValueType RangeValueType
Definition: integrands.hh:146
Definition: integrands.hh:368
RangeValueType boundary(const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u) const
Definition: integrands.hh:533
DomainValue DomainValueType
Definition: integrands.hh:373
VirtualizedIntegrands & operator=(const This &other)
Definition: integrands.hh:492
RangeValueType interior(const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u) const
Definition: integrands.hh:507
auto linearizedBoundary(const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u) const
Definition: integrands.hh:545
GridPart GridPartType
Definition: integrands.hh:372
VirtualizedIntegrands(Integrands integrands)
Definition: integrands.hh:478
VirtualizedIntegrands(This &&)=default
std::array< int, RRangeType::dimension > DirichletComponentType
Definition: integrands.hh:380
typename detail::GetDimRange< std::tuple_element_t< 0, RangeValueType > >::type RRangeType
Definition: integrands.hh:379
VirtualizedIntegrands(const This &other)
Definition: integrands.hh:487
bool hasSkeleton() const
Definition: integrands.hh:556
auto linearizedSkeleton(const Fem::QuadraturePointWrapper< Quadrature > &xIn, const DomainValueType &uIn, const Fem::QuadraturePointWrapper< Quadrature > &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:571
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: integrands.hh:376
EntityType::Geometry::LocalCoordinate DomainType
Definition: integrands.hh:381
bool init(const IntersectionType &intersection)
Definition: integrands.hh:502
bool isDirichletIntersection(const IntersectionType &inter, DirichletComponentType &dirichletComponent) const
Definition: integrands.hh:586
RangeValue RangeValueType
Definition: integrands.hh:374
bool hasInterior() const
Definition: integrands.hh:504
GridPartType::IntersectionType IntersectionType
Definition: integrands.hh:377
bool hasBoundary() const
Definition: integrands.hh:530
auto linearizedInterior(const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u) const
Definition: integrands.hh:519
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:559
bool hasDirichletBoundary() const
Definition: integrands.hh:582
bool init(const EntityType &entity)
Definition: integrands.hh:501
void dirichlet(int bndId, const DomainType &x, RRangeType &value) const
Definition: integrands.hh:590
Definition: integrands.hh:609
GridPartType::IntersectionType IntersectionType
Definition: integrands.hh:614
Model::RangeType RangeType
Definition: integrands.hh:616
Model ModelType
Definition: integrands.hh:610
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: integrands.hh:613
std::tuple< RangeType, JacobianRangeType > RangeValueType
Definition: integrands.hh:620
Model::JacobianRangeType JacobianRangeType
Definition: integrands.hh:617
auto linearizedBoundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:662
const Model & model() const
Definition: integrands.hh:671
bool init(const IntersectionType &intersection)
Definition: integrands.hh:626
std::tuple< RangeType, JacobianRangeType > DomainValueType
Definition: integrands.hh:619
bool init(const EntityType &entity)
Definition: integrands.hh:624
DiffusionModelIntegrands(const Model &model)
Definition: integrands.hh:622
RangeValueType boundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:654
RangeValueType interior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:632
auto linearizedInterior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:642
Model::GridPartType GridPartType
Definition: integrands.hh:611
Definition: integrands.hh:684
std::pair< RangeValueType, RangeValueType > skeleton(const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:768
RangeValueType boundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:750
bool init(const IntersectionType &intersection)
Definition: integrands.hh:709
auto linearizedInterior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:738
Model::GridPartType GridPartType
Definition: integrands.hh:686
FieldTraits< RangeType >::field_type RangeFieldType
Definition: integrands.hh:694
RangeValueType interior(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:728
std::tuple< RangeType, JacobianRangeType > RangeValueType
Definition: integrands.hh:697
Model::JacobianRangeType JacobianRangeType
Definition: integrands.hh:692
Model ModelType
Definition: integrands.hh:685
bool init(const EntityType &entity)
Definition: integrands.hh:703
GridPartType::IntersectionType IntersectionType
Definition: integrands.hh:689
DGDiffusionModelIntegrands(const Model &model, RangeFieldType penalty)
Definition: integrands.hh:699
std::tuple< RangeType, JacobianRangeType > DomainValueType
Definition: integrands.hh:696
auto linearizedBoundary(const Point &x, const DomainValueType &u) const
Definition: integrands.hh:758
const Model & model() const
Definition: integrands.hh:887
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: integrands.hh:688
auto linearizedSkeleton(const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut) const
Definition: integrands.hh:809
Model::RangeType RangeType
Definition: integrands.hh:691