benchmark

benchmark.py

@@ -123,8 +123,8 @@ def make_mesh_quad():
     yarray = np.concatenate([np.linspace(y0, y1, int(np.ceil((y1 - y0) / h)), endpoint=False),
                              np.linspace(y1, y2, int(np.ceil((y2 - y1) / h)), endpoint=False),
                              np.array([y2, 0.19, 0.21]),
-                             np.linspace(y3, y4, int(np.ceil((x4 - x3) / h)), endpoint=False),
-                             np.linspace(y4, y5, int(np.ceil((x5 - x4) / h)) + 1, endpoint=True)])
+                             np.linspace(y3, y4, int(np.ceil((y4 - y3) / h)), endpoint=False),
+                             np.linspace(y4, y5, int(np.ceil((y5 - y4) / h)) + 1, endpoint=True)])
     mesh = TensorRectangleMesh(xarray, yarray, quadrilateral=True)
     x, y = SpatialCoordinate(mesh)
     V = FunctionSpace(mesh, "DQ", 0)
@@ -169,7 +169,8 @@ def _elevate_degree(mesh, degree, label):
     f.dat.data_with_halos[bc.nodes, 1] = (f.dat.data_with_halos[bc.nodes, 1] - Cy) * (r / r_) + Cy
     return Mesh(f)
 
-
+
+quad = True 
 T = 20 # 10.0 # 12.0
 dt_float = 0.0025  #.002
 dt = Constant(dt_float)  #0.001
@@ -180,7 +181,7 @@ def _elevate_degree(mesh, degree, label):
 nref = 4 #  # 2 - 5 tested for CSM 1 and 2
 h = 0.1 / 2 ** nref
 degree = 3 #3  # 2 - 4
-if True:
+if not quad:
     mesh, label_left, label_right, label_bottom, label_top, label_interface_orig, label_interface, label_struct_base, label_circle = make_mesh(h)
     mesh = _elevate_degree(mesh, degree, label_circle + label_struct_base)
 else:
@@ -204,7 +205,10 @@ def _elevate_degree(mesh, degree, label):
 x_f, y_f = SpatialCoordinate(mesh_f)
 n_f = FacetNormal(mesh_f)
 mesh_s = Submesh(mesh, dmcommon.CELL_SETS_LABEL, label_struct, mesh.topological_dimension())
-mesh_s = _elevate_degree(mesh_s, degree, label_struct_base)
+if quad:
+    mesh_s = _elevate_degree(mesh_s, degree, label_circle)
+else:
+    mesh_s = _elevate_degree(mesh_s, degree, label_struct_base)
 x_s, y_s = SpatialCoordinate(mesh_s)
 n_s = FacetNormal(mesh_s)
 dx = Measure("dx", domain=mesh)
@@ -241,6 +245,7 @@ def _elevate_degree(mesh, degree, label):
     triplot(mesh_s, axes=axes)
     axes.legend()
     plt.savefig('mesh_s.pdf')
+    raise RuntimeError("not error")
 
 
 case = "FSI3_2"