correcting bugs

scontact/scontact.c

@@ -197,7 +197,7 @@ void particle_problem_get_delassus(ParticleProblem *p, double *delassus){
     double w[DIMENSION][DIMENSION] = {0};
 #if DIMENSION==2
     Body *body = &p->bodies[p->particles[i].body];
-    body_get_point(body,&p->particles[i].x,r);
+    body_get_point(body,p->particles[i].x,r);
     r[0] -= body->position[0];
     r[1] -= body->position[1];
     w[0][1] = w[1][0] = -r[0]*r[1]*body->invertI;
@@ -223,20 +223,17 @@ void particle_problem_get_particles_position(ParticleProblem *p, double *x){
 void particle_problem_get_particles_velocity(ParticleProblem *p, double *v){
   for(int i = 0;i<vectorSize(p->particles);i++){
     const Body *body = &p->bodies[p->particles[i].body];
-    double r[3], res[3], omega[3];
-    for (int d = 0;d<DIMENSION;d++){
-#if DIMENSION==3
-      omega[d] = p->bodies[p->particles[i].body].omega[DIMENSION+d];
-#else
-      omega[d] = 0;
-#endif
-      r[d] = p->particles[i].x[d];
-    }
+    double res[DIMENSION] = {0.};
 #if DIMENSION==2
-    r[2] = 0; // TODO FIXME !!!! cross does not work in 2D
-    omega[2] =  p->bodies[p->particles[i].body].omega;
+    double r[2];
+    body_get_point(body,p->particles[i].x,r);
+    r[0] -= body->position[0];
+    r[1] -= body->position[1];
+    res[0] = -body->omega*r[1];
+    res[1] = body->omega*r[0];
+#else
+//TODO
 #endif
-    _cross(omega,r,res);
     for(int d = 0;d<DIMENSION;d++){
       v[i*DIMENSION + d] = body->velocity[d] + res[d];
     }
@@ -251,7 +248,7 @@ static double contact_update_particle_particle(const Contact *c, const ParticleP
   double nnorm = 0;
   for (int k = 0; k < DIMENSION; ++k) {
     basis[0][k] = x1[k] - x0[k];
-    nnorm += basis[0][k] * basis[0][k];
+    nnorm += basis[0][k] * basis[0][k]; 
   }
   nnorm = sqrt(nnorm);
   double onorm = 1./nnorm;
@@ -1195,9 +1192,7 @@ static int particle_problem_solve_contact_queue(ParticleProblem *p, double dt, d
     activeContact[i] = 1;
   }
   int iter =0;
-  static size_t itertot = 0;
   for (size_t ic = fifo_pop(queue); ic != FIFO_EMPTY; ic=fifo_pop(queue)){
-    itertot ++;
     Contact *c = &p->contacts[ic];
     if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){ 
       printf("Warning : Convergence not reached\n");
@@ -1230,7 +1225,6 @@ static int particle_problem_solve_contact_queue(ParticleProblem *p, double dt, d
     }
     iter++;
   }
-  printf("iter = %d / %e\n", iter, (double)itertot);
   fifo_free(queue);
   free(activeContact);
   free(contactByBodyP);

testcases/depot-2d/depot.py

@@ -48,18 +48,27 @@ def genInitialPosition(filename, r, H, ly, lx, rhop) :
     """
     # Particles structure builder
     p = scontact.ParticleProblem(2)
+    #p2 = scontact.ParticleProblem(2)
     # Load mesh.msh file specifying physical boundaries names
-    p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])
-
+    #p2.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])
+    #for i in range(p2.segments().shape[0]):
+    #  p.add_boundary_segment((0,0),(1,1),0)
+    #p.segments()[:] = p2.segments()[:]
+    #print(p2.segments()[:],p.segments()[:])
+    p.add_boundary_segment((-0.2,-0.3),(-0.2,0.3),0)
+    p.add_boundary_segment((-0.2,0.3),(0.2,0.3),0)
+    p.add_boundary_segment((0.2,0.3),(0.2,-0.3),0)
+    p.add_boundary_segment((0.2,-0.3),(-0.2,-0.3),0)
     #Definition of the points where the particles are located
-    x = np.arange(-lx/2+r-1e-8, lx/2-r+1e-8, 2*r)
-    y = np.arange(H/2-r, H/2-ly+r, -2*r)
+    x = np.arange(-lx/2+r-1e-8, lx/2-r+1e-8, 2.05*r)
+    y = np.arange(H/2-r, H/2-ly+r, -2.05*r)
     x, y = np.meshgrid(x, y)
     x = x.flat
     y = y.flat
-    # Add a grain at each centre position
+    #Add a grain at each centre position
     for i in range(len(x)) :
         p.add_particle(np.array([x[i], y[i]]), r, r**2 * np.pi * rhop)
+    #p.add_particle_body(np.array([[0,0],[6*r,0],[12*r,0]]),np.array([3*r,3*r,3*r]),np.pi*rhop*np.array([3*r,3*r,3*r])**2)
     p.write_vtk(filename,0,0)
 
 # Define output directory
@@ -75,7 +84,7 @@ rho = 1000                                      # fluid density
 nu = 1e-6                                       # kinematic viscosity
 
 # Numerical parameters
-outf = 5                                        # number of iterations between output files
+outf = 10                                        # number of iterations between output files
 dt = 1e-3                                     # time step
 tEnd = 100                                      # final time
 
@@ -91,7 +100,6 @@ H = 0.6                                         # domain height
 genInitialPosition(outputdir, r, H, ly, lx, rhop)
 p = scontact.ParticleProblem(2)
 p.read_vtk(outputdir,0)
-print(p.body_position())
 
 
 # Initial time and iteration