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