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 
25 #include <dune/fem/misc/l2norm.hh>
26 
32 
34 
35 // fempy includes
36 #include <dune/fempy/quadrature/fempyquadratures.hh>
37 
38 namespace 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
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
static constexpr T max(T a)
Definition: utility.hh:77
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 bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
static ParameterContainer & container()
Definition: io/parameter.hh:193
operator providing a Jacobian through automatic differentiation
Definition: automaticdifferenceoperator.hh:86
abstract differentiable operator
Definition: differentiableoperator.hh:29
BaseType::RangeFunctionType RangeFunctionType
type of discrete function in the operator's range
Definition: differentiableoperator.hh:40
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
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
const GridPartType & gridPart() const
Definition: galerkin.hh:982
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
const RangeDiscreteFunctionSpaceType & rangeSpace() const
Definition: galerkin.hh:1037
void jacobian(const GridFunction &u, JacobianOperatorType &jOp) const
Definition: galerkin.hh:1028
const DomainDiscreteFunctionSpaceType & domainSpace() const
Definition: galerkin.hh:1033
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
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
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