dune-fem  2.8-git
dgelliptic.hh
Go to the documentation of this file.
1 /**************************************************************************
2 
3  The dune-fem module is a module of DUNE (see www.dune-project.org).
4  It is based on the dune-grid interface library
5  extending the grid interface by a number of discretization algorithms
6  for solving non-linear systems of partial differential equations.
7 
8  Copyright (C) 2003 - 2015 Robert Kloefkorn
9  Copyright (C) 2003 - 2010 Mario Ohlberger
10  Copyright (C) 2004 - 2015 Andreas Dedner
11  Copyright (C) 2005 Adrian Burri
12  Copyright (C) 2005 - 2015 Mirko Kraenkel
13  Copyright (C) 2006 - 2015 Christoph Gersbacher
14  Copyright (C) 2006 - 2015 Martin Nolte
15  Copyright (C) 2011 - 2015 Tobias Malkmus
16  Copyright (C) 2012 - 2015 Stefan Girke
17  Copyright (C) 2013 - 2015 Claus-Justus Heine
18  Copyright (C) 2013 - 2014 Janick Gerstenberger
19  Copyright (C) 2013 Sven Kaulman
20  Copyright (C) 2013 Tom Ranner
21  Copyright (C) 2015 Marco Agnese
22  Copyright (C) 2015 Martin Alkaemper
23 
24 
25  The dune-fem module is free software; you can redistribute it and/or
26  modify it under the terms of the GNU General Public License as
27  published by the Free Software Foundation; either version 2 of
28  the License, or (at your option) any later version.
29 
30  The dune-fem module is distributed in the hope that it will be useful,
31  but WITHOUT ANY WARRANTY; without even the implied warranty of
32  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33  GNU General Public License for more details.
34 
35  You should have received a copy of the GNU General Public License along
36  with this program; if not, write to the Free Software Foundation, Inc.,
37  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
38 
39 **************************************************************************/
40 #ifndef DGELLIPTIC_HH
41 #define DGELLIPTIC_HH
42 
43 #include <dune/common/fmatrix.hh>
44 
50 
52 
53 // EllipticOperator
54 // ----------------
55 
56 template <class DFSpace>
58 {
59  DefaultPenalty(const DFSpace &space, double penalty)
60  : space_(space)
61  , penalty_(penalty)
62  {}
63  template <class Intersection>
64  double operator()(const Intersection &intersection,
65  double intersectionArea, double area, double nbArea) const
66  {
67  const double hInv = intersectionArea / std::min( area, nbArea );
68  return penalty_ * hInv;
69  }
70  const double &factor() const { return penalty_; }
71  private:
72  const DFSpace &space_;
73  double penalty_;
74 };
75 
76 template< class DiscreteFunction, class Model,
79 : public virtual Dune::Fem::Operator< DiscreteFunction >
80 {
81  typedef DiscreteFunction DiscreteFunctionType;
82  typedef Model ModelType;
83 
84  typedef DiscreteFunction DomainDiscreteFunctionType;
85  typedef DiscreteFunction RangeDiscreteFunctionType;
86 protected:
87  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
88  typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
89  typedef typename LocalFunctionType::RangeType RangeType;
90  typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
91 
92  typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
93  typedef typename IteratorType::Entity EntityType;
94  typedef typename EntityType::Geometry GeometryType;
95 
96  typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
97 
98  typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
99  typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
100  typedef typename IntersectionIteratorType::Intersection IntersectionType;
101  typedef typename IntersectionType::Geometry IntersectionGeometryType;
102 
105 
106  static const int dimDomain = LocalFunctionType::dimDomain;
107  static const int dimRange = LocalFunctionType::dimRange;
108 
109 public:
112  const ModelType &model,
114  : model_( model ),
115  penalty_( space, parameter.getValue< double >( "penalty", 40 ) )
116  {
117  // std::cout << "dg operator with penalty:" << penalty_.factor() << std::endl;
118  }
119 
121  const DiscreteFunctionSpaceType &rSpace,
122  const ModelType &model,
124  : model_( model ),
125  penalty_( dSpace, parameter.getValue< double >( "penalty", 40 ) )
126  {
127  // std::cout << "dg operator with penalty:" << penalty_.factor() << std::endl;
128  }
129 
132  { apply(u,w); }
133  template <class GF>
134  void apply( const GF &u, RangeDiscreteFunctionType &w ) const;
135 
136 protected:
137  const ModelType &model () const { return model_; }
138  Penalty penalty() const { return penalty_; }
139 
140 private:
141  const ModelType &model_;
142  Penalty penalty_;
143 };
144 
145 // DifferentiableDGEllipticOperator
146 // ------------------------------
147 
148 template< class JacobianOperator, class Model,
151 : public DGEllipticOperator< typename JacobianOperator::DomainFunctionType, Model, Penalty >,
152  public Dune::Fem::DifferentiableOperator< JacobianOperator >
153 {
155 
156  typedef JacobianOperator JacobianOperatorType;
157 
159  typedef typename BaseType::ModelType ModelType;
162 
163 protected:
164  typedef typename DomainDiscreteFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType;
165  typedef typename DomainDiscreteFunctionType::LocalFunctionType DomainLocalFunctionType;
166  typedef typename DomainLocalFunctionType::RangeType DomainRangeType;
167  typedef typename DomainLocalFunctionType::JacobianRangeType DomainJacobianRangeType;
168  typedef typename RangeDiscreteFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType;
169  typedef typename RangeDiscreteFunctionType::LocalFunctionType RangeLocalFunctionType;
170  typedef typename RangeLocalFunctionType::RangeType RangeRangeType;
171  typedef typename RangeLocalFunctionType::JacobianRangeType RangeJacobianRangeType;
172 
173  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
174  typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
175  typedef typename LocalFunctionType::RangeType RangeType;
176  typedef typename LocalFunctionType::RangeFieldType RangeFieldType;
177  typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
178 
179  typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
180  typedef typename IteratorType::Entity EntityType;
181  typedef typename EntityType::Geometry GeometryType;
182 
183  typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
184 
185  typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
186  typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
187  typedef typename IntersectionIteratorType::Intersection IntersectionType;
188  typedef typename IntersectionType::Geometry IntersectionGeometryType;
189 
192 
193  static const int dimDomain = LocalFunctionType::dimDomain;
194  static const int dimRange = LocalFunctionType::dimRange;
195 
196 public:
199  const ModelType &model,
201  : BaseType( space, model, parameter )
202  {}
204  const DiscreteFunctionSpaceType &rSpace,
205  const ModelType &model,
207  : BaseType( dSpace, rSpace, model, parameter )
208  {}
209 
212  { apply(u,jOp); }
213  template <class GridFunctionType>
214  void apply ( const GridFunctionType &u, JacobianOperatorType &jOp ) const;
215  using BaseType::apply;
216 
217 protected:
218  using BaseType::model;
219  using BaseType::penalty;
220 };
221 
222 // Implementation of DGEllipticOperator
223 // ----------------------------------
224 
225 template< class RangeDiscreteFunction, class Model, class Penalty >
226 template<class GF>
228  ::apply( const GF &u, RangeDiscreteFunctionType &w ) const
229 {
230  // clear destination
231  w.clear();
232 
233  // get discrete function space
234  const DiscreteFunctionSpaceType &dfSpace = w.space();
235 
236  // iterate over grid
237  const IteratorType end = dfSpace.end();
238  for( IteratorType it = dfSpace.begin(); it != end; ++it )
239  {
240  // get entity (here element)
241  const EntityType &entity = *it;
242 
243  bool needsCalculation = model().init( entity );
244  if (! needsCalculation )
245  continue;
246 
247  // get elements geometry
248  const GeometryType &geometry = entity.geometry();
249 
250  // get local representation of the discrete functions
251  const auto uLocal = u.localFunction( entity );
252  LocalFunctionType wLocal = w.localFunction( entity );
253 
254  // obtain quadrature order
255  const int quadOrder = uLocal.order() + wLocal.order();
256 
257  { // element integral
258  QuadratureType quadrature( entity, quadOrder );
259  const size_t numQuadraturePoints = quadrature.nop();
260  for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
261  {
262  const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
263  const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );
264 
265  RangeType vu;
266  uLocal.evaluate( quadrature[ pt ], vu );
268  uLocal.jacobian( quadrature[ pt ], du );
269 
270  // compute mass contribution (studying linear case so linearizing around zero)
271  RangeType avu( 0 );
272  model().source( quadrature[ pt ], vu, du, avu );
273  avu *= weight;
274  // add to local functional wLocal.axpy( quadrature[ pt ], avu );
275 
276  JacobianRangeType adu( 0 );
277  // apply diffusive flux
278  model().diffusiveFlux( quadrature[ pt ], vu, du, adu );
279  adu *= weight;
280 
281  // add to local function
282  wLocal.axpy( quadrature[ pt ], avu, adu );
283  }
284  }
285  if ( ! dfSpace.continuous() )
286  {
287  const double area = entity.geometry().volume();
288  const IntersectionIteratorType iitend = dfSpace.gridPart().iend( entity );
289  for( IntersectionIteratorType iit = dfSpace.gridPart().ibegin( entity ); iit != iitend; ++iit ) // looping over intersections
290  {
292  const IntersectionType &intersection = *iit;
293  if ( intersection.neighbor() )
294  {
295  // const EntityType outside = Dune::Fem::make_entity( intersection.outside() );
296  const EntityType outside = intersection.outside() ;
297 
298  typedef typename IntersectionType::Geometry IntersectionGeometryType;
299  const IntersectionGeometryType &intersectionGeometry = intersection.geometry();
300 
301  // compute penalty factor
302  const double intersectionArea = intersectionGeometry.volume();
303  const double beta = penalty()(intersection,intersectionArea,area,outside.geometry().volume());
304 
305  auto uOutLocal = u.localFunction( outside ); // local u on outisde element
306 
307  FaceQuadratureType quadInside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
308  FaceQuadratureType quadOutside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::OUTSIDE );
309  const size_t numQuadraturePoints = quadInside.nop();
311 
312  for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
313  {
315  // get coordinate of quadrature point on the reference element of the intersection
316  const typename FaceQuadratureType::LocalCoordinateType &x = quadInside.localPoint( pt );
317  DomainType normal = intersection.integrationOuterNormal( x );
318  double faceVol = normal.two_norm();
319  normal /= faceVol; // make it into a unit normal
320  const double weight = quadInside.weight( pt ) * faceVol;
321 
322  RangeType vuIn,vuOut;
323  JacobianRangeType duIn, duOut;
324  uLocal.evaluate( quadInside[ pt ], vuIn );
325  uLocal.jacobian( quadInside[ pt ], duIn );
326  uOutLocal.evaluate( quadOutside[ pt ], vuOut );
327  uOutLocal.jacobian( quadOutside[ pt ], duOut );
328  RangeType jump = vuIn - vuOut;
329  // compute -0.5 * [u] x normal
330  JacobianRangeType dvalue;
331  for (int r=0;r<dimRange;++r)
332  for (int d=0;d<dimDomain;++d)
333  dvalue[r][d] = -0.5 * normal[d] * jump[r];
334  JacobianRangeType aduIn,aduOut;
335  model().init( outside );
336  model().diffusiveFlux( quadOutside[ pt ], vuOut, duOut, aduOut );
337  auto pFactor = model().penalty( quadOutside[ pt ], vuOut );
338  model().init( entity );
339  model().diffusiveFlux( quadInside[ pt ], vuIn, duIn, aduIn );
340  JacobianRangeType affine;
341  model().diffusiveFlux( quadInside[ pt ], RangeType(0), JacobianRangeType(0), affine);
342  pFactor += model().penalty( quadInside[ pt ], vuIn );
344 
346  RangeType value;
347  JacobianRangeType advalue;
348  // penalty term : beta [u] [phi] = beta (u+ - u-)(phi+ - phi-)=beta (u+ - u-)phi+
349  value = jump;
350  for (unsigned int r=0;r<dimRange;++r)
351  value[r] *= beta * pFactor[r]/2.;
352  // {A grad u}.[phi] = {A grad u}.phi+ n_+ = 0.5*(grad u+ + grad u-).n_+ phi+
353  aduIn += aduOut;
354  aduIn *= -0.5;
355  aduIn.umv(normal,value);
356  // [ u ] * { {A} grad phi_en } = -normal(u+ - u-) * 0.5 {A} grad phi_en
357  // we actually need G(u)tau with G(x,u) = d/sigma_j D_i(x,u,sigma)
358  // - we might need to assume D(x,u,sigma) = G(x,u)sigma + affine(x)
359  model().diffusiveFlux( quadInside[ pt ], vuIn, dvalue, advalue );
360  // advalue -= affine;
361 
362  value *= weight;
363  advalue *= weight;
364  wLocal.axpy( quadInside[ pt ], value, advalue );
366  }
367  }
368  else if( intersection.boundary() )
369  {
370  std::array<int,dimRange> components;
371  components.fill(0);
372  model().isDirichletIntersection( intersection, components);
373 
374  typedef typename IntersectionType::Geometry IntersectionGeometryType;
375  const IntersectionGeometryType &intersectionGeometry = intersection.geometry();
376 
377  // compute penalty factor
378  const double intersectionArea = intersectionGeometry.volume();
379  const double beta = penalty()(intersection,intersectionArea,area,area);
380 
381  FaceQuadratureType quadInside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
382  const size_t numQuadraturePoints = quadInside.nop();
383 
384  for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
385  {
386  const typename FaceQuadratureType::LocalCoordinateType &x = quadInside.localPoint( pt );
387  const DomainType normal = intersection.integrationOuterNormal( x );
388  const double weight = quadInside.weight( pt );
389 
390  RangeType bndValue;
391  model().dirichlet(1, quadInside[pt], bndValue);
392 
393  RangeType value;
394  JacobianRangeType dvalue,advalue;
395 
396  RangeType vuIn,jump;
397  JacobianRangeType duIn, aduIn;
398  uLocal.evaluate( quadInside[ pt ], vuIn );
399  uLocal.jacobian( quadInside[ pt ], duIn );
400  jump = vuIn;
401  jump -= bndValue;
402 
403  // penalty term : beta [u] [phi] = beta (u+ - u-)(phi+ - phi-)=beta (u+ - u-)phi+
404  auto pFactor = model().penalty( quadInside[ pt ], vuIn );
405  value = jump;
406  for (unsigned int r=0;r<dimRange;++r)
407  value[r] *= beta * pFactor[r] * intersectionGeometry.integrationElement( x );
408 
409  // [ u ] * { grad phi_en } = -normal(u+ - u-) * 0.5 grad phi_en
410  // here we need a diadic product of u x n
411  for (int r=0;r<dimRange;++r)
412  for (int d=0;d<dimDomain;++d)
413  dvalue[r][d] = -0.5 * normal[d] * jump[r];
414 
415  model().diffusiveFlux( quadInside[ pt ], jump, dvalue, advalue );
416 
417  // consistency term
418  // {A grad u}.[phi] = {A grad u}.phi+ n_+ = 0.5*(grad u+ + grad u-).n_+ phi+
419  model().diffusiveFlux( quadInside[ pt ], vuIn, duIn, aduIn );
420  aduIn.mmv(normal,value);
421 
422  for (int r=0;r<dimRange;++r)
423  if (!components[r]) // do not use dirichlet constraints here
424  {
425  value[r] = 0;
426  advalue[r] = 0;
427  }
428 
429  value *= weight;
430  advalue *= weight;
431  wLocal.axpy( quadInside[ pt ], value, advalue );
432  }
433  }
434  }
435  }
436 
437  }
438 
439  // communicate data (in parallel runs)
440  w.communicate();
441 }
442 
443 // Implementation of DifferentiableDGEllipticOperator
444 // ------------------------------------------------
445 template< class JacobianOperator, class Model, class Penalty >
446 template<class GF>
448  ::apply ( const GF &u, JacobianOperator &jOp ) const
449 {
450  // std::cout << "starting assembly\n";
451  // Dune::Timer timer;
452  typedef typename JacobianOperator::LocalMatrixType LocalMatrixType;
453  typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType BasisFunctionSetType;
454 
455  const DomainDiscreteFunctionSpaceType &domainSpace = jOp.domainSpace();
456  const RangeDiscreteFunctionSpaceType &rangeSpace = jOp.rangeSpace();
457 
459  jOp.reserve(stencil);
460  jOp.clear();
461 
462  const GridPartType& gridPart = rangeSpace.gridPart();
463 
464  const unsigned int numDofs = rangeSpace.blockMapper().maxNumDofs() *
465  DiscreteFunctionSpaceType :: localBlockSize ;
466 
467  std::vector< RangeType > phi( numDofs );
468  std::vector< JacobianRangeType > dphi( numDofs );
469 
470  std::vector< RangeType > phiNb( numDofs );
471  std::vector< JacobianRangeType > dphiNb( numDofs );
472 
473  // std::cout << " in assembly: start element loop size=" << rangeSpace.gridPart().grid().size(0) << " time= " << timer.elapsed() << std::endl;;
474  const IteratorType end = rangeSpace.end();
475  for( IteratorType it = rangeSpace.begin(); it != end; ++it )
476  {
477  const EntityType &entity = *it;
478 
479  bool needsCalculation = model().init( entity );
480  if (! needsCalculation )
481  continue;
482 
483  const GeometryType geometry = entity.geometry();
484 
485  const auto uLocal = u.localFunction( entity );
486 
487  LocalMatrixType jLocal = jOp.localMatrix( entity, entity );
488 
489  const BasisFunctionSetType &baseSet = jLocal.domainBasisFunctionSet();
490  const unsigned int numBaseFunctions = baseSet.size();
491 
492  QuadratureType quadrature( entity, 2*rangeSpace.order() );
493  const size_t numQuadraturePoints = quadrature.nop();
494  for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
495  {
496  const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
497  const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );
498 
499  // evaluate all basis functions at given quadrature point
500  baseSet.evaluateAll( quadrature[ pt ], phi );
501 
502  // evaluate jacobians of all basis functions at given quadrature point
503  baseSet.jacobianAll( quadrature[ pt ], dphi );
504 
505  // get value for linearization
506  RangeType u0;
507  JacobianRangeType jacU0;
508  uLocal.evaluate( quadrature[ pt ], u0 );
509  uLocal.jacobian( quadrature[ pt ], jacU0 );
510 
511  RangeType aphi( 0 );
512  JacobianRangeType adphi( 0 );
513  for( unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol )
514  {
515  // if mass terms or right hand side is present
516  model().linSource( u0, jacU0, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], aphi );
517 
518  // if gradient term is present
519  model().linDiffusiveFlux( u0, jacU0, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], adphi );
520 
521  // get column object and call axpy method
522  jLocal.column( localCol ).axpy( phi, dphi, aphi, adphi, weight );
523  }
524  }
525  if ( rangeSpace.continuous() )
526  continue;
527 
528  double area = geometry.volume();
529  const IntersectionIteratorType endiit = gridPart.iend( entity );
530  for ( IntersectionIteratorType iit = gridPart.ibegin( entity );
531  iit != endiit ; ++ iit )
532  {
533  const IntersectionType& intersection = *iit ;
534 
535  if( intersection.neighbor() )
536  {
537  const EntityType neighbor = Dune::Fem::make_entity( intersection.outside() );
538 
539  typedef typename IntersectionType::Geometry IntersectionGeometryType;
540  const IntersectionGeometryType &intersectionGeometry = intersection.geometry();
541 
543  // get local matrix for face entries
544  LocalMatrixType localOpNb = jOp.localMatrix( neighbor, entity );
545  // get neighbor's base function set
546  const BasisFunctionSetType &baseSetNb = localOpNb.domainBasisFunctionSet();
548 
549  // compute penalty factor
550  const double intersectionArea = intersectionGeometry.volume();
551  const double beta = penalty()(intersection,intersectionArea,area,neighbor.geometry().volume());
552 
553  // here we assume that the intersection is conforming
554  FaceQuadratureType faceQuadInside(gridPart, intersection, 2*rangeSpace.order() + 1,
555  FaceQuadratureType::INSIDE);
556  FaceQuadratureType faceQuadOutside(gridPart, intersection, 2*rangeSpace.order() + 1,
557  FaceQuadratureType::OUTSIDE);
558 
559  const size_t numFaceQuadPoints = faceQuadInside.nop();
560  for( size_t pt = 0; pt < numFaceQuadPoints; ++pt )
561  {
562  const typename FaceQuadratureType::LocalCoordinateType &x = faceQuadInside.localPoint( pt );
563  DomainType normal = intersection.integrationOuterNormal( x );
564  double faceVol = normal.two_norm();
565  normal /= faceVol; // make it into a unit normal
566 
567  const double quadWeight = faceQuadInside.weight( pt );
568  const double weight = quadWeight * faceVol;
569 
571  RangeType u0En;
572  JacobianRangeType u0EnJac;
573  uLocal.evaluate( faceQuadInside[ pt ], u0En );
574  uLocal.jacobian( faceQuadInside[ pt ], u0EnJac );
575 
577  // evaluate basis function of face inside E^- (entity)
579 
580  // evaluate all basis functions for quadrature point pt
581  baseSet.evaluateAll( faceQuadInside[ pt ], phi );
582 
583  // evaluate the jacobians of all basis functions
584  baseSet.jacobianAll( faceQuadInside[ pt ], dphi );
585 
587  // evaluate basis function of face inside E^+ (neighbor)
589 
590  // evaluate all basis functions for quadrature point pt on neighbor
591  baseSetNb.evaluateAll( faceQuadOutside[ pt ], phiNb );
592 
593  // evaluate the jacobians of all basis functions on neighbor
594  baseSetNb.jacobianAll( faceQuadOutside[ pt ], dphiNb );
595 
596  model().init( entity );
597  for( unsigned int i = 0; i < numBaseFunctions; ++i )
598  {
599  JacobianRangeType adphiEn = dphi[ i ];
600  model().linDiffusiveFlux( u0En, u0EnJac, faceQuadInside[ pt ], phi[i], adphiEn, dphi[ i ] );
601  }
602  auto pFactor = model().penalty( faceQuadInside[ pt ], u0En );
603  model().init( neighbor );
604  for( unsigned int i = 0; i < numBaseFunctions; ++i )
605  {
606  JacobianRangeType adphiNb = dphiNb[ i ];
607  model().linDiffusiveFlux( u0En, u0EnJac, faceQuadOutside[ pt ], phiNb[i], adphiNb, dphiNb[ i ] );
608  }
609  pFactor += model().penalty( faceQuadOutside[ pt ], u0En );
610  model().init( entity );
612 
614  for( unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol )
615  {
616  RangeType valueEn(0), valueNb(0);
617  JacobianRangeType dvalueEn(0), dvalueNb(0);
618 
619  // -{ A grad u } * [ phi_en ]
620  dphi[localCol].usmv( -0.5, normal, valueEn );
621 
622  // -{ A grad u } * [ phi_en ]
623  dphiNb[localCol].usmv( -0.5, normal, valueNb );
624 
625  // [ u ] * [ phi_en ] = u^- * phi_en^-
626  for (unsigned int r=0;r<dimRange;++r)
627  {
628  valueEn[r] += beta*pFactor[r]/2.*phi[localCol][r];
629  valueNb[r] -= beta*pFactor[r]/2.*phiNb[localCol][r];
630  }
631  // here we need a diadic product of u x n
632  for ( int r=0; r< dimRange; ++r )
633  for ( int d=0; d< dimDomain; ++d )
634  {
635  // [ u ] * { grad phi_en }
636  dvalueEn[r][d] = - 0.5 * normal[d] * phi[localCol][r];
637 
638  // [ u ] * { grad phi_en }
639  dvalueNb[r][d] = 0.5 * normal[d] * phiNb[localCol][r];
640  }
641 
642  jLocal.column( localCol ).axpy( phi, dphi, valueEn, dvalueEn, weight );
643  localOpNb.column( localCol ).axpy( phi, dphi, valueNb, dvalueNb, weight );
644  }
646  }
647  }
648  else if( intersection.boundary() )
649  {
650  std::array<int,dimRange> components;
651  components.fill(0);
652  model().isDirichletIntersection( intersection, components);
653 
654  typedef typename IntersectionType::Geometry IntersectionGeometryType;
655  const IntersectionGeometryType &intersectionGeometry = intersection.geometry();
656 
657  // compute penalty factor
658  const double intersectionArea = intersectionGeometry.volume();
659  const double beta = penalty()(intersection,intersectionArea,area,area);
660 
661  // here we assume that the intersection is conforming
662  FaceQuadratureType faceQuadInside(gridPart, intersection, 2*rangeSpace.order() + 1,
663  FaceQuadratureType::INSIDE);
664 
665  const size_t numFaceQuadPoints = faceQuadInside.nop();
666  for( size_t pt = 0; pt < numFaceQuadPoints; ++pt )
667  {
668  const typename FaceQuadratureType::LocalCoordinateType &x = faceQuadInside.localPoint( pt );
669  DomainType normal = intersection.integrationOuterNormal( x );
670  double faceVol = normal.two_norm();
671  normal /= faceVol; // make it into a unit normal
672 
673  const double quadWeight = faceQuadInside.weight( pt );
674  const double weight = quadWeight * faceVol;
675 
676  RangeType u0En;
677  JacobianRangeType u0EnJac;
678  uLocal.evaluate( faceQuadInside[ pt ], u0En );
679  uLocal.jacobian( faceQuadInside[ pt ], u0EnJac );
680 
682  // evaluate basis function of face inside E^- (entity)
684 
685  // evaluate all basis functions for quadrature point pt
686  baseSet.evaluateAll( faceQuadInside[ pt ], phi );
687 
688  // evaluate the jacobians of all basis functions
689  baseSet.jacobianAll( faceQuadInside[ pt ], dphi );
690 
691  for( unsigned int i = 0; i < numBaseFunctions; ++i )
692  {
693  JacobianRangeType adphiEn = dphi[ i ];
694  model().linDiffusiveFlux( u0En, u0EnJac, faceQuadInside[ pt ], phi[i], adphiEn, dphi[ i ] );
695  }
696 
697  auto pFactor = model().penalty( faceQuadInside[ pt ], u0En );
698 
699  for( unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol )
700  {
701  RangeType valueEn(0);
702  JacobianRangeType dvalueEn(0);
703 
704  // -{ A grad u } * [ phi_en ]
705  dphi[localCol].usmv( -1.0, normal, valueEn );
706 
707  // [ u ] * [ phi_en ] = u^- * phi_en^-
708  for (unsigned int r=0;r<dimRange;++r)
709  valueEn[r] += beta*pFactor[r]*phi[ localCol ][r];
710 
711  // here we need a diadic product of u x n
712  for ( int r=0; r< dimRange; ++r )
713  for ( int d=0; d< dimDomain; ++d )
714  {
715  // [ u ] * { grad phi_en }
716  dvalueEn[r][d] = - 0.5 * normal[d] * phi[localCol][r];
717  }
718 
719  for (int r=0;r<dimRange;++r)
720  if (!components[r]) // do not use dirichlet constraints here
721  {
722  valueEn[r] = 0;
723  dvalueEn[r] = 0;
724  }
725 
726  jLocal.column( localCol ).axpy( phi, dphi, valueEn, dvalueEn, weight );
727  }
728  }
729  }
730  }
731  } // end grid traversal
732  jOp.flushAssembly();
733  // std::cout << " in assembly: final " << timer.elapsed() << std::endl;;
734 }
735 
736 #endif // ELLIPTIC_HH
Dune::Entity< codim, dim, Grid, Implementation > make_entity(Dune::Entity< codim, dim, Grid, Implementation > entity)
Definition: compatibility.hh:21
static constexpr T min(T a)
Definition: utility.hh:93
static ParameterContainer & container()
Definition: io/parameter.hh:193
abstract differentiable operator
Definition: differentiableoperator.hh:29
JacobianOperator JacobianOperatorType
type of linear operator modelling the operator's Jacobian
Definition: differentiableoperator.hh:35
abstract operator
Definition: operator.hh:34
Stencil contaning the entries (en,en) and (en,nb) for all entities en in the space and neighbors nb o...
Definition: stencil.hh:336
quadrature on the codim-0 reference element
Definition: elementquadrature.hh:18
Definition: dgelliptic.hh:58
const double & factor() const
Definition: dgelliptic.hh:70
DefaultPenalty(const DFSpace &space, double penalty)
Definition: dgelliptic.hh:59
double operator()(const Intersection &intersection, double intersectionArea, double area, double nbArea) const
Definition: dgelliptic.hh:64
Definition: dgelliptic.hh:80
DiscreteFunctionSpaceType::IteratorType IteratorType
Definition: dgelliptic.hh:92
const ModelType & model() const
Definition: dgelliptic.hh:137
DGEllipticOperator(const DiscreteFunctionSpaceType &space, const ModelType &model, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
constructor
Definition: dgelliptic.hh:111
DiscreteFunctionSpaceType::DomainType DomainType
Definition: dgelliptic.hh:96
GridPartType::IntersectionIteratorType IntersectionIteratorType
Definition: dgelliptic.hh:99
DGEllipticOperator(const DiscreteFunctionSpaceType &dSpace, const DiscreteFunctionSpaceType &rSpace, const ModelType &model, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: dgelliptic.hh:120
static const int dimDomain
Definition: dgelliptic.hh:106
IteratorType::Entity EntityType
Definition: dgelliptic.hh:93
EntityType::Geometry GeometryType
Definition: dgelliptic.hh:94
DiscreteFunctionSpaceType::GridPartType GridPartType
Definition: dgelliptic.hh:98
void apply(const GF &u, RangeDiscreteFunctionType &w) const
Definition: dgelliptic.hh:228
Penalty penalty() const
Definition: dgelliptic.hh:138
IntersectionType::Geometry IntersectionGeometryType
Definition: dgelliptic.hh:101
IntersectionIteratorType::Intersection IntersectionType
Definition: dgelliptic.hh:100
LocalFunctionType::JacobianRangeType JacobianRangeType
Definition: dgelliptic.hh:90
virtual void operator()(const DomainDiscreteFunctionType &u, RangeDiscreteFunctionType &w) const
application operator
Definition: dgelliptic.hh:131
static const int dimRange
Definition: dgelliptic.hh:107
DiscreteFunction DomainDiscreteFunctionType
Definition: dgelliptic.hh:84
Dune::Fem::ElementQuadrature< GridPartType, 1 > FaceQuadratureType
Definition: dgelliptic.hh:103
LocalFunctionType::RangeType RangeType
Definition: dgelliptic.hh:89
Dune::Fem::CachingQuadrature< GridPartType, 0 > QuadratureType
Definition: dgelliptic.hh:104
DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
Definition: dgelliptic.hh:87
DiscreteFunction RangeDiscreteFunctionType
Definition: dgelliptic.hh:85
DiscreteFunctionType::LocalFunctionType LocalFunctionType
Definition: dgelliptic.hh:88
DiscreteFunction DiscreteFunctionType
Definition: dgelliptic.hh:81
Model ModelType
Definition: dgelliptic.hh:82
Definition: dgelliptic.hh:153
RangeLocalFunctionType::JacobianRangeType RangeJacobianRangeType
Definition: dgelliptic.hh:171
static const int dimRange
Definition: dgelliptic.hh:194
DiscreteFunctionSpaceType::IteratorType IteratorType
Definition: dgelliptic.hh:179
static const int dimDomain
Definition: dgelliptic.hh:193
DiscreteFunctionSpaceType::DomainType DomainType
Definition: dgelliptic.hh:183
DomainLocalFunctionType::JacobianRangeType DomainJacobianRangeType
Definition: dgelliptic.hh:167
BaseType::ModelType ModelType
Definition: dgelliptic.hh:159
DiscreteFunctionType::LocalFunctionType LocalFunctionType
Definition: dgelliptic.hh:174
DifferentiableDGEllipticOperator(const DiscreteFunctionSpaceType &dSpace, const DiscreteFunctionSpaceType &rSpace, const ModelType &model, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
Definition: dgelliptic.hh:203
DomainDiscreteFunctionType::LocalFunctionType DomainLocalFunctionType
Definition: dgelliptic.hh:165
IntersectionType::Geometry IntersectionGeometryType
Definition: dgelliptic.hh:188
Dune::Fem::ElementQuadrature< GridPartType, 1 > FaceQuadratureType
Definition: dgelliptic.hh:190
void apply(const GridFunctionType &u, JacobianOperatorType &jOp) const
IteratorType::Entity EntityType
Definition: dgelliptic.hh:180
DiscreteFunctionSpaceType::GridPartType GridPartType
Definition: dgelliptic.hh:185
DomainDiscreteFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: dgelliptic.hh:164
LocalFunctionType::JacobianRangeType JacobianRangeType
Definition: dgelliptic.hh:177
DGEllipticOperator< typename JacobianOperator::DomainFunctionType, Model, Penalty > BaseType
Definition: dgelliptic.hh:154
LocalFunctionType::RangeType RangeType
Definition: dgelliptic.hh:175
DifferentiableDGEllipticOperator(const DiscreteFunctionSpaceType &space, const ModelType &model, const Dune::Fem::ParameterReader &parameter=Dune::Fem::Parameter::container())
constructor
Definition: dgelliptic.hh:198
RangeDiscreteFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType
Definition: dgelliptic.hh:168
GridPartType::IntersectionIteratorType IntersectionIteratorType
Definition: dgelliptic.hh:186
DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
Definition: dgelliptic.hh:173
RangeDiscreteFunctionType::LocalFunctionType RangeLocalFunctionType
Definition: dgelliptic.hh:169
BaseType::DiscreteFunctionType DiscreteFunctionType
Definition: dgelliptic.hh:158
BaseType::DomainDiscreteFunctionType DomainDiscreteFunctionType
Definition: dgelliptic.hh:160
Dune::Fem::CachingQuadrature< GridPartType, 0 > QuadratureType
Definition: dgelliptic.hh:191
RangeLocalFunctionType::RangeType RangeRangeType
Definition: dgelliptic.hh:170
void jacobian(const DomainDiscreteFunctionType &u, JacobianOperatorType &jOp) const
method to setup the jacobian of the operator for storage in a matrix
Definition: dgelliptic.hh:211
IntersectionIteratorType::Intersection IntersectionType
Definition: dgelliptic.hh:187
JacobianOperator JacobianOperatorType
Definition: dgelliptic.hh:156
BaseType::RangeDiscreteFunctionType RangeDiscreteFunctionType
Definition: dgelliptic.hh:161
LocalFunctionType::RangeFieldType RangeFieldType
Definition: dgelliptic.hh:176
DomainLocalFunctionType::RangeType DomainRangeType
Definition: dgelliptic.hh:166
EntityType::Geometry GeometryType
Definition: dgelliptic.hh:181