1#ifndef DUNE_FEM_SCHEMES_INTEGRANDS_HH
2#define DUNE_FEM_SCHEMES_INTEGRANDS_HH
14#include <dune/fempy/quadrature/cachingpoint.hh>
15#include <dune/fempy/quadrature/elementpoint.hh>
16#include <dune/fem/quadrature/cachingquadrature.hh>
17#include <dune/fem/common/explicitfieldvector.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 >
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 );
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 ( ... );
75 template< class Integrands, std::enable_if_t< std::is_same< decltype( std::declval< const Integrands & >().nonlinear() ),
bool >::value,
int > = 0 >
76 std::true_type hasNonLinear (
const Integrands & );
78 std::false_type hasNonLinear ( ... );
81 template<
class Integrands >
82 std::true_type boundary (
const Integrands &,
decltype( std::declval< const Integrands & >().boundary( std::declval<
const SurfaceQuadraturePointType< Integrands > & >(), std::declval<
const DomainValueType< Integrands > & >() ) ) * =
nullptr );
84 std::false_type boundary ( ... );
86 template< class Integrands, std::enable_if_t< std::is_same< decltype( std::declval< const Integrands & >().hasBoundary() ),
bool >::value,
int > = 0 >
87 std::true_type hasBoundary (
const Integrands & );
89 std::false_type hasBoundary ( ... );
92 template<
class Integrands >
93 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 );
95 std::false_type skeleton ( ... );
97 template< class Integrands, std::enable_if_t< std::is_same< decltype( std::declval< const Integrands & >().hasSkeleton() ),
bool >::value,
int > = 0 >
98 std::true_type hasSkeleton (
const Integrands & );
100 std::false_type hasSkeleton ( ... );
102 template<
class Integrands >
103 using Get = std::decay_t< decltype( std::ref( std::declval< const Integrands & >() ).
get() ) >;
105 template<
class Integrands >
106 using RRangeType =
typename Get<Integrands>::RRangeType;
107 template<
class Integrands >
108 using RJacobianRangeType =
typename Get<Integrands>::RJacobianRangeType;
109 template<
class Integrands >
110 using DirichletComponentType =
typename Get<Integrands>::DirichletComponentType;
115 template<
class Integrands >
116 struct IntegrandsTraits
118 typedef Impl::IntegrandsTraits::DomainValueType< Integrands > DomainValueType;
119 typedef Impl::IntegrandsTraits::RangeValueType< Integrands > RangeValueType;
120 typedef Impl::IntegrandsTraits::GridPartType< Integrands > GridPartType;
122 typedef Impl::IntegrandsTraits::EntityType< Integrands > EntityType;
123 typedef Impl::IntegrandsTraits::IntersectionType< Integrands > IntersectionType;
125 static const bool hasNonLinear =
decltype( Impl::IntegrandsTraits::hasNonLinear( std::declval< const Integrands & >() ) )::value;
127 static const bool interior =
decltype( Impl::IntegrandsTraits::interior( std::declval< const Integrands & >() ) )::value;
128 static const bool hasInterior =
decltype( Impl::IntegrandsTraits::hasInterior( std::declval< const Integrands & >() ) )::value;
130 static const bool boundary =
decltype( Impl::IntegrandsTraits::boundary( std::declval< const Integrands & >() ) )::value;
131 static const bool hasBoundary =
decltype( Impl::IntegrandsTraits::hasBoundary( std::declval< const Integrands & >() ) )::value;
133 static const bool skeleton =
decltype( Impl::IntegrandsTraits::skeleton( std::declval< const Integrands & >() ) )::value;
134 static const bool hasSkeleton =
decltype( Impl::IntegrandsTraits::hasSkeleton( std::declval< const Integrands & >() ) )::value;
136 static_assert( (!hasInterior ||
interior),
"Existence of method 'hasInterior' implies existence of method interior." );
137 static_assert( (!hasBoundary || boundary),
"Existence of method 'hasBoundary' implies existence of method boundary." );
138 static_assert( (!hasSkeleton || skeleton),
"Existence of method 'hasSkeleton' implies existence of method skeleton." );
140 static const bool isFull = hasInterior && hasBoundary && hasSkeleton;
142 typedef Impl::IntegrandsTraits::RRangeType< Integrands > RRangeType;
143 typedef Impl::IntegrandsTraits::RJacobianRangeType< Integrands > RJacobianRangeType;
144 typedef Impl::IntegrandsTraits::DirichletComponentType< Integrands > DirichletComponentType;
152 template<
class Integrands >
153 struct FullIntegrands
155 typedef typename IntegrandsTraits< Integrands >::DomainValueType DomainValueType;
156 typedef typename IntegrandsTraits< Integrands >::RangeValueType RangeValueType;
157 typedef typename IntegrandsTraits< Integrands >::GridPartType GridPartType;
159 typedef typename IntegrandsTraits< Integrands >::EntityType EntityType;
160 typedef typename IntegrandsTraits< Integrands >::IntersectionType IntersectionType;
163 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasNonLinear,
int > = 0 >
164 static bool nonlinear (
const T &integrands )
166 return integrands.nonlinear();
169 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasNonLinear,
int > = 0 >
170 static bool nonlinear (
const T &integrands )
176 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasInterior,
int > = 0 >
177 static bool hasInterior (
const T &integrands )
179 return integrands.hasInterior();
182 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasInterior,
int > = 0 >
183 static bool hasInterior (
const T &integrands )
185 return IntegrandsTraits< T >::interior;
188 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::interior,
int > = 0 >
189 static RangeValueType
interior (
const T &integrands,
const Point &x,
const DomainValueType &u )
191 return integrands.interior( x, u );
194 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::interior,
int > = 0 >
195 static RangeValueType
interior (
const T &integrands,
const Point &x,
const DomainValueType &u )
197 return RangeValueType();
200 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::interior,
int > = 0 >
201 static auto linearizedInterior (
const T &integrands,
const Point &x,
const DomainValueType &u )
203 return integrands.linearizedInterior( x, u );
206 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::interior,
int > = 0 >
207 static auto linearizedInterior (
const T &integrands,
const Point &x,
const DomainValueType &u )
209 return [] (
const DomainValueType & ) {
return RangeValueType(); };
212 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasBoundary,
int > = 0 >
213 static bool hasBoundary (
const T &integrands )
215 return integrands.hasBoundary();
218 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasBoundary,
int > = 0 >
219 static bool hasBoundary (
const T &integrands )
221 return IntegrandsTraits< T >::boundary;
224 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::boundary,
int > = 0 >
225 static RangeValueType boundary (
const T &integrands,
const Point &x,
const DomainValueType &u )
227 return integrands.boundary( x, u );
230 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::boundary,
int > = 0 >
231 static RangeValueType boundary (
const T &integrands,
const Point &x,
const DomainValueType &u )
233 return RangeValueType();
236 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::boundary,
int > = 0 >
237 static auto linearizedBoundary (
const T &integrands,
const Point &x,
const DomainValueType &u )
239 return integrands.linearizedBoundary( x, u );
242 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::boundary,
int > = 0 >
243 static auto linearizedBoundary (
const T &integrands,
const Point &x,
const DomainValueType &u )
245 return [] (
const DomainValueType & ) {
return RangeValueType(); };
248 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasSkeleton,
int > = 0 >
249 static bool hasSkeleton (
const T &integrands )
251 return integrands.hasSkeleton();
254 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasSkeleton,
int > = 0 >
255 static bool hasSkeleton (
const T &integrands )
257 return IntegrandsTraits< T >::skeleton;
260 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::skeleton,
int > = 0 >
261 static std::pair< RangeValueType, RangeValueType > skeleton (
const T &integrands,
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
263 return integrands.skeleton( xIn, uIn, xOut, uOut );
266 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::skeleton,
int > = 0 >
267 static std::pair< RangeValueType, RangeValueType > skeleton (
const T &integrands,
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
269 return std::make_pair( RangeValueType(), RangeValueType() );
272 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::skeleton,
int > = 0 >
273 static auto linearizedSkeleton (
const T &integrands,
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
275 return integrands.linearizedSkeleton( xIn, uIn, xOut, uOut );
278 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::skeleton,
int > = 0 >
279 static auto linearizedSkeleton (
const T &integrands,
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
281 auto zero = [] (
const DomainValueType & ) {
return std::make_pair( RangeValueType(), RangeValueType() ); };
282 return std::make_pair( zero, zero );
286 template<
class... Args >
287 explicit FullIntegrands ( Args &&... args )
288 : integrands_(
std::forward< Args >( args )... ),
289 rInt_(
std::ref( integrands_ ).
get() )
293 bool init (
const EntityType &entity ) {
return integrands().init( entity ); }
294 bool init (
const IntersectionType &intersection ) {
return integrands().init( intersection ); }
295 void unbind ( ) { integrands().unbind( ); }
297 bool nonlinear()
const {
return nonlinear( integrands() ); }
299 bool hasInterior ()
const {
return hasInterior( integrands() ); }
301 template<
class Po
int >
302 RangeValueType
interior (
const Point &x,
const DomainValueType &u )
const
304 return interior( integrands(), x, u );
307 template<
class Po
int >
308 auto linearizedInterior (
const Point &x,
const DomainValueType &u )
const
310 return linearizedInterior( integrands(), x, u );
313 bool hasBoundary ()
const {
return hasBoundary( integrands() ); }
315 template<
class Po
int >
316 RangeValueType boundary (
const Point &x,
const DomainValueType &u )
const
318 return boundary( integrands(), x, u );
321 template<
class Po
int >
322 auto linearizedBoundary (
const Point &x,
const DomainValueType &u )
const
324 return linearizedBoundary( integrands(), x, u );
327 bool hasSkeleton ()
const {
return hasSkeleton( integrands() ); }
329 template<
class Po
int >
330 std::pair< RangeValueType, RangeValueType > skeleton (
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
const
332 return skeleton( integrands(), xIn, uIn, xOut, uOut );
335 template<
class Po
int >
336 auto linearizedSkeleton (
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
const
338 return linearizedSkeleton( integrands(), xIn, uIn, xOut, uOut );
342 typedef typename Integrands::type RealIntegrands;
345 RealIntegrands& integrands () {
return rInt_; }
346 const RealIntegrands& integrands ()
const {
return rInt_; }
350 Integrands integrands_;
351 RealIntegrands rInt_;
354 typedef typename IntegrandsTraits< Integrands >::RRangeType RRangeType;
355 typedef typename IntegrandsTraits< Integrands >::RJacobianRangeType RJacobianRangeType;
356 typedef typename IntegrandsTraits< Integrands >::DirichletComponentType DirichletComponentType;
357 bool hasDirichletBoundary ()
const
359 return integrands().hasDirichletBoundary();
361 bool isDirichletIntersection(
const IntersectionType& inter, DirichletComponentType &dirichletComponent )
const
363 return integrands().isDirichletIntersection(inter,dirichletComponent);
365 template <
class Po
int>
366 void dirichlet(
int bndId,
const Point &x, RRangeType &value)
const
368 return integrands().dirichlet(bndId,x,value);
370 template <
class Po
int>
371 void dDirichlet(
int bndId,
const Point &x, RJacobianRangeType &value)
const
373 return integrands().dDirichlet(bndId,x,value);
386 template <
class FT,
int r>
387 struct GetDimRange<
Dune::FieldVector<FT,r>>
392 static const int value = r;
394 template <
class FT,
int r,
int c>
395 struct GetDimRange<
Dune::FieldMatrix<FT,r,c>>
400 static const int value = r;
402 template <
class FT,
int r,
int c>
403 struct GetDimRange<
Dune::Fem::ExplicitFieldVector<Dune::FieldMatrix<FT,c,c>,r>>
408 static const int value = r;
412 template<
class Gr
idPart,
class DomainValue,
class RangeValue = DomainValue >
413 class VirtualizedIntegrands
415 typedef VirtualizedIntegrands< GridPart, DomainValue, RangeValue > This;
418 typedef GridPart GridPartType;
419 typedef DomainValue DomainValueType;
420 typedef RangeValue RangeValueType;
422 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
425 using RRangeType =
typename detail::GetDimRange<std::tuple_element_t<0,RangeValueType>>::type;
426 using RJacobianRangeType =
typename detail::GetDimRange<std::tuple_element_t<0,RangeValueType>>::template dtype<GridPart::dimension>;
427 typedef std::array<int,RRangeType::dimension> DirichletComponentType;
428 typedef typename EntityType::Geometry::LocalCoordinate DomainType;
431 typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
433 typedef FemPy::CachingPoint< GridPart, LocalCoordinateType, 0 > InteriorCachingPointType;
434 typedef FemPy::ElementPoint< GridPart, LocalCoordinateType, 0 > InteriorElementPointType;
435 typedef FemPy::CachingPoint< GridPart, LocalCoordinateType, 1 > SurfaceCachingPointType;
436 typedef FemPy::ElementPoint< GridPart, LocalCoordinateType, 1 > SurfaceElementPointType;
439 static Fem::QuadraturePointWrapper< QP > asQP (
const QP &qp )
441 return static_cast< Fem::QuadraturePointWrapper< QP >
>( qp );
445 using Linearization = std::function< R(
const DomainValueType & ) >;
448 using LinearizationPair = std::pair< Linearization< std::pair< R, R > >, Linearization< std::pair< R, R > > >;
452 virtual ~Interface () {}
453 virtual Interface *clone ()
const = 0;
455 virtual bool init (
const EntityType &entity ) = 0;
456 virtual bool init (
const IntersectionType &intersection ) = 0;
457 virtual void unbind ( ) = 0;
459 virtual bool nonlinear()
const = 0;
461 virtual bool hasInterior ()
const = 0;
462 virtual RangeValueType
interior (
const InteriorCachingPointType &x,
const DomainValueType &u )
const = 0;
463 virtual RangeValueType
interior (
const InteriorElementPointType &x,
const DomainValueType &u )
const = 0;
464 virtual Linearization< RangeValueType > linearizedInterior (
const InteriorCachingPointType &x,
const DomainValueType &u )
const = 0;
465 virtual Linearization< RangeValueType > linearizedInterior (
const InteriorElementPointType &x,
const DomainValueType &u )
const = 0;
467 virtual bool hasBoundary ()
const = 0;
468 virtual RangeValueType boundary (
const SurfaceCachingPointType &x,
const DomainValueType &u )
const = 0;
469 virtual RangeValueType boundary (
const SurfaceElementPointType &x,
const DomainValueType &u )
const = 0;
470 virtual Linearization< RangeValueType > linearizedBoundary (
const SurfaceCachingPointType &x,
const DomainValueType &u )
const = 0;
471 virtual Linearization< RangeValueType > linearizedBoundary (
const SurfaceElementPointType &x,
const DomainValueType &u )
const = 0;
473 virtual bool hasSkeleton ()
const = 0;
474 virtual std::pair< RangeValueType, RangeValueType > skeleton (
const SurfaceCachingPointType &xIn,
const DomainValueType &uIn,
const SurfaceCachingPointType &xOut,
const DomainValueType &uOut )
const = 0;
475 virtual std::pair< RangeValueType, RangeValueType > skeleton (
const SurfaceElementPointType &xIn,
const DomainValueType &uIn,
const SurfaceElementPointType &xOut,
const DomainValueType &uOut )
const = 0;
476 virtual LinearizationPair< RangeValueType > linearizedSkeleton (
const SurfaceCachingPointType &xIn,
const DomainValueType &uIn,
const SurfaceCachingPointType &xOut,
const DomainValueType &uOut )
const = 0;
477 virtual LinearizationPair< RangeValueType > linearizedSkeleton (
const SurfaceElementPointType &xIn,
const DomainValueType &uIn,
const SurfaceElementPointType &xOut,
const DomainValueType &uOut )
const = 0;
479 virtual bool hasDirichletBoundary ()
const = 0;
480 virtual bool isDirichletIntersection(
const IntersectionType& inter, DirichletComponentType &dirichletComponent )
const = 0;
481 virtual void dirichlet(
int bndId,
const DomainType &x,RRangeType &value)
const = 0;
482 virtual void dDirichlet(
int bndId,
const DomainType &x,RJacobianRangeType &value)
const = 0;
485 template<
class Impl >
489 Implementation ( Impl impl ) : impl_(
std::move( impl ) )
492 virtual Interface *clone ()
const override {
return new Implementation( *
this ); }
494 virtual bool init (
const EntityType &entity )
override {
return impl().init( entity ); }
495 virtual bool init (
const IntersectionType &intersection )
override {
return impl().init( intersection ); }
496 virtual void unbind ( )
override { impl().unbind( ); }
498 virtual bool nonlinear ()
const override {
return impl().nonlinear(); }
500 virtual bool hasInterior ()
const override {
return impl().hasInterior(); }
501 virtual RangeValueType
interior (
const InteriorCachingPointType &x,
const DomainValueType &u )
const override {
return impl().interior( asQP( x ), u ); }
502 virtual RangeValueType
interior (
const InteriorElementPointType &x,
const DomainValueType &u )
const override {
return impl().interior( asQP( x ), u ); }
503 virtual Linearization< RangeValueType > linearizedInterior (
const InteriorCachingPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedInterior( asQP( x ), u ); }
504 virtual Linearization< RangeValueType > linearizedInterior (
const InteriorElementPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedInterior( asQP( x ), u ); }
506 virtual bool hasBoundary ()
const override {
return impl().hasBoundary(); }
507 virtual RangeValueType boundary (
const SurfaceCachingPointType &x,
const DomainValueType &u )
const override {
return impl().boundary( asQP( x ), u ); }
508 virtual RangeValueType boundary (
const SurfaceElementPointType &x,
const DomainValueType &u )
const override {
return impl().boundary( asQP( x ), u ); }
509 virtual Linearization< RangeValueType > linearizedBoundary (
const SurfaceCachingPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedBoundary( asQP( x ), u ); }
510 virtual Linearization< RangeValueType > linearizedBoundary (
const SurfaceElementPointType &x,
const DomainValueType &u )
const override {
return impl().linearizedBoundary( asQP( x ), u ); }
512 virtual bool hasSkeleton ()
const override {
return impl().hasSkeleton(); }
513 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 ); }
514 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 ); }
515 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 ); }
516 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 ); }
518 virtual bool hasDirichletBoundary ()
const override {
return impl().hasDirichletBoundary(); }
519 virtual bool isDirichletIntersection(
const IntersectionType& inter, DirichletComponentType &dirichletComponent )
const override {
return impl().isDirichletIntersection(inter,dirichletComponent); }
520 virtual void dirichlet(
int bndId,
const DomainType &x,RRangeType &value)
const override { impl().dirichlet(bndId,x,value); }
521 virtual void dDirichlet(
int bndId,
const DomainType &x,RJacobianRangeType &value)
const override { impl().dDirichlet(bndId,x,value); }
524 const auto &impl ()
const {
return std::cref( impl_ ).get(); }
525 auto &impl () {
return std::ref( impl_ ).get(); }
530 template<
class Integrands >
531 using isVirtualized = std::is_same< std::decay_t< decltype( std::ref( std::declval< Integrands & >() ).
get() ) >, This >;
534 template<
class Integrands, std::enable_if_t< IntegrandsTraits< std::decay_t< Integrands > >::isFull && !isVirtualized< Integrands >::value,
int > = 0 >
535 explicit VirtualizedIntegrands ( Integrands integrands )
536 : impl_( new Implementation< Integrands >(
std::move( integrands ) ) )
539 template< class Integrands, std::enable_if_t< !IntegrandsTraits< Integrands >::isFull && !isVirtualized< Integrands >::value,
int > = 0 >
540 explicit VirtualizedIntegrands ( Integrands integrands )
541 : VirtualizedIntegrands( FullIntegrands<
std::decay_t< Integrands > >(
std::move( integrands ) ) )
544 VirtualizedIntegrands (
const This &other ) : impl_( other ? other.impl().clone() : nullptr )
547 VirtualizedIntegrands ( This && ) =
default;
549 VirtualizedIntegrands &operator= (
const This &other )
551 impl_.reset( other ? other.impl().clone() :
nullptr );
554 VirtualizedIntegrands &operator= ( This && ) =
default;
556 explicit operator bool ()
const {
return static_cast< bool >( impl_ ); }
558 bool init (
const EntityType &entity ) {
return impl().init( entity ); }
559 bool init (
const IntersectionType &intersection ) {
return impl().init( intersection ); }
560 void unbind ( ) { impl().unbind( ); }
562 bool nonlinear()
const {
return impl().nonlinear(); }
564 bool hasInterior ()
const {
return impl().hasInterior(); }
566 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
567 RangeValueType
interior (
const Fem::QuadraturePointWrapper< Quadrature > &x,
const DomainValueType &u )
const
569 return impl().interior( InteriorCachingPointType( x ), u );
572 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
573 RangeValueType
interior (
const Fem::QuadraturePointWrapper< Quadrature > &x,
const DomainValueType &u )
const
575 return impl().interior( InteriorElementPointType( x ), u );
578 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
579 auto linearizedInterior (
const Fem::QuadraturePointWrapper< Quadrature > &x,
const DomainValueType &u )
const
581 return impl().linearizedInterior( InteriorCachingPointType( x ), u );
584 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
585 auto linearizedInterior (
const Fem::QuadraturePointWrapper< Quadrature > &x,
const DomainValueType &u )
const
587 return impl().linearizedInterior( InteriorElementPointType( x ), u );
590 bool hasBoundary ()
const {
return impl().hasBoundary(); }
592 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
593 RangeValueType boundary (
const Fem::QuadraturePointWrapper< Quadrature > &x,
const DomainValueType &u )
const
595 return impl().boundary( SurfaceCachingPointType( x ), u );
598 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
599 RangeValueType boundary (
const Fem::QuadraturePointWrapper< Quadrature > &x,
const DomainValueType &u )
const
601 return impl().boundary( SurfaceElementPointType( x ), u );
604 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
605 auto linearizedBoundary (
const Fem::QuadraturePointWrapper< Quadrature > &x,
const DomainValueType &u )
const
607 return impl().linearizedBoundary( SurfaceCachingPointType( x ), u );
610 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
611 auto linearizedBoundary (
const Fem::QuadraturePointWrapper< Quadrature > &x,
const DomainValueType &u )
const
613 return impl().linearizedBoundary( SurfaceElementPointType( x ), u );
616 bool hasSkeleton ()
const {
return impl().hasSkeleton(); }
618 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
619 std::pair< RangeValueType, RangeValueType > skeleton (
const Fem::QuadraturePointWrapper< Quadrature > &xIn,
const DomainValueType &uIn,
const Fem::QuadraturePointWrapper< Quadrature > &xOut,
const DomainValueType &uOut )
const
621 return impl().skeleton( SurfaceCachingPointType( xIn ), uIn, SurfaceCachingPointType( xOut ), uOut );
624 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
625 std::pair< RangeValueType, RangeValueType > skeleton (
const Fem::QuadraturePointWrapper< Quadrature > &xIn,
const DomainValueType &uIn,
const Fem::QuadraturePointWrapper< Quadrature > &xOut,
const DomainValueType &uOut )
const
627 return impl().skeleton( SurfaceElementPointType( xIn ), uIn, SurfaceElementPointType( xOut ), uOut );
630 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
631 auto linearizedSkeleton (
const Fem::QuadraturePointWrapper< Quadrature > &xIn,
const DomainValueType &uIn,
const Fem::QuadraturePointWrapper< Quadrature > &xOut,
const DomainValueType &uOut )
const
633 return impl().linearizedSkeleton( SurfaceCachingPointType( xIn ), uIn, SurfaceCachingPointType( xOut ), uOut );
636 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value,
int > = 0 >
637 auto linearizedSkeleton (
const Fem::QuadraturePointWrapper< Quadrature > &xIn,
const DomainValueType &uIn,
const Fem::QuadraturePointWrapper< Quadrature > &xOut,
const DomainValueType &uOut )
const
639 return impl().linearizedSkeleton( SurfaceElementPointType( xIn ), uIn, SurfaceElementPointType( xOut ), uOut );
642 bool hasDirichletBoundary ()
const
644 return impl().hasDirichletBoundary();
646 bool isDirichletIntersection(
const IntersectionType& inter, DirichletComponentType &dirichletComponent )
const
648 return impl().isDirichletIntersection(inter,dirichletComponent);
650 void dirichlet(
int bndId,
const DomainType &x,RRangeType &value)
const
652 return impl().dirichlet(bndId,x,value);
654 void dDirichlet(
int bndId,
const DomainType &x,RJacobianRangeType &value)
const
656 return impl().dDirichlet(bndId,x,value);
660 const Interface &impl ()
const { assert( impl_ );
return *impl_; }
661 Interface &impl () { assert( impl_ );
return *impl_; }
663 std::unique_ptr< Interface > impl_;
671 template<
class Model >
672 struct ConservationLawModelIntegrands
674 typedef Model ModelType;
675 typedef typename Model::GridPartType GridPartType;
677 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
680 typedef typename Model::RangeType RangeType;
681 typedef typename Model::JacobianRangeType JacobianRangeType;
683 typedef std::tuple< RangeType, JacobianRangeType > DomainValueType;
684 typedef std::tuple< RangeType, JacobianRangeType > RangeValueType;
686 explicit ConservationLawModelIntegrands (
const Model &model ) : model_( &model ) {}
688 bool init (
const EntityType &entity ) {
return model().init( entity ); }
690 bool init (
const IntersectionType &intersection )
692 return (intersection.boundary() && model().hasNeumanBoundary() && model().init( intersection.inside() ));
695 void unbind ( ) { model().unbind( ); }
697 template<
class Po
int >
698 RangeValueType
interior (
const Point &x,
const DomainValueType &u )
const
700 RangeType source( 0 );
701 model().source( x, std::get< 0 >( u ), std::get< 1 >( u ), source );
702 JacobianRangeType dFlux( 0 );
703 model().flux( x, std::get< 0 >( u ), std::get< 1 >( u ), dFlux );
704 return std::make_tuple( source, dFlux );
707 template<
class Po
int >
708 auto linearizedInterior (
const Point &x,
const DomainValueType &u )
const
710 return [
this, x, u ] (
const DomainValueType &phi ) {
711 RangeType source( 0 );
712 model().linSource( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), source );
713 JacobianRangeType dFlux( 0 );
714 model().linFlux( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), dFlux );
715 return std::make_tuple( source, dFlux );
719 template<
class Po
int >
720 RangeValueType boundary (
const Point &x,
const DomainValueType &u )
const
722 RangeType alpha( 0 );
723 model().alpha( x, std::get< 0 >( u ), alpha );
724 return std::make_tuple( alpha, 0 );
727 template<
class Po
int >
728 auto linearizedBoundary (
const Point &x,
const DomainValueType &u )
const
730 return [
this, x, u ] (
const DomainValueType &phi ) {
731 RangeType alpha( 0 );
732 model().linAlpha( std::get< 0 >( u ), x, std::get< 0 >( phi ), alpha );
733 return std::make_tuple( alpha, 0 );
737 const Model &model ()
const {
return *model_; }
740 const Model *model_ =
nullptr;
748 template<
class Model >
749 struct DGConservationLawModelIntegrands
751 typedef Model ModelType;
752 typedef typename Model::GridPartType GridPartType;
754 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
757 typedef typename Model::RangeType RangeType;
758 typedef typename Model::JacobianRangeType JacobianRangeType;
760 typedef typename FieldTraits< RangeType >::field_type RangeFieldType;
762 typedef std::tuple< RangeType, JacobianRangeType > DomainValueType;
763 typedef std::tuple< RangeType, JacobianRangeType > RangeValueType;
765 DGConservationLawModelIntegrands (
const Model &model, RangeFieldType penalty )
766 : model_( &model ), penalty_( penalty )
769 bool init (
const EntityType &entity )
771 intersection_ =
nullptr;
772 return model().init( entity );
775 bool init (
const IntersectionType &intersection )
777 intersection_ = &intersection;
778 if( intersection.boundary() )
780 const EntityType inside = intersection.inside();
781 beta_ = penalty_ * intersection.geometry().volume() / inside.geometry().volume();
782 return (model().hasNeumanBoundary() && model().init( inside ));
784 else if( intersection.neighbor() )
786 const auto volIn = intersection.inside().geometry().volume();
787 const auto volOut = intersection.outside().geometry().volume();
788 beta_ = penalty_ * intersection.geometry().volume() / std::min( volIn, volOut );
798 template<
class Po
int >
799 RangeValueType
interior (
const Point &x,
const DomainValueType &u )
const
801 RangeType source( 0 );
802 model().source( x, std::get< 0 >( u ), std::get< 1 >( u ), source );
803 JacobianRangeType dFlux( 0 );
804 model().flux( x, std::get< 0 >( u ), std::get< 1 >( u ), dFlux );
805 return RangeValueType( source, dFlux );
808 template<
class Po
int >
809 auto linearizedInterior (
const Point &x,
const DomainValueType &u )
const
811 return [
this, x, u ] (
const DomainValueType &phi ) {
812 RangeType source( 0 );
813 model().linSource( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), source );
814 JacobianRangeType dFlux( 0 );
815 model().linFlux( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), dFlux );
816 return RangeValueType( source, dFlux );
820 template<
class Po
int >
821 RangeValueType boundary (
const Point &x,
const DomainValueType &u )
const
823 RangeType alpha( 0 );
824 model().alpha( x, std::get< 0 >( u ), alpha );
825 return RangeValueType( alpha, 0 );
828 template<
class Po
int >
829 auto linearizedBoundary (
const Point &x,
const DomainValueType &u )
const
831 return [
this, x, u ] (
const DomainValueType &phi ) {
832 RangeType alpha( 0 );
833 model().linAlpha( std::get< 0 >( u ), x, std::get< 0 >( phi ), alpha );
834 return RangeValueType( alpha, 0 );
838 template<
class Po
int >
839 std::pair< RangeValueType, RangeValueType > skeleton (
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
const
841 const EntityType inside = intersection().inside();
842 const EntityType outside = intersection().outside();
844 const RangeFieldType half = RangeFieldType( 1 ) / RangeFieldType( 2 );
845 const auto normal = intersection().unitOuterNormal( xIn.localPosition() );
847 DomainValueType uJump( 0, 0 );
848 std::get< 0 >( uJump ) = std::get< 0 >( uOut ) - std::get< 0 >( uIn );
849 for(
int i = 0; i < RangeType::dimension; ++i )
850 std::get< 1 >( uJump )[ i ].axpy( std::get< 0 >( uJump )[ i ], normal );
852 model().init( outside );
853 JacobianRangeType dFluxOut( 0 ), dFluxPrimeOut( 0 );
854 model().flux( xOut, std::get< 0 >( uOut ), std::get< 1 >( uJump ), dFluxPrimeOut );
855 model().flux( xOut, std::get< 0 >( uOut ), 0, dFluxOut );
856 dFluxPrimeOut -= dFluxOut;
858 model().flux( xOut, std::get< 0 >( uOut ), std::get< 1 >( uOut ), dFluxOut );
860 model().init( inside );
861 JacobianRangeType dFluxIn( 0 ), dFluxPrimeIn( 0 );
862 model().flux( xIn, std::get< 0 >( uIn ), std::get< 1 >( uJump ), dFluxPrimeIn );
863 model().flux( xIn, std::get< 0 >( uIn ), 0, dFluxIn );
864 dFluxPrimeIn -= dFluxIn;
866 model().flux( xIn, std::get< 0 >( uIn ), std::get< 1 >( uIn ), dFluxIn );
868 RangeType int0 = std::get< 0 >( uJump );
871 dFluxIn.usmv( -half, normal, int0 );
873 dFluxPrimeIn *= -half;
874 dFluxPrimeOut *= -half;
876 return std::make_pair( RangeValueType( -int0, dFluxPrimeIn ), RangeValueType( int0, dFluxPrimeOut ) );
879 template<
class Po
int >
880 auto linearizedSkeleton (
const Point &xIn,
const DomainValueType &uIn,
const Point &xOut,
const DomainValueType &uOut )
const
882 const auto normal = intersection().unitOuterNormal( xIn.localPosition() );
884 DomainValueType uJump( 0, 0 );
885 std::get< 0 >( uJump ) = std::get< 0 >( uOut ) - std::get< 0 >( uIn );
886 for(
int i = 0; i < RangeType::dimension; ++i )
887 std::get< 1 >( uJump )[ i ].axpy( std::get< 0 >( uJump )[ i ], normal );
889 auto intIn = [
this, xIn, uIn, xOut, uOut, normal, uJump ] (
const DomainValueType &phiIn ) {
890 const EntityType inside = intersection().inside();
891 const EntityType outside = intersection().outside();
893 const RangeFieldType half = RangeFieldType( 1 ) / RangeFieldType( 2 );
895 DomainValueType phiJump( 0, 0 );
896 std::get< 0 >( phiJump ) -= std::get< 0 >( phiIn );
897 for(
int i = 0; i < RangeType::dimension; ++i )
898 std::get< 1 >( phiJump )[ i ].axpy( std::get< 0 >( phiJump )[ i ], normal );
900 model().init( outside );
901 JacobianRangeType dFluxPrimeOut( 0 );
902 model().linFlux( std::get< 0 >( uOut ), std::get< 1 >( uJump ), xOut, 0, std::get< 1 >( phiJump ), dFluxPrimeOut );
904 model().init( inside );
905 JacobianRangeType dFluxIn( 0 ), dFluxPrimeIn( 0 );
906 model().linFlux( std::get< 0 >( uIn ), std::get< 1 >( uJump ), xIn, std::get< 0 >( phiIn ), std::get< 1 >( phiJump ), dFluxPrimeIn );
907 model().linFlux( std::get< 0 >( uIn ), 0, xIn, std::get< 0 >( phiIn ), 0, dFluxIn );
908 dFluxPrimeIn -= dFluxIn;
910 model().linFlux( std::get< 0 >( uIn ), std::get< 1 >( uIn ), xIn, std::get< 0 >( phiIn ), std::get< 1 >( phiIn ), dFluxIn );
912 RangeType int0 = std::get< 0 >( phiJump );
914 dFluxIn.usmv( -half, normal, int0 );
916 dFluxPrimeIn *= -half;
917 dFluxPrimeOut *= -half;
919 return std::make_pair( RangeValueType( -int0, dFluxPrimeIn ), RangeValueType( int0, dFluxPrimeOut ) );
922 auto intOut = [
this, xIn, uIn, xOut, uOut, normal, uJump ] (
const DomainValueType &phiOut ) {
923 const EntityType inside = intersection().inside();
924 const EntityType outside = intersection().outside();
926 const RangeFieldType half = RangeFieldType( 1 ) / RangeFieldType( 2 );
928 DomainValueType phiJump( 0, 0 );
929 std::get< 0 >( phiJump ) = std::get< 0 >( phiOut );
930 for(
int i = 0; i < RangeType::dimension; ++i )
931 std::get< 1 >( phiJump )[ i ].axpy( std::get< 0 >( phiJump )[ i ], normal );
933 model().init( outside );
934 JacobianRangeType dFluxOut( 0 ), dFluxPrimeOut( 0 );
935 model().linFlux( std::get< 0 >( uOut ), std::get< 1 >( uJump ), xOut, std::get< 0 >( phiOut ), std::get< 1 >( phiJump ), dFluxPrimeOut );
936 model().linFlux( std::get< 0 >( uOut ), 0, xOut, std::get< 0 >( phiOut ), 0, dFluxOut );
937 dFluxPrimeOut -= dFluxOut;
939 model().linFlux( std::get< 0 >( uOut ), std::get< 1 >( uOut ), xOut, std::get< 0 >( phiOut ), std::get< 1 >( phiOut ), dFluxOut );
941 model().init( inside );
942 JacobianRangeType dFluxPrimeIn( 0 );
943 model().linFlux( std::get< 0 >( uIn ), std::get< 1 >( uJump ), xIn, 0, std::get< 1 >( phiJump ), dFluxPrimeIn );
945 RangeType int0 = std::get< 0 >( phiJump );
947 dFluxOut.usmv( -half, normal, int0 );
949 dFluxPrimeIn *= -half;
950 dFluxPrimeOut *= -half;
952 return std::make_pair( RangeValueType( -int0, dFluxPrimeIn ), RangeValueType( int0, dFluxPrimeOut ) );
955 return std::make_pair( intIn, intOut );
958 const Model &model ()
const {
return *model_; }
961 const IntersectionType &intersection ()
const { assert( intersection_ );
return *intersection_; }
963 const Model *model_ =
nullptr;
964 RangeFieldType penalty_;
965 const IntersectionType *intersection_ =
nullptr;
966 RangeFieldType beta_;
IntersectionIteratorType::Intersection IntersectionType
type of intersection
Definition: adaptiveleafgridpart.hh:95
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Type traits to determine the type of reals (when working with complex numbers)
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:271
Dune namespace.
Definition: alignedallocator.hh:13
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
#define DUNE_PRIVATE
Mark a symbol as being for internal use within the current DSO only.
Definition: visibility.hh:28