fix 3d contacts visu

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 = []
points = []

#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
vals_t.append(in_contacts.FieldData["particle_particle_t"])
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"])
    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"])
    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"])
    vals.append(in_contacts.FieldData["particle_triangle_t"])
    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)
#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")
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_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

@@ -43,11 +43,14 @@ def RequestData():
 
     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
@@ -55,6 +58,8 @@ def RequestData():
         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],:]
@@ -67,6 +72,8 @@ def RequestData():
         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],:]
@@ -80,10 +87,11 @@ def RequestData():
 
     #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:]
+        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"])
-        vals.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],:]
@@ -95,6 +103,8 @@ def RequestData():
     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)
@@ -110,6 +120,8 @@ def RequestData():
     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) :