Skip to content

GitLab

Commit 2a0e2eeb authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

2d seems ok

parent d7e5c563
Pipeline #9308 failed with stages
in 2 minutes and 23 seconds
Hide whitespace changes
Inline Side-by-side
scontact/scontact.c
View file @ 2a0e2eeb
......@@ -110,6 +110,7 @@ struct _Contact {
double dv[DIMENSION];
double basis[DIMENSION][DIMENSION];
ContactType type;
double D0;
};
struct _PeriodicEntity{
......@@ -221,7 +222,8 @@ static int contact_init_particle_particle (Contact *c, const ParticleProblem *p,
for (int i = 0; i < DIMENSION; ++i)
c->basis[0][i] /= nnorm;
contactBuildBasis(c);
c->D = fmax(0., nnorm - (p0->r + p1->r));
//c->D = fmax(0., nnorm - (p0->r + p1->r));
c->D = nnorm - (p0->r + p1->r);
c->imass0 = p->ForcedFlag[p0->body] ? 0. : 1/p->bodies[p0->body].m;
c->imass1 = p->ForcedFlag[p1->body] ? 0. : 1/p->bodies[p1->body].m;
#if ROTATION_ENABLED
......@@ -828,11 +830,14 @@ void contact_tree_add_particle(ContactTree *tree, int id, const Contact *old_con
||(p->ForcedFlag[p->particles[c->o0].body]&&p->ForcedFlag[p->particles[c->o1].body])){
vectorPop(p->contacts);
}
else if ((cold = findContactSorted(c, old_contacts, iold))) {
for (int d = 0; d < DIMENSION; ++d) {
c->dv[d] = cold->dv[d];
else {
if ((cold = findContactSorted(c, old_contacts, iold))) {
for (int d = 0; d < DIMENSION; ++d) {
c->dv[d] = cold->dv[d];
}
contact_apply(c, p);
}
contact_apply(c, p);
c->D0 = fmin(c->D, 0);
}
}
}
......@@ -876,11 +881,14 @@ static void contact_tree_gen_boundary_contact(ContactTree *tree, ContactType typ
const Contact *cold;
if(!r || p->ForcedFlag[p->particles[c->o1].body])
vectorPop(p->contacts);
else if ((cold = findContactSorted(c, old_contacts, iold))) {
for (int d = 0; d < DIMENSION; ++d) {
c->dv[d] = cold->dv[d];
else {
if ((cold = findContactSorted(c, old_contacts, iold))) {
for (int d = 0; d < DIMENSION; ++d) {
c->dv[d] = cold->dv[d];
}
contact_apply(c, p);
}
contact_apply(c, p);
c->D0 = fmin(c->D, 0);
}
}
}
......@@ -956,6 +964,18 @@ static void fifoFree(Fifo *f) {
free(f->data);
free(f);
}
inline static int projection_is_on_segment(const Segment *s, double x[DIMENSION]) {
double alpha = 0, beta = 0;
for (int i = 0; i < DIMENSION; ++i) {
const double d = s->p[1][i] - s->p[0][i];
alpha += (x[i] - s->p[0][i]) * d;
beta += (x[i] - s->p[1][i]) * d;
}
int r = alpha >=0 && beta <=0;
return r;
}
inline static void contact_avoid(Contact *c, ParticleProblem *p, double vnfree)
{
......@@ -1034,15 +1054,89 @@ inline static void contact_avoid(Contact *c, ParticleProblem *p, double vnfree)
static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol) {
if (c->type == PARTICLE_PARTICLE && p->particles[c->o0].body==p->particles[c->o1].body)
return 1;
if (c->type == PARTICLE_PARTICLE) {
const Particle *p0 = &p->particles[c->o0];
const Particle *p1 = &p->particles[c->o1];
double x0[DIMENSION], x1[DIMENSION];
Body *body0 = &p->bodies[p0->body];
Body *body1 = &p->bodies[p1->body];
#if ROTATION_ENABLED
body0->theta += dt*p->omega[p0->body];
body1->theta += dt*p->omega[p1->body];
for (int d = 0; d < DIMENSION; ++d) {
p->position[p0->body*DIMENSION+d] += dt*p->velocity[p0->body*DIMENSION+d];
p->position[p1->body*DIMENSION+d] += dt*p->velocity[p1->body*DIMENSION+d];
}
particle_get_position(p, p0, x0);
particle_get_position(p, p1, x1);
double dv = c->dv[0];
contact_init_particle_particle(c, p, c->o0, p0, x0, c->o1, p1, x1, 0);
c->D -= c->D0;
c->dv[0] = dv;
body0->theta -= dt*p->omega[p0->body];
body1->theta -= dt*p->omega[p1->body];
for (int d = 0; d < DIMENSION; ++d) {
p->position[p0->body*DIMENSION+d] -= dt*p->velocity[p0->body*DIMENSION+d];
p->position[p1->body*DIMENSION+d] -= dt*p->velocity[p1->body*DIMENSION+d];
}
#endif
}
int contact_avoided = 0;
if (c->type == PARTICLE_SEGMENT) {
const Particle *p1 = &p->particles[c->o1];
double x1[DIMENSION];
Body *body1 = &p->bodies[p1->body];
#if ROTATION_ENABLED
body1->theta += dt*p->omega[p1->body];
for (int d = 0; d < DIMENSION; ++d) {
p->position[p1->body*DIMENSION+d] += dt*p->velocity[p1->body*DIMENSION+d];
}
particle_get_position(p, p1, x1);
double dv = c->dv[0];
contact_init_segment_particle(c, p, c->o0, c->o1, x1, 0);
if (!projection_is_on_segment(&p->segments[c->o0], x1)) contact_avoided = 1;
c->D -= c->D0;
c->dv[0] = dv;
body1->theta -= dt*p->omega[p1->body];
for (int d = 0; d < DIMENSION; ++d) {
p->position[p1->body*DIMENSION+d] -= dt*p->velocity[p1->body*DIMENSION+d];
}
#endif
}
if (c->type == PARTICLE_DISK) {
const Particle *p1 = &p->particles[c->o1];
double x1[DIMENSION];
Body *body1 = &p->bodies[p1->body];
#if ROTATION_ENABLED
body1->theta += dt*p->omega[p1->body];
for (int d = 0; d < DIMENSION; ++d) {
p->position[p1->body*DIMENSION+d] += dt*p->velocity[p1->body*DIMENSION+d];
}
particle_get_position(p, p1, x1);
double dv = c->dv[0];
contact_init_disk_particle(c, p, c->o0, c->o1, x1, 0);
c->D -= c->D0;
c->dv[0] = dv;
body1->theta -= dt*p->omega[p1->body];
for (int d = 0; d < DIMENSION; ++d) {
p->position[p1->body*DIMENSION+d] -= dt*p->velocity[p1->body*DIMENSION+d];
}
#endif
}
double *v0, *v1;
contact_get_velocity_pointers(c, p, &v0, &v1);
double v[DIMENSION]; // local free velocity in the contact basis
double dvold[DIMENSION] = {0};
coordAdd(dvold,1,c->dv);
contact_get_local_free_velocity_normal(c, p, v0, v1, v);
c->dv[0] = fmax(0., v[0] - c->D/dt);
contact_avoid(c, p, v[0]);
double dvmax = fabs(c->dv[0]-dvold[0]);
//contact_get_local_free_velocity_normal(c, p, v0, v1, v);
double incdv0 = -c->D/dt;
double newdv0 = fmax(c->dv[0] + incdv0, 0.);
if (contact_avoided) newdv0 = 0;
double dvmax = fabs(c->dv[0]-newdv0);
c->dv[0] = newdv0;
//c->dv[0] = fmax(0., v[0] - c->D/dt);
//contact_avoid(c, p, v[0]);
//double dvmax = fabs(c->dv[0]-dvold[0]);
if(c->mu<1e-10 || c->dv[0] == 0.){
for (int d = 1; d < DIMENSION; ++d) {
c->dv[d] = 0.;
......@@ -1317,9 +1411,12 @@ static int _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, dou
fifoPush(queue, contactByBody[j]);
}
}
fifoPush(queue,ic);
}
else {
activeContact[ic] = 0;
}
iter++;
activeContact[ic] = 0;
}
fifoFree(queue);
free(activeContact);
......
testcases/bodies/tetris.py 0 → 100644
View file @ 2a0e2eeb
from migflow import scontact
import random
import os
import numpy as np
os.makedirs("output",exist_ok=True)
pieces = np.array([
[[1,0,0,0],[1,1,1,0]],
[[0,0,0,1],[0,1,1,1]],
[[1,1,1,1],[0,0,0,0]],
[[0,1,1,0],[0,1,1,0]],
[[0,0,1,1],[0,1,1,0]],
[[0,1,0,0],[1,1,1,0]],
[[0,1,1,0],[0,0,1,1]]
])
def add_random_piece(p):
pmodel = random.choice(range(len(pieces)))
print(pmodel)
piece = pieces[pmodel]
y, x = np.where(piece)
R = np.repeat(r,y.size)
coord = r*2*np.column_stack((x-2.5+2*random.random(),y+20))
p.add_body(coord, R, np.pi*R**2*rho,material=str(pmodel),forced=0)
p.omega()[-1,0] = 10*np.pi*(-1+2*np.random.random())
p.velocity()[-1,1] = -0.1
n_inserted_pieces = 0
p = scontact.ParticleProblem(2,friction_enabled=False,rotation_enabled=True)
a = 1
g = np.array([0,-9.81])
r = 0.05
rho = 1000
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],[0,-0.7],"bnd",material="Steel")
p.add_boundary_segment([0,-0.7],[-a,-a],"bnd",material="Steel")
p.add_boundary_disk([0,-0.7],0,"bnd","Steel")
tend=300
#p.set_friction_coefficient(.2,"Sand","Sand")
#p.set_friction_coefficient(.5,"Sand","Steel")
i = 0
p.write_vtk("output",0,0)
dt = 0.003
outf=3
t = 0
while t<tend :
if t*4 > n_inserted_pieces:
add_random_piece(p)
n_inserted_pieces += 1
t += dt
f = g*(np.pi*p.r()**2*rho)
p.iterate(dt,f,tol=1e-8,force_motion=1)
if i %outf == 0 :
ioutput = int(i/outf) + 1
p.write_vtk("output", ioutput, t)
i += 1
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