DUNE-FEM (unstable)

integrands.hh
1#ifndef DUNE_FEM_SCHEMES_INTEGRANDS_HH
2#define DUNE_FEM_SCHEMES_INTEGRANDS_HH
3
4#include <cassert>
5
6#include <algorithm>
7#include <functional>
8#include <tuple>
9#include <type_traits>
10#include <utility>
11
12#include <dune/common/ftraits.hh>
13
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>
18
19namespace Dune
20{
21
22 namespace Fem
23 {
24
25 // IntegrandsTraits
26 // ----------------
27
28 namespace Impl
29 {
30
31 namespace IntegrandsTraits
32 {
33
34 template< class Integrands >
35 using DomainValueType = typename std::decay_t< decltype( std::ref( std::declval< const Integrands & >() ).get() ) >::DomainValueType;
36
37 template< class Integrands >
38 using RangeValueType = typename std::decay_t< decltype( std::ref( std::declval< const Integrands & >() ).get() ) >::RangeValueType;
39
40 template< class Integrands >
41 using GridPartType = typename std::decay_t< decltype( std::ref( std::declval< const Integrands & >() ).get() ) >::GridPartType;
42
43
44 template< class Integrands >
45 using EntityType = typename GridPartType< Integrands >::template Codim< 0 >::EntityType;
46
47 template< class Integrands >
48 using IntersectionType = typename GridPartType< Integrands >::IntersectionType;
49
50
51 template< class Integrands >
52 using InteriorQuadratureType = CachingQuadrature< GridPartType< Integrands >, 0 >;
53
54 template< class Integrands >
55 using SurfaceQuadratureType = CachingQuadrature< GridPartType< Integrands >, 1 >;
56
57
58 template< class Integrands >
59 using InteriorQuadraturePointType = QuadraturePointWrapper< InteriorQuadratureType< Integrands > >;
60
61 template< class Integrands >
62 using SurfaceQuadraturePointType = QuadraturePointWrapper< SurfaceQuadratureType< Integrands > >;
63
64
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 );
67
68 std::false_type interior ( ... );
69
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 & );
72
73 std::false_type hasInterior ( ... );
74
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 & );
77
78 std::false_type hasNonLinear ( ... );
79
80
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 );
83
84 std::false_type boundary ( ... );
85
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 & );
88
89 std::false_type hasBoundary ( ... );
90
91
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 );
94
95 std::false_type skeleton ( ... );
96
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 & );
99
100 std::false_type hasSkeleton ( ... );
101
102 template< class Integrands >
103 using Get = std::decay_t< decltype( std::ref( std::declval< const Integrands & >() ).get() ) >;
104
105 template< class Integrands >
106 using RRangeType = typename Get<Integrands>::RRangeType;
107 template< class Integrands >
108 using DirichletComponentType = typename Get<Integrands>::DirichletComponentType;
109 } // namespace IntegrandsTraits
110
111 } // namespace Impl
112
113 template< class Integrands >
114 struct IntegrandsTraits
115 {
116 typedef Impl::IntegrandsTraits::DomainValueType< Integrands > DomainValueType;
117 typedef Impl::IntegrandsTraits::RangeValueType< Integrands > RangeValueType;
118 typedef Impl::IntegrandsTraits::GridPartType< Integrands > GridPartType;
119
120 typedef Impl::IntegrandsTraits::EntityType< Integrands > EntityType;
121 typedef Impl::IntegrandsTraits::IntersectionType< Integrands > IntersectionType;
122
123 static const bool hasNonLinear = decltype( Impl::IntegrandsTraits::hasNonLinear( std::declval< const Integrands & >() ) )::value;
124
125 static const bool interior = decltype( Impl::IntegrandsTraits::interior( std::declval< const Integrands & >() ) )::value;
126 static const bool hasInterior = decltype( Impl::IntegrandsTraits::hasInterior( std::declval< const Integrands & >() ) )::value;
127
128 static const bool boundary = decltype( Impl::IntegrandsTraits::boundary( std::declval< const Integrands & >() ) )::value;
129 static const bool hasBoundary = decltype( Impl::IntegrandsTraits::hasBoundary( std::declval< const Integrands & >() ) )::value;
130
131 static const bool skeleton = decltype( Impl::IntegrandsTraits::skeleton( std::declval< const Integrands & >() ) )::value;
132 static const bool hasSkeleton = decltype( Impl::IntegrandsTraits::hasSkeleton( std::declval< const Integrands & >() ) )::value;
133
134 static_assert( (!hasInterior || interior), "Existence of method 'hasInterior' implies existence of method interior." );
135 static_assert( (!hasBoundary || boundary), "Existence of method 'hasBoundary' implies existence of method boundary." );
136 static_assert( (!hasSkeleton || skeleton), "Existence of method 'hasSkeleton' implies existence of method skeleton." );
137
138 static const bool isFull = hasInterior && hasBoundary && hasSkeleton;
139
140 typedef Impl::IntegrandsTraits::RRangeType< Integrands > RRangeType;
141 typedef Impl::IntegrandsTraits::DirichletComponentType< Integrands > DirichletComponentType;
142 };
143
144
145
146 // FullIntegrands
147 // --------------
148
149 template< class Integrands >
150 struct FullIntegrands
151 {
152 typedef typename IntegrandsTraits< Integrands >::DomainValueType DomainValueType;
153 typedef typename IntegrandsTraits< Integrands >::RangeValueType RangeValueType;
154 typedef typename IntegrandsTraits< Integrands >::GridPartType GridPartType;
155
156 typedef typename IntegrandsTraits< Integrands >::EntityType EntityType;
157 typedef typename IntegrandsTraits< Integrands >::IntersectionType IntersectionType;
158
159 private:
160 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasNonLinear, int > = 0 >
161 static bool nonlinear ( const T &integrands )
162 {
163 return integrands.nonlinear();
164 }
165
166 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasNonLinear, int > = 0 >
167 static bool nonlinear ( const T &integrands )
168 {
169 // default for nonlinear is true assuming that the model needs nonlinear solver
170 return true;
171 }
172
173 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasInterior, int > = 0 >
174 static bool hasInterior ( const T &integrands )
175 {
176 return integrands.hasInterior();
177 }
178
179 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasInterior, int > = 0 >
180 static bool hasInterior ( const T &integrands )
181 {
182 return IntegrandsTraits< T >::interior;
183 }
184
185 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::interior, int > = 0 >
186 static RangeValueType interior ( const T &integrands, const Point &x, const DomainValueType &u )
187 {
188 return integrands.interior( x, u );
189 }
190
191 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::interior, int > = 0 >
192 static RangeValueType interior ( const T &integrands, const Point &x, const DomainValueType &u )
193 {
194 return RangeValueType();
195 }
196
197 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::interior, int > = 0 >
198 static auto linearizedInterior ( const T &integrands, const Point &x, const DomainValueType &u )
199 {
200 return integrands.linearizedInterior( x, u );
201 }
202
203 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::interior, int > = 0 >
204 static auto linearizedInterior ( const T &integrands, const Point &x, const DomainValueType &u )
205 {
206 return [] ( const DomainValueType & ) { return RangeValueType(); };
207 }
208
209 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasBoundary, int > = 0 >
210 static bool hasBoundary ( const T &integrands )
211 {
212 return integrands.hasBoundary();
213 }
214
215 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasBoundary, int > = 0 >
216 static bool hasBoundary ( const T &integrands )
217 {
218 return IntegrandsTraits< T >::boundary;
219 }
220
221 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::boundary, int > = 0 >
222 static RangeValueType boundary ( const T &integrands, const Point &x, const DomainValueType &u )
223 {
224 return integrands.boundary( x, u );
225 }
226
227 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::boundary, int > = 0 >
228 static RangeValueType boundary ( const T &integrands, const Point &x, const DomainValueType &u )
229 {
230 return RangeValueType();
231 }
232
233 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::boundary, int > = 0 >
234 static auto linearizedBoundary ( const T &integrands, const Point &x, const DomainValueType &u )
235 {
236 return integrands.linearizedBoundary( x, u );
237 }
238
239 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::boundary, int > = 0 >
240 static auto linearizedBoundary ( const T &integrands, const Point &x, const DomainValueType &u )
241 {
242 return [] ( const DomainValueType & ) { return RangeValueType(); };
243 }
244
245 template< class T, std::enable_if_t< IntegrandsTraits< T >::hasSkeleton, int > = 0 >
246 static bool hasSkeleton ( const T &integrands )
247 {
248 return integrands.hasSkeleton();
249 }
250
251 template< class T, std::enable_if_t< !IntegrandsTraits< T >::hasSkeleton, int > = 0 >
252 static bool hasSkeleton ( const T &integrands )
253 {
254 return IntegrandsTraits< T >::skeleton;
255 }
256
257 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::skeleton, int > = 0 >
258 static std::pair< RangeValueType, RangeValueType > skeleton ( const T &integrands, const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut )
259 {
260 return integrands.skeleton( xIn, uIn, xOut, uOut );
261 }
262
263 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::skeleton, int > = 0 >
264 static std::pair< RangeValueType, RangeValueType > skeleton ( const T &integrands, const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut )
265 {
266 return std::make_pair( RangeValueType(), RangeValueType() );
267 }
268
269 template< class T, class Point, std::enable_if_t< IntegrandsTraits< T >::skeleton, int > = 0 >
270 static auto linearizedSkeleton ( const T &integrands, const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut )
271 {
272 return integrands.linearizedSkeleton( xIn, uIn, xOut, uOut );
273 }
274
275 template< class T, class Point, std::enable_if_t< !IntegrandsTraits< T >::skeleton, int > = 0 >
276 static auto linearizedSkeleton ( const T &integrands, const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut )
277 {
278 auto zero = [] ( const DomainValueType & ) { return std::make_pair( RangeValueType(), RangeValueType() ); };
279 return std::make_pair( zero, zero );
280 }
281
282 public:
283 template< class... Args >
284 explicit FullIntegrands ( Args &&... args )
285 : integrands_( std::forward< Args >( args )... ),
286 rInt_( std::ref( integrands_ ).get() )
287 {
288 }
289
290 bool init ( const EntityType &entity ) { return integrands().init( entity ); }
291 bool init ( const IntersectionType &intersection ) { return integrands().init( intersection ); }
292 void unbind ( ) { integrands().unbind( ); }
293
294 bool nonlinear() const { return nonlinear( integrands() ); }
295
296 bool hasInterior () const { return hasInterior( integrands() ); }
297
298 template< class Point >
299 RangeValueType interior ( const Point &x, const DomainValueType &u ) const
300 {
301 return interior( integrands(), x, u );
302 }
303
304 template< class Point >
305 auto linearizedInterior ( const Point &x, const DomainValueType &u ) const
306 {
307 return linearizedInterior( integrands(), x, u );
308 }
309
310 bool hasBoundary () const { return hasBoundary( integrands() ); }
311
312 template< class Point >
313 RangeValueType boundary ( const Point &x, const DomainValueType &u ) const
314 {
315 return boundary( integrands(), x, u );
316 }
317
318 template< class Point >
319 auto linearizedBoundary ( const Point &x, const DomainValueType &u ) const
320 {
321 return linearizedBoundary( integrands(), x, u );
322 }
323
324 bool hasSkeleton () const { return hasSkeleton( integrands() ); }
325
326 template< class Point >
327 std::pair< RangeValueType, RangeValueType > skeleton ( const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut ) const
328 {
329 return skeleton( integrands(), xIn, uIn, xOut, uOut );
330 }
331
332 template< class Point >
333 auto linearizedSkeleton ( const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut ) const
334 {
335 return linearizedSkeleton( integrands(), xIn, uIn, xOut, uOut );
336 }
337
338 protected:
339 typedef typename Integrands::type RealIntegrands;
340 //decltype( auto ) integrands () { return std::ref( integrands_ ).get(); }
341 //decltype( auto ) integrands () const { return std::ref( integrands_ ).get(); }
342 RealIntegrands& integrands () { return rInt_; }
343 const RealIntegrands& integrands () const { return rInt_; }
344
345
346
347 Integrands integrands_;
348 RealIntegrands rInt_;
349
350 public:
351 typedef typename IntegrandsTraits< Integrands >::RRangeType RRangeType;
352 typedef typename IntegrandsTraits< Integrands >::DirichletComponentType DirichletComponentType;
353 bool hasDirichletBoundary () const
354 {
355 return integrands().hasDirichletBoundary();
356 }
357 bool isDirichletIntersection( const IntersectionType& inter, DirichletComponentType &dirichletComponent ) const
358 {
359 return integrands().isDirichletIntersection(inter,dirichletComponent);
360 }
361 template <class Point>
362 void dirichlet( int bndId, const Point &x, RRangeType &value) const
363 {
364 return integrands().dirichlet(bndId,x,value);
365 }
366 };
367
368
369
370 // VirtualizedIntegrands
371 // ---------------------
372
373 namespace detail
374 {
375 template <class T>
376 struct GetDimRange;
377 template <class FT,int r>
378 struct GetDimRange<Dune::FieldVector<FT,r>>
379 {
380 typedef Dune::FieldVector<FT,r> type;
381 static const int value = r;
382 };
383 template <class FT,int r,int c>
384 struct GetDimRange<Dune::FieldMatrix<FT,r,c>>
385 {
386 typedef Dune::FieldVector<FT,r> type;
387 static const int value = r;
388 };
389 template <class FT,int r,int c>
390 struct GetDimRange<Dune::Fem::ExplicitFieldVector<Dune::FieldMatrix<FT,c,c>,r>>
391 {
392 typedef Dune::FieldVector<FT,r> type;
393 static const int value = r;
394 };
395 }
396
397 template< class GridPart, class DomainValue, class RangeValue = DomainValue >
398 class VirtualizedIntegrands
399 {
400 typedef VirtualizedIntegrands< GridPart, DomainValue, RangeValue > This;
401
402 public:
403 typedef GridPart GridPartType;
404 typedef DomainValue DomainValueType;
405 typedef RangeValue RangeValueType;
406
407 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
408 typedef typename GridPartType::IntersectionType IntersectionType;
409
410 using RRangeType = typename detail::GetDimRange<std::tuple_element_t<0,RangeValueType>>::type;
411 typedef std::array<int,RRangeType::dimension> DirichletComponentType;
412 typedef typename EntityType::Geometry::LocalCoordinate DomainType;
413
414 private:
415 typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
416
417 typedef FemPy::CachingPoint< GridPart, LocalCoordinateType, 0 > InteriorCachingPointType;
418 typedef FemPy::ElementPoint< GridPart, LocalCoordinateType, 0 > InteriorElementPointType;
419 typedef FemPy::CachingPoint< GridPart, LocalCoordinateType, 1 > SurfaceCachingPointType;
420 typedef FemPy::ElementPoint< GridPart, LocalCoordinateType, 1 > SurfaceElementPointType;
421
422 template< class QP >
423 static Fem::QuadraturePointWrapper< QP > asQP ( const QP &qp )
424 {
425 return static_cast< Fem::QuadraturePointWrapper< QP > >( qp );
426 }
427
428 template< class R >
429 using Linearization = std::function< R( const DomainValueType & ) >;
430
431 template< class R >
432 using LinearizationPair = std::pair< Linearization< std::pair< R, R > >, Linearization< std::pair< R, R > > >;
433
434 struct Interface
435 {
436 virtual ~Interface () {}
437 virtual Interface *clone () const = 0;
438
439 virtual bool init ( const EntityType &entity ) = 0;
440 virtual bool init ( const IntersectionType &intersection ) = 0;
441 virtual void unbind ( ) = 0;
442
443 virtual bool nonlinear() const = 0;
444
445 virtual bool hasInterior () const = 0;
446 virtual RangeValueType interior ( const InteriorCachingPointType &x, const DomainValueType &u ) const = 0;
447 virtual RangeValueType interior ( const InteriorElementPointType &x, const DomainValueType &u ) const = 0;
448 virtual Linearization< RangeValueType > linearizedInterior ( const InteriorCachingPointType &x, const DomainValueType &u ) const = 0;
449 virtual Linearization< RangeValueType > linearizedInterior ( const InteriorElementPointType &x, const DomainValueType &u ) const = 0;
450
451 virtual bool hasBoundary () const = 0;
452 virtual RangeValueType boundary ( const SurfaceCachingPointType &x, const DomainValueType &u ) const = 0;
453 virtual RangeValueType boundary ( const SurfaceElementPointType &x, const DomainValueType &u ) const = 0;
454 virtual Linearization< RangeValueType > linearizedBoundary ( const SurfaceCachingPointType &x, const DomainValueType &u ) const = 0;
455 virtual Linearization< RangeValueType > linearizedBoundary ( const SurfaceElementPointType &x, const DomainValueType &u ) const = 0;
456
457 virtual bool hasSkeleton () const = 0;
458 virtual std::pair< RangeValueType, RangeValueType > skeleton ( const SurfaceCachingPointType &xIn, const DomainValueType &uIn, const SurfaceCachingPointType &xOut, const DomainValueType &uOut ) const = 0;
459 virtual std::pair< RangeValueType, RangeValueType > skeleton ( const SurfaceElementPointType &xIn, const DomainValueType &uIn, const SurfaceElementPointType &xOut, const DomainValueType &uOut ) const = 0;
460 virtual LinearizationPair< RangeValueType > linearizedSkeleton ( const SurfaceCachingPointType &xIn, const DomainValueType &uIn, const SurfaceCachingPointType &xOut, const DomainValueType &uOut ) const = 0;
461 virtual LinearizationPair< RangeValueType > linearizedSkeleton ( const SurfaceElementPointType &xIn, const DomainValueType &uIn, const SurfaceElementPointType &xOut, const DomainValueType &uOut ) const = 0;
462
463 virtual bool hasDirichletBoundary () const = 0;
464 virtual bool isDirichletIntersection( const IntersectionType& inter, DirichletComponentType &dirichletComponent ) const = 0;
465 virtual void dirichlet( int bndId, const DomainType &x,RRangeType &value) const = 0;
466 };
467
468 template< class Impl >
469 struct DUNE_PRIVATE Implementation final
470 : public Interface
471 {
472 Implementation ( Impl impl ) : impl_( std::move( impl ) )
473 {
474 }
475 virtual Interface *clone () const override { return new Implementation( *this ); }
476
477 virtual bool init ( const EntityType &entity ) override { return impl().init( entity ); }
478 virtual bool init ( const IntersectionType &intersection ) override { return impl().init( intersection ); }
479 virtual void unbind ( ) override { impl().unbind( ); }
480
481 virtual bool nonlinear () const override { return impl().nonlinear(); }
482
483 virtual bool hasInterior () const override { return impl().hasInterior(); }
484 virtual RangeValueType interior ( const InteriorCachingPointType &x, const DomainValueType &u ) const override { return impl().interior( asQP( x ), u ); }
485 virtual RangeValueType interior ( const InteriorElementPointType &x, const DomainValueType &u ) const override { return impl().interior( asQP( x ), u ); }
486 virtual Linearization< RangeValueType > linearizedInterior ( const InteriorCachingPointType &x, const DomainValueType &u ) const override { return impl().linearizedInterior( asQP( x ), u ); }
487 virtual Linearization< RangeValueType > linearizedInterior ( const InteriorElementPointType &x, const DomainValueType &u ) const override { return impl().linearizedInterior( asQP( x ), u ); }
488
489 virtual bool hasBoundary () const override { return impl().hasBoundary(); }
490 virtual RangeValueType boundary ( const SurfaceCachingPointType &x, const DomainValueType &u ) const override { return impl().boundary( asQP( x ), u ); }
491 virtual RangeValueType boundary ( const SurfaceElementPointType &x, const DomainValueType &u ) const override { return impl().boundary( asQP( x ), u ); }
492 virtual Linearization< RangeValueType > linearizedBoundary ( const SurfaceCachingPointType &x, const DomainValueType &u ) const override { return impl().linearizedBoundary( asQP( x ), u ); }
493 virtual Linearization< RangeValueType > linearizedBoundary ( const SurfaceElementPointType &x, const DomainValueType &u ) const override { return impl().linearizedBoundary( asQP( x ), u ); }
494
495 virtual bool hasSkeleton () const override { return impl().hasSkeleton(); }
496 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 ); }
497 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 ); }
498 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 ); }
499 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 ); }
500
501 virtual bool hasDirichletBoundary () const override { return impl().hasDirichletBoundary(); }
502 virtual bool isDirichletIntersection( const IntersectionType& inter, DirichletComponentType &dirichletComponent ) const override { return impl().isDirichletIntersection(inter,dirichletComponent); }
503 virtual void dirichlet( int bndId, const DomainType &x,RRangeType &value) const override { impl().dirichlet(bndId,x,value); }
504
505 private:
506 const auto &impl () const { return std::cref( impl_ ).get(); }
507 auto &impl () { return std::ref( impl_ ).get(); }
508
509 Impl impl_;
510 };
511
512 template< class Integrands >
513 using isVirtualized = std::is_same< std::decay_t< decltype( std::ref( std::declval< Integrands & >() ).get() ) >, This >;
514
515 public:
516 template< class Integrands, std::enable_if_t< IntegrandsTraits< std::decay_t< Integrands > >::isFull && !isVirtualized< Integrands >::value, int > = 0 >
517 explicit VirtualizedIntegrands ( Integrands integrands )
518 : impl_( new Implementation< Integrands >( std::move( integrands ) ) )
519 {}
520
521 template< class Integrands, std::enable_if_t< !IntegrandsTraits< Integrands >::isFull && !isVirtualized< Integrands >::value, int > = 0 >
522 explicit VirtualizedIntegrands ( Integrands integrands )
523 : VirtualizedIntegrands( FullIntegrands< std::decay_t< Integrands > >( std::move( integrands ) ) )
524 {}
525
526 VirtualizedIntegrands ( const This &other ) : impl_( other ? other.impl().clone() : nullptr )
527 {}
528
529 VirtualizedIntegrands ( This && ) = default;
530
531 VirtualizedIntegrands &operator= ( const This &other )
532 {
533 impl_.reset( other ? other.impl().clone() : nullptr );
534 }
535
536 VirtualizedIntegrands &operator= ( This && ) = default;
537
538 explicit operator bool () const { return static_cast< bool >( impl_ ); }
539
540 bool init ( const EntityType &entity ) { return impl().init( entity ); }
541 bool init ( const IntersectionType &intersection ) { return impl().init( intersection ); }
542 void unbind ( ) { impl().unbind( ); }
543
544 bool nonlinear() const { return impl().nonlinear(); }
545
546 bool hasInterior () const { return impl().hasInterior(); }
547
548 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
549 RangeValueType interior ( const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u ) const
550 {
551 return impl().interior( InteriorCachingPointType( x ), u );
552 }
553
554 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
555 RangeValueType interior ( const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u ) const
556 {
557 return impl().interior( InteriorElementPointType( x ), u );
558 }
559
560 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
561 auto linearizedInterior ( const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u ) const
562 {
563 return impl().linearizedInterior( InteriorCachingPointType( x ), u );
564 }
565
566 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
567 auto linearizedInterior ( const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u ) const
568 {
569 return impl().linearizedInterior( InteriorElementPointType( x ), u );
570 }
571
572 bool hasBoundary () const { return impl().hasBoundary(); }
573
574 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
575 RangeValueType boundary ( const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u ) const
576 {
577 return impl().boundary( SurfaceCachingPointType( x ), u );
578 }
579
580 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
581 RangeValueType boundary ( const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u ) const
582 {
583 return impl().boundary( SurfaceElementPointType( x ), u );
584 }
585
586 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
587 auto linearizedBoundary ( const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u ) const
588 {
589 return impl().linearizedBoundary( SurfaceCachingPointType( x ), u );
590 }
591
592 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
593 auto linearizedBoundary ( const Fem::QuadraturePointWrapper< Quadrature > &x, const DomainValueType &u ) const
594 {
595 return impl().linearizedBoundary( SurfaceElementPointType( x ), u );
596 }
597
598 bool hasSkeleton () const { return impl().hasSkeleton(); }
599
600 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
601 std::pair< RangeValueType, RangeValueType > skeleton ( const Fem::QuadraturePointWrapper< Quadrature > &xIn, const DomainValueType &uIn, const Fem::QuadraturePointWrapper< Quadrature > &xOut, const DomainValueType &uOut ) const
602 {
603 return impl().skeleton( SurfaceCachingPointType( xIn ), uIn, SurfaceCachingPointType( xOut ), uOut );
604 }
605
606 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
607 std::pair< RangeValueType, RangeValueType > skeleton ( const Fem::QuadraturePointWrapper< Quadrature > &xIn, const DomainValueType &uIn, const Fem::QuadraturePointWrapper< Quadrature > &xOut, const DomainValueType &uOut ) const
608 {
609 return impl().skeleton( SurfaceElementPointType( xIn ), uIn, SurfaceElementPointType( xOut ), uOut );
610 }
611
612 template< class Quadrature, std::enable_if_t< std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
613 auto linearizedSkeleton ( const Fem::QuadraturePointWrapper< Quadrature > &xIn, const DomainValueType &uIn, const Fem::QuadraturePointWrapper< Quadrature > &xOut, const DomainValueType &uOut ) const
614 {
615 return impl().linearizedSkeleton( SurfaceCachingPointType( xIn ), uIn, SurfaceCachingPointType( xOut ), uOut );
616 }
617
618 template< class Quadrature, std::enable_if_t< !std::is_convertible< Quadrature, Fem::CachingInterface >::value, int > = 0 >
619 auto linearizedSkeleton ( const Fem::QuadraturePointWrapper< Quadrature > &xIn, const DomainValueType &uIn, const Fem::QuadraturePointWrapper< Quadrature > &xOut, const DomainValueType &uOut ) const
620 {
621 return impl().linearizedSkeleton( SurfaceElementPointType( xIn ), uIn, SurfaceElementPointType( xOut ), uOut );
622 }
623
624 bool hasDirichletBoundary () const
625 {
626 return impl().hasDirichletBoundary();
627 }
628 bool isDirichletIntersection( const IntersectionType& inter, DirichletComponentType &dirichletComponent ) const
629 {
630 return impl().isDirichletIntersection(inter,dirichletComponent);
631 }
632 void dirichlet( int bndId, const DomainType &x,RRangeType &value) const
633 {
634 return impl().dirichlet(bndId,x,value);
635 }
636
637 private:
638 const Interface &impl () const { assert( impl_ ); return *impl_; }
639 Interface &impl () { assert( impl_ ); return *impl_; }
640
641 std::unique_ptr< Interface > impl_;
642 };
643
644
645
646 // ConservationLawModelIntegrands
647 // ------------------------------
648
649 template< class Model >
650 struct ConservationLawModelIntegrands
651 {
652 typedef Model ModelType;
653 typedef typename Model::GridPartType GridPartType;
654
655 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
656 typedef typename GridPartType::IntersectionType IntersectionType;
657
658 typedef typename Model::RangeType RangeType;
659 typedef typename Model::JacobianRangeType JacobianRangeType;
660
661 typedef std::tuple< RangeType, JacobianRangeType > DomainValueType;
662 typedef std::tuple< RangeType, JacobianRangeType > RangeValueType;
663
664 explicit ConservationLawModelIntegrands ( const Model &model ) : model_( &model ) {}
665
666 bool init ( const EntityType &entity ) { return model().init( entity ); }
667
668 bool init ( const IntersectionType &intersection )
669 {
670 return (intersection.boundary() && model().hasNeumanBoundary() && model().init( intersection.inside() ));
671 }
672
673 void unbind ( ) { model().unbind( ); }
674
675 template< class Point >
676 RangeValueType interior ( const Point &x, const DomainValueType &u ) const
677 {
678 RangeType source( 0 );
679 model().source( x, std::get< 0 >( u ), std::get< 1 >( u ), source );
680 JacobianRangeType dFlux( 0 );
681 model().flux( x, std::get< 0 >( u ), std::get< 1 >( u ), dFlux );
682 return std::make_tuple( source, dFlux );
683 }
684
685 template< class Point >
686 auto linearizedInterior ( const Point &x, const DomainValueType &u ) const
687 {
688 return [ this, x, u ] ( const DomainValueType &phi ) {
689 RangeType source( 0 );
690 model().linSource( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), source );
691 JacobianRangeType dFlux( 0 );
692 model().linFlux( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), dFlux );
693 return std::make_tuple( source, dFlux );
694 };
695 }
696
697 template< class Point >
698 RangeValueType boundary ( const Point &x, const DomainValueType &u ) const
699 {
700 RangeType alpha( 0 );
701 model().alpha( x, std::get< 0 >( u ), alpha );
702 return std::make_tuple( alpha, 0 );
703 }
704
705 template< class Point >
706 auto linearizedBoundary ( const Point &x, const DomainValueType &u ) const
707 {
708 return [ this, x, u ] ( const DomainValueType &phi ) {
709 RangeType alpha( 0 );
710 model().linAlpha( std::get< 0 >( u ), x, std::get< 0 >( phi ), alpha );
711 return std::make_tuple( alpha, 0 );
712 };
713 }
714
715 const Model &model () const { return *model_; }
716
717 private:
718 const Model *model_ = nullptr;
719 };
720
721
722
723 // DGConservationLawModelIntegrands
724 // --------------------------------
725
726 template< class Model >
727 struct DGConservationLawModelIntegrands
728 {
729 typedef Model ModelType;
730 typedef typename Model::GridPartType GridPartType;
731
732 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
733 typedef typename GridPartType::IntersectionType IntersectionType;
734
735 typedef typename Model::RangeType RangeType;
736 typedef typename Model::JacobianRangeType JacobianRangeType;
737
738 typedef typename FieldTraits< RangeType >::field_type RangeFieldType;
739
740 typedef std::tuple< RangeType, JacobianRangeType > DomainValueType;
741 typedef std::tuple< RangeType, JacobianRangeType > RangeValueType;
742
743 DGConservationLawModelIntegrands ( const Model &model, RangeFieldType penalty )
744 : model_( &model ), penalty_( penalty )
745 {}
746
747 bool init ( const EntityType &entity )
748 {
749 intersection_ = nullptr;
750 return model().init( entity );
751 }
752
753 bool init ( const IntersectionType &intersection )
754 {
755 intersection_ = &intersection;
756 if( intersection.boundary() )
757 {
758 const EntityType inside = intersection.inside();
759 beta_ = penalty_ * intersection.geometry().volume() / inside.geometry().volume();
760 return (model().hasNeumanBoundary() && model().init( inside ));
761 }
762 else if( intersection.neighbor() )
763 {
764 const auto volIn = intersection.inside().geometry().volume();
765 const auto volOut = intersection.outside().geometry().volume();
766 beta_ = penalty_ * intersection.geometry().volume() / std::min( volIn, volOut );
767 return true;
768 }
769 }
770
771 void unbind ( )
772 {
773 model().unbind( );
774 }
775
776 template< class Point >
777 RangeValueType interior ( const Point &x, const DomainValueType &u ) const
778 {
779 RangeType source( 0 );
780 model().source( x, std::get< 0 >( u ), std::get< 1 >( u ), source );
781 JacobianRangeType dFlux( 0 );
782 model().flux( x, std::get< 0 >( u ), std::get< 1 >( u ), dFlux );
783 return RangeValueType( source, dFlux );
784 }
785
786 template< class Point >
787 auto linearizedInterior ( const Point &x, const DomainValueType &u ) const
788 {
789 return [ this, x, u ] ( const DomainValueType &phi ) {
790 RangeType source( 0 );
791 model().linSource( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), source );
792 JacobianRangeType dFlux( 0 );
793 model().linFlux( std::get< 0 >( u ), std::get< 1 >( u ), x, std::get< 0 >( phi ), std::get< 1 >( phi ), dFlux );
794 return RangeValueType( source, dFlux );
795 };
796 }
797
798 template< class Point >
799 RangeValueType boundary ( const Point &x, const DomainValueType &u ) const
800 {
801 RangeType alpha( 0 );
802 model().alpha( x, std::get< 0 >( u ), alpha );
803 return RangeValueType( alpha, 0 );
804 }
805
806 template< class Point >
807 auto linearizedBoundary ( const Point &x, const DomainValueType &u ) const
808 {
809 return [ this, x, u ] ( const DomainValueType &phi ) {
810 RangeType alpha( 0 );
811 model().linAlpha( std::get< 0 >( u ), x, std::get< 0 >( phi ), alpha );
812 return RangeValueType( alpha, 0 );
813 };
814 }
815
816 template< class Point >
817 std::pair< RangeValueType, RangeValueType > skeleton ( const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut ) const
818 {
819 const EntityType inside = intersection().inside();
820 const EntityType outside = intersection().outside();
821
822 const RangeFieldType half = RangeFieldType( 1 ) / RangeFieldType( 2 );
823 const auto normal = intersection().unitOuterNormal( xIn.localPosition() );
824
825 DomainValueType uJump( 0, 0 );
826 std::get< 0 >( uJump ) = std::get< 0 >( uOut ) - std::get< 0 >( uIn );
827 for( int i = 0; i < RangeType::dimension; ++i )
828 std::get< 1 >( uJump )[ i ].axpy( std::get< 0 >( uJump )[ i ], normal );
829
830 model().init( outside );
831 JacobianRangeType dFluxOut( 0 ), dFluxPrimeOut( 0 );
832 model().flux( xOut, std::get< 0 >( uOut ), std::get< 1 >( uJump ), dFluxPrimeOut );
833 model().flux( xOut, std::get< 0 >( uOut ), 0, dFluxOut );
834 dFluxPrimeOut -= dFluxOut;
835 dFluxOut = 0;
836 model().flux( xOut, std::get< 0 >( uOut ), std::get< 1 >( uOut ), dFluxOut );
837
838 model().init( inside );
839 JacobianRangeType dFluxIn( 0 ), dFluxPrimeIn( 0 );
840 model().flux( xIn, std::get< 0 >( uIn ), std::get< 1 >( uJump ), dFluxPrimeIn );
841 model().flux( xIn, std::get< 0 >( uIn ), 0, dFluxIn );
842 dFluxPrimeIn -= dFluxIn;
843 dFluxIn = 0;
844 model().flux( xIn, std::get< 0 >( uIn ), std::get< 1 >( uIn ), dFluxIn );
845
846 RangeType int0 = std::get< 0 >( uJump );
847 int0 *= beta_;
848 dFluxIn += dFluxOut;
849 dFluxIn.usmv( -half, normal, int0 );
850
851 dFluxPrimeIn *= -half;
852 dFluxPrimeOut *= -half;
853
854 return std::make_pair( RangeValueType( -int0, dFluxPrimeIn ), RangeValueType( int0, dFluxPrimeOut ) );
855 }
856
857 template< class Point >
858 auto linearizedSkeleton ( const Point &xIn, const DomainValueType &uIn, const Point &xOut, const DomainValueType &uOut ) const
859 {
860 const auto normal = intersection().unitOuterNormal( xIn.localPosition() );
861
862 DomainValueType uJump( 0, 0 );
863 std::get< 0 >( uJump ) = std::get< 0 >( uOut ) - std::get< 0 >( uIn );
864 for( int i = 0; i < RangeType::dimension; ++i )
865 std::get< 1 >( uJump )[ i ].axpy( std::get< 0 >( uJump )[ i ], normal );
866
867 auto intIn = [ this, xIn, uIn, xOut, uOut, normal, uJump ] ( const DomainValueType &phiIn ) {
868 const EntityType inside = intersection().inside();
869 const EntityType outside = intersection().outside();
870
871 const RangeFieldType half = RangeFieldType( 1 ) / RangeFieldType( 2 );
872
873 DomainValueType phiJump( 0, 0 );
874 std::get< 0 >( phiJump ) -= std::get< 0 >( phiIn );
875 for( int i = 0; i < RangeType::dimension; ++i )
876 std::get< 1 >( phiJump )[ i ].axpy( std::get< 0 >( phiJump )[ i ], normal );
877
878 model().init( outside );
879 JacobianRangeType dFluxPrimeOut( 0 );
880 model().linFlux( std::get< 0 >( uOut ), std::get< 1 >( uJump ), xOut, 0, std::get< 1 >( phiJump ), dFluxPrimeOut );
881
882 model().init( inside );
883 JacobianRangeType dFluxIn( 0 ), dFluxPrimeIn( 0 );
884 model().linFlux( std::get< 0 >( uIn ), std::get< 1 >( uJump ), xIn, std::get< 0 >( phiIn ), std::get< 1 >( phiJump ), dFluxPrimeIn );
885 model().linFlux( std::get< 0 >( uIn ), 0, xIn, std::get< 0 >( phiIn ), 0, dFluxIn );
886 dFluxPrimeIn -= dFluxIn;
887 dFluxIn = 0;
888 model().linFlux( std::get< 0 >( uIn ), std::get< 1 >( uIn ), xIn, std::get< 0 >( phiIn ), std::get< 1 >( phiIn ), dFluxIn );
889
890 RangeType int0 = std::get< 0 >( phiJump );
891 int0 *= beta_;
892 dFluxIn.usmv( -half, normal, int0 );
893
894 dFluxPrimeIn *= -half;
895 dFluxPrimeOut *= -half;
896
897 return std::make_pair( RangeValueType( -int0, dFluxPrimeIn ), RangeValueType( int0, dFluxPrimeOut ) );
898 };
899
900 auto intOut = [ this, xIn, uIn, xOut, uOut, normal, uJump ] ( const DomainValueType &phiOut ) {
901 const EntityType inside = intersection().inside();
902 const EntityType outside = intersection().outside();
903
904 const RangeFieldType half = RangeFieldType( 1 ) / RangeFieldType( 2 );
905
906 DomainValueType phiJump( 0, 0 );
907 std::get< 0 >( phiJump ) = std::get< 0 >( phiOut );
908 for( int i = 0; i < RangeType::dimension; ++i )
909 std::get< 1 >( phiJump )[ i ].axpy( std::get< 0 >( phiJump )[ i ], normal );
910
911 model().init( outside );
912 JacobianRangeType dFluxOut( 0 ), dFluxPrimeOut( 0 );
913 model().linFlux( std::get< 0 >( uOut ), std::get< 1 >( uJump ), xOut, std::get< 0 >( phiOut ), std::get< 1 >( phiJump ), dFluxPrimeOut );
914 model().linFlux( std::get< 0 >( uOut ), 0, xOut, std::get< 0 >( phiOut ), 0, dFluxOut );
915 dFluxPrimeOut -= dFluxOut;
916 dFluxOut = 0;
917 model().linFlux( std::get< 0 >( uOut ), std::get< 1 >( uOut ), xOut, std::get< 0 >( phiOut ), std::get< 1 >( phiOut ), dFluxOut );
918
919 model().init( inside );
920 JacobianRangeType dFluxPrimeIn( 0 );
921 model().linFlux( std::get< 0 >( uIn ), std::get< 1 >( uJump ), xIn, 0, std::get< 1 >( phiJump ), dFluxPrimeIn );
922
923 RangeType int0 = std::get< 0 >( phiJump );
924 int0 *= beta_;
925 dFluxOut.usmv( -half, normal, int0 );
926
927 dFluxPrimeIn *= -half;
928 dFluxPrimeOut *= -half;
929
930 return std::make_pair( RangeValueType( -int0, dFluxPrimeIn ), RangeValueType( int0, dFluxPrimeOut ) );
931 };
932
933 return std::make_pair( intIn, intOut );
934 }
935
936 const Model &model () const { return *model_; }
937
938 private:
939 const IntersectionType &intersection () const { assert( intersection_ ); return *intersection_; }
940
941 const Model *model_ = nullptr;
942 RangeFieldType penalty_;
943 const IntersectionType *intersection_ = nullptr;
944 RangeFieldType beta_;
945 };
946
947 } // namespace Fem
948
949} // namespace Dune
950
951#endif // #ifndef DUNE_FEM_SCHEMES_INTEGRANDS_HH
Dune::Fem::AdaptiveGridPartBase< AdaptiveLeafGridPartTraits< Grid, idxpitype, onlyCodimensionZero > >::IntersectionType
IntersectionIteratorType::Intersection IntersectionType
type of intersection
Definition: adaptiveleafgridpart.hh:95
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:95
ftraits.hh
Type traits to determine the type of reals (when working with complex numbers)
Dune::Partitions::interior
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:271
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::get
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
std
STL namespace.
DUNE_PRIVATE
#define DUNE_PRIVATE
Mark a symbol as being for internal use within the current DSO only.
Definition: visibility.hh:28
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 27, 22:29, 2024)