From 20853255fbcfa71ed2a9637a2975f251d750c9c7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Thu, 22 Feb 2018 13:46:21 +0100
Subject: [PATCH] Adapt DG tests (except for hyperbolic stuff) to new
 self/neighbor

I also made the order of the terms more consistent and it is now much
easier to compare the skeletons and the boundaries.
---
 test/heatequation/heatequation_dg.ufl         | 25 +++++++++-------
 .../navierstokes_2d_dg_quadrilateral.ufl      |  8 ++---
 .../navierstokes_3d_dg_quadrilateral.ufl      | 29 ++++---------------
 test/nonlinear/diffusivewave.ufl              | 22 +++++++++-----
 test/nonlinear/nonlinear_dg.ufl               | 28 +++++++++++-------
 .../poisson_1d_dg_interval.ufl                | 24 +++++++++------
 .../poisson_2d_dg_quadrilateral.ufl           | 24 +++++++++------
 .../poisson_2d_dg_triangle.ufl                | 24 +++++++++------
 .../poisson_3d_dg_hexahedron.ufl              | 26 ++++++++++-------
 .../poisson_3d_dg_tetrahedron.ufl             | 26 ++++++++++-------
 test/poisson/opcount_poisson_dg.ufl           | 16 +++++-----
 test/poisson/poisson_dg.ufl                   | 10 +++----
 test/poisson/poisson_dg_neumann.ufl           | 12 ++++----
 test/poisson/poisson_dg_quadrilateral.ufl     | 24 +++++++++------
 test/poisson/poisson_dg_tensor.ufl            | 14 ++++-----
 test/stokes/stokes_3d_dg_quadrilateral.ufl    |  8 ++---
 test/stokes/stokes_dg.ufl                     |  8 ++---
 test/stokes/stokes_dg_quadrilateral.ufl       |  8 ++---
 test/sumfact/poisson/poisson_dg_2d.ufl        | 21 +++++++++-----
 test/sumfact/poisson/poisson_dg_3d.ufl        | 23 +++++++++------
 test/sumfact/poisson/poisson_dg_tensor.ufl    | 28 ++++++++++--------
 test/sumfact/stokes/stokes_3d_dg.ufl          |  8 ++---
 test/sumfact/stokes/stokes_dg.ufl             |  8 ++---
 23 files changed, 236 insertions(+), 188 deletions(-)

diff --git a/test/heatequation/heatequation_dg.ufl b/test/heatequation/heatequation_dg.ufl
index 4f5c2da5..1d113a88 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/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl b/test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl
index c1aaf70a..957d715e 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 f6e8752f..96038a8b 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 fc92f13d..5dddfe3c 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 47b85075..41dfb430 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 15a60efc..51776b3a 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 5b4cf2d7..87fd9b59 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 2121ff1b..250b6900 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 3925b02e..810bfa46 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 5cc027cf..f3af6ada 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 5db12e14..4962ce0c 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 fa263e98..75b536b8 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 0e3a8995..bdcea093 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 c8e1107d..73878943 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 437e2f70..ba8eb74d 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 3fefa9b7..8193a893 100644
--- a/test/stokes/stokes_3d_dg_quadrilateral.ufl
+++ b/test/stokes/stokes_3d_dg_quadrilateral.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/stokes/stokes_dg.ufl b/test/stokes/stokes_dg.ufl
index 22d239c8..2398c6f6 100644
--- a/test/stokes/stokes_dg.ufl
+++ b/test/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 \
diff --git a/test/stokes/stokes_dg_quadrilateral.ufl b/test/stokes/stokes_dg_quadrilateral.ufl
index 3342fd9a..0062259b 100644
--- a/test/stokes/stokes_dg_quadrilateral.ufl
+++ b/test/stokes/stokes_dg_quadrilateral.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/poisson/poisson_dg_2d.ufl b/test/sumfact/poisson/poisson_dg_2d.ufl
index 1da5733c..3c2cf876 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 20705f2f..0f7f0399 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 8388ed75..2734383f 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 3fefa9b7..8193a893 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 ddefd870..7f873537 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 \
-- 
GitLab