Merge branch 'VTK_periodic' into nathan_master_release

fluid/fluid_problem.c

@@ -1742,19 +1742,6 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_p
   free(old_particle_uvw);
 }
 
-
-void fluid_problem_after_import(FluidProblem *problem)
-{
-  if(problem->mesh_tree)
-    mesh_tree_free(problem->mesh_tree);
-  problem->mesh_tree = mesh_tree_create(problem->mesh);
-  for (int i = 0; i < problem->n_particles; ++i)
-    problem->particle_element_id[i] = -1;
-  mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
-  fluid_problem_compute_node_volume(problem);
-  compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw, NULL);
-}
-
 void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact) {
   if (problem->n_particles != n) {
     problem->n_particles = n;

fluid/fluid_problem.h

@@ -138,7 +138,6 @@ void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_e
 int fluid_problem_n_mesh_boundaries(const FluidProblem *p);
 void fluid_problem_mesh_boundary_info(const FluidProblem *p, int bid, char **bname, int *bsize);
 void fluid_problem_mesh_boundary_interfaces(const FluidProblem *p, int bid, int *binterfaces);
-int *fluid_problem_parent_tag(const FluidProblem *p);
 
 
 void fluid_problem_set_stab_param(FluidProblem *p, double  stab_param);

paraview_utils/MigFlow.xml

@@ -63,7 +63,7 @@
        name="Script"
        command="SetScript"
        number_of_elements="1"
-       default_values="import numpy as np

in_particles,in_bnd,in_contacts = list(inputs[0])

n_disks = np.count_nonzero(in_bnd.CellTypes == 1)
n_segments = np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3)
n_triangles = np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)

vals = []
vals_t = []
vals_s = [] if ("particle_particle_s" in in_contacts.FieldData.keys()) else None
points = []

#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
vals_t.append(in_contacts.FieldData["particle_particle_t"])
if vals_s is not None:
    vals_s.append(in_contacts.FieldData["particle_particle_s"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
if n_disks :
    disks = in_bnd.Cells[1:2*n_disks:2]
    vals.append(in_contacts.FieldData["particle_disk"])
    vals_t.append(in_contacts.FieldData["particle_disk_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_disk_s"])
    contacts = in_contacts.FieldData["particle_disk_idx"]
    points_d = np.ndarray((contacts.shape[0],2,3))
    points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]
    points.append(points_d)

#particle-segments contacts

if n_segments :
    segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
    vals.append(in_contacts.FieldData["particle_segment"])
    vals_t.append(in_contacts.FieldData["particle_segment_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_segment_s"])
    contacts = in_contacts.FieldData["particle_segment_idx"]
    points_s = np.ndarray((contacts.shape[0],2,3))
    points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
    s = in_bnd.Points[segments[contacts[:,0],:]]
    t = s[:,1,:]-s[:,0,:]
    t /= ((t[:,0]**2+t[:,1]**2+t[:,2]**2)**0.5)[:,None]
    d = points_s[:,0,:]-s[:,0,:]
    l = d[:,0]*t[:,0]+d[:,1]*t[:,1]+d[:,2]*t[:,2]
    points_s[:,1,:] = s[:,0,:]+l[:,None]*t
    points.append(points_s)

#particle-triangle contacts
if n_triangles :
    triangles = in_bnd.Cells[2*n_disks+3*n_segments:2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
    vals.append(in_contacts.FieldData["particle_triangle"])
    vals_t.append(in_contacts.FieldData["particle_triangle_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_triangle_s"])
    contacts = in_contacts.FieldData["particle_triangle_idx"]
    points_t = np.ndarray((contacts.shape[0],2,3))
    points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)
    points.append(points_t)

#merge everything

points = np.vstack(points)
vals = np.hstack(vals)
vals_t = np.hstack(vals_t)
if vals_s is not None :
    vals_s = np.hstack(vals_s)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
ez = np.where(np.abs(t[:,1,None])>np.abs(t[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
n1 = np.cross(t,ez)
n2 = np.cross(t,n1)
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals**(1./2)*1

opoints = points[:,None,:,:] \
    +n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
    +n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
if vals_s is not None :
    output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
for i in range(nf) :
    j = (i+1)%nf
    pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4
cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
output.SetCells(types, locations, cells.reshape([-1]))
"
+       default_values="import numpy as np

in_particles,in_bnd,in_periodic,in_contacts = list(inputs[0])

n_disks = np.count_nonzero(in_bnd.CellTypes == 1)
n_segments = np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3)
n_triangles = np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)

vals = []
vals_t = []
vals_s = [] if ("particle_particle_s" in in_contacts.FieldData.keys()) else None
points = []

#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
vals_t.append(in_contacts.FieldData["particle_particle_t"])
if vals_s is not None:
    vals_s.append(in_contacts.FieldData["particle_particle_s"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
if n_disks :
    disks = in_bnd.Cells[1:2*n_disks:2]
    vals.append(in_contacts.FieldData["particle_disk"])
    vals_t.append(in_contacts.FieldData["particle_disk_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_disk_s"])
    contacts = in_contacts.FieldData["particle_disk_idx"]
    points_d = np.ndarray((contacts.shape[0],2,3))
    points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]
    points.append(points_d)

#particle-segments contacts

if n_segments :
    segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
    vals.append(in_contacts.FieldData["particle_segment"])
    vals_t.append(in_contacts.FieldData["particle_segment_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_segment_s"])
    contacts = in_contacts.FieldData["particle_segment_idx"]
    points_s = np.ndarray((contacts.shape[0],2,3))
    points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
    s = in_bnd.Points[segments[contacts[:,0],:]]
    t = s[:,1,:]-s[:,0,:]
    t /= ((t[:,0]**2+t[:,1]**2+t[:,2]**2)**0.5)[:,None]
    d = points_s[:,0,:]-s[:,0,:]
    l = d[:,0]*t[:,0]+d[:,1]*t[:,1]+d[:,2]*t[:,2]
    points_s[:,1,:] = s[:,0,:]+l[:,None]*t
    points.append(points_s)

#particle-triangle contacts
if n_triangles :
    triangles = in_bnd.Cells[2*n_disks+3*n_segments:2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
    vals.append(in_contacts.FieldData["particle_triangle"])
    vals_t.append(in_contacts.FieldData["particle_triangle_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_triangle_s"])
    contacts = in_contacts.FieldData["particle_triangle_idx"]
    points_t = np.ndarray((contacts.shape[0],2,3))
    points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)
    points.append(points_t)

#merge everything

points = np.vstack(points)
vals = np.hstack(vals)
vals_t = np.hstack(vals_t)
if vals_s is not None :
    vals_s = np.hstack(vals_s)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
ez = np.where(np.abs(t[:,1,None])>np.abs(t[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
n1 = np.cross(t,ez)
n2 = np.cross(t,n1)
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals**(1./2)*1

opoints = points[:,None,:,:] \
    +n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
    +n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
if vals_s is not None :
    output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
for i in range(nf) :
    j = (i+1)%nf
    pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4
cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
output.SetCells(types, locations, cells.reshape([-1]))
"
        panel_visibility="advanced">
        <Hints>
         <Widget type="multi_line"/>

paraview_utils/contacts-extract.py

@@ -34,12 +34,13 @@ Properties = dict(nf = 5, rf = 2.5e-4)
 
 def RequestData():
     import numpy as np
+    in_particles,in_bnd,in_periodic,in_contacts = list(inputs[0])
 
-    in_particles,in_bnd,in_contacts = list(inputs[0])
-
-    n_disks = np.count_nonzero(in_bnd.CellTypes == 1)
-    n_segments = np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3)
-    n_triangles = np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)
+    n_disks = np.count_nonzero(in_bnd.CellTypes == 1) + np.count_nonzero(in_periodic.CellTypes == 1)
+    n_segments  = (np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3) 
+                + np.count_nonzero(in_periodic.CellTypes == 7) + np.count_nonzero(in_periodic.CellTypes == 3))
+    n_triangles = (np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)
+                + np.count_nonzero(in_periodic.CellTypes == 13)+np.count_nonzero(in_periodic.CellTypes==5))
 
     vals = []
     vals_t = []
@@ -138,4 +139,4 @@ def RequestData():
     cells = np.ndarray([n,nf,5],np.uint32)
     cells[:,:,0] = 4
     cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
-    output.SetCells(types, locations, cells.reshape([-1]))
+    output.SetCells(types, locations, cells.reshape([-1]))
\ No newline at end of file

paraview_utils/show_problem.py

@@ -15,8 +15,6 @@ def showFluid(fluid, renderView, animationScene):
     fluidDisplay.Representation = 'Surface'
     ColorBy(fluidDisplay, ('POINTS', 'velocity', 'Magnitude'))
     fluidDisplay.SetScalarBarVisibility(renderView, True)
-    animationScene.GoToLast()
-    fluidDisplay.RescaleTransferFunctionToDataRange(False, True)
     RenameSource('Fluid', fluid)
     return fluidDisplay
 
@@ -26,8 +24,8 @@ def showParticles(particleProblem, renderView):
     particles.BlockIndices = [1]
 
     glyph = Glyph(Input=particles, GlyphType='Sphere')
-    glyph.GlyphType.PhiResolution = 1
-    glyph.GlyphType.ThetaResolution = 20
+    glyph.GlyphType.PhiResolution = 10
+    glyph.GlyphType.ThetaResolution = 10
     glyph.ScaleArray = ['POINTS', 'Radius']
     glyph.ScaleFactor = 2.0
     glyph.GlyphMode = 'All Points'
@@ -44,6 +42,15 @@ def showBoundaries(particleProblem, renderView):
     RenameSource('Boundaries', boundaries)
     return boundaries
 
+def showPeriodic(particleProblem, renderView):
+    periodic = ExtractBlock(particleProblem)
+    periodic.BlockIndices = [3]
+    periodicDisplay = Show(periodic, renderView)
+    periodicDisplay.AmbientColor = [1.0,0.0,0.0]
+    periodicDisplay.DiffuseColor = [1.0,0.0,0.0]
+    RenameSource('Periodic Boundaries', periodic)
+    return periodic
+
 #Show the contacts
 def showContacts(particleProblem, renderView):
     try:
@@ -61,11 +68,24 @@ def showContacts(particleProblem, renderView):
         print("MigFlow extract contacts Filter has not been charged !\nPlease see the wiki to add it in your Paraview filter.") 
         return None
 
+def showFluidVelocity(fluid, renderView):
+    fluidVelocity = Calculator(Input = fluid)
+    fluidVelocity.ResultArrayName = 'FluidVelocity'
+    fluidVelocity.Function = 'velocity/porosity'
+    RenameSource('Fluid Velocity', fluidVelocity)
+    fluidVelocityDisplay = Show(fluidVelocity, renderView)
+    ColorBy(fluidVelocityDisplay, ('POINTS', 'FluidVelocity', 'Magnitude'))
+    fluidVelocityDisplay.SetScalarBarVisibility(renderView, True)
+    Hide(fluidVelocity)
+    return fluidVelocity
+
+
+
 # Animation Scene
 animationScene = GetAnimationScene()
 animationScene.PlayMode = 'Snap To TimeSteps'
 renderView = GetActiveViewOrCreate('RenderView')
-renderView.InteractionMode = '2D'
+renderView.InteractionMode = '3D'
 
 # Read the pvd files
 particleProblemFile = str(dirname)+"/particle_problem.pvd"
@@ -76,11 +96,20 @@ if isfile(particleProblemFile):
     showContacts(particleProblem, renderView)
     showParticles(particleProblem, renderView)
     showBoundaries(particleProblem, renderView)
+    showPeriodic(particleProblem, renderView)
     RenameSource('Particle Problem', particleProblem)
 
 if isfile(fluidFile):
     fluid = PVDReader(FileName=fluidFile)
-    showFluid(fluid, renderView ,animationScene)
+    fluidDisplay = showFluid(fluid, renderView ,animationScene)
+    if isfile(particleProblemFile):
+        Hide(fluid, renderView)
+        showFluidVelocity(fluid, renderView)
+    else:
+        animationScene.GoToLast()
+        fluidDisplay.RescaleTransferFunctionToDataRange(False, True)
+
+
 
 # Show the first time step
 animationScene.GoToFirst()

python/_tools.py

@@ -21,4 +21,8 @@ def timeit(func):
         return r
     return wrapper
 
-atexit.register(print, timers)
+def timeprint(timers) :
+    if len(timers) != 0 :
+        print(timers)
+
+atexit.register(timeprint, timers)

python/fluid.py

@@ -33,7 +33,7 @@ import numpy as np
 import os
 import sys
 from . import VTK
-from ._tools import gmsh
+from ._tools import gmsh,timeit
 try :
     from .petsclsys import LinearSystem
 except :
@@ -182,7 +182,8 @@ class FluidProblem :
         Keyword argument:
         mesh_file_name -- Name of the mesh.msh file containing information about the domain
         """
-        self._lib.fluid_problem_set_mesh(self._fp, _load_msh(mesh_file_name, self._lib, self.dimension()))
+        mesh = _load_msh(mesh_file_name, self._lib, self.dimension())
+        self._lib.fluid_problem_set_mesh(self._fp, mesh)
         self.sys = None
         gmsh.model.remove()
 
@@ -325,13 +326,12 @@ class FluidProblem :
             bnds[bname.value] = nodes
         return bnds
 
-    def export_vtk(self, output_dir, t, it,stab=False) :
+    def write_vtk(self, output_dir, it ,t, stab=False):
         """Writes output file for post-visualization.
-
         Keyword arguments:
         output_dir -- Output directory
-        t -- Computation time
         it -- Computation iteration
+        t -- Computation time
         stab -- If True exports the stabilisation parametres in the output files
         """
         v = self.solution()[:,:self._dim]
@@ -345,10 +345,11 @@ class FluidProblem :
             ("velocity",v),
             ("porosity",self.porosity()),
             ("old_porosity",self.old_porosity()),
-            ("grad",da)
+            ("grad",da),
+            ("parent_node_id", self.parent_nodes())
             ]
-        cell_data.append(("curvature",self.curvature()))
         if self._n_fluids == 2 :
+            cell_data.append(("curvature",self.curvature()))
             cell_data.append(("concentration",self.concentration_dg()))
             cell_data.append(("stab",self._get_matrix("stab_param",self.n_elements(),1)))
         field_data = [(b"Boundary %s"%(name), nodes) for name,nodes in self._mesh_boundaries().items()]
@@ -357,12 +358,17 @@ class FluidProblem :
         offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
         VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
 
-    def import_vtk(self, filename) :
+    def export_vtk(self, output_dir, t, it, stab=False) :
+        print("Careful this function is deprecated... \n\tUse write_vtk instead !")
+        self.write_vtk(output_dir, it, t, stab)
+
+    def read_vtk(self, dirname, i):
         """Reads output file to reload computation.
 
         Keyword arguments:
         filename -- Name of the file to read
         """
+        filename = dirname + "/fluid_%05d.vtu"%i
         x,cells,data,cdata,fdata = VTK.read(filename)
         mesh_boundaries = {name[9:]:nodes for name,nodes in fdata.items() if name.startswith("Boundary ")}
         cbnd_names = (c_char_p*len(mesh_boundaries))(*(i.encode() for i in mesh_boundaries.keys()))
@@ -371,11 +377,14 @@ class FluidProblem :
         bnds = np.vstack(list(mesh_boundaries.values()))
         bnd_tags = np.repeat(list(range(nbnd)),list([v.shape[0] for v in mesh_boundaries.values()]))
         bnd_tags = np.require(bnd_tags,np.int32,"C")
-        self._lib.fluid_problem_set_elements(self._fp,
+        self._lib.mesh_new_from_elements.restype = c_void_p
+        _mesh = c_void_p(self._lib.mesh_new_from_elements(
                 c_int(x.shape[0]),_np2c(x,np.float64),
                 c_int(el.shape[0]),_np2c(el,np.int32),
-                c_int(bnds.shape[0]),c_void_p(bnds.ctypes.data),c_void_p(bnd_tags.ctypes.data),c_int(len(cbnd_names)),cbnd_names
-                )
+                c_int(bnds.shape[0]),_np2c(bnds,np.int32),
+                _np2c(bnd_tags,np.int32),c_int(len(cbnd_names)),cbnd_names,
+                _np2c(data["parent_node_id"],np.int32) if "parent_node_id"  in data else None))
+        self._lib.fluid_problem_set_mesh(self._fp, _mesh)
         sol = self.solution()
         sol[:,:self._dim] = data["velocity"][:,:self._dim]
         sol[:,[self._dim]] = data["pressure"]
@@ -383,9 +392,36 @@ class FluidProblem :
             self.concentration_dg()[:] = cdata["concentration"]
         self.porosity()[:] = data["porosity"]
         self.old_porosity()[:] = data["old_porosity"]
-        self._lib.fluid_problem_after_import(self._fp)
         self.sys = None
 
+
+    def import_vtk(self, filename) :
+        print("Careful this function is deprecated... \n\tUse read_vtk instead !")
+        x,cells,data,cdata,fdata = VTK.read(filename)
+        mesh_boundaries = {name[9:]:nodes for name,nodes in fdata.items() if name.startswith("Boundary ")}
+        cbnd_names = (c_char_p*len(mesh_boundaries))(*(i.encode() for i in mesh_boundaries.keys()))
+        el = cells["connectivity"].reshape([-1,self._dim+1])
+        nbnd = len(mesh_boundaries)
+        bnds = np.vstack(list(mesh_boundaries.values()))
+        bnd_tags = np.repeat(list(range(nbnd)),list([v.shape[0] for v in mesh_boundaries.values()]))
+        bnd_tags = np.require(bnd_tags,np.int32,"C")
+        self._lib.mesh_new_from_elements.restype = c_void_p
+        _mesh = c_void_p(self._lib.mesh_new_from_elements(
+                c_int(x.shape[0]),_np2c(x,np.float64),
+                c_int(el.shape[0]),_np2c(el,np.int32),
+                c_int(bnds.shape[0]),_np2c(bnds,np.int32),
+                _np2c(bnd_tags,np.int32),c_int(len(cbnd_names)),cbnd_names,
+                _np2c(data["parent_node_id"],np.int32) if "parent_node_id"  in data else None))
+        self._lib.fluid_problem_set_mesh(self._fp, _mesh)
+        sol = self.solution()
+        sol[:,:self._dim] = data["velocity"][:,:self._dim]
+        sol[:,[self._dim]] = data["pressure"]
+        if self._n_fluids == 2 :
+            self.concentration_dg()[:] = cdata["concentration"]
+        self.porosity()[:] = data["porosity"]
+        self.old_porosity()[:] = data["old_porosity"]
+        self.sys = None
+        
     def compute_node_force(self, dt) :
         """Computes the forces to apply on each grain resulting from the fluid motion.
 
@@ -407,7 +443,6 @@ class FluidProblem :
         self._lib.fluid_problem_compute_node_particle_torque(self._fp, c_double(dt), c_void_p(torques.ctypes.data))
         return torques
 
-
     def implicit_euler(self, dt, check_residual_norm=-1, reduced_gravity=0, stab_param=0.) :
         """Solves the fluid equations.
 
@@ -523,6 +558,10 @@ class FluidProblem :
         """Gives access to the coordinate of the mesh nodes."""
         return self._get_matrix("coordinates",self.n_nodes(), 3)
 
+    def parent_nodes(self):
+        """Gives access to the parent nodes of each node."""
+        return self._get_matrix("periodic_mapping", self.n_nodes(),1, typec=c_int32)
+
     def n_fluids(self) :
         """Returns the number of fluids."""
         return self._n_fluids
@@ -534,17 +573,17 @@ class FluidProblem :
     def n_elements(self):
         """Returns the number of mesh nodes."""
         self._lib.fluid_problem_n_elements.restype = c_int
-        return self._lib.fluid_problem_n_elements(self._fp);
+        return self._lib.fluid_problem_n_elements(self._fp)
 
     def n_fields(self):
         """Returns the number of fluid fields."""
         self._lib.fluid_problem_n_fields.restype = c_int
-        return self._lib.fluid_problem_n_fields(self._fp);
+        return self._lib.fluid_problem_n_fields(self._fp)
 
     def n_nodes(self):
         """Returns the number of mesh nodes."""
         self._lib.fluid_problem_n_nodes.restype = c_int
-        return self._lib.fluid_problem_n_nodes(self._fp);
+        return self._lib.fluid_problem_n_nodes(self._fp)
 
     def porosity(self):
         """Returns the porosity at nodes"""

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([('entity_id', np.int64), ('p',np.float64,(2,dim))])
+        self._periodicTriangletype = np.dtype([('entity_id', np.int64),('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("PeriodicEntity", self._periodicEntitytype)
+
+    def periodicSegments(self):
+        """Returns the list of periodic segments."""
+        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("PeriodicTriangle", self._periodicTriangletype)
+
     def triangles(self):
         """Returns the list of boundary triangles."""
         if self._dim == 2:
@@ -346,6 +363,7 @@ class ParticleProblem :
         self._lib.particleProblemGetMaterialTagName.restype = c_char_p
         tags = list([self._lib.particleProblemGetMaterialTagName(self._p,i) for i in range(nmat)])
         VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp",field_data=[(b"MaterialNames",VTK.string_array_encode(tags))])
+
         disks = self.disks()
         segments = self.segments()
         triangles = self.triangles()
@@ -361,6 +379,23 @@ 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()
+        periodicX = np.vstack([periodicSegments["p"].reshape([-1,self._dim]), periodicTriangles["p"].reshape([-1,self._dim])])
+        if self._dim == 2 :
+            periodicX = np.insert(periodicX,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]))], field_data=[(b"Entity_transformation", entityTransformation)])
+
+        #Contacts
         self._lib.particleProblemContactN.restype = c_size_t
         fdata = []
         for name in ("particle_particle","particle_disk","particle_segment","particle_triangle") :
@@ -374,7 +409,7 @@ class ParticleProblem :
             fdata.append((nameb+b"_idx",o))
         x = np.array([[0.,0.,0.]])
         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)
+        VTK.write_multipart(odir+"/particle_problem",["particles_%05d.vtp"%i,"boundaries_%05d.vtu"%i,"periodicBoundaries_%05d.vtu"%i,"contacts_%05d.vtp"%i],i,t)
 
     def read_vtk(self,dirname,i,contact_factor=1) :
         """Reads an existing output file to set the particle data.
@@ -394,14 +429,17 @@ class ParticleProblem :
         self._lib.particleProblemClearParticles(self._p)
         self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
                 _np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(partmat,dtype=np.int32),_np2c(forced,dtype=np.int32),_np2c(d["Contacts"][:,:self._dim]*contact_factor))
+       
         x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
         tagnames =  VTK.string_array_decode(fdata["TagNames"])
         tagmap = {}
+        print(fdata)
+        print(tagnames)
         for j,n in enumerate(tagnames) :
             tagmap[j] = self._lib.particleProblemGetTagId(self._p,n)
-
         offsets = np.hstack([[0],cells["offsets"]])
         el = cells["connectivity"]
+        print(tagmap)
         tags = np.vectorize(tagmap.get)(cdata["Tag"])
         materials = cdata["Material"]
         for idx in np.where(cells["types"] == 1)[0] :
@@ -416,6 +454,32 @@ class ParticleProblem :
             x1 = x[el[offsets[idx]+1],:self._dim]
             x2 = x[el[offsets[idx]+2],:self._dim]
             self._lib.particleProblemAddBoundaryTriangleTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(tags[idx,0]), c_int(materials[idx,0]))
+        
+        periodicX, periodicCells,_, periodicCellsData, periodicFieldsData = VTK.read(dirname+"/periodicBoundaries_%05d.vtu"%i)
+        entityTransformation = periodicFieldsData["Entity_transformation"]
+        entity_ids = periodicCellsData["Entity_id"].ravel()
+        periodicOffsets = np.hstack([[0], periodicCells["offsets"]])
+        periodicEl = periodicCells["connectivity"]
+        # add entity
+        for entity_id in range(np.shape(entityTransformation)[0]):
+            trans = tuple(entityTransformation[entity_id])
+            self.add_boundary_periodic_entity(entity_id, self._dim-1, trans)
+        # add periodic segment
+        for idx in np.where(periodicCells["types"] == 3)[0]:
+            pidx = periodicEl[periodicOffsets[idx]]
+            entity_id = entity_ids[pidx//2]
+            x0 = periodicX[pidx,:self._dim]
+            x1 = periodicX[pidx+1,:self._dim]
+            self.add_boundary_periodic_segment(tuple(x0), tuple(x1),entity_id)
+        # add periodic triangles
+        for idx in np.where(periodicCells["types"] == 5)[0]:
+            pidx = periodicEl[periodicOffsets[idx]]
+            entity_id = entity_ids[pidx//3]
+            x0 = periodicX[pidx,:self._dim]
+            x1 = periodicX[pidx+1,:self._dim]
+            x2 = periodicX[pidx+2,:self._dim]
+            self.add_boundary_periodic_triangle(tuple(x0), tuple(x1), tuple(x2), entity_id)
+
         _,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
         ks = ("particle_particle","particle_disk","particle_segment","particle_triangle")
         v = []
@@ -450,9 +514,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") :
@@ -482,6 +560,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
@@ -626,7 +712,6 @@ class ParticleProblem :
                         self.add_boundary_periodic_triangle(
                                 x[pmap[l[0]]]+shift, x[pmap[l[1]]]+shift,
                                 x[pmap[l[2]]]+shift, ptag)
-
         for dim, tag in gmsh.model.getEntities(0) :
             stag = get_entity_name(dim,tag)
             if (dim,tag) in periodic_entities : continue

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;
+  size_t entity_id;
   double p[2][DIMENSION];
 };
 
@@ -331,8 +330,7 @@ struct _Triangle {
   double v[3][DIMENSION];
 };
 struct _PeriodicTriangle {
-  int etag;
-  double* transform;
+  size_t 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;
-  for(size_t  ie = 0; ie < vectorSize(p->periodicEntities); ++ie){
+  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];