diff --git a/applications/convection_diffusion/conv_diff_dg.ufl b/applications/convection_diffusion/conv_diff_dg.ufl
index 38e8dc7334af54b167a51955186f63660af12e4c..80e9205514aa0d7c1d813b2abb825fb654c6443f 100644
--- a/applications/convection_diffusion/conv_diff_dg.ufl
+++ b/applications/convection_diffusion/conv_diff_dg.ufl
@@ -1,3 +1,4 @@
+dim = 3
x = SpatialCoordinate(cell)
g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
@@ -13,17 +14,22 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+
theta = -1.0
r = (inner(A*grad(u), grad(v)) + (c*u-f)*v)*dx \
+ inner(n, A*avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
+ theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
- inner(n, A*grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
- theta*u*inner(A*grad(v), n)*ds \
+ theta*g*inner(A*grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]
diff --git a/applications/poisson_dg/poisson_dg.ufl b/applications/poisson_dg/poisson_dg.ufl
index 536360533793fe795d9b5bd303e023f9ecb31c19..7d532811958c1abc947c71cdfacdec57e0bc3ff0 100644
--- a/applications/poisson_dg/poisson_dg.ufl
+++ b/applications/poisson_dg/poisson_dg.ufl
@@ -1,3 +1,4 @@
+dim = 3
x = SpatialCoordinate(cell)
f = -6.
g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
@@ -9,18 +10,23 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+
theta = 1.0
r = inner(grad(u), grad(v))*dx \
+ inner(n, avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
+ theta*jump(u)*inner(avg(grad(v)), n)*dS \
- inner(n, grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
- theta*u*inner(grad(v), n)*ds \
- f*v*dx \
+ theta*g*inner(grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]
diff --git a/python/dune/perftool/loopy/target.py b/python/dune/perftool/loopy/target.py
index 9427a401e8b0a51f8b43ddf52abfd2cc3b9bc753..6ef27648aa9d659ab63838d84798f19eec3a91f3 100644
--- a/python/dune/perftool/loopy/target.py
+++ b/python/dune/perftool/loopy/target.py
@@ -118,7 +118,7 @@ class DuneCExpressionToCodeMapper(CExpressionToCodeMapper):
result = self.rec(children.pop(), PREC_NONE)
while children:
- result = "std::%s(%s, %s)" % (what,
+ result = "%s(%s, %s)" % (what,
self.rec(children.pop(), PREC_NONE),
result)
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 8529eb1e55186f52c379337be47165800cf21cee..dcc17068303501fc01d96e842089454907342c3b 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -11,6 +11,8 @@ from dune.perftool.pdelab.basis import (pymbolic_basis,
pymbolic_reference_gradient,
)
from dune.perftool.pdelab.geometry import (dimension_iname,
+ name_cell_volume,
+ name_facet_area,
name_facet_jacobian_determinant,
name_jacobian_determinant,
name_jacobian_inverse_transposed,
@@ -112,6 +114,12 @@ class PDELabInterface(object):
def name_unit_outer_normal(self):
return name_unit_outer_normal()
+ def name_cell_volume(self, restriction):
+ return name_cell_volume(restriction)
+
+ def name_facet_area(self):
+ return name_facet_area()
+
#
# Quadrature related generator functions
#
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index b66f58fc875370737ec9a9d267b53acedcb1efe6..ce080eba080be2cf75e5d4ee6993232dfb237079 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -405,3 +405,33 @@ def to_global(local):
name = get_pymbolic_basename(local) + "_global"
apply_to_global_transformation(name, local)
return prim.Variable(name)
+
+
+def define_cell_volume(name, restriction):
+ geo = name_cell_geometry(restriction)
+ temporary_variable(name, shape=())
+ code = "{} = {}.volume();".format(name, geo)
+ return quadrature_preamble(code,
+ assignees=frozenset({name}),
+ )
+
+
+def name_cell_volume(restriction):
+ name = restricted_name("volume", restriction)
+ define_cell_volume(name, restriction)
+ return name
+
+
+def define_facet_area(name):
+ geo = name_intersection_geometry()
+ temporary_variable(name, shape=())
+ code = "{} = {}.volume();".format(name, geo)
+ return quadrature_preamble(code,
+ assignees=frozenset({name}),
+ )
+
+
+def name_facet_area():
+ name = "area"
+ define_facet_area(name)
+ return name
diff --git a/python/dune/perftool/ufl/preprocess.py b/python/dune/perftool/ufl/preprocess.py
index 540300bb3ee61abaa89b6b958df65e879af69e46..913bf5dad763483e56d619f9be6226b92ee9ffd8 100644
--- a/python/dune/perftool/ufl/preprocess.py
+++ b/python/dune/perftool/ufl/preprocess.py
@@ -1,22 +1,20 @@
""" Preprocessing algorithms for UFL forms """
+import ufl.classes as uc
-def preprocess_form(form):
- from ufl.classes import (FacetJacobianDeterminant,
- FacetNormal,
- JacobianDeterminant,
- JacobianInverse,
- )
+def preprocess_form(form):
from ufl.algorithms import compute_form_data
formdata = compute_form_data(form,
do_apply_function_pullbacks=True,
do_apply_geometry_lowering=True,
do_apply_integral_scaling=True,
- preserve_geometry_types=(FacetJacobianDeterminant,
- FacetNormal,
- JacobianDeterminant,
- JacobianInverse,
+ preserve_geometry_types=(uc.CellVolume,
+ uc.FacetArea,
+ uc.FacetJacobianDeterminant,
+ uc.FacetNormal,
+ uc.JacobianDeterminant,
+ uc.JacobianInverse,
)
)
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index 1ba817c856c7ba76fe36b40bf9ec4b94552a5b89..de92f3a49f23346324a280c0e5b1565563adaf3a 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -265,6 +265,12 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
raise NotImplementedError("Power function not really implemented")
+ def max_value(self, o):
+ return prim.Max(tuple(self.call(op) for op in o.ufl_operands))
+
+ def min_value(self, o):
+ return prim.Min(tuple(self.call(op) for op in o.ufl_operands))
+
#
# Handler for conditionals, use pymbolic base implementation
#
@@ -360,6 +366,16 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
def facet_jacobian_determinant(self, o):
return Variable(self.interface.name_facet_jacobian_determinant())
+ def cell_volume(self, o):
+ restriction = self.restriction
+ if self.measure == 'exterior_facet':
+ restriction = Restriction.NEGATIVE
+
+ return prim.Variable(self.interface.name_cell_volume(restriction))
+
+ def facet_area(self, o):
+ return prim.Variable(self.interface.name_facet_area())
+
#
# Equality/Hashability of the visitor
# In order to pass this visitor into (cached) generator functions, it is beneficial
diff --git a/test/poisson/opcount_poisson_dg.ufl b/test/poisson/opcount_poisson_dg.ufl
index 3ced1b8263f75d01feeb90acae9a64ae1bcba482..ed537bc936a03fb1fa21a87419d4c62f67246d4f 100644
--- a/test/poisson/opcount_poisson_dg.ufl
+++ b/test/poisson/opcount_poisson_dg.ufl
@@ -11,20 +11,24 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
# penalty factor
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
theta = 1.0
r = inner(grad(u), grad(v))*dx \
+ inner(n, avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
- theta*jump(u)*inner(avg(grad(v)), n)*dS \
- inner(n, grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
+ theta*u*inner(grad(v), n)*ds \
- f*v*dx \
- theta*g*inner(grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]
diff --git a/test/poisson/poisson_dg.ufl b/test/poisson/poisson_dg.ufl
index daf2793ee5c327dd389d05fffda0c7c8dc9dd244..b8ad04251f4fb94d483c7c01b7430fb3fc440e88 100644
--- a/test/poisson/poisson_dg.ufl
+++ b/test/poisson/poisson_dg.ufl
@@ -1,11 +1,13 @@
cell = triangle
+degree = 1
+dim = 2
x = SpatialCoordinate(cell)
c = (0.5-x[0])**2 + (0.5-x[1])**2
g = exp(-1.*c)
f = 4*(1.-c)*g
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
u = TrialFunction(V)
v = TestFunction(V)
@@ -13,20 +15,24 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
# penalty factor
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
theta = 1.0
r = inner(grad(u), grad(v))*dx \
+ inner(n, avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
- theta*jump(u)*inner(avg(grad(v)), n)*dS \
- inner(n, grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
+ theta*u*inner(grad(v), n)*ds \
- f*v*dx \
- theta*g*inner(grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]
diff --git a/test/poisson/poisson_dg_neumann.ufl b/test/poisson/poisson_dg_neumann.ufl
index 1dd15d4f39c18d2b07567c3d8d4afe46dbd72a18..884787e6374a742ab5276a792f6975019a65f747 100644
--- a/test/poisson/poisson_dg_neumann.ufl
+++ b/test/poisson/poisson_dg_neumann.ufl
@@ -1,3 +1,5 @@
+dim = 3
+degree = 1
cell = triangle
x = SpatialCoordinate(cell)
@@ -8,7 +10,7 @@ sgn = conditional(x[1] > 0.5, 1., -1.)
j = -2.*sgn*(x[1]-0.5)*g
bctype = conditional(Or(x[1]<1e-8, x[1]>1.-1e-8), 0, 1)
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
u = TrialFunction(V)
v = TestFunction(V)
@@ -18,21 +20,25 @@ n = FacetNormal(cell)('+')
ds = ds(subdomain_data=bctype)
# penalty factor
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
theta = 1.0
r = inner(grad(u), grad(v))*dx \
+ inner(n, avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
- theta*jump(u)*inner(avg(grad(v)), n)*dS \
- inner(n, grad(u))*v*ds(1) \
- + gamma*u*v*ds(1) \
+ + gamma_ext*u*v*ds(1) \
+ theta*u*inner(grad(v), n)*ds(1) \
- f*v*dx \
- theta*g*inner(grad(v), n)*ds(1) \
- - gamma*g*v*ds(1) \
+ - gamma_ext*g*v*ds(1) \
- j*v*ds(0)
forms = [r]
diff --git a/test/poisson/poisson_dg_tensor.ufl b/test/poisson/poisson_dg_tensor.ufl
index d84cebf6d1420b8ec3fd35086d40833937e55f60..1ae0139b3d63fe588b69ce221e17a488b0e1535a 100644
--- a/test/poisson/poisson_dg_tensor.ufl
+++ b/test/poisson/poisson_dg_tensor.ufl
@@ -1,3 +1,5 @@
+dim = 2
+degree = 1
cell = quadrilateral
x = SpatialCoordinate(cell)
@@ -8,7 +10,7 @@ g = x[0]**2 + x[1]**2
c = 8.
f = -4.
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
u = TrialFunction(V)
v = TestFunction(V)
@@ -16,20 +18,24 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
# penalty factor
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
theta = 1.0
r = (inner(A*grad(u), grad(v)) + c*u*v)*dx \
+ inner(n, A*avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
- theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
- inner(n, A*grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
+ theta*u*inner(A*grad(v), n)*ds \
- f*v*dx \
- theta*g*inner(A*grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]