benchmark

Author: ksagiyam

benchmark.py

@@ -357,20 +357,29 @@ def _elevate_degree(mesh, degree):
     E_s = mu_s * 2 * (1 + nu_s)
     lambda_s = nu_s * E_s / (1 + nu_s) / (1 - 2 * nu_s)
     V_0 = VectorFunctionSpace(mesh_f, "P", degree)
-    V_1 = FunctionSpace(mesh_f, "P", degree - 2)
-    V = V_0 * V_1
+    V_1 = FunctionSpace(mesh_f, "P", degree - 1)
+    V = V_0 * V_0 * V_1
     solution = Function(V)
     solution_0 = Function(V)
-    v_f, p_f = split(solution)
-    u_f = Function(V_0)
-    v_f_0, p_f_0 = split(solution_0)
-    dv_f, dp_f = split(TestFunction(V))
+    v_f, u_f, p_f = split(solution)
+    v_f_0, u_f_0, p_f_0 = split(solution_0)
+    dv_f, du_f, dp_f = split(TestFunction(V))
+    F = Identity(dim) + grad(u_f)
+    J = det(F)
+    E = 1. / 2. * (dot(transpose(F), F) - Identity(dim))
+    S = lambda_s * tr(E) * Identity(dim) + 2.0 * mu_s * E
+    #residual = (
+    #    rho_f * J * inner(v_f - v_f_0, dv_f) / dt +
+    #    rho_f * J * inner(dot(grad(v_f), dot(inv(F), v_f - (u_f - u_f_0) / 2)), dv_f) +
+    #    rho_f * J * nu_f * inner(dot(2 * sym(dot(grad(v_f), inv(F))), transpose(inv(F))), grad(dv_f)) -
+    #    J * inner(p_f * transpose(inv(F)), grad(dv_f)) +
+    #    J * inner(tr(dot(grad(v_f), inv(F))), dp_f)) * dx
     residual = (
-                rho_f * inner(v_f - v_f_0, dv_f) / dt +
-                rho_f * inner(dot(grad(v_f), v_f), dv_f) +
-                rho_f * nu_f * inner(2 * sym(grad(v_f)), grad(dv_f)) -
-                inner(p_f, div(dv_f)) +
-                inner(div(v_f), dp_f)) * dx
+        rho_f * J * inner(v_f - v_f_0, dv_f) / dt +
+        rho_f * J * inner(dot(grad(v_f), dot(inv(F), v_f - (u_f - u_f_0) / 2)), dv_f) +
+        rho_f * J * nu_f * inner(2 * sym(dot(grad(v_f), inv(F))), dot(grad(dv_f), inv(F))) -
+        J * inner(p_f, tr(dot(grad(dv_f), inv(F)))) +
+        J * inner(tr(dot(grad(v_f), inv(F))), dp_f) + J * inner(u_f, du_f)) * dx
     v_f_left = 1.5 * Ubar * y_f * (H - y_f) / ((H / 2) ** 2) * conditional(t < 2.0, (1 - cos(pi / 2 * t)) / 2., 1.)
     bc_inflow = DirichletBC(V.sub(0), as_vector([v_f_left, 0.]), label_left)
     bc_noslip = DirichletBC(V.sub(0), Constant((0, 0)), label_bottom + label_top + label_circle + label_interface)