Commit 36abc17f authored by Nathan Coppin's avatar Nathan Coppin
Browse files

added fric in 2D

parent a061413c
Pipeline #8660 failed with stages
in 3 minutes and 4 seconds
......@@ -44,11 +44,11 @@ assert(np.lib.NumpyVersion(np.__version__) >= "1.17.0")
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = np.ctypeslib.load_library("libscontact2",dir_path)
#lib2f = np.ctypeslib.load_library("libscontact2f",dir_path)
#lib2fwr = np.ctypeslib.load_library("libscontact2fwr",dir_path)
#lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
#lib3f = np.ctypeslib.load_library("libscontact3f",dir_path)
#lib3fwr = np.ctypeslib.load_library("libscontact3fwr",dir_path)
lib2f = np.ctypeslib.load_library("libscontact2f",dir_path)
lib2fwr = np.ctypeslib.load_library("libscontact2fwr",dir_path)
lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
lib3f = np.ctypeslib.load_library("libscontact3f",dir_path)
lib3fwr = np.ctypeslib.load_library("libscontact3fwr",dir_path)
is_64bits = sys.maxsize > 2**32
......
......@@ -26,23 +26,23 @@ set(SRC
include_directories(. ../tools)
#add_library(scontact3f SHARED ${SRC})
#target_compile_definitions(scontact3f PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
add_library(scontact3f SHARED ${SRC})
target_compile_definitions(scontact3f PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
add_library(scontact2 SHARED ${SRC})
target_compile_definitions(scontact2 PUBLIC "-DDIMENSION=2")
#add_library(scontact2f SHARED ${SRC})
#target_compile_definitions(scontact2f PUBLIC -DDIMENSION=2 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
add_library(scontact2f SHARED ${SRC})
target_compile_definitions(scontact2f PUBLIC -DDIMENSION=2 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
#add_library(scontact2fwr SHARED ${SRC})
#target_compile_definitions(scontact2fwr PUBLIC -DDIMENSION=2 -DFRICTION_ENABLED=1)
add_library(scontact2fwr SHARED ${SRC})
target_compile_definitions(scontact2fwr PUBLIC -DDIMENSION=2 -DFRICTION_ENABLED=1)
#add_library(scontact3 SHARED ${SRC})
#target_compile_definitions(scontact3 PUBLIC "-DDIMENSION=3")
add_library(scontact3 SHARED ${SRC})
target_compile_definitions(scontact3 PUBLIC "-DDIMENSION=3")
#add_library(scontact3fwr SHARED ${SRC})
#target_compile_definitions(scontact3fwr PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1)
add_library(scontact3fwr SHARED ${SRC})
target_compile_definitions(scontact3fwr PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1)
......
......@@ -131,6 +131,7 @@ static void contactBuildBasis(Contact *c) {
}
static int particleInitContact(ParticleProblem *p, size_t id0, Particle *p0, double *x0, size_t id1, Particle *p1, double *x1, double alert, Contact *c) {
if(p->ForcedFlag[id0]*p->ForcedFlag[id1]) return 0;
double nnorm = 0;
c->dv = 0;
c->ct = 0;
......@@ -206,7 +207,6 @@ 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;
for (int i = 0; i < DIMENSION; ++i) {
c->n[i] /= nnorm * (d->r < 0 ? -1 : 1);
}
......@@ -238,17 +238,17 @@ 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) {
static double diskIsInside(const void *pd, const double *x, const double *n, double vnfree, double *vbnd, double dv) {
const Disk *d = (const Disk*)pd;
double vs[DIMENSION];
for (int i = 0; i < DIMENSION; ++i) {
vnfree -= d->v[i] * n[i];
vbnd[i] = d->v[i];
vnfree -= vbnd[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) {
static double segmentProjectionIsInside(const void *ps, const double *x, const double *n, double vnfree,double *vbnd, double dv) {
const Segment *s = (const Segment*)ps;
double alpha = 0, beta = 0;
for (int i = 0; i < DIMENSION; ++i) {
......@@ -257,11 +257,10 @@ static double segmentProjectionIsInside(const void *ps, const double *x, const d
beta += (x[i] - s->p[1][i]) * d;
}
if (alpha < 0 || beta > 0) return 0;
double vs[DIMENSION];
double xi = alpha/(alpha-beta);
for (int i = 0; i < DIMENSION; ++i) {
vs[i] = (1-xi)*s->v[0][i] + xi*s->v[1][i];
vnfree -= vs[i] * n[i];
vbnd[i] = (1-xi)*s->v[0][i] + xi*s->v[1][i];
vnfree -= vbnd[i] * n[i];
}
if (dv > fabs(vnfree)) dv = fabs(vnfree);
return dv;
......@@ -294,12 +293,10 @@ static int segmentInitContact(size_t id, const Segment *s, size_t particle, cons
}
contactBuildBasis(c);
c->D = nnorm - p->r;
//if (c->D < 0 && segmentProjectionIsInside(s, x)) c->D = 0;
c->a0 = 0;
c->a1 = 1;
c->type = PARTICLE_SEGMENT;
return c->D < alert;
//return c->D >= 0 && c->D < alert;
}
......@@ -397,64 +394,58 @@ static double particleProblemGetMu(const ParticleProblem *p, int mat0, int mat1)
#endif
static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double dt, double tol) {
if(p->ForcedFlag[c->o0] == 1 && p->ForcedFlag[c->o1] == 1){
return 1;
}
double vn = c->dv;
for (int i = 0; i<DIMENSION; ++i)
vn += (p->velocity[c->o0*DIMENSION+i] - p->velocity[c->o1*DIMENSION+i]) * c->n[i];
double dp = fmax(0., vn - c->D/dt);
#if FRICTION_ENABLED
#if DIMENSION == 2
double dv = fmax(0., vn - c->D/dt);
dv -= c->dv;
int criterion = fabs(dv) <= tol/dt;
coordAdd(&p->velocity[c->o0*DIMENSION], -dv*c->a0, c->n);
coordAdd(&p->velocity[c->o1*DIMENSION], dv*c->a1, c->n);
c->dv += dv;
#if !FRICTION_ENABLED
return criterion;
#else
dv = c->dv;
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o0],p->particleMaterial[c->o1]);
if(mu<1e-10) return criterion;
double ct = 0;
double vt = c->ct;
for(int i = 0;i<2;i++){
vt += (p->velocity[c->o0*2+i]-p->velocity[c->o1*2+i])*c->t[i];
#if DIMENSION==3
double cs = 0;
double vs = c->cs;
#endif
for(int i = 0;i<DIMENSION;i++){
vt += (p->velocity[c->o0*DIMENSION+i]-p->velocity[c->o1*DIMENSION+i])*c->t[i];
#if DIMENSION==3
vs += (p->velocity[c->o0*3+i]-p->velocity[c->o1*3+i])*c->s[i];
#endif
}
#if ROTATION_ENABLED
#if ROTATION_ENABLED && DIMENSION==2
vt += p->omega[c->o0]*p->particles[c->o0].r+p->omega[c->o1]*p->particles[c->o1].r;
#endif
#else
double ct, cs = 0;
double vt = c->ct;
double vs = c->cs;
#if ROTATION_ENABLED
#elif ROTATION_ENABLED && DIMENSION==3
double res0[3], res1[3];
_cross(&p->omega[c->o0*3], c->n, res0);
_cross(&p->omega[c->o1*3], c->n, res1);
#endif
for(int i = 0;i<DIMENSION;i++){
vt += (p->velocity[c->o0*3+i]-p->velocity[c->o1*3+i])*c->t[i];
vs += (p->velocity[c->o0*3+i]-p->velocity[c->o1*3+i])*c->s[i];
#if ROTATION_ENABLED
vt += (res1[i]*p->particles[c->o1].r+res0[i]*p->particles[c->o0].r)*c->t[i];
vs += (res1[i]*p->particles[c->o1].r+res0[i]*p->particles[c->o0].r)*c->s[i];
#endif
}
#endif
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o0],p->particleMaterial[c->o1]);
#if DIMENSION==2
if (fabs(vt)> 3*mu*dp) vt *= 3*mu*dp/fabs(vt);
ct = vt;
ct = vt*fmin(fabs(vt),3*mu*dv)/(fabs(vt) == 0 ? 1 : fabs(vt));
#else
double norm = sqrt(vt*vt + vs*vs);
if(norm>3.5*mu*dp){
vt *= 3.5*mu*dp/norm;
vs *= 3.5*mu*dp/norm;
}
ct = vt;
cs = vs;
ct = vt*fmin(norm,3.5*mu*dv)/(norm == 0 ? 1 : norm);
cs = vs*fmin(norm,3.5*mu*dv)/(norm == 0 ? 1 : norm);
#endif
#if ROTATION_ENABLED
ct -= c->ct;
#if DIMENSION==2
#if ROTATION_ENABLED && DIMENSION==2
coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*c->a0/3, c->t);
coordAdd(&p->velocity[c->o1 * DIMENSION], ct*c->a1/3, c->t);
p->omega[c->o0] -= (1/p->particles[c->o0].r)*ct*c->a0*2/3;
p->omega[c->o1] -= (1/p->particles[c->o1].r)*ct*c->a1*2/3;
#else
#elif ROTATION_ENABLED && DIMENSION==3
cs -= c->cs;
coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*c->a0*2/7, c->t);
coordAdd(&p->velocity[c->o1 * DIMENSION], ct*c->a1*2/7, c->t);
......@@ -465,78 +456,63 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
p->omega[c->o1*3 + i] += (1/p->particles[c->o1].r)*c->a1*(cs*c->t[i]-ct*c->s[i])*5/7;
}
c->cs += cs;
#endif
#else
ct -= c->ct;
#elif DIMENSION==2
coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*c->a0, c->t);
coordAdd(&p->velocity[c->o1 * DIMENSION], ct*c->a1, c->t);
#if DIMENSION==3
#elif DIMENSION==3
cs -= c->cs;
coordAdd(&p->velocity[c->o0 * DIMENSION], -cs*c->a0, c->s);
coordAdd(&p->velocity[c->o1 * DIMENSION], cs*c->a1, c->s);
c->cs += cs;
#endif
#endif
c->ct += ct;
#endif
dp -= c->dv;
coordAdd(&p->velocity[c->o0*DIMENSION], -dp*c->a0, c->n);
coordAdd(&p->velocity[c->o1*DIMENSION], dp*c->a1, c->n);
c->dv += dp;
int rfriction = 1;
#if FRICTION_ENABLED
rfriction &= fabs(ct) <= tol/dt;
int rfriction = fabs(ct) <= tol/dt;
#if DIMENSION == 3
rfriction &= fabs(cs) <= tol/dt;
#endif
return (criterion && rfriction);
#endif
return (fabs(dp) <= tol/dt && rfriction);
}
#if FRICTION_ENABLED
static int contactParticleBndFrictionSolve(ParticleProblem *p, Contact *c, Boundary *bnd, double dp, double toldt) {
#if DIMENSION == 2
static int contactParticleBndFrictionSolve(ParticleProblem *p, Contact *c, Boundary *bnd, double dv, double toldt, double *vbnd) {
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],bnd->material);
if(mu<1e-10) return 1;
double ct = 0;
double vt = c->ct + (p->velocity[c->o1*2+0]-bnd->v[0])*c->t[0]+ (p->velocity[c->o1*2+1]-bnd->v[1])*c->t[1]+bnd->vt;
#if ROTATION_ENABLED
vt += p->omega[c->o1]*p->particles[c->o1].r;
#endif
#else
double ct, cs = 0;
double vt = c->ct ;
double vt = c->ct;
#if DIMENSION == 3
double cs = 0;
double vs = c->cs ;
#if ROTATION_ENABLED
#endif
for(int i=0;i<DIMENSION;i++){
vt += (p->velocity[c->o1*DIMENSION+i]- vbnd[i])*c->t[i];
#if DIMENSION==3
vs += (p->velocity[c->o1*3+i]- vbnd[i])*c->s[i];
#endif
}
#if ROTATION_ENABLED && DIMENSION==2
vt += p->omega[c->o1]*p->particles[c->o1].r;
#elif ROTATION_ENABLED && DIMENSION==3
double res1[3];
_cross(&p->omega[c->o1*3], c->n, res1);
#endif
for(int i=0;i<DIMENSION;i++){
vt += (p->velocity[c->o1*3+i]- bnd->v[i])*c->t[i];
vs += (p->velocity[c->o1*3+i]- bnd->v[i])*c->s[i];
#if ROTATION_ENABLED
vt += (res1[i]*p->particles[c->o1].r)*c->t[i];
vs += (res1[i]*p->particles[c->o1].r)*c->s[i];
#endif
}
#endif
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],bnd->material);
#if DIMENSION==2
if (fabs(vt)> 3*mu*dp) vt *= 3*mu*dp/fabs(vt);
ct = vt;
ct = vt*fmin(fabs(vt),3*mu*dv)/(fabs(vt) == 0 ? 1 : fabs(vt));
#else
double norm = sqrt(vt*vt + vs*vs);
if(norm>3.5*mu*dp){
vt *= 3.5*mu*dp/norm;
vs *= 3.5*mu*dp/norm;
}
ct = vt;
cs = vs;
ct = vt*fmin(norm,3.5*mu*dv)/(norm == 0 ? 1 : norm);
cs = vs*fmin(norm,3.5*mu*dv)/(norm == 0 ? 1 : norm);
#endif
#if ROTATION_ENABLED
ct -= c->ct;
#if DIMENSION==2
#if ROTATION_ENABLED && DIMENSION==2
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct/3, c->t);
p->omega[c->o1] -= (1/p->particles[c->o1].r)*ct*2/3;
#else
#elif ROTATION_ENABLED && DIMENSION==3
cs -= c->cs;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*2/7, c->t);
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*2/7, c->s);
......@@ -544,9 +520,7 @@ static int contactParticleBndFrictionSolve(ParticleProblem *p, Contact *c, Bound
p->omega[c->o1*3 + i] += (1/p->particles[c->o1].r)*(cs*c->t[i]-ct*c->s[i])*5/7;
}
c->cs += cs;
#endif
#else
ct -= c->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct, c->t);
#if DIMENSION==3
cs -= c->cs;
......@@ -555,40 +529,38 @@ static int contactParticleBndFrictionSolve(ParticleProblem *p, Contact *c, Bound
#endif
#endif
c->ct += ct;
int criterion3D = 1;
int rfriction = fabs(ct) <= toldt;
#if DIMENSION == 3
criterion3D &= fabs(cs) <= toldt;
rfriction &= fabs(cs) <= toldt;
#endif
return (fabs(ct) <= toldt && criterion3D);
return rfriction;
}
#endif
typedef double IsInsideFunc(const void*, const double *x, const double *n, double vnfree, double dv);
typedef double IsInsideFunc(const void*, const double *x, const double *n, double vnfree, double *vbnd, double dv);
static int contactParticleBoundarySolve(Contact *c, ParticleProblem *p, double dt, double tol, Boundary *bnd, IsInsideFunc *isInsideFunc) {
if(p->ForcedFlag[c->o1])
return 1;
double vlocfree[DIMENSION];
double vn = 0;
double *vbnd = (double *)malloc(sizeof(double)*DIMENSION);
for (int i= 0; i < DIMENSION; ++i) {
vlocfree[i] = p->velocity[c->o1*DIMENSION + i] + c->n[i] * c->dv;
vn += vlocfree[i] * c->n[i];
}
double dv = fmax(0, vn - c->D/dt);
double r = p->particles[c->o1].r;
if (isInsideFunc) {
double xn[DIMENSION];
for (int i = 0; i < DIMENSION; ++i)
xn[i] = p->position[c->o1*DIMENSION + i] + r*c->n[i] + vlocfree[i]*dt;//vlocfree[i]*c->D / fmax(vn,1e-16);
dv = isInsideFunc(bnd, xn, c->n, vn, dv);
}
double xn[DIMENSION];
for (int i = 0; i < DIMENSION; ++i)
xn[i] = p->position[c->o1*DIMENSION + i] + r*c->n[i] + vlocfree[i]*dt;
dv = isInsideFunc(bnd, xn, c->n, vn, vbnd, dv);
int rfriction = 1;
#if FRICTION_ENABLED
rfriction = contactParticleBndFrictionSolve(p,c,bnd,dv,tol/dt);
rfriction = contactParticleBndFrictionSolve(p,c,bnd,dv,tol/dt,vbnd);
#endif
dv -= c->dv;
coordAdd(&p->velocity[c->o1*DIMENSION], -dv, c->n);
c->dv += dv;
free(vbnd);
return (fabs(dv) <= tol/dt && rfriction);
}
......@@ -669,7 +641,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
for (size_t j = 0; j < vectorSize(found); ++j) {
ntot++;
Contact *c= vectorPush(&p->contacts);
if(!particleInitContact(p, i, &p->particles[i], &p->position[i * DIMENSION], found[j], &p->particles[found[j]], &p->position[found[j] * DIMENSION], alert, c)||(p->ForcedFlag[c->o0]&&p->ForcedFlag[c->o1]))
if(!particleInitContact(p, i, &p->particles[i], &p->position[i * DIMENSION], found[j], &p->particles[found[j]], &p->position[found[j] * DIMENSION], alert, c))
vectorPop(p->contacts);
else if ((cold = findContactSorted(c, oldContacts, &iold))) {
coordAdd(&p->velocity[c->o0 * DIMENSION], -cold->dv * c->a0, c->n);
......@@ -697,14 +669,11 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
}
#endif
#else
#if DIMENSION==2
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0, c->t);
coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*c->a1, c->t);
#else
#if DIMENSION==3
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->cs*c->a0, c->s);
coordAdd(&p->velocity[c->o1 * DIMENSION], c->cs*c->a1, c->s);
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0, c->t);
coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*c->a1, c->t);
#endif
#endif
#endif
......@@ -719,17 +688,15 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
cellSearch(tree, amin, amax, &found);
for (size_t j = 0; j < vectorSize(found); ++j) {
Contact *c = vectorPush(&p->contacts);
if(!diskInitContact(i, &p->disks[i], found[j], &p->particles[found[j]], &p->position[found[j] * DIMENSION], alert, c) || p->ForcedFlag[c->o1])
if(p->ForcedFlag[found[j]] ||!diskInitContact(i, &p->disks[i], found[j], &p->particles[found[j]], &p->position[found[j] * DIMENSION], alert, c))
vectorPop(p->contacts);
else if ((cold = findContactSorted(c, oldContacts, &iold))) {
coordAdd(&p->velocity[c->o1 * DIMENSION], -cold->dv, c->n);
c->dv = cold->dv;
#if FRICTION_ENABLED
#if DIMENSION == 2
c->ct = cold->ct;
#else
#if DIMENSION == 3
c->cs = cold->cs;
c->ct = cold->ct;
#endif
#if ROTATION_ENABLED
#if DIMENSION==2
......@@ -743,11 +710,9 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
}
#endif
#else
#if DIMENSION==2
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct, c->t);
#else
#if DIMENSION==3
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs, c->s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct, c->t);
#endif
#endif
#endif
......@@ -762,17 +727,15 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
cellSearch(tree, amin, amax, &found);
for (size_t j = 0; j < vectorSize(found); ++j) {
Contact *c = vectorPush(&p->contacts);
if (!segmentInitContact(i, &p->segments[i], found[j], &p->particles[found[j]], &p->position[found[j] * DIMENSION], alert, c) || p->ForcedFlag[c->o1])
if (p->ForcedFlag[found[j]]||!segmentInitContact(i, &p->segments[i], found[j], &p->particles[found[j]], &p->position[found[j] * DIMENSION], alert, c))
vectorPop(p->contacts);
else if ((cold = findContactSorted(c, oldContacts, &iold))) {
coordAdd(&p->velocity[c->o1 * DIMENSION], -cold->dv, c->n);
c->dv = cold->dv;
#if FRICTION_ENABLED
#if DIMENSION == 2
c->ct = cold->ct;
#else
#if DIMENSION == 3
c->cs = cold->cs;
c->ct = cold->ct;
#endif
#if ROTATION_ENABLED
#if DIMENSION==2
......@@ -786,11 +749,9 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
}
#endif
#else
#if DIMENSION==2
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct, c->t);
#else
#if DIMENSION==3
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs, c->s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct, c->t);
#endif
#endif
#endif
......@@ -806,7 +767,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
cellSearch(tree, amin, amax, &found);
for (int j = 0; j < vectorSize(found); ++j) {
Contact *c = vectorPush(&p->contacts);
if (!triangleInitContact(i, &p->triangles[i], found[j], &p->particles[found[j]], &p->position[found[j] * DIMENSION], alert, c) || p->ForcedFlag[c->o1])
if (p->ForcedFlag[found[j]]||!triangleInitContact(i, &p->triangles[i], found[j], &p->particles[found[j]], &p->position[found[j] * DIMENSION], alert, c))
vectorPop(p->contacts);
else if ((cold = findContactSorted(c, oldContacts, &iold))) {
coordAdd(&p->velocity[c->o1 * DIMENSION], -cold->dv, c->n);
......@@ -901,29 +862,23 @@ void randomize(int arr[], int n) {
static int _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol) {
int converged = 0;
//#if FRICTION_ENABLED
int iter = 0;
//#endif
int *order = malloc(vectorSize(p->contacts)*sizeof(int));
for (int i = 0; i < vectorSize(p->contacts); ++i) {
order[i] = i;
}
while (! converged) {
randomize(order, vectorSize(p->contacts));
//#if FRICTION_ENABLED
//randomize(order, vectorSize(p->contacts));
if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){
printf("Warning : Convergence not reached\n");
return 0;
}
//#endif
converged = 1;
for (int i = 0; i < vectorSize(p->contacts); ++i) {
converged &= contactSolve(p,&p->contacts[order[i]],dt,tol);
}
//#if FRICTION_ENABLED
iter++;
//#endif
}
return 1;
}
......@@ -1090,12 +1045,9 @@ static int _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, dou
fifoPush(queue, i);
activeContact[i] = 1;
}
#if FRICTION_ENABLED
int iter =0;
#endif
for (size_t ic = fifoPop(queue); ic != FIFO_EMPTY; ic=fifoPop(queue)){
Contact *c = &p->contacts[ic];
#if FRICTION_ENABLED
if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){
printf("Warning : Convergence not reached\n");
fifoFree(queue);
......@@ -1104,7 +1056,6 @@ static int _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, dou
free(contactByParticle);
return 0;
}
#endif
int converged = contactSolve(p,c,dt,tol);
if (!converged) {
if(c->type == PARTICLE_PARTICLE) {
......@@ -1122,9 +1073,7 @@ static int _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, dou
}
}
}
#if FRICTION_ENABLED
iter++;
#endif
activeContact[ic] = 0;
}
fifoFree(queue);
......
......@@ -4,25 +4,27 @@ import numpy as np
os.makedirs("output",exist_ok=True)
p = scontact.ParticleProblem(2,False,False)
p = scontact.ParticleProblem(2,True)
a = 1
r = 0.05
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,0],[-a/3-r,0],"bnd",material="Steel")
p.add_boundary_segment([a/3+r,0],[a,0],"bnd",material="Steel")
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):
m = np.random.random()
p.add_particle((x+m*r/4,y),r-m*r/4,1)
#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",material="Steel")
#p.add_boundary_disk([-a/3-r,0],0,"bnd",material="Steel")
p.add_boundary_segment([-a,-a],[-a,a],"bnd",material="Steel")
p.add_boundary_segment([-a,a],[a,a],"bnd",material="Steel")
p.add_boundary_segment([a,a],[a,-a],"bnd",material="Steel")
p.add_boundary_segment([a,-a],[-a,-a],"bnd",material="Steel")
for x in np.arange(-a+r,a-r,3*r):
for y in np.arange(-a+r,-r,3*r):
m = 0.0001#np.random.random()
p.add_particle((x+m*r/4,y),r-m*r/4,1,"Sand")
p.set_friction_coefficient(0.2,"Sand","Sand")
p.set_friction_coefficient(0.2,"Sand","Steel")
i = 0
f = np.zeros_like(p.velocity())