update/Still Bugs

python/scontact.py

@@ -66,6 +66,7 @@ class ParticleProblem :
         if ptr is None :
             return np.ndarray((0),dtype)
         size = cast(ptr-(8 if is_64bits else 4), POINTER(c_size_t)).contents.value//dtype.itemsize
+        print(size)
         buf = (size*np.ctypeslib.as_ctypes_type(dtype)).from_address(ptr)
         return np.ctypeslib.as_array(buf)
 
@@ -124,8 +125,8 @@ class ParticleProblem :
         self._segmenttype = np.dtype(bndtype+[('p',np.float64,(2,dim)),('v',np.float64,(2,dim))])
         self._triangletype = np.dtype(bndtype+[('p',np.float64,(3,dim)),('v',np.float64,(3,dim))])
         self._periodicEntitytype = np.dtype([('etag', np.int32), ('edim', np.int32), ('periodic_transformation', np.float64, dim)])
-        self._periodicSegmenttype = np.dtype([('etag', np.int32), ('periodic_transformation', np.float64, dim), ('p', np.float64, (2,dim))])
-        self._periodicTriangletype = np.dtype([('etag', np.int32), ('periodic_transformation', np.float64, dim), ('p', np.float64, (3,dim))])
+        self._periodicSegmenttype = np.dtype([('entity_id', np.int32),('p',np.float64,(2,dim))])
+        self._periodicTriangletype = np.dtype([('entity_id', np.int32),('p', np.float64,(3,dim))])
 
     def __del__(self):
         """Deletes the particle solver structure."""
@@ -216,17 +217,18 @@ class ParticleProblem :
 
     def periodicEntity(self):
         """Returns the list of periodic entity."""
-        return self._get_array("Entity", self._periodicEntitytype)
+        return self._get_array("PeriodicEntity", self._periodicEntitytype)
 
     def periodicSegments(self):
         """Returns the list of periodic segments."""
-        return self._get_array("Segment", self._periodicSegmenttype)
+        print("Periodic Segment")
+        return self._get_array("PeriodicSegment", self._periodicSegmenttype)
 
     def periodicTriangles(self):
         """Returns the list of periodic triangles."""
         if self._dim == 2:
             return np.ndarray((0),self._periodicTriangletype)
-        return self._get_array("Triangle", self._periodicTriangletype)
+        return self._get_array("PeriodicTriangle", self._periodicTriangletype)
 
     def triangles(self):
         """Returns the list of boundary triangles."""
@@ -379,9 +381,24 @@ class ParticleProblem :
         tagnames = list([self._lib.particleProblemGetTagName(self._p,i) for i in range(ntagname)])
         VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,cell_data=[("Tag",tags.reshape([-1,1])),("Material",materials.reshape([-1,1]))],field_data=[(b"TagNames",VTK.string_array_encode(tagnames))])
         
-        # periodicEntity = self.periodicEntity()
-        # periodicSegments = self.periodicSegments()
-        # periodicTriangles = self.periodicTriangles()
+        periodicEntity = self.periodicEntity()
+        print("PEntity : \n", periodicEntity)
+        periodicSegments = self.periodicSegments()
+        print("periodicSegment itemsize :", self._periodicSegmenttype.itemsize)
+        print("PSegments : \n",periodicSegments)
+        periodicTriangles = self.periodicTriangles()
+        # print(periodicTriangles)
+        periodicX = np.vstack([periodicSegments["p"].reshape([-1,self._dim]), triangles["p"].reshape([-1,self._dim])])
+        if self._dim == 2 :
+            periodicX = np.insert(x,self._dim, 0, axis=1)
+        periodicConnectivity = np.arange(0,periodicX.shape[0], dtype = np.int32)
+        periodicTypes = np.repeat(np.array([3,5],dtype=np.int32),[periodicSegments.shape[0],periodicTriangles.shape[0]])
+        periodicOffsets = np.cumsum(np.repeat([2,3],[periodicSegments.shape[0],periodicTriangles.shape[0]]),dtype=np.int32)
+        periodicEntityIds = np.array(np.hstack([periodicSegments["entity_id"],periodicTriangles["entity_id"]]))
+        # entityTransformation = np.array(periodicEntity["periodic_transformation"])
+        VTK.write(odir+"/periodicBoundaries", i, t, (periodicTypes, periodicConnectivity, periodicOffsets),
+                    periodicX, cell_data=[("Entity_id", periodicEntityIds.reshape([-1,1]))])
+
         #Contacts
         self._lib.particleProblemContactN.restype = c_size_t
         fdata = []
@@ -472,9 +489,23 @@ class ParticleProblem :
         self._lib.particleProblemAddBoundaryDisk.restype = c_size_t
         return self._lib.particleProblemAddBoundaryDisk(self._p, self._coord_type(*x0), c_double(r), tag.encode(),material.encode())
     def add_boundary_periodic_entity(self, etag, edim, transform):
+        """Adds a periodic entity.
+
+        Keyword arguments:
+        etag -- tag of the entity
+        edim -- dimension of the entity
+        transform -- tuple of the transformation to applied to the periodic entity
+        """
         self._lib.particleProblemAddBoundaryPeriodicEntity.restype = c_size_t
         return self._lib.particleProblemAddBoundaryPeriodicEntity(self._p, c_int(etag), c_int(edim), self._coord_type(*transform))
     def add_boundary_periodic_segment(self, x0, x1, etag) :
+        """Adds a boundary periodic segment.
+
+        Keyword arguments:
+        x0 -- Tuple of the coordinates of the first endpoint
+        x1 -- Tuple of the coordinates of the second endpoint
+        tag -- entity tag
+        """
         self._lib.particleProblemAddBoundaryPeriodicSegment.restype = c_size_t
         return self._lib.particleProblemAddBoundaryPeriodicSegment(self._p, self._coord_type(*x0), self._coord_type(*x1), c_int(etag))
     def add_boundary_segment(self, x0, x1, tag, material="default") :
@@ -504,6 +535,14 @@ class ParticleProblem :
         self._lib.particleProblemAddBoundaryTriangle.restype = c_size_t
         return self._lib.particleProblemAddBoundaryTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1),  self._coord_type(*x2), tag.encode(),material.encode())
     def add_boundary_periodic_triangle(self, x0, x1, x2, etag):
+        """Adds a boundary periodic triangle.
+
+        Keyword arguments:
+        x0 -- Tuple of the coordinates of the first summit 
+        x1 -- Tuple of the coordinates of the second summit 
+        x2 -- Tuple of the coordinates of the third summit 
+        etag -- tag of the periodic entity 
+        """
         if self._dim != 3 :
             raise NameError("Triangle boundaries only available in 3D")
         self._lib.particleProblemAddBoundaryPeriodicTriangle.restype = c_size_t
@@ -645,7 +684,7 @@ class ParticleProblem :
                         self.add_boundary_periodic_triangle(
                                 x[pmap[l[0]]]+shift, x[pmap[l[1]]]+shift,
                                 x[pmap[l[2]]]+shift, ptag)
-
+        # print(self.periodicSegments())
         for dim, tag in gmsh.model.getEntities(0) :
             stag = get_entity_name(dim,tag)
             if (dim,tag) in periodic_entities : continue

scontact/scontact.c

@@ -354,14 +354,14 @@ void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3],
 size_t particleProblemAddBoundaryPeriodicTriangle(ParticleProblem  *p, const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION], const int etag){
   PeriodicTriangle *t = vectorPush(&p->periodicTriangles);
   t->entity_id = -1;
-  for(size_t  ie = 0; ie < vectorSize(p->periodicEntities); ++ie){
+  for(size_t ie = 0; ie < vectorSize(p->periodicEntities); ++ie){
     PeriodicEntity *pe = &p->periodicEntities[ie];
     if(etag == pe->etag){
-      t->entity_id = ie;
+      t->entity_id = (int) ie;
       break;
     }
   }
-  if (t->entity_id == -1){printf("Entity not found in periodic triangles !\n");
+  if (t->entity_id == -1) printf("Entity not found in periodic triangles !\n");
   for(int i = 0; i < DIMENSION; ++i){
     t->p[0][i] = x0[i];
     t->p[1][i] = x1[i];
@@ -1387,7 +1387,7 @@ size_t particleProblemAddBoundaryPeriodicSegment(ParticleProblem *p, const doubl
   for(size_t ie = 0; ie < vectorSize(p->periodicEntities); ++ie){
     PeriodicEntity *pe = &p->periodicEntities[ie];
     if (etag == pe->etag){
-      ps->entity_id = ie;
+      ps->entity_id = (int)ie;
       break;
     }
   }
@@ -1563,13 +1563,19 @@ void particleProblemSetContacts(ParticleProblem *p, size_t n, const size_t *o, c
 int *particleProblemForcedFlag(ParticleProblem *p) {return p->ForcedFlag;}
 #if DIMENSION == 3
 double *particleProblemTriangle(ParticleProblem *p) {return (double*)&p->triangles[0];}
-double *particleProblemPeriodicTriangle(ParticleProblem *p) {return (double*)&p->triangles[0];}
+double *particleProblemPeriodicTriangle(ParticleProblem *p) {return (double*)&p->periodicTriangles[0];}
 #endif
 double *particleProblemVelocity(ParticleProblem *p) {return &p->velocity[0];}
 double *particleProblemDisk(ParticleProblem *p) { return (double*)&p->disks[0];}
 double *particleProblemPeriodicEntity(ParticleProblem *p) {return (double*)&p->periodicEntities[0];}
 double *particleProblemSegment(ParticleProblem *p) {return (double*)&p->segments[0];}
-double *particleProblemPeriodicSegment(ParticleProblem *p){return (double*) &p->periodicSegments[0];}
+double *particleProblemPeriodicSegment(ParticleProblem *p){
+  for(size_t i = 0; i < vectorSize(p->periodicSegments); ++i){
+    const PeriodicSegment pe = p->periodicSegments[i];
+    printf("\ni : %lu\tetag : %d\tp1 = [%f %f],\tp2 = [%f %f]",i, pe.entity_id, pe.p[0][0], pe.p[0][1], pe.p[1][0], pe.p[1][1]);
+  }
+  printf("\n");
+  return (double*)&p->periodicSegments[0];}
 double *particleProblemParticle(ParticleProblem *p) {return (double*)&p->particles[0];}
 double *particleProblemPosition(ParticleProblem *p){return &p->position[0];}
 double *particleProblemContactForces(ParticleProblem *p){return &p->contactForces[0];}

testcases/periodicTest/2D/grains.py

@@ -65,6 +65,7 @@ p.add_particle((L - 2*r,0.5),r,r**2*np.pi*rhop)
 p.velocity()[1,0] = 0.5
 
 p.write_vtk(outputdir,0,0)
+p.read_vtk(outputdir,0)
 
 #Initial time and iteration
 t = 0
@@ -75,7 +76,7 @@ tic = time.time()
 #
 # COMPUTATION LOOP
 #
-while t < tEnd :
+while t < 0 :
     # Time integration
     time_integration.iterate(None,p,dt)