Implement handlers to allow generation of correct penalty term

applications/convection_diffusion/conv_diff_dg.ufl

+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]

applications/poisson_dg/poisson_dg.ufl

+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]

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)
 

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
     #

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

python/dune/perftool/ufl/preprocess.py

 """ 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,
                                                           )
                                  )
 

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

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]

test/poisson/poisson_dg.ufl

 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]

test/poisson/poisson_dg_neumann.ufl

+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]

test/poisson/poisson_dg_tensor.ufl

+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]