diff --git a/dune/fem-dg/models/modelwrapper.hh b/dune/fem-dg/models/modelwrapper.hh
index 7009d6e46274901c8c0fd67ffbd0c9937b437e1e..2cabcf41971ec88485b869ea31d63316a2921c06 100644
--- a/dune/fem-dg/models/modelwrapper.hh
+++ b/dune/fem-dg/models/modelwrapper.hh
@@ -240,8 +240,8 @@ namespace Fem
void obtainBounds( RangeType& globalMin, RangeType& globalMax) const
{
- globalMin = 0;
- globalMax = std::numeric_limits<double>::max();
+ assert( hasAdvection );
+ advection().obtainBounds(globalMin, globalMax);
}
bool isConstant( const RangeType& min, const RangeType& max ) const
diff --git a/dune/fem-dg/operator/limiter/limitpass.hh b/dune/fem-dg/operator/limiter/limitpass.hh
index 6caab8ed03a9aea60244d01fb001cfe39e77c6d8..b83ef10637cef4d2c4bb8a07e03265683ab6d525 100644
--- a/dune/fem-dg/operator/limiter/limitpass.hh
+++ b/dune/fem-dg/operator/limiter/limitpass.hh
@@ -418,7 +418,7 @@ namespace Fem
storedComboSets_(),
tolFactor_( getTolFactor() ),
indicator_( getIndicator(parameter) ),
- tol_1_( 1.0/ parameter.getValue("femdg.limiter.tolerance", double(0.5) ) ),
+ tol_1_( 2./ parameter.getValue("femdg.limiter.tolerance", double(1.0) ) ),
geoInfo_( gridPart_.indexSet() ),
faceGeoInfo_( geoInfo_.geomTypes(1) ),
phi0_( 0 ),
diff --git a/python/dune/femdg/_operators.py b/python/dune/femdg/_operators.py
index 56afc95edeb8f66f3ae42dfc383a605eac84e5e0..6fc8e6e0dd9d6c99320e857814c6245412bb425a 100644
--- a/python/dune/femdg/_operators.py
+++ b/python/dune/femdg/_operators.py
@@ -294,8 +294,7 @@ def femDGOperator(Model, space,
limitedDimRange = "FunctionSpace :: dimRange"
else:
limiterModified = {}
- count = 0
- for id,f in limiterModifiedDict.items(): count += 1
+ count = len( limiterModifiedDict.items() )
limitedDimRange = str(count)
struct.append([Declaration(
Variable("const int", "limitedDimRange = " + limitedDimRange),
diff --git a/python/dune/femdg/patch.py b/python/dune/femdg/patch.py
index f99bc3fbdb91d7580dc0f70bc84d8b48f5dd47c5..1279a99a91522dacec6e799ae3691e9bb3e75abd 100644
--- a/python/dune/femdg/patch.py
+++ b/python/dune/femdg/patch.py
@@ -1,9 +1,12 @@
-from __future__ import division, print_function, unicode_literals
-# from dune.ufl.codegen import generateMethod
-from ufl import grad, TrialFunction, SpatialCoordinate, FacetNormal, Coefficient, replace, diff, as_vector
+from functools import reduce
+
+from ufl import grad, TrialFunction, SpatialCoordinate, FacetNormal,\
+ Coefficient, replace, diff, as_vector,\
+ conditional
from ufl.core.expr import Expr
-from dune.source.cplusplus import Variable, UnformattedExpression, AccessModifier, Declaration
+from dune.source.cplusplus import Variable, UnformattedExpression,\
+ AccessModifier, Declaration, Method
from ufl.algorithms import expand_compounds, expand_derivatives, expand_indices, expand_derivatives
def uflExpr(Model,space,t):
@@ -11,6 +14,32 @@ def uflExpr(Model,space,t):
n = FacetNormal(space.cell())
x = SpatialCoordinate(space.cell())
+ physicalBound = None
+ lowerBound = getattr(Model,"lowerBound",None)
+ if lowerBound is None:
+ lowerBound = space.dimRange*["std::numeric_limits<double>::min()"]
+ else:
+ physicalBound = reduce( (lambda x,y: x*y),
+ [conditional(u[i]>=lowerBound[i],1,0)
+ for i in range(len(lowerBound))
+ if lowerBound[i] is not None
+ ], None )
+ if physicalBound == None: lowerBound=space.dimRange*["std::numeric_limits<double>::min()"]
+ upperBound = getattr(Model,"upperBound",None)
+ if upperBound is None:
+ upperBound = space.dimRange*["std::numeric_limits<double>::max()"]
+ else:
+ physicalBound_ = reduce( (lambda x,y: x*y),
+ [conditional(u[i]>=upperBound[i],1,0)
+ for i in range(len(upperBound))
+ if upperBound[i] is not None
+ ], None )
+ if physicalBound_ == None: lowerBound=space.dimRange*["std::numeric_limits<double>::max()"]
+ elif physicalBound_ is not None:
+ physicalBound = physicalBound_ if physicalBound is None else\
+ physicalBound*physicalBound_
+
+ if physicalBound == [None]: physicalBound = None
maxSpeed = getattr(Model,"maxLambda",None)
if maxSpeed is not None:
maxSpeed = maxSpeed(t,x,u,n)
@@ -20,9 +49,14 @@ def uflExpr(Model,space,t):
diffusionTimeStep = getattr(Model,"maxDiffusion",None)
if diffusionTimeStep is not None:
diffusionTimeStep = diffusionTimeStep(t,x,u)
- physical = getattr(Model,"physical",None)
- if physical is not None:
+ physical = getattr(Model,"physical",True)
+ if not isinstance(physical,bool):
physical = physical(t,x,u)
+ if physicalBound is not None:
+ physical = physical*physicalBound
+ else:
+ if physicalBound is not None:
+ physical = physical*physicalBound
# TODO: jump is problematic since the coefficient 'w' causes issues -
# no predefined when extracting required Coefficients so needs fixing
# So jump is not returned by this method and is constructed again in
@@ -78,12 +112,13 @@ def uflExpr(Model,space,t):
limitedDimRange = "DFunctionSpace :: dimRange"
else:
limiterModified = {}
- count = 0
- for id,f in limiterModifiedDict.items(): count += 1
+ count = len(limiterModifiedDict.items())
+ # for id,f in limiterModifiedDict.items(): count += 1
limitedDimRange = str(count)
# jump = None # TODO: see comment above
return maxSpeed, velocity, diffusionTimeStep, physical, jump,\
- boundaryAFlux, boundaryDFlux, boundaryValue, hasBoundaryValue
+ boundaryAFlux, boundaryDFlux, boundaryValue, hasBoundaryValue,\
+ physicalBound
def codeFemDg(self):
code = self._code()
@@ -100,6 +135,31 @@ def codeFemDg(self):
n = FacetNormal(space.cell())
x = SpatialCoordinate(space.cell())
+ physicalBound = None
+ lowerBound = getattr(Model,"lowerBound",None)
+ if lowerBound is None:
+ lowerBound = space.dimRange*["std::numeric_limits<double>::min()"]
+ else:
+ physicalBound = reduce( (lambda x,y: x*y),
+ [conditional(u[i]>=lowerBound[i],1,0)
+ for i in range(len(lowerBound))
+ if lowerBound[i] is not None
+ ], None )
+ if physicalBound == None: lowerBound=space.dimRange*["std::numeric_limits<double>::min()"]
+ upperBound = getattr(Model,"upperBound",None)
+ if upperBound is None:
+ upperBound = space.dimRange*["std::numeric_limits<double>::max()"]
+ else:
+ physicalBound_ = reduce( (lambda x,y: x*y),
+ [conditional(u[i]>=upperBound[i],1,0)
+ for i in range(len(upperBound))
+ if upperBound[i] is not None
+ ], None )
+ if physicalBound_ == None: upperBound=space.dimRange*["std::numeric_limits<double>::max()"]
+ elif physicalBound_ is not None:
+ physicalBound = physicalBound_ if physicalBound is None else\
+ physicalBound*physicalBound_
+
# TODO come up with something better!
hasGamma = getattr(Model,"gamma",None)
code.append([Declaration(
@@ -158,10 +218,15 @@ def codeFemDg(self):
targs=['class Entity, class Point, class T'], const=True,
predefined=predefined)
- # QUESTION: should `physical` actually depend on x? Perhaps even t?
physical = getattr(Model,"physical",True)
if not isinstance(physical,bool):
physical = physical(t,x,u)
+ if physicalBound is not None:
+ physical = physical*physicalBound
+ else:
+ if physicalBound is not None:
+ physical = physical*physicalBound
+
self.generateMethod(code, physical,
'double', 'physical',
args=['const Entity &entity', 'const Point &x',
@@ -285,15 +350,24 @@ def codeFemDg(self):
limitedDimRange = "DFunctionSpace :: dimRange"
else:
limiterModified = {}
- count = 0
- for id,f in limiterModifiedDict.items(): count += 1
+ count = len(limiterModifiedDict.items())
+ # for id,f in limiterModifiedDict.items(): count += 1
limitedDimRange = str(count)
self.generateMethod(code, limiterModified,
'void', 'limitedRange',
args=['LimitedRange& limRange'],
targs=['class LimitedRange'], const=True, evalSwitch=False,
predefined=predefined)
-
+ obtainBounds = Method("void","obtainBounds",
+ targs=['class LimitedRange'],
+ args=["LimitedRange& globalMin","LimitedRange& globalMax"],const=True)
+ obtainBounds.append("globalMin = {"
+ + ', '.join(map(str, lowerBound))
+ +"};")
+ obtainBounds.append("globalMax = {"
+ + ', '.join(map(str, upperBound))
+ +"};")
+ code.append( [obtainBounds] )
return code
def transform(Model,space,t, name="" ):