1#ifndef DUNE_FEM_SCHEMES_MASSLUMPING_HH
2#define DUNE_FEM_SCHEMES_MASSLUMPING_HH
5#include <dune/fem/schemes/galerkin.hh>
6#include <dune/fem/function/localfunction/temporary.hh>
7#include <dune/fem/common/bindguard.hh>
8#include <dune/fem/quadrature/lumpingquadrature.hh>
18 template <
typename DomainFunction,
typename RangeFunction >
19 class PetscLinearOperator ;
25 template<
class Integrands,
class MassIntegrands,
class DomainFunction,
class RangeFunction = DomainFunction >
26 struct MassLumpingOperator
27 :
public virtual Operator< DomainFunction, RangeFunction >
29 typedef DomainFunction DomainFunctionType;
30 typedef RangeFunction RangeFunctionType;
33 static_assert( std::is_same< typename DomainFunctionType::GridPartType, typename RangeFunctionType::GridPartType >::value,
"DomainFunction and RangeFunction must be defined on the same grid part." );
35 typedef typename RangeFunctionType::GridPartType GridPartType;
38 typedef Impl::GalerkinOperator< GridPartType > GalerkinOperatorImplType;
39 typedef Impl::LocalGalerkinOperator< Integrands > LocalGalerkinOperatorImplType;
41 template <
class Space>
42 struct LumpingQuadratureSelector
44 typedef typename Space :: GridPartType GridPartType;
45 typedef CachingLumpingQuadrature< GridPartType, 0 > InteriorQuadratureType;
48 typedef CachingQuadrature< GridPartType, 1, Capabilities::DefaultQuadrature< Space > :: template DefaultQuadratureTraits > SurfaceQuadratureType;
52 typedef Impl::LocalGalerkinOperator< MassIntegrands, LumpingQuadratureSelector > MassOperatorImplType;
54 typedef typename RangeFunctionType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
56 template<
class... Args >
57 explicit MassLumpingOperator (
const GridPartType &gridPart,
58 const Integrands& integrands,
59 const MassIntegrands& massIntegrands,
61 : iterators_( gridPart ),
63 localOp_( gridPart,
std::move( integrands ),
std::forward< Args >( args )... ),
64 mass_( gridPart,
std::move( massIntegrands ),
std::forward< Args >( args )... ),
65 gridSizeInterior_( 0 ),
68 if( mass().hasSkeleton() )
69 DUNE_THROW(InvalidStateException,
"MassLumpingOperator: Mass model cannot have a skeleton term");
72 void setCommunicate(
const bool communicate )
74 communicate_ = communicate;
77 std::cout <<
"MassLumpingOperator::setCommunicate: communicate was disabled!" << std::endl;
81 void setQuadratureOrders(
unsigned int interior,
unsigned int surface)
84 size_t size = localOp_.size();
85 for(
size_t i=0; i<
size; ++i )
86 localOp_[ i ].setQuadratureOrders(interior,surface);
89 virtual void operator() (
const DomainFunctionType &u, RangeFunctionType &w )
const final override
94 template<
class Gr
idFunction >
95 void operator() (
const GridFunction &u, RangeFunctionType &w )
const
100 const GridPartType &gridPart ()
const {
return op().gridPart(); }
102 typedef Integrands ModelType;
103 typedef Integrands DirichletModelType;
104 ModelType &model()
const {
return localOperator().model(); }
106 [[deprecated(
"Use localOperator instead!")]]
107 const LocalGalerkinOperatorImplType& impl()
const {
return localOperator(); }
110 const LocalGalerkinOperatorImplType& localOperator()
const {
return (*localOp_); }
111 const MassOperatorImplType& mass()
const {
return (*mass_); }
113 std::size_t gridSizeInterior ()
const {
return gridSizeInterior_; }
116 const GalerkinOperatorImplType& op()
const {
return (*opImpl_); }
118 template<
class Gr
idFunction >
119 void evaluate(
const GridFunction &u, RangeFunctionType &w )
const
124 std::shared_mutex mutex;
126 auto doEval = [
this, &u, &w, &mutex] ()
129 std::tuple<
const LocalGalerkinOperatorImplType&,
130 const MassOperatorImplType& > integrands( localOperator(), mass() );
134 this->op().evaluate( u, w, this->iterators_, integrands, mutex );
137 bool singleThreadModeError = false ;
141 MPIManager :: run ( doEval );
144 gridSizeInterior_ = Impl::accumulateGridSize( opImpl_ );
146 catch (
const SingleThreadModeError& e )
148 singleThreadModeError =
true;
152 if( singleThreadModeError )
158 std::tuple<
const LocalGalerkinOperatorImplType&,
159 const MassOperatorImplType& > integrands( localOperator(), mass() );
161 op().evaluate( u, w, iterators_, integrands );
164 gridSizeInterior_ = op().gridSizeInterior();
174 mutable ThreadIteratorType iterators_;
177 ThreadSafeValue< MassOperatorImplType > mass_;
179 mutable std::size_t gridSizeInterior_;
188 template<
class Integrands,
class MassIntegrands,
class JacobianOperator >
189 class MassLumpingDifferentiableOperator
190 :
public MassLumpingOperator< Integrands, MassIntegrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType >,
191 public DifferentiableOperator< JacobianOperator >
193 typedef MassLumpingOperator< Integrands, MassIntegrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType > BaseType;
196 typedef typename BaseType :: LocalGalerkinOperatorImplType LocalGalerkinOperatorImplType;
197 typedef typename BaseType :: MassOperatorImplType MassOperatorImplType;
199 typedef JacobianOperator JacobianOperatorType;
201 typedef typename BaseType::DomainFunctionType DomainFunctionType;
202 typedef typename BaseType::RangeFunctionType RangeFunctionType;
203 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType;
204 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType;
206 typedef DiagonalAndNeighborStencil< DomainDiscreteFunctionSpaceType, RangeDiscreteFunctionSpaceType > DiagonalAndNeighborStencilType;
207 typedef DiagonalStencil< DomainDiscreteFunctionSpaceType, RangeDiscreteFunctionSpaceType > DiagonalStencilType;
209 typedef typename BaseType::GridPartType GridPartType;
211 template<
class... Args >
212 explicit MassLumpingDifferentiableOperator (
const DomainDiscreteFunctionSpaceType &dSpace,
213 const RangeDiscreteFunctionSpaceType &rSpace,
215 : BaseType( rSpace.gridPart(),
std::forward< Args >( args )... ),
218 domainSpaceSequence_(dSpace.sequence()),
219 rangeSpaceSequence_(rSpace.sequence()),
220 stencilDAN_(), stencilD_()
223 stencilDAN_.reset(
new DiagonalAndNeighborStencilType( dSpace_, rSpace_ ) );
225 stencilD_.reset(
new DiagonalStencilType( dSpace_, rSpace_ ) );
228 virtual void jacobian (
const DomainFunctionType &u, JacobianOperatorType &jOp )
const final override
234 template<
class Gr
idFunction >
235 void jacobian (
const GridFunction &u, JacobianOperatorType &jOp )
const
240 const DomainDiscreteFunctionSpaceType& domainSpace()
const
244 const RangeDiscreteFunctionSpaceType& rangeSpace()
const
249 const RangeDiscreteFunctionSpaceType& space()
const
254 using BaseType::gridPart;
255 using BaseType::localOperator;
258 bool hasSkeleton()
const
260 std::tuple< const LocalGalerkinOperatorImplType&, const MassOperatorImplType& > integrands( localOperator(), mass() );
261 return op().hasSkeleton( integrands );
265 using BaseType::mass;
266 using BaseType::iterators_;
267 using BaseType::gridSizeInterior_;
269 void prepare( JacobianOperatorType& jOp )
const
271 if ( domainSpaceSequence_ != domainSpace().sequence()
272 || rangeSpaceSequence_ != rangeSpace().sequence() )
274 domainSpaceSequence_ = domainSpace().sequence();
275 rangeSpaceSequence_ = rangeSpace().sequence();
278 assert( stencilDAN_ );
279 stencilDAN_->update();
289 jOp.reserve( *stencilDAN_ );
291 jOp.reserve( *stencilD_ );
296 template <
class Gr
idFunction >
297 void assemble(
const GridFunction &u, JacobianOperatorType &jOp )
const
302 bool singleThreadModeError =
false;
303 std::shared_mutex mutex;
305 auto doAssemble = [
this, &u, &jOp, &mutex] ()
308 std::tuple< const LocalGalerkinOperatorImplType&, const MassOperatorImplType& > integrands( localOperator(), mass() );
309 this->op().assemble( u, jOp, this->iterators_, integrands, mutex );
314 MPIManager :: run ( doAssemble );
317 gridSizeInterior_ = Impl::accumulateGridSize( this->opImpl_ );
319 catch (
const SingleThreadModeError& e )
321 singleThreadModeError =
true;
324 if (singleThreadModeError)
329 std::tuple< const LocalGalerkinOperatorImplType&, const MassOperatorImplType& > integrands( localOperator(), mass() );
330 op().assemble( u, jOp, iterators_, integrands );
333 gridSizeInterior_ = op().gridSizeInterior();
341 const DomainDiscreteFunctionSpaceType &dSpace_;
342 const RangeDiscreteFunctionSpaceType &rSpace_;
343 mutable int domainSpaceSequence_, rangeSpaceSequence_;
345 mutable std::unique_ptr< DiagonalAndNeighborStencilType > stencilDAN_;
346 mutable std::unique_ptr< DiagonalStencilType > stencilD_;
354 template<
class Integrands,
class DomainFunction,
class RangeFunction >
355 class MassLumpingAutomaticDifferenceGalerkinOperator
356 :
public MassLumpingOperator< Integrands, DomainFunction, RangeFunction >,
357 public AutomaticDifferenceOperator< DomainFunction, RangeFunction >
359 typedef MassLumpingOperator< Integrands, DomainFunction, RangeFunction > BaseType;
363 typedef typename BaseType::GridPartType GridPartType;
365 template<
class... Args >
366 explicit MassLumpingAutomaticDifferenceGalerkinOperator (
const GridPartType &gridPart, Args &&... args )
367 : BaseType( gridPart,
std::forward< Args >( args )... ), AutomaticDifferenceOperatorType()
378 template<
class Integrands,
class MassIntegrands,
379 class LinearOperator,
bool addDirichletBC,
380 template <
class,
class,
class>
class DifferentiableGalerkinOperatorImpl >
381 struct MassLumpingSchemeTraits
383 template <
class O,
bool addDBC>
384 struct DirichletBlockSelector {
using type = void; };
386 struct DirichletBlockSelector<O,true> {
using type =
typename O::DirichletBlockVector; };
388 using DifferentiableOperatorType = std::conditional_t< addDirichletBC,
389 DirichletWrapperOperator< DifferentiableGalerkinOperatorImpl< Integrands, MassIntegrands, LinearOperator >>,
390 DifferentiableGalerkinOperatorImpl< Integrands, MassIntegrands, LinearOperator > >;
392 using DirichletBlockVector =
typename DirichletBlockSelector<
393 DirichletWrapperOperator<
394 DifferentiableGalerkinOperatorImpl< Integrands, MassIntegrands, LinearOperator >>,
395 addDirichletBC>::type;
397 typedef DifferentiableOperatorType type;
403 template<
class Integrands,
class MassIntegrands,
404 class LinearOperator,
class LinearInverseOperator,
bool addDirichletBC,
405 template <
class,
class,
class>
class DifferentiableGalerkinOperatorImpl = MassLumpingDifferentiableOperator >
406 struct MassLumpingSchemeImpl :
public FemScheme< typename MassLumpingSchemeTraits< Integrands, MassIntegrands, LinearOperator,
407 addDirichletBC, DifferentiableGalerkinOperatorImpl>::type,
408 LinearInverseOperator >
410 typedef FemScheme<
typename MassLumpingSchemeTraits< Integrands, MassIntegrands, LinearOperator,
411 addDirichletBC, DifferentiableGalerkinOperatorImpl>::type,
412 LinearInverseOperator >
416 typedef typename BaseType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
419 typedef MassIntegrands MassModelType;
421 MassLumpingSchemeImpl (
const DiscreteFunctionSpaceType &dfSpace,
422 const Integrands &integrands,
423 const MassIntegrands& massIntegrands,
424 const ParameterReader& parameter = Parameter::container() )
427 std::move(integrands),
428 std::move(massIntegrands) )
437 template<
class Integrands,
class MassIntegrands,
class LinearOperator,
class InverseOperator,
bool addDirichletBC >
438 using MassLumpingScheme = Impl::MassLumpingSchemeImpl< Integrands, MassIntegrands, LinearOperator, InverseOperator, addDirichletBC,
439 MassLumpingDifferentiableOperator >;
static bool verbose()
obtain the cached value for fem.verbose with default verbosity level 2
Definition: parameter.hh:466
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75