From b46a974d9b135e24c54ead798a93c7a315677495 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Mon, 5 Mar 2018 16:30:27 +0100
Subject: [PATCH] Adjust penalty term in stokes dg examples
---
applications/stokes_dg/stokes_dg.ufl | 28 ++++++++++++++--------
test/stokes/stokes_3d_dg_quadrilateral.ufl | 27 ++++++++++++++-------
test/stokes/stokes_dg.ufl | 27 ++++++++++++++-------
test/stokes/stokes_dg_quadrilateral.ufl | 27 ++++++++++++++-------
4 files changed, 72 insertions(+), 37 deletions(-)
diff --git a/applications/stokes_dg/stokes_dg.ufl b/applications/stokes_dg/stokes_dg.ufl
index 6189ae1b..d36556db 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/test/stokes/stokes_3d_dg_quadrilateral.ufl b/test/stokes/stokes_3d_dg_quadrilateral.ufl
index 8193a893..c47773f5 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 \
+ + 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 2398c6f6..d4f4225b 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 \
+ + 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 0062259b..0b37429b 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 \
+ + 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
--
GitLab