Commit 5a645419 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

fix

parent 1b5a5e23
Pipeline #6791 passed with stage
in 1 minute and 21 seconds
......@@ -108,12 +108,12 @@ static int particleInitContact(size_t id0, Particle *p0, double *x0, size_t id1,
c->o0 = id0;
c->o1 = id1;
for (int k = 0; k < DIMENSION; ++k) {
c->n[k] = x0[k] - x1[k];
c->n[k] = x1[k] - x0[k];
nnorm += c->n[k] * c->n[k];
}
nnorm = sqrt(nnorm);
for (int i = 0; i < DIMENSION; ++i)
c->n[i] /= -nnorm;
c->n[i] /= nnorm;
c->D = fmax(0., nnorm - (p0->r + p1->r));
c->a0 = p1->m / (p0->m + p1->m);
c->a1 = p0->m / (p0->m + p1->m);
......@@ -1034,85 +1034,72 @@ static void _genContactByParticle(ParticleProblem *p, size_t **contactByParticle
void particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes, const double radius, const int n_nodes)
{
size_t *contactByParticle, *contactByParticleP;
printf("RADIUS = %g\n",radius);
vectorClear(p->stressTensor);
vectorPushN(&p->stressTensor,DIMENSION*DIMENSION*n_nodes);
for (int i = 0; i < n_nodes*DIMENSION*DIMENSION; ++i) {
p->stressTensor[i] = 0;
}
#if DIMENSION == 2
double V = M_PI*radius*radius;
#else
double V = 4./3.*M_PI*radius*radius*radius;
#endif
_genContactByParticle(p, &contactByParticleP, &contactByParticle);
double bbmin[DIMENSION], bbmax[DIMENSION];
_particleProblemBoundingBox(p, bbmin, bbmax);
Cell *tree = cellNew(bbmin, bbmax);
double amin[DIMENSION], amax[DIMENSION];
int *found = NULL;
vectorPushN(&found, 100);
vectorClear(found);
int ntot = 0;
for (size_t i = 0; i < vectorSize(p->particles); ++i) {
particleBoundingBox(&p->particles[i], &p->position[i * DIMENSION], amin, amax);
vectorClear(found);
cellAdd(tree, amin, amax, i, &found);
for (size_t ic = 0; ic < vectorSize(p->contacts); ++ic) {
const Contact *c = p->contacts + ic;
double *xp = p->position + c->o1*DIMENSION;
double x[DIMENSION];
for (int d = 0; d < DIMENSION; ++d) {
x[d] = xp[d] - c->n[d]*p->particles[c->o1].r;
}
cellAdd(tree,x,x,ic,NULL);
}
for (int i = 0; i < n_nodes; ++i){
double F[DIMENSION][DIMENSION] = {0};
for (int j = 0; j < DIMENSION; ++j){
amin[j] = nodes[i*DIMENSION+j]-radius;
amax[j] = nodes[i*DIMENSION+j]+radius;
for (int inode = 0; inode < n_nodes; ++inode) {
for (int d = 0; d < DIMENSION; ++d) {
amin[d] = nodes[inode*DIMENSION+d] -radius;
amax[d] = nodes[inode*DIMENSION+d] +radius;
}
vectorClear(found);
cellSearch(tree,amin,amax,&found); //Use tree to find all the grains in the square of side 2*radius
//if(vectorSize(found)!=0){
//printf("x=[%g,%g]\n",nodes[i*DIMENSION],nodes[i*DIMENSION+1]);
//}
for (size_t j = 0; j < vectorSize(found); ++j){
for (int k = contactByParticleP[found[j]]; k < contactByParticleP[found[j]+1]; ++k){ //Loop over all the contacts of the grains in the square
Contact *c = &p->contacts[contactByParticle[k]];
double *xg0 = &p->position[c->o0*DIMENSION];
double *xg1 = &p->position[c->o1*DIMENSION];
double d0 = 0;
double d1 = 0;
for (int l = 0; l < DIMENSION; ++l){
d0 += (xg0[l]-nodes[i*DIMENSION+l])*(xg0[l]-nodes[i*DIMENSION+l]);
d1 += (xg1[l]-nodes[i*DIMENSION+l])*(xg1[l]-nodes[i*DIMENSION+l]);
}
if (j!=c->o0 && sqrt(d0)<radius && sqrt(d1)<radius){
continue; //If grains of the contacts are both in the reference area, just consider the contact when looping on the first object of the contact. Avoid to take twice a contact into account
}
double r = p->particles[c->o1].r;
double xc[DIMENSION];
cellSearch(tree,amin,amax,&found);
for (size_t ifound = 0; ifound < vectorSize(found); ++ifound) {
const Contact *c = p->contacts+found[ifound];
double x[DIMENSION];
double *xp = p->position + c->o1*DIMENSION;
for (int d = 0; d < DIMENSION; ++d) {
x[d] = xp[d] - c->n[d]*p->particles[c->o1].r;
}
double d2 = 0;
for (int d = 0; d < DIMENSION; ++d) d2 += (x[d]-nodes[inode*DIMENSION+d])*(x[d]-nodes[inode*DIMENSION+d]);
if (d2 > radius*radius) continue;
double *fj = p->stressTensor + DIMENSION*DIMENSION*inode;
double m = p->particles[c->o1].m*c->a1;
double alpha = 0, d = 0;
if (c->type == PARTICLE_PARTICLE) {
alpha = 1/(1/p->particles[c->o0].m+1/p->particles[c->o1].m);
d = p->particles[c->o0].r + p->particles[c->o1].r;
}
else {
alpha = p->particles[c->o1].m;
d = p->particles[c->o1].r;
}
#if DIMENSION == 2
double t[DIMENSION] = {-c->n[1],c->n[0]};
double t[DIMENSION] = {-c->n[1],c->n[0]};
#else
double t[DIMENSION] = {0,0,0}; //TODO when doing 3D friction -> Nathan
double t[DIMENSION] = {0,0,0}; //TODO when doing 3D friction -> Nathan
#endif
double d = 0;
for (int l = 0; l < DIMENSION; ++l){
xc[l] = xg1[l]-r*c->n[l];
d += (xc[l]-nodes[i*DIMENSION+l])*(xc[l]-nodes[i*DIMENSION+l]); //Compute the contact point position
}
d = sqrt(d);
if (d<radius){//If contact is in the reference volume compute the stress tensor
for (int iD = 0; iD < DIMENSION; iD++){
for (int jD = 0; jD < DIMENSION; jD++){
F[iD][jD] += c->n[iD]*c->dv*(p->particles[c->o0].m+p->particles[c->o1].m)*(xg1[jD]-xg0[jD])/V;
for (int iD = 0; iD < DIMENSION; iD++){
double fi = c->n[iD]*alpha*c->dv;
#if FRICTION_ENABLED
F[iD][jD] += c->n[iD]*c->ct*(p->particles[c->o0].m+p->particles[c->o1].m)*(xg1[jD]-xg0[jD])/V;
fi += t[iD]*alpha*c->ct;
#endif
}
}
for (int jD = 0; jD < DIMENSION; jD++){
fj[iD*DIMENSION+jD] += fi*(-d*c->n[jD])/V;
}
else{
continue;
}
}
}
for (int iD = 0; iD < DIMENSION; iD++){
for (int jD = 0; jD < DIMENSION; jD++){
p->stressTensor[i*DIMENSION*DIMENSION+iD*DIMENSION+jD] = F[iD][jD];
//printf("node %d, iD=%d, jD=%d, F=%g; value=%g\n",i,iD,jD,F[iD][jD],p->stressTensor[i*DIMENSION*DIMENSION+iD*DIMENSION+jD]);
}
}
}
......
......@@ -111,7 +111,7 @@ fluid.set_wall_boundary("Top",pressure=0)
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.export_vtk(outputdir,0,0)
p.computeStressTensor(fluid.coordinates(),100*r)
p.computeStressTensor(fluid.coordinates(),0.01)
p.write_vtk(outputdir, 0, 0, fluid.coordinates(),fluid.elements())
tic = time.time()
......@@ -127,7 +127,7 @@ while t < tEnd :
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.computeStressTensor(fluid.coordinates(),5e-2)
p.computeStressTensor(fluid.coordinates(),0.01)
p.write_vtk(outputdir, ioutput, t, fluid.coordinates(),fluid.elements())
fluid.export_vtk(outputdir, t, ioutput)
ii += 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