Commit 2fb603b2 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

correct bug in triangle ctct and in force computation

parent d3feed44
Pipeline #6870 failed with stage
in 2 minutes and 9 seconds
......@@ -24,7 +24,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,momentum=None,iter
nsub = max(min_nsub, int(np.ceil((vmax * dt * 8)/min(particles.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(particles.r()))
#If time step was split too many times, ignore the convergence and move the grains
if iteration == 5:
if iteration == 3:
for i in range(nsub) :
particles.iterate(dt/nsub, f, tol=contact_tol,force_motion=1)
return
......@@ -35,7 +35,9 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,momentum=None,iter
_advance_particles(particles,f,dt/nsub,2,contact_tol/2,momentum=momentum,iteration=iteration+1)
print("Convergence reached for subinvervals of level %d"%(iteration+1))
if momentum is not None:
momentum += particles.momentum()*nsub/dt
#print(particles.momentum()*nsub/dt)
momentum += particles.momentum()
momentum /= dt
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,momentum=None) :
......
......@@ -308,6 +308,10 @@ void particleProblemAddBoundaryTriangleTagId(ParticleProblem *p, const double x0
t->p[2][i] = x2[i];
t->v[i] = 0.;
}
#if FRICTION_ENABLED
t->vt = 0;
t->vs = 0;
#endif
}
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], const char *tag, const char *material) {
......@@ -388,8 +392,19 @@ static void build_basis(const double *a, double *b
, double *c){
b[0] = -a[2];
b[1] = a[0];
b[2] = -a[1];
b[2] = -a[1];
double provi[3];
double ddot = dot(a,b);
provi[0] = b[0] - ddot*a[0];
provi[1] = b[1] - ddot*a[1];
provi[2] = b[2] - ddot*a[2];
double norm = sqrt(provi[0]*provi[0] + provi[1]*provi[1] + provi[2]*provi[2]);
if(norm==0) norm=1.;
b[0] = provi[0]/norm;
b[1] = provi[1]/norm;
b[2] = provi[2]/norm;
_cross(a, b, c);
#else
){
b[0] = -a[1];
......@@ -454,7 +469,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
if (dp == 0.){
ct = 0;
}
//if(ct > 1e-3*7e7*(p->particles[c->o1].r+p->particles[c->o0].r)*(p->particles[c->o1].r+p->particles[c->o0].r)*(1/4)*M_PI/(p->particles[c->o1].m+p->particles[c->o0].m)) ct = 1e-3*7e7*(p->particles[c->o1].r+p->particles[c->o0].r)*(p->particles[c->o1].r+p->particles[c->o0].r)*(1/4)*M_PI/(p->particles[c->o1].m+p->particles[c->o0].m);
#else
double norme = sqrt(vt*vt + vs*vs);
if(norme>3.5*mu*dp){
......@@ -467,14 +482,16 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
cs = 0;
ct = 0;
}
//if(ct > 1e-3*7e7*(p->particles[c->o1].r+p->particles[c->o0].r)*(p->particles[c->o1].r+p->particles[c->o0].r)*(1/4)*M_PI/(p->particles[c->o1].m+p->particles[c->o0].m)) ct = 1e-3*7e7*(p->particles[c->o1].r+p->particles[c->o0].r)*(p->particles[c->o1].r+p->particles[c->o0].r)*(1/4)*M_PI/(p->particles[c->o1].m+p->particles[c->o0].m);
//if(cs > 1e-3*7e7*(p->particles[c->o1].r+p->particles[c->o0].r)*(p->particles[c->o1].r+p->particles[c->o0].r)*(1/4)*M_PI/(p->particles[c->o1].m+p->particles[c->o0].m)) cs = 1e-3*7e7*(p->particles[c->o1].r+p->particles[c->o0].r)*(p->particles[c->o1].r+p->particles[c->o0].r)*(1/4)*M_PI/(p->particles[c->o1].m+p->particles[c->o0].m);
#endif
#if ROTATION_ENABLED
ct -= c->ct;
#if DIMENSION==2
coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*(p->a/(p->a+1))*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], ct*(p->a/(p->a+1))*c->a1, t);
p->omega[c->o0] -= (1/p->particles[c->o0].r)*(1/((p->a+1)))*ct;
p->omega[c->o1] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct;
p->omega[c->o0] -= (1/p->particles[c->o0].r)*(1/((p->a+1)))*ct*c->a0;
p->omega[c->o1] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*c->a1;
#else
cs -= c->cs;
......@@ -483,10 +500,10 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
coordAdd(&p->velocity[c->o0 * DIMENSION], -cs*(p->a/(p->a+1))*c->a0, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], cs*(p->a/(p->a+1))*c->a1, s);
for (int i = 0; i<3; i++){
p->omega[c->o0*3 + i] += -(1/p->particles[c->o0].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o0].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o0*3 + i] += -(1/p->particles[c->o0].r)*(1/((p->a+1)))*ct*s[i]*c->a0
+ (1/p->particles[c->o0].r)*(1/((p->a+1)))*cs*t[i]*c->a0;
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]*c->a1
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i]*c->a1;
}
c->cs += cs;
#endif
......@@ -569,12 +586,13 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
#endif
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->diskMaterial[c->o0]);
#if DIMENSION==2
if(isnan(dp)){printf("hello\n");}
if (fabs(vt)> ((p->a+1)/p->a)*mu*dp) vt *= ((p->a+1)/p->a)*mu*dp/fabs(vt);
ct = vt;
if (dp == 0.){
ct = 0;
}
if(ct > 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m) ct = 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m;
#else
double norme = sqrt(vt*vt + vs*vs);
if(norme>3.5*mu*dp){
......@@ -587,19 +605,21 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
cs = 0;
ct = 0;
}
//if(ct > 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m) ct = 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m;
//if(cs > 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m) cs = 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m;
#endif
#if ROTATION_ENABLED
ct -= c->ct;
#if DIMENSION==2
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
p->omega[c->o1] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct;
p->omega[c->o1] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*c->a1;
#else
cs -= c->cs;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*(p->a/(p->a+1))*c->a1, s);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]*c->a1
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i]*c->a1;
}
c->cs += cs;
#endif
......@@ -656,8 +676,8 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double dt
double t[3], s[3];
build_basis(c->n,t,s);
double ct, cs = 0;
double vt = c->ct + p->segments[c->o0].vt;
double vs = c->cs + p->segments[c->o0].vs;
double vt = c->ct ;//+ p->segments[c->o0].vt;
double vs = c->cs ;//+ p->segments[c->o0].vs;
#if ROTATION_ENABLED
double res1[3];
_cross(&p->omega[c->o1*3], c->n, res1);
......@@ -673,12 +693,13 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double dt
#endif
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->segmentMaterial[c->o0]);
#if DIMENSION==2
if(isnan(dp)){printf("hello\n");}
if (fabs(vt)> ((p->a+1)/p->a)*mu*dp) vt *= ((p->a+1)/p->a)*mu*dp/fabs(vt);
ct = vt;
if (dp == 0.){
ct = 0;
}
if(ct > 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m) ct = 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m;
#else
double norme = sqrt(vt*vt + vs*vs);
if(norme>3.5*mu*dp){
......@@ -691,19 +712,21 @@ if (fabs(vt)> ((p->a+1)/p->a)*mu*dp) vt *= ((p->a+1)/p->a)*mu*dp/fabs(vt);
cs = 0;
ct = 0;
}
//if(ct > 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m) ct = 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m;
//if(cs > 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m) cs = 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m;
#endif
#if ROTATION_ENABLED
ct -= c->ct;
#if DIMENSION==2
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
p->omega[c->o1] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct;
p->omega[c->o1] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*c->a1;
#else
cs -= c->cs;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*(p->a/(p->a+1))*c->a1, s);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]*c->a1
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i]*c->a1;
}
c->cs += cs;
#endif
......@@ -750,13 +773,15 @@ static int contactParticleTriangleSolve(Contact *c, ParticleProblem *p,double dt
double ct, cs = 0;
double vt = c->ct + p->triangles[c->o0].vt;
double vs = c->cs + p->triangles[c->o0].vs;
#if ROTATION_ENABLED
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]- p->triangles[c->o0].v[i])*t[i];
vs += (p->velocity[c->o1*3+i]- p->triangles[c->o0].v[i])*s[i];
vt += (p->velocity[c->o1*3+i])*t[i];
vs += (p->velocity[c->o1*3+i])*s[i];
#if ROTATION_ENABLED
vt += (res1[i]*p->particles[c->o1].r)*t[i];
vs += (res1[i]*p->particles[c->o1].r)*s[i];
......@@ -768,20 +793,23 @@ static int contactParticleTriangleSolve(Contact *c, ParticleProblem *p,double dt
vt *= 3.5*mu*dp/norme;
vs *= 3.5*mu*dp/norme;
}
ct = vt;
cs = vs;
if (dp == 0.){
cs = 0;
ct = 0;
}
//if(ct > 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m) ct = 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m;
//if(cs > 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m) cs = 1e-3*7e7*p->particles[c->o1].r*p->particles[c->o1].r*M_PI/p->particles[c->o1].m;
cs -= c->cs;
ct -= c->ct;
#if ROTATION_ENABLED
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]*c->a1
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i]*c->a1;
}
#else
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*c->a1, s);
......@@ -965,8 +993,8 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*(p->a/(p->a+1))*c->a1, t);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]*c->a1
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i]*c->a1;
}
#endif
#else
......@@ -1014,8 +1042,8 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*(p->a/(p->a+1))*c->a1, t);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]*c->a1
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i]*c->a1;
}
#endif
#else
......@@ -1054,8 +1082,8 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*(p->a/(p->a+1))*c->a1, t);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]*c->a1
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i]*c->a1;
}
#else
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs*c->a1, s);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment