1#ifndef DUNE_FEM_ConservationLaw_MODEL_HH
2#define DUNE_FEM_ConservationLaw_MODEL_HH
9#include <dune/fem/solver/timeprovider.hh>
10#include <dune/fem/io/parameter.hh>
11#include <dune/fem/space/common/functionspace.hh>
12#include <dune/fem/function/common/gridfunctionadapter.hh>
13#include <dune/fem/quadrature/cachingquadrature.hh>
14#include <dune/fempy/quadrature/fempyquadratures.hh>
16#define VirtualConservationLawModelMethods(POINT) \
17 virtual void source ( const POINT &x,\
18 const DRangeType &value,\
19 const DJacobianRangeType &gradient,\
20 RRangeType &flux ) const = 0;\
21 virtual void linSource ( const DRangeType& uBar,\
22 const DJacobianRangeType &gradientBar,\
24 const DRangeType &value,\
25 const DJacobianRangeType &gradient,\
26 RRangeType &flux ) const = 0;\
27 virtual void flux ( const POINT &x,\
28 const DRangeType &value,\
29 const DJacobianRangeType &gradient,\
30 RJacobianRangeType &result ) const = 0;\
32 virtual void diffusiveFlux ( const POINT &x,\
33 const DRangeType &value,\
34 const DJacobianRangeType &gradient,\
35 RJacobianRangeType &result ) const \
36 { flux(x, value, gradient, result ); } \
37 virtual void linFlux ( const DRangeType& uBar,\
38 const DJacobianRangeType& gradientBar,\
40 const DRangeType &value,\
41 const DJacobianRangeType &gradient,\
42 RJacobianRangeType &flux ) const = 0;\
44 virtual void linDiffusiveFlux ( const DRangeType& uBar,\
45 const DJacobianRangeType& gradientBar,\
47 const DRangeType &value,\
48 const DJacobianRangeType &gradient,\
49 RJacobianRangeType &flux ) const \
50 { linFlux(uBar, gradientBar, x, value, gradient, flux ); } \
51 virtual void fluxDivergence( const POINT &x,\
52 const DRangeType &value,\
53 const DJacobianRangeType &jacobian,\
54 const DHessianRangeType &hessian,\
55 RRangeType &flux) const = 0;\
56 virtual void alpha(const POINT &x,\
57 const DRangeType &value,\
58 RRangeType &val) const = 0;\
59 virtual void linAlpha(const DRangeType &uBar,\
61 const DRangeType &value,\
62 RRangeType &val) const = 0;\
63 virtual void dirichlet( int bndId, const POINT &x,\
64 RRangeType &value) const = 0;
66#define WrapperConservationLawModelMethods(POINT) \
67 virtual void source ( const POINT &x,\
68 const DRangeType &value,\
69 const DJacobianRangeType &gradient,\
70 RRangeType &flux ) const \
71 { impl().source(x, value, gradient, flux); } \
72 virtual void linSource ( const DRangeType& uBar,\
73 const DJacobianRangeType &gradientBar,\
75 const DRangeType &value,\
76 const DJacobianRangeType &gradient,\
77 RRangeType &flux ) const \
78 { impl().linSource(uBar, gradientBar, x, value, gradient, flux); } \
79 virtual void flux ( const POINT &x,\
80 const DRangeType &value,\
81 const DJacobianRangeType &gradient,\
82 RJacobianRangeType &flux ) const \
83 { impl().flux(x, value, gradient, flux); } \
84 virtual void linFlux ( const DRangeType& uBar,\
85 const DJacobianRangeType& gradientBar,\
87 const DRangeType &value,\
88 const DJacobianRangeType &gradient,\
89 RJacobianRangeType &flux ) const \
90 { impl().linFlux(uBar, gradientBar, x, value, gradient, flux); } \
91 virtual void fluxDivergence( const POINT &x,\
92 const DRangeType &value,\
93 const DJacobianRangeType &jacobian,\
94 const DHessianRangeType &hessian,\
95 RRangeType &flux) const \
96 { impl().fluxDivergence(x, value, jacobian, hessian, flux); } \
97 virtual void alpha(const POINT &x,\
98 const DRangeType &value,\
99 RRangeType &val) const \
100 { impl().alpha(x, value, val); } \
101 virtual void linAlpha(const DRangeType &uBar,\
103 const DRangeType &value,\
104 RRangeType &val) const \
105 { impl().linAlpha(uBar, x, value, val); } \
106 virtual void dirichlet( int bndId, const POINT &x,\
107 RRangeType &value) const \
108 { impl().dirichlet(bndId, x, value); }
110template<
class Gr
idPart,
int dimDomain,
int dimRange=dimDomain,
class RangeField =
double >
111struct ConservationLawModel
113 typedef GridPart GridPartType;
114 static const int dimD = dimDomain;
115 static const int dimR = dimRange;
116 typedef ConservationLawModel<GridPartType, dimD, dimR, RangeField> ModelType;
117 typedef RangeField RangeFieldType;
120 GridPart::dimensionworld, dimD > DFunctionSpaceType;
122 GridPart::dimensionworld, dimR > RFunctionSpaceType;
133 typedef typename GridPartType::template Codim<0>::EntityType EntityType;
135 typedef typename EntityType::Geometry::LocalCoordinate LocalDomainType;
138 template <
class F,
int d>
139 using Traits = Dune::FemPy::FempyQuadratureTraits<F,d>;
141 QuadraturePointWrapperType Point;
143 QuadraturePointWrapperType IntersectionPoint;
144 typedef typename Dune::Fem::ElementQuadrature< GridPartType, 0, Traits >::
145 QuadraturePointWrapperType ElementPoint;
146 typedef typename Dune::Fem::ElementQuadrature< GridPartType, 1, Traits >::
147 QuadraturePointWrapperType ElementIntersectionPoint;
151 QuadraturePointWrapperType OriginalPoint;
153 QuadraturePointWrapperType OriginalIntersectionPoint;
154 typedef typename Dune::Fem::ElementQuadrature< GridPartType, 0 >::
155 QuadraturePointWrapperType OriginalElementPoint;
156 typedef typename Dune::Fem::ElementQuadrature< GridPartType, 1 >::
157 QuadraturePointWrapperType OriginalElementIntersectionPoint;
165 ConservationLawModel( )
167 virtual ~ConservationLawModel() {}
169 virtual std::string name()
const = 0;
171 virtual bool init(
const EntityType &entity)
const = 0;
172 virtual void unbind( )
const = 0;
176 virtual void setTime(
const double t )
const = 0;
180 virtual double time()
const = 0;
183 VirtualConservationLawModelMethods(Point)
184 VirtualConservationLawModelMethods(ElementPoint)
185 VirtualConservationLawModelMethods(IntersectionPoint)
186 VirtualConservationLawModelMethods(ElementIntersectionPoint)
189 VirtualConservationLawModelMethods(OriginalPoint)
190 VirtualConservationLawModelMethods(OriginalElementPoint)
191 VirtualConservationLawModelMethods(OriginalIntersectionPoint)
192 VirtualConservationLawModelMethods(OriginalElementIntersectionPoint)
194 VirtualConservationLawModelMethods(LocalDomainType)
196 typedef std::array<int, dimRange> DirichletComponentType;
197 virtual bool hasDirichletBoundary ()
const = 0;
198 virtual bool hasNeumanBoundary ()
const = 0;
199 virtual bool isDirichletIntersection(
const IntersectionType& inter, DirichletComponentType &dirichletComponent )
const = 0;
205 class CheckTimeMethod
208 static std::true_type testSignature(
double& (T::*)());
211 static decltype(testSignature(&T::time)) test(std::nullptr_t);
214 static std::false_type test(...);
216 using type =
decltype(test<M>(
nullptr));
218 static const bool value = type::value;
222 template <
class Model,
bool>
225 static void setTime( Model&,
const double ) {}
226 static double time(
const Model& ) {
return -1.0; }
229 template <
class Model>
230 struct CallSetTime< Model, true >
232 static void setTime( Model& model,
const double t ) { model.time() = t; }
233 static double time (
const Model& model ) {
return model.time(); }
238template <
class ModelImpl >
239struct ConservationLawModelWrapper
240 :
public ConservationLawModel<typename ModelImpl::GridPartType,
241 ModelImpl::dimD, ModelImpl::dimR, typename ModelImpl::RRangeFieldType>
243 typedef ModelImpl Impl;
244 typedef typename ModelImpl::GridPartType GridPartType;
245 static const int dimD = ModelImpl::dimD;
246 static const int dimR = ModelImpl::dimR;
247 typedef ConservationLawModel<GridPartType, dimD, dimR, typename ModelImpl::RRangeFieldType> Base;
249 typedef typename Base::Point Point;
250 typedef typename Base::IntersectionPoint IntersectionPoint;
251 typedef typename Base::ElementPoint ElementPoint;
252 typedef typename Base::ElementIntersectionPoint ElementIntersectionPoint;
253 typedef typename Base::OriginalPoint OriginalPoint;
254 typedef typename Base::OriginalIntersectionPoint OriginalIntersectionPoint;
255 typedef typename Base::OriginalElementPoint OriginalElementPoint;
256 typedef typename Base::OriginalElementIntersectionPoint OriginalElementIntersectionPoint;
257 typedef typename Base::LocalDomainType LocalDomainType;
258 typedef typename Base::DomainType DomainType;
259 typedef typename Base::DRangeType DRangeType;
260 typedef typename Base::DJacobianRangeType DJacobianRangeType;
261 typedef typename Base::DHessianRangeType DHessianRangeType;
262 typedef typename Base::RRangeType RRangeType;
263 typedef typename Base::RJacobianRangeType RJacobianRangeType;
264 typedef typename Base::RHessianRangeType RHessianRangeType;
265 typedef typename Base::EntityType EntityType;
266 typedef typename Base::IntersectionType IntersectionType;
268 template<
class... Args, std::enable_if_t< std::is_constructible< ModelImpl, Args &&... >::value,
int > = 0 >
269 explicit ConservationLawModelWrapper ( Args &&... args )
270 : impl_(
std::forward< Args >( args )... )
273 ~ConservationLawModelWrapper()
278 WrapperConservationLawModelMethods(Point);
279 WrapperConservationLawModelMethods(ElementPoint);
280 WrapperConservationLawModelMethods(IntersectionPoint);
281 WrapperConservationLawModelMethods(ElementIntersectionPoint);
284 WrapperConservationLawModelMethods(OriginalPoint);
285 WrapperConservationLawModelMethods(OriginalElementPoint);
286 WrapperConservationLawModelMethods(OriginalIntersectionPoint);
287 WrapperConservationLawModelMethods(OriginalElementIntersectionPoint);
289 WrapperConservationLawModelMethods(LocalDomainType);
292 virtual std::string name()
const
294 return impl().name();
299 virtual void setTime(
const double t )
const
301 detail::CallSetTime< ModelImpl, detail::CheckTimeMethod< ModelImpl >::value >
302 ::setTime(
const_cast< ModelImpl&
> (impl()), t );
307 virtual double time()
const
309 return detail::CallSetTime< ModelImpl, detail::CheckTimeMethod< ModelImpl >::value >::time( impl() );
312 typedef std::array<int, dimR> DirichletComponentType;
313 virtual bool hasDirichletBoundary ()
const
315 return impl().hasDirichletBoundary();
317 virtual bool hasNeumanBoundary ()
const
319 return impl().hasNeumanBoundary();
321 virtual bool isDirichletIntersection(
const IntersectionType& inter, DirichletComponentType &dirichletComponent )
const
323 return impl().isDirichletIntersection(inter, dirichletComponent);
325 virtual bool init(
const EntityType &entity)
const
327 return impl().init(entity);
329 virtual void unbind()
const
331 return impl().unbind();
333 const ModelImpl &impl()
const
345#define VirtualPenaltyMethods(POINT) \
346 virtual RRangeType penalty ( const POINT &x,\
347 const DRangeType &value ) const = 0;\
348 virtual RRangeType linPenalty ( const POINT &x,\
349 const DRangeType &value ) const = 0;
350#define WrapperPenaltyMethods(POINT) \
351 virtual RRangeType penalty ( const POINT &x,\
352 const DRangeType &value ) const \
353 { return impl().penalty(x, value); } \
354 virtual RRangeType linPenalty ( const POINT &x,\
355 const DRangeType &value ) const \
356 { return impl().linPenalty( x, value); }
358template<
class Gr
idPart,
int dimDomain,
int dimRange=dimDomain,
class RangeField =
double >
359struct DGConservationLawModel
361 typedef GridPart GridPartType;
362 static const int dimD = dimDomain;
363 static const int dimR = dimRange;
364 typedef DGConservationLawModel<GridPartType, dimD, dimR, RangeField> ModelType;
365 typedef RangeField RangeFieldType;
368 GridPart::dimensionworld, dimD > DFunctionSpaceType;
370 GridPart::dimensionworld, dimR > RFunctionSpaceType;
381 typedef typename GridPartType::template Codim<0>::EntityType EntityType;
383 typedef typename EntityType::Geometry::LocalCoordinate LocalDomainType;
386 template <
class F,
int d>
387 using Traits = Dune::FemPy::FempyQuadratureTraits<F,d>;
389 QuadraturePointWrapperType Point;
391 QuadraturePointWrapperType IntersectionPoint;
392 typedef typename Dune::Fem::ElementQuadrature< GridPartType, 0, Traits >::
393 QuadraturePointWrapperType ElementPoint;
394 typedef typename Dune::Fem::ElementQuadrature< GridPartType, 1, Traits >::
395 QuadraturePointWrapperType ElementIntersectionPoint;
399 QuadraturePointWrapperType OriginalPoint;
401 QuadraturePointWrapperType OriginalIntersectionPoint;
402 typedef typename Dune::Fem::ElementQuadrature< GridPartType, 0 >::
403 QuadraturePointWrapperType OriginalElementPoint;
404 typedef typename Dune::Fem::ElementQuadrature< GridPartType, 1 >::
405 QuadraturePointWrapperType OriginalElementIntersectionPoint;
413 DGConservationLawModel( )
415 virtual ~DGConservationLawModel() {}
417 virtual std::string name()
const = 0;
419 virtual bool init(
const EntityType &entity)
const = 0;
420 virtual void unbind()
const = 0;
424 virtual void setTime(
const double t )
const = 0;
428 virtual double time()
const = 0;
431 VirtualConservationLawModelMethods(Point)
432 VirtualConservationLawModelMethods(ElementPoint)
433 VirtualConservationLawModelMethods(IntersectionPoint)
434 VirtualConservationLawModelMethods(ElementIntersectionPoint)
437 VirtualConservationLawModelMethods(OriginalPoint)
438 VirtualConservationLawModelMethods(OriginalElementPoint)
439 VirtualConservationLawModelMethods(OriginalIntersectionPoint)
440 VirtualConservationLawModelMethods(OriginalElementIntersectionPoint)
442 VirtualConservationLawModelMethods(LocalDomainType)
444 typedef std::array<int, dimRange> DirichletComponentType;
445 virtual bool hasDirichletBoundary ()
const = 0;
446 virtual bool hasNeumanBoundary ()
const = 0;
447 virtual bool isDirichletIntersection(
const IntersectionType& inter, DirichletComponentType &dirichletComponent )
const = 0;
449 VirtualPenaltyMethods(Point)
450 VirtualPenaltyMethods(ElementPoint)
451 VirtualPenaltyMethods(IntersectionPoint)
452 VirtualPenaltyMethods(ElementIntersectionPoint)
455 VirtualPenaltyMethods(OriginalPoint)
456 VirtualPenaltyMethods(OriginalElementPoint)
457 VirtualPenaltyMethods(OriginalIntersectionPoint)
458 VirtualPenaltyMethods(OriginalElementIntersectionPoint)
460 VirtualPenaltyMethods(LocalDomainType)
464template <
class ModelImpl >
465struct DGConservationLawModelWrapper
466:
public DGConservationLawModel<typename ModelImpl::GridPartType, ModelImpl::dimD, ModelImpl::dimR, typename ModelImpl::RRangeFieldType>
468 typedef ModelImpl Impl;
469 typedef typename ModelImpl::GridPartType GridPartType;
470 static const int dimD = ModelImpl::dimD;
471 static const int dimR = ModelImpl::dimR;
472 typedef DGConservationLawModel<GridPartType, dimD, dimR, typename ModelImpl::RRangeFieldType> Base;
474 typedef typename Base::Point Point;
475 typedef typename Base::IntersectionPoint IntersectionPoint;
476 typedef typename Base::ElementPoint ElementPoint;
477 typedef typename Base::ElementIntersectionPoint ElementIntersectionPoint;
478 typedef typename Base::OriginalPoint OriginalPoint;
479 typedef typename Base::OriginalIntersectionPoint OriginalIntersectionPoint;
480 typedef typename Base::OriginalElementPoint OriginalElementPoint;
481 typedef typename Base::OriginalElementIntersectionPoint OriginalElementIntersectionPoint;
482 typedef typename Base::LocalDomainType LocalDomainType;
483 typedef typename Base::DomainType DomainType;
484 typedef typename Base::DRangeType DRangeType;
485 typedef typename Base::DJacobianRangeType DJacobianRangeType;
486 typedef typename Base::DHessianRangeType DHessianRangeType;
487 typedef typename Base::RRangeType RRangeType;
488 typedef typename Base::RJacobianRangeType RJacobianRangeType;
489 typedef typename Base::RHessianRangeType RHessianRangeType;
490 typedef typename Base::EntityType EntityType;
491 typedef typename Base::IntersectionType IntersectionType;
493 template<
class... Args, std::enable_if_t< std::is_constructible< ModelImpl, Args &&... >::value,
int > = 0 >
494 explicit DGConservationLawModelWrapper ( Args &&... args )
495 : impl_(
std::forward< Args >( args )... )
498 ~DGConservationLawModelWrapper()
503 WrapperConservationLawModelMethods(Point);
504 WrapperConservationLawModelMethods(ElementPoint);
505 WrapperConservationLawModelMethods(IntersectionPoint);
506 WrapperConservationLawModelMethods(ElementIntersectionPoint);
509 WrapperConservationLawModelMethods(OriginalPoint);
510 WrapperConservationLawModelMethods(OriginalElementPoint);
511 WrapperConservationLawModelMethods(OriginalIntersectionPoint);
512 WrapperConservationLawModelMethods(OriginalElementIntersectionPoint);
514 WrapperConservationLawModelMethods(LocalDomainType);
516 WrapperPenaltyMethods(Point)
517 WrapperPenaltyMethods(ElementPoint)
518 WrapperPenaltyMethods(IntersectionPoint)
519 WrapperPenaltyMethods(ElementIntersectionPoint)
522 WrapperPenaltyMethods(OriginalPoint);
523 WrapperPenaltyMethods(OriginalElementPoint);
524 WrapperPenaltyMethods(OriginalIntersectionPoint);
525 WrapperPenaltyMethods(OriginalElementIntersectionPoint);
527 WrapperPenaltyMethods(LocalDomainType);
529 virtual std::string name()
const
531 return impl().name();
536 virtual void setTime(
const double t )
const
538 detail::CallSetTime< ModelImpl, detail::CheckTimeMethod< ModelImpl >::value >
539 ::setTime(
const_cast< ModelImpl&
> (impl()), t );
544 virtual double time()
const
546 return detail::CallSetTime< ModelImpl, detail::CheckTimeMethod< ModelImpl >::value >::time( impl() );
549 typedef std::array<int, dimR> DirichletComponentType;
550 virtual bool hasDirichletBoundary ()
const
552 return impl().hasDirichletBoundary();
554 virtual bool hasNeumanBoundary ()
const
556 return impl().hasNeumanBoundary();
558 virtual bool isDirichletIntersection(
const IntersectionType& inter, DirichletComponentType &dirichletComponent )
const
560 return impl().isDirichletIntersection(inter, dirichletComponent);
562 virtual bool init(
const EntityType &entity)
const
564 return impl().init(entity);
566 virtual void unbind()
const
568 return impl().unbind();
570 const ModelImpl &impl()
const
IntersectionIteratorType::Intersection IntersectionType
type of intersection
Definition: adaptiveleafgridpart.hh:95
quadrature class supporting base function caching
Definition: cachingquadrature.hh:101
ExplicitFieldVector< FieldMatrix< RangeFieldType, dimDomain, dimDomain >, dimRange > HessianRangeType
Intrinsic type used for the hessian values has a Dune::FieldMatrix type interface.
Definition: functionspaceinterface.hh:79
FunctionSpaceTraits::DomainFieldType DomainFieldType
Intrinsic type used for values in the domain field (usually a double)
Definition: functionspaceinterface.hh:60
FunctionSpaceTraits::RangeType RangeType
Type of range vector (using type of range field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:71
FunctionSpaceTraits::LinearMappingType JacobianRangeType
Intrinsic type used for the jacobian values has a Dune::FieldMatrix type interface.
Definition: functionspaceinterface.hh:75
FunctionSpaceTraits::DomainType DomainType
Type of domain vector (using type of domain field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:67
A vector valued function space.
Definition: functionspace.hh:60
Definition of macros controlling symbol visibility at the ABI level.