adding test for file in loadmsh

python/scontact.py

@@ -562,6 +562,9 @@ class ParticleProblem :
         shift -- Optional argument to shift the numerical domain
         material -- Material of the boundary
         """
+        if not os.path.isfile(filename):
+          print("Error : no such file as " + filename)
+          exit(1)
         if shift is None :
             shift = [0]*self._dim
         shift = np.array(shift)
@@ -663,4 +666,4 @@ class ParticleProblem :
             
     def dim(self) :
         """Returns the dimension of the particle problem"""
-        return self._dim
\ No newline at end of file
+        return self._dim

testcases/drag-2d/drag.py

@@ -23,6 +23,7 @@
 
 # TESTCASE DESCRIPTION
 # Study of the drag on a fixed disk in a 2D immersed granular flow
+# Requires a pre-computed deposit
 from migflow import scontact
 from migflow import fluid
 from migflow import time_integration
@@ -42,9 +43,9 @@ dt = 1e-3                                      #time step
 t = 0                                          #initial time
 tEnd = 10                                      #final time
 ii = 0                                         #iteration number
-outf = 1000                                      #iterations between data frames  
+outf = 1000                                    #iterations between data frames  
 # Define output directory
-outputdir = "notOkTrial/" + vitesse[index] 
+outputdir = "output/" + vitesse[index] 
 if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
 filename = outputdir + "/Drag.csv"
@@ -53,14 +54,14 @@ filename = outputdir + "/Drag.csv"
 #
 #Injected particles
 p = scontact.ParticleProblem(2,True,True)
-p.read_vtk("depot/outputFluid",10)
+p.read_vtk("depot/output",10)
 p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.61))
 #Flow particles
 p2 = scontact.ParticleProblem(2,True,True)
-p2.read_vtk("depot/outputFluid",10)
+p2.read_vtk("depot/output",10)
 p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
 p2.clear_boundaries()
-p2.load_msh_boundaries("meshP.msh", ["Inner","Left","Right"],material = "Steel") 
+p2.load_msh_boundaries("mesh.msh", ["Inner","Left","Right"],material = "Steel") 
 p2.set_friction_coefficient(0.3,"Sand","Sand")
 p2.set_friction_coefficient(0.3,"Sand","Steel")
 p2.write_vtk(outputdir, 0, 0)
@@ -79,22 +80,18 @@ p2.write_vtk(outputdir, 0, 0)
 # 
 # COMPUTATION LOOP
 #
-def accumulate(F,nsub):
-  F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)/nsub
+def accumulate(F,n_divide):
+  F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)/n_divide
 while t < tEnd : 
-     print(vit)
      p2.forced_flag()[(p2.position()[:,1] + p2.r()[:,0] <= -0.5)] = 1
-     if t > 2:
-       p2.velocity()[p2.forced_flag()==1,:] = [0, vit]
-     else :
-       p2.velocity()[p2.forced_flag()==1,:] = [0, 0]
+     p2.velocity()[p2.forced_flag()==1,:] = [0, vit]
      p2.omega()[p2.forced_flag()==0] = 1
      if  max(p2.position()[:,1]+p2.r()[:,0]) < 0.5 and t < 9 :
        p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
      p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.55))
      F = np.zeros(2)
      #time_integration.particle_changed(fluid,p2)
-     time_integration.iterate(None,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda nsub : accumulate(F,nsub),max_nsub=4)
+     time_integration.iterate(None,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda n_divide : accumulate(F,n_divide))
      with open(filename,"a") as file1:
         #F += fluid.boundary_force()[0,:]
         file1.write( str(F[0])+";"+str(F[1])+";"+str(t)+"\n")

testcases/drag-2d/mesh.geo

 rin = 0.01905;
-wout = 0.56;
+wout = 0.7;
 lout = 0.15;
 lcin = wout/2; 
 lcout = wout/2;

testcases/drag-3d/drag.py

@@ -88,8 +88,8 @@ p2.write_vtk(outputdir, 0, 0)
 # 
 # COMPUTATION LOOP
 #
-def accumulate(F,nsub):
-  F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)/nsub
+def accumulate(F,n_divide):
+  F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)/n_divide
 while t < tEnd : 
      p2.forced_flag()[(p2.position()[:,1] + p2.r()[:,0] <= -0.5)] = 1
      p2.velocity()[p2.forced_flag()==1,:] = [0, vit, 0]
@@ -99,7 +99,7 @@ while t < tEnd :
      p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.52))
      F = np.zeros(3)
      #time_integration.particle_changed(fluid,p2)
-     time_integration.iterate(None,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda nsub : accumulate(F,nsub))
+     time_integration.iterate(None,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda n_divide : accumulate(F,n_divide))
      with open(filename,"a") as file1:
         #F += fluid.boundary_force()[0,:]
         file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+"\n")