tenseur des contraintes wip

python/scontact.py

@@ -133,6 +133,10 @@ class ParticleProblem :
         assert self._friction_enabled
         return self._get_matrix("Omega", 1)
 
+    def stressTensor(self):
+        """Return the stress tensor for post visualization."""
+        return self._get_matrix("StressTensor",self._dim*self._dim)
+
     def disk_tag(self):
         return self._get_idx("DiskTag")
 
@@ -232,7 +236,12 @@ class ParticleProblem :
         self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
         self._fc = (self.velocity()-vn)*self.mass()/dt - forces
 
-    def write_vtk(self, odir, i, t) :
+    def computeStressTensor(self, nodes, radius):
+        """Compute the stres tensor of the granular material"""
+        n_nodes = len(nodes[:,0])
+        self._lib.particleProblemComputeStressTensor(self._p,_np2c(nodes[:,:self._dim]),c_double(radius),c_int(n_nodes))
+
+    def write_vtk(self, odir, i, t, nodes=None, elements=None) :
         """Write output files for post-visualization."""
 
         v = self.velocity()
@@ -284,6 +293,17 @@ class ParticleProblem :
         VTK.write(odir+"/contacts",i,t,None,x,field_data=fdata,vtktype="vtp")
 
         VTK.write_multipart(odir+"/particle_problem",["particles_%05d.vtp"%i,"boundaries_%05d.vtu"%i,"contacts_%05d.vtp"%i],i,t)
+
+        if nodes is not None:
+            stressTensor = self.stressTensor()
+            data = [("stressTensor",stressTensor)]
+            nmat = self._lib.particleProblemNMaterial(self._p)
+            self._lib.particleProblemGetMaterialTagName.restype = c_char_p
+            tags = list([self._lib.particleProblemGetMaterialTagName(self._p,i) for i in range(nmat)])
+            connectivity = elements
+            types = np.repeat([5 if self._dim == 2 else 10],connectivity.shape[0])
+            offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0]))
+            VTK.write(odir+"/post",i,t,(types,connectivity,offsets),nodes,data,vtktype="vtu")
         
     def read_vtk(self,dirname,i,contact_factor=1) :
         """Read an existing output file to set the particle data."""

scontact/scontact.c

@@ -42,6 +42,7 @@ struct _ParticleProblem{
   Contact *contacts;
   Contact *savedContacts;
   double *position, *velocity, *externalForces;
+  double *stressTensor;
   Disk *disks;
   char **tagname, **materialName;
   int *diskTag, *segmentTag;  
@@ -1031,9 +1032,11 @@ static void _genContactByParticle(ParticleProblem *p, size_t **contactByParticle
   *contactByParticlePPtr = contactByParticleP;
 }
 
-static void _particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes, const double radius, const int n_nodes)
+void particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes, const double radius, const int n_nodes)
 {
   size_t *contactByParticle, *contactByParticleP;
+  vectorClear(p->stressTensor);
+  vectorPushN(&p->stressTensor,DIMENSION*DIMENSION*n_nodes);
 #if DIMENSION == 2
   double V = M_PI*radius*radius;
 #else
@@ -1047,9 +1050,6 @@ static void _particleProblemComputeStressTensor(ParticleProblem *p, const double
   int *found = NULL;
   vectorPushN(&found, 100);
   vectorClear(found);
-  // Particles
-  Contact *oldContacts = vectorDup(p->contacts), *cold;
-  vectorClear(p->contacts);
   int ntot = 0;
   for (size_t i = 0; i < vectorSize(p->particles); ++i) {
     particleBoundingBox(&p->particles[i], &p->position[i * DIMENSION], amin, amax);
@@ -1057,32 +1057,44 @@ static void _particleProblemComputeStressTensor(ParticleProblem *p, const double
     cellAdd(tree, amin, amax, i, &found);
   }
   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;
     }
     vectorClear(found);
-    cellSearch(tree,amin,amax,&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){
+      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];
 #if DIMENSION == 2
         double t[DIMENSION] = {-c->n[1],c->n[0]};
 #else
-        double t[DIMENSION] = {0,0,0}; //TODO when doing 3D friction
+        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]);
+          d += (xc[l]-nodes[i*DIMENSION+l])*(xc[l]-nodes[i*DIMENSION+l]); //Compute the contact point position
         }
         d = sqrt(d);
-        double F[DIMENSION][DIMENSION] = {0};
-        if (d<radius){
+        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;
@@ -1097,6 +1109,12 @@ static void _particleProblemComputeStressTensor(ParticleProblem *p, const double
         }
       }
     }
+    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]);
+      }
+    }
   }
   vectorFree(found);
   cellFree(tree);
@@ -1359,6 +1377,7 @@ int *particleProblemTriangleMaterial(ParticleProblem *p) {return p->triangleMate
 double *particleProblemTriangle(ParticleProblem *p) {return (double*)&p->triangles[0];}
 #endif
 double *particleProblemVelocity(ParticleProblem *p) {return &p->velocity[0];}
+double *particleProblemStressTensor(ParticleProblem *p) {return &p->stressTensor[0];}
 double *particleProblemDisk(ParticleProblem *p) {
   return (double*)&p->disks[0];
 }
@@ -1381,6 +1400,7 @@ void particleProblemDelete(ParticleProblem *p) {
   vectorFree(p->particleMaterial);
   vectorFree(p->velocity);
   vectorFree(p->externalForces);
+  vectorFree(p->stressTensor);
 #if FRICTION_ENABLED
   vectorFree(p->mu);
   vectorFree(p->omega);
@@ -1460,6 +1480,7 @@ ParticleProblem *particleProblemNew() {
   p->particleMaterial = NULL;
   p->position = NULL;
   p->velocity = NULL;
+  p->stressTensor = NULL;
 #if FRICTION_ENABLED
   p->friction_relaxation = 1.0;
   p->omega = NULL;
@@ -1482,6 +1503,7 @@ size_t particleProblemContactN(ParticleProblem *p, int ctype){
   size_t n = 0;
   for (size_t i = 0; i < vectorSize(p->contacts); ++i)
     if (p->contacts[i].type == ctype) n+=1;
+  printf("type=%i,n=%i,sizeContacts=%d\n",ctype,n,vectorSize(p->contacts));
   return n;
 }
 

testcases/depot-2d/depot.py

@@ -111,6 +111,8 @@ 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.write_vtk(outputdir, 0, 0, fluid.coordinates(),fluid.elements())
 
 tic = time.time()
 fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
@@ -125,7 +127,8 @@ while t < tEnd :
     # Output files writting
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
-        p.write_vtk(outputdir, ioutput, t)
+        p.computeStressTensor(fluid.coordinates(),5e-2)
+        p.write_vtk(outputdir, ioutput, t, fluid.coordinates(),fluid.elements())
         fluid.export_vtk(outputdir, t, ioutput)
     ii += 1
     print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))