Commit c643f52b authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

remove Contact.D

parent e78f4e83
Pipeline #9330 failed with stages
in 2 minutes and 24 seconds
......@@ -98,14 +98,15 @@ struct _Body{
struct _Contact {
size_t o0, o1;
ContactType type;
double D0;
//
double r0[DIMENSION], r1[DIMENSION]; // vector from body center to contact point
double impulse[DIMENSION]; // impulse on o1 expressed in global basis
//
double imass0, imass1;
double mu;
double D;
double iinertmom0, iinertmom1;
double r0[DIMENSION], r1[DIMENSION]; // vector from body center to contact point
double impulse[DIMENSION]; // impulse on o1 expressed in global basis
ContactType type;
double D0;
};
struct _PeriodicEntity{
......@@ -196,7 +197,7 @@ void getParticleVelocity(ParticleProblem *p, double *v){
}
}
static void contact_update_particle_particle(Contact *c, const ParticleProblem *p, const double *x0, const double *x1, double basis[DIMENSION][DIMENSION]) {
static double contact_update_particle_particle(Contact *c, const ParticleProblem *p, const double *x0, const double *x1, double basis[DIMENSION][DIMENSION]) {
const Particle *p0 = &p->particles[c->o0];
const Particle *p1 = &p->particles[c->o1];
double nnorm = 0;
......@@ -208,14 +209,13 @@ static void contact_update_particle_particle(Contact *c, const ParticleProblem *
for (int i = 0; i < DIMENSION; ++i)
basis[0][i] /= nnorm;
basis_from_normal_vector(basis);
//c->D = fmax(0., nnorm - (p0->r + p1->r));
c->D = nnorm - (p0->r + p1->r);
const double *xbody0 = &p->position[p0->body*DIMENSION];
const double *xbody1 = &p->position[p1->body*DIMENSION];
for (int i = 0;i<DIMENSION;i++){
c->r1[i] = x1[i]-xbody1[i] - basis[0][i]*p1->r;
c->r0[i] = x0[i]-xbody0[i] + basis[0][i]*p0->r;
}
return nnorm - (p0->r + p1->r);
}
......@@ -236,8 +236,9 @@ static int contact_init_particle_particle(Contact *c, const ParticleProblem *p,
c->type = PARTICLE_PARTICLE;
c->mu = particleProblemGetMu(p,p->bodyMaterial[p->particles[id0].body],p->bodyMaterial[p->particles[id1].body]);;
double basis[DIMENSION][DIMENSION];
contact_update_particle_particle(c, p, x0, x1, basis);
return c->D < alert;
double D = contact_update_particle_particle(c, p, x0, x1, basis);
c->D0 = fmin(D, 0);
return D < alert;
}
typedef struct {
......@@ -277,7 +278,7 @@ size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x[DIMENSI
return particleProblemAddBoundaryDiskTagId(p,x,r,particleProblemGetTagId(p,tag),particleProblemGetMaterialTagId(p,materialTag));
}
static void contact_update_disk_particle(Contact *c, const ParticleProblem *p, double *x, double basis[DIMENSION][DIMENSION]) {
static double contact_update_disk_particle(Contact *c, const ParticleProblem *p, double *x, double basis[DIMENSION][DIMENSION]) {
double nnorm = 0;
const Particle *p1 = &p->particles[c->o1];
const Disk *d = &p->disks[c->o0];
......@@ -286,7 +287,6 @@ static void contact_update_disk_particle(Contact *c, const ParticleProblem *p, d
nnorm += basis[0][i] * basis[0][i];
}
nnorm = sqrt(nnorm);
c->D = (nnorm - fabs(p1->r + d->r)) * (d->r < 0 ? -1 : 1);
for (int i = 0; i < DIMENSION; ++i) {
basis[0][i] /= (nnorm == 0 ? 1 : nnorm) * (d->r < 0 ? -1 : 1);
}
......@@ -296,6 +296,7 @@ static void contact_update_disk_particle(Contact *c, const ParticleProblem *p, d
c->r0[i] = 0; // todo change when disks move
c->r1[i] = (x[i]-xbody[i] -basis[0][i]*p1->r);
}
return (nnorm - fabs(p1->r + d->r)) * (d->r < 0 ? -1 : 1);
}
static int contact_init_disk_particle(Contact *c, ParticleProblem *p, size_t disk_id, size_t particle_id, double *x, double alert) {
......@@ -313,8 +314,9 @@ static int contact_init_disk_particle(Contact *c, ParticleProblem *p, size_t dis
c->type = PARTICLE_DISK;
c->mu = particleProblemGetMu(p,d->b.material,p->bodyMaterial[p->particles[particle_id].body]);
double basis[DIMENSION][DIMENSION];
contact_update_disk_particle(c, p, x, basis);
return c->D < alert;
double D = contact_update_disk_particle(c, p, x, basis);
c->D0 = fmin(D, 0);
return D < alert;
}
struct _PeriodicSegment{
......@@ -351,7 +353,7 @@ static void segmentBoundingBox(const Segment *s, double *pmin, double *pmax) {
}
}
void contact_update_segment_particle(Contact *c, const ParticleProblem *p, double *x1, double basis[DIMENSION][DIMENSION]) {
static double contact_update_segment_particle(Contact *c, const ParticleProblem *p, double *x1, double basis[DIMENSION][DIMENSION]) {
const Segment *s = &p->segments[c->o0];
const Particle *p1 = &p->particles[c->o1];
double nt2 = 0;
......@@ -373,12 +375,12 @@ void contact_update_segment_particle(Contact *c, const ParticleProblem *p, doubl
basis[0][i] /= nnorm;
}
basis_from_normal_vector(basis);
c->D = nnorm - p1->r;
const double *xbody = &p->position[p1->body*DIMENSION];
for (int i = 0;i<DIMENSION;i++){
c->r0[i] = 0; // todo change if segment moves
c->r1[i] = x1[i] - basis[0][i]*p1->r -xbody[i];
}
return nnorm - p1->r;
}
static int contact_init_segment_particle(Contact *c, ParticleProblem *p, size_t segment_id, size_t particle_id, double *x, double alert) {
......@@ -398,8 +400,9 @@ static int contact_init_segment_particle(Contact *c, ParticleProblem *p, size_t
c->type = PARTICLE_SEGMENT;
c->mu = particleProblemGetMu(p,s->b.material,p->bodyMaterial[p->particles[particle_id].body]);
double basis[DIMENSION][DIMENSION];
contact_update_segment_particle(c, p, x, basis);
return c->D < alert;
double D = contact_update_segment_particle(c, p, x, basis);
c->D0 = fmin(D, 0);
return D < alert;
}
......@@ -480,14 +483,14 @@ static int contact_init_triangle_particle(Contact *c, const ParticleProblem *p,
basis[0][i] = N[i] / (nn == 0 ? 1 : nn);
}
double dd[3] = {x[0]+particle->x[0] - t->p[0][0],x[1] + particle->x[1]- t->p[0][1],x[2] + particle->x[2]- t->p[0][2]};
c->D = dot(basis[0], dd);
if (c->D < 0) {
double D = dot(basis[0], dd);
if (D < 0) {
for (int i = 0; i <3; ++i)
basis[0][i] = -basis[0][i];
c->D = -c->D;
D = -D;
}
basis_from_normal_vector(basis);
c->D -= particle->r;
D -= particle->r;
c->imass0 = 0;
c->imass1 = 1/p->bodies[particle->body].m;
c->iinertmom0 = 0;
......@@ -498,7 +501,8 @@ static int contact_init_triangle_particle(Contact *c, const ParticleProblem *p,
}
c->type = PARTICLE_TRIANGLE;
c->mu = particleProblemGetMu(p,t->b.material,p->bodyMaterial[p->particles[particle_id].body]);
return c->D < alert;
c->D0 = fmin(D, 0);
return D < alert;
}
#endif
......@@ -787,7 +791,6 @@ void contact_tree_add_particle(ContactTree *tree, int id, const Contact *old_con
}
contact_apply(c, p);
}
c->D0 = fmin(c->D, 0);
}
}
}
......@@ -838,7 +841,6 @@ static void contact_tree_gen_boundary_contact(ContactTree *tree, ContactType typ
}
contact_apply(c, p);
}
c->D0 = fmin(c->D, 0);
}
}
}
......@@ -935,6 +937,7 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
// position prediction with the velocity taking into account the previous contact
// then initialize the contact with those positions and revert to the position in n
double D;
{
const Particle *p1 = &p->particles[c->o1];
p->bodies[p1->body].theta += dt*p->omega[p1->body];
......@@ -952,18 +955,18 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
p->position[p0->body*DIMENSION+d] += dt*p->velocity[p0->body*DIMENSION+d];
}
particle_get_position(p, p0, x0);
contact_update_particle_particle(c, p, x0, x1, basis);
D = contact_update_particle_particle(c, p, x0, x1, basis);
body0->theta -= dt*p->omega[p0->body];
for (int d = 0; d < DIMENSION; ++d) {
p->position[p0->body*DIMENSION+d] -= dt*p->velocity[p0->body*DIMENSION+d];
}
}
if (c->type == PARTICLE_SEGMENT) {
contact_update_segment_particle(c, p, x1, basis);
D = contact_update_segment_particle(c, p, x1, basis);
if (!projection_is_on_segment(&p->segments[c->o0], x1)) contact_avoided = 1;
}
if (c->type == PARTICLE_DISK) {
contact_update_disk_particle(c, p, x1, basis);
D = contact_update_disk_particle(c, p, x1, basis);
}
p->bodies[p1->body].theta -= dt*p->omega[p1->body];
for (int d = 0; d < DIMENSION; ++d) {
......@@ -971,7 +974,7 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
}
}
// allow the pre-existing (from previous time steps) interpenetrations
c->D -= c->D0;
D -= c->D0;
// un-apply the old contact
for (int d=0; d<DIMENSION; ++d) {
oldcontact.impulse[d] = -oldcontact.impulse[d];
......@@ -994,7 +997,7 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
contact_get_local_velocity(c, p, v, basis); // v[0] < 0 if particles get closer
double oldpn = c->impulse[0]*basis[0][0] + c->impulse[1]*basis[0][1];
double oldpt = c->impulse[0]*basis[1][0] + c->impulse[1]*basis[1][1];
v[0] = fmin(c->D/dt-oldpn*wnn-oldpt*wnt, 0.);
v[0] = fmin(D/dt-oldpn*wnn-oldpt*wnt, 0.);
double pstick[2] = {
-(+wtt*v[0]-wnt*v[1])/wdet,
-(-wnt*v[0]+wnn*v[1])/wdet
......
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