new mesh + bubble inside

testcases/xmesh/bubble.py

@@ -33,7 +33,7 @@ class Poiseuille(unittest.TestCase) :
         if not os.path.isdir(outputdir) :
             os.makedirs(outputdir)
 
-        subprocess.call(["gmsh", "-2", "mesh.geo","-clscale", "2"])
+        subprocess.call(["gmsh", "-2", "mesh2.geo","-clscale", "2"])
         t = 0
         ii = 0
 
@@ -50,29 +50,53 @@ class Poiseuille(unittest.TestCase) :
                 return Uin
             else:
                 return 0
+        
+        def setInitialConcentration(r):
+            concentration = np.full(fluid.n_nodes(),0)
+            #tags, x, _ = gmsh.model.mesh.get_nodes(includeBoundary=False)
+            x = fluid.coordinates()
+            print(x.shape)
+            #x = x.reshape([-1,3])
+            # circle centered at (0,0.15), radius r
+            for i in range(fluid.n_nodes()):
+                if (x[i,0]-0)**2 + (x[i,1]-0.15)**2 < r**2 :
+                    concentration[i] = 1
+            fluid.set_concentration_cg(concentration)
+        
 
         #numerical parameters
         lcmin = .1                                  # mesh size
-        dt = 0.0001                                    # time step
+        dt = 0.001                                    # time step
         alpha = 1e-4                                    # stabilization coefficient
 
-        shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+        shutil.copy("mesh2.msh", outputdir +"/mesh2.msh")
 
         outf = 1                                       # number of iterations between output files
 
         #Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
         fluid = mbfluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1])
-        fluid.load_msh("mesh.msh")
+        fluid.load_msh("mesh2.msh")
 
-        fluid.set_open_boundary("BottomMiddle", velocity=[0, boundaryV], concentration=1)
+        #fluid.set_open_boundary("BottomMiddle", velocity=[0, boundaryV], concentration=1)
+        #fluid.set_wall_boundary("BottomLeft",velocity=[0,0])
+        #fluid.set_wall_boundary("BottomRight",velocity=[0,0])
+        
+        fluid.set_wall_boundary("BottomMiddle", velocity=[0,0])
         fluid.set_wall_boundary("BottomLeft",velocity=[0,0])
         fluid.set_wall_boundary("BottomRight",velocity=[0,0])
+        
+        #fluid.set_open_boundary("BottomMiddle", velocity=[0, boundaryV], concentration=1)
+        #fluid.set_open_boundary("BottomLeft",velocity=[0, boundaryV], concentration=1)
+        #fluid.set_open_boundary("BottomRight",velocity=[0, boundaryV], concentration=1)
+        
         # fluid.set_open_boundary("Insert", velocity=[0, boundaryV], concentration=1)
         fluid.set_open_boundary("Top",pressure=0, concentration=0)
         fluid.set_wall_boundary("Right",velocity=[0,0])
         fluid.set_wall_boundary("Left",velocity=[0,0])
-
-
+        
+        setInitialConcentration(0.04)
+        
+        
         ii = 0
         t = 0
 
@@ -80,22 +104,24 @@ class Poiseuille(unittest.TestCase) :
 
         fluid.write_vtk(outputdir,0,0)
 
-        mesh,x = xmesh2d.Mesh.import_from_file("mesh.msh")
+        mesh,x = xmesh2d.Mesh.import_from_file("mesh2.msh")
         tracker = xmesh2d.FrontTracker(mesh,x)
 
         ii = 0
         tic = time.time()
-        while ii <  1000:
+        while ii <  200:
             #Fluid solver
             time_integration.iterate(fluid,None,dt)
             t += dt
 
             concentration = fluid.concentration_cg().reshape((fluid.n_nodes(),))
+            print(concentration)
             level = concentration - 0.5;
+            #level = np.where(concentration < 0.5, -1, 1)
 
             newx = np.copy(fluid.coordinates()[:])
             newx = tracker.move_front(newx, level[:,])
-            newx = tracker.relax(newx, dt*5)
+            newx = tracker.relax(newx, 1)
             print(x[:,:2].shape)
             print(np.sum(newx-fluid.coordinates()))
             fluid.mesh_velocity()[:] = (newx[:,:2] - fluid.coordinates()[:,:2])/dt

testcases/xmesh/mesh2.geo

+L = 0.3;
+H = 0.4;
+y = 0;
+lc = .01;
+
+Point(1) = {-L, H, 0, lc};
+Point(2) = {-L, 0, 0, lc};
+Point(3) = {L, 0, 0, lc};
+Point(4) = {L, H, 0, lc};
+Point(5) = {-L/10, 0, 0, lc};
+Point(6) = {L/10, 0, 0, lc};
+Line(1) = {1, 2};
+Line(2) = {2, 5};
+Line(3) = {5, 6};
+Line(4) = {6, 3};
+Line(5) = {3, 4};
+Line(6) = {4, 1};
+
+Line Loop(1) = {1:6};
+
+Plane Surface(1) = {1};
+
+Physical Line("Left") = {1};
+Physical Line("BottomLeft") = {2};
+Physical Line("BottomMiddle") = {3};
+Physical Line("BottomRight") = {4};
+Physical Line("Right") = {5};
+Physical Line("Top") = {6};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+Mesh.MshFileVersion = 2;

testcases/xmesh/xmesh2d.py

@@ -117,7 +117,7 @@ class Mesh():
         model = gmsh.model.get_current()
         self.tritags, tri = gmsh.model.mesh.get_elements_by_type(2)
         self.tri = tri.reshape([-1,3])
-        tags, x, _ = gmsh.model.mesh.get_nodes()
+        tags, x, _ = gmsh.model.mesh.get_nodes(includeBoundary=False)
         order = np.argsort(tags)
         self.nodetags = tags[order]
         x = x.reshape([-1,3])[order]
@@ -147,6 +147,7 @@ class FrontTracker():
         self.previous_sign = np.full(x.shape[0], 1)
 
     def move_front(self, x, levelset):
+        print("in xmesh")
 
         def move_node(tn) :
             all_neighbours = neigh[neigh_start[tn]:neigh_start[tn+1]]
@@ -170,6 +171,7 @@ class FrontTracker():
                 return False
 
         neigh = self.mesh.neigh
+        boundary_nodes = self.mesh.boundary_nodes
         neigh_start = self.mesh.neigh_start
         front = self.front
         sign = np.where(levelset > 0, 1, -1)
@@ -181,6 +183,7 @@ class FrontTracker():
             nsign = sign[all_neighbours]
             nsign[front[all_neighbours]] = 0
             nsign = np.r_[sign[tn],nsign]
+            print(nsign)
             if  nsign.max() - nsign.min() != 2 :
                 front[tn] = False
         # move the nodes already in the front