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