Commit 4a96fd99 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

added fric 3D

parent 36abc17f
Pipeline #8661 passed with stages
in 3 minutes and 8 seconds
......@@ -131,7 +131,6 @@ 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;
......@@ -208,7 +207,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);
for (int i = 0; i < DIMENSION; ++i) {
c->n[i] /= nnorm * (d->r < 0 ? -1 : 1);
c->n[i] /= (nnorm == 0 ? 1 : nnorm) * (d->r < 0 ? -1 : 1);
}
contactBuildBasis(c);
c->a0 = 0;
......@@ -287,7 +286,7 @@ static int segmentInitContact(size_t id, const Segment *s, size_t particle, cons
c->n[i] = s->p[0][i] - x[i] - t[i] / nt2 * dt;
nn2 += c->n[i] * c->n[i];
}
const double nnorm = sqrt(nn2);
const double nnorm = nn2 == 0 ? 1 : sqrt(nn2);
for (int i = 0; i < DIMENSION; ++i) {
c->n[i] /= nnorm;
}
......@@ -332,17 +331,17 @@ static void triangleBoundingBox(const Triangle *t, double *pmin, double *pmax) {
}
}
static int triangleProjectionIsInside(const void *pt, const double *x) {
static double triangleProjectionIsInside(const void *pt, const double *x, const double *n, double vnfree,double *vbnd, double dv) {
const Triangle *t = (const Triangle*)pt;
double d0[3] = {t->p[1][0] - t->p[0][0], t->p[1][1] - t->p[0][1], t->p[1][2] - t->p[0][2]};
double d1[3] = {t->p[2][0] - t->p[0][0], t->p[2][1] - t->p[0][1], t->p[2][2] - t->p[0][2]};
double n[3];
_cross(d0, d1, n);
double tn[3];
_cross(d0, d1, tn);
double xp[3];
double dp[3] = {t->p[0][0] - x[0], t->p[0][1] - x[1], t->p[0][2] - x[2]};
const double nd = dot(n, dp);
const double nd = dot(tn, dp);
for (int i = 0; i < 3; ++i)
xp[i] = x[i] + n[i] * nd;
xp[i] = x[i] + tn[i] * nd;
double dx[3][3] =
{{xp[0] - t->p[0][0], xp[1] - t->p[0][1], xp[2] - t->p[0][2]},
{xp[0] - t->p[1][0], xp[1] - t->p[1][1], xp[2] - t->p[1][2]},
......@@ -351,9 +350,15 @@ static int triangleProjectionIsInside(const void *pt, const double *x) {
for (int i = 0; i < 3; ++i) {
double d[3];
_cross(dx[(i + 1) % 3], dx[(i + 2) %3], d);
alpha[i] = 0.5 * dot(n, d);
alpha[i] = dot(tn, d)/(tn[0]*tn[0] + tn[1]*tn[1] + tn[2]*tn[2]);
}
return ((alpha[0] <= 0 && alpha[1] <= 0 && alpha[2]<= 0) || (alpha[0] >= 0 && alpha[1] >= 0 && alpha[2] >= 0));
if (!((alpha[0] <= 0 && alpha[1] <= 0 && alpha[2]<= 0) || (alpha[0] >= 0 && alpha[1] >= 0 && alpha[2] >= 0))) return 0;
for (int i = 0; i < DIMENSION; ++i) {
vbnd[i] = alpha[0]*t->v[0][i] + alpha[1]*t->v[1][i] + alpha[2]*t->v[2][i];
vnfree -= vbnd[i] * n[i];
}
if (dv > fabs(vnfree)) dv = fabs(vnfree);
return dv;
}
static int triangleInitContact(size_t id, const Triangle *t, size_t particle, const Particle *p, double *x, double alert, Contact *c) {
......@@ -368,7 +373,7 @@ static int triangleInitContact(size_t id, const Triangle *t, size_t particle, co
_cross(d0, d1, N);
const double nn = sqrt(dot(N, N));
for (int i = 0; i < 3; ++i)
c->n[i] = N[i] / nn;
c->n[i] = N[i] / (nn == 0 ? 1 : nn);
double dd[3] = {t->p[0][0] - x[0], t->p[0][1] - x[1], t->p[0][2] - x[2]};
c->D = dot(c->n, dd);
if (c->D < 0) {
......@@ -380,10 +385,8 @@ static int triangleInitContact(size_t id, const Triangle *t, size_t particle, co
c->D -= p->r;
c->a0 = 0;
c->a1 = 1;
if (c->D < 0)
c->D = 0;
c->type = PARTICLE_TRIANGLE;
return (c->D > 0 || triangleProjectionIsInside(t, x)) && c->D < alert;
return c->D < alert;
}
#endif
......
......@@ -10,10 +10,10 @@ r = 0.05
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",material="Steel")
#p.add_boundary_disk([-a/3-r,0],0,"bnd",material="Steel")
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")
......@@ -23,8 +23,8 @@ 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")
p.set_friction_coefficient(1.5,"Sand","Sand")
p.set_friction_coefficient(1.5,"Sand","Steel")
i = 0
f = np.zeros_like(p.velocity())
f[:,1] -= 0.0002
......
Markdown is supported
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