entity periodic

python/scontact.py

@@ -123,6 +123,9 @@ class ParticleProblem :
         self._disktype = np.dtype(bndtype+[('x',np.float64,dim),('v',np.float64,dim),('r',np.float64)])
         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))])
 
     def __del__(self):
         """Deletes the particle solver structure."""
@@ -211,6 +214,20 @@ class ParticleProblem :
         """Returns the list of boundary segments."""
         return self._get_array("Segment",self._segmenttype)
 
+    def periodicEntity(self):
+        """Returns the list of periodic entity."""
+        return self._get_array("Entity", self._periodicEntitytype)
+
+    def periodicSegments(self):
+        """Returns the list of periodic segments."""
+        return self._get_array("Segment", 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)
+
     def triangles(self):
         """Returns the list of boundary triangles."""
         if self._dim == 2:
@@ -361,6 +378,11 @@ class ParticleProblem :
         self._lib.particleProblemGetTagName.restype = c_char_p
         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()
+        #Contacts
         self._lib.particleProblemContactN.restype = c_size_t
         fdata = []
         for name in ("particle_particle","particle_disk","particle_segment","particle_triangle") :

scontact/scontact.c

@@ -51,7 +51,7 @@ struct _ParticleProblem{
   int *particleMaterial;
   int *ForcedFlag;
   Segment *segments;
-  PeriodicEntity* periodicEntities;
+  PeriodicEntity *periodicEntities;
   PeriodicSegment *periodicSegments;
 #if FRICTION_ENABLED
 #if ROTATION_ENABLED
@@ -228,8 +228,7 @@ static int diskInitContact(size_t id, const Disk *d, size_t particle, const Part
   return c->D < alert;
 }
 struct _PeriodicSegment{
-  int etag;
-  double* transform;
+  int entity_id;
   double p[2][DIMENSION];
 };
 
@@ -331,8 +330,7 @@ struct _Triangle {
   double v[3][DIMENSION];
 };
 struct _PeriodicTriangle {
-  int etag;
-  double* transform;
+  int entity_id;
   double p[3][DIMENSION];
 };
 
@@ -355,14 +353,15 @@ 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->etag = etag;
+  t->entity_id = -1;
   for(size_t  ie = 0; ie < vectorSize(p->periodicEntities); ++ie){
     PeriodicEntity *pe = &p->periodicEntities[ie];
-    if(t->etag == pe->etag){
-      t->transform = &pe->transform[0];
+    if(etag == pe->etag){
+      t->entity_id = ie;
       break;
     }
   }
+  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];
@@ -801,9 +800,9 @@ void contact_tree_add_periodic_segment(ContactTree *tree, PeriodicSegment *segme
   *vectorPush(&tree->id) = vectorSize(tree->periodic_tag);
   *vectorPush(&tree->type) = CPERIODIC;
   for (int d = 0; d < DIMENSION; ++d) {
-    *vectorPush(&tree->periodic_translate) = *(segment->transform + d);
+    *vectorPush(&tree->periodic_translate) = (tree->problem->periodicEntities[segment->entity_id]).transform[d];
   }
-  *vectorPush(&tree->periodic_tag) = segment->etag;
+  *vectorPush(&tree->periodic_tag) = segment->entity_id;
   double amin[DIMENSION], amax[DIMENSION];
   periodicSegmentBoundingBox(segment, amin, amax);
   addAlert(tree->alert/2, amin, amax);
@@ -815,9 +814,9 @@ void contact_tree_add_periodic_triangle(ContactTree *tree, PeriodicTriangle *tri
   *vectorPush(&tree->id) = vectorSize(tree->periodic_tag);
   *vectorPush(&tree->type) = CPERIODIC;
   for (int d = 0; d < DIMENSION; ++d) {
-    *vectorPush(&tree->periodic_translate) = *(triangle->transform + d);
+    *vectorPush(&tree->periodic_translate) = (tree->problem->periodicEntities[triangle->entity_id]).transform[d];
   }
-  *vectorPush(&tree->periodic_tag) = triangle->etag;
+  *vectorPush(&tree->periodic_tag) = triangle->entity_id;
   double amin[DIMENSION], amax[DIMENSION];
   periodicTriangleBoundingBox(triangle, amin, amax);
   addAlert(tree->alert/2, amin, amax);
@@ -1384,14 +1383,15 @@ size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIM
 }
 size_t particleProblemAddBoundaryPeriodicSegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const int etag){
   PeriodicSegment *ps = vectorPush(&p->periodicSegments);
-  ps->etag = etag;
+  ps->entity_id = -1;
   for(size_t ie = 0; ie < vectorSize(p->periodicEntities); ++ie){
     PeriodicEntity *pe = &p->periodicEntities[ie];
-    if (ps->etag == pe->etag){
-      ps->transform = &pe->transform[0];
+    if (etag == pe->etag){
+      ps->entity_id = ie;
       break;
     }
   }
+  if (ps->entity_id == -1) printf("There is no entity found in periodic segment!");
   for(int i = 0; i < DIMENSION; ++i){
     ps->p[0][i] = x0[i];
     ps->p[1][i] = x1[i];
@@ -1563,9 +1563,11 @@ 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];}
 #endif
 double *particleProblemVelocity(ParticleProblem *p) {return &p->velocity[0];}
-double *particleProblemDisk(ParticleProblem *p) { return (double*)&p->disks[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 *particleProblemParticle(ParticleProblem *p) {return (double*)&p->particles[0];}

scontact/scontact.h

@@ -46,7 +46,10 @@ const char* particleProblemGetMaterialTagName(ParticleProblem *p, size_t tag);
 void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION], const char *tag, const char *material);
 #endif
 double *particleProblemDisk(ParticleProblem *p);
+double *particleProblemPeriodicEntity(ParticleProblem *p);
+double *particleProblemTriangle(ParticleProblem *p);
 double *particleProblemSegment(ParticleProblem *p);
+double *particleProblemPeriodicSegment(ParticleProblem *p);
 double *particleProblemParticle(ParticleProblem *p);
 double *particleProblemVelocity(ParticleProblem *p);
 double *particleProblemPosition(ParticleProblem *p);
@@ -65,6 +68,8 @@ int *particleProblemForcedFlag(ParticleProblem *p);
 int *particleProblemTriangleTag(ParticleProblem *p);
 int *particleProblemTriangleMaterial(ParticleProblem *p);
 double *particleProblemTriangle(ParticleProblem *p);
+double *particleProblemPeriodicTriangle(ParticleProblem *p);
+
 #endif
 #if FRICTION_ENABLED
 void particleProblemSetFrictionCoefficient(ParticleProblem *p, double mu, const char *mat0, const char *mat1);

validation/periodic3d/periodic3d.py

@@ -102,7 +102,7 @@ class Periodic2d(unittest.TestCase):
         # COMPUTATION LOOP
         #
         while t < tEnd :
-            time_integration.iterate(fluid, p, dt,min_nsub=5, contact_tol = 1e-8,external_particles_forces=G)
+            time_integration.iterate(None, p, dt,min_nsub=5, contact_tol = 1e-8,external_particles_forces=G)
 
             p.position()[:,0] = np.fmod(p.position()[:,0],L)
             ind = np.where(p.position()[:,0] <= 0)