DUNE-FEM (unstable)

galerkin.hh
1#ifndef DUNE_FEM_SCHEMES_GALERKIN_HH
2#define DUNE_FEM_SCHEMES_GALERKIN_HH
3
4#include <cstddef>
5
6#include <tuple>
7#include <type_traits>
8#include <utility>
9#include <shared_mutex>
10#include <vector>
11#include <memory>
12
13#include <dune/common/hybridutilities.hh>
14#include <dune/common/timer.hh>
15
16#include <dune/grid/common/rangegenerators.hh>
17
18#include <dune/fem/function/localfunction/temporary.hh>
19#include <dune/fem/io/parameter/reader.hh>
20#include <dune/fem/operator/common/automaticdifferenceoperator.hh>
21#include <dune/fem/operator/common/differentiableoperator.hh>
22#include <dune/fem/operator/common/operator.hh>
23#include <dune/fem/operator/common/stencil.hh>
24#include <dune/fem/operator/common/temporarylocalmatrix.hh>
25#include <dune/fem/quadrature/cachingquadrature.hh>
26#include <dune/fem/quadrature/intersectionquadrature.hh>
27#include <dune/fem/common/bindguard.hh>
28
29#include <dune/fem/misc/threads/threaditerator.hh>
30#include <dune/fem/misc/threads/threadsafevalue.hh>
31#include <dune/fem/misc/hasboundaryintersection.hh>
32#include <dune/fem/misc/griddeclaration.hh>
33
34#include <dune/fem/operator/common/localmatrixcolumn.hh>
35#include <dune/fem/operator/common/localcontribution.hh>
36#include <dune/fem/operator/1order/localmassmatrix.hh>
37#include <dune/fem/schemes/integrands.hh>
38#include <dune/fem/schemes/dirichletwrapper.hh>
39#include <dune/fem/schemes/femscheme.hh>
40
41#include <dune/fem/space/common/capabilities.hh>
42
43// fempy includes
44#include <dune/fempy/quadrature/fempyquadratures.hh>
45
46namespace Dune
47{
48
49 namespace Fem
50 {
51
52 namespace Impl
53 {
54 template <class M>
55 class CallOrder
56 {
57
58 template <class F>
59 static int callOrder(const F& f, char)
60 {
61#ifndef NDEBUG
62 std::cerr << "WARNING: no order method available on " << typeid(F).name() << ", defaulting to 1!" << std::endl;
63#endif
64 return 1;
65 }
66
67 template <class F>
68 static auto callOrder(const F& f, int) -> decltype( f.order() )
69 {
70 return f.order();
71 }
72
73 public:
74 template <class F>
75 static int order (const F& f ) { return callOrder(f, 0); }
76 };
77
78 // GalerkinOperator
79 // ----------------
80
81 template <class Space>
82 struct DefaultGalerkinOperatorQuadratureSelector
83 {
84 typedef typename Space :: GridPartType GridPartType;
85 typedef typename GridPartType :: GridType GridType;
86
87 template <class Grid>
88 struct AnotherQuadSelector
89 {
90 typedef CachingQuadrature< GridPartType, 0, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > InteriorQuadratureType;
91 typedef CachingQuadrature< GridPartType, 1, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > SurfaceQuadratureType;
92 };
93
94
95#if HAVE_DUNE_POLYGONGRID
96 template < class ct >
97 struct AnotherQuadSelector< Dune::PolygonGrid< ct > >
98 {
99 typedef ElementQuadrature< GridPartType, 0, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > InteriorQuadratureType;
100 typedef ElementQuadrature< GridPartType, 1, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > SurfaceQuadratureType;
101 };
102#endif
103
104 typedef typename AnotherQuadSelector< GridType >::InteriorQuadratureType InteriorQuadratureType;
105 typedef typename AnotherQuadSelector< GridType >::SurfaceQuadratureType SurfaceQuadratureType;
106 };
107
108 // LocalGalerkinOperator
109 // ---------------------
110
111 template< class Integrands, template <class> class QuadSelector = DefaultGalerkinOperatorQuadratureSelector >
112 struct LocalGalerkinOperator
113 {
114 typedef LocalGalerkinOperator<Integrands> ThisType;
115 typedef std::conditional_t< Fem::IntegrandsTraits< Integrands >::isFull, Integrands, FullIntegrands< Integrands > > IntegrandsType;
116
117 typedef typename IntegrandsType::GridPartType GridPartType;
118
119 typedef typename GridPartType::ctype ctype;
120 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
121
122 // typedef QuadratureSelector
123 template <class Space>
124 using QuadratureSelector = QuadSelector< Space >;
125
126 // constructor
127 template< class... Args >
128 explicit LocalGalerkinOperator ( const GridPartType &gridPart, Args &&... args )
129 : gridPart_( gridPart ),
130 integrands_( std::forward< Args >( args )... ),
131 defaultInteriorOrder_( [] (const int order) { return 2 * order; } ),
132 defaultSurfaceOrder_ ( [] (const int order) { return 2 * order + 1; } ),
133 interiorQuadOrder_(0), surfaceQuadOrder_(0)
134 {
135 }
136
137 protected:
138 typedef typename IntegrandsType::DomainValueType DomainValueType;
139 typedef typename IntegrandsType::RangeValueType RangeValueType;
140 typedef std::make_index_sequence< std::tuple_size< DomainValueType >::value > DomainValueIndices;
141 typedef std::make_index_sequence< std::tuple_size< RangeValueType >::value > RangeValueIndices;
142
143
144 template< std::size_t... i >
145 static auto makeDomainValueVector ( std::size_t maxNumLocalDofs, std::index_sequence< i... > )
146 {
147 return std::make_tuple( std::vector< std::tuple_element_t< i, DomainValueType > >( maxNumLocalDofs )... );
148 }
149
150 static auto makeDomainValueVector ( std::size_t maxNumLocalDofs )
151 {
152 return makeDomainValueVector( maxNumLocalDofs, DomainValueIndices() );
153 }
154
155 template< std::size_t... i >
156 static auto makeRangeValueVector ( std::size_t maxNumLocalDofs, std::index_sequence< i... > )
157 {
158 return std::make_tuple( std::vector< std::tuple_element_t< i, RangeValueType > >( maxNumLocalDofs )... );
159 }
160
161 static auto makeRangeValueVector ( std::size_t maxNumLocalDofs )
162 {
163 return makeRangeValueVector( maxNumLocalDofs, RangeValueIndices() );
164 }
165
166 typedef decltype( makeDomainValueVector( 0u ) ) DomainValueVectorType;
167 typedef decltype( makeRangeValueVector( 0u ) ) RangeValueVectorType;
168
169 static void resizeDomainValueVector ( DomainValueVectorType& vec, const std::size_t size )
170 {
171 Hybrid::forEach( DomainValueIndices(), [ &vec, &size ] ( auto i ) {
172 std::get< i >( vec ).resize( size );
173 } );
174 }
175
176 static void resizeRangeValueVector ( RangeValueVectorType& vec, const std::size_t size )
177 {
178 Hybrid::forEach( RangeValueIndices(), [ &vec, &size ] ( auto i ) {
179 std::get< i >( vec ).resize( size );
180 } );
181 }
182
183 public:
184 void prepare( const std::size_t size ) const
185 {
186 resizeDomainValueVector( phiIn_, size );
187 resizeDomainValueVector( phiOut_, size );
188 resizeDomainValueVector( basisValues_, size );
189 resizeDomainValueVector( domainValues_, size );
190 resizeDomainValueVector( domainValuesOut_, size );
191 }
192
193 template< class LocalFunction, class Quadrature >
194 static void evaluateQuadrature ( const LocalFunction &u, const Quadrature &quad, std::vector< typename LocalFunction::RangeType > &phi )
195 {
196 u.evaluateQuadrature( quad, phi );
197 }
198
199 template< class LocalFunction, class Quadrature>
200 static void evaluateQuadrature ( const LocalFunction &u, const Quadrature &quad, std::vector< typename LocalFunction::JacobianRangeType > &phi )
201 {
202 u.jacobianQuadrature( quad, phi );
203 }
204
205 template< class LocalFunction, class Quadrature >
206 static void evaluateQuadrature ( const LocalFunction &u, const Quadrature &quad, std::vector< typename LocalFunction::HessianRangeType > &phi )
207 {
208 u.hessianQuadrature( quad, phi );
209 }
210
211 protected:
212 template< class LocalFunction, class Point >
213 static void value ( const LocalFunction &u, const Point &x, typename LocalFunction::RangeType &phi )
214 {
215 u.evaluate( x, phi );
216 }
217
218 template< class LocalFunction, class Point >
219 static void value ( const LocalFunction &u, const Point &x, typename LocalFunction::JacobianRangeType &phi )
220 {
221 u.jacobian( x, phi );
222 }
223
224 template< class LocalFunction, class Point >
225 static void value ( const LocalFunction &u, const Point &x, typename LocalFunction::HessianRangeType &phi )
226 {
227 u.hessian( x, phi );
228 }
229
230 template< class LocalFunction, class Point, class... T >
231 static void value ( const LocalFunction &u, const Point &x, std::tuple< T... > &phi )
232 {
233 Hybrid::forEach( std::index_sequence_for< T... >(), [ &u, &x, &phi ] ( auto i ) { LocalGalerkinOperator::value( u, x, std::get< i >( phi ) ); } );
234 }
235
236 template< class Basis, class Point >
237 static void values ( const Basis &basis, const Point &x, std::vector< typename Basis::RangeType > &phi )
238 {
239 basis.evaluateAll( x, phi );
240 }
241
242 template< class Basis, class Point >
243 static void values ( const Basis &basis, const Point &x, std::vector< typename Basis::JacobianRangeType > &phi )
244 {
245 basis.jacobianAll( x, phi );
246 }
247
248 template< class Basis, class Point >
249 static void values ( const Basis &basis, const Point &x, std::vector< typename Basis::HessianRangeType > &phi )
250 {
251 basis.hessianAll( x, phi );
252 }
253
254 template< class Basis, class Point, class... T >
255 static void values ( const Basis &basis, const Point &x, std::tuple< std::vector< T >... > &phi )
256 {
257 Hybrid::forEach( std::index_sequence_for< T... >(), [ &basis, &x, &phi ] ( auto i ) { LocalGalerkinOperator::values( basis, x, std::get< i >( phi ) ); } );
258 }
259
260 template< class LocalFunction, class Point >
261 static DomainValueType domainValue ( const LocalFunction &u, const Point &x )
262 {
263 DomainValueType phi;
264 value( u, x, phi );
265 return phi;
266 }
267
268 static DomainValueType domainValue ( const unsigned int qpIdx, DomainValueVectorType& vec)
269 {
270 DomainValueType phi;
271 Hybrid::forEach( DomainValueIndices(), [ &qpIdx, &vec, &phi ] ( auto i ) {
272 std::get< i > ( phi ) = std::get< i >( vec )[ qpIdx ];
273 } );
274 return phi;
275 }
276
277 template< class LocalFunction, class Quadrature >
278 static void domainValue ( const LocalFunction &u, const Quadrature& quadrature, DomainValueVectorType &result )
279 {
280 Hybrid::forEach( DomainValueIndices(), [ &u, &quadrature, &result ] ( auto i ) {
281 auto& vec = std::get< i >( result );
282 vec.resize( quadrature.nop() );
283 ThisType::evaluateQuadrature( u, quadrature, vec );
284 } );
285 }
286
287 template< class Phi, std::size_t... i >
288 static auto value ( const Phi &phi, std::size_t col, std::index_sequence< i... > )
289 {
290 return std::make_tuple( std::get< i >( phi )[ col ]... );
291 }
292
293 template< class... T >
294 static auto value ( const std::tuple< std::vector< T >... > &phi, std::size_t col )
295 {
296 return value( phi, col, std::index_sequence_for< T... >() );
297 }
298
299 static void assignRange( RangeValueVectorType& ranges, const std::size_t idx, const RangeValueType& range )
300 {
301 Hybrid::forEach( RangeValueIndices(), [ &ranges, &idx, &range ] ( auto i ) {
302 std::get< i >( ranges )[ idx ] = std::get< i >( range );
303 });
304 }
305 template <class W>
306 static void assignRange( RangeValueVectorType& ranges, const std::size_t idx, const RangeValueType& range, const W &weight )
307 {
308 Hybrid::forEach( RangeValueIndices(), [ &ranges, &idx, &range, &weight ] ( auto i ) {
309 std::get< i >( ranges )[ idx ] = std::get< i >( range );
310 std::get< i >( ranges )[ idx ] *= weight;
311 });
312 }
313
314 static void assignDomain( DomainValueVectorType& domains, const std::size_t idx, const DomainValueType& domain )
315 {
316 Hybrid::forEach( DomainValueIndices(), [ &domains, &idx, &domain ] ( auto i ) {
317 std::get< i >( domains )[ idx ] = std::get< i >( domain );
318 });
319 }
320
321 template <class W, class Quadrature>
322 static void axpyQuadrature( W& w, const Quadrature& quadrature, RangeValueVectorType& ranges )
323 {
324 Hybrid::forEach( RangeValueIndices(), [ &w, &quadrature, &ranges ] ( auto i ) {
325 w.axpyQuadrature( quadrature, std::get< i >( ranges ) );
326 } );
327 }
328
329 public:
330 // interior integral
331
332 template< class U, class W >
333 void addInteriorIntegral ( const U &u, W &w ) const
334 {
335 if( !integrands().init( u.entity() ) )
336 return;
337
338 const auto& geometry = u.geometry();
339
340 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: InteriorQuadratureType InteriorQuadratureType;
341 const InteriorQuadratureType quadrature( u.entity(), interiorQuadratureOrder(maxOrder(u, w)) );
342
343 // evaluate u for all quadrature points
344 DomainValueVectorType& domains = domainValues_;
345 domainValue( u, quadrature, domains );
346
347 auto& ranges = values_;
348 resizeRangeValueVector( ranges, quadrature.nop() );
349
350 // evaluate integrands for all quadrature points
351 for( const auto qp : quadrature )
352 {
353 const ctype weight = qp.weight() * geometry.integrationElement( qp.position() );
354 assignRange( ranges, qp.index(), integrands().interior( qp, domainValue( qp.index(), domains ) ), weight );
355 }
356
357 // add to w for all quadrature points
358 axpyQuadrature( w, quadrature, ranges );
359 integrands().unbind();
360 }
361
362 template< class U, class J >
363 void addLinearizedInteriorIntegral ( const U &u, J &j ) const
364 {
365 if( !integrands().init( u.entity() ) )
366 return;
367
368 const auto &geometry = u.geometry();
369 const auto &domainBasis = j.domainBasisFunctionSet();
370 const auto &rangeBasis = j.rangeBasisFunctionSet();
371
372 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: InteriorQuadratureType InteriorQuadratureType;
373 const InteriorQuadratureType quadrature( u.entity(), interiorQuadratureOrder( maxOrder( u, domainBasis, rangeBasis )) );
374 const size_t domainSize = domainBasis.size();
375 const size_t quadNop = quadrature.nop();
376
377 auto& basisValues = basisValues_;
378 resizeDomainValueVector( basisValues, domainSize );
379
380 // evaluate u for all quadrature points
381 auto& rangeValues = rangeValues_;
382 DomainValueVectorType& domains = domainValues_;
383 domainValue( u, quadrature, domains );
384
385 rangeValues.resize( domainSize );
386 for( std::size_t col = 0; col < domainSize; ++col )
387 {
388 resizeRangeValueVector( rangeValues[ col ], quadNop );
389 }
390
391 // evaluate all basis functions and integrands
392 for( const auto qp : quadrature )
393 {
394 values( domainBasis, qp, basisValues );
395 const auto weight = qp.weight() * geometry.integrationElement( qp.position() );
396 auto integrand = integrands().linearizedInterior( qp, domainValue( qp.index(), domains ) );
397 for( std::size_t col = 0; col < domainSize; ++col )
398 {
399 assignRange( rangeValues[ col ], qp.index(), integrand( value( basisValues, col ) ), weight );
400 }
401 }
402
403 // add to local matrix for all quadrature points and basis functions
404 for( std::size_t col = 0; col < domainSize; ++col )
405 {
406 LocalMatrixColumn< J > jCol( j, col );
407 axpyQuadrature( jCol, quadrature, rangeValues[ col ] );
408 }
409 integrands().unbind();
410 }
411
412 // boundary integral
413
414 template< class Intersection, class U, class W >
415 void addBoundaryIntegral ( const Intersection &intersection, const U &u, W &w ) const
416 {
417 if( !integrands().init( intersection ) )
418 return;
419
420 const auto geometry = intersection.geometry();
421 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
422 const SurfaceQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( u, w )), SurfaceQuadratureType::INSIDE );
423 for( const auto qp : quadrature )
424 {
425 const ctype weight = qp.weight() * geometry.integrationElement( qp.localPosition() );
426
427 RangeValueType integrand = integrands().boundary( qp, domainValue( u, qp ) );
428
429 Hybrid::forEach( RangeValueIndices(), [ &qp, &w, &integrand, weight ] ( auto i ) {
430 std::get< i >( integrand ) *= weight;
431 w.axpy( qp, std::get< i >( integrand ) );
432 } );
433 }
434 integrands().unbind();
435 }
436
437 template< class Intersection, class U, class J >
438 void addLinearizedBoundaryIntegral ( const Intersection &intersection, const U &u, J &j ) const
439 {
440 if( !integrands().init( intersection ) )
441 return;
442
443 DomainValueVectorType &phi = phiIn_;
444
445 const auto geometry = intersection.geometry();
446 const auto &domainBasis = j.domainBasisFunctionSet();
447 const auto &rangeBasis = j.rangeBasisFunctionSet();
448
449 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
450 const SurfaceQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder(u, domainBasis, rangeBasis )), SurfaceQuadratureType::INSIDE );
451 for( const auto qp : quadrature )
452 {
453 const ctype weight = qp.weight() * geometry.integrationElement( qp.localPosition() );
454
455 values( domainBasis, qp, phi );
456 auto integrand = integrands().linearizedBoundary( qp, domainValue( u, qp ) );
457
458 for( std::size_t col = 0, cols = domainBasis.size(); col < cols; ++col )
459 {
460 LocalMatrixColumn< J > jCol( j, col );
461 RangeValueType intPhi = integrand( value( phi, col ) );
462
463 Hybrid::forEach( RangeValueIndices(), [ &qp, &jCol, &intPhi, weight ] ( auto i ) {
464 std::get< i >( intPhi ) *= weight;
465 jCol.axpy( qp, std::get< i >( intPhi ) );
466 } );
467 }
468 }
469 integrands().unbind();
470 }
471
472 // addSkeletonIntegral
473
474 protected:
475 template< bool conforming, class Intersection, class U, class W >
476 void addSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut, W &wIn ) const
477 {
478 const auto geometry = intersection.geometry();
479
480 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
481 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
482 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( uIn, uOut, wIn)), false );
483
484 DomainValueVectorType& domainsIn = domainValues_;
485 domainValue( uIn, quadrature.inside(), domainsIn );
486
487 DomainValueVectorType& domainsOut = domainValuesOut_;
488 domainValue( uOut, quadrature.outside(), domainsOut );
489
490 for( unsigned int qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
491 {
492 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
493
494 const auto qpIn = quadrature.inside()[ qp ];
495 const auto qpOut = quadrature.outside()[ qp ];
496 std::pair< RangeValueType, RangeValueType > integrand
497 = integrands().skeleton( qpIn, domainValue( qp, domainsIn ), qpOut, domainValue( qp, domainsOut ) );
498
499 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &wIn, &integrand, weight ] ( auto i ) {
500 std::get< i >( integrand.first ) *= weight;
501 wIn.axpy( qpIn, std::get< i >( integrand.first ) );
502 } );
503 }
504 }
505
506 template< bool conforming, class Intersection, class U, class W >
507 void addSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut, W &wIn, W &wOut ) const
508 {
509 const auto geometry = intersection.geometry();
510 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
511 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
512 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( uIn, uOut, wIn, wOut)), false );
513
514 DomainValueVectorType& domainsIn = domainValues_;
515 domainValue( uIn, quadrature.inside(), domainsIn );
516
517 DomainValueVectorType& domainsOut = domainValuesOut_;
518 domainValue( uOut, quadrature.outside(), domainsOut );
519
520 for( unsigned int qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
521 {
522 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
523
524 const auto qpIn = quadrature.inside()[ qp ];
525 const auto qpOut = quadrature.outside()[ qp ];
526 std::pair< RangeValueType, RangeValueType > integrand
527 = integrands().skeleton( qpIn, domainValue( qp, domainsIn ), qpOut, domainValue( qp, domainsOut ) );
528
529 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &wIn, &qpOut, &wOut, &integrand, weight ] ( auto i ) {
530 std::get< i >( integrand.first ) *= weight;
531 wIn.axpy( qpIn, std::get< i >( integrand.first ) );
532
533 std::get< i >( integrand.second ) *= weight;
534 wOut.axpy( qpOut, std::get< i >( integrand.second ) );
535 } );
536 }
537 }
538
539 template< bool conforming, class Intersection, class U, class J >
540 void addLinearizedSkeletonIntegral ( const Intersection &intersection,
541 const U &uIn, const U &uOut, J &jInIn, J &jOutIn ) const
542 {
543 DomainValueVectorType &phiIn = phiIn_;
544 DomainValueVectorType &phiOut = phiOut_;
545
546 const auto &domainBasisIn = jInIn.domainBasisFunctionSet();
547 const auto &domainBasisOut = jOutIn.domainBasisFunctionSet();
548
549 const auto &rangeBasisIn = jInIn.rangeBasisFunctionSet();
550
551 const int order = std::max( maxOrder(uIn, uOut), maxOrder( domainBasisIn, domainBasisOut, rangeBasisIn ));
552
553 const auto geometry = intersection.geometry();
554 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
555 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
556 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(order), false );
557
558 DomainValueVectorType& domainsIn = domainValues_;
559 domainValue( uIn, quadrature.inside(), domainsIn );
560
561 DomainValueVectorType& domainsOut = domainValuesOut_;
562 domainValue( uOut, quadrature.outside(), domainsOut );
563
564 for( unsigned int qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
565 {
566 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
567
568 const auto qpIn = quadrature.inside()[ qp ];
569 const auto qpOut = quadrature.outside()[ qp ];
570
571 values( domainBasisIn, qpIn, phiIn );
572 values( domainBasisOut, qpOut, phiOut );
573
574 auto integrand = integrands().linearizedSkeleton( qpIn, domainValue( qp, domainsIn ), qpOut, domainValue( qp, domainsOut ) );
575 for( std::size_t col = 0, cols = domainBasisIn.size(); col < cols; ++col )
576 {
577 LocalMatrixColumn< J > jInInCol( jInIn, col );
578 std::pair< RangeValueType, RangeValueType > intPhi = integrand.first( value( phiIn, col ) );
579
580 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jInInCol, &intPhi, weight ] ( auto i ) {
581 std::get< i >( intPhi.first ) *= weight;
582 jInInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
583 } );
584 }
585 for( std::size_t col = 0, cols = domainBasisOut.size(); col < cols; ++col )
586 {
587 LocalMatrixColumn< J > jOutInCol( jOutIn, col );
588 std::pair< RangeValueType, RangeValueType > intPhi = integrand.second( value( phiOut, col ) );
589
590 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jOutInCol, &intPhi, weight ] ( auto i ) {
591 std::get< i >( intPhi.first ) *= weight;
592 jOutInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
593 } );
594 }
595 }
596 }
597
598 template< bool conforming, class Intersection, class U, class J >
599 void addLinearizedSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut,
600 J &jInIn, J &jOutIn, J &jInOut, J &jOutOut ) const
601 {
602 DomainValueVectorType &phiIn = phiIn_;
603 DomainValueVectorType &phiOut = phiOut_;
604
605 const auto &domainBasisIn = jInIn.domainBasisFunctionSet();
606 const auto &domainBasisOut = jOutIn.domainBasisFunctionSet();
607
608 const auto &rangeBasisIn = jInIn.rangeBasisFunctionSet();
609 const auto &rangeBasisOut = jInOut.rangeBasisFunctionSet();
610
611 const int order = std::max( maxOrder(uIn, uOut), maxOrder( domainBasisIn, domainBasisOut, rangeBasisIn, rangeBasisOut ));
612
613 const auto geometry = intersection.geometry();
614 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
615 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
616 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(order), false );
617
618 DomainValueVectorType& domainsIn = domainValues_;
619 domainValue( uIn, quadrature.inside(), domainsIn );
620
621 DomainValueVectorType& domainsOut = domainValuesOut_;
622 domainValue( uOut, quadrature.outside(), domainsOut );
623
624 for( unsigned int qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
625 {
626 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
627
628 const auto qpIn = quadrature.inside()[ qp ];
629 const auto qpOut = quadrature.outside()[ qp ];
630
631 values( domainBasisIn, qpIn, phiIn );
632 values( domainBasisOut, qpOut, phiOut );
633
634 auto integrand = integrands().linearizedSkeleton( qpIn, domainValue( qp, domainsIn ), qpOut, domainValue( qp, domainsOut ) );
635 for( std::size_t col = 0, cols = domainBasisIn.size(); col < cols; ++col )
636 {
637 LocalMatrixColumn< J > jInInCol( jInIn, col );
638 LocalMatrixColumn< J > jInOutCol( jInOut, col );
639 std::pair< RangeValueType, RangeValueType > intPhi = integrand.first( value( phiIn, col ) );
640
641 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jInInCol, &qpOut, &jInOutCol, &intPhi, weight ] ( auto i ) {
642 std::get< i >( intPhi.first ) *= weight;
643 jInInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
644
645 std::get< i >( intPhi.second ) *= weight;
646 jInOutCol.axpy( qpOut, std::get< i >( intPhi.second ) );
647 } );
648 }
649 for( std::size_t col = 0, cols = domainBasisOut.size(); col < cols; ++col )
650 {
651 LocalMatrixColumn< J > jOutInCol( jOutIn, col );
652 LocalMatrixColumn< J > jOutOutCol( jOutOut, col );
653 std::pair< RangeValueType, RangeValueType > intPhi = integrand.second( value( phiOut, col ) );
654
655 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jOutInCol, &qpOut, &jOutOutCol, &intPhi, weight ] ( auto i ) {
656 std::get< i >( intPhi.first ) *= weight;
657 jOutInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
658
659 std::get< i >( intPhi.second ) *= weight;
660 jOutOutCol.axpy( qpOut, std::get< i >( intPhi.second ) );
661 } );
662 }
663 }
664 }
665
666 public:
667 template< class Intersection, class U, class... W >
668 void addSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut, W &... w ) const
669 {
670 if( !integrands().init( intersection ) )
671 return;
672
673 if( intersection.conforming() )
674 addSkeletonIntegral< true >( intersection, uIn, uOut, w... );
675 else
676 addSkeletonIntegral< false >( intersection, uIn, uOut, w... );
677 integrands().unbind();
678 }
679
680 template< class Intersection, class U, class... J >
681 void addLinearizedSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut, J &... j ) const
682 {
683 if( !integrands().init( intersection ) )
684 return;
685
686 if( intersection.conforming() )
687 addLinearizedSkeletonIntegral< true >( intersection, uIn, uOut, j... );
688 else
689 addLinearizedSkeletonIntegral< false >( intersection, uIn, uOut, j... );
690 integrands().unbind();
691 }
692
693 void setQuadratureOrders(unsigned int interior, unsigned int surface)
694 {
695 interiorQuadOrder_ = interior;
696 surfaceQuadOrder_ = surface;
697 }
698
699 void setQuadratureOrderFunctions( std::function<int(const int)> interiorOrder,
700 std::function<int(const int)> surfaceOrder )
701 {
702 defaultInteriorOrder_ = interiorOrder;
703 defaultSurfaceOrder_ = surfaceOrder;
704 }
705
706 IntegrandsType& model() const
707 {
708 return integrands();
709 }
710 bool nonlinear() const { return model().nonlinear(); }
711 bool hasInterior() const { return model().hasInterior(); }
712 bool hasSkeleton() const { return model().hasSkeleton(); }
713 bool hasBoundary() const { return model().hasBoundary(); }
714
715 private:
716 IntegrandsType& integrands() const
717 {
718 return integrands_;
719 }
720
721 public:
722 // accessors
723 const GridPartType &gridPart () const { return gridPart_; }
724
725 unsigned int interiorQuadratureOrder(unsigned int order) const { return interiorQuadOrder_ == 0 ? defaultInteriorOrder_(order) : interiorQuadOrder_; }
726 unsigned int surfaceQuadratureOrder(unsigned int order) const { return surfaceQuadOrder_ == 0 ? defaultSurfaceOrder_ (order) : surfaceQuadOrder_; }
727
728 protected:
729 template <class U>
730 int maxOrder( const U& u ) const
731 {
732 return CallOrder< U > :: order( u );
733 }
734
735 template< class U, class W >
736 int maxOrder( const U& u, const W& w ) const
737 {
738 return std::max( maxOrder( u ), maxOrder( w ) );
739 }
740
741 template< class U, class V, class W >
742 int maxOrder( const U& u, const V& v, const W& w ) const
743 {
744 return std::max( maxOrder( u, v ), maxOrder( w ) );
745 }
746
747 template< class U, class V, class W, class X >
748 int maxOrder( const U& u, const V& v, const W& w, const X& x ) const
749 {
750 return std::max( maxOrder( u, v ), maxOrder( w, x) );
751 }
752
753 protected:
754 const GridPartType &gridPart_;
755
756 mutable IntegrandsType integrands_;
757
758 mutable std::function<int(const int)> defaultInteriorOrder_;
759 mutable std::function<int(const int)> defaultSurfaceOrder_;
760
761 unsigned int interiorQuadOrder_;
762 unsigned int surfaceQuadOrder_;
763
764 mutable std::vector< RangeValueVectorType > rangeValues_;
765 mutable RangeValueVectorType values_;
766 mutable DomainValueVectorType phiIn_;
767 mutable DomainValueVectorType phiOut_;
768 mutable DomainValueVectorType basisValues_;
769 mutable DomainValueVectorType domainValues_;
770 mutable DomainValueVectorType domainValuesOut_;
771 };
772
773
774
775 // GalerkinOperator
776 // ----------------
777
778 template< class GridPart >
779 // Integrands, template <class> class QuadSelector = DefaultGalerkinOperatorQuadratureSelector >
780 struct GalerkinOperator
781 {
782 typedef GridPart GridPartType;
783 typedef GalerkinOperator< GridPartType > ThisType;
784
785 typedef typename GridPartType::ctype ctype;
786 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
787
788 // constructor
789 explicit GalerkinOperator ( const GridPartType &gridPart )
790 : gridPart_( gridPart ),
791 gridSizeInterior_( 0 )
792 {
793 }
794
795 protected:
796 template <class IntegrandsTuple>
797 bool hasBoundary( const IntegrandsTuple& integrandsTuple ) const
798 {
799 typedef std::make_index_sequence< std::tuple_size< IntegrandsTuple >::value > Indices;
800 bool hasBoundary = false ;
801 Hybrid::forEach( Indices(), [&integrandsTuple, &hasBoundary]( auto i ) {
802 if( std::get< i > (integrandsTuple).hasBoundary() )
803 {
804 hasBoundary = true ;
805 return ;
806 }
807 });
808 return hasBoundary;
809 }
810
811 template< class GridFunction, class DiscreteFunction, class Iterators, class IntegrandsTuple, class Functor, bool hasSkeleton >
812 void evaluateImpl ( const GridFunction &u, DiscreteFunction &w, const Iterators& iterators,
813 const IntegrandsTuple& integrandsTuple, Functor& addLocalDofs, std::integral_constant<bool, hasSkeleton> ) const
814 {
815 Dune::Fem::ConstLocalFunction< GridFunction > uInside( u );
816 Dune::Fem::ConstLocalFunction< GridFunction > uOutside( u );
817
818 typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
819 TemporaryLocalFunction< DiscreteFunctionSpaceType > wInside( w.space() ), wOutside( w.space() );
820
821 // element counter
822 gridSizeInterior_ = 0;
823
824 typedef std::make_index_sequence< std::tuple_size< IntegrandsTuple >::value > Indices;
825
826 // true if one of the integrands has a boundary term
827 const bool hasBnd = hasBoundary( integrandsTuple );
828
829 const auto &indexSet = gridPart().indexSet();
830 const auto end = iterators.end();
831 for( auto it = iterators.begin(); it != end; ++it )
832 {
833 // assert( iterators.thread( *it ) == MPIManager::thread() );
834 const EntityType inside = *it ;
835
836 // increase counter for interior elements
837 ++gridSizeInterior_;
838
839 auto uGuard = bindGuard( uInside, inside );
840 auto wGuard = bindGuard( wInside, inside );
841 wInside.clear();
842
843 auto addInteriorIntegral = [&integrandsTuple, &uInside, &wInside]( auto i )
844 {
845 const auto& integrands = std::get< i >( integrandsTuple );
846 if( integrands.hasInterior() )
847 integrands.addInteriorIntegral( uInside, wInside );
848 };
849 // add interior integral of any integrands
850 Hybrid::forEach( Indices(), addInteriorIntegral );
851
852 if( hasSkeleton || (hasBnd && HasBoundaryIntersection<GridPartType>::apply(inside) ) )
853 {
854 for( const auto &intersection : intersections( gridPart(), inside ) )
855 {
856 bool neighbor = false;
857 if constexpr ( hasSkeleton )
858 {
859 // check neighbor first since on periodic boundaries both,
860 // neighbor and boundary are true, so we treat neighbor first
861 if( intersection.neighbor() )
862 {
863 neighbor = true;
864 const EntityType outside = intersection.outside();
865
866 if( outside.partitionType() != InteriorEntity )
867 {
868 auto uOutGuard = bindGuard( uOutside, outside );
869
870 auto addSkeletonIntegral = [&integrandsTuple, &intersection, &uInside, &uOutside, &wInside] ( auto i )
871 {
872 const auto& integrands = std::get< i >( integrandsTuple );
873 if( integrands.hasSkeleton() )
874 integrands.addSkeletonIntegral( intersection, uInside, uOutside, wInside );
875 };
876 // add skeleton integral of any integrands
877 Hybrid::forEach( Indices(), addSkeletonIntegral );
878 }
879 else if( indexSet.index( inside ) < indexSet.index( outside ) )
880 {
881 auto uOutGuard = bindGuard( uOutside, outside );
882 auto wOutGuard = bindGuard( wOutside, outside );
883 wOutside.clear();
884
885 auto addSkeletonIntegral = [&integrandsTuple, &intersection, &uInside, &uOutside, &wInside, &wOutside] ( auto i )
886 {
887 const auto& integrands = std::get< i >( integrandsTuple );
888 if( integrands.hasSkeleton() )
889 integrands.addSkeletonIntegral( intersection, uInside, uOutside, wInside, wOutside );
890 };
891 // add skeleton integral of any integrands
892 Hybrid::forEach( Indices(), addSkeletonIntegral );
893
894 // addLocalDofs calls w.addLocalDofs but also
895 // prevents race condition for thread parallel runs
896 addLocalDofs( outside, wOutside );
897 }
898 }
899 } // end skeleton
900
901 if( ! neighbor && intersection.boundary() )
902 {
903 auto addBoundaryIntegral = [&integrandsTuple, &intersection, &uInside, &wInside]( auto i )
904 {
905 const auto& integrands = std::get< i >( integrandsTuple );
906 if( integrands.hasBoundary() )
907 integrands.addBoundaryIntegral( intersection, uInside, wInside );
908 };
909 // add boundary integral of any integrands
910 Hybrid::forEach( Indices(), addBoundaryIntegral );
911 } // end boundary
912 }
913 } // end intersections
914
915 addLocalDofs( inside, wInside );
916 }
917 }
918
919 template <class Space>
920 struct InsideEntity
921 {
922 typedef typename Space::EntityType EntityType;
923 template <class Iterators>
924 InsideEntity(const Space &space, const Iterators& iterators)
925 : space_(space), dofThread_(space.size(),-1)
926 , thread_( iterators.thread() )
927 {
928 const auto& mapper = space_.blockMapper();
929 for (const auto &entity : space_)
930 {
931 int t=iterators.threadParallel(entity);
932 mapper.mapEach(entity, [ this, t ] ( int local, auto global )
933 { dofThread_[global] = (dofThread_[global]==t || dofThread_[global]==-1)?
934 t : -2 ; } ); // -2: shared dof
935 }
936 }
937 bool operator()(const EntityType &entity) const
938 {
939 bool needsLocking = false;
940 space_.blockMapper().mapEach(entity,
941 [ this, &needsLocking ] ( int local, auto global )
942 { needsLocking = (needsLocking || dofThread_[global]!=thread_); });
943 return !needsLocking;
944 }
945 const Space &space_;
946 std::vector<int> dofThread_;
947 int thread_;
948 };
949
950 template <class DiscreteFunction>
951 struct AddLocalEvaluate
952 {
953 AddLocalEvaluate(DiscreteFunction &w)
954 : w_(w) {}
955 template <class LocalDofs>
956 void operator () (const EntityType& entity, const LocalDofs& wLocal ) const
957 {
958 w_.addLocalDofs( entity, wLocal.localDofVector() );
959 }
960 DiscreteFunction &w_;
961 };
962
963 template <class DiscreteFunction>
964 struct AddLocalEvaluateLocked : public AddLocalEvaluate<DiscreteFunction>
965 {
966 typedef AddLocalEvaluate<DiscreteFunction> BaseType;
967
968 std::shared_mutex& mutex_;
969 InsideEntity<typename DiscreteFunction::DiscreteFunctionSpaceType> inside_;
970
971 template <class Iterators>
972 AddLocalEvaluateLocked(DiscreteFunction &w, std::shared_mutex& mtx, const Iterators &iterators)
973 : BaseType(w), mutex_(mtx), inside_(w.space(),iterators) {}
974
975 template <class LocalDofs>
976 void operator () (const EntityType& entity, const LocalDofs& wLocal ) const
977 {
978 // call addLocalDofs on w
979 if (inside_(entity))
980 {
981 std::shared_lock<std::shared_mutex> guard ( mutex_ );
982 BaseType::operator()( entity, wLocal );
983 }
984 else
985 {
986 // lock mutex (unlock on destruction)
987 std::lock_guard<std::shared_mutex> guard ( mutex_ );
988 BaseType::operator()( entity, wLocal );
989 }
990 }
991 };
992
993 template< class GridFunction, class DiscreteFunction, class Iterators, class IntegrandsTuple, class Functor >
994 void evaluate ( const GridFunction &u, DiscreteFunction &w, const Iterators& iterators,
995 const IntegrandsTuple& integrandsTuple, Functor& addLocalDofs ) const
996 {
997 static_assert( std::is_same< typename GridFunction::GridPartType, GridPartType >::value, "Argument 'u' and Integrands must be defined on the same grid part." );
998 static_assert( std::is_same< typename DiscreteFunction::GridPartType, GridPartType >::value, "Argument 'w' and Integrands must be defined on the same grid part." );
999
1000 if( hasSkeleton( integrandsTuple ) )
1001 evaluateImpl( u, w, iterators, integrandsTuple, addLocalDofs, std::true_type() );
1002 else
1003 evaluateImpl( u, w, iterators, integrandsTuple, addLocalDofs, std::false_type() );
1004 }
1005
1006 public:
1007 template <class IntegrandsTuple>
1008 bool hasSkeleton( const IntegrandsTuple& integrandsTuple ) const
1009 {
1010 typedef std::make_index_sequence< std::tuple_size< IntegrandsTuple >::value > Indices;
1011 bool hasSkeleton = false ;
1012 Hybrid::forEach( Indices(), [&integrandsTuple, &hasSkeleton] ( auto i ) {
1013 if( std::get< i >( integrandsTuple ).hasSkeleton() )
1014 {
1015 hasSkeleton = true;
1016 return ;
1017 }
1018 });
1019 return hasSkeleton ;
1020 }
1021
1022 template< class GridFunction, class DiscreteFunction, class Iterators, class IntegrandsTuple >
1023 void evaluate ( const GridFunction &u, DiscreteFunction &w, const Iterators& iterators,
1024 const IntegrandsTuple& integrandsTuple, std::shared_mutex& mtx ) const
1025 {
1026 AddLocalEvaluateLocked<DiscreteFunction> addLocalEvaluate(w,mtx,iterators);
1027 evaluate( u, w, iterators, integrandsTuple, addLocalEvaluate );
1028 }
1029
1030 template< class GridFunction, class DiscreteFunction, class Iterators, class IntegrandsTuple >
1031 void evaluate ( const GridFunction &u, DiscreteFunction &w, const Iterators& iterators, const IntegrandsTuple& integrandsTuple ) const
1032 {
1033 AddLocalEvaluate<DiscreteFunction> addLocalEvaluate(w);
1034 evaluate( u, w, iterators, integrandsTuple, addLocalEvaluate );
1035 }
1036
1037 protected:
1038 template<class T, int length>
1039 class FiniteStack
1040 {
1041 public :
1042 // Makes empty stack
1043 FiniteStack () : _f(0) {}
1044
1045 // Returns true if the stack is empty
1046 bool empty () const { return _f <= 0; }
1047
1048 // Returns true if the stack is full
1049 bool full () const { return (_f >= length); }
1050
1051 // clear stack
1052 void clear() { _f = 0; }
1053
1054 // Puts a new object onto the stack
1055 void push (const T& t)
1056 {
1057 assert ( _f < length );
1058 _s[_f++] = t;
1059 }
1060
1061 // Removes and returns the uppermost object from the stack
1062 T pop () {
1063 assert ( _f > 0 );
1064 return _s[--_f];
1065 }
1066
1067 // Returns the uppermost object on the stack
1068 T top () const {
1069 assert ( _f > 0 );
1070 return _s[_f-1];
1071 }
1072
1073 // stacksize
1074 int size () const { return _f; }
1075
1076 private:
1077 T _s[length]; // the stack
1078 int _f; // actual position in stack
1079 };
1080
1081
1082 template <class JacobianOperator>
1083 struct AddLocalAssemble
1084 {
1085 typedef typename JacobianOperator::DomainSpaceType DomainSpaceType;
1086 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
1087 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
1088 JacobianOperator &jOp_;
1089 std::vector< TemporaryLocalMatrixType > jOpLocal_;
1090
1091 FiniteStack< TemporaryLocalMatrixType*, 12 > jOpLocalFinalized_;
1092 FiniteStack< TemporaryLocalMatrixType*, 12 > jOpLocalFree_;
1093
1094 std::size_t locked, notLocked, timesLocked;
1095 AddLocalAssemble(JacobianOperator& jOp)
1096 : jOp_(jOp)
1097 , jOpLocal_(12, TemporaryLocalMatrixType(jOp_.domainSpace(), jOp_.rangeSpace()))
1098 , jOpLocalFinalized_()
1099 , jOpLocalFree_()
1100 , locked(0), notLocked(0), timesLocked(0)
1101 {
1102 for( auto& jOpLocal : jOpLocal_ )
1103 jOpLocalFree_.push( &jOpLocal );
1104 }
1105
1106 TemporaryLocalMatrixType& bind(const EntityType& dE, const EntityType& rE)
1107 {
1108 assert( ! jOpLocalFree_.empty() );
1109 TemporaryLocalMatrixType& lop = *(jOpLocalFree_.pop());
1110 lop.bind(dE,rE);
1111 lop.clear();
1112 return lop;
1113 }
1114
1115 void unbind(TemporaryLocalMatrixType &lop)
1116 {
1117 notLocked += 1;
1118 jOp_.addLocalMatrix( lop.domainEntity(), lop.rangeEntity(), lop );
1119 lop.unbind();
1120 jOpLocalFree_.push( &lop );
1121 }
1122
1123 void finalize()
1124 {
1125 locked += jOpLocalFinalized_.size();
1126 while ( ! jOpLocalFinalized_.empty() )
1127 {
1128 TemporaryLocalMatrixType &lop = *(jOpLocalFinalized_.pop());
1129 jOp_.addLocalMatrix( lop.domainEntity(), lop.rangeEntity(), lop );
1130 lop.unbind();
1131 jOpLocalFree_.push( &lop );
1132 }
1133 }
1134 };
1135
1136 template <class JacobianOperator>
1137 struct AddLocalAssembleLocked : public AddLocalAssemble<JacobianOperator>
1138 {
1139 typedef AddLocalAssemble<JacobianOperator> BaseType;
1140 typedef typename BaseType::TemporaryLocalMatrixType TemporaryLocalMatrixType;
1141 using BaseType::jOpLocalFinalized_;
1142 using BaseType::jOpLocalFree_;
1143
1144 std::shared_mutex& mutex_;
1145 InsideEntity<typename JacobianOperator::DomainSpaceType> insideDomain_;
1146 InsideEntity<typename JacobianOperator::RangeSpaceType> insideRange_;
1147
1148 template <class Iterators>
1149 AddLocalAssembleLocked(JacobianOperator &jOp, std::shared_mutex &mtx, const Iterators &iterators)
1150 : BaseType(jOp)
1151 , mutex_(mtx)
1152 , insideDomain_(jOp.domainSpace(),iterators)
1153 , insideRange_(jOp.rangeSpace(),iterators)
1154 {}
1155
1156 void finalize()
1157 {
1158 // lock mutex (unlock on destruction)
1159 ++BaseType::timesLocked;
1160 std::lock_guard<std::shared_mutex> guard ( mutex_ );
1161 BaseType::finalize();
1162 }
1163
1164 TemporaryLocalMatrixType& bind(const EntityType& dE, const EntityType& rE)
1165 {
1166 if ( jOpLocalFree_.empty() )
1167 {
1168 finalize();
1169 }
1170 return BaseType::bind(dE,rE);
1171 }
1172
1173 void unbind(TemporaryLocalMatrixType &lop)
1174 {
1175 /* // always lock
1176 ++BaseType::timesLocked;
1177 ++BaseType::locked;
1178 std::lock_guard guard ( mutex_ );
1179 BaseType::unbind(lop);
1180 return;
1181 */
1182 if ( insideDomain_(lop.domainEntity()) &&
1183 insideRange_(lop.rangeEntity()) )
1184 {
1185 std::shared_lock<std::shared_mutex> guard ( mutex_ );
1186 BaseType::unbind(lop);
1187 }
1188 else
1189 {
1190 jOpLocalFinalized_.push( &lop );
1191 }
1192 }
1193 };
1194
1195 template< class GridFunction, class JacobianOperator, class Iterators, class IntegrandsTuple, class Functor, bool hasSkeleton >
1196 void assembleImpl ( const GridFunction &u, JacobianOperator &jOp, const Iterators& iterators, const IntegrandsTuple& integrandsTuple,
1197 Functor& addLocalMatrix, std::integral_constant<bool, hasSkeleton> ) const
1198 {
1199 typedef typename JacobianOperator::DomainSpaceType DomainSpaceType;
1200 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
1201
1202 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
1203
1204 Dune::Fem::ConstLocalFunction< GridFunction > uIn( u );
1205 Dune::Fem::ConstLocalFunction< GridFunction > uOut( u );
1206
1207 typedef std::make_index_sequence< std::tuple_size< IntegrandsTuple >::value > Indices;
1208 const std::size_t maxNumLocalDofs = jOp.domainSpace().blockMapper().maxNumDofs() * jOp.domainSpace().localBlockSize;
1209
1210 // initialize local temporary data
1211 Hybrid::forEach( Indices(), [&integrandsTuple, &maxNumLocalDofs] ( auto i ) {
1212 const auto& integrands = std::get< i >( integrandsTuple );
1213 integrands.prepare( maxNumLocalDofs );
1214 });
1215
1216 // element counter
1217 gridSizeInterior_ = 0;
1218
1219 // true if one of the integrands has a boundary term
1220 const bool hasBnd = hasBoundary( integrandsTuple );
1221
1222 const auto &indexSet = gridPart().indexSet();
1223 // threaded iterators provide from outside
1224 const auto end = iterators.end();
1225 for( auto it = iterators.begin(); it != end; ++it )
1226 {
1227 // increase counter for interior elements
1228 ++gridSizeInterior_;
1229
1230 const EntityType inside = *it;
1231
1232 auto uiGuard = bindGuard( uIn, inside );
1233
1234 TemporaryLocalMatrixType& jOpInIn = addLocalMatrix.bind( inside, inside );
1235 auto addLinearizedInteriorIntegral = [&integrandsTuple, &uIn, &jOpInIn]( auto i )
1236 {
1237 const auto& integrands = std::get< i >( integrandsTuple );
1238 if( integrands.hasInterior() )
1239 integrands.addLinearizedInteriorIntegral( uIn, jOpInIn );
1240 };
1241 // add interior integral of any integrands
1242 Hybrid::forEach( Indices(), addLinearizedInteriorIntegral );
1243
1244 if( hasSkeleton || (hasBnd && HasBoundaryIntersection<GridPartType>::apply(inside) ) )
1245 {
1246 for( const auto &intersection : intersections( gridPart(), inside ) )
1247 {
1248 bool neighbor = false ;
1249 // check neighbor first since on periodic boundaries both,
1250 // neighbor and boundary are true, so we treat neighbor first
1251 if constexpr ( hasSkeleton )
1252 {
1253 if( intersection.neighbor() )
1254 {
1255 neighbor = true ;
1256 const EntityType &outside = intersection.outside();
1257
1258 TemporaryLocalMatrixType &jOpOutIn = addLocalMatrix.bind( outside, inside );
1259
1260 auto uoGuard = bindGuard( uOut, outside );
1261
1262 if( outside.partitionType() != InteriorEntity )
1263 {
1264 auto addLinearizedSkeletonIntegral = [&integrandsTuple, &intersection, &uIn, &uOut, &jOpInIn, &jOpOutIn]( auto i )
1265 {
1266 const auto& integrands = std::get< i >( integrandsTuple );
1267 if( integrands.hasSkeleton() )
1268 integrands.addLinearizedSkeletonIntegral( intersection, uIn, uOut, jOpInIn, jOpOutIn );
1269 };
1270 // add skeleton integral of any integrands
1271 Hybrid::forEach( Indices(), addLinearizedSkeletonIntegral );
1272 }
1273 else if( indexSet.index( inside ) < indexSet.index( outside ) )
1274 {
1275 TemporaryLocalMatrixType &jOpInOut = addLocalMatrix.bind( inside, outside );
1276 TemporaryLocalMatrixType &jOpOutOut = addLocalMatrix.bind( outside, outside );
1277
1278 auto addLinearizedSkeletonIntegral = [&integrandsTuple, &intersection, &uIn, &uOut, &jOpInIn, &jOpOutIn, &jOpInOut, &jOpOutOut]( auto i )
1279 {
1280 const auto& integrands = std::get< i >( integrandsTuple );
1281 if( integrands.hasSkeleton() )
1282 integrands.addLinearizedSkeletonIntegral( intersection, uIn, uOut, jOpInIn, jOpOutIn, jOpInOut, jOpOutOut );
1283 };
1284 // add skeleton integral of any integrands
1285 Hybrid::forEach( Indices(), addLinearizedSkeletonIntegral );
1286
1287 addLocalMatrix.unbind(jOpInOut);
1288 addLocalMatrix.unbind(jOpOutOut);
1289 }
1290
1291 addLocalMatrix.unbind(jOpOutIn);
1292 }
1293 } // end skeleton
1294
1295 if( !neighbor && intersection.boundary() )
1296 {
1297 auto addLinearizedBoundaryIntegral = [&integrandsTuple, &intersection, &uIn, &jOpInIn]( auto i )
1298 {
1299 const auto& integrands = std::get< i >( integrandsTuple );
1300 if( integrands.hasBoundary() )
1301 integrands.addLinearizedBoundaryIntegral( intersection, uIn, jOpInIn );
1302 };
1303 // add skeleton integral of any integrands
1304 Hybrid::forEach( Indices(), addLinearizedBoundaryIntegral );
1305
1306 } // end boundary
1307 }
1308 } // end intersection
1309 addLocalMatrix.unbind(jOpInIn);
1310 }
1311
1312 // complete the matrix build
1313 addLocalMatrix.finalize();
1314 }
1315
1316
1317 template< class GridFunction, class JacobianOperator, class Iterators, class IntegrandsTuple, class Functor >
1318 void assemble ( const GridFunction &u, JacobianOperator &jOp, const Iterators& iterators,
1319 const IntegrandsTuple& integrandsTuple, Functor& addLocalMatrix, int ) const
1320 {
1321 static_assert( std::is_same< typename GridFunction::GridPartType, GridPartType >::value, "Argument 'u' and Integrands must be defined on the same grid part." );
1322 static_assert( std::is_same< typename JacobianOperator::DomainSpaceType::GridPartType, GridPartType >::value, "Argument 'jOp' and Integrands must be defined on the same grid part." );
1323 static_assert( std::is_same< typename JacobianOperator::RangeSpaceType::GridPartType, GridPartType >::value, "Argument 'jOp' and Integrands must be defined on the same grid part." );
1324
1325 if( hasSkeleton( integrandsTuple ) )
1326 assembleImpl( u, jOp, iterators, integrandsTuple ,addLocalMatrix, std::true_type() );
1327 else
1328 assembleImpl( u, jOp, iterators, integrandsTuple, addLocalMatrix, std::false_type() );
1329 }
1330
1331 public:
1332 template< class GridFunction, class JacobianOperator, class Iterators, class IntegrandsTuple>
1333 void assemble ( const GridFunction &u, JacobianOperator &jOp, const Iterators& iterators,
1334 const IntegrandsTuple& integrandsTuple, std::shared_mutex& mtx) const
1335 {
1336 AddLocalAssembleLocked<JacobianOperator> addLocalAssemble( jOp, mtx, iterators);
1337 assemble( u, jOp, iterators, integrandsTuple, addLocalAssemble, 10 );
1338 #if 0 // print information about how many times a lock was used during assemble
1339 std::lock_guard guard ( mtx );
1340 std::cout << MPIManager::thread() << " : "
1341 << addLocalAssemble.locked << " " << addLocalAssemble.notLocked << " "
1342 << addLocalAssemble.timesLocked << std::endl;
1343 #endif
1344 }
1345
1346 template< class GridFunction, class JacobianOperator, class Iterators, class IntegrandsTuple>
1347 void assemble ( const GridFunction &u, JacobianOperator &jOp, const Iterators& iterators, const IntegrandsTuple& integrandsTuple ) const
1348 {
1349 AddLocalAssemble<JacobianOperator> addLocalAssemble(jOp);
1350 assemble( u, jOp, iterators, integrandsTuple, addLocalAssemble, 10 );
1351 }
1352
1353 // accessors
1354 const GridPartType &gridPart () const { return gridPart_; }
1355
1356 std::size_t gridSizeInterior () const { return gridSizeInterior_; }
1357
1358 protected:
1359 const GridPartType &gridPart_;
1360 mutable std::size_t gridSizeInterior_;
1361 };
1362
1363
1364 template <class GalerkinOperator >
1365 static std::size_t accumulateGridSize( const ThreadSafeValue< GalerkinOperator >& ops )
1366 {
1367 std::size_t s = ops.size();
1368 std::size_t sum = 0;
1369 for( std::size_t i=0; i<s; ++i )
1370 sum += ops[ i ].gridSizeInterior();
1371 return sum;
1372 }
1373
1374 } // namespace Impl
1375
1378
1379
1380 // GalerkinOperator
1381 // ----------------
1382
1383 template< class Integrands, class DomainFunction, class RangeFunction = DomainFunction >
1384 struct GalerkinOperator
1385 : public virtual Operator< DomainFunction, RangeFunction >
1386 {
1387 typedef DomainFunction DomainFunctionType;
1388 typedef RangeFunction RangeFunctionType;
1389
1390 typedef typename RangeFunctionType::GridPartType GridPartType;
1391
1392 typedef Impl::LocalGalerkinOperator< Integrands > LocalGalerkinOperatorImplType;
1393 typedef Impl::GalerkinOperator< GridPartType > GalerkinOperatorImplType;
1394
1395 static_assert( std::is_same< typename DomainFunctionType::GridPartType, typename RangeFunctionType::GridPartType >::value, "DomainFunction and RangeFunction must be defined on the same grid part." );
1396
1397 typedef ThreadIterator< GridPartType > ThreadIteratorType;
1398
1399 template< class... Args >
1400 explicit GalerkinOperator ( const GridPartType &gridPart, Args &&... args )
1401 : iterators_( gridPart ),
1402 opImpl_( gridPart ),
1403 localOp_( gridPart, std::forward< Args >( args )... ),
1404 gridSizeInterior_( 0 ),
1405 communicate_( true )
1406 {
1407 }
1408
1409 void setCommunicate( const bool communicate )
1410 {
1411 communicate_ = communicate;
1412 if( ! communicate_ && Dune::Fem::Parameter::verbose() )
1413 {
1414 std::cout << "GalerkinOperator::setCommunicate: communicate was disabled!" << std::endl;
1415 }
1416 }
1417
1418 void setQuadratureOrders(unsigned int interior, unsigned int surface)
1419 {
1420 size_t size = localOp_.size();
1421 for( size_t i=0; i<size; ++i )
1422 localOp_[ i ].setQuadratureOrders(interior,surface);
1423 }
1424
1425 virtual bool nonlinear() const final override
1426 {
1427 return localOperator().nonlinear();
1428 }
1429
1430 virtual void operator() ( const DomainFunctionType &u, RangeFunctionType &w ) const final override
1431 {
1432 evaluate( u, w );
1433 }
1434
1435 template< class GridFunction >
1436 void operator() ( const GridFunction &u, RangeFunctionType &w ) const
1437 {
1438 evaluate( u, w );
1439 }
1440
1441 const GridPartType &gridPart () const { return op().gridPart(); }
1442
1443 typedef Integrands ModelType;
1444 typedef Integrands DirichletModelType;
1445 ModelType &model() const { return localOperator().model(); }
1446
1447 [[deprecated("Use localOperator instead!")]]
1448 const LocalGalerkinOperatorImplType& impl() const { return localOperator(); }
1449
1451 const LocalGalerkinOperatorImplType& localOperator() const { return *localOp_; }
1452
1453 std::size_t gridSizeInterior () const { return gridSizeInterior_; }
1454
1455 protected:
1457 const GalerkinOperatorImplType& op() const { return *opImpl_; }
1458
1459 template < class GridFunction >
1460 void evaluate( const GridFunction &u, RangeFunctionType &w ) const
1461 {
1462 iterators_.update();
1463 w.clear();
1464
1465 std::shared_mutex mutex;
1466
1467 auto doEval = [this, &u, &w, &mutex] ()
1468 {
1469 // TODO: Move this to be a class variable
1470 std::tuple< const LocalGalerkinOperatorImplType& > integrands( localOperator() );
1471 this->op().evaluate( u, w, this->iterators_, integrands, mutex );
1472 };
1473
1474 try {
1475 // execute in parallel
1476 MPIManager :: run ( doEval );
1477
1478 // update number of interior elements as sum over threads
1479 gridSizeInterior_ = Impl::accumulateGridSize( opImpl_ );
1480 }
1481 catch ( const SingleThreadModeError& e )
1482 {
1483 // reset w from previous entries
1484 w.clear();
1485 // re-run in single thread mode if previous attempt failed
1486 std::tuple< const LocalGalerkinOperatorImplType& > integrands( localOperator() );
1487 op().evaluate( u, w, iterators_, integrands );
1488
1489 // update number of interior elements as sum over threads
1490 gridSizeInterior_ = op().gridSizeInterior();
1491 }
1492
1493 // synchronize result
1494 if( communicate_ )
1495 w.communicate();
1496 }
1497
1498 mutable ThreadIteratorType iterators_;
1499 ThreadSafeValue< GalerkinOperatorImplType > opImpl_;
1500 ThreadSafeValue< LocalGalerkinOperatorImplType > localOp_;
1501
1502 mutable std::size_t gridSizeInterior_;
1503 bool communicate_;
1504 };
1505
1506
1507
1508 // DifferentiableGalerkinOperator
1509 // ------------------------------
1510
1511 template< class Integrands, class JacobianOperator >
1512 class DifferentiableGalerkinOperator
1513 : public GalerkinOperator< Integrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType >,
1514 public DifferentiableOperator< JacobianOperator >
1515 {
1516 typedef GalerkinOperator< Integrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType > BaseType;
1517
1518 typedef typename BaseType :: LocalGalerkinOperatorImplType LocalGalerkinOperatorImplType;
1519 public:
1520 typedef JacobianOperator JacobianOperatorType;
1521
1522 typedef typename BaseType::DomainFunctionType DomainFunctionType;
1523 typedef typename BaseType::RangeFunctionType RangeFunctionType;
1524 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType;
1525 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType;
1526
1527 typedef DiagonalAndNeighborStencil< DomainDiscreteFunctionSpaceType, RangeDiscreteFunctionSpaceType > DiagonalAndNeighborStencilType;
1528 typedef DiagonalStencil< DomainDiscreteFunctionSpaceType, RangeDiscreteFunctionSpaceType > DiagonalStencilType;
1529
1530 typedef typename BaseType::GridPartType GridPartType;
1531
1532 template< class... Args >
1533 explicit DifferentiableGalerkinOperator ( const DomainDiscreteFunctionSpaceType &dSpace,
1534 const RangeDiscreteFunctionSpaceType &rSpace,
1535 Args &&... args )
1536 : BaseType( rSpace.gridPart(), std::forward< Args >( args )... ),
1537 dSpace_(dSpace), rSpace_(rSpace),
1538 domainSpaceSequence_(dSpace.sequence()),
1539 rangeSpaceSequence_(rSpace.sequence()),
1540 stencilDAN_(), stencilD_()
1541 {
1542 if( hasSkeleton() )
1543 stencilDAN_.reset( new DiagonalAndNeighborStencilType( dSpace_, rSpace_ ) );
1544 else
1545 stencilD_.reset( new DiagonalStencilType( dSpace_, rSpace_ ) );
1546 }
1547
1548 virtual void jacobian ( const DomainFunctionType &u, JacobianOperatorType &jOp ) const final override
1549 {
1550 assemble( u, jOp );
1551 }
1552
1553 template< class GridFunction >
1554 void jacobian ( const GridFunction &u, JacobianOperatorType &jOp ) const
1555 {
1556 assemble( u, jOp );
1557 }
1558
1559 const DomainDiscreteFunctionSpaceType& domainSpace() const
1560 {
1561 return dSpace_;
1562 }
1563 const RangeDiscreteFunctionSpaceType& rangeSpace() const
1564 {
1565 return rSpace_;
1566 }
1567
1568 using BaseType::localOperator;
1569 using BaseType::nonlinear;
1570
1571 protected:
1572 using BaseType::op;
1573
1574 bool hasSkeleton() const
1575 {
1576 std::tuple< const LocalGalerkinOperatorImplType& > integrands( localOperator() );
1577 return op().hasSkeleton( integrands );
1578 }
1579
1580 void prepare( JacobianOperatorType& jOp ) const
1581 {
1582 if ( domainSpaceSequence_ != domainSpace().sequence()
1583 || rangeSpaceSequence_ != rangeSpace().sequence() )
1584 {
1585 domainSpaceSequence_ = domainSpace().sequence();
1586 rangeSpaceSequence_ = rangeSpace().sequence();
1587 if( hasSkeleton() )
1588 {
1589 assert( stencilDAN_ );
1590 stencilDAN_->update();
1591 }
1592 else
1593 {
1594 assert( stencilD_ );
1595 stencilD_->update();
1596 }
1597 }
1598 if( hasSkeleton() )
1599 jOp.reserve( *stencilDAN_ );
1600 else
1601 jOp.reserve( *stencilD_ );
1602 // set all entries to zero
1603 jOp.clear();
1604 }
1605
1606 template < class GridFunction >
1607 void assemble( const GridFunction &u, JacobianOperatorType &jOp ) const
1608 {
1609 // reserve memory and clear entries
1610 {
1611 prepare( jOp );
1612 iterators_.update();
1613 }
1614
1615 std::shared_mutex mutex;
1616
1617 auto doAssemble = [this, &u, &jOp, &mutex] ()
1618 {
1619 std::tuple< const LocalGalerkinOperatorImplType& > integrands( localOperator() );
1620 this->op().assemble( u, jOp, this->iterators_, integrands, mutex );
1621 };
1622
1623 try {
1624 // execute in parallel
1625 MPIManager :: run ( doAssemble );
1626
1627 // update number of interior elements as sum over threads
1628 gridSizeInterior_ = Impl::accumulateGridSize( this->opImpl_ );
1629 }
1630 catch ( const SingleThreadModeError& e )
1631 {
1632 // redo assemble since it failed previously
1633 jOp.clear();
1634 std::tuple< const LocalGalerkinOperatorImplType& > integrands( localOperator() );
1635 op().assemble( u, jOp, iterators_, integrands );
1636 // update number of interior elements as sum over threads
1637 gridSizeInterior_ = op().gridSizeInterior();
1638 }
1639
1640 // note: assembly done without local contributions so need
1641 // to call flush assembly
1642 jOp.flushAssembly();
1643 }
1644
1645 using BaseType::iterators_;
1646 using BaseType::gridSizeInterior_;
1647
1648 const DomainDiscreteFunctionSpaceType &dSpace_;
1649 const RangeDiscreteFunctionSpaceType &rSpace_;
1650
1651 mutable int domainSpaceSequence_, rangeSpaceSequence_;
1652
1653 mutable std::unique_ptr< DiagonalAndNeighborStencilType > stencilDAN_;
1654 mutable std::unique_ptr< DiagonalStencilType > stencilD_;
1655 };
1656
1657
1658
1659 // AutomaticDifferenceGalerkinOperator
1660 // -----------------------------------
1661
1662 template< class Integrands, class DomainFunction, class RangeFunction >
1663 class AutomaticDifferenceGalerkinOperator
1664 : public GalerkinOperator< Integrands, DomainFunction, RangeFunction >,
1665 public AutomaticDifferenceOperator< DomainFunction, RangeFunction >
1666 {
1667 typedef GalerkinOperator< Integrands, DomainFunction, RangeFunction > BaseType;
1668 typedef AutomaticDifferenceOperator< DomainFunction, RangeFunction > AutomaticDifferenceOperatorType;
1669
1670 public:
1671 typedef typename BaseType::GridPartType GridPartType;
1672
1673 template< class... Args >
1674 explicit AutomaticDifferenceGalerkinOperator ( const GridPartType &gridPart, Args &&... args )
1675 : BaseType( gridPart, std::forward< Args >( args )... ), AutomaticDifferenceOperatorType()
1676 {}
1677 };
1678
1679
1680
1681 // ModelDifferentiableGalerkinOperator
1682 // -----------------------------------
1683
1684 template < class LinearOperator, class ModelIntegrands >
1685 struct ModelDifferentiableGalerkinOperator
1686 : public DifferentiableGalerkinOperator< ModelIntegrands, LinearOperator >
1687 {
1688 typedef DifferentiableGalerkinOperator< ModelIntegrands, LinearOperator > BaseType;
1689
1690 typedef typename ModelIntegrands::ModelType ModelType;
1691
1692 typedef typename LinearOperator::DomainFunctionType RangeFunctionType;
1693 typedef typename LinearOperator::RangeSpaceType DiscreteFunctionSpaceType;
1694
1695 ModelDifferentiableGalerkinOperator ( ModelType &model, const DiscreteFunctionSpaceType &dfSpace )
1696 : BaseType( dfSpace.gridPart(), model )
1697 {}
1698
1699 template< class GridFunction >
1700 void apply ( const GridFunction &u, RangeFunctionType &w ) const
1701 {
1702 (*this)( u, w );
1703 }
1704
1705 template< class GridFunction >
1706 void apply ( const GridFunction &u, LinearOperator &jOp ) const
1707 {
1708 (*this).jacobian( u, jOp );
1709 }
1710 };
1711
1712 namespace Impl
1713 {
1714
1715 // GalerkinSchemeImpl
1716 // ------------------
1717 template< class Integrands, class LinearOperator, bool addDirichletBC,
1718 template <class,class> class DifferentiableGalerkinOperatorImpl >
1719 struct GalerkinSchemeTraits
1720 {
1721 template <class O, bool addDBC>
1722 struct DirichletBlockSelector { using type = void; };
1723 template <class O>
1724 struct DirichletBlockSelector<O,true> { using type = typename O::DirichletBlockVector; };
1725
1726 using DifferentiableOperatorType = std::conditional_t< addDirichletBC,
1727 DirichletWrapperOperator< DifferentiableGalerkinOperatorImpl< Integrands, LinearOperator >>,
1728 DifferentiableGalerkinOperatorImpl< Integrands, LinearOperator > >;
1729 using DirichletBlockVector = typename DirichletBlockSelector<
1730 DirichletWrapperOperator<
1731 DifferentiableGalerkinOperatorImpl< Integrands, LinearOperator >>,
1732 addDirichletBC>::type;
1733
1734 typedef DifferentiableOperatorType type;
1735 };
1736
1737 template< class Integrands, class LinearOperator, class LinearInverseOperator, bool addDirichletBC,
1738 template <class,class> class DifferentiableGalerkinOperatorImpl = DifferentiableGalerkinOperator >
1739 struct GalerkinSchemeImpl : public FemScheme< typename
1740 GalerkinSchemeTraits< Integrands, LinearOperator,
1741 addDirichletBC, DifferentiableGalerkinOperatorImpl>::type, // Operator
1742 LinearInverseOperator > // LinearInverseOperator
1743 {
1744 typedef FemScheme< typename GalerkinSchemeTraits< Integrands, LinearOperator,
1745 addDirichletBC, DifferentiableGalerkinOperatorImpl>::type, // Operator
1746 LinearInverseOperator > // LinearInverseOperator
1747 BaseType;
1748
1749 typedef typename BaseType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
1750
1751 GalerkinSchemeImpl ( const DiscreteFunctionSpaceType &dfSpace,
1752 const Integrands &integrands,
1753 const ParameterReader& parameter = Parameter::container() )
1754 : BaseType(dfSpace,
1755 parameter,
1756 std::move(integrands))
1757 {}
1758 };
1759
1760 } // end namespace Impl
1761
1762 // GalerkinScheme
1763 // --------------
1764
1765 template< class Integrands, class LinearOperator, class InverseOperator, bool addDirichletBC >
1766 using GalerkinScheme = Impl::GalerkinSchemeImpl< Integrands, LinearOperator, InverseOperator, addDirichletBC,
1767 DifferentiableGalerkinOperator >;
1768
1769 } // namespace Fem
1770
1771} // namespace Dune
1772
1773#endif // #ifndef DUNE_FEM_SCHEMES_GALERKIN_HH
Dune::Fem::LocalFunction::RangeType
FunctionSpaceType::RangeType RangeType
type of range vectors, i.e., type of function values
Definition: localfunction.hh:110
Dune::Fem::LocalFunction::HessianRangeType
FunctionSpaceType::HessianRangeType HessianRangeType
type of the Hessian
Definition: localfunction.hh:114
Dune::Fem::LocalFunction::JacobianRangeType
FunctionSpaceType::JacobianRangeType JacobianRangeType
type of the Jacobian, i.e., type of evaluated Jacobian matrix
Definition: localfunction.hh:112
Dune::Fem::Parameter::verbose
static bool verbose()
obtain the cached value for fem.verbose with default verbosity level 2
Definition: parameter.hh:466
ElementQuadrature
quadrature on the codim-0 reference element
Quadrature
actual interface class for quadratures
Dune::InteriorEntity
@ InteriorEntity
all interior entities
Definition: gridenums.hh:31
Dune::Hybrid::forEach
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:261
Dune::range
static constexpr IntegralRange< std::decay_t< T > > range(T &&from, U &&to) noexcept
free standing function for setting up a range based for loop over an integer range for (auto i: range...
Definition: rangeutilities.hh:288
Dune::Partitions::interior
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:271
Dune
Dune namespace
Definition: alignedallocator.hh:13
Dune::size
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Dune::empty
constexpr std::bool_constant<(sizeof...(II)==0)> empty(std::integer_sequence< T, II... >)
Checks whether the sequence is empty.
Definition: integersequence.hh:80
std
STL namespace.
Dune::Fem::Operator< DomainFunction, DomainFunction >::DomainFunctionType
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:36
timer.hh
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Mar 8, 23:37, 2026)