Commit 61b668a0 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

wip

parent fc3f2adf
Pipeline #10141 failed with stages
in 61 minutes and 54 seconds
......@@ -63,7 +63,7 @@
name="Script"
command="SetScript"
number_of_elements="1"
default_values="import numpy as np

in_particles,in_bnd,in_periodic,in_contacts,_ = list(inputs[0])
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_triangles = 0
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"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-segments contacts

if n_segments :
 segments = in_bnd.Cells[:3*n_segments].reshape([-1,3])[:,1:]
 vals.append(in_contacts.FieldData["particle_segment"])
 print(in_contacts.FieldData["particle_segment"].shape)
 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[3*n_segments:3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
 vals.append(in_contacts.FieldData["particle_triangle"])
 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.vstack(vals)


#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*np.linalg.norm(vals[:],axis=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,nf*2,axis=0).reshape(-1,vals.shape[1]),"Reaction")
"
default_values="import numpy as np

in_particles,in_bnd,in_periodic,in_contacts,_ = list(inputs[0])
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_triangles = 0
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"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-segments contacts

if n_segments :
 segments = in_bnd.Cells[:3*n_segments].reshape([-1,3])[:,1:]
 vals.append(in_contacts.FieldData["particle_segment"])
 print(in_contacts.FieldData["particle_segment"].shape)
 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[3*n_segments:3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
 vals.append(in_contacts.FieldData["particle_triangle"])
 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.vstack(vals)


#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*np.linalg.norm(vals[:],axis=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,nf*2,axis=0).reshape(-1,vals.shape[1]),"Reaction")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int32)
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"/>
......
......@@ -97,3 +97,15 @@ def RequestData():
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,nf*2,axis=0).reshape(-1,vals.shape[1]),"Reaction")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int32)
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]))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment