Commit b3b27b8c authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

lmgc particle-segments contacts

parent 7fdc442b
Pipeline #4711 passed with stage
in 22 seconds
......@@ -63,7 +63,7 @@
name="Script"
command="SetScript"
number_of_elements="1"
default_values="import numpy as np

#input
#nf = 5
#rf = 2.5e-4

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 = []
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-disk contacts
if n_disks :
 disks = in_bnd.Cells[1:2*n_disks:2]
 vals.append(in_contacts.FieldData["particle_disk"])
 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"])
 contacts = in_contacts.FieldData["particle_segment_idx"]
 points_s = np.ndarray((contacts.shape[0],2,3))
 points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
 points_s[:,1,:] = np.mean(in_bnd.Points[segments[contacts[:,0],:]],axis=1)
 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"])
 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)

#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")
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 = []
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-disk contacts
if n_disks :
 disks = in_bnd.Cells[1:2*n_disks:2]
 vals.append(in_contacts.FieldData["particle_disk"])
 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 :
 loc = in_bnd.CellLocations[in_bnd.CellTypes == 3]
 segments = np.column_stack([in_bnd.Cells[loc+1],in_bnd.Cells[loc+2]])
 loc = in_bnd.CellLocations[in_bnd.CellTypes == 7]
 segments_lmgc = np.column_stack([in_bnd.Cells[loc+2],in_bnd.Cells[loc+3]])
 segments = np.vstack([segments,segments_lmgc])
 vals.append(in_contacts.FieldData["particle_segment"])
 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"])
 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)
#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")
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"/>
......
......@@ -35,10 +35,6 @@ Properties = dict(nf = 5, rf = 2.5e-4)
def RequestData():
import numpy as np
#input
#nf = 5
#rf = 2.5e-4
in_particles,in_bnd,in_contacts = list(inputs[0])
n_disks = np.count_nonzero(in_bnd.CellTypes == 1)
......@@ -63,13 +59,23 @@ def RequestData():
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:]
loc = in_bnd.CellLocations[in_bnd.CellTypes == 3]
segments = np.column_stack([in_bnd.Cells[loc+1],in_bnd.Cells[loc+2]])
loc = in_bnd.CellLocations[in_bnd.CellTypes == 7]
segments_lmgc = np.column_stack([in_bnd.Cells[loc+2],in_bnd.Cells[loc+3]])
segments = np.vstack([segments,segments_lmgc])
vals.append(in_contacts.FieldData["particle_segment"])
contacts = in_contacts.FieldData["particle_segment_idx"]
points_s = np.ndarray((contacts.shape[0],2,3))
points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
points_s[:,1,:] = np.mean(in_bnd.Points[segments[contacts[:,0],:]],axis=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
......@@ -83,9 +89,9 @@ def RequestData():
points.append(points_t)
#merge everything
points = np.vstack(points)
vals = np.hstack(vals)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
......@@ -94,6 +100,7 @@ def RequestData():
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]
......
......@@ -62,16 +62,17 @@ def stdout_redirector(tfile):
#tfile.close()
os.close(saved_stdout_fd)
def get_contact_ids() :
def get_contact_ids(lmgctype) :
tfile = tempfile.TemporaryFile(mode='w+b')
with stdout_redirector(tfile):
chipy.DKDKx_DisplayProxTactors()
getattr(chipy,lmgctype+"_DisplayProxTactors")()
tfile.seek(0,io.SEEK_SET)
contactids = []
for l in tfile :
if l.startswith(b" RBDY2") :
if l.startswith(b" RBDY") :
w = l.split()
contactids.append((int(w[1]),int(w[6])))
p = [i for i,x in enumerate(w) if x==b"RBDY2"]
contactids.append((int(w[p[0]+1]),int(w[p[1]+1])))
tfile.close()
return contactids
###########################
......@@ -241,25 +242,28 @@ class ParticleProblem(object):
VTK.write(odir+"/boundaries",iiter,t,(types,connectivity,offsets),x)
# 2D: DKDKx_ID, DKJCx_ID, PLPLx_ID, PLJCx_ID, DKPLx_ID, DKKDx_ID
# 3D: SPSPx_ID, PRPRx_ID, PRPLx_ID, SPPLx_ID ...
if chipy.inter_handler_2D_getNb( chipy.DKDKx_ID ) > 0:
mmap = {}
for i,v in enumerate(self._d2bMap):
mmap[v[0]] = i
# coorx,coory,nx,ny,fn,ft,g
contacts = chipy.inter_handler_2D_getAll(chipy.DKDKx_ID)
contactids = get_contact_ids()
contactids = np.array(list((mmap[i[0]],mmap[i[1]]) for i in contactids))
else :
contacts = np.ndarray((0,1),np.float64)
contactids = np.ndarray((0,2),np.uint64)
print(contactids)
def get_contacts(lmgctype,map0,map1) :
if chipy.inter_handler_2D_getNb(getattr(chipy,lmgctype+"_ID")) > 0:
# coorx,coory,nx,ny,fn,ft,g
contacts = chipy.inter_handler_2D_getAll(getattr(chipy,lmgctype+"_ID"))
contactids = get_contact_ids(lmgctype)
contactids = np.array(list((map0[i[0]],map1[i[1]]) for i in contactids))
#contactids = np.array(contactids)-1
else :
contacts = np.ndarray((0,1),np.float64)
contactids = np.ndarray((0,2),np.uint64)
return contacts,contactids
fdata = []
fdata.append((b"particle_particle",contacts[:,[4]]))
fdata.append((b"particle_particle_idx",contactids))
dmap = dict((v[0],i) for i,v in enumerate(self._d2bMap))
pmap = dict((v[0],i) for i,v in enumerate(self._p2bMap))
dkdk,dkdkid = get_contacts("DKDKx",dmap,dmap)
dkpl,dkplid = get_contacts("DKPLx",pmap,dmap)
fdata.append((b"particle_particle",dkdk[:,[4]]))
fdata.append((b"particle_particle_idx",dkdkid))
fdata.append((b"particle_disk",np.ndarray((0,1),np.float64)))
fdata.append((b"particle_segment",np.ndarray((0,1),np.float64)))
fdata.append((b"particle_segment",dkpl[:,[4]]))
fdata.append((b"particle_segment_idx",dkplid))
fdata.append((b"particle_disk_idx",np.ndarray((0,2),np.uint64)))
fdata.append((b"particle_segment_idx",np.ndarray((0,2),np.uint64)))
x = np.array([[0.,0.,0.]])
i = iiter
VTK.write(odir+"/contacts",i,t,None,x,field_data=fdata,vtktype="vtp")
......
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