Commit 760da152 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

remove disks

parent 0fac4d43
Pipeline #9341 failed with stages
in 2 minutes and 21 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_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],:]
 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)

#generate tubes
points = np.vstack(points)
b0 = points[:,0,:]-points[:,1,:]
b0 /= np.sqrt(b0[:,0,None]**2+b0[:,1,None]**2+b0[:,2,None]**2)
ez = np.where(np.abs(b0[:,1,None])>np.abs(b0[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
b1 = np.cross(b0,ez)
b2 = np.cross(b0,b1)
#merge everything
vals = np.vstack(vals)
if vals.shape[1] == 3:
 vals_n = np.abs(vals[:,0]*b0[:,0]+vals[:,1]*b0[:,1]+vals[:,2]*b0[:,2])
 vals_t = np.hypot(vals[:,0]*b1[:,0]+vals[:,1]*b1[:,1]+vals[:,2]*b1[:,2],
 vals[:,0]*b2[:,0]+vals[:,1]*b2[:,1]+vals[:,2]*b2[:,2])
else:
 vals_n = np.abs(vals[:,0]*b0[:,0]+vals[:,1]*b0[:,1])
 vals_t = np.abs(vals[:,0]*b1[:,0]+vals[:,1]*b1[:,1])
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals_n**(1./2)*1

opoints = points[:,None,:,:] \
 +b1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
 +b2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals_n,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
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]))
"
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)

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-segments contacts
if n_segments :
 segments = in_bnd.Cells[0: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],:]
 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)

#generate tubes
points = np.vstack(points)
b0 = points[:,0,:]-points[:,1,:]
b0 /= np.sqrt(b0[:,0,None]**2+b0[:,1,None]**2+b0[:,2,None]**2)
ez = np.where(np.abs(b0[:,1,None])>np.abs(b0[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
b1 = np.cross(b0,ez)
b2 = np.cross(b0,b1)
#merge everything
vals = np.vstack(vals)
if vals.shape[1] == 3:
 vals_n = np.abs(vals[:,0]*b0[:,0]+vals[:,1]*b0[:,1]+vals[:,2]*b0[:,2])
 vals_t = np.hypot(vals[:,0]*b1[:,0]+vals[:,1]*b1[:,1]+vals[:,2]*b1[:,2],
 vals[:,0]*b2[:,0]+vals[:,1]*b2[:,1]+vals[:,2]*b2[:,2])
else:
 vals_n = np.abs(vals[:,0]*b0[:,0]+vals[:,1]*b0[:,1])
 vals_t = np.abs(vals[:,0]*b1[:,0]+vals[:,1]*b1[:,1])
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals_n**(1./2)*1

opoints = points[:,None,:,:] \
 +b1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
 +b2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals_n,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
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"/>
......
......@@ -37,7 +37,6 @@ def RequestData():
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)
......@@ -48,20 +47,10 @@ def RequestData():
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:]
segments = in_bnd.Cells[0: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))
......@@ -76,7 +65,7 @@ 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[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))
......
......@@ -102,7 +102,6 @@ class ParticleProblem :
raise NameError("Cannot create particle problem.")
bndtype =[('material',np.int32),('tag',np.int32),('body',np.uint64)]
self._particletype = np.dtype([('r',np.float64),('body',np.uint64),('x',np.float64,dim),('material',np.int32)], align=True)
self._disktype = np.dtype(bndtype+[('x',np.float64,(1,dim)),('r',np.float64)])
self._segmenttype = np.dtype(bndtype+[('x',np.float64,(2,dim))])
self._triangletype = np.dtype(bndtype+[('x',np.float64,(3,dim))])
self._periodicEntitytype = np.dtype([('etag', np.int32),('edim', np.int32),('periodic_transformation', np.float64, dim)])
......@@ -168,7 +167,6 @@ class ParticleProblem :
self._saved_velocity = np.copy(self.body_velocity())
self._saved_position = np.copy(self.body_position())
self._saved_segments = np.copy(self.segments())
self._saved_disk = np.copy(self.disks())
if self.dim() == 3 :
self._saved_triangles = np.copy(self.triangles())
self._saved_omega = np.copy(self.body_omega())
......@@ -180,7 +178,6 @@ class ParticleProblem :
self.body_omega()[:] = self._saved_omega
self._bodies()["theta"] = _self.saved_theta
self.segments()[:] = self._saved_segments
self.disks()[:] = self._saved_disk
if self.dim() == 3 :
self.triangles()[:] = self._saved_triangles
......@@ -188,10 +185,6 @@ class ParticleProblem :
"""Returns the list of particle materials."""
return self._get_array("particle", self._particletype)["material"].reshape(-1,1)
def disks(self) :
"""Returns the list of boundary disks."""
return self._get_array("disk",self._disktype)
def get_tag_id(self, tag="default"):
"""Returns the number associated to a string tag."""
if tag not in self.tag2id.keys():
......@@ -383,7 +376,6 @@ class ParticleProblem :
tags = list (s.encode() for s in self.id2material)
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()
def absolute_coord(o):
......@@ -394,19 +386,18 @@ class ParticleProblem :
xo[:,:,0] = xb[:,:,0]+o["x"][:,:,0]*ct-o["x"][:,:,1]*st
xo[:,:,1] = xb[:,:,1]+o["x"][:,:,0]*st+o["x"][:,:,1]*ct
return xo
xdisks = absolute_coord(disks).reshape([-1,self._dim])
xsegs = absolute_coord(segments).reshape([-1,self._dim])
xtris = absolute_coord(triangles).reshape([-1,self._dim])
x = np.vstack([xdisks,xsegs,xtris])
xrel = np.vstack([disks["x"].reshape(-1,self._dim),segments["x"].reshape(-1,self._dim),triangles["x"].reshape(-1,self._dim)])
x = np.vstack([xsegs,xtris])
xrel = np.vstack([segments["x"].reshape(-1,self._dim),triangles["x"].reshape(-1,self._dim)])
if self._dim == 2 :
x = np.insert(x,self._dim,0,axis=1)
xrel = np.insert(xrel,self._dim,0,axis=1)
connectivity = np.arange(0,x.shape[0],dtype=np.int32)
types = np.repeat(np.array([1,3,5],dtype=np.int32),[disks.shape[0],segments.shape[0],triangles.shape[0]])
offsets = np.cumsum(np.repeat([1,2,3],[disks.shape[0],segments.shape[0],triangles.shape[0]]),dtype=np.int32)
tags = np.array(np.hstack([disks["tag"],segments["tag"],triangles["tag"]]))
materials = np.array(np.hstack([disks["material"],segments["material"],triangles["material"]]))
types = np.repeat(np.array([3,5],dtype=np.int32),[segments.shape[0],triangles.shape[0]])
offsets = np.cumsum(np.repeat([2,3],[segments.shape[0],triangles.shape[0]]),dtype=np.int32)
tags = np.array(np.hstack([segments["tag"],triangles["tag"]]))
materials = np.array(np.hstack([segments["material"],triangles["material"]]))
tagnames = list(s.encode() for s in self.id2tag)
VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,data=[("RelativePosition",xrel)],cell_data=[("Tag",tags.reshape([-1,1])),("Material",materials.reshape([-1,1]))],field_data=[(b"TagNames",VTK.string_array_encode(tagnames))])
periodicEntity = self.periodic_entity()
......@@ -427,7 +418,7 @@ class ParticleProblem :
##Contacts
fdata = []
contacts = self.contacts()
for tid,name in enumerate(("particle_particle","particle_disk","particle_segment","particle_triangle")) :
for tid,name in enumerate(("particle_particle","particle_segment","particle_triangle")) :
o = contacts[contacts["type"]==tid]["o"]
v = contacts[contacts["type"]==tid]["reaction"]
nameb = name.encode()
......@@ -467,14 +458,6 @@ class ParticleProblem :
self.tag2id = dict((name,i) for i,name in enumerate(self.id2tag))
offsets = np.hstack([[0],cells["offsets"]])
el = cells["connectivity"]
diskidx = np.where(cells["types"] ==1)[0]
self._lib.particle_problem_allocate_disks(self._p, c_size_t(diskidx.shape[0]))
disks = self.disks()
disks["x"] = x[el[offsets[diskidx]],:self._dim]
disks["v"] = 0
disks["r"] = 0
disks["material"] = cdata["Material"][diskidx][:,0]
disks["tag"] = cdata["Tag"][diskidx][:,0]
segidx = np.where(cells["types"] == 3)[0]
self._lib.particle_problem_allocate_segments(self._p, c_size_t(segidx.shape[0]))
segxidx = np.column_stack([el[offsets[segidx]], el[offsets[segidx]+1]])
......@@ -517,7 +500,7 @@ class ParticleProblem :
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")
ks = ("particle_particle","particle_segment","particle_triangle")
v = np.vstack(list([fdata[k] for k in ks]))
o = np.vstack([fdata[k+"_idx"] for k in ks])
t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks])
......@@ -526,20 +509,6 @@ class ParticleProblem :
self.contacts()["o"] = o
self.contacts()["type"] = t
def add_boundary_disk(self, x0, r, tag, material="default") :
"""Adds a boundary disk.
Keyword arguments:
x0 -- Tuple of the coordinates of the centre
r -- Disk radius
tag -- Disk tag
material -- Disk material
"""
self._lib.particle_problem_add_boundary_disk.restype = c_size_t
imat = self.get_material_id(material)
itag = self.get_tag_id(material)
return self._lib.particle_problem_add_boundary_disk(self._p, self._coord_type(*x0), c_double(r), c_int(itag), c_int(imat))
def add_boundary_periodic_entity(self, etag, edim, transform):
"""Adds a periodic entity.
......@@ -607,6 +576,9 @@ class ParticleProblem :
self._lib.particle_problem_add_boundary_periodic_triangle.restype = c_size_t
return self._lib.particle_problem_add_boundary_periodic_triangle(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(etag))
def add_body(self, x, invertM, invertI):
self._lib.particle_problem_add_body.restype = c_size_t
return self._lib.particle_problem_add_body(self._p, self._coord_type(*x), c_double(invertM), c_double(invertI))
def add_particle_body(self, positions,radii,masses,material="default",forced=False):
"""Adds a body in the particle problem.
......@@ -625,8 +597,7 @@ class ParticleProblem :
I = np.sum(masses.reshape((-1,1))*radii.reshape((-1,1))**2/2 + np.sum((positions-x)**2,axis=1)*masses.reshape((-1,1)))
invertI = 0 if forced else 1/I
invertM = 0 if forced else 1/np.sum(masses)
self._lib.particle_problem_add_body.restype = c_size_t
idd = self._lib.particle_problem_add_body(self._p, self._coord_type(*x), c_double(invertM), c_double(invertI))
idd = self.add_body(x, invertM, invertI)
for i in range(positions.shape[0]):
self.add_particle(positions[i,:]-x,radii[i],idd, material)
......@@ -687,6 +658,8 @@ class ParticleProblem :
periodic_entities = set()
_, x, _ = gmsh.model.mesh.getNodes()
x = x.reshape([-1,3])[:,:self._dim]
self._lib.particle_problem_add_body.restype = c_size_t
bndbody = self._lib.particle_problem_add_body(self._p, self._coord_type(0,0,0), c_double(0), c_double(0))
def get_entity_name(edim, etag) :
for tag in gmsh.model.getPhysicalGroupsForEntity(edim, etag):
......@@ -697,8 +670,8 @@ class ParticleProblem :
def add_disk(t, stag) :
if t in addv : return
self.add_particle(x[t]+shift, 0, bndbody, material)
addv.add(t)
self.add_boundary_disk(x[t]+shift, 0., stag, material)
def add_segment(t0, t1, stag) :
key = (min(t0,t1),max(t0,t1))
......
......@@ -24,6 +24,7 @@ a = 1
g = np.array([0,-9.81])
r = 0.05
rho = 1000
bndbody = p.add_body((0,0), 0, 0)
p.add_boundary_segment([-a,-a],[-a,a],"bnd",material="Steel")
#p.add_boundary_segment([-a,a],[a,a],"bnd",material="Steel")
p.add_boundary_segment([a,a],[a,-a],"bnd",material="Steel")
......@@ -31,8 +32,8 @@ p.add_boundary_segment([a,-a],[-a,-a],"bnd",material="Steel")
#p.add_boundary_segment([a,-a],[0,-0.7],"bnd",material="Steel")
#p.add_boundary_segment([0,-0.7],[-a,-a],"bnd",material="Steel")
#p.add_boundary_disk([0,-0.7],0,"bnd","Steel")
p.add_boundary_disk([-a,a],0,"bnd","Steel")
p.add_boundary_disk([a,a],0,"bnd","Steel")
p.add_particle([-a,a],0,bndbody,"Steel")
p.add_particle([a,a],0,bndbody,"Steel")
tend=30
if friction:
p.set_friction_coefficient(.2,"Sand","Sand")
......
Supports Markdown
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