diff --git a/benchmark.py b/benchmark.py
index eba7b61a88..64543607dc 100644
--- a/benchmark.py
+++ b/benchmark.py
@@ -173,7 +173,7 @@ def _elevate_degree(mesh, degree):
     return Mesh(f)
 
 
-use_netgen = True
+use_netgen = False
 quadrilateral = True
 
 T = 20 # 10.0 # 12.0
@@ -183,7 +183,7 @@ def _elevate_degree(mesh, degree):
 ntimesteps = int(T / dt_float)
 t = Constant(0.0)
 dim = 2
-degree = 2  # 2 - 4
+degree = 3  # 2 - 4
 if use_netgen:
     nref = 1 #  # 2 - 5 tested for CSM 1 and 2
     mesh  = make_mesh_netgen(0.1 / 2 ** nref)
@@ -683,10 +683,10 @@ def compute_elast_tensors(dim, u, lambda_s, mu_s):
         inner(dot(- p_mid * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
         #inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
     else:  # CN
-        #theta_p = Constant(1. / 2. + 100 * float(dt))
-        #theta_m = Constant(1. / 2. - 100 * float(dt))
-        theta_p = Constant(1.)
-        theta_m = Constant(0.)
+        theta_p = Constant(1. / 2. + 100 * float(dt))
+        theta_m = Constant(1. / 2. - 100 * float(dt))
+        #theta_p = Constant(1.)
+        #theta_m = Constant(0.)
         v_f_dot = (v_f - v_f_0) / dt
         u_f_dot = (u_f - u_f_0) / dt
         v_s_dot = (v_s - v_s_0) / dt
@@ -695,8 +695,8 @@ def _fluid(v_f, u_f, p):
             F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
             return (inner(rho_f * J_f * v_f_dot, dv_f) +
                     inner(rho_f * J_f * dot(dot(grad(v_f), inv(F_f)), v_f - u_f_dot), dv_f) +
-                    inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f), inv(F_f))), dot(grad(dv_f), inv(F_f))) +
-                    J_f * inner(p, tr(dot(grad(dv_f), inv(F_f)))) -
+                    inner(rho_f * J_f * nu_f * 2 * sym(dot(grad(v_f), inv(F_f))), dot(grad(dv_f), inv(F_f))) -
+                    J_f * inner(p, tr(dot(grad(dv_f), inv(F_f)))) +
                     J_f * inner(tr(dot(grad(v_f), inv(F_f))), dp) +
                     J_f * inner(dot(grad(u_f), inv(F_f)), dot(grad(du_f), inv(F_f)))) * dx_f
         def _struct(v_f, u_f, p, v_s, u_s):
@@ -707,17 +707,17 @@ def _struct(v_f, u_f, p, v_s, u_s):
                     inner(rho_s * J_s * as_vector([0., - g_s]), dv_s) +
                     inner(J_s * (u_s_dot - v_s), du_s)) * dx_s + \
                    inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
-        #residual_f = theta_p * _fluid(v_f, u_f, p) + \
-        #             theta_m * _fluid(v_f_0, u_f_0, p_0)
-        #residual_s = theta_p * _struct(v_f, u_f, p, v_s, u_s) + \
-        #             theta_m * _struct(v_f_0, u_f_0, p_0, v_s_0, u_s_0)
+        residual_f = theta_p * _fluid(v_f, u_f, p) + \
+                     theta_m * _fluid(v_f_0, u_f_0, p_0)
+        residual_s = theta_p * _struct(v_f, u_f, p, v_s, u_s) + \
+                     theta_m * _struct(v_f_0, u_f_0, p_0, v_s_0, u_s_0)
+        """
         F_f, J_f, E_f, S_f = compute_elast_tensors(dim, u_f, lambda_s, mu_s)
         F_s, J_s, E_s, S_s = compute_elast_tensors(dim, u_s, lambda_s, mu_s)
         v_f_mid = v_f
         v_s_mid = v_s
         u_f_mid = u_f
         p_mid = p
-        """
         residual_f = (
             inner(rho_f * J_f * (v_f - v_f_0) / dt, dv_f) +
             inner(rho_f * J_f * dot(dot(grad(v_f_mid), inv(F_f)), v_f_mid - (u_f - u_f_0) / dt), dv_f) +
@@ -735,8 +735,6 @@ def _struct(v_f, u_f, p, v_s, u_s):
         inner(dot(- p('|') * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f('|')), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
         #inner(dot(- p_mid * Identity(dim) + rho_f * nu_f * 2 * sym(dot(grad(v_f_mid), inv(F_f))), dot(J_f * transpose(inv(F_f)), n_f)), dv_s('|')) * ds_s(label_interface)
         """
-        residual_f = _fluid(v_f, u_f, p)
-        residual_s = _struct(v_f, u_f, p, v_s, u_s)
     residual = residual_f + residual_s
     def v_f_left(t_):
         return 1.5 * Ubar * y_f * (H - y_f) / ((H / 2) ** 2) * conditional(t_ < 2.0 + dt / 10., (1 - cos(pi / 2 * t_)) / 2., 1.)
diff --git a/mesh_f.pdf b/mesh_f.pdf
index cadef13fc9..a08635e8b2 100644
Binary files a/mesh_f.pdf and b/mesh_f.pdf differ
diff --git a/mesh_orig.pdf b/mesh_orig.pdf
index 805b768bd0..3238d469c0 100644
Binary files a/mesh_orig.pdf and b/mesh_orig.pdf differ
diff --git a/mesh_s.pdf b/mesh_s.pdf
index af3b90378b..5fd87cdfee 100644
Binary files a/mesh_s.pdf and b/mesh_s.pdf differ
diff --git a/time_series.pdf b/time_series.pdf
index 7eba2d278d..2b41b76c1d 100644
Binary files a/time_series.pdf and b/time_series.pdf differ