43#include <dune/common/fmatrix.hh>
56template <
class DFSpace>
63 template <
class Intersection>
65 double intersectionArea,
double area,
double nbArea)
const
67 const double hInv = intersectionArea /
std::min( area, nbArea );
68 return penalty_ * hInv;
70 const double &
factor()
const {
return penalty_; }
72 const DFSpace &space_;
76template<
class DiscreteFunction,
class Model,
89 typedef typename LocalFunctionType::RangeType
RangeType;
92 typedef typename DiscreteFunctionSpaceType::IteratorType
IteratorType;
96 typedef typename DiscreteFunctionSpaceType::DomainType
DomainType;
98 typedef typename DiscreteFunctionSpaceType::GridPartType
GridPartType;
106 static const int dimDomain = LocalFunctionType::dimDomain;
107 static const int dimRange = LocalFunctionType::dimRange;
115 penalty_( space, parameter.getValue< double >(
"penalty", 40 ) )
125 penalty_( dSpace, parameter.getValue< double >(
"penalty", 40 ) )
148template<
class JacobianOperator,
class Model,
151:
public DGEllipticOperator< typename JacobianOperator::DomainFunctionType, Model, Penalty >,
175 typedef typename LocalFunctionType::RangeType
RangeType;
183 typedef typename DiscreteFunctionSpaceType::DomainType
DomainType;
193 static const int dimDomain = LocalFunctionType::dimDomain;
194 static const int dimRange = LocalFunctionType::dimRange;
213 template <
class Gr
idFunctionType>
225template<
class RangeDiscreteFunction,
class Model,
class Penalty >
238 for(
IteratorType it = dfSpace.begin(); it != end; ++it )
243 bool needsCalculation = model().init( entity );
244 if (! needsCalculation )
251 const auto uLocal = u.localFunction( entity );
255 const int quadOrder = uLocal.order() + wLocal.order();
259 const size_t numQuadraturePoints = quadrature.nop();
260 for(
size_t pt = 0; pt < numQuadraturePoints; ++pt )
262 const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
263 const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );
266 uLocal.evaluate( quadrature[ pt ], vu );
268 uLocal.jacobian( quadrature[ pt ], du );
272 model().source( quadrature[ pt ], vu, du, avu );
278 model().flux( quadrature[ pt ], vu, du, adu );
282 wLocal.axpy( quadrature[ pt ], avu, adu );
285 if ( ! dfSpace.continuous() )
287 const double area = entity.geometry().volume();
293 if ( intersection.neighbor() )
296 const EntityType outside = intersection.outside() ;
302 const double intersectionArea = intersectionGeometry.volume();
303 const double beta = penalty()(intersection,intersectionArea,area,outside.geometry().volume());
305 auto uOutLocal = u.localFunction( outside );
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();
312 for(
size_t pt = 0; pt < numQuadraturePoints; ++pt )
316 const typename FaceQuadratureType::LocalCoordinateType &x = quadInside.localPoint( pt );
317 DomainType normal = intersection.integrationOuterNormal( x );
318 double faceVol = normal.two_norm();
320 const double weight = quadInside.weight( pt ) * faceVol;
324 uLocal.evaluate( quadInside[ pt ], vuIn );
325 uLocal.jacobian( quadInside[ pt ], duIn );
326 uOutLocal.evaluate( quadOutside[ pt ], vuOut );
327 uOutLocal.jacobian( quadOutside[ pt ], duOut );
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];
335 model().init( outside );
336 model().flux( quadOutside[ pt ], vuOut, duOut, aduOut );
337 auto pFactor = model().penalty( quadOutside[ pt ], vuOut );
338 model().init( entity );
339 model().flux( quadInside[ pt ], vuIn, duIn, aduIn );
342 pFactor += model().penalty( quadInside[ pt ], vuIn );
350 for (
unsigned int r=0;r<dimRange;++r)
351 value[r] *= beta * pFactor[r]/2.;
355 aduIn.umv(normal,value);
359 model().flux( quadInside[ pt ], vuIn, dvalue, advalue );
364 wLocal.axpy( quadInside[ pt ], value, advalue );
368 else if( intersection.boundary() )
370 std::array<int,dimRange> components;
372 model().isDirichletIntersection( intersection, components);
378 const double intersectionArea = intersectionGeometry.volume();
379 const double beta = penalty()(intersection,intersectionArea,area,area);
381 FaceQuadratureType quadInside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
382 const size_t numQuadraturePoints = quadInside.nop();
384 for(
size_t pt = 0; pt < numQuadraturePoints; ++pt )
386 const typename FaceQuadratureType::LocalCoordinateType &x = quadInside.localPoint( pt );
387 const DomainType normal = intersection.integrationOuterNormal( x );
388 const double weight = quadInside.weight( pt );
391 model().dirichlet(1, quadInside[pt], bndValue);
398 uLocal.evaluate( quadInside[ pt ], vuIn );
399 uLocal.jacobian( quadInside[ pt ], duIn );
404 auto pFactor = model().penalty( quadInside[ pt ], vuIn );
406 for (
unsigned int r=0;r<dimRange;++r)
407 value[r] *= beta * pFactor[r] * intersectionGeometry.integrationElement( x );
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];
415 model().flux( quadInside[ pt ], jump, dvalue, advalue );
419 model().flux( quadInside[ pt ], vuIn, duIn, aduIn );
420 aduIn.mmv(normal,value);
422 for (
int r=0;r<dimRange;++r)
431 wLocal.axpy( quadInside[ pt ], value, advalue );
445template<
class JacobianOperator,
class Model,
class Penalty >
448 ::apply (
const GF &u, JacobianOperator &jOp )
const
452 typedef typename JacobianOperator::LocalMatrixType LocalMatrixType;
453 typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType BasisFunctionSetType;
459 jOp.reserve(stencil);
464 const unsigned int numDofs = rangeSpace.blockMapper().maxNumDofs() *
465 DiscreteFunctionSpaceType :: localBlockSize ;
467 std::vector< RangeType > phi( numDofs );
468 std::vector< JacobianRangeType > dphi( numDofs );
470 std::vector< RangeType > phiNb( numDofs );
471 std::vector< JacobianRangeType > dphiNb( numDofs );
475 for(
IteratorType it = rangeSpace.begin(); it != end; ++it )
479 bool needsCalculation = model().init( entity );
480 if (! needsCalculation )
485 const auto uLocal = u.localFunction( entity );
487 LocalMatrixType jLocal = jOp.localMatrix( entity, entity );
489 const BasisFunctionSetType &baseSet = jLocal.domainBasisFunctionSet();
490 const unsigned int numBaseFunctions = baseSet.size();
493 const size_t numQuadraturePoints = quadrature.nop();
494 for(
size_t pt = 0; pt < numQuadraturePoints; ++pt )
496 const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
497 const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );
500 baseSet.evaluateAll( quadrature[ pt ], phi );
503 baseSet.jacobianAll( quadrature[ pt ], dphi );
508 uLocal.evaluate( quadrature[ pt ], u0 );
509 uLocal.jacobian( quadrature[ pt ], jacU0 );
513 for(
unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol )
516 model().linSource( u0, jacU0, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], aphi );
519 model().linFlux( u0, jacU0, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], adphi );
522 jLocal.column( localCol ).axpy( phi, dphi, aphi, adphi, weight );
525 if ( rangeSpace.continuous() )
528 double area = geometry.volume();
531 iit != endiit ; ++ iit )
535 if( intersection.neighbor() )
544 LocalMatrixType localOpNb = jOp.localMatrix( neighbor, entity );
546 const BasisFunctionSetType &baseSetNb = localOpNb.domainBasisFunctionSet();
550 const double intersectionArea = intersectionGeometry.volume();
551 const double beta = penalty()(intersection,intersectionArea,area,neighbor.geometry().volume());
555 FaceQuadratureType::INSIDE);
557 FaceQuadratureType::OUTSIDE);
559 const size_t numFaceQuadPoints = faceQuadInside.nop();
560 for(
size_t pt = 0; pt < numFaceQuadPoints; ++pt )
562 const typename FaceQuadratureType::LocalCoordinateType &x = faceQuadInside.localPoint( pt );
563 DomainType normal = intersection.integrationOuterNormal( x );
564 double faceVol = normal.two_norm();
567 const double quadWeight = faceQuadInside.weight( pt );
568 const double weight = quadWeight * faceVol;
573 uLocal.evaluate( faceQuadInside[ pt ], u0En );
574 uLocal.jacobian( faceQuadInside[ pt ], u0EnJac );
581 baseSet.evaluateAll( faceQuadInside[ pt ], phi );
584 baseSet.jacobianAll( faceQuadInside[ pt ], dphi );
591 baseSetNb.evaluateAll( faceQuadOutside[ pt ], phiNb );
594 baseSetNb.jacobianAll( faceQuadOutside[ pt ], dphiNb );
596 model().init( entity );
597 for(
unsigned int i = 0; i < numBaseFunctions; ++i )
600 model().linFlux( u0En, u0EnJac, faceQuadInside[ pt ], phi[i], adphiEn, dphi[ i ] );
602 auto pFactor = model().penalty( faceQuadInside[ pt ], u0En );
603 model().init( neighbor );
604 for(
unsigned int i = 0; i < numBaseFunctions; ++i )
607 model().linFlux( u0En, u0EnJac, faceQuadOutside[ pt ], phiNb[i], adphiNb, dphiNb[ i ] );
609 pFactor += model().penalty( faceQuadOutside[ pt ], u0En );
610 model().init( entity );
614 for(
unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol )
620 dphi[localCol].usmv( -0.5, normal, valueEn );
623 dphiNb[localCol].usmv( -0.5, normal, valueNb );
626 for (
unsigned int r=0;r<dimRange;++r)
628 valueEn[r] += beta*pFactor[r]/2.*phi[localCol][r];
629 valueNb[r] -= beta*pFactor[r]/2.*phiNb[localCol][r];
632 for (
int r=0; r< dimRange; ++r )
633 for (
int d=0; d< dimDomain; ++d )
636 dvalueEn[r][d] = - 0.5 * normal[d] * phi[localCol][r];
639 dvalueNb[r][d] = 0.5 * normal[d] * phiNb[localCol][r];
642 jLocal.column( localCol ).axpy( phi, dphi, valueEn, dvalueEn, weight );
643 localOpNb.column( localCol ).axpy( phi, dphi, valueNb, dvalueNb, weight );
648 else if( intersection.boundary() )
650 std::array<int,dimRange> components;
652 model().isDirichletIntersection( intersection, components);
658 const double intersectionArea = intersectionGeometry.volume();
659 const double beta = penalty()(intersection,intersectionArea,area,area);
663 FaceQuadratureType::INSIDE);
665 const size_t numFaceQuadPoints = faceQuadInside.nop();
666 for(
size_t pt = 0; pt < numFaceQuadPoints; ++pt )
668 const typename FaceQuadratureType::LocalCoordinateType &x = faceQuadInside.localPoint( pt );
669 DomainType normal = intersection.integrationOuterNormal( x );
670 double faceVol = normal.two_norm();
673 const double quadWeight = faceQuadInside.weight( pt );
674 const double weight = quadWeight * faceVol;
678 uLocal.evaluate( faceQuadInside[ pt ], u0En );
679 uLocal.jacobian( faceQuadInside[ pt ], u0EnJac );
686 baseSet.evaluateAll( faceQuadInside[ pt ], phi );
689 baseSet.jacobianAll( faceQuadInside[ pt ], dphi );
691 for(
unsigned int i = 0; i < numBaseFunctions; ++i )
694 model().linFlux( u0En, u0EnJac, faceQuadInside[ pt ], phi[i], adphiEn, dphi[ i ] );
697 auto pFactor = model().penalty( faceQuadInside[ pt ], u0En );
699 for(
unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol )
705 dphi[localCol].usmv( -1.0, normal, valueEn );
708 for (
unsigned int r=0;r<dimRange;++r)
709 valueEn[r] += beta*pFactor[r]*phi[ localCol ][r];
712 for (
int r=0; r< dimRange; ++r )
713 for (
int d=0; d< dimDomain; ++d )
716 dvalueEn[r][d] = - 0.5 * normal[d] * phi[localCol][r];
719 for (
int r=0;r<dimRange;++r)
726 jLocal.column( localCol ).axpy( phi, dphi, valueEn, dvalueEn, weight );
double min(const Dune::Fem::Double &v, const double p)
Definition: double.hh:953
Dune::Entity< codim, dim, Grid, Implementation > make_entity(Dune::Entity< codim, dim, Grid, Implementation > entity)
Definition: compatibility.hh:21
static ParameterContainer & container()
Definition: io/parameter.hh:193
abstract differentiable operator
Definition: differentiableoperator.hh:29
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:385
quadrature on the codim-0 reference element
Definition: elementquadrature.hh:58
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
DGEllipticOperator(const DiscreteFunctionSpaceType &space, const ModelType &model, const Dune::Fem::ParameterReader ¶meter=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 ¶meter=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
const ModelType & model() const
Definition: dgelliptic.hh:137
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 ¶meter=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
const ModelType & model() const
Definition: dgelliptic.hh:137
DifferentiableDGEllipticOperator(const DiscreteFunctionSpaceType &space, const ModelType &model, const Dune::Fem::ParameterReader ¶meter=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