Update Shear

fluid/fluid_problem.c

@@ -193,7 +193,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
     dp[iD] = ds[P*D+iD];
     unold += sold[U+iD]*n[iD];
     unmesh += mesh_velocity[iD]*n[iD];
-    uext[iD] = (vid<0?s[U+iD]:data[vid+iD]);
+    uext[iD] = (vid<0?s[U+iD]:data[vid+iD]*c);
     unext += uext[iD]*n[iD];
     dcn += dc[iD]*n[iD];
     s_c += (u[iD]-uext[iD])*n[iD];

testcases/shear-2d/mesh.geo

@@ -4,10 +4,10 @@ H = a;
 y = 0;
 lc = 0.1*a;
 
-Point(1) = {0,0,0,lc};
-Point(2) = {L,0,0,lc};
-Point(3) = {L,H,0,lc};
-Point(4) = {0,H,0,lc};
+Point(1) = {-L/2,-H/2,0,lc};
+Point(2) = {L/2,-H/2,0,lc};
+Point(3) = {L/2,H/2,0,lc};
+Point(4) = {-L/2,H/2,0,lc};
 Line(1) = {1,2};
 Line(2) = {2,3};
 Line(3) = {4,3};

testcases/shear-2d/shear.py

@@ -29,6 +29,7 @@ import time
 import shutil
 import random
 
+np.random.seed(0)
 dir_path = os.path.dirname(os.path.realpath(__file__))
 os.chdir(dir_path)
 outputdir = "outputShear2D_2"
@@ -36,6 +37,8 @@ if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
 
 subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
+subprocess.call(["gmsh", "-2", "meshFluid.geo","-clscale","1"])
+
 t = 0
 ii = 0
 
@@ -59,13 +62,15 @@ dt = 1e-3                                                # time step
 # dt = 10                                             # time step
 
 shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+shutil.copy("meshFluid.msh", outputdir +"/mesh.msh")
+
 outf = 20                                           # 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,mu,rho,petsc_solver_type="-pc_type lu")
-fluid.load_msh("mesh.msh")
-fluid.set_wall_boundary("Bottom",velocity=[0,0], pressure = 0)
-fluid.set_wall_boundary("Top",velocity=[1.0,0])
+fluid.load_msh("meshFluid.msh")
+fluid.set_wall_boundary("Bottom", velocity=[0,0], pressure = 0)
+fluid.set_wall_boundary("Top", velocity=[1.0,0])
 time_integration.iterate(fluid,None,10)
 
 # if strong boundary on periodic line, it should be forced on both sides
@@ -76,9 +81,10 @@ time_integration.iterate(fluid,None,10)
 # Particle Problem
 p = scontact.ParticleProblem(2)
 p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Left", "Right"])
-for x in np.arange(r,L-r,2.1*r):
-    for y in np.arange(r,H-r,2.1*r):
+for x in np.arange(-L/2+r,L/2-r,2.1*r):
+    for y in np.arange(-H/2+r,H/2-r,2.1*r):
         R = r*(1-np.random.random()/4)
+        # R = r
         p.add_particle((x+R,y),R,R**2*np.pi*rhop)
 
 
@@ -88,24 +94,40 @@ t = 0
 
 #set initial_condition
 fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
-fluid.export_vtk(outputdir,0,0)
+fluid.write_vtk(outputdir,0,0)
 p.write_vtk(outputdir,0,0)
 
 ii = 0
 tic = time.time()
+
+xTop = np.array([0.0, H/2])
+vTop = np.array([0.0, 0.0])
+k = 5.0
+Fext = lambda x: np.array([0.0, k*(H/2 - x[1])])
 # while t < tEnd :
 while t < tEnd:
     #Fluid solver
     time_integration.iterate(fluid,p,dt,min_nsub=5)
 
+    fTop = p.get_boundary_forces("Top")
+    mTop = np.mean(p.mass())
+    print(f"F Top normal = {fTop[0]}")
+    print(f"F ext = {Fext(xTop)}")
+
+    vTop += (fTop[0] + Fext(xTop))*dt/mTop
+    xTop += vTop*dt
+    p.set_boundary_velocity("Top", vTop)
+
     #Ensure particle teleportation
-    p.position()[:,0] = np.fmod(p.position()[:,0],L)
+    p.position()[:,0] = np.fmod(p.position()[:,0] + L/2,L) - L/2
+    ind = np.where(p.position()[:,0] < -L/2)
+    p.position()[ind,0] += L
 
     t += dt
     #Output files writting
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
-        fluid.export_vtk(outputdir, t, ioutput)
+        fluid.write_vtk(outputdir, ioutput, t)
         p.write_vtk(outputdir, ioutput, t)
     ii += 1
     print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
\ No newline at end of file