dune-fem 2.8-git
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 <vector>
10
11#include <dune/common/hybridutilities.hh>
12
13#include <dune/grid/common/rangegenerators.hh>
14
26
32
34
35// fempy includes
36#include <dune/fempy/quadrature/fempyquadratures.hh>
37
38namespace Dune
39{
40
41 namespace Fem
42 {
43
44 namespace Impl
45 {
46 template <class M>
47 class CallOrder
48 {
49 // check for 'int order () const' or 'int order()'
50 // and variants returning unsigned int or size_t
51 template <class T> static std::true_type testSignature(int (T::*)() const);
52 template <class T> static std::true_type testSignature(int (T::*)());
53 template <class T> static std::true_type testSignature(unsigned int (T::*)() const);
54 template <class T> static std::true_type testSignature(unsigned int (T::*)());
55 template <class T> static std::true_type testSignature(std::size_t (T::*)() const);
56 template <class T> static std::true_type testSignature(std::size_t (T::*)());
57
58 template <class T>
59 static decltype(testSignature(&T::order)) test(std::nullptr_t);
60
61 template <class T>
62 static std::false_type test(...);
63
64 using type = decltype(test<M>(nullptr));
65
66 template <class F>
67 static int callOrder(const F& f, std::false_type)
68 {
69#ifndef NDEBUG
70 std::cerr << "WARNING: not order method available on " << typeid(F).name() << ", defaulting to 1!" << std::endl;
71#endif
72 return 1;
73 }
74
75 template <class F>
76 static int callOrder(const F& f, std::true_type)
77 {
78 return f.order();
79 }
80
81 public:
82 template <class F>
83 static int order (const F& f ) { return callOrder(f, type() ); }
84 };
85
86 // GalerkinOperator
87 // ----------------
88
89 template< class Integrands >
90 struct GalerkinOperator
91 {
92 typedef GalerkinOperator<Integrands> ThisType;
93 typedef std::conditional_t< Fem::IntegrandsTraits< Integrands >::isFull, Integrands, FullIntegrands< Integrands > > IntegrandsType;
94
95 typedef typename IntegrandsType::GridPartType GridPartType;
96
97 typedef typename GridPartType::ctype ctype;
98 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
99
100 template <class Space>
101 struct QuadratureSelector
102 {
103 typedef CachingQuadrature< GridPartType, 0, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > InteriorQuadratureType;
104 typedef CachingQuadrature< GridPartType, 1, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > SurfaceQuadratureType;
105 // typedef CachingQuadrature< GridPartType, 0, Dune::FemPy::FempyQuadratureTraits > InteriorQuadratureType;
106 // typedef CachingQuadrature< GridPartType, 1, Dune::FemPy::FempyQuadratureTraits > SurfaceQuadratureType;
107 };
108
109 private:
110 typedef typename IntegrandsType::DomainValueType DomainValueType;
111 typedef typename IntegrandsType::RangeValueType RangeValueType;
112 typedef std::make_index_sequence< std::tuple_size< DomainValueType >::value > DomainValueIndices;
113 typedef std::make_index_sequence< std::tuple_size< RangeValueType >::value > RangeValueIndices;
114
115
116 template< std::size_t... i >
117 static auto makeDomainValueVector ( std::size_t maxNumLocalDofs, std::index_sequence< i... > )
118 {
119 return std::make_tuple( std::vector< std::tuple_element_t< i, DomainValueType > >( maxNumLocalDofs )... );
120 }
121
122 static auto makeDomainValueVector ( std::size_t maxNumLocalDofs )
123 {
124 return makeDomainValueVector( maxNumLocalDofs, DomainValueIndices() );
125 }
126
127 template< std::size_t... i >
128 static auto makeRangeValueVector ( std::size_t maxNumLocalDofs, std::index_sequence< i... > )
129 {
130 return std::make_tuple( std::vector< std::tuple_element_t< i, RangeValueType > >( maxNumLocalDofs )... );
131 }
132
133 static auto makeRangeValueVector ( std::size_t maxNumLocalDofs )
134 {
135 return makeRangeValueVector( maxNumLocalDofs, RangeValueIndices() );
136 }
137
138 typedef decltype( makeDomainValueVector( 0u ) ) DomainValueVectorType;
139 typedef decltype( makeRangeValueVector( 0u ) ) RangeValueVectorType;
140
141 static void resizeDomainValueVector ( DomainValueVectorType& vec, const std::size_t size )
142 {
143 Hybrid::forEach( DomainValueIndices(), [ &vec, &size ] ( auto i ) {
144 std::get< i >( vec ).resize( size );
145 } );
146 }
147
148 static void resizeRangeValueVector ( RangeValueVectorType& vec, const std::size_t size )
149 {
150 Hybrid::forEach( RangeValueIndices(), [ &vec, &size ] ( auto i ) {
151 std::get< i >( vec ).resize( size );
152 } );
153 }
154
155 template< class LocalFunction, class Quadrature >
156 static void evaluateQuadrature ( const LocalFunction &u, const Quadrature &quad, std::vector< typename LocalFunction::RangeType > &phi )
157 {
158 u.evaluateQuadrature( quad, phi );
159 }
160
161 template< class LocalFunction, class Quadrature>
162 static void evaluateQuadrature ( const LocalFunction &u, const Quadrature &quad, std::vector< typename LocalFunction::JacobianRangeType > &phi )
163 {
164 u.jacobianQuadrature( quad, phi );
165 }
166
167 template< class LocalFunction, class Quadrature >
168 static void evaluateQuadrature ( const LocalFunction &u, const Quadrature &quad, std::vector< typename LocalFunction::HessianRangeType > &phi )
169 {
170 u.hessianQuadrature( quad, phi );
171 }
172
173 template< class LocalFunction, class Point >
174 static void value ( const LocalFunction &u, const Point &x, typename LocalFunction::RangeType &phi )
175 {
176 u.evaluate( x, phi );
177 }
178
179 template< class LocalFunction, class Point >
180 static void value ( const LocalFunction &u, const Point &x, typename LocalFunction::JacobianRangeType &phi )
181 {
182 u.jacobian( x, phi );
183 }
184
185 template< class LocalFunction, class Point >
186 static void value ( const LocalFunction &u, const Point &x, typename LocalFunction::HessianRangeType &phi )
187 {
188 u.hessian( x, phi );
189 }
190
191 template< class LocalFunction, class Point, class... T >
192 static void value ( const LocalFunction &u, const Point &x, std::tuple< T... > &phi )
193 {
194 Hybrid::forEach( std::index_sequence_for< T... >(), [ &u, &x, &phi ] ( auto i ) { GalerkinOperator::value( u, x, std::get< i >( phi ) ); } );
195 }
196
197 template< class Basis, class Point >
198 static void values ( const Basis &basis, const Point &x, std::vector< typename Basis::RangeType > &phi )
199 {
200 basis.evaluateAll( x, phi );
201 }
202
203 template< class Basis, class Point >
204 static void values ( const Basis &basis, const Point &x, std::vector< typename Basis::JacobianRangeType > &phi )
205 {
206 basis.jacobianAll( x, phi );
207 }
208
209 template< class Basis, class Point >
210 static void values ( const Basis &basis, const Point &x, std::vector< typename Basis::HessianRangeType > &phi )
211 {
212 basis.hessianAll( x, phi );
213 }
214
215 template< class Basis, class Point, class... T >
216 static void values ( const Basis &basis, const Point &x, std::tuple< std::vector< T >... > &phi )
217 {
218 Hybrid::forEach( std::index_sequence_for< T... >(), [ &basis, &x, &phi ] ( auto i ) { GalerkinOperator::values( basis, x, std::get< i >( phi ) ); } );
219 }
220
221 template< class LocalFunction, class Point >
222 static DomainValueType domainValue ( const LocalFunction &u, const Point &x )
223 {
224 DomainValueType phi;
225 value( u, x, phi );
226 return phi;
227 }
228
229 static DomainValueType domainValue ( const unsigned int qpIdx, DomainValueVectorType& vec)
230 {
231 DomainValueType phi;
232 Hybrid::forEach( DomainValueIndices(), [ &qpIdx, &vec, &phi ] ( auto i ) {
233 std::get< i > ( phi ) = std::get< i >( vec )[ qpIdx ];
234 } );
235 return phi;
236 }
237
238 template< class LocalFunction, class Quadrature >
239 static void domainValue ( const LocalFunction &u, const Quadrature& quadrature, DomainValueVectorType &result )
240 {
241 Hybrid::forEach( DomainValueIndices(), [ &u, &quadrature, &result ] ( auto i ) {
242 auto& vec = std::get< i >( result );
243 vec.resize( quadrature.nop() );
244 ThisType::evaluateQuadrature( u, quadrature, vec );
245 } );
246 }
247
248 template< class Phi, std::size_t... i >
249 static auto value ( const Phi &phi, std::size_t col, std::index_sequence< i... > )
250 {
251 return std::make_tuple( std::get< i >( phi )[ col ]... );
252 }
253
254 template< class... T >
255 static auto value ( const std::tuple< std::vector< T >... > &phi, std::size_t col )
256 {
257 return value( phi, col, std::index_sequence_for< T... >() );
258 }
259
260 static void assignRange( RangeValueVectorType& ranges, const std::size_t idx, const RangeValueType& range )
261 {
262 Hybrid::forEach( RangeValueIndices(), [ &ranges, &idx, &range ] ( auto i ) {
263 std::get< i >( ranges )[ idx ] = std::get< i >( range );
264 });
265 }
266 template <class W>
267 static void assignRange( RangeValueVectorType& ranges, const std::size_t idx, const RangeValueType& range, const W &weight )
268 {
269 Hybrid::forEach( RangeValueIndices(), [ &ranges, &idx, &range, &weight ] ( auto i ) {
270 std::get< i >( ranges )[ idx ] = std::get< i >( range );
271 std::get< i >( ranges )[ idx ] *= weight;
272 });
273 }
274
275 static void assignDomain( DomainValueVectorType& domains, const std::size_t idx, const DomainValueType& domain )
276 {
277 Hybrid::forEach( DomainValueIndices(), [ &domains, &idx, &domain ] ( auto i ) {
278 std::get< i >( domains )[ idx ] = std::get< i >( domain );
279 });
280 }
281
282 template <class W, class Quadrature>
283 static void axpyQuadrature( W& w, const Quadrature& quadrature, RangeValueVectorType& ranges )
284 {
285 Hybrid::forEach( RangeValueIndices(), [ &w, &quadrature, &ranges ] ( auto i ) {
286 w.axpyQuadrature( quadrature, std::get< i >( ranges ) );
287 } );
288 }
289
290 public:
291 // interior integral
292
293 template< class U, class W >
294 void addInteriorIntegral ( const U &u, W &w ) const
295 {
296 if( !integrands_.init( u.entity() ) )
297 return;
298
299 const auto geometry = u.entity().geometry();
300
301 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: InteriorQuadratureType InteriorQuadratureType;
302 InteriorQuadratureType quadrature( u.entity(), interiorQuadratureOrder(maxOrder(u, w)) );
303
304 // evaluate u for all quadrature points
305 DomainValueVectorType& domains = domainValues_;
306 domainValue( u, quadrature, domains );
307
308 auto& ranges = values_;
309 resizeRangeValueVector( ranges, quadrature.nop() );
310
311 // evaluate integrands for all quadrature points
312 for( const auto qp : quadrature )
313 {
314 const ctype weight = qp.weight() * geometry.integrationElement( qp.position() );
315 assignRange( ranges, qp.index(), integrands_.interior( qp, domainValue( qp.index(), domains ) ), weight );
316 }
317
318 // add to w for all quadrature points
319 axpyQuadrature( w, quadrature, ranges );
320 }
321
322 template< class U, class J >
323 void addLinearizedInteriorIntegral ( const U &u, DomainValueVectorType &phi, J &j ) const
324 {
325 if( !integrands_.init( u.entity() ) )
326 return;
327
328 const auto geometry = u.entity().geometry();
329 const auto &domainBasis = j.domainBasisFunctionSet();
330 const auto &rangeBasis = j.rangeBasisFunctionSet();
331
332 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: InteriorQuadratureType InteriorQuadratureType;
333 InteriorQuadratureType quadrature( u.entity(), interiorQuadratureOrder( maxOrder( u, domainBasis, rangeBasis )) );
334 const size_t domainSize = domainBasis.size();
335 const size_t quadNop = quadrature.nop();
336
337 auto& basisValues = basisValues_;
338 resizeDomainValueVector( basisValues_, domainSize );
339
340 // evaluate u for all quadrature points
341 DomainValueVectorType& domains = domainValues_;
342 domainValue( u, quadrature, domains );
343
344 rangeValues_.resize( domainSize );
345 for( std::size_t col = 0; col < domainSize; ++col )
346 {
347 resizeRangeValueVector( rangeValues_[ col ], quadNop );
348 }
349
350 // evaluate all basis functions and integrands
351 for( const auto qp : quadrature )
352 {
353 values( domainBasis, qp, basisValues );
354 const auto weight = qp.weight() * geometry.integrationElement( qp.position() );
355 auto integrand = integrands_.linearizedInterior( qp, domainValue( qp.index(), domains ) );
356 for( std::size_t col = 0; col < domainSize; ++col )
357 {
358 assignRange( rangeValues_[ col ], qp.index(), integrand( value( basisValues, col ) ), weight );
359 }
360 }
361
362 // add to local matrix for all quadrature points and basis functions
363 for( std::size_t col = 0; col < domainSize; ++col )
364 {
365 LocalMatrixColumn< J > jCol( j, col );
366 axpyQuadrature( jCol, quadrature, rangeValues_[ col ] );
367 }
368 }
369
370 // boundary integral
371
372 template< class Intersection, class U, class W >
373 void addBoundaryIntegral ( const Intersection &intersection, const U &u, W &w ) const
374 {
375 if( !integrands_.init( intersection ) )
376 return;
377
378 const auto geometry = intersection.geometry();
379 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
380 SurfaceQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( u, w )), SurfaceQuadratureType::INSIDE );
381 for( const auto qp : quadrature )
382 {
383 const ctype weight = qp.weight() * geometry.integrationElement( qp.localPosition() );
384
385 RangeValueType integrand = integrands_.boundary( qp, domainValue( u, qp ) );
386
387 Hybrid::forEach( RangeValueIndices(), [ &qp, &w, &integrand, weight ] ( auto i ) {
388 std::get< i >( integrand ) *= weight;
389 w.axpy( qp, std::get< i >( integrand ) );
390 } );
391 }
392 }
393
394 template< class Intersection, class U, class J >
395 void addLinearizedBoundaryIntegral ( const Intersection &intersection, const U &u, DomainValueVectorType &phi, J &j ) const
396 {
397 if( !integrands_.init( intersection ) )
398 return;
399
400 const auto geometry = intersection.geometry();
401 const auto &domainBasis = j.domainBasisFunctionSet();
402 const auto &rangeBasis = j.rangeBasisFunctionSet();
403
404 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
405 SurfaceQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder(u, domainBasis, rangeBasis )), SurfaceQuadratureType::INSIDE );
406 for( const auto qp : quadrature )
407 {
408 const ctype weight = qp.weight() * geometry.integrationElement( qp.localPosition() );
409
410 values( domainBasis, qp, phi );
411 auto integrand = integrands_.linearizedBoundary( qp, domainValue( u, qp ) );
412
413 for( std::size_t col = 0, cols = domainBasis.size(); col < cols; ++col )
414 {
415 LocalMatrixColumn< J > jCol( j, col );
416 RangeValueType intPhi = integrand( value( phi, col ) );
417
418 Hybrid::forEach( RangeValueIndices(), [ &qp, &jCol, &intPhi, weight ] ( auto i ) {
419 std::get< i >( intPhi ) *= weight;
420 jCol.axpy( qp, std::get< i >( intPhi ) );
421 } );
422 }
423 }
424 }
425
426 // addSkeletonIntegral
427
428 private:
429 template< bool conforming, class Intersection, class U, class W >
430 void addSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut, W &wIn ) const
431 {
432 const auto geometry = intersection.geometry();
433
434 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
435 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
436 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( uIn, uOut, wIn)), false );
437 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
438 {
439 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
440
441 const auto qpIn = quadrature.inside()[ qp ];
442 const auto qpOut = quadrature.outside()[ qp ];
443 std::pair< RangeValueType, RangeValueType > integrand = integrands_.skeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
444
445 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &wIn, &integrand, weight ] ( auto i ) {
446 std::get< i >( integrand.first ) *= weight;
447 wIn.axpy( qpIn, std::get< i >( integrand.first ) );
448 } );
449 }
450 }
451
452 template< bool conforming, class Intersection, class U, class W >
453 void addSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut, W &wIn, W &wOut ) const
454 {
455 const auto geometry = intersection.geometry();
456 typedef typename QuadratureSelector< typename W::DiscreteFunctionSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
457 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
458 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(maxOrder( uIn, uOut, wIn, wOut)), false );
459 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
460 {
461 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
462
463 const auto qpIn = quadrature.inside()[ qp ];
464 const auto qpOut = quadrature.outside()[ qp ];
465 std::pair< RangeValueType, RangeValueType > integrand = integrands_.skeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
466
467 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &wIn, &qpOut, &wOut, &integrand, weight ] ( auto i ) {
468 std::get< i >( integrand.first ) *= weight;
469 wIn.axpy( qpIn, std::get< i >( integrand.first ) );
470
471 std::get< i >( integrand.second ) *= weight;
472 wOut.axpy( qpOut, std::get< i >( integrand.second ) );
473 } );
474 }
475 }
476
477 template< bool conforming, class Intersection, class U, class J >
478 void addLinearizedSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut,
479 DomainValueVectorType &phiIn, DomainValueVectorType &phiOut, J &jInIn, J &jOutIn ) const
480 {
481 const auto &domainBasisIn = jInIn.domainBasisFunctionSet();
482 const auto &domainBasisOut = jOutIn.domainBasisFunctionSet();
483
484 const auto &rangeBasisIn = jInIn.rangeBasisFunctionSet();
485
486 const int order = std::max( maxOrder(uIn, uOut), maxOrder( domainBasisIn, domainBasisOut, rangeBasisIn ));
487
488 const auto geometry = intersection.geometry();
489 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
490 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
491 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(order), false );
492 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
493 {
494 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
495
496 const auto qpIn = quadrature.inside()[ qp ];
497 const auto qpOut = quadrature.outside()[ qp ];
498
499 values( domainBasisIn, qpIn, phiIn );
500 values( domainBasisOut, qpOut, phiOut );
501
502 auto integrand = integrands_.linearizedSkeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
503 for( std::size_t col = 0, cols = domainBasisIn.size(); col < cols; ++col )
504 {
505 LocalMatrixColumn< J > jInInCol( jInIn, col );
506 std::pair< RangeValueType, RangeValueType > intPhi = integrand.first( value( phiIn, col ) );
507
508 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jInInCol, &intPhi, weight ] ( auto i ) {
509 std::get< i >( intPhi.first ) *= weight;
510 jInInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
511 } );
512 }
513 for( std::size_t col = 0, cols = domainBasisOut.size(); col < cols; ++col )
514 {
515 LocalMatrixColumn< J > jOutInCol( jOutIn, col );
516 std::pair< RangeValueType, RangeValueType > intPhi = integrand.second( value( phiOut, col ) );
517
518 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jOutInCol, &intPhi, weight ] ( auto i ) {
519 std::get< i >( intPhi.first ) *= weight;
520 jOutInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
521 } );
522 }
523 }
524 }
525
526 template< bool conforming, class Intersection, class U, class J >
527 void addLinearizedSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut,
528 DomainValueVectorType &phiIn, DomainValueVectorType &phiOut, J &jInIn, J &jOutIn, J &jInOut, J &jOutOut ) const
529 {
530 const auto &domainBasisIn = jInIn.domainBasisFunctionSet();
531 const auto &domainBasisOut = jOutIn.domainBasisFunctionSet();
532
533 const auto &rangeBasisIn = jInIn.rangeBasisFunctionSet();
534 const auto &rangeBasisOut = jInOut.rangeBasisFunctionSet();
535
536 const int order = std::max( maxOrder(uIn, uOut), maxOrder( domainBasisIn, domainBasisOut, rangeBasisIn, rangeBasisOut ));
537
538 const auto geometry = intersection.geometry();
539 typedef typename QuadratureSelector< typename J::RangeSpaceType > :: SurfaceQuadratureType SurfaceQuadratureType;
540 typedef IntersectionQuadrature< SurfaceQuadratureType, conforming > IntersectionQuadratureType;
541 const IntersectionQuadratureType quadrature( gridPart(), intersection, surfaceQuadratureOrder(order), false );
542 for( std::size_t qp = 0, nop = quadrature.nop(); qp != nop; ++qp )
543 {
544 const ctype weight = quadrature.weight( qp ) * geometry.integrationElement( quadrature.localPoint( qp ) );
545
546 const auto qpIn = quadrature.inside()[ qp ];
547 const auto qpOut = quadrature.outside()[ qp ];
548
549 values( domainBasisIn, qpIn, phiIn );
550 values( domainBasisOut, qpOut, phiOut );
551
552 auto integrand = integrands_.linearizedSkeleton( qpIn, domainValue( uIn, qpIn ), qpOut, domainValue( uOut, qpOut ) );
553 for( std::size_t col = 0, cols = domainBasisIn.size(); col < cols; ++col )
554 {
555 LocalMatrixColumn< J > jInInCol( jInIn, col );
556 LocalMatrixColumn< J > jInOutCol( jInOut, col );
557 std::pair< RangeValueType, RangeValueType > intPhi = integrand.first( value( phiIn, col ) );
558
559 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jInInCol, &qpOut, &jInOutCol, &intPhi, weight ] ( auto i ) {
560 std::get< i >( intPhi.first ) *= weight;
561 jInInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
562
563 std::get< i >( intPhi.second ) *= weight;
564 jInOutCol.axpy( qpOut, std::get< i >( intPhi.second ) );
565 } );
566 }
567 for( std::size_t col = 0, cols = domainBasisOut.size(); col < cols; ++col )
568 {
569 LocalMatrixColumn< J > jOutInCol( jOutIn, col );
570 LocalMatrixColumn< J > jOutOutCol( jOutOut, col );
571 std::pair< RangeValueType, RangeValueType > intPhi = integrand.second( value( phiOut, col ) );
572
573 Hybrid::forEach( RangeValueIndices(), [ &qpIn, &jOutInCol, &qpOut, &jOutOutCol, &intPhi, weight ] ( auto i ) {
574 std::get< i >( intPhi.first ) *= weight;
575 jOutInCol.axpy( qpIn, std::get< i >( intPhi.first ) );
576
577 std::get< i >( intPhi.second ) *= weight;
578 jOutOutCol.axpy( qpOut, std::get< i >( intPhi.second ) );
579 } );
580 }
581 }
582 }
583
584 public:
585 template< class Intersection, class U, class... W >
586 void addSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut, W &... w ) const
587 {
588 if( !integrands_.init( intersection ) )
589 return;
590
591 if( intersection.conforming() )
592 addSkeletonIntegral< true >( intersection, uIn, uOut, w... );
593 else
594 addSkeletonIntegral< false >( intersection, uIn, uOut, w... );
595 }
596
597 template< class Intersection, class U, class... J >
598 void addLinearizedSkeletonIntegral ( const Intersection &intersection, const U &uIn, const U &uOut, DomainValueVectorType &phiIn, DomainValueVectorType &phiOut, J &... j ) const
599 {
600 if( !integrands_.init( intersection ) )
601 return;
602
603 if( intersection.conforming() )
604 addLinearizedSkeletonIntegral< true >( intersection, uIn, uOut, phiIn, phiOut, j... );
605 else
606 addLinearizedSkeletonIntegral< false >( intersection, uIn, uOut, phiIn, phiOut, j... );
607 }
608
609 // constructor
610
611 template< class... Args >
612 explicit GalerkinOperator ( const GridPartType &gridPart, Args &&... args )
613 : gridPart_( gridPart ),
614 integrands_( std::forward< Args >( args )... ),
615 defaultInteriorOrder_( [] (const int order) { return 2 * order; } ),
616 defaultSurfaceOrder_ ( [] (const int order) { return 2 * order + 1; } ),
617 interiorQuadOrder_(0), surfaceQuadOrder_(0),
618 communicate_( true )
619 {
620 }
621
622 void setCommunicate( const bool communicate )
623 {
624 communicate_ = communicate;
625 if( ! communicate_ && Dune::Fem::Parameter::verbose() )
626 {
627 std::cout << "Impl::GalerkinOperator::setCommunicate: communicate was disabled!" << std::endl;
628 }
629 }
630
631 void setQuadratureOrders(unsigned int interior, unsigned int surface)
632 {
633 interiorQuadOrder_ = interior;
634 surfaceQuadOrder_ = surface;
635 }
636
637 IntegrandsType &model() const
638 {
639 return integrands_;
640 }
641
642 private:
643 template< class GridFunction, class DiscreteFunction >
644 void evaluate ( const GridFunction &u, DiscreteFunction &w, std::false_type ) const
645 {
646 w.clear();
647
648 typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
649
651
652 defaultInteriorOrder_ = [] (const int order) { return Capabilities::DefaultQuadrature< DiscreteFunctionSpaceType >::volumeOrder(order); };
653 defaultSurfaceOrder_ = [] (const int order) { return Capabilities::DefaultQuadrature< DiscreteFunctionSpaceType >::surfaceOrder(order); };
654
656 for( const EntityType &entity : elements( gridPart(), Partitions::interiorBorder ) )
657 {
658 uLocal.bind( entity );
659 wLocal.bind( entity );
660 wLocal.clear();
661
662 if( integrands_.hasInterior() )
663 addInteriorIntegral( uLocal, wLocal );
664
665 if( integrands_.hasBoundary() && entity.hasBoundaryIntersections() )
666 {
667 for( const auto &intersection : intersections( gridPart(), entity ) )
668 {
669 if( intersection.boundary() )
670 addBoundaryIntegral( intersection, uLocal, wLocal );
671 }
672 }
673
674 w.addLocalDofs( entity, wLocal.localDofVector() );
675 }
676 }
677
678 template< class GridFunction, class DiscreteFunction >
679 void evaluate ( const GridFunction &u, DiscreteFunction &w, std::true_type ) const
680 {
681 w.clear();
682
683 typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
684
685 TemporaryLocalFunction< DiscreteFunctionSpaceType > wInside( w.space() ), wOutside( w.space() );
686
687 defaultInteriorOrder_ = [] (const int order) { return Capabilities::DefaultQuadrature< DiscreteFunctionSpaceType >::volumeOrder(order); };
688 defaultSurfaceOrder_ = [] (const int order) { return Capabilities::DefaultQuadrature< DiscreteFunctionSpaceType >::surfaceOrder(order); };
689
691
692 const auto &indexSet = gridPart().indexSet();
693 for( const EntityType &inside : elements( gridPart(), Partitions::interiorBorder ) )
694 {
695 uInside.bind( inside );
696 wInside.bind( inside );
697 wInside.clear();
698
699 if( integrands_.hasInterior() )
700 addInteriorIntegral( uInside, wInside );
701
702 for( const auto &intersection : intersections( gridPart(), inside ) )
703 {
704 // check neighbor first since on periodic boundaries both,
705 // neighbor and boundary are true, so we treat neighbor first
706 if( intersection.neighbor() )
707 {
708 const EntityType &outside = intersection.outside();
709
711 uOutside.bind( outside );
712 if( outside.partitionType() != InteriorEntity )
713 addSkeletonIntegral( intersection, uInside, uOutside, wInside );
714 else if( indexSet.index( inside ) < indexSet.index( outside ) )
715 {
716 wOutside.bind( outside );
717 wOutside.clear();
718 addSkeletonIntegral( intersection, uInside, uOutside, wInside, wOutside );
719 w.addLocalDofs( outside, wOutside.localDofVector() );
720 }
721 }
722 else if( intersection.boundary() )
723 {
724 if( integrands_.hasBoundary() )
725 addBoundaryIntegral( intersection, uInside, wInside );
726 }
727 }
728
729 w.addLocalDofs( inside, wInside.localDofVector() );
730 }
731 }
732
733 public:
734 template< class GridFunction, class DiscreteFunction >
735 void evaluate ( const GridFunction &u, DiscreteFunction &w ) const
736 {
737 static_assert( std::is_same< typename GridFunction::GridPartType, GridPartType >::value, "Argument 'u' and Integrands must be defined on the same grid part." );
738 static_assert( std::is_same< typename DiscreteFunction::GridPartType, GridPartType >::value, "Argument 'w' and Integrands must be defined on the same grid part." );
739
740 if( integrands_.hasSkeleton() )
741 evaluate( u, w, std::true_type() );
742 else
743 evaluate( u, w, std::false_type() );
744
745 // synchronize result
746 if( communicate_ )
747 w.communicate();
748 }
749
750 // assemble
751
752 private:
753 template< class GridFunction, class JacobianOperator >
754 void assemble ( const GridFunction &u, JacobianOperator &jOp, std::false_type ) const
755 {
756 typedef typename JacobianOperator::DomainSpaceType DomainSpaceType;
757 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
758
759 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
760
761 // select correct default quadrature orders
762 defaultInteriorOrder_ = [] (const int order) { return Capabilities::DefaultQuadrature< RangeSpaceType >::volumeOrder(order); };
763 defaultSurfaceOrder_ = [] (const int order) { return Capabilities::DefaultQuadrature< RangeSpaceType >::surfaceOrder(order); };
764
765 DiagonalStencil< DomainSpaceType, RangeSpaceType > stencil( jOp.domainSpace(), jOp.rangeSpace() );
766 jOp.reserve( stencil );
767 jOp.clear();
768
769 const std::size_t maxNumLocalDofs = jOp.domainSpace().blockMapper().maxNumDofs() * jOp.domainSpace().localBlockSize;
770 DomainValueVectorType phi = makeDomainValueVector( maxNumLocalDofs );
771
772 TemporaryLocalMatrixType jOpLocal( jOp.domainSpace(), jOp.rangeSpace() );
774
775 for( const EntityType &entity : elements( gridPart(), Partitions::interiorBorder ) )
776 {
777 uLocal.bind( entity );
778
779 jOpLocal.init( entity, entity );
780 jOpLocal.clear();
781
782 if( integrands_.hasInterior() )
783 addLinearizedInteriorIntegral( uLocal, phi, jOpLocal );
784
785 if( integrands_.hasBoundary() && entity.hasBoundaryIntersections() )
786 {
787 for( const auto &intersection : intersections( gridPart(), entity ) )
788 {
789 if( intersection.boundary() )
790 addLinearizedBoundaryIntegral( intersection, uLocal, phi, jOpLocal );
791 }
792 }
793
794 jOp.addLocalMatrix( entity, entity, jOpLocal );
795 }
796 }
797
798 template< class GridFunction, class JacobianOperator >
799 void assemble ( const GridFunction &u, JacobianOperator &jOp, std::true_type ) const
800 {
801 typedef typename JacobianOperator::DomainSpaceType DomainSpaceType;
802 typedef typename JacobianOperator::RangeSpaceType RangeSpaceType;
803
804 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
805
806 // select correct default quadrature orders
807 defaultInteriorOrder_ = [] (const int order) { return Capabilities::DefaultQuadrature< RangeSpaceType >::volumeOrder(order); };
808 defaultSurfaceOrder_ = [] (const int order) { return Capabilities::DefaultQuadrature< RangeSpaceType >::surfaceOrder(order); };
809
810 DiagonalAndNeighborStencil< DomainSpaceType, RangeSpaceType > stencil( jOp.domainSpace(), jOp.rangeSpace() );
811 jOp.reserve( stencil );
812 jOp.clear();
813
814 const std::size_t maxNumLocalDofs = jOp.domainSpace().blockMapper().maxNumDofs() * jOp.domainSpace().localBlockSize;
815 DomainValueVectorType phiIn = makeDomainValueVector( maxNumLocalDofs );
816 DomainValueVectorType phiOut = makeDomainValueVector( maxNumLocalDofs );
817
818 TemporaryLocalMatrixType jOpInIn( jOp.domainSpace(), jOp.rangeSpace() ), jOpOutIn( jOp.domainSpace(), jOp.rangeSpace() );
819 TemporaryLocalMatrixType jOpInOut( jOp.domainSpace(), jOp.rangeSpace() ), jOpOutOut( jOp.domainSpace(), jOp.rangeSpace() );
821
822 const auto &indexSet = gridPart().indexSet();
823 for( const EntityType &inside : elements( gridPart(), Partitions::interiorBorder ) )
824 {
825 uIn.bind( inside );
826
827 jOpInIn.init( inside, inside );
828 jOpInIn.clear();
829
830 if( integrands_.hasInterior() )
831 addLinearizedInteriorIntegral( uIn, phiIn, jOpInIn );
832
833 for( const auto &intersection : intersections( gridPart(), inside ) )
834 {
835 // check neighbor first since on periodic boundaries both,
836 // neighbor and boundary are true, so we treat neighbor first
837 if( intersection.neighbor() )
838 {
839 const EntityType &outside = intersection.outside();
840
841 jOpOutIn.init( outside, inside );
842 jOpOutIn.clear();
843
845 uOut.bind( outside );
846
847 if( outside.partitionType() != InteriorEntity )
848 addLinearizedSkeletonIntegral( intersection, uIn, uOut, phiIn, phiOut, jOpInIn, jOpOutIn );
849 else if( indexSet.index( inside ) < indexSet.index( outside ) )
850 {
851 jOpInOut.init( inside, outside );
852 jOpInOut.clear();
853 jOpOutOut.init( outside, outside );
854 jOpOutOut.clear();
855
856 addLinearizedSkeletonIntegral( intersection, uIn, uOut, phiIn, phiOut, jOpInIn, jOpOutIn, jOpInOut, jOpOutOut );
857
858 jOp.addLocalMatrix( inside, outside, jOpInOut );
859 jOp.addLocalMatrix( outside, outside, jOpOutOut );
860 }
861
862 jOp.addLocalMatrix( outside, inside, jOpOutIn );
863 }
864 else if( intersection.boundary() )
865 {
866 if( integrands_.hasBoundary() )
867 addLinearizedBoundaryIntegral( intersection, uIn, phiIn, jOpInIn );
868 }
869 }
870
871 jOp.addLocalMatrix( inside, inside, jOpInIn );
872 }
873 }
874
875 public:
876 template< class GridFunction, class JacobianOperator >
877 void assemble ( const GridFunction &u, JacobianOperator &jOp ) const
878 {
879 static_assert( std::is_same< typename GridFunction::GridPartType, GridPartType >::value, "Argument 'u' and Integrands must be defined on the same grid part." );
880 static_assert( std::is_same< typename JacobianOperator::DomainSpaceType::GridPartType, GridPartType >::value, "Argument 'jOp' and Integrands must be defined on the same grid part." );
881 static_assert( std::is_same< typename JacobianOperator::RangeSpaceType::GridPartType, GridPartType >::value, "Argument 'jOp' and Integrands must be defined on the same grid part." );
882
883 if( integrands_.hasSkeleton() )
884 assemble( u, jOp, std::true_type() );
885 else
886 assemble( u, jOp, std::false_type() );
887
888 // note: assembly done without local contributions so need
889 // to call flush assembly
890 jOp.flushAssembly();
891 }
892
893 // accessors
894
895 const GridPartType &gridPart () const { return gridPart_; }
896
897 unsigned int interiorQuadratureOrder(unsigned int order) const { return interiorQuadOrder_ == 0 ? defaultInteriorOrder_(order) : interiorQuadOrder_; }
898 unsigned int surfaceQuadratureOrder(unsigned int order) const { return surfaceQuadOrder_ == 0 ? defaultSurfaceOrder_ (order) : surfaceQuadOrder_; }
899
900 protected:
901 template <class U>
902 int maxOrder( const U& u ) const
903 {
904 return CallOrder< U > :: order( u );
905 }
906
907 template< class U, class W >
908 int maxOrder( const U& u, const W& w ) const
909 {
910 return std::max( maxOrder( u ), maxOrder( w ) );
911 }
912
913 template< class U, class V, class W >
914 int maxOrder( const U& u, const V& v, const W& w ) const
915 {
916 return std::max( maxOrder( u, v ), maxOrder( w ) );
917 }
918
919 template< class U, class V, class W, class X >
920 int maxOrder( const U& u, const V& v, const W& w, const X& x ) const
921 {
922 return std::max( maxOrder( u, v ), maxOrder( w, x) );
923 }
924
925 protected:
926 const GridPartType &gridPart_;
927 mutable IntegrandsType integrands_;
928
929 mutable std::function<int(const int)> defaultInteriorOrder_;
930 mutable std::function<int(const int)> defaultSurfaceOrder_;
931
932 unsigned int interiorQuadOrder_;
933 unsigned int surfaceQuadOrder_;
934
935 mutable std::vector< RangeValueVectorType > rangeValues_;
936 mutable RangeValueVectorType values_;
937 mutable DomainValueVectorType basisValues_;
938 mutable DomainValueVectorType domainValues_;
939
940 bool communicate_;
941 };
942
943 } // namespace Impl
944
945
946
947
948 // GalerkinOperator
949 // ----------------
950
951 template< class Integrands, class DomainFunction, class RangeFunction = DomainFunction >
953 : public virtual Operator< DomainFunction, RangeFunction >
954 {
955 typedef DomainFunction DomainFunctionType;
956 typedef RangeFunction RangeFunctionType;
957
958 static_assert( std::is_same< typename DomainFunctionType::GridPartType, typename RangeFunctionType::GridPartType >::value, "DomainFunction and RangeFunction must be defined on the same grid part." );
959
960 typedef typename RangeFunctionType::GridPartType GridPartType;
961
962 template< class... Args >
963 explicit GalerkinOperator ( const GridPartType &gridPart, Args &&... args )
964 : impl_( gridPart, std::forward< Args >( args )... )
965 {}
966
967 void setCommunicate( const bool communicate ) { impl_.setCommunicate( communicate ); }
968
969 void setQuadratureOrders(unsigned int interior, unsigned int surface) { impl_.setQuadratureOrders(interior,surface); }
970
971 virtual void operator() ( const DomainFunctionType &u, RangeFunctionType &w ) const final override
972 {
973 impl_.evaluate( u, w );
974 }
975
976 template< class GridFunction >
977 void operator() ( const GridFunction &u, RangeFunctionType &w ) const
978 {
979 return impl_.evaluate( u, w );
980 }
981
982 const GridPartType &gridPart () const { return impl_.gridPart(); }
983
984 typedef Integrands ModelType;
985 typedef Integrands DirichletModelType;
986 ModelType &model() const { return impl_.model(); }
987
988 protected:
989 Impl::GalerkinOperator< Integrands > impl_;
990 };
991
992
993
994 // DifferentiableGalerkinOperator
995 // ------------------------------
996
997 template< class Integrands, class JacobianOperator >
999 : public GalerkinOperator< Integrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType >,
1000 public DifferentiableOperator< JacobianOperator >
1001 {
1003
1004 public:
1005 typedef JacobianOperator JacobianOperatorType;
1006
1009 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType;
1010 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType;
1011
1013
1014 template< class... Args >
1016 const RangeDiscreteFunctionSpaceType &rSpace,
1017 Args &&... args )
1018 : BaseType( rSpace.gridPart(), std::forward< Args >( args )... ),
1019 dSpace_(dSpace), rSpace_(rSpace)
1020 {}
1021
1022 virtual void jacobian ( const DomainFunctionType &u, JacobianOperatorType &jOp ) const final override
1023 {
1024 impl_.assemble( u, jOp );
1025 }
1026
1027 template< class GridFunction >
1028 void jacobian ( const GridFunction &u, JacobianOperatorType &jOp ) const
1029 {
1030 impl_.assemble( u, jOp );
1031 }
1032
1034 {
1035 return dSpace_;
1036 }
1038 {
1039 return rSpace_;
1040 }
1041
1042 protected:
1043 using BaseType::impl_;
1046 };
1047
1048
1049
1050 // AutomaticDifferenceGalerkinOperator
1051 // -----------------------------------
1052
1053 template< class Integrands, class DomainFunction, class RangeFunction >
1055 : public GalerkinOperator< Integrands, DomainFunction, RangeFunction >,
1056 public AutomaticDifferenceOperator< DomainFunction, RangeFunction >
1057 {
1060
1061 public:
1063
1064 template< class... Args >
1065 explicit AutomaticDifferenceGalerkinOperator ( const GridPartType &gridPart, Args &&... args )
1066 : BaseType( gridPart, std::forward< Args >( args )... ), AutomaticDifferenceOperatorType()
1067 {}
1068 };
1069
1070
1071
1072 // ModelDifferentiableGalerkinOperator
1073 // -----------------------------------
1074
1075 template < class LinearOperator, class ModelIntegrands >
1077 : public DifferentiableGalerkinOperator< ModelIntegrands, LinearOperator >
1078 {
1080
1081 typedef typename ModelIntegrands::ModelType ModelType;
1082
1084 typedef typename LinearOperator::RangeSpaceType DiscreteFunctionSpaceType;
1085
1087 : BaseType( dfSpace.gridPart(), model )
1088 {}
1089
1090 template< class GridFunction >
1091 void apply ( const GridFunction &u, RangeFunctionType &w ) const
1092 {
1093 (*this)( u, w );
1094 }
1095
1096 template< class GridFunction >
1097 void apply ( const GridFunction &u, LinearOperator &jOp ) const
1098 {
1099 (*this).jacobian( u, jOp );
1100 }
1101 };
1102
1103 namespace Impl
1104 {
1105
1106 // GalerkinSchemeImpl
1107 // ------------------
1108 template <class O, bool addDirichletBC>
1109 struct DirichletBlockSelector { using type = void; };
1110 template <class O>
1111 struct DirichletBlockSelector<O,true> { using type = typename O::DirichletBlockVector; };
1112 template< class Integrands, class LinearOperator, class InverseOperator, bool addDirichletBC,
1113 template <class,class> class DifferentiableGalerkinOperatorImpl = DifferentiableGalerkinOperator >
1114 struct GalerkinSchemeImpl
1115 {
1116 typedef InverseOperator InverseOperatorType;
1117 typedef Integrands ModelType;
1118 using DifferentiableOperatorType = std::conditional_t< addDirichletBC,
1120 DifferentiableGalerkinOperatorImpl< Integrands, LinearOperator > >;
1121 using DirichletBlockVector = typename DirichletBlockSelector<
1123 DifferentiableGalerkinOperatorImpl< Integrands, LinearOperator >>,
1124 addDirichletBC>::type;
1125
1126 typedef typename DifferentiableOperatorType::DomainFunctionType DomainFunctionType;
1127 typedef typename DifferentiableOperatorType::RangeFunctionType RangeFunctionType;
1128 typedef typename DifferentiableOperatorType::JacobianOperatorType LinearOperatorType;
1129 typedef typename DifferentiableOperatorType::JacobianOperatorType JacobianOperatorType;
1130
1131 typedef RangeFunctionType DiscreteFunctionType;
1132 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeFunctionSpaceType;
1133 typedef typename RangeFunctionType::DiscreteFunctionSpaceType DomainFunctionSpaceType;
1134 typedef typename RangeFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
1135
1136 typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
1137 typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
1138
1140 typedef InverseOperator LinearInverseOperatorType;
1141 typedef typename NewtonOperatorType::ErrorMeasureType ErrorMeasureType;
1142
1143 struct SolverInfo
1144 {
1145 SolverInfo ( bool converged, int linearIterations, int nonlinearIterations )
1146 : converged( converged ), linearIterations( linearIterations ), nonlinearIterations( nonlinearIterations )
1147 {}
1148
1149 bool converged;
1150 int linearIterations, nonlinearIterations;
1151 };
1152
1153 GalerkinSchemeImpl ( const DiscreteFunctionSpaceType &dfSpace,
1154 const Integrands &integrands,
1155 const ParameterReader& parameter = Parameter::container() )
1156 : dfSpace_( dfSpace ),
1157 fullOperator_( dfSpace, dfSpace, std::move( integrands ) ),
1158 invOp_(parameter)
1159 {}
1160
1161 void setQuadratureOrders(unsigned int interior, unsigned int surface) { fullOperator().setQuadratureOrders(interior,surface); }
1162
1163 const DifferentiableOperatorType &fullOperator() const { return fullOperator_; }
1164 DifferentiableOperatorType &fullOperator() { return fullOperator_; }
1165
1166 void constraint ( DiscreteFunctionType &u ) const {}
1167
1168 template< class GridFunction >
1169 void operator() ( const GridFunction &u, DiscreteFunctionType &w ) const
1170 {
1171 fullOperator()( u, w );
1172 }
1173
1174 void setErrorMeasure(ErrorMeasureType &errorMeasure) const
1175 {
1176 invOp_.setErrorMeasure(errorMeasure);
1177 }
1178
1179 SolverInfo solve ( const DiscreteFunctionType &rhs, DiscreteFunctionType &solution ) const
1180 {
1181 DiscreteFunctionType rhs0 = rhs;
1182 setZeroConstraints( rhs0 );
1183 setModelConstraints( solution );
1184
1185 invOp_.bind(fullOperator());
1186 invOp_( rhs0, solution );
1187 invOp_.unbind();
1188 return SolverInfo( invOp_.converged(), invOp_.linearIterations(), invOp_.iterations() );
1189 }
1190
1191 SolverInfo solve ( DiscreteFunctionType &solution ) const
1192 {
1193 DiscreteFunctionType bnd( solution );
1194 bnd.clear();
1195 setModelConstraints( solution );
1196 invOp_.bind(fullOperator());
1197 invOp_( bnd, solution );
1198 invOp_.unbind();
1199 return SolverInfo( invOp_.converged(), invOp_.linearIterations(), invOp_.iterations() );
1200 }
1201
1202 template< class GridFunction >
1203 void jacobian( const GridFunction &ubar, LinearOperatorType &linearOp) const
1204 {
1205 fullOperator().jacobian( ubar, linearOp );
1206 }
1207
1208 const DiscreteFunctionSpaceType &space () const { return dfSpace_; }
1209 const GridPartType &gridPart () const { return space().gridPart(); }
1210 ModelType &model() const { return fullOperator().model(); }
1211
1212 void setConstraints( DomainFunctionType &u ) const
1213 {
1214 if constexpr (addDirichletBC)
1215 fullOperator().setConstraints( u );
1216 }
1217 void setConstraints( const typename DiscreteFunctionType::RangeType &value, DiscreteFunctionType &u ) const
1218 {
1219 if constexpr (addDirichletBC)
1220 fullOperator().setConstraints( value, u );
1221 }
1222 void setConstraints( const DiscreteFunctionType &u, DiscreteFunctionType &v ) const
1223 {
1224 if constexpr (addDirichletBC)
1225 fullOperator().setConstraints( u, v );
1226 }
1227 void subConstraints( const DiscreteFunctionType &u, DiscreteFunctionType &v ) const
1228 {
1229 if constexpr (addDirichletBC)
1230 fullOperator().subConstraints( u, v );
1231 }
1232 const auto& dirichletBlocks() const
1233 {
1234 if constexpr (addDirichletBC)
1235 return fullOperator().dirichletBlocks();
1236 }
1237
1238 protected:
1239 void setZeroConstraints( DiscreteFunctionType &u ) const
1240 {
1241 if constexpr (addDirichletBC)
1242 fullOperator().setConstraints( typename DiscreteFunctionType::RangeType(0), u );
1243 }
1244 void setModelConstraints( DiscreteFunctionType &u ) const
1245 {
1246 if constexpr (addDirichletBC)
1247 fullOperator().setConstraints( u );
1248 }
1249 const DiscreteFunctionSpaceType &dfSpace_;
1250 DifferentiableOperatorType fullOperator_;
1251 mutable NewtonOperatorType invOp_;
1252 };
1253
1254 } // end namespace Impl
1255
1256 // GalerkinScheme
1257 // --------------
1258
1259 template< class Integrands, class LinearOperator, class InverseOperator, bool addDirichletBC >
1260 using GalerkinScheme = Impl::GalerkinSchemeImpl< Integrands, LinearOperator, InverseOperator, addDirichletBC,
1262
1263 } // namespace Fem
1264
1265} // namespace Dune
1266
1267#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:1261
typename Impl::ConstLocalFunction< GridFunction >::Type ConstLocalFunction
Definition: const.hh:517
BasicParameterReader< std::function< const std::string *(const std::string &, const std::string *) > > ParameterReader
Definition: reader.hh:315
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
operator providing a Jacobian through automatic differentiation
Definition: automaticdifferenceoperator.hh:86
abstract differentiable operator
Definition: differentiableoperator.hh:29
JacobianOperator JacobianOperatorType
type of linear operator modelling the operator's Jacobian
Definition: differentiableoperator.hh:35
BaseType::DomainFunctionType DomainFunctionType
type of discrete function in the operator's domain
Definition: differentiableoperator.hh:38
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
Definition: dirichletwrapper.hh:27
Definition: galerkin.hh:954
GalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: galerkin.hh:963
const GridPartType & gridPart() const
Definition: galerkin.hh:982
ModelType & model() const
Definition: galerkin.hh:986
Impl::GalerkinOperator< Integrands > impl_
Definition: galerkin.hh:989
Integrands DirichletModelType
Definition: galerkin.hh:985
DomainFunction DomainFunctionType
Definition: galerkin.hh:955
RangeFunction RangeFunctionType
Definition: galerkin.hh:956
virtual void operator()(const DomainFunctionType &u, RangeFunctionType &w) const final override
Definition: galerkin.hh:971
Integrands ModelType
Definition: galerkin.hh:984
RangeFunctionType::GridPartType GridPartType
Definition: galerkin.hh:958
void setQuadratureOrders(unsigned int interior, unsigned int surface)
Definition: galerkin.hh:969
void setCommunicate(const bool communicate)
Definition: galerkin.hh:967
Definition: galerkin.hh:1001
BaseType::DomainFunctionType DomainFunctionType
Definition: galerkin.hh:1007
BaseType::GridPartType GridPartType
Definition: galerkin.hh:1012
void jacobian(const GridFunction &u, JacobianOperatorType &jOp) const
Definition: galerkin.hh:1028
RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType
Definition: galerkin.hh:1010
virtual void jacobian(const DomainFunctionType &u, JacobianOperatorType &jOp) const final override
obtain linearization
Definition: galerkin.hh:1022
BaseType::RangeFunctionType RangeFunctionType
Definition: galerkin.hh:1008
const RangeDiscreteFunctionSpaceType & rangeSpace() const
Definition: galerkin.hh:1037
DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: galerkin.hh:1009
const RangeDiscreteFunctionSpaceType & rSpace_
Definition: galerkin.hh:1045
JacobianOperator JacobianOperatorType
Definition: galerkin.hh:1005
const DomainDiscreteFunctionSpaceType & dSpace_
Definition: galerkin.hh:1044
const DomainDiscreteFunctionSpaceType & domainSpace() const
Definition: galerkin.hh:1033
DifferentiableGalerkinOperator(const DomainDiscreteFunctionSpaceType &dSpace, const RangeDiscreteFunctionSpaceType &rSpace, Args &&... args)
Definition: galerkin.hh:1015
BaseType::GridPartType GridPartType
Definition: galerkin.hh:1062
AutomaticDifferenceGalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: galerkin.hh:1065
ModelDifferentiableGalerkinOperator(ModelType &model, const DiscreteFunctionSpaceType &dfSpace)
Definition: galerkin.hh:1086
ModelIntegrands::ModelType ModelType
Definition: galerkin.hh:1081
LinearOperator::DomainFunctionType RangeFunctionType
Definition: galerkin.hh:1083
void apply(const GridFunction &u, RangeFunctionType &w) const
Definition: galerkin.hh:1091
DifferentiableGalerkinOperator< ModelIntegrands, LinearOperator > BaseType
Definition: galerkin.hh:1079
LinearOperator::RangeSpaceType DiscreteFunctionSpaceType
Definition: galerkin.hh:1084
void apply(const GridFunction &u, LinearOperator &jOp) const
Definition: galerkin.hh:1097
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