2d no friction ok

scontact/scontact.c

@@ -206,7 +206,7 @@ static int diskInitContact(size_t id, const Disk *d, size_t particle, const Part
   }
   nnorm = sqrt(nnorm);
   c->D = (nnorm - fabs(p->r + d->r)) * (d->r < 0 ? -1 : 1);
-  if (c->D < 0) c->D = 0;
+  //if (c->D < 0) c->D = 0;
   for (int i = 0; i < DIMENSION; ++i) {
     c->n[i] /= nnorm * (d->r < 0 ? -1 : 1);
   }
@@ -238,6 +238,16 @@ static void segmentBoundingBox(const Segment *s, double *pmin, double *pmax) {
   }
 }
 
+static double diskIsInside(const void *pd, const double *x, const double *n, double vnfree, double dv) {
+  const Disk *d = (const Disk*)pd;
+  double vs[DIMENSION];
+  for (int i = 0; i < DIMENSION; ++i) {
+    vnfree -= d->v[i] * n[i];
+  }
+  if (dv > fabs(vnfree)) dv = fabs(vnfree);
+  return dv;
+}
+
 static double segmentProjectionIsInside(const void *ps, const double *x, const double *n, double vnfree, double dv) {
   const Segment *s = (const Segment*)ps;
   double alpha = 0, beta = 0;
@@ -862,7 +872,7 @@ static int contactSolve(ParticleProblem *p, Contact *c, double dt, double tol) {
       return contactParticleParticleSolve(c, p, dt,tol);
     case PARTICLE_DISK :
       return contactParticleBoundarySolve(c, p, dt, tol,
-          (Boundary*)&p->disks[c->o0],NULL);
+          (Boundary*)&p->disks[c->o0],diskIsInside);
     case PARTICLE_SEGMENT :
       return contactParticleBoundarySolve(c, p, dt, tol,
           (Boundary*)&p->segments[c->o0],segmentProjectionIsInside);
@@ -874,25 +884,46 @@ static int contactSolve(ParticleProblem *p, Contact *c, double dt, double tol) {
   }
 }
 
+void swap(int *a, int *b) {
+    int temp = *a;
+    *a = *b;
+    *b = temp;
+}
+
+void randomize(int arr[], int n) {
+    srand(time(NULL));
+    int i;
+    for(i = n-1; i > 0; i--) {
+        int j = rand() % (i+1);
+        swap(&arr[i], &arr[j]);
+    }
+}
+
 static int _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol) {
   int converged = 0;
-#if FRICTION_ENABLED
+//#if FRICTION_ENABLED
   int iter = 0;
-#endif
+//#endif
+  int *order = malloc(vectorSize(p->contacts)*sizeof(int));
+  for (int i = 0; i < vectorSize(p->contacts); ++i) {
+    order[i] = i;
+  }
+
   while (! converged) {
-#if FRICTION_ENABLED
+    randomize(order, vectorSize(p->contacts));
+//#if FRICTION_ENABLED
     if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){ 
       printf("Warning : Convergence not reached\n");
       return 0;
     }
-#endif
+//#endif
     converged = 1;
-    for (int ic = vectorSize(p->contacts)-1; ic >= 0;  --ic) {
-      converged &= contactSolve(p,&p->contacts[ic],dt,tol);
+    for (int i = 0; i < vectorSize(p->contacts); ++i) {
+      converged &= contactSolve(p,&p->contacts[order[i]],dt,tol);
     }
-#if FRICTION_ENABLED
+//#if FRICTION_ENABLED
     iter++;
-#endif
+//#endif
   }
   return 1;
 }

testcases/bndrotation/test.py

@@ -7,38 +7,57 @@ os.makedirs("output",exist_ok=True)
 p = scontact.ParticleProblem(2,False,False)
 a = 1
 r = 0.05
-p.add_boundary_segment([0,0],[0,a],"bnd")
-#p.add_boundary_segment([0,a],[-a/3,a],"bnd")
-p.add_boundary_disk([-a/3,a],0,"bnd")
-p.add_boundary_segment([-a,-a],[-a,a],"wall")
-p.add_boundary_segment([-a,a],[a,a],"wall")
-p.add_boundary_segment([a,a],[a,-a],"wall")
-p.add_boundary_segment([a,-a],[-a,-a],"wall")
-#p.add_boundary_segment([-a+3*r,a-2.1*r],[a-3*r,a-2.1*r],"wall")
+p.add_boundary_segment([-a,0],[-a/3-r,0],"bnd")
+p.add_boundary_segment([a/3+r,0],[a,0],"bnd")
+p.add_boundary_segment([-a/3+3*r,0],[a/3-3*r,0],"wall")
+p.add_boundary_disk([a/3-3*r,0],0,"wall")
+p.add_boundary_disk([-a/3+3*r,0],0,"wall")
+p.add_boundary_disk([a/3+r,0],0,"bnd")
+p.add_boundary_disk([-a/3-r,0],0,"bnd")
+p.add_boundary_segment([-a,-a],[-a,a],"bnd")
+p.add_boundary_segment([-a,a],[a,a],"bnd")
+p.add_boundary_segment([a,a],[a,-a],"bnd")
+p.add_boundary_segment([a,-a],[-a,-a],"bnd")
+
 for x in np.arange(-a+r,a-r,2*r):
     for y in np.arange(-a+r,-r,2*r):
-        p.add_particle((x,y),r-np.random.random()*r/2,1)
+        m = np.random.random()
+        p.add_particle((x+m*r/4,y),r-m*r/4,1)
 i = 0
 f = np.zeros_like(p.velocity())
 f[:,1] -= 0.0002
 
 p.write_vtk("output",0,0)
+
+dt = 1
+def rotate(x, o) :
+    xn = np.zeros_like(x)
+    r = np.hypot(x[...,0],x[...,1])
+    a = np.arctan2(x[...,1],x[...,0])+o
+    xn[...,0] = r*np.cos(a)
+    xn[...,1] = r*np.sin(a)
+    return xn
+
+wallseg0 = np.copy(p.segments()['p'][p.segments()['tag']==1])
+walldisk0 = np.copy(p.segments()['p'][p.segments()['tag']==1])
+
+t = 0
 while True :
 
+    t += dt
     x = p.segments()['p'][p.segments()['tag']==0]
-    v = np.zeros_like(x)
-    v[:,:,0] = -x[:,:,1]*0.01
-    v[:,:,1] = x[:,:,0]*0.01
-    p.segments()['v'][p.segments()['tag']==0] = v
+    p.segments()['v'][p.segments()['tag']==0] = (rotate(x,1e-2*dt)-x)/dt
 
     x = p.disks()['x'][p.disks()['tag']==0]
-    v = np.zeros_like(x)
-    v[:,0] = -x[:,1]*0.01
-    v[:,1] = x[:,0]*0.01
-    p.disks()['v'][p.disks()['tag']==0] = v
+    p.disks()['v'][p.disks()['tag']==0] = (rotate(x,1e-2*dt)-x)/dt
+
+    x = p.segments()['p'][p.segments()['tag']==1]
+    p.segments()['v'][p.segments()['tag']==1] = (wallseg0*(1+np.abs(2*np.sin(t*1e-2)))-x)/dt
+    x = p.disks()['x'][p.disks()['tag']==1]
+    p.disks()['v'][p.disks()['tag']==1] = (walldisk0*(1+np.abs(2*np.sin(t*1e-2)))-x)/dt
 
     p.set_use_queue(0)
-    p.iterate(1,f,tol=1e-8,force_motion=1)
+    p.iterate(dt,f,tol=1e-8,force_motion=1)
     i += 1
     print(i)
     p.write_vtk("output",i,i)