diff --git a/applications/stokes_dg/stokes_dg.ufl b/applications/stokes_dg/stokes_dg.ufl
index 6189ae1be86f2a1e582e9bb6ad5c9acef83c47e2..d36556dbaf38b76dbc565437a28de38889b62d0f 100644
--- a/applications/stokes_dg/stokes_dg.ufl
+++ b/applications/stokes_dg/stokes_dg.ufl
@@ -1,4 +1,5 @@
 cell = hexahedron
+dim = 3
 
 x = SpatialCoordinate(cell)
 g_v = as_vector((4.*x[1]*(1.-x[1]), 0.0, 0.0))
@@ -14,21 +15,28 @@ u, p = TrialFunctions(TH)
 ds = ds(subdomain_id=1, subdomain_data=bctype)
 
 n = FacetNormal(cell)('+')
-eps = -1.0
-sigma = 1.0
-h_e = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+
+# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
+theta = -1.0
+
+# penalty factor
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * v_degree * (v_degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * v_degree * (v_degree + dim - 1)) / h_int
 
 r = inner(grad(u), grad(v))*dx \
   - p*div(v)*dx \
   - q*div(u)*dx \
-  + inner(avg(grad(u))*n, jump(v))*dS \
-  + sigma / h_e * inner(jump(u), jump(v))*dS \
-  - eps * inner(avg(grad(v))*n, jump(u))*dS \
-  - avg(p)*inner(jump(v), n)*dS \
-  - avg(q)*inner(jump(u), n)*dS \
+  - inner(avg(grad(u))*n, jump(v))*dS \
+  + gamma_int * inner(jump(u), jump(v))*dS \
+  + theta * inner(avg(grad(v))*n, jump(u))*dS \
+  + avg(p)*inner(jump(v), n)*dS \
+  + avg(q)*inner(jump(u), n)*dS \
   - inner(grad(u)*n, v)*ds \
-  + sigma / h_e * inner(u-g_v, v)*ds \
-  + eps * inner(grad(v)*n, u-g_v)*ds \
+  + gamma_ext * inner(u-g_v, v)*ds \
+  + theta * inner(grad(v)*n, u-g_v)*ds \
   + p*inner(v, n)*ds \
   + q*inner(u-g_v, n)*ds
 
diff --git a/python/dune/perftool/blockstructured/tools.py b/python/dune/perftool/blockstructured/tools.py
index 31f57e70f7bfbdafc061a0b248e6614a7ee6480c..029b0f752b6c5f37854550588070f6072a7de977 100644
--- a/python/dune/perftool/blockstructured/tools.py
+++ b/python/dune/perftool/blockstructured/tools.py
@@ -38,7 +38,7 @@ def sub_element_inames():
 def sub_facet_inames():
     subelem_inames = sub_element_inames()
 
-    center = pymbolic_in_cell_coordinates(prim.Variable(name_localcenter()), Restriction.NEGATIVE)
+    center = pymbolic_in_cell_coordinates(prim.Variable(name_localcenter()), Restriction.POSITIVE)
 
     # check if iname[index] must be constant or not
     def predicate(index):
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index e4dfb972a48b34428b54bef25f142eafd2c86866..30449edea105a0a962c5883559e61ced1fc4cdb0 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -177,14 +177,14 @@ def name_accumulation_variable(restrictions=None):
             if measure == "cell":
                 restrictions = (Restriction.NONE,)
             else:
-                restrictions = (Restriction.NEGATIVE,)
+                restrictions = (Restriction.POSITIVE,)
         return name_residual(*restrictions)
     if ft == 'jacobian':
         if restrictions is None:
             if measure == "cell":
                 restrictions = (Restriction.NONE, Restriction.NONE)
             else:
-                restrictions = (Restriction.NEGATIVE, Restriction.NEGATIVE)
+                restrictions = (Restriction.POSITIVE, Restriction.POSITIVE)
         return name_jacobian(*restrictions)
     assert False
 
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index d536027362ee2ef347ad1847934c50fb94cf3fc4..c382faeea09e13208aa7155b120fb070789c80e6 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -112,7 +112,7 @@ def type_geometry_wrapper():
 @preamble
 def define_restricted_cell(name, restriction):
     ig = name_intersection_geometry_wrapper()
-    which = "inside" if restriction == Restriction.NEGATIVE else "outside"
+    which = "inside" if restriction == Restriction.POSITIVE else "outside"
     return "const auto& {} = {}.{}();".format(name,
                                               ig,
                                               which,
@@ -124,7 +124,7 @@ def name_cell(restriction):
         eg = name_element_geometry_wrapper()
         return "{}.entity()".format(eg)
     else:
-        which = "inside" if restriction == Restriction.NEGATIVE else "outside"
+        which = "inside" if restriction == Restriction.POSITIVE else "outside"
         name = "{}_cell".format(which)
         define_restricted_cell(name, restriction)
         return name
@@ -186,7 +186,7 @@ def name_geometry():
 @preamble
 def define_in_cell_geometry(restriction, name):
     ig = name_intersection_geometry_wrapper()
-    which = "In" if restriction == Restriction.NEGATIVE else "Out"
+    which = "In" if restriction == Restriction.POSITIVE else "Out"
     return "auto {} = {}.geometryIn{}side();".format(name,
                                                      ig,
                                                      which
@@ -196,7 +196,7 @@ def define_in_cell_geometry(restriction, name):
 def name_in_cell_geometry(restriction):
     assert restriction is not Restriction.NONE
 
-    name = "geo_in_{}side".format("in" if restriction is Restriction.NEGATIVE else "out")
+    name = "geo_in_{}side".format("in" if restriction is Restriction.POSITIVE else "out")
     define_in_cell_geometry(restriction, name)
     return name
 
@@ -216,7 +216,7 @@ def apply_in_cell_transformation(name, local, restriction):
 
 def pymbolic_in_cell_coordinates(local, restriction):
     basename = get_pymbolic_basename(local)
-    name = "{}_in_{}side".format(basename, "in" if restriction is Restriction.NEGATIVE else "out")
+    name = "{}_in_{}side".format(basename, "in" if restriction is Restriction.POSITIVE else "out")
     temporary_variable(name, shape=(world_dimension(),), shape_impl=("fv",))
     apply_in_cell_transformation(name, local, restriction)
     return Variable(name)
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 03496a5cfbc24b0191e29df5bf83bac13fa99bc9..1a9d9c4661af5e4454e81098a22cbe1797b9eaa0 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -161,9 +161,9 @@ def localoperator_basename(form_ident):
 
 def name_gridfunction_member(coeff, restriction, diffOrder=0):
     # We reuse the grid function for volume integrals in skeleton integrals
-    if restriction == Restriction.NEGATIVE:
+    if restriction == Restriction.POSITIVE:
         restriction = Restriction.NONE
-    restr = "_n" if restriction == Restriction.POSITIVE else ""
+    restr = "_n" if restriction == Restriction.NEGATIVE else ""
     name = "local_gridfunction_coeff{}_diff{}{}".format(coeff.count(), diffOrder, restr)
     define_gridfunction_member(name, coeff, restriction, diffOrder)
     return name
@@ -347,13 +347,12 @@ def _list_infos(expr, number, visitor):
         return
     element = ma[0].argexpr.ufl_element()
 
-    from dune.perftool.ufl.modified_terminals import Restriction
     if visitor.measure == "cell":
         restrictions = (Restriction.NONE,)
     elif visitor.measure == "exterior_facet":
-        restrictions = (Restriction.NEGATIVE,)
+        restrictions = (Restriction.POSITIVE,)
     elif visitor.measure == "interior_facet":
-        restrictions = (Restriction.NEGATIVE, Restriction.POSITIVE)
+        restrictions = (Restriction.POSITIVE, Restriction.NEGATIVE)
     for res in restrictions:
         for ei in range(element.value_size()):
             yield PDELabAccumulationInfo(element_index=ei, restriction=res)
@@ -378,8 +377,7 @@ def get_accumulation_info(expr, visitor):
 
     restriction = visitor.restriction
     if visitor.measure == 'exterior_facet':
-        from dune.perftool.pdelab.restriction import Restriction
-        restriction = Restriction.NEGATIVE
+        restriction = Restriction.POSITIVE
 
     inames = visitor.interface.lfs_inames(leaf_element,
                                           restriction,
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 19339acfd315038f8e7d03ed8409265084fad2f4..031d97b019df0cab8355d4295b5229584e077b9b 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -15,7 +15,6 @@ from dune.perftool.generation import (backend,
                                       )
 from dune.perftool.pdelab.localoperator import lop_template_range_field
 from dune.perftool.options import get_form_option
-from dune.perftool.ufl.modified_terminals import Restriction
 
 from pymbolic.primitives import Variable, Subscript
 
diff --git a/python/dune/perftool/pdelab/restriction.py b/python/dune/perftool/pdelab/restriction.py
index 7d77a17b6c7e14107a359f711d110584b1a83cb2..03c55eaec4874c2611d892d74d8a14476cf8435d 100644
--- a/python/dune/perftool/pdelab/restriction.py
+++ b/python/dune/perftool/pdelab/restriction.py
@@ -2,9 +2,24 @@ from dune.perftool.ufl.modified_terminals import Restriction
 
 
 def restricted_name(name, restriction):
+    """Adapt name according to the restictrion
+
+    Some remarks:
+
+    - UFL defines the jump the following: jump(v) = v('+') - v('-').
+
+    - The corresponding outer normal vector is n =
+      FacetNormal(cell)('+'). The user needs to make the right choice
+      in the UFL file.
+
+    - In the literature this convention is sometimes swapped. In order
+      to be consistent with UFL we choose ('+') as self and ('-') as
+      neighbor and choose the outer unit normal vector accordingly.
+
+    """
     if restriction == Restriction.NONE:
         return name
     if restriction == Restriction.POSITIVE:
-        return name + '_n'
-    if restriction == Restriction.NEGATIVE:
         return name + '_s'
+    if restriction == Restriction.NEGATIVE:
+        return name + '_n'
diff --git a/python/dune/perftool/pdelab/signatures.py b/python/dune/perftool/pdelab/signatures.py
index 6eb30cd11b977e4b6aa930dd86558219e9da3280..2990c8b24b9c89bcca243dce78c9144914c26b44 100644
--- a/python/dune/perftool/pdelab/signatures.py
+++ b/python/dune/perftool/pdelab/signatures.py
@@ -141,10 +141,10 @@ def alpha_boundary_templates():
 
 def alpha_boundary_args():
     geo = name_geometry_wrapper()
-    lfsu = name_trialfunctionspace(Restriction.NEGATIVE)
-    lfsv = name_testfunctionspace(Restriction.NEGATIVE)
-    cc = name_coefficientcontainer(Restriction.NEGATIVE)
-    av = name_accumulation_variable((Restriction.NEGATIVE,))
+    lfsu = name_trialfunctionspace(Restriction.POSITIVE)
+    lfsv = name_testfunctionspace(Restriction.POSITIVE)
+    cc = name_coefficientcontainer(Restriction.POSITIVE)
+    av = name_accumulation_variable((Restriction.POSITIVE,))
     return ((True, geo), (True, lfsu), (True, cc), (True, lfsv), (False, av))
 
 
@@ -159,14 +159,14 @@ def alpha_skeleton_templates():
 
 def alpha_skeleton_args():
     geo = name_geometry_wrapper()
-    lfsu_s = name_trialfunctionspace(Restriction.NEGATIVE)
-    lfsu_n = name_trialfunctionspace(Restriction.POSITIVE)
-    lfsv_s = name_testfunctionspace(Restriction.NEGATIVE)
-    lfsv_n = name_testfunctionspace(Restriction.POSITIVE)
-    cc_s = name_coefficientcontainer(Restriction.NEGATIVE)
-    cc_n = name_coefficientcontainer(Restriction.POSITIVE)
-    av_s = name_accumulation_variable((Restriction.NEGATIVE,))
-    av_n = name_accumulation_variable((Restriction.POSITIVE,))
+    lfsu_s = name_trialfunctionspace(Restriction.POSITIVE)
+    lfsu_n = name_trialfunctionspace(Restriction.NEGATIVE)
+    lfsv_s = name_testfunctionspace(Restriction.POSITIVE)
+    lfsv_n = name_testfunctionspace(Restriction.NEGATIVE)
+    cc_s = name_coefficientcontainer(Restriction.POSITIVE)
+    cc_n = name_coefficientcontainer(Restriction.NEGATIVE)
+    av_s = name_accumulation_variable((Restriction.POSITIVE,))
+    av_n = name_accumulation_variable((Restriction.NEGATIVE,))
     return ((True, geo), (True, lfsu_s), (True, cc_s), (True, lfsv_s), (True, lfsu_n), (True, cc_n), (True, lfsv_n), (False, av_s), (False, av_n))
 
 
@@ -199,10 +199,10 @@ def jacobian_boundary_templates():
 
 def jacobian_boundary_args():
     geo = name_geometry_wrapper()
-    lfsu = name_trialfunctionspace(Restriction.NEGATIVE)
-    lfsv = name_testfunctionspace(Restriction.NEGATIVE)
-    cc = name_coefficientcontainer(Restriction.NEGATIVE)
-    av = name_accumulation_variable((Restriction.NEGATIVE, Restriction.NEGATIVE))
+    lfsu = name_trialfunctionspace(Restriction.POSITIVE)
+    lfsv = name_testfunctionspace(Restriction.POSITIVE)
+    cc = name_coefficientcontainer(Restriction.POSITIVE)
+    av = name_accumulation_variable((Restriction.POSITIVE, Restriction.POSITIVE))
     return ((True, geo), (True, lfsu), (True, cc), (True, lfsv), (False, av))
 
 
@@ -217,16 +217,16 @@ def jacobian_skeleton_templates():
 
 def jacobian_skeleton_args():
     geo = name_geometry_wrapper()
-    lfsu_s = name_trialfunctionspace(Restriction.NEGATIVE)
-    lfsu_n = name_trialfunctionspace(Restriction.POSITIVE)
-    lfsv_s = name_testfunctionspace(Restriction.NEGATIVE)
-    lfsv_n = name_testfunctionspace(Restriction.POSITIVE)
-    cc_s = name_coefficientcontainer(Restriction.NEGATIVE)
-    cc_n = name_coefficientcontainer(Restriction.POSITIVE)
-    av_ss = name_accumulation_variable((Restriction.NEGATIVE, Restriction.NEGATIVE))
-    av_sn = name_accumulation_variable((Restriction.NEGATIVE, Restriction.POSITIVE))
-    av_ns = name_accumulation_variable((Restriction.POSITIVE, Restriction.NEGATIVE))
-    av_nn = name_accumulation_variable((Restriction.POSITIVE, Restriction.POSITIVE))
+    lfsu_s = name_trialfunctionspace(Restriction.POSITIVE)
+    lfsu_n = name_trialfunctionspace(Restriction.NEGATIVE)
+    lfsv_s = name_testfunctionspace(Restriction.POSITIVE)
+    lfsv_n = name_testfunctionspace(Restriction.NEGATIVE)
+    cc_s = name_coefficientcontainer(Restriction.POSITIVE)
+    cc_n = name_coefficientcontainer(Restriction.NEGATIVE)
+    av_ss = name_accumulation_variable((Restriction.POSITIVE, Restriction.POSITIVE))
+    av_sn = name_accumulation_variable((Restriction.POSITIVE, Restriction.NEGATIVE))
+    av_ns = name_accumulation_variable((Restriction.NEGATIVE, Restriction.POSITIVE))
+    av_nn = name_accumulation_variable((Restriction.NEGATIVE, Restriction.NEGATIVE))
     return ((True, geo), (True, lfsu_s), (True, cc_s), (True, lfsv_s), (True, lfsu_n), (True, cc_n), (True, lfsv_n), (False, av_ss), (False, av_sn), (False, av_ns), (False, av_nn))
 
 
@@ -259,10 +259,10 @@ def jacobian_apply_boundary_templates():
 
 def jacobian_apply_boundary_args():
     geo = name_geometry_wrapper()
-    lfsu = name_trialfunctionspace(Restriction.NEGATIVE)
-    lfsv = name_testfunctionspace(Restriction.NEGATIVE)
-    ac = name_applycontainer(Restriction.NEGATIVE)
-    av = name_accumulation_variable((Restriction.NEGATIVE,))
+    lfsu = name_trialfunctionspace(Restriction.POSITIVE)
+    lfsv = name_testfunctionspace(Restriction.POSITIVE)
+    ac = name_applycontainer(Restriction.POSITIVE)
+    av = name_accumulation_variable((Restriction.POSITIVE,))
     return ((True, geo), (True, lfsu), (True, ac), (True, lfsv), (False, av))
 
 
@@ -277,14 +277,14 @@ def jacobian_apply_skeleton_templates():
 
 def jacobian_apply_skeleton_args():
     geo = name_geometry_wrapper()
-    lfsu_s = name_trialfunctionspace(Restriction.NEGATIVE)
-    lfsu_n = name_trialfunctionspace(Restriction.POSITIVE)
-    lfsv_s = name_testfunctionspace(Restriction.NEGATIVE)
-    lfsv_n = name_testfunctionspace(Restriction.POSITIVE)
-    ac_s = name_applycontainer(Restriction.NEGATIVE)
-    ac_n = name_applycontainer(Restriction.POSITIVE)
-    av_s = name_accumulation_variable((Restriction.NEGATIVE,))
-    av_n = name_accumulation_variable((Restriction.POSITIVE,))
+    lfsu_s = name_trialfunctionspace(Restriction.POSITIVE)
+    lfsu_n = name_trialfunctionspace(Restriction.NEGATIVE)
+    lfsv_s = name_testfunctionspace(Restriction.POSITIVE)
+    lfsv_n = name_testfunctionspace(Restriction.NEGATIVE)
+    ac_s = name_applycontainer(Restriction.POSITIVE)
+    ac_n = name_applycontainer(Restriction.NEGATIVE)
+    av_s = name_accumulation_variable((Restriction.POSITIVE,))
+    av_n = name_accumulation_variable((Restriction.NEGATIVE,))
     return ((True, geo), (True, lfsu_s), (True, ac_s), (True, lfsv_s), (True, lfsu_n), (True, ac_n), (True, lfsv_n), (False, av_s), (False, av_n))
 
 
@@ -318,11 +318,11 @@ def nonlinear_jacobian_apply_boundary_templates():
 
 def nonlinear_jacobian_apply_boundary_args():
     geo = name_geometry_wrapper()
-    lfsu = name_trialfunctionspace(Restriction.NEGATIVE)
-    lfsv = name_testfunctionspace(Restriction.NEGATIVE)
-    cc = name_coefficientcontainer(Restriction.NEGATIVE)
-    ac = name_applycontainer(Restriction.NEGATIVE)
-    av = name_accumulation_variable((Restriction.NEGATIVE,))
+    lfsu = name_trialfunctionspace(Restriction.POSITIVE)
+    lfsv = name_testfunctionspace(Restriction.POSITIVE)
+    cc = name_coefficientcontainer(Restriction.POSITIVE)
+    ac = name_applycontainer(Restriction.POSITIVE)
+    av = name_accumulation_variable((Restriction.POSITIVE,))
     return ((True, geo), (True, lfsu), (True, cc), (True, ac), (True, lfsv), (False, av))
 
 
@@ -337,14 +337,14 @@ def nonlinear_jacobian_apply_skeleton_templates():
 
 def nonlinear_jacobian_apply_skeleton_args():
     geo = name_geometry_wrapper()
-    lfsu_s = name_trialfunctionspace(Restriction.NEGATIVE)
-    lfsu_n = name_trialfunctionspace(Restriction.POSITIVE)
-    lfsv_s = name_testfunctionspace(Restriction.NEGATIVE)
-    lfsv_n = name_testfunctionspace(Restriction.POSITIVE)
-    cc_s = name_coefficientcontainer(Restriction.NEGATIVE)
-    cc_n = name_coefficientcontainer(Restriction.POSITIVE)
-    ac_s = name_applycontainer(Restriction.NEGATIVE)
-    ac_n = name_applycontainer(Restriction.POSITIVE)
-    av_s = name_accumulation_variable((Restriction.NEGATIVE,))
-    av_n = name_accumulation_variable((Restriction.POSITIVE,))
+    lfsu_s = name_trialfunctionspace(Restriction.POSITIVE)
+    lfsu_n = name_trialfunctionspace(Restriction.NEGATIVE)
+    lfsv_s = name_testfunctionspace(Restriction.POSITIVE)
+    lfsv_n = name_testfunctionspace(Restriction.NEGATIVE)
+    cc_s = name_coefficientcontainer(Restriction.POSITIVE)
+    cc_n = name_coefficientcontainer(Restriction.NEGATIVE)
+    ac_s = name_applycontainer(Restriction.POSITIVE)
+    ac_n = name_applycontainer(Restriction.NEGATIVE)
+    av_s = name_accumulation_variable((Restriction.POSITIVE,))
+    av_n = name_accumulation_variable((Restriction.NEGATIVE,))
     return ((True, geo), (True, lfsu_s), (True, cc_s), (True, ac_s), (True, lfsv_s), (True, lfsu_n), (True, cc_n), (True, ac_n), (True, lfsv_n), (False, av_s), (False, av_n))
diff --git a/python/dune/perftool/pdelab/spaces.py b/python/dune/perftool/pdelab/spaces.py
index 706be5b91a186b990719227fe61f1336d01313a0..d5f62735ea0b1d98c3496f49f0d062fc62bd1e5a 100644
--- a/python/dune/perftool/pdelab/spaces.py
+++ b/python/dune/perftool/pdelab/spaces.py
@@ -125,7 +125,7 @@ type_gfs = partial(_function_space_traversal, defaultname=available_gfs_names, r
 def initialize_function_spaces(expr, visitor):
     restriction = visitor.restriction
     if visitor.measure == 'exterior_facet':
-        restriction = Restriction.NEGATIVE
+        restriction = Restriction.POSITIVE
 
     index = None
     from ufl import MixedElement
diff --git a/python/dune/perftool/sumfact/accumulation.py b/python/dune/perftool/sumfact/accumulation.py
index 02eb474754f1233e8cfbe86f57c474db4d2ed9d2..e928b817eaf7b5f0eb9492f64880f554578c6fe0 100644
--- a/python/dune/perftool/sumfact/accumulation.py
+++ b/python/dune/perftool/sumfact/accumulation.py
@@ -269,7 +269,7 @@ def get_accumulation_info(expr, visitor):
     restriction = visitor.restriction
     if visitor.measure == 'exterior_facet':
         from dune.perftool.pdelab.restriction import Restriction
-        restriction = Restriction.NEGATIVE
+        restriction = Restriction.POSITIVE
 
     inames = visitor.interface.lfs_inames(leaf_element,
                                           restriction,
@@ -311,9 +311,9 @@ def _test_generator(expr, visitor):
     if visitor.measure == "cell":
         restrictions = (Restriction.NONE,)
     elif visitor.measure == "exterior_facet":
-        restrictions = (Restriction.NEGATIVE,)
+        restrictions = (Restriction.POSITIVE,)
     elif visitor.measure == "interior_facet":
-        restrictions = (Restriction.NEGATIVE, Restriction.POSITIVE)
+        restrictions = (Restriction.POSITIVE, Restriction.NEGATIVE)
     for res in restrictions:
         for ei, e in _get_childs(element):
             for grad in (None,) + tuple(range(dim)):
@@ -334,9 +334,9 @@ def _trial_generator(expr, visitor):
     if visitor.measure == "cell":
         restrictions = (Restriction.NONE,)
     elif visitor.measure == "exterior_facet":
-        restrictions = (Restriction.NEGATIVE,)
+        restrictions = (Restriction.POSITIVE,)
     elif visitor.measure == "interior_facet":
-        restrictions = (Restriction.NEGATIVE, Restriction.POSITIVE)
+        restrictions = (Restriction.POSITIVE, Restriction.NEGATIVE)
     for res in restrictions:
         for ei, e in _get_childs(element):
             yield SumfactAccumulationInfo(element_index=ei, restriction=res, element=e)
diff --git a/python/dune/perftool/sumfact/geometry.py b/python/dune/perftool/sumfact/geometry.py
index 0fd86eec08e519960d344cf13cf894f7243f287b..d0763db759836bcd49150c9a983dcd357a49006d 100644
--- a/python/dune/perftool/sumfact/geometry.py
+++ b/python/dune/perftool/sumfact/geometry.py
@@ -152,7 +152,7 @@ def pymbolic_spatial_coordinate_axiparallel(do_predicates, visitor):
     restriction = Restriction.NONE
     from dune.perftool.generation import get_global_context_value
     if get_global_context_value("integral_type") == "interior_facet":
-        restriction = Restriction.NEGATIVE
+        restriction = Restriction.POSITIVE
     from dune.perftool.sumfact.switch import get_facedir
     face = get_facedir(restriction)
 
@@ -183,8 +183,8 @@ def pymbolic_unit_outer_normal(visitor_indices):
     assert isinstance(index, int)
     if get_form_option("diagonal_transformation_matrix"):
         from dune.perftool.sumfact.switch import get_facedir, get_facemod
-        if index == get_facedir(Restriction.NEGATIVE):
-            if get_facemod(Restriction.NEGATIVE):
+        if index == get_facedir(Restriction.POSITIVE):
+            if get_facemod(Restriction.POSITIVE):
                 return 1, None
             else:
                 return -1, None
@@ -200,8 +200,8 @@ def pymbolic_unit_inner_normal(visitor_indices):
     assert isinstance(index, int)
     if get_form_option("diagonal_transformation_matrix"):
         from dune.perftool.sumfact.switch import get_facedir, get_facemod
-        if index == get_facedir(Restriction.NEGATIVE):
-            if get_facemod(Restriction.NEGATIVE):
+        if index == get_facedir(Restriction.POSITIVE):
+            if get_facemod(Restriction.POSITIVE):
                 return -1, None
             else:
                 return 1, None
@@ -221,7 +221,7 @@ def pymbolic_facet_jacobian_determinant():
 
 
 def pymbolic_constant_facet_jacobian_determinant():
-    facedir = get_facedir(Restriction.NEGATIVE)
+    facedir = get_facedir(Restriction.POSITIVE)
     assert isinstance(facedir, int)
 
     name = "fdetjac"
diff --git a/python/dune/perftool/sumfact/switch.py b/python/dune/perftool/sumfact/switch.py
index 873d7d462f7cd560d0c46a1ea811c34d9cec2ef5..8c6a0f13ace27b6e3081030607149b7b003aa82c 100644
--- a/python/dune/perftool/sumfact/switch.py
+++ b/python/dune/perftool/sumfact/switch.py
@@ -138,9 +138,9 @@ def generate_interior_facet_switch():
 
 def get_facedir(restriction):
     from dune.perftool.pdelab.restriction import Restriction
-    if restriction == Restriction.NEGATIVE or get_global_context_value("integral_type") == "exterior_facet":
+    if restriction == Restriction.POSITIVE or get_global_context_value("integral_type") == "exterior_facet":
         return get_global_context_value("facedir_s")
-    if restriction == Restriction.POSITIVE:
+    if restriction == Restriction.NEGATIVE:
         return get_global_context_value("facedir_n")
     if restriction == Restriction.NONE:
         return None
@@ -149,9 +149,9 @@ def get_facedir(restriction):
 
 def get_facemod(restriction):
     from dune.perftool.pdelab.restriction import Restriction
-    if restriction == Restriction.NEGATIVE or get_global_context_value("integral_type") == "exterior_facet":
+    if restriction == Restriction.POSITIVE or get_global_context_value("integral_type") == "exterior_facet":
         return get_global_context_value("facemod_s")
-    if restriction == Restriction.POSITIVE:
+    if restriction == Restriction.NEGATIVE:
         return get_global_context_value("facemod_n")
     if restriction == Restriction.NONE:
         return None
diff --git a/python/dune/perftool/ufl/modified_terminals.py b/python/dune/perftool/ufl/modified_terminals.py
index ac372d7fde8ba533fc52232f663a292026c2b7b4..bf3a6df939ddf787967b16c316cd2aaab308ce07 100644
--- a/python/dune/perftool/ufl/modified_terminals.py
+++ b/python/dune/perftool/ufl/modified_terminals.py
@@ -9,8 +9,8 @@ import ufl.classes as uc
 
 class Restriction:
     NONE = 0
-    NEGATIVE = 1
-    POSITIVE = 2
+    POSITIVE = 1
+    NEGATIVE = 2
 
 
 class ModifiedArgument(Record):
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index fb2d779316cc58eba19fd2ebf45c2ffa854aed01..3c4562c96dc1532461f46d53cd6ef8cc9edd7e8f 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -104,7 +104,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         # Correct the restriction on boundary integrals
         restriction = self.restriction
         if self.measure == 'exterior_facet':
-            restriction = Restriction.NEGATIVE
+            restriction = Restriction.POSITIVE
         leaf_element = o.ufl_element()
 
         # Select the correct leaf element in the case of this being a mixed finite element
@@ -130,7 +130,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         # Correct the restriction on boundary integrals
         restriction = self.restriction
         if self.measure == 'exterior_facet':
-            restriction = Restriction.NEGATIVE
+            restriction = Restriction.POSITIVE
 
         # Do something different for trial function and coefficients from jacobian apply
         if o.count() == 0 or o.count() == 1:
@@ -414,6 +414,11 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         # The normal must be restricted to be well-defined
         assert self.restriction is not Restriction.NONE
 
+        # Note: In UFL the jump is defined as: jump(v) = v('+') -
+        # v('-'). The corresponding outer unit normal is
+        # n=FacetNormal(cell)('+'). In order to be consisten with UFL
+        # we need to create the outer unit normal if the restriction
+        # is positive.
         if self.restriction == Restriction.POSITIVE:
             return self.interface.pymbolic_unit_outer_normal()
         if self.restriction == Restriction.NEGATIVE:
@@ -432,7 +437,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
     def jacobian_inverse(self, o):
         restriction = self.restriction
         if self.measure == 'exterior_facet':
-            restriction = Restriction.NEGATIVE
+            restriction = Restriction.POSITIVE
 
         assert(len(self.indices) == 2)
         i, j = self.indices
@@ -454,7 +459,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
     def cell_volume(self, o):
         restriction = self.restriction
         if self.measure == 'exterior_facet':
-            restriction = Restriction.NEGATIVE
+            restriction = Restriction.POSITIVE
 
         return self.interface.pymbolic_cell_volume(restriction)
 
diff --git a/test/heatequation/heatequation_dg.ufl b/test/heatequation/heatequation_dg.ufl
index 4f5c2da5f393ab31c8084405fe0d74daca9f2394..1d113a88f35b8404d7129a2fe5a2494f060e394b 100644
--- a/test/heatequation/heatequation_dg.ufl
+++ b/test/heatequation/heatequation_dg.ufl
@@ -1,12 +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)
@@ -14,21 +15,25 @@ 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
 
 poisson = inner(grad(u), grad(v))*dx \
-  + inner(n, avg(grad(u)))*jump(v)*dS \
-  + gamma*jump(u)*jump(v)*dS \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - f*v*dx \
+  - inner(n, avg(grad(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 \
+  - theta*g*inner(grad(v), n)*ds
 
 mass = (u*v)*dx
 
diff --git a/test/hyperbolic/linearacoustics.ufl b/test/hyperbolic/linearacoustics.ufl
index 844a3078bc6fb4994f9379a85aa8fa879c27df07..5a78e7848578053dbb8d7f75a94c2f901f5831d3 100644
--- a/test/hyperbolic/linearacoustics.ufl
+++ b/test/hyperbolic/linearacoustics.ufl
@@ -21,11 +21,11 @@ flux = as_matrix([[q0,  q1],
                   [0., rho]])
 
 # Define numerical fluxes to choose from
-llf_flux = dot(avg(flux), n) - 0.5*jump(u)
+llf_flux = dot(avg(flux), n) + 0.5*jump(u)
 numerical_flux = llf_flux
 
 r = -1. * inner(flux, grad(v))*dx \
-  - inner(numerical_flux, jump(v))*dS \
+  + inner(numerical_flux, jump(v))*dS \
   + inner(u, v)*ds
 
 interpolate_expression = f, 0.0, 0.0
diff --git a/test/hyperbolic/lineartransport.ufl b/test/hyperbolic/lineartransport.ufl
index 04497a40d8d85b24ec44d31d86a1328add3d2623..eaf8d33efe67a07a9a688efd3949e709a0a8122b 100644
--- a/test/hyperbolic/lineartransport.ufl
+++ b/test/hyperbolic/lineartransport.ufl
@@ -15,15 +15,15 @@ v = TestFunction(V)
 beta = as_vector((1., 1.))
 n = FacetNormal(cell)('+')
 
-def numerical_flux(normal, outside, inside):
+def numerical_flux(normal, inside, outside):
 	return conditional(inner(beta, n) > 0, inside, outside)*inner(beta, n)
 
 mass = u*v*dx
 
 r = -1.*u*inner(beta, grad(v))*dx \
-  - numerical_flux(n, u('+'), u('-'))*jump(v)*dS \
+  + numerical_flux(n, u('+'), u('-'))*jump(v)*dS \
   + inner(beta, n)*u*v*dso \
-  + numerical_flux(n, 0.0, u('-'))*v*dsd
+  + numerical_flux(n, u('+'), 0.0)*v*dsd
 
 is_dirichlet = dirichlet
 interpolate_expression = initial
diff --git a/test/hyperbolic/shallowwater.ufl b/test/hyperbolic/shallowwater.ufl
index 56267c99a725c7357bc923c2bdd844e81d141af9..c58429633b4ec11f5a69fc87fdbf0ed82c034a1b 100644
--- a/test/hyperbolic/shallowwater.ufl
+++ b/test/hyperbolic/shallowwater.ufl
@@ -23,13 +23,13 @@ b_flux = as_matrix([[-1.* q], [q*q/h + 0.5*g*h*h]])
 
 # Define numerical fluxes to choose from
 alpha = Max(abs(n[0]*q('+')) / h('+') + sqrt(g*h('+')), abs(n[0]*q('-')) / h('-') + sqrt(g*h('-')))
-llf_flux = dot(avg(flux), n) - 0.5*alpha*jump(u)
+llf_flux = dot(avg(flux), n) + 0.5*alpha*jump(u)
 alpha_b = abs(n[0]*q) / h + sqrt(g*h)
 boundary_flux = 0.5*dot(flux + b_flux, n) + alpha_b * as_vector([0., q])
 numerical_flux = llf_flux
 
 r = -1. * inner(flux, grad(v))*dx \
-  - inner(numerical_flux, jump(v))*dS \
+  + inner(numerical_flux, jump(v))*dS \
   + inner(boundary_flux, v)*ds
 
 interpolate_expression = f, 0.0
diff --git a/test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl b/test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl
index c1aaf70a3a0def29377613974547b8df88f57d3d..957d715e82acf50a936b9977957b3ae16c1d3e3b 100644
--- a/test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl
+++ b/test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl
@@ -37,11 +37,11 @@ r = mu * inner(grad(u), grad(v))*dx \
   - p*div(v)*dx \
   - q*div(u)*dx \
   + rho * inner(grad(u)*u,v)*dx \
-  + mu * inner(avg(grad(u))*n, jump(v))*dS \
+  - mu * inner(avg(grad(u))*n, jump(v))*dS \
   + mu * gamma_int * inner(jump(u), jump(v))*dS \
-  - mu * theta * inner(avg(grad(v))*n, jump(u))*dS \
-  - avg(p)*inner(jump(v), n)*dS \
-  - avg(q)*inner(jump(u), n)*dS \
+  + mu * theta * inner(avg(grad(v))*n, jump(u))*dS \
+  + avg(p)*inner(jump(v), n)*dS \
+  + avg(q)*inner(jump(u), n)*dS \
 
 interpolate_expression = g_v, g_p
 exact_solution = g_v, g_p
diff --git a/test/navier-stokes/navierstokes_3d_dg_quadrilateral.ufl b/test/navier-stokes/navierstokes_3d_dg_quadrilateral.ufl
index f6e8752f841256413fdd9dc5c7604e3543d3023c..96038a8b4ad4cc0b2ec8254f5a0bfcbcd4922373 100644
--- a/test/navier-stokes/navierstokes_3d_dg_quadrilateral.ufl
+++ b/test/navier-stokes/navierstokes_3d_dg_quadrilateral.ufl
@@ -2,7 +2,7 @@
 
 cell = hexahedron
 degree = 2
-dim = 2
+dim = 3
 
 x = SpatialCoordinate(cell)
 time = get_time(cell)
@@ -14,10 +14,6 @@ TH = P2 * P1
 v, q = TestFunctions(TH)
 u, p = TrialFunctions(TH)
 
-# g_v = as_vector((4.*x[1]*(1.-x[1]), 0.0, 0.0))
-# bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
-# ds = ds(subdomain_id=1, subdomain_data=bctype)
-
 n = FacetNormal(cell)('+')
 
 rho = 1.0
@@ -47,31 +43,16 @@ r = mu * inner(grad(u), grad(v))*dx \
   - p*div(v)*dx \
   - q*div(u)*dx \
   + rho * inner(grad(u)*u,v)*dx \
-  + mu * inner(avg(grad(u))*n, jump(v))*dS \
+  - mu * inner(avg(grad(u))*n, jump(v))*dS \
   + mu * gamma_int * inner(jump(u), jump(v))*dS \
-  - mu * theta * inner(avg(grad(v))*n, jump(u))*dS \
-  - avg(p)*inner(jump(v), n)*dS \
-  - avg(q)*inner(jump(u), n)*dS \
+  + mu * theta * inner(avg(grad(v))*n, jump(u))*dS \
+  + avg(p)*inner(jump(v), n)*dS \
+  + avg(q)*inner(jump(u), n)*dS \
   - mu * inner(grad(u)*n, v)*ds \
   + mu * gamma_ext * inner(u-g_v, v)*ds \
   + mu * theta * inner(grad(v)*n, u-g_v)*ds \
   + p*inner(v, n)*ds \
   + q*inner(u-g_v, n)*ds
 
-# r = mu * inner(grad(u), grad(v))*dx \
-#   - p*div(v)*dx \
-#   - q*div(u)*dx \
-#   + rho * inner(grad(u)*u,v)*dx \
-#   + mu * inner(avg(grad(u))*n, jump(v))*dS \
-#   + mu / h_e * sigma * inner(jump(u), jump(v))*dS \
-#   - mu * theta * inner(avg(grad(v))*n, jump(u))*dS \
-#   - avg(p)*inner(jump(v), n)*dS \
-#   - avg(q)*inner(jump(u), n)*dS \
-#   - mu * inner(grad(u)*n, v)*ds \
-#   + mu / h_e * sigma * inner(u-g_v, v)*ds \
-#   + mu * theta * inner(grad(v)*n, u-g_v)*ds \
-#   + p*inner(v, n)*ds \
-#   + q*inner(u-g_v, n)*ds
-
 interpolate_expression = g_v, g_p
 exact_solution = g_v, g_p
\ No newline at end of file
diff --git a/test/nonlinear/diffusivewave.ufl b/test/nonlinear/diffusivewave.ufl
index fc92f13d4efd9dcaa0826ae47892148b6e6bd234..5dddfe3cb719343a82fd74f9a109460a97ecb71e 100644
--- a/test/nonlinear/diffusivewave.ufl
+++ b/test/nonlinear/diffusivewave.ufl
@@ -1,8 +1,10 @@
 cell = quadrilateral
+degree = 1
+dim = 2
 
 x = SpatialCoordinate(cell)
 
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
 
 u = TrialFunction(V)
 v = TestFunction(V)
@@ -10,7 +12,11 @@ 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
@@ -20,14 +26,14 @@ K = u**(5./3.)
 # / Max(1e-8, norm)
 
 poisson = inner(K*grad(u), grad(v))*dx \
-  + inner(n, avg(K*grad(u)))*jump(v)*dS \
-  + gamma*jump(u)*jump(v)*dS \
-  - theta*jump(u)*inner(avg(K*grad(v)), n)*dS
+  - inner(n, avg(K*grad(u)))*jump(v)*dS \
+  + gamma_int*jump(u)*jump(v)*dS \
+  + theta*jump(u)*inner(avg(K*grad(v)), n)*dS
 #  - inner(n, K*grad(u))*v*ds \
-#  + gamma*u*v*ds \
+#  + gamma_ext*u*v*ds \
 #  + theta*u*inner(K*grad(v), n)*ds \
-#  - theta*g*inner(K*grad(v), n)*ds \
-#  - gamma*g*v*ds
+#  - gamma_ext*g*v*ds \
+#  - theta*g*inner(K*grad(v), n)*ds
 
 mass = (u*v)*dx
 
diff --git a/test/nonlinear/nonlinear_dg.ufl b/test/nonlinear/nonlinear_dg.ufl
index 47b850759e24461a394c230b112a7e5929ae1ea4..41dfb430c1048e011c04843f496b35947695877f 100644
--- a/test/nonlinear/nonlinear_dg.ufl
+++ b/test/nonlinear/nonlinear_dg.ufl
@@ -1,10 +1,12 @@
 cell = "triangle"
-x = SpatialCoordinate(cell)
+degree = 1
+dim = 2
 
+x = SpatialCoordinate(cell)
 f = -4.
 g = x[0]*x[0] + x[1]*x[1]
 
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
 
 u = TrialFunction(V)
 v = TestFunction(V)
@@ -15,21 +17,25 @@ def q(u):
     return u*u
 
 # 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 \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
-  - inner(n, grad(u))*v*ds \
-  + gamma*u*v*ds \
-  + theta*u*inner(grad(v), n)*ds \
   - f*v*dx \
   + q(u)*v*dx \
-  - theta*g*inner(grad(v), n)*ds \
-  - gamma*g*v*ds
+  - inner(n, avg(grad(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_ext*u*v*ds \
+  + theta*u*inner(grad(v), n)*ds \
+  - gamma_ext*g*v*ds \
+  - theta*g*inner(grad(v), n)*ds
 
 exact_solution = g
diff --git a/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl b/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl
index 15a60efc606078c91fcd91b77c3b21180d73808e..51776b3acb57f87519bb2d493a1e833034c29b3d 100644
--- a/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_1d_dg_interval.ufl
@@ -1,11 +1,13 @@
 cell = "interval"
+degree = 1
+dim = 1
 
 x = SpatialCoordinate(cell)
 c = (0.5-x[0])**2
 g = exp(-1.*c)
 f = 2*(1.-2*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 \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - f*v*dx \
+  - inner(n, avg(grad(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 \
+  - theta*g*inner(grad(v), n)*ds
 
 exact_solution = g
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl b/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl
index 5b4cf2d71b09181e307ce9dbe835a0d030e05935..87fd9b591b7d6ce406b655d52926d7117218a90d 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.ufl
@@ -1,11 +1,13 @@
 cell = "quadrilateral"
+degree = 1
+dim = 2
 
 x = SpatialCoordinate(cell)
 c = (0.5-x[0])**2 + (0.5-x[1])**2
 g = exp(-1.*c)
 f = 2*(2.-2*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 \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - f*v*dx \
+  - inner(n, avg(grad(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 \
+  - theta*g*inner(grad(v), n)*ds
 
 exact_solution = g
diff --git a/test/poisson/dimension-grid-variations/poisson_2d_dg_triangle.ufl b/test/poisson/dimension-grid-variations/poisson_2d_dg_triangle.ufl
index 2121ff1b2a635ec4128d521d143f0456c046911a..250b6900d28fa305534334c37d6e07563ad36aa5 100644
--- a/test/poisson/dimension-grid-variations/poisson_2d_dg_triangle.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_2d_dg_triangle.ufl
@@ -1,10 +1,12 @@
 cell = triangle
+degree = 1
+dim = 2
 
 x = SpatialCoordinate(cell)
 f = -4.
 g = x[0]*x[0] + x[1]*x[1]
 
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
 
 u = TrialFunction(V)
 v = TestFunction(V)
@@ -12,20 +14,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 \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - f*v*dx \
+  - inner(n, avg(grad(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 \
+  - theta*g*inner(grad(v), n)*ds
 
 exact_solution = g
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl
index 3925b02e2f18c7e477aa6109dd4b56fd4ce9b740..810bfa46234b7a3951173cb36e9825a3dd7338c0 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_dg_hexahedron.ufl
@@ -1,11 +1,13 @@
 cell = "hexahedron"
+degree = 1
+dim = 3
 
 x = SpatialCoordinate(cell)
 c = (0.5 - x[0])**2 + (0.5 - x[1])**2 + (0.5 - x[2])**2
 g = exp(-1.*c)
 f = 2*(3.-2*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 \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - f*v*dx \
+  - inner(n, avg(grad(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 \
+  - theta*g*inner(grad(v), n)*ds
 
-exact_solution = g
+exact_solution = g
\ No newline at end of file
diff --git a/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl b/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl
index 5cc027cf1fc4ddf0eca4109982c15cb9e340c19c..f3af6ada266ea8b3099e37ebc544688bf6bfca18 100644
--- a/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl
+++ b/test/poisson/dimension-grid-variations/poisson_3d_dg_tetrahedron.ufl
@@ -1,10 +1,12 @@
 cell = "tetrahedron"
+degree = 1
+dim = 3
 
 x = SpatialCoordinate(cell)
 f = -6.
 g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
 
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
 
 u = TrialFunction(V)
 v = TestFunction(V)
@@ -12,20 +14,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 \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - f*v*dx \
+  - inner(n, avg(grad(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 \
+  - theta*g*inner(grad(v), n)*ds
 
-exact_solution = g
+exact_solution = g
\ No newline at end of file
diff --git a/test/poisson/opcount_poisson_dg.ufl b/test/poisson/opcount_poisson_dg.ufl
index 5db12e14b8a85c7c417534269c18d52e7a3e3596..4962ce0cefce2cd5ed901994eb8a3b8b75f187a9 100644
--- a/test/poisson/opcount_poisson_dg.ufl
+++ b/test/poisson/opcount_poisson_dg.ufl
@@ -1,9 +1,9 @@
-degree = 1
 cell = "quadrilateral"
 
+degree = 1
 dim = 2
-x = SpatialCoordinate(cell)
 
+x = SpatialCoordinate(cell)
 f = -4.
 g = x[0]*x[0] + x[1]*x[1]
 
@@ -25,14 +25,14 @@ 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 \
+  - f*v*dx \
+  - inner(n, avg(grad(u)))*jump(v)*dS \
   + gamma_int*jump(u)*jump(v)*dS \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  + theta*jump(u)*inner(avg(grad(v)), n)*dS \
   - inner(n, grad(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_ext*g*v*ds
+  - gamma_ext*g*v*ds \
+  - theta*g*inner(grad(v), n)*ds
 
-exact_solution = g
+exact_solution = g
\ No newline at end of file
diff --git a/test/poisson/poisson_dg.ufl b/test/poisson/poisson_dg.ufl
index fa263e98c6732ba53ca58ab4fd7daa8422d86f36..75b536b851aa6752ea3c7505087f61f7f85177bd 100644
--- a/test/poisson/poisson_dg.ufl
+++ b/test/poisson/poisson_dg.ufl
@@ -25,14 +25,14 @@ 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 \
+  - f*v*dx \
+  - inner(n, avg(grad(u)))*jump(v)*dS \
   + gamma_int*jump(u)*jump(v)*dS \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  + theta*jump(u)*inner(avg(grad(v)), n)*dS \
   - inner(n, grad(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_ext*g*v*ds
+  - gamma_ext*g*v*ds \
+  - theta*g*inner(grad(v), n)*ds
 
 exact_solution = g
\ No newline at end of file
diff --git a/test/poisson/poisson_dg_neumann.ufl b/test/poisson/poisson_dg_neumann.ufl
index 0e3a899544f9c2d5941e6d954d2a85e3f65f51ac..bdcea093ef8d0d0b53f4a965513d5d8ae1b058e2 100644
--- a/test/poisson/poisson_dg_neumann.ufl
+++ b/test/poisson/poisson_dg_neumann.ufl
@@ -1,6 +1,6 @@
-dim = 3
-degree = 1
 cell = triangle
+degree = 1
+dim = 2
 
 x = SpatialCoordinate(cell)
 c = (0.5-x[0])**2 + (0.5-x[1])**2
@@ -30,15 +30,15 @@ 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 \
+  - f*v*dx \
+  - inner(n, avg(grad(u)))*jump(v)*dS \
   + gamma_int*jump(u)*jump(v)*dS \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  + theta*jump(u)*inner(avg(grad(v)), n)*dS \
   - inner(n, grad(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_ext*g*v*ds(1) \
+  - theta*g*inner(grad(v), n)*ds(1) \
   - j*v*ds(0)
 
 exact_solution = g
diff --git a/test/poisson/poisson_dg_quadrilateral.ufl b/test/poisson/poisson_dg_quadrilateral.ufl
index c8e1107df230b5af44cc2262c5e5c20f262ada3f..7387894341f25221d412ae967717a3bfca8788f4 100644
--- a/test/poisson/poisson_dg_quadrilateral.ufl
+++ b/test/poisson/poisson_dg_quadrilateral.ufl
@@ -1,11 +1,13 @@
 cell = "quadrilateral"
+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 \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - f*v*dx \
+  - inner(n, avg(grad(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 \
+  - theta*g*inner(grad(v), n)*ds
 
 exact_solution = g
diff --git a/test/poisson/poisson_dg_tensor.ufl b/test/poisson/poisson_dg_tensor.ufl
index 437e2f708bdcc4a6b455920c0864c20ccd406a8c..ba8eb74d084f7e4b3efc9e940beb616e305b609c 100644
--- a/test/poisson/poisson_dg_tensor.ufl
+++ b/test/poisson/poisson_dg_tensor.ufl
@@ -1,6 +1,6 @@
-dim = 2
-degree = 1
 cell = quadrilateral
+degree = 1
+dim = 2
 
 x = SpatialCoordinate(cell)
 
@@ -28,14 +28,14 @@ gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
 theta = 1.0
 
 r = (inner(A*grad(u), grad(v)) + c*u*v)*dx \
-  + inner(n, A*avg(grad(u)))*jump(v)*dS \
+  - f*v*dx \
+  - inner(n, A*avg(grad(u)))*jump(v)*dS \
   + gamma_int*jump(u)*jump(v)*dS \
-  - theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
+  + theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
   - inner(n, A*grad(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_ext*g*v*ds
+  - gamma_ext*g*v*ds \
+  - theta*g*inner(A*grad(v), n)*ds
 
 exact_solution = g
diff --git a/test/stokes/stokes_3d_dg_quadrilateral.ufl b/test/stokes/stokes_3d_dg_quadrilateral.ufl
index 3fefa9b763714ca27b7fa7d741bfe94059dae836..c47773f5821ce09ecb6182e81eeaff6bf17a7fc7 100644
--- a/test/stokes/stokes_3d_dg_quadrilateral.ufl
+++ b/test/stokes/stokes_3d_dg_quadrilateral.ufl
@@ -1,11 +1,13 @@
 cell = hexahedron
+degree = 2
+dim = 3
 
 x = SpatialCoordinate(cell)
 g_v = as_vector((4.*x[1]*(1.-x[1]), 0.0, 0.0))
 bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
 
-P2 = VectorElement("DG", cell, 2)
-P1 = FiniteElement("DG", cell, 1)
+P2 = VectorElement("DG", cell, degree)
+P1 = FiniteElement("DG", cell, degree-1)
 TH = P2 * P1
 
 v, q = TestFunctions(TH)
@@ -14,21 +16,28 @@ u, p = TrialFunctions(TH)
 ds = ds(subdomain_id=1, subdomain_data=bctype)
 
 n = FacetNormal(cell)('+')
-eps = -1.0
-sigma = 1.0
-h_e = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+
+# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
+theta = -1.0
+
+# penalty factor
+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
 
 r = inner(grad(u), grad(v))*dx \
   - p*div(v)*dx \
   - q*div(u)*dx \
-  + inner(avg(grad(u))*n, jump(v))*dS \
-  + sigma / h_e * inner(jump(u), jump(v))*dS \
-  - eps * inner(avg(grad(v))*n, jump(u))*dS \
-  - avg(p)*inner(jump(v), n)*dS \
-  - avg(q)*inner(jump(u), n)*dS \
+  - inner(avg(grad(u))*n, jump(v))*dS \
+  + gamma_int * inner(jump(u), jump(v))*dS \
+  + theta * inner(avg(grad(v))*n, jump(u))*dS \
+  + avg(p)*inner(jump(v), n)*dS \
+  + avg(q)*inner(jump(u), n)*dS \
   - inner(grad(u)*n, v)*ds \
-  + sigma / h_e * inner(u-g_v, v)*ds \
-  + eps * inner(grad(v)*n, u-g_v)*ds \
+  + gamma_ext * inner(u-g_v, v)*ds \
+  + theta * inner(grad(v)*n, u-g_v)*ds \
   + p*inner(v, n)*ds \
   + q*inner(u-g_v, n)*ds
 
diff --git a/test/stokes/stokes_dg.ufl b/test/stokes/stokes_dg.ufl
index 22d239c8166311e74b9dbb489f7b4499a5e9aaa3..d4f4225b939d228fc4436d482a3cc989a97ad78f 100644
--- a/test/stokes/stokes_dg.ufl
+++ b/test/stokes/stokes_dg.ufl
@@ -1,11 +1,13 @@
 cell = triangle
+degree = 2
+dim = 2
 
 x = SpatialCoordinate(cell)
 g_v = as_vector((4*x[1]*(1.-x[1]), 0.0))
 bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
 
-P2 = VectorElement("DG", cell, 2)
-P1 = FiniteElement("DG", cell, 1)
+P2 = VectorElement("DG", cell, degree)
+P1 = FiniteElement("DG", cell, degree-1)
 TH = P2 * P1
 
 v, q = TestFunctions(TH)
@@ -14,21 +16,28 @@ u, p = TrialFunctions(TH)
 ds = ds(subdomain_id=1, subdomain_data=bctype)
 
 n = FacetNormal(cell)('+')
-eps = -1.0
-sigma = 1.0
-h_e = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+
+# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
+theta = -1.0
+
+# penalty factor
+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
 
 r = inner(grad(u), grad(v))*dx \
   - p*div(v)*dx \
   - q*div(u)*dx \
-  + inner(avg(grad(u))*n, jump(v))*dS \
-  + sigma / h_e * inner(jump(u), jump(v))*dS \
-  - eps * inner(avg(grad(v))*n, jump(u))*dS \
-  - avg(p)*inner(jump(v), n)*dS \
-  - avg(q)*inner(jump(u), n)*dS \
+  - inner(avg(grad(u))*n, jump(v))*dS \
+  + gamma_int * inner(jump(u), jump(v))*dS \
+  + theta * inner(avg(grad(v))*n, jump(u))*dS \
+  + avg(p)*inner(jump(v), n)*dS \
+  + avg(q)*inner(jump(u), n)*dS \
   - inner(grad(u)*n, v)*ds \
-  + sigma / h_e * inner(u-g_v, v)*ds \
-  + eps * inner(grad(v)*n, u-g_v)*ds \
+  + gamma_ext * inner(u-g_v, v)*ds \
+  + theta * inner(grad(v)*n, u-g_v)*ds \
   + p*inner(v, n)*ds \
   + q*inner(u-g_v, n)*ds
 
diff --git a/test/stokes/stokes_dg_quadrilateral.ufl b/test/stokes/stokes_dg_quadrilateral.ufl
index 3342fd9a67cc1961817cc9537e50623b012f6ac0..0b37429b54cbe4e77fad3eb0abd88d9e345b7dcf 100644
--- a/test/stokes/stokes_dg_quadrilateral.ufl
+++ b/test/stokes/stokes_dg_quadrilateral.ufl
@@ -1,11 +1,13 @@
 cell = quadrilateral
+degree = 2
+dim = 2
 
 x = SpatialCoordinate(cell)
 g_v = as_vector((4*x[1]*(1.-x[1]), 0.0))
 bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
 
-P2 = VectorElement("DG", cell, 2)
-P1 = FiniteElement("DG", cell, 1)
+P2 = VectorElement("DG", cell, degree)
+P1 = FiniteElement("DG", cell, degree-1)
 TH = P2 * P1
 
 v, q = TestFunctions(TH)
@@ -14,21 +16,28 @@ u, p = TrialFunctions(TH)
 ds = ds(subdomain_id=1, subdomain_data=bctype)
 
 n = FacetNormal(cell)('+')
-eps = -1.0
-sigma = 1.0
-h_e = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+
+# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
+theta = -1.0
+
+# penalty factor
+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
 
 r = inner(grad(u), grad(v))*dx \
   - p*div(v)*dx \
   - q*div(u)*dx \
-  + inner(avg(grad(u))*n, jump(v))*dS \
-  + sigma / h_e * inner(jump(u), jump(v))*dS \
-  - eps * inner(avg(grad(v))*n, jump(u))*dS \
-  - avg(p)*inner(jump(v), n)*dS \
-  - avg(q)*inner(jump(u), n)*dS \
+  - inner(avg(grad(u))*n, jump(v))*dS \
+  + gamma_int * inner(jump(u), jump(v))*dS \
+  + theta * inner(avg(grad(v))*n, jump(u))*dS \
+  + avg(p)*inner(jump(v), n)*dS \
+  + avg(q)*inner(jump(u), n)*dS \
   - inner(grad(u)*n, v)*ds \
-  + sigma / h_e * inner(u-g_v, v)*ds \
-  + eps * inner(grad(v)*n, u-g_v)*ds \
+  + gamma_ext * inner(u-g_v, v)*ds \
+  + theta * inner(grad(v)*n, u-g_v)*ds \
   + p*inner(v, n)*ds \
   + q*inner(u-g_v, n)*ds
 
diff --git a/test/sumfact/hyperbolic/linearacoustics.ufl b/test/sumfact/hyperbolic/linearacoustics.ufl
index 844a3078bc6fb4994f9379a85aa8fa879c27df07..5a78e7848578053dbb8d7f75a94c2f901f5831d3 100644
--- a/test/sumfact/hyperbolic/linearacoustics.ufl
+++ b/test/sumfact/hyperbolic/linearacoustics.ufl
@@ -21,11 +21,11 @@ flux = as_matrix([[q0,  q1],
                   [0., rho]])
 
 # Define numerical fluxes to choose from
-llf_flux = dot(avg(flux), n) - 0.5*jump(u)
+llf_flux = dot(avg(flux), n) + 0.5*jump(u)
 numerical_flux = llf_flux
 
 r = -1. * inner(flux, grad(v))*dx \
-  - inner(numerical_flux, jump(v))*dS \
+  + inner(numerical_flux, jump(v))*dS \
   + inner(u, v)*ds
 
 interpolate_expression = f, 0.0, 0.0
diff --git a/test/sumfact/hyperbolic/lineartransport.ufl b/test/sumfact/hyperbolic/lineartransport.ufl
index 0afaef57f46bb9c8a8f9d56fe05c9695c4ff6a8f..15a8021854222931b8dbef707082e90b406a47dc 100644
--- a/test/sumfact/hyperbolic/lineartransport.ufl
+++ b/test/sumfact/hyperbolic/lineartransport.ufl
@@ -15,15 +15,15 @@ v = TestFunction(V)
 beta = as_vector((1., 1.))
 n = FacetNormal(cell)('+')
 
-def numerical_flux(normal, outside, inside):
+def numerical_flux(normal, inside, outside):
 	return conditional(inner(beta, n) > 0, inside, outside)*inner(beta, n)
 
 mass = u*v*dx
 
 r = -1.*u*inner(beta, grad(v))*dx \
-  - numerical_flux(n, u('+'), u('-'))*jump(v)*dS \
+  + numerical_flux(n, u('+'), u('-'))*jump(v)*dS \
   + inner(beta, n)*u*v*dso \
-  + numerical_flux(n, 0.0, u('-'))*v*dsd
+  + numerical_flux(n, u('+'), 0.0)*v*dsd
 
 exact_solution = 0.0
 is_dirichlet = dirichlet
diff --git a/test/sumfact/hyperbolic/shallowwater.ufl b/test/sumfact/hyperbolic/shallowwater.ufl
index 0f5f24443a63d041163ead2ad70a8401a099d5b8..d0bfc147cf83062855836606edaf4a337f5d21c4 100644
--- a/test/sumfact/hyperbolic/shallowwater.ufl
+++ b/test/sumfact/hyperbolic/shallowwater.ufl
@@ -27,12 +27,12 @@ bflux = as_matrix([[-q0,                  -q1],
 
 
 # Define numerical fluxes to choose from
-llf_flux = dot(avg(flux), n) - 0.5*jump(u)
+llf_flux = dot(avg(flux), n) + 0.5*jump(u)
 boundary_flux = 0.5*dot(flux + bflux, n) + as_vector([0., q0, q1])
 numerical_flux = llf_flux
 
 r = -1. * inner(flux, grad(v))*dx \
-  - inner(numerical_flux, jump(v))*dS \
+  + inner(numerical_flux, jump(v))*dS \
   + inner(boundary_flux, v)*ds
 
 interpolate_expression = f, 0.0, 0.0
diff --git a/test/sumfact/poisson/poisson_dg_2d.ufl b/test/sumfact/poisson/poisson_dg_2d.ufl
index 1da5733c19e48a96e931948c2ed5e46214c068d7..3c2cf8767cdc395643366936529fe44c2a39be5b 100644
--- a/test/sumfact/poisson/poisson_dg_2d.ufl
+++ b/test/sumfact/poisson/poisson_dg_2d.ufl
@@ -1,4 +1,5 @@
 cell = "quadrilateral"
+dim = 2
 
 x = SpatialCoordinate(cell)
 c = (0.5-x[0])**2 + (0.5-x[1])**2
@@ -13,20 +14,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 \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - f*v*dx \
+  - inner(n, avg(grad(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 \
+  - theta*g*inner(grad(v), n)*ds
 
 exact_solution = g
diff --git a/test/sumfact/poisson/poisson_dg_3d.ufl b/test/sumfact/poisson/poisson_dg_3d.ufl
index 20705f2f58a0ba1cdb66dfdcf01f1f0835612ebf..0f7f0399b8e07ad85e45a20071304e469b1b6133 100644
--- a/test/sumfact/poisson/poisson_dg_3d.ufl
+++ b/test/sumfact/poisson/poisson_dg_3d.ufl
@@ -1,4 +1,5 @@
-cell = "hexahedron"
+cell = hexahedron
+dim = 3
 
 x = SpatialCoordinate(cell)
 c = (0.5-x[0])**2 + (0.5-x[1])**2 + (0.5-x[2])**2
@@ -13,20 +14,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 \
-  - theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - f*v*dx \
+  - inner(n, avg(grad(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 \
+  - theta*g*inner(grad(v), n)*ds
 
 exact_solution = g
diff --git a/test/sumfact/poisson/poisson_dg_tensor.ufl b/test/sumfact/poisson/poisson_dg_tensor.ufl
index 8388ed75b40f2108ba56076e974404267d2c7310..2734383f71d3afddc39c6c72f41ba5545def3b54 100644
--- a/test/sumfact/poisson/poisson_dg_tensor.ufl
+++ b/test/sumfact/poisson/poisson_dg_tensor.ufl
@@ -1,14 +1,14 @@
 cell = hexahedron
+dim = 3
 
 x = SpatialCoordinate(cell)
-
 I = Identity(3)
 A = as_matrix([[x[i]*x[j] + I[i,j] for j in range(3)] for i in range(3)])
 g = x[0]**2 + x[1]**2 + x[2]**2
 c = 10.
 f = -6.
 
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
 
 u = TrialFunction(V)
 v = TestFunction(V)
@@ -16,20 +16,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 \
-  - theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
-  - inner(n, A*grad(u))*v*ds \
-  + gamma*u*v*ds \
-  + theta*u*inner(A*grad(v), n)*ds \
+r = inner(grad(u), grad(v))*dx \
   - f*v*dx \
-  - theta*g*inner(A*grad(v), n)*ds \
-  - gamma*g*v*ds
+  - inner(n, avg(grad(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_ext*u*v*ds \
+  + theta*u*inner(grad(v), n)*ds \
+  - gamma_ext*g*v*ds \
+  - theta*g*inner(grad(v), n)*ds
 
 exact_solution = g
diff --git a/test/sumfact/stokes/stokes_3d_dg.ufl b/test/sumfact/stokes/stokes_3d_dg.ufl
index 3fefa9b763714ca27b7fa7d741bfe94059dae836..8193a8933e9a3f16722737802e93199b2739ec30 100644
--- a/test/sumfact/stokes/stokes_3d_dg.ufl
+++ b/test/sumfact/stokes/stokes_3d_dg.ufl
@@ -21,11 +21,11 @@ h_e = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
 r = inner(grad(u), grad(v))*dx \
   - p*div(v)*dx \
   - q*div(u)*dx \
-  + inner(avg(grad(u))*n, jump(v))*dS \
+  - inner(avg(grad(u))*n, jump(v))*dS \
   + sigma / h_e * inner(jump(u), jump(v))*dS \
-  - eps * inner(avg(grad(v))*n, jump(u))*dS \
-  - avg(p)*inner(jump(v), n)*dS \
-  - avg(q)*inner(jump(u), n)*dS \
+  + eps * inner(avg(grad(v))*n, jump(u))*dS \
+  + avg(p)*inner(jump(v), n)*dS \
+  + avg(q)*inner(jump(u), n)*dS \
   - inner(grad(u)*n, v)*ds \
   + sigma / h_e * inner(u-g_v, v)*ds \
   + eps * inner(grad(v)*n, u-g_v)*ds \
diff --git a/test/sumfact/stokes/stokes_dg.ufl b/test/sumfact/stokes/stokes_dg.ufl
index ddefd870a413fe7099150e619d336c6d089da5d6..7f873537ee0c6b5f1015091abef27d660823b331 100644
--- a/test/sumfact/stokes/stokes_dg.ufl
+++ b/test/sumfact/stokes/stokes_dg.ufl
@@ -21,11 +21,11 @@ h_e = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
 r = inner(grad(u), grad(v))*dx \
   - p*div(v)*dx \
   - q*div(u)*dx \
-  + inner(avg(grad(u))*n, jump(v))*dS \
+  - inner(avg(grad(u))*n, jump(v))*dS \
   + sigma / h_e * inner(jump(u), jump(v))*dS \
-  - eps * inner(avg(grad(v))*n, jump(u))*dS \
-  - avg(p)*inner(jump(v), n)*dS \
-  - avg(q)*inner(jump(u), n)*dS \
+  + eps * inner(avg(grad(v))*n, jump(u))*dS \
+  + avg(p)*inner(jump(v), n)*dS \
+  + avg(q)*inner(jump(u), n)*dS \
   - inner(grad(u)*n, v)*ds \
   + sigma / h_e * inner(u-g_v, v)*ds \
   + eps * inner(grad(v)*n, u-g_v)*ds \