writing contact forces instead of impulses

python/scontact.py

@@ -273,6 +273,8 @@ class ParticleProblem :
         return self._lib.particleProblemIterate(self._p, c_double(np.min(self.r()) if self.r().size != 0 else tol), c_double(dt), c_double(tol), c_int(-1),c_int(force_motion),_np2c(forces))
 
     def get_contacts(self,contact_type) :
+        """Gives the contact forces. Attention : during the resolution fo the contacts, 
+           the considered quantities are impulses, while forces are written and read."""
         ctype = {"particle_particle":0,"particle_disk":1,"particle_segment":2,"particle_triangle":3}[contact_type]
         n = self._lib.particleProblemContactN(self._p,c_int(ctype))
         basis = np.ndarray((n,self._dim*self._dim),dtype=np.float64,order="C")

scontact/scontact.c

@@ -1121,6 +1121,13 @@ int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol
 
 int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit,int force_motion, double *externalForces)
 {
+  for (size_t j = 0; j<vectorSize(p->contacts); ++j){
+    p->contacts[j].dv *= dt;
+    p->contacts[j].ct *= dt;
+#if DIMENSION==3
+    p->contacts[j].cs *= dt;
+#endif
+  }
   double *oldVelocity = malloc(sizeof(double)*vectorSize(p->velocity));
   for (size_t j = 0; j < vectorSize(p->particles); ++j){
     for (size_t i = 0; i < DIMENSION; ++i){
@@ -1142,6 +1149,13 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
       }
     }
     free(oldVelocity);
+    for (size_t j = 0; j<vectorSize(p->contacts); ++j){
+      p->contacts[j].dv /= dt;
+      p->contacts[j].ct /= dt;
+#if DIMENSION==3
+      p->contacts[j].cs /= dt;
+#endif
+    }
     return 0;
   }
   for (size_t i = 0; i < vectorSize(p->position); ++i){
@@ -1169,6 +1183,13 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
       p->contactForces[j*DIMENSION+i] = (p->velocity[j*DIMENSION+i]-p->contactForces[j*DIMENSION+i])*p->particles[j].m/dt;
     }
   }
+  for (size_t j = 0; j<vectorSize(p->contacts); ++j){
+    p->contacts[j].dv /= dt;
+    p->contacts[j].ct /= dt;
+#if DIMENSION==3
+    p->contacts[j].cs /= dt;
+#endif
+  }
   free(oldVelocity);
   return 1;
 }

testcases/drag-3d/drag.py

@@ -46,9 +46,9 @@ g = np.array([0,-9.81,0])                      #gravity
 # Numerical parameters
 dt = 1e-3                                      #time step
 t = 0                                          #initial time
-tEnd = 35                                      #final time
+tEnd = 10                                      #final time
 ii = 0                                         #iteration number
-outf = 10                                      #iterations between data frames  
+outf = 100                                      #iterations between data frames  
 # Define output directory
 outputdir = "drag" + str(int(listeVit[index])) + "/" + str(numDep)
 if not os.path.isdir(outputdir) :
@@ -64,12 +64,12 @@ if os.path.exists(filename):
 #
 #Injected particles
 p = scontact.ParticleProblem(3,True,True)
-p.read_vtk("depot/" + "output3mm",5)
+p.read_vtk("outputhahalol",5)
 p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.52))
 p.velocity()[:,1] = vit
 #Flow particles
 p2 = scontact.ParticleProblem(3,True,True)
-p2.read_vtk("depot/" + "output3mm",5)
+p2.read_vtk("outputhahalol",5)
 p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
 p2.clear_boundaries()
 p2.load_msh_boundaries("mesh.msh", ["Inner","Left", "Right","Front","Rear"],material = "Steel")
@@ -97,10 +97,7 @@ def accumulate(F):
   F+=np.sum(p2.get_boundary_impulsion("Inner"),axis=0)
 while t < tEnd : 
      p2.forced_flag()[(p2.position()[:,1] + p2.r()[:,0] <= -0.5)] = 1
-     #if t > 5 or t<-1*dt: 
      p2.velocity()[p2.forced_flag()==1,:] = [0, vit, 0]
-     #else : 
-     #  p2.velocity()[p2.forced_flag()==1,:] = [0, 0, 0]
      p2.omega()[p2.forced_flag()==1,:] = [0,0,0]
      if  max(p2.position()[:,1]+p2.r()[:,0]) < 0.5 and t < 21 :
        p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
@@ -109,7 +106,6 @@ while t < tEnd :
      time_integration.particle_changed(fluid,p2)
      time_integration.iterate(fluid,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
      with open(filename,"a") as file1:
-        F /= dt
         F += fluid.boundary_force()[0,:]
         meanVel =  abs((np.mean(np.linalg.norm(p2.velocity()[(p2.position()[:,1] > -0.45)*(p2.position()[:,1]<-0.3),:2],axis=1)) + np.mean(np.linalg.norm(p2.velocity()[(p2.position()[:,1] < 0.45)*(p2.position()[:,1]>0.3),:2],axis=1)))/2)
         file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+";"+str(meanVel)+"\n")

testcases/hopper-2d/hopper.py

@@ -117,7 +117,7 @@ while t < tEnd :
     speedAbove = np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1] > 0.3),1])
     speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
     with open(filename,"a") as file1:
-        file1.write(str(F[0]/dt)+";"+str(F[1]/dt)+";"
+        file1.write(str(F[0])+";"+str(F[1])+";"
                     +str(t)+";"+str(0 if np.isnan(speedAbove) else speedAbove)+";"
                     +str(0 if np.isnan(speedBelow) else speedBelow)+"\n")
     t += dt

testcases/hopper-2d/hopperFluid.py

@@ -124,7 +124,6 @@ while t < tEnd :
     speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
     #print("Speed Above",speedAbove,"Speed Below",speedBelow)
     with open(filename,"a") as file1:
-        F /=dt
         F += fluid.boundary_force()[0,:];
         file1.write(str(F[0])+";"+str(F[1])+";"+str(t)+";"
                    +str(0 if np.isnan(speedAbove) else speedAbove)+";"+str(0 if np.isnan(speedBelow) else speedBelow)+"\n")

testcases/hopper-3d/hopper.py

@@ -35,7 +35,7 @@ import sys
 #index = int(inter/5)
 #numM = liste[index]
 #num = inter%5
-use_lmgc = True
+use_lmgc = False
 num = int(sys.argv[1])
 friction= [[0.0,0.0],[0.0,0.1],[0.0,0.2],[0.0,0.3],[0.0,0.4],[0.0,0.5],
 [0.1,0.0],[0.1,0.1],[0.1,0.2],[0.1,0.3],[0.1,0.4],[0.1,0.5],
@@ -109,8 +109,8 @@ p2.remove_particles_flag((p2.position()[:,1] + p2.r()[:,0] >1.1))
 p2.remove_particles_flag((p2.position()[:,0] + p2.r()[:,0] <0.18))
 p2.remove_particles_flag((p2.position()[:,0] + p2.r()[:,0] >-0.18))
 p2.position()[:,1] -= 0.5
-#p.clear_boundaries()
-#p.load_msh_boundaries("mesh04" + ".msh", ["Inner","Right","Left","RightUp","RightDown","LeftUp","LeftDown","Front","Rear","LockDown"], material = "Steel")
+p.clear_boundaries()
+p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","RightUp","RightDown","LeftUp","LeftDown","Front","Rear","LockDown"], material = "Steel")
 
 p.write_vtk(outputdir, 0, 0)
 # 
@@ -123,20 +123,20 @@ while t < tEnd :
     if t>=10 and opened == 0:
         opened = 1
         p.clear_boundaries()
-        p.load_msh_boundaries("mesh04" + ".msh", ["Inner","Right","Left","RightUp","RightDown","LeftUp","LeftDown","Front","Rear"],material="Steel")
+        p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","RightUp","RightDown","LeftUp","LeftDown","Front","Rear"],material="Steel")
 #Adding new particles
-    #if (max(p.position()[:,1])+max(p.r()) <0.6) and t > 10 and t < 25:
-    #  p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
+    if (max(p.position()[:,1])+max(p.r()) <0.6) and t > 10 and t < 25:
+      p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
 #Solving the contacts
     F = np.zeros(3)
-    time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=1e-6,external_particles_forces=g*p.mass())#,after_sub_iter= lambda subdt : accumulate(F))
+    time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=1e-6,external_particles_forces=g*p.mass(),after_sub_iter= lambda subdt : accumulate(F))
 #Removing particles outside the hopper
-    #p.remove_particles_flag( (p.position()[:,1] >-0.6))
+    p.remove_particles_flag( (p.position()[:,1] >-0.6))
 #Computing mean velocity and writing to file
     speedAbove = np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1] > 0.3),1])
     speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
     with open(filename,"a") as file1:
-        file1.write(str(F[0]/dt)+";"+str(F[1]/dt)+";"+str(F[2]/dt)+";"
+        file1.write(str(F[0])+";"+str(F[1])+";"+str(F[2])+";"
                     +str(t)+";"+str(0 if np.isnan(speedAbove) else speedAbove)+";"
                     +str(0 if np.isnan(speedBelow) else speedBelow)+"\n")
     t += dt

testcases/hopper-3d/hopperFluid.py

@@ -95,7 +95,6 @@ while t > tEnd :
     time_integration.particle_changed(fluid,p)
     time_integration.iterate(fluid,p,dt,min_nsub=8,contact_tol=1e-5,external_particles_forces=g*p.mass(),after_sub_iter=lambda subdt : accumulate(F))
     with open(filename,"a") as file1:
-        F /= dt
         F +=  fluid.boundary_force()[0,:]
         file1.write(str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+";"+
                    str(np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1]>0.25),1]))