@@ -177,7 +177,7 @@ def _elevate_degree(mesh, degree, label):
177177ntimesteps = int (T / dt_float )
178178t = Constant (0.0 )
179179dim = 2
180- nref = 0 # # 2 - 5 tested for CSM 1 and 2
180+ nref = 4 # # 2 - 5 tested for CSM 1 and 2
181181h = 0.1 / 2 ** nref
182182degree = 3 #3 # 2 - 4
183183if True :
@@ -187,9 +187,9 @@ def _elevate_degree(mesh, degree, label):
187187 # quad mesh
188188 mesh , label_left , label_right , label_bottom , label_top , label_interface_orig , label_interface , label_struct_base , label_circle = make_mesh_quad ()
189189 mesh = _elevate_degree (mesh , degree , label_circle + label_struct_base )
190+ x , y = SpatialCoordinate (mesh )
190191"""
191192#mesh.topology_dm.viewFromOptions("-dm_view")
192- x, y = SpatialCoordinate(mesh)
193193v = assemble(Constant(1.0, domain=mesh) * ds(label_circle + label_struct_base))
194194print(v - 2 * pi * r)
195195print(assemble(x * dx(label_struct)))
@@ -217,6 +217,11 @@ def _elevate_degree(mesh, degree, label):
217217dS_f = Measure ("dS" , domain = mesh_f )
218218dS_s = Measure ("dS" , domain = mesh_s )
219219
220+ #print(repr(mesh))
221+ #print(repr(mesh_f))
222+ #print(repr(mesh_s), flush=True)
223+
224+
220225if mesh .comm .size == 1 :
221226 fig , axes = plt .subplots ()
222227 axes .axis ('equal' )
@@ -496,7 +501,8 @@ def _elevate_degree(mesh, degree, label):
496501 J * inner ((u - u_0 ) / dt - v , du )
497502 ) * dx_s # dx(label_struct)
498503 residual = residual_f + residual_s
499- 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. )
504+ #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.)
505+ v_f_left = 1.5 * Ubar * y * (H - y ) / ((H / 2 ) ** 2 ) * conditional (t < 2.0 , (1 - cos (pi / 2 * t )) / 2. , 1. )
500506 bc_inflow = DirichletBC (V .sub (0 ), as_vector ([v_f_left , 0. ]), label_left )
501507 bc_noslip_v = DirichletBC (V .sub (0 ), Constant ((0 , 0 )), label_bottom + label_top + label_circle + label_struct_base )
502508 bc_noslip_u = DirichletBC (V .sub (1 ), Constant ((0 , 0 )), label_left + label_right + label_bottom + label_top + label_circle + label_struct_base )
0 commit comments