Mise à jour Paraview utils to add periodic and fluidVelocity

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,23 @@ 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)
+    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 +95,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/scontact.py

@@ -409,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.

testcases/avalanch/avalanch1fluid/avalanch1fluidfriction.py

@@ -88,7 +88,7 @@ fluid.export_vtk(outputdir,0,0)
 
 tic = time.time()
 
-G = np.zeros((p.n_particles(),p.dim())
+G = np.zeros((p.n_particles(),p.dim()))
 G[:,1] = g[1]*p.mass()[:,0]
 # 
 # COMPUTATION LOOP

testcases/avalanch/avalanch1fluid/avalanch1fluidnofriction.py

@@ -113,7 +113,7 @@ fluid.set_strong_boundary("Top",1,0)
 fluid.set_strong_boundary("Left",0,0)
 fluid.set_strong_boundary("Right",0,0)
 # Set locations of the grains 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(), init=True)
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
 # Pressure initialization
 fluid.solution()[:,2] =  (H-fluid.coordinates()[:,1])*(-g[1])*rhof
 # Write output files for post-visualization.
@@ -121,7 +121,7 @@ fluid.export_vtk(outputdir,0,0)
 
 tic = time.time()
 
-G = np.zeros((p.n_particles(),p.dim())
+G = np.zeros((p.n_particles(),p.dim()))
 G[:,1] = g[1]*p.mass()[:,0]
 #
 # COMPUTATION LOOP

testcases/avalanch/avalanchnofluid/avalanchnofluidfriction.py

@@ -68,7 +68,7 @@ p.set_friction_coefficient(0.9,"Sand","Sand")#Particle-Particle
 t = 0
 ii = 0
 outf = 5
-G = np.zeros((p.n_particles(),p.dim())
+G = np.zeros((p.n_particles(),p.dim()))
 G[:,1] = g[1]*p.mass()[:,0]
 # 
 # COMPUTATION LOOP