Works with relaxation

scontact/scontact.c

@@ -375,13 +375,13 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
   double muk = p->muk[index];
   if(p->frictionFlag && mus!=0){
     if(vt>(7./2.)*mus*(*dp)){
-      *ct = mus*(*dp);//printf("Slip\n");
+      *ct = muk*(*dp);//printf("Slip\n");
     }
     else if(vt<(-7./2.)*mus*(*dp)){
-      *ct = -mus*(*dp);//printf("Slip\n");
+      *ct = -muk*(*dp);//printf("Slip\n");
     }
     else{
-      *ct = (2./7.)*vt;//printf("Stick\n");
+      *ct = 0.8*(2./7.)*vt;//printf("Stick\n");
     }
     *ct -= c->ct;
     coordAdd(&p->velocity[c->o0 * DIMENSION], -*ct*c->a0, t);
@@ -396,18 +396,13 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
   //printf("ct : %.8f dp : %.8f\n",*ct,*dp);
   coordAdd(&p->velocity[c->o0 * DIMENSION], -*dp*c->a0, c->n);
   coordAdd(&p->velocity[c->o1 * DIMENSION], *dp*c->a1, c->n);
- vt = (p->velocity[c->o0*2+0]-p->velocity[c->o1*2+0])*t[0]+
-    (p->velocity[c->o0*2+1]-p->velocity[c->o1*2+1])*t[1]+
-    p->omega[c->o0]*p->particles[c->o0].r +
-    p->omega[c->o1]*p->particles[c->o1].r;
-  //printf("vt après : %.8f\n",vt);
   #endif  
   if(fabs(*dp) < tol/dt && fabs(*ct) <tol/dt){
-    printf("b : %d b : %d ok\n",c->o0,c->o1);
+    //printf("b : %d b : %d ok\n",c->o0,c->o1);
     return 1;
   }
   else{
-    printf("b : %d b : %d not ok\n",c->o0,c->o1);
+    //printf("b : %d b : %d not ok\n",c->o0,c->o1);
     return 0;
   }
 }
@@ -442,13 +437,13 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double *ct,
   double muk = p->muk[index];
   if(p->frictionFlag && mus!=0){
     if(vt>(7./2.)*mus*(*dp)){
-      *ct = mus*(*dp);//printf("Slip\n");
+      *ct = muk*(*dp);//printf("Slip\n");
     }
     else if(vt<(-7./2.)*mus*(*dp)){
-      *ct = -mus*(*dp);//printf("Slip\n");
+      *ct = -muk*(*dp);//printf("Slip\n");
     }
     else{
-      *ct = (2./7.)*vt;//printf("Stick\n");
+      *ct = 0.8*(2./7.)*vt;//printf("Stick\n");
     }
     *ct -= c->ct;
     coordAdd(&p->velocity[c->o1 * DIMENSION], -*ct*c->a1, t);
@@ -460,8 +455,7 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double *ct,
   *dp -=c->dv;
   coordAdd(&p->velocity[c->o1 * DIMENSION], -*dp, c->n);
 //printf("dis : %d ct : %f dp : %f\n",c->o0,*ct,*dp);
-  vt = p->velocity[c->o1 * DIMENSION]*t[0] +p->velocity[c->o1*DIMENSION +1]*t[1] +  p->omega[c->o1]*p->particles[c->o1].r + p->disks[c->o0].vt;
-  //printf("vt après : %f\n",vt);
+ 
 
   #else
 
@@ -507,13 +501,13 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double *d
   double muk = p->muk[index];
   if(p->frictionFlag && mus!=0){
     if(vt>(7./2.)*mus*(*dp)){
-      *ct = mus*(*dp);//printf("Slip\n");
+      *ct = muk*(*dp);//printf("Slip\n");
     }
     else if(vt<(-7./2.)*mus*(*dp)){
-      *ct = -mus*(*dp);//printf("Slip\n");
+      *ct = -muk*(*dp);//printf("Slip\n");
     }
     else{
-      *ct = (2./7.)*vt;//printf("Stick\n");
+      *ct = 0.8*(2./7.)*vt;//printf("Stick\n");
     }
     *ct -= c->ct;
     coordAdd(&p->velocity[c->o1 * DIMENSION], -*ct*c->a1, t);
@@ -525,8 +519,6 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double *d
   *dp -=c->dv;
   //printf("seg : %d ct : %f dp : %f\n",c->o0,*ct,*dp);
   coordAdd(&p->velocity[c->o1 * DIMENSION], -*dp, c->n);
-  vt = p->velocity[c->o1 * DIMENSION]*t[0] +p->velocity[c->o1*DIMENSION +1]*t[1] +  p->omega[c->o1]*p->particles[c->o1].r + p->segments[c->o0].vt;
-  //printf("vt après : %f\n",vt);
 
   #else
   
@@ -536,11 +528,11 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double *d
   
   #endif
   if(fabs(*dp) < tol/dt &&fabs(*ct) <tol/dt){
-    printf("seg : %d b : %d ok\n",c->o0,c->o1);
+    //printf("seg : %d b : %d ok\n",c->o0,c->o1);
     return 1;
   }
   else{
-    printf("seg : %d b : %d not ok\n",c->o0,c->o1);
+    //printf("seg : %d b : %d not ok\n",c->o0,c->o1);
     return 0;
   }
 }

testcases/inclined-2d/mesh.geo

-rout = 0.2;
+rout = 1000;
 lcout = rout/20;
 
 Point(1) = {0, 0, 0};

testcases/mill-2d/mill.py

@@ -21,6 +21,7 @@
 
 #!/usr/bin/env python
 from migflow import scontact
+from migflow import lmgc90Interface
 
 import numpy as np
 import os
@@ -69,18 +70,25 @@ dt = 1e-3                               # time step
 
 #geometry parameters
 rout = 0.11                           # outer radius
-r  = 1e-2                              # grains radius
+r  = 0.1e-2                              # grains radius
 Fr = 1.3				#Froude number
 shutil.copy("mesh.msh", outputdir +"/mesh.msh")
 v = rout*np.sqrt(Fr*(-g)/rout)
 print('v = %s\n'%(v))
 #Object particles creation
 genInitialPosition(outputdir, r, rout, rhop)
-p = scontact.ParticleProblem(2)
-p.read_vtk(outputdir,0)
+use_lmgc90 = True
+if use_lmgc90 :
+    friction=1.0                                # friction coefficient
+    lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
+    p = lmgc90Interface.ParticleProblem(2)
+else :
+    p = scontact.ParticleProblem(2)
+    p.read_vtk(outputdir,0)
+
 outf = 5                                # number of iterations between output files
 #Taking friction into account
-p.add_friction([1.0, 1.0],[0.5, 0.5])
+p.add_friction([1.0, 1.0],[0.8, 0.8])
 v = 1.0;
 p.set_tangent_boundary_velocity("",v)
 print(p.velocity().shape[0])

testcases/sablier-2d/sablier.py

@@ -45,7 +45,7 @@ def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
        for j in range(y.shape[0]) :
           rr=r*(0.95+random.random()*0.05)
           #Addition of an particle object at each point
-          p.add_particle((x[i], y[j]), rr, rr**2 * np.pi * rhop);
+          p.add_particle((x[i], y[j]), rr, rr**2 * np.pi * rhop,0);
     p.write_vtk(filename,0,0)