Skip to content

GitLab

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

wip

parent 0253319b
Hide whitespace changes
Inline Side-by-side
python/scontact.py
View file @ b02feb06
......@@ -108,7 +108,7 @@ class ParticleProblem :
self._periodicEntitytype = np.dtype([('etag', np.int32),('edim', np.int32),('periodic_transformation', np.float64, dim)])
self._periodicSegmenttype = np.dtype([('entity_id', np.int64), ('p',np.float64,(2,dim))])
self._periodicTriangletype = np.dtype([('entity_id', np.int64),('p', np.float64,(3,dim))])
self._bodytype = np.dtype([('inverseI', np.float64),('inverseM', np.float64), ('theta', np.float64),('omega',np.float64)])
self._bodytype = np.dtype([('inverseI', np.float64),('inverseM', np.float64), ('theta', np.float64),('omega',np.float64),('position',np.float64,dim),('velocity',np.float64,dim)], align=True)
self._contacttype = np.dtype([('o', np.uint64,(2,)), ('type', np.int32), ('space',np.float64),('r', np.float64,(2,self._dim)) , ('impulse', np.float64,(self._dim,))],align=True)
self.material2id = {}
self.tag2id = {}
......@@ -140,7 +140,7 @@ class ParticleProblem :
return self._get_array("body", self._bodytype)["inverseI"][:,None]
def body_position(self):
return self._get_matrix("position", self._dim)
return self._get_array("body", self._bodytype)["position"]
def position(self):
ppos = np.zeros((self.n_particles(),self._dim))
......@@ -155,7 +155,7 @@ class ParticleProblem :
return pvel
def body_velocity(self):
return self._get_matrix("velocity", self._dim)
return self._get_array("body", self._bodytype)["velocity"]
def body_omega(self) :
d = 1 if self._dim==2 else 3
......@@ -351,6 +351,7 @@ class ParticleProblem :
i -- Number of the fiel to write
t -- Time at which the simulation is
"""
bodies = self._get_array("body",self._bodytype)
v = self.body_velocity()
x = self.body_position()
if self._dim == 2 :
......@@ -453,24 +454,28 @@ class ParticleProblem :
el = cells["connectivity"]
diskidx = np.where(cells["types"] ==1)[0]
self._lib.particle_problem_allocate_disks(self._p, c_size_t(diskidx.shape[0]))
disk = self._get_array("disk",self._disktype)
disk["x"] = x[el[offsets[diskidx]],:self._dim]
disk["v"] = 0
disk["r"] = 0
disk["material"] = cdata["Material"][diskidx][:,0]
disk["tag"] = cdata["Tag"][diskidx][:,0]
#for idx in np.where(cells["types"] == 3)[0] :
# x0 = x[el[offsets[idx]],:self._dim]
# x1 = x[el[offsets[idx]+1],:self._dim]
# self._lib.particleProblemAddBoundarySegmentTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), c_int(tags[idx,0]), c_int(materials[idx,0]))
#for idx in np.where(cells["types"] == 5)[0] :
# x0 = x[el[offsets[idx]],:self._dim]
# x1 = x[el[offsets[idx]+1],:self._dim]
# x2 = x[el[offsets[idx]+2],:self._dim]
# self._lib.particleProblemAddBoundaryTriangleTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(tags[idx,0]), c_int(materials[idx,0]))
disks = self.disks()
disks["x"] = x[el[offsets[diskidx]],:self._dim]
disks["v"] = 0
disks["r"] = 0
disks["material"] = cdata["Material"][diskidx][:,0]
disks["tag"] = cdata["Tag"][diskidx][:,0]
segidx = np.where(cells["types"] == 3)[0]
self._lib.particle_problem_allocate_segments(self._p, c_size_t(segidx.shape[0]))
segxidx = np.column_stack([el[offsets[segidx]], el[offsets[segidx]+1]])
segments = self.segments()
segments["p"] = x[el[segxidx],:self._dim]
segments["v"] = 0
segments["material"] = cdata["Material"][segidx][:,0]
segments["tag"] = cdata["Tag"][segidx][:,0]
if self._dim == 3:
trixidx = np.column_stack([el[offsets[segidx]], el[offsets[segidx]+1],el[offsets[segidx]+3]])
triangles = self.triangles()
triangles["p"] = x[el[trixidx],:self._dim]
triangles["v"] = 0
triangles["material"] = cdata["Material"][triidx][:,0]
triangles["tag"] = cdata["Tag"][triidx][:,0]
return
periodicX, periodicCells,_, periodicCellsData, periodicFieldsData = VTK.read(dirname+"/periodicBoundaries_%05d.vtu"%i)
entityTransformation = periodicFieldsData["Entity_transformation"]
entity_ids = periodicCellsData["Entity_id"].ravel()
......
scontact/scontact.c
View file @ b02feb06
......@@ -46,7 +46,6 @@ struct _ParticleProblem{
Particle *particles;
Body *bodies;
Contact *contacts;
double *position, *velocity;
Disk *disks;
Segment *segments;
PeriodicEntity *periodicEntities;
......@@ -60,21 +59,6 @@ struct _ParticleProblem{
int use_queue;
};
#if DIMENSION == 3
typedef enum {PARTICLE_PARTICLE=0, PARTICLE_DISK, PARTICLE_SEGMENT, PARTICLE_TRIANGLE} ContactType;
#else
typedef enum {PARTICLE_PARTICLE=0, PARTICLE_DISK, PARTICLE_SEGMENT} ContactType;
#endif
static double dot(const double *x, const double *y) {
#if DIMENSION == 3
return x[0] * y[0] + x[1] * y[1] + x[2] * y[2];
#else
return x[0] * y[0] + x[1] * y[1];
#endif
}
/* Particle */
struct _Particle{
double r;
size_t body;
......@@ -92,8 +76,40 @@ struct _Body{
#else
double omega;
#endif
double position[DIMENSION];
double velocity[DIMENSION];
};
typedef struct {
int material;
int tag;
}Boundary;
struct _Disk {
Boundary b;
double x[DIMENSION];
double v[DIMENSION];
double r;
};
struct _Segment{
Boundary b;
double p[2][DIMENSION];
double v[DIMENSION];
};
struct _Triangle {
Boundary b;
double p[3][DIMENSION];
double v[DIMENSION];
};
#if DIMENSION == 3
typedef enum {PARTICLE_PARTICLE=0, PARTICLE_DISK, PARTICLE_SEGMENT, PARTICLE_TRIANGLE} ContactType;
#else
typedef enum {PARTICLE_PARTICLE=0, PARTICLE_DISK, PARTICLE_SEGMENT} ContactType;
#endif
struct _Contact {
// constant over one time step
size_t o0, o1;
......@@ -110,6 +126,30 @@ struct _PeriodicEntity{
double transform[DIMENSION];
};
struct _PeriodicSegment{
size_t entity_id;
double p[2][DIMENSION];
};
#if DIMENSION == 3
struct _PeriodicTriangle {
size_t entity_id;
double p[3][DIMENSION];
};
#endif
static double dot(const double *x, const double *y) {
#if DIMENSION == 3
return x[0] * y[0] + x[1] * y[1] + x[2] * y[2];
#else
return x[0] * y[0] + x[1] * y[1];
#endif
}
/* Particle */
static void particleBoundingBox(const Particle *p, const double x[DIMENSION], double pmin[DIMENSION], double pmax[DIMENSION]){
for (int i = 0; i < DIMENSION; ++i) {
pmin[i] = x[i] - p->r;
......@@ -153,12 +193,12 @@ static void basis_from_normal_vector(double basis[DIMENSION][DIMENSION]) {
}
static void particle_get_position(const ParticleProblem *p, const Particle *particle, double x[DIMENSION]) {
size_t body = particle->body;
Body *body = &p->bodies[particle->body];
const double *xp = particle->x;
double theta = p->bodies[body].theta;
double theta = body->theta;
#if DIMENSION == 2
x[0] = p->position[body*DIMENSION+0] + cos(theta)*xp[0] - sin(theta)*xp[1];
x[1] = p->position[body*DIMENSION+1] + sin(theta)*xp[0] + cos(theta)*xp[1];
x[0] = body->position[0] + cos(theta)*xp[0] - sin(theta)*xp[1];
x[1] = body->position[1] + sin(theta)*xp[0] + cos(theta)*xp[1];
#else
// TODO
#endif
......@@ -172,6 +212,7 @@ void particle_problem_get_particles_position(ParticleProblem *p, double *x){
void particle_problem_get_particles_velocity(ParticleProblem *p, double *v){
for(int i = 0;i<vectorSize(p->particles);i++){
const Body *body = &p->bodies[p->particles[i].body];
double r[3], res[3], omega[3];
for (int d = 0;d<DIMENSION;d++){
#if DIMENSION==3
......@@ -187,7 +228,7 @@ void particle_problem_get_particles_velocity(ParticleProblem *p, double *v){
#endif
_cross(omega,r,res);
for(int d = 0;d<DIMENSION;d++){
v[i*DIMENSION + d] = p->velocity[p->particles[i].body*DIMENSION + d] +res[d];
v[i*DIMENSION + d] = body->velocity[d] + res[d];
}
}
}
......@@ -206,11 +247,11 @@ static double contact_update_particle_particle(const Contact *c, const ParticleP
for (int i = 0; i < DIMENSION; ++i)
basis[0][i] /= nnorm;
basis_from_normal_vector(basis);
const double *xbody0 = &p->position[p0->body*DIMENSION];
const double *xbody1 = &p->position[p1->body*DIMENSION];
const double *xbody0 = p->bodies[p0->body].position;
const double *xbody1 = p->bodies[p1->body].position;
for (int i = 0;i<DIMENSION;i++){
r1[i] = x1[i]-xbody1[i] - basis[0][i]*p1->r;
r0[i] = x0[i]-xbody0[i] + basis[0][i]*p0->r;
r1[i] = x1[i]-xbody1[i] - basis[0][i]*p1->r;
}
return nnorm - (p0->r + p1->r);
}
......@@ -223,7 +264,7 @@ static int contact_init_particle_particle(Contact *c, const ParticleProblem *p,
const Body *body1 = &p->bodies[p1->body];
if(p0->body==p1->body)
return 0;
if(body0->invertM == 0 && body1->invertM == 0 && body1->invertI == 0 && body1->invertI == 0)
if(body0->invertM == 0 && body1->invertM == 0 && body0->invertI == 0 && body1->invertI == 0)
return 0;
c->o0 = id0;
c->o1 = id1;
......@@ -238,18 +279,6 @@ static int contact_init_particle_particle(Contact *c, const ParticleProblem *p,
return D < alert;
}
typedef struct {
int material;
int tag;
}Boundary;
struct _Disk {
Boundary b;
double x[DIMENSION];
double v[DIMENSION];
double r;
};
static void diskBoundingBox(const Disk *d, double *pmin, double *pmax) {
const double r = fabs(d->r);
for (int i = 0; i < DIMENSION; ++i) {
......@@ -284,7 +313,7 @@ static double contact_update_disk_particle(const Contact *c, const ParticleProbl
basis[0][i] /= (nnorm == 0 ? 1 : nnorm) * (d->r < 0 ? -1 : 1);
}
basis_from_normal_vector(basis);
const double *xbody = &p->position[p1->body*DIMENSION];
const double *xbody = p->bodies[p1->body].position;
for (int i = 0;i<DIMENSION;i++){
r0[i] = 0; // todo change when disks move
r1[i] = (x[i]-xbody[i] -basis[0][i]*p1->r);
......@@ -295,8 +324,8 @@ static double contact_update_disk_particle(const Contact *c, const ParticleProbl
static int contact_init_disk_particle(Contact *c, ParticleProblem *p, size_t disk_id, size_t particle_id, double *x, double alert) {
const Disk *d = &p->disks[disk_id];
const Particle *particle = &p->particles[particle_id];
const Body *body0 = &p->bodies[particle->body];
if (body0->invertI == 0 && body0->invertM == 0)
const Body *body = &p->bodies[particle->body];
if (body->invertI == 0 && body->invertM == 0)
return 0;
c->o0 = disk_id;
c->o1 = particle_id;
......@@ -310,11 +339,6 @@ static int contact_init_disk_particle(Contact *c, ParticleProblem *p, size_t dis
return D < alert;
}
struct _PeriodicSegment{
size_t entity_id;
double p[2][DIMENSION];
};
static void periodicSegmentBoundingBox(const PeriodicSegment *s, double *pmin, double *pmax) {
for (int i = 0; i < DIMENSION; ++i) {
pmin[i] = fmin(s->p[0][i], s->p[1][i]);
......@@ -322,13 +346,6 @@ static void periodicSegmentBoundingBox(const PeriodicSegment *s, double *pmin, d
}
}
struct _Segment{
Boundary b;
double p[2][DIMENSION];
double v[DIMENSION];
};
static void coordAdd(double *x, double a, const double *y) {
x[0] += a * y[0];
x[1] += a * y[1];
......@@ -366,7 +383,7 @@ static double contact_update_segment_particle(const Contact *c, const ParticlePr
basis[0][i] /= nnorm;
}
basis_from_normal_vector(basis);
const double *xbody = &p->position[p1->body*DIMENSION];
const double *xbody = p->bodies[p1->body].position;
for (int i = 0;i<DIMENSION;i++){
r0[i] = 0; // todo change if segment moves
r1[i] = x1[i] - basis[0][i]*p1->r -xbody[i];
......@@ -379,8 +396,8 @@ static int contact_init_segment_particle(Contact *c, ParticleProblem *p, size_t
double xcheck[3];
particle_get_position(p, &p->particles[particle_id], xcheck);
const Particle *particle = &p->particles[particle_id];
const Body *body0 = &p->bodies[particle->body];
if (body0->invertI == 0 && body0->invertM == 0)
const Body *body= &p->bodies[particle->body];
if (body->invertI == 0 && body->invertM == 0)
return 0;
c->o0 = segment_id;
c->o1 = particle_id;
......@@ -396,16 +413,6 @@ static int contact_init_segment_particle(Contact *c, ParticleProblem *p, size_t
#if DIMENSION == 3
struct _Triangle {
Boundary b;
double p[3][DIMENSION];
double v[DIMENSION];
};
struct _PeriodicTriangle {
size_t entity_id;
double p[3][DIMENSION];
};
void particle_problem_add_boundary_triangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], int tag, int materialTag) {
Triangle *t = vectorPush(&p->triangles);
t->b.tag = tag;
......@@ -456,8 +463,8 @@ static int contact_init_triangle_particle(Contact *c, const ParticleProblem *p,
size_t triangle_id, size_t particle_id, double *x, double alert) {
double basis[DIMENSION][DIMENSION];
const Particle *particle = &p->particles[particle_id];
const Body *body0 = &p->bodies[particle->body];
if (body0->invertI == 0 && body0->invertM == 0)
const Body *body = &p->bodies[particle->body];
if (body->invertI == 0 && body->invertM == 0)
return 0;
const Triangle *t = &p->triangles[triangle_id];
c->o0 = triangle_id;
......@@ -510,9 +517,9 @@ static void contact_get_omega_pointers(const Contact *c, ParticleProblem *p, dou
static void contact_get_velocity_pointers(const Contact *c, ParticleProblem *p, double **v0, double **v1) {
*v1 = &p->velocity[p->particles[c->o1].body*(DIMENSION)];
*v1 = p->bodies[p->particles[c->o1].body].velocity;
if (c->type == PARTICLE_PARTICLE) {
*v0 = &p->velocity[p->particles[c->o0].body*(DIMENSION)];
*v0 = p->bodies[p->particles[c->o0].body].velocity;
}
else if (c->type == PARTICLE_DISK) {
*v0 = p->disks[c->o0].v;
......@@ -940,13 +947,14 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
double r0[DIMENSION], r1[DIMENSION];
{
const Particle *p1 = &p->particles[c->o1];
Body *body1 = &p->bodies[p1->body];
#if DIMENSION == 2
p->bodies[p1->body].theta += dt*p->bodies[p1->body].omega;
body1->theta += dt*body1->omega;
#else
//todo
#endif
for (int d = 0; d < DIMENSION; ++d) {
p->position[p1->body*DIMENSION+d] += dt*p->velocity[p1->body*DIMENSION+d];
body1->position[d] += dt*body1->velocity[d];
}
double x1[DIMENSION];
particle_get_position(p, p1, x1);
......@@ -960,22 +968,22 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
double x0[DIMENSION];
Body *body0 = &p->bodies[p0->body];
#if DIMENSION == 2
body0->theta += dt*p->bodies[p0->body].omega;
body0->theta += dt*body0->omega;
#else
//todo
#endif
for (int d = 0; d < DIMENSION; ++d) {
p->position[p0->body*DIMENSION+d] += dt*p->velocity[p0->body*DIMENSION+d];
body0->position[d] += dt*body0->velocity[d];
}
particle_get_position(p, p0, x0);
D = contact_update_particle_particle(c, p, x0, x1, basis, r0, r1);
#if DIMENSION == 2
body0->theta -= dt*p->bodies[p0->body].omega;
body0->theta -= dt*body0->omega;
#else
//todo
#endif
for (int d = 0; d < DIMENSION; ++d) {
p->position[p0->body*DIMENSION+d] -= dt*p->velocity[p0->body*DIMENSION+d];
body0->position[d] -= dt*body0->velocity[d];
}
}
else {
......@@ -994,12 +1002,12 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
if (!projection_is_on_segment(&p->segments[c->o0], x1)) contact_avoided = 1;
}
#if DIMENSION == 2
p->bodies[p1->body].theta -= dt*p->bodies[p1->body].omega;
body1->theta -= dt*body1->omega;
#else
//todo
#endif
for (int d = 0; d < DIMENSION; ++d) {
p->position[p1->body*DIMENSION+d] -= dt*p->velocity[p1->body*DIMENSION+d];
body1->position[d] -= dt*body1->velocity[d];
}
}
// allow the pre-existing (from previous time steps) interpenetrations
......@@ -1154,8 +1162,9 @@ void particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes,
for (size_t ic = 0; ic < vectorSize(p->contacts); ++ic) {
const Contact *c = p->contacts + ic;
double x[DIMENSION];
const Body *body1 = &p->bodies[p->particles[c->o1].body];
for (int d = 0; d < DIMENSION; ++d) {
x[d] = p->position[p->particles[c->o1].body*DIMENSION+d] + c->r1[d];
x[d] = body1->position[d] + c->r1[d];
}
cellAdd(tree,x,x,ic,NULL);
}
......@@ -1220,8 +1229,9 @@ void particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes,
for (size_t ifound = 0; ifound < vectorSize(found); ++ifound) {
const Contact *c = p->contacts+found[ifound];
double x[DIMENSION];
const Body *body1 = &p->bodies[p->particles[c->o1].body];
for (int d = 0; d < DIMENSION; ++d) {
x[d] = p->position[p->particles[c->o1].body*DIMENSION+d] + c->r1[d];
x[d] = body1->position[d] + c->r1[d];
}
double d2 = 0;
for (int d = 0; d < DIMENSION; ++d) d2 += (x[d]-nodes[inode*DIMENSION+d])*(x[d]-nodes[inode*DIMENSION+d]);
......@@ -1340,18 +1350,24 @@ int particle_problem_iterate(ParticleProblem *p, double alert, double dt, double
}
#endif
double *oldVelocity = malloc(sizeof(double)*vectorSize(p->velocity));
double *oldVelocity = malloc(sizeof(double)*vectorSize(p->bodies)*DIMENSION);
for (size_t j = 0; j < vectorSize(p->bodies); ++j){
for (size_t i = 0; i < DIMENSION; ++i){
oldVelocity[j*DIMENSION+i] = p->bodies[j].velocity[i];
}
}
for (size_t j = 0; j < vectorSize(p->particles); ++j){
Body *body = &p->bodies[p->particles[j].body];
for (size_t i = 0; i < DIMENSION; ++i){
oldVelocity[p->particles[j].body*DIMENSION + i] = p->velocity[p->particles[j].body * DIMENSION + i];
p->velocity[p->particles[j].body*DIMENSION+i] += externalForces[j * DIMENSION + i] * dt * p->bodies[p->particles[j].body].invertM;
body->velocity[i] += externalForces[j*DIMENSION+i]*dt*body->invertM;
//Ajouter torques
}
}
if(!particleProblemSolve(p, alert, dt, tol, maxit) && !force_motion){
for (size_t j = 0; j < vectorSize(p->bodies); ++j){
Body *body = &p->bodies[p->particles[j].body];
for (size_t i = 0; i < DIMENSION; ++i){
p->velocity[j * DIMENSION + i] = oldVelocity[j*DIMENSION + i];
body->velocity[i] = oldVelocity[j*DIMENSION + i];
}
}
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
......@@ -1363,12 +1379,13 @@ int particle_problem_iterate(ParticleProblem *p, double alert, double dt, double
return 0;
}
for (size_t i = 0; i < vectorSize(p->position); ++i){
p->position[i] += p->velocity[i] * dt;
}
for (size_t i = 0; i < vectorSize(p->bodies); ++i){
Body *body = &p->bodies[i];
for (int d = 0; d < DIMENSION; ++d) {
body->position[d] += body->velocity[d]*dt;
}
#if DIMENSION ==2
p->bodies[i].theta += p->bodies[i].omega * dt;
body->theta += body->omega*dt;
#else
//TODO
#endif
......@@ -1423,8 +1440,6 @@ size_t particle_problem_add_boundary_periodic_segment(ParticleProblem *p, const
}
void particleProblemRemoveBodies(ParticleProblem *p, const int *keep_flag) {
vectorRemoveFlag(p->position,keep_flag,DIMENSION);
vectorRemoveFlag(p->velocity,keep_flag,DIMENSION);
int *keepParticles = malloc(sizeof(size_t)*vectorSize(p->particles));
for (int i = 0;i<vectorSize(p->particles);i++){
if (keep_flag[p->particles[i].body])
......@@ -1482,11 +1497,9 @@ size_t particle_problem_add_body(ParticleProblem *p, const double x[DIMENSION],
for (int d = 0; d < 3; ++d)
body->omega[d] = 0;
#endif
vectorPushN(&p->position, DIMENSION);
vectorPushN(&p->velocity, DIMENSION);
for (int i = 0; i < DIMENSION; ++i) {
p->position[n * DIMENSION + i] = x[i];
p->velocity[n * DIMENSION + i] = 0;
body->position[i] = x[i];
body->velocity[i] = 0;
}
return n;
}
......@@ -1516,15 +1529,14 @@ void particle_problem_set_contacts(ParticleProblem *p, size_t n, const size_t *o
double *particle_problem_triangle(ParticleProblem *p) {return (double*)&p->triangles[0];}
double *particle_problem_periodic_triangle(ParticleProblem *p) {return (double*)&p->periodicTriangles[0];}
#endif
double *particle_problem_velocity(ParticleProblem *p) {return &p->velocity[0];}
double *particle_problem_disk(ParticleProblem *p) { return (double*)&p->disks[0];}
double *particle_problem_(ParticleProblem *p) { return (double*)&p->disks[0];}
double *particle_problem_periodic_entity(ParticleProblem *p) {return (double*)&p->periodicEntities[0];}
double *particle_problem_segment(ParticleProblem *p) {return (double*)&p->segments[0];}
double *particle_problem_periodic_segment(ParticleProblem *p){return (double*)&p->periodicSegments[0];}
void *particle_problem_particle(ParticleProblem *p) {return (double*)&p->particles[0];}
void *particle_problem_body(ParticleProblem *p) {return (void*)&p->bodies[0];}
void *particle_problem_contact(ParticleProblem *p) {return (void*)&p->contacts[0];}
double *particle_problem_position(ParticleProblem *p){return &p->position[0];}
double *particle_problem_mu(ParticleProblem *p){return p->mu;}
void particle_problem_allocate_disks(ParticleProblem *p, size_t n) {
vectorClear(p->disks);
......@@ -1556,10 +1568,6 @@ void particle_problem_allocate_periodic_triangles(ParticleProblem *p, size_t n)
void particle_problem_allocate_bodies(ParticleProblem *p, size_t n) {
vectorClear(p->bodies);
vectorPushN(&p->bodies,n);
vectorClear(p->velocity);
vectorPushN(&p->velocity,n*DIMENSION);
vectorClear(p->position);
vectorPushN(&p->position,n*DIMENSION);
}
void particle_problem_allocate_particles(ParticleProblem *p, size_t n) {
......@@ -1574,8 +1582,6 @@ void particle_problem_delete(ParticleProblem *p) {
vectorFree(p->periodicEntities);
vectorFree(p->periodicSegments);
vectorFree(p->contacts);
vectorFree(p->position);
vectorFree(p->velocity);
vectorFree(p->mu);
#if DIMENSION == 3
vectorFree(p->periodicTriangles);
......@@ -1600,8 +1606,6 @@ ParticleProblem *particle_problem_new() {
p->segments = NULL;
p->periodicEntities = NULL;
p->periodicSegments = NULL;
p->position = NULL;
p->velocity = NULL;
p->mu = NULL;
p->n_material = 0;
#if DIMENSION == 3
......
testcases/bodies/tetris.py
View file @ b02feb06
......@@ -53,5 +53,5 @@ while t<tend :
if i %outf == 0 :
ioutput = int(i/outf) + 1
p.write_vtk(odir, ioutput, t)
p.read_vtk(odir, ioutput, t)
#p.read_vtk(odir, ioutput, t)
i += 1
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