40#ifndef DUNE_FEM_SCHEMES_ELLIPTIC_HH
41#define DUNE_FEM_SCHEMES_ELLIPTIC_HH
45#include <dune/common/timer.hh>
46#include <dune/common/fmatrix.hh>
64#include <dune/fempy/quadrature/fempyquadratures.hh>
70template<
class DomainDiscreteFunction,
class RangeDiscreteFunction,
class Model>
90 typedef typename RangeDiscreteFunctionSpaceType::IteratorType
IteratorType;
93 typedef typename RangeDiscreteFunctionSpaceType::DomainType
DomainType;
94 typedef typename RangeDiscreteFunctionSpaceType::GridPartType
GridPartType;
112 dSpace_(dSpace), rSpace_(rSpace),
156template<
class JacobianOperator,
class Model >
158:
public EllipticOperator< typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType, Model >,
180 typedef typename RangeDiscreteFunctionSpaceType::IteratorType
IteratorType;
183 typedef typename RangeDiscreteFunctionSpaceType::DomainType
DomainType;
184 typedef typename RangeDiscreteFunctionSpaceType::GridPartType
GridPartType;
209 template <
class Gr
idFunctionType>
212 template <
class Gr
idFunctionType>
214 using BaseType::operator();
226template<
class DomainDiscreteFunction,
class RangeDiscreteFunction,
class Model >
242 if( !model().init( entity ) )
252 const int quadOrder = interiorOrder_==-1?
253 uLocal.order() + wLocal.order() : interiorOrder_;
257 const size_t numQuadraturePoints = quadrature.nop();
258 for(
size_t pt = 0; pt < numQuadraturePoints; ++pt )
261 const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
262 const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );
265 uLocal.evaluate( quadrature[ pt ], vu );
267 uLocal.jacobian( quadrature[ pt ], du );
271 model().source( quadrature[ pt ], vu, du, avu );
277 model().flux( quadrature[ pt ], vu, du, adu );
281 wLocal.axpy( quadrature[ pt ], avu, adu );
286 if( model().hasNeumanBoundary() && entity.hasBoundaryIntersections() )
292 if( !intersection.boundary() )
295 std::array<int,RangeRangeType::dimension> components(0);
297 const bool hasDirichletComponent = model().isDirichletIntersection( intersection, components );
299 const auto &intersectionGeometry = intersection.geometry();
300 const int quadOrder = surfaceOrder_==-1?
301 uLocal.order() + wLocal.order() : surfaceOrder_;
302 FaceQuadratureType quadInside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
303 const std::size_t numQuadraturePoints = quadInside.nop();
304 for( std::size_t pt = 0; pt < numQuadraturePoints; ++pt )
306 const auto &x = quadInside.localPoint( pt );
307 double weight = quadInside.weight( pt ) * intersectionGeometry.integrationElement( x );
309 uLocal.evaluate( quadInside[ pt ], vu );
311 model().alpha( quadInside[ pt ], vu, alpha );
313 for(
int k = 0; k < RangeRangeType::dimension; ++k )
314 if( hasDirichletComponent && components[ k ] )
316 wLocal.axpy( quadInside[ pt ], alpha );
330template<
class JacobianOperator,
class Model >
337 typedef typename JacobianOperator::LocalMatrixType LocalMatrixType;
338 typedef typename DomainDiscreteFunctionSpaceType::BasisFunctionSetType DomainBasisFunctionSetType;
339 typedef typename RangeDiscreteFunctionSpaceType::BasisFunctionSetType RangeBasisFunctionSetType;
345 jOp.reserve(stencil);
348 const int domainBlockSize = domainSpace.localBlockSize;
349 std::vector< typename DomainLocalFunctionType::RangeType > phi( domainSpace.blockMapper().maxNumDofs()*domainBlockSize );
350 std::vector< typename DomainLocalFunctionType::JacobianRangeType > dphi( domainSpace.blockMapper().maxNumDofs()*domainBlockSize );
351 const int rangeBlockSize = rangeSpace.localBlockSize;
352 std::vector< typename RangeLocalFunctionType::RangeType > rphi( rangeSpace.blockMapper().maxNumDofs()*rangeBlockSize );
353 std::vector< typename RangeLocalFunctionType::JacobianRangeType > rdphi( rangeSpace.blockMapper().maxNumDofs()*rangeBlockSize );
359 if( !model().init( entity ) )
365 LocalMatrixType jLocal = jOp.localMatrix( entity, entity );
367 const DomainBasisFunctionSetType &domainBaseSet = jLocal.domainBasisFunctionSet();
368 const RangeBasisFunctionSetType &rangeBaseSet = jLocal.rangeBasisFunctionSet();
369 const std::size_t domainNumBasisFunctions = domainBaseSet.size();
371 const int quadOrder = interiorOrder_==-1?
372 domainSpace.order() + rangeSpace.order() : interiorOrder_;
374 const size_t numQuadraturePoints = quadrature.nop();
375 for(
size_t pt = 0; pt < numQuadraturePoints; ++pt )
378 const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
379 const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );
382 domainBaseSet.evaluateAll( quadrature[ pt ], phi );
383 rangeBaseSet.evaluateAll( quadrature[ pt ], rphi );
386 domainBaseSet.jacobianAll( quadrature[ pt ], dphi );
387 rangeBaseSet.jacobianAll( quadrature[ pt ], rdphi );
392 uLocal.evaluate( quadrature[ pt ], u0 );
393 uLocal.jacobian( quadrature[ pt ], jacU0 );
397 for( std::size_t localCol = 0; localCol < domainNumBasisFunctions; ++localCol )
400 model().linSource( u0, jacU0, quadrature[ pt ], phi[ localCol ], dphi[localCol], aphi );
403 model().linFlux( u0, jacU0, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], adphi );
406 jLocal.column( localCol ).axpy( rphi, rdphi, aphi, adphi, weight );
411 if( model().hasNeumanBoundary() && entity.hasBoundaryIntersections() )
417 if( !intersection.boundary() )
420 std::array<int,RangeRangeType::dimension> components(0);
421 bool hasDirichletComponent = model().isDirichletIntersection( intersection, components );
423 const auto &intersectionGeometry = intersection.geometry();
424 const int quadOrder = surfaceOrder_==-1?
425 domainSpace.order() + rangeSpace.order() : surfaceOrder_;
426 FaceQuadratureType quadInside( rangeSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
427 const std::size_t numQuadraturePoints = quadInside.nop();
428 for( std::size_t pt = 0; pt < numQuadraturePoints; ++pt )
430 const auto &x = quadInside.localPoint( pt );
431 const double weight = quadInside.weight( pt ) * intersectionGeometry.integrationElement( x );
433 uLocal.evaluate( quadInside[ pt ], u0 );
434 domainBaseSet.evaluateAll( quadInside[ pt ], phi );
435 for( std::size_t localCol = 0; localCol < domainNumBasisFunctions; ++localCol )
438 model().linAlpha( u0,quadInside[ pt ], phi[ localCol ], alpha );
439 for(
int k = 0; k < RangeRangeType::dimension; ++k )
440 if( hasDirichletComponent && components[ k ] )
442 jLocal.column( localCol ).axpy( phi, alpha, weight );
typename Impl::ConstLocalFunction< GridFunction >::Type ConstLocalFunction
Definition: const.hh:604
static auto bindGuard(Object &object, Args &&... args) -> std::enable_if_t< isBindable< Object, Args... >::value, BindGuard< Object > >
Definition: bindguard.hh:67
Definition: common/localcontribution.hh:14
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) for all entities in the space. Defailt for an operator over a L...
Definition: stencil.hh:348
quadrature class supporting base function caching
Definition: cachingquadrature.hh:41
[Class for elliptic operator]
Definition: elliptic.hh:74
DomainLocalFunctionType::JacobianRangeType DomainJacobianRangeType
Definition: elliptic.hh:83
RangeLocalFunctionType::JacobianRangeType RangeJacobianRangeType
Definition: elliptic.hh:87
Model ModelType
Definition: elliptic.hh:77
DomainFunctionType::LocalFunctionType DomainLocalFunctionType
Definition: elliptic.hh:81
int interiorOrder_
Definition: elliptic.hh:143
ModelType & model()
Definition: elliptic.hh:134
IteratorType::Entity EntityType
Definition: elliptic.hh:91
RangeDiscreteFunctionSpaceType::GridPartType GridPartType
Definition: elliptic.hh:94
RangeLocalFunctionType::RangeType RangeRangeType
Definition: elliptic.hh:86
RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType
Definition: elliptic.hh:84
const DomainDiscreteFunctionSpaceType & domainSpace() const
Definition: elliptic.hh:125
RangeDiscreteFunctionSpaceType::IteratorType IteratorType
Definition: elliptic.hh:90
Dune::Fem::CachingQuadrature< GridPartType, 0, Dune::FemPy::FempyQuadratureTraits > QuadratureType
Definition: elliptic.hh:98
EllipticOperator(const RangeDiscreteFunctionSpaceType &rangeSpace, ModelType &model, const Dune::Fem::ParameterReader ¶meter=Dune::Fem::Parameter::container())
Definition: elliptic.hh:101
void apply(const GF &u, RangeFunctionType &w) const
Definition: elliptic.hh:229
GridPartType::IntersectionIteratorType IntersectionIteratorType
Definition: elliptic.hh:95
void setQuadratureOrders(unsigned int interior, unsigned int surface)
Definition: elliptic.hh:136
Model DirichletModelType
Definition: elliptic.hh:78
void operator()(const GF &u, RangeFunctionType &w) const
Definition: elliptic.hh:120
DomainDiscreteFunction DomainFunctionType
Definition: elliptic.hh:75
RangeFunctionType::LocalFunctionType RangeLocalFunctionType
Definition: elliptic.hh:85
virtual void operator()(const DomainFunctionType &u, RangeFunctionType &w) const
application operator
Definition: elliptic.hh:117
RangeDiscreteFunctionSpaceType::DomainType DomainType
Definition: elliptic.hh:93
DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: elliptic.hh:80
Dune::Fem::CachingQuadrature< GridPartType, 1, Dune::FemPy::FempyQuadratureTraits > FaceQuadratureType
Definition: elliptic.hh:99
IntersectionIteratorType::Intersection IntersectionType
Definition: elliptic.hh:96
EntityType::Geometry GeometryType
Definition: elliptic.hh:92
const ModelType & model() const
Definition: elliptic.hh:133
int surfaceOrder_
Definition: elliptic.hh:143
const RangeDiscreteFunctionSpaceType & rangeSpace() const
Definition: elliptic.hh:129
RangeDiscreteFunction RangeFunctionType
Definition: elliptic.hh:76
DomainLocalFunctionType::RangeType DomainRangeType
Definition: elliptic.hh:82
EllipticOperator(const DomainDiscreteFunctionSpaceType &dSpace, const RangeDiscreteFunctionSpaceType &rSpace, ModelType &model, const Dune::Fem::ParameterReader ¶meter=Dune::Fem::Parameter::container())
Definition: elliptic.hh:107
[Class for linearizable elliptic operator]
Definition: elliptic.hh:161
DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: elliptic.hh:170
DifferentiableEllipticOperator(const DomainDiscreteFunctionSpaceType &dSpace, const RangeDiscreteFunctionSpaceType &rSpace, ModelType &model, const Dune::Fem::ParameterReader ¶meter=Dune::Fem::Parameter::container())
contructor
Definition: elliptic.hh:199
EllipticOperator< typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType, Model > BaseType
Definition: elliptic.hh:162
RangeDiscreteFunctionSpaceType::IteratorType IteratorType
Definition: elliptic.hh:180
DifferentiableEllipticOperator(const RangeDiscreteFunctionSpaceType &space, ModelType &model, const Dune::Fem::ParameterReader ¶meter=Dune::Fem::Parameter::container())
contructor
Definition: elliptic.hh:193
BaseType::DomainFunctionType DomainFunctionType
Definition: elliptic.hh:166
void assemble(const GridFunctionType &u, JacobianOperatorType &jOp) const
IntersectionIteratorType::Intersection IntersectionType
Definition: elliptic.hh:186
DomainLocalFunctionType::JacobianRangeType DomainJacobianRangeType
Definition: elliptic.hh:173
BaseType::RangeFunctionType RangeFunctionType
Definition: elliptic.hh:167
RangeLocalFunctionType::JacobianRangeType RangeJacobianRangeType
Definition: elliptic.hh:177
void jacobian(const GridFunctionType &u, JacobianOperatorType &jOp) const
Definition: elliptic.hh:210
RangeDiscreteFunctionSpaceType::GridPartType GridPartType
Definition: elliptic.hh:184
RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType
Definition: elliptic.hh:174
RangeDiscreteFunctionSpaceType::DomainType DomainType
Definition: elliptic.hh:183
BaseType::QuadratureType QuadratureType
Definition: elliptic.hh:188
BaseType::FaceQuadratureType FaceQuadratureType
Definition: elliptic.hh:190
DomainFunctionType::LocalFunctionType DomainLocalFunctionType
Definition: elliptic.hh:171
RangeFunctionType::LocalFunctionType RangeLocalFunctionType
Definition: elliptic.hh:175
JacobianOperator JacobianOperatorType
Definition: elliptic.hh:164
const ModelType & model() const
Definition: elliptic.hh:133
DomainLocalFunctionType::RangeType DomainRangeType
Definition: elliptic.hh:172
GridPartType::IntersectionIteratorType IntersectionIteratorType
Definition: elliptic.hh:185
EntityType::Geometry GeometryType
Definition: elliptic.hh:182
BaseType::ModelType ModelType
Definition: elliptic.hh:168
IteratorType::Entity EntityType
Definition: elliptic.hh:181
void jacobian(const DomainFunctionType &u, JacobianOperatorType &jOp) const
method to setup the jacobian of the operator for storage in a matrix
Definition: elliptic.hh:207
RangeLocalFunctionType::RangeType RangeRangeType
Definition: elliptic.hh:176