Commit 61e1b247 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

filter points near boundary for stressTensor 3D

parent a94df624
Pipeline #7362 passed with stage
in 1 minute and 50 seconds
......@@ -23,7 +23,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
nsub = max(min_nsub, int(np.ceil((vmax * dt * 8)/min(particles.r()))))
else:
nsub = 1
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(particles.r()))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(particles.r()) if particles.r() is not None else 0 )
#If time step was split too many times, ignore the convergence and move the grains
if iteration == max_split:
for i in range(nsub) :
......
......@@ -977,6 +977,13 @@ void particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes,
segmentBoundingBox(&p->segments[is], amin, amax);
cellAdd(stree,amin,amax,is,NULL);
}
#if DIMENSION == 3
Cell *ttree = cellNew(bbmin, bbmax);
for (size_t it = 0; it < vectorSize(p->triangles); ++it) {
triangleBoundingBox(&p->triangles[it], amin, amax);
cellAdd(ttree,amin,amax,it,NULL);
}
#endif
for (int inode = 0; inode < n_nodes; ++inode) {
const double *x = nodes+inode*DIMENSION;
for (int d = 0; d < DIMENSION; ++d) {
......@@ -1000,6 +1007,26 @@ void particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes,
}
}
vectorClear(found);
#if DIMENSION==3
cellSearch(ttree,amin,amax,&found);
for (size_t ifound = 0; ifound < vectorSize(found); ++ifound) {
const Triangle *t = &p->triangles[found[ifound]];
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);
const double nn = sqrt(dot(N, N));
for (int i = 0; i < 3; ++i)
N[i] /= nn;
double dd[3] = {t->p[0][0] - x[0], t->p[0][1] - x[1], t->p[0][2] - x[2]};
double dist = fabs(dot(N, dd));
if (dist < radius*0.999){
close_to_bnd = 1;
break;
}
}
vectorClear(found);
#endif
if (close_to_bnd == 1) continue;
cellSearch(tree,amin,amax,&found);
for (size_t ifound = 0; ifound < vectorSize(found); ++ifound) {
......@@ -1039,6 +1066,9 @@ void particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes,
vectorFree(found);
cellFree(tree);
cellFree(stree);
#if DIMENSION == 3
cellFree(ttree);
#endif
}
static int _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, double tol) {
......
......@@ -15,8 +15,6 @@ p.read_vtk(directory,0)
x = p.disks()["x"]
xmin,ymin,zmin = np.min(x,axis=0)
xmax,ymax,zmax = np.max(x,axis=0)
#Time between output files
timeStep = 1
z = np.arange(zmin,zmax,step)
x = np.arange(xmin,xmax,step)
y = np.arange(ymin,ymax,step)
......@@ -28,7 +26,7 @@ for filename,it,t in VTK.read_index(directory+"/particles.pvd") :
p.read_vtk(directory,it)
stress = p.computeStressTensor(points,radius)
tensor = stress[:,[0,4,8,1,5,2]]
lmax = np.amax(abs(np.linalg.eigvalsh(stress.reshape((-1,3,3))),axis=1)
lmax = np.amax(abs(np.linalg.eigvalsh(stress.reshape((-1,3,3)))),axis=1)
types = np.full([points.shape[0]],1)
location = np.arange(0,points.shape[0])
cells = np.zeros((points.shape[0],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