yade-users team mailing list archive

Thread
Date
Re: [Question #680139]: rigid (gts)Pfacets?

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: Bettina Suhr <question680139@xxxxxxxxxxxxxxxxxxxxx>
Date: Wed, 10 Apr 2019 12:57:39 -0000
Reply-to: question680139@xxxxxxxxxxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx
Question #680139 on Yade changed:
https://answers.launchpad.net/yade/+question/680139

    Status: Answered => Open

Bettina Suhr is still having a problem:
Dear Bruno,

thank you very much for your advice. It’s a good idea to add
nodes+edges+pfacets to the clump.

In my script, I tried to do this but faced the following problem. Normally, the gtsPFacet() function calls a couple of subfunctions, which create the gridNodes, gridConnections and pfacets and directly append them to the scene with O.bodies.append(). Objects which are already registered to the scene cannot be part of a  clump (at least that is how I understand the error message of O.bodies.appendClumped(): IndexError: Body already has id 1538 set; appending such body (for the second time) is not allowed.)
 
Therefore,  I copied all functions called by  gtsPFacet() in my MWE and tried to adapt them with the aim: create all needed gridNodes, gridConnections and pfacets and at the very last line use  O.bodies.appendClumped() to add them all to the scene. 
The problem is that in my adapted function gridConnection an interaction should be created and parameters of i.phys and i.geom are set. My script crushes, because the createInteraction() function expects the the partners of the interaction to be already part of the scene (which is what I try to avoid). I don’t know how to proceed. At the end, the clump should be rigid (i.e. ignoring internal interactions). But this interaction is also needed for the contacts between the gridConnetion and other bodies(spheres or other clumped of gtsPfacets)?
 Is this the right way at all?  

Below is my script, any help would be appreciated.

Thanks and best regards,
Bettina


from yade import qt
from yade.gridpfacet import *
import gts, os.path, locale

locale.setlocale(locale.LC_ALL, 'en_US.UTF-8')   # Note: gts is locale-
dependent. If, for example, german locale is used, gts.read()-function
does not import floats normally

'''
if you get "Error: unsupported locale setting"
-> type as root: "dpkg-reconfigure locales"
-> choose "en_US.UTF-8" (press space to choose)
'''

#-------------------------------------------------------------------------------------
# adapted functions from trunk/py/gridpfacet.py to make clump of pfacets
#-------------------------------------------------------------------------------------

def ClumpedGtsPFacet(meshfile,shift=Vector3.Zero,scale=1.0,
							radius=1,wire=True,fixed=True,materialNodes=-1,material=-1,color=None):
	"""doc deleted
	"""
	import gts,yade.pack
	surf=gts.read(open(meshfile))
	surf.scale(scale,scale,scale)
	surf.translate(shift[0],shift[1],shift[2]) 
	nodesList=[]
	connList=[]
	pfList=[]
	
	for face in surf.faces():
		a=face.vertices()[0].coords()
		b=face.vertices()[1].coords()
		c=face.vertices()[2].coords()
		ClumpedPfacetCreator1([a,b,c],radius=radius,nodesList=nodesList,connList=connList,pfList=pfList,wire=wire,fixed=fixed,materialNodes=materialNodes,material=material,color=color)
		#print a,b,c
	pfIds=O.bodies.appendClumped(pfList+connList+nodesList)
	#return nodesIds,cylIds,pfIds


def ClumpedPfacetCreator1(vertices,radius,nodesList=[],connList=[],pfList=[],wire=False, fixed=True,materialNodes=-1,material=-1,color=None):
	"""doc deleted
	"""
	n=len(nodesList)
	k=[0,0,0]
	f=[0,0,0]
	u=0
	nod=0
	for i in vertices:
		u=0
		for (j,b) in enumerate(nodesList):
			if(i==nodesList[j].state.pos):
				f[nod]=j
				k[nod]=1
				u+=1
		nod+=1
		test=True
		#if(u==0):
		for (GN,b) in enumerate(nodesList):
			if(i==nodesList[GN].state.pos):
				  u=1
		if(u==0):
			#nodesIds.append( O.bodies.append(gridNode(i,radius,wire=wire,fixed=fixed,material=materialNodes,color=color)) )
			nodesList.append( gridNode(i,radius,wire=wire,fixed=fixed,material=materialNodes,color=color)) 

	if(k==[0,0,0]):
		pfObj=ClumpedPfacetCreator3(n,n+1,n+2,nodesList, connList=connList, pfList=pfList,wire=wire,material=material,color=color,fixed=fixed )
	if(k==[1,0,0]):
		pfObj=ClumpedPfacetCreator3(f[0],n,n+1,nodesList, connList=connList, pfList=pfList,wire=wire,material=material,color=color,fixed=fixed )
	if(k==[0,1,0]):
		pfObj=ClumpedPfacetCreator3(n,f[1],n+1, nodesList, connList=connList, pfList=pfList,wire=wire,material=material,color=color,fixed=fixed )
	if(k==[0,0,1]):
		pfObj=ClumpedPfacetCreator3(n, n+1,f[2],nodesList, connList=connList, pfList=pfList,wire=wire,material=material,color=color,fixed=fixed )
	if(k==[1,1,0]):
		pfObj=ClumpedPfacetCreator3(f[0],f[1],n,nodesList, connList=connList, pfList=pfList,wire=wire,material=material,color=color,fixed=fixed )
	if(k==[0,1,1]):
		pfObj=ClumpedPfacetCreator3(n,f[1],f[2],nodesList,connList=connList, pfList=pfList,wire=wire,material=material,color=color,fixed=fixed )
	if(k==[1,0,1]):
		pfObj=ClumpedPfacetCreator3(f[0],n,f[2],nodesList,connList=connList, pfList=pfList,wire=wire,material=material,color=color,fixed=fixed )
	if(k==[1,1,1]):
		pfObj=ClumpedPfacetCreator3(f[0],f[1],f[2],nodesList, connList=connList, pfList=pfList,wire=wire,material=material,color=color,fixed=fixed )
	return pfObj


def ClumpedPfacetCreator3(id1,id2,id3,nodesList=[], connList=[],pfList=[], wire=True, material=-1, color=None, fixed=True,mask=-1):
	"""doc deleted
	"""
	IdList=[]
	radius=nodesList[id1].shape.radius
	try:
        #cylIds.append(O.bodies.append(gridConnection(id1,id2,radius=radius,material=material,color=color,wire=wire,mask=mask)))
		connList.append(ClumpedGridConnection(id1,id2,nodesList=nodesList, radius=radius, material=material, color=color,wire=wire,mask=mask))
		IdList.append([id1,id2])
	except:
		print "skipped a ", id1, id2
		pass
	try:
        #cylIds.append(O.bodies.append(gridConnection(id2,id3,radius=radius,material=material,color=color,wire=wire,mask=mask)))
		connList.append(ClumpedGridConnection(id2,id3,nodesList,radius=radius,material=material, color=color,wire=wire,mask=mask))
		IdList.append([id2,id3])
	except:
		print "skipped b ", id2, id3
		pass	      
	try:
        #cylIds.append(O.bodies.append(gridConnection(id3,id1,radius=radius,material=material,color=color,wire=wire,mask=mask)))
		connList.append(ClumpedGridConnection(id3,id1,nodesList,radius=radius,material=material, color=color,wire=wire,mask=mask))
		IdList.append([id3,id1])
	except:
		print "skipped c ", id3, id1
		pass	      
	#pfIds.append(O.bodies.appendClumped())
	pfList.append(ClumpedPfacet(id1,id2,id3,nodesList,connList,IdList,wire=wire,material=material,color=color,mask=mask))



def ClumpedGridConnection(id1,id2,nodesList,radius,wire=False,color=None,highlight=False, material=-1,mask=1,cellDist=None): 
	"""doc deleted
	"""
	b=Body()
	b.shape=GridConnection(radius=radius,color=color if color else utils.randomColor(),wire=wire,highlight=highlight)
	#sph1=O.bodies[id1] ; sph2=O.bodies[id2]
	sph1=nodesList[id1] ; sph2=nodesList[id2]
	print "Big problem here: createInteraction"
	i=createInteraction(id1,id2)#(sph1, sph2)#
	nodeMat=sph1.material
	b.shape.node1=sph1 ; b.shape.node2=sph2
	sph1.shape.addConnection(b) ; sph2.shape.addConnection(b)
	#if(O.periodic):##UNTESTED!!
		#if(cellDist!=None):
			#i.cellDist=cellDist
		#segt=sph2.state.pos + O.cell.hSize*i.cellDist - sph1.state.pos
	#else: 
	if 1:
		segt=sph2.state.pos - sph1.state.pos
	L=segt.norm()
	V=0.5*L*math.pi*radius**2
	geomInert=(2./5.)*V*radius**2
	utils._commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert),material,pos=sph1.state.pos,dynamic=False,fixed=True)
	sph1.state.mass = sph1.state.mass + V*nodeMat.density
	sph2.state.mass = sph2.state.mass + V*nodeMat.density
	for k in [0,1,2]:
		sph1.state.inertia[k] = sph1.state.inertia[k] + geomInert*nodeMat.density
		sph2.state.inertia[k] = sph2.state.inertia[k] + geomInert*nodeMat.density
	b.aspherical=False
	#if O.periodic:##UNTESTED!!
		#i.phys.unp= -(sph2.state.pos + O.cell.hSize*i.cellDist - sph1.state.pos).norm() + sph1.shape.radius + sph2.shape.radius
		#b.shape.periodic=True
		#b.shape.cellDist=i.cellDist
	#else:
	if 1:
		i.phys.unp= -(sph2.state.pos - sph1.state.pos).norm() + sph1.shape.radius + sph2.shape.radius	
	i.geom.connectionBody=b
	I=math.pi*(2.*radius)**4/64.
	E=nodeMat.young
	i.phys.kn=E*math.pi*(radius**2)/L
	i.phys.kr=E*I/L
	i.phys.ks=12.*E*I/(L**3)
	G=E/(2.*(1+nodeMat.poisson))
	i.phys.ktw=2.*I*G/L
	b.mask=mask
	return b


def ClumpedPfacet(id1,id2,id3,nodesList,connList, IdList,wire=True,color=None,highlight=False, material=-1,mask=1, cellDist=None):
	"""doc deleted
	"""
	b=Body()
	#GridN1=O.bodies[id1]; GridN2=O.bodies[id2]; GridN3=O.bodies[id3] 
	GridN1=nodesList[id1]; GridN2=nodesList[id2]; GridN3=nodesList[id3] 
	b.shape=PFacet(color=color if color else randomColor(),wire=wire,highlight=highlight,node1=GridN1,node2=GridN2,node3=GridN3)
	GridN1.bounded=False; GridN2.bounded=False; GridN3.bounded=False
	#GridC1=O.bodies[O.interactions[id1,id2].geom.connectionBody.id]
	#GridC2=O.bodies[O.interactions[id2,id3].geom.connectionBody.id]
	#GridC3=O.bodies[O.interactions[id1,id3].geom.connectionBody.id]
	#print connList, id1, id2, id3
	GridC1=connList[-3]
	GridC2=connList[-2]
	GridC3=connList[-1]
	GridC1.bounded=False
	GridC2.bounded=False
	GridC3.bounded=False
	
	b.shape.conn1=GridC1
	b.shape.conn2=GridC2 
	b.shape.conn3=GridC3
	
	b.shape.radius=GridN1.shape.radius
	GridC1.shape.addPFacet(b) 
	GridC2.shape.addPFacet(b)
	GridC3.shape.addPFacet(b)
	GridN1.shape.addPFacet(b); GridN2.shape.addPFacet(b); GridN3.shape.addPFacet(b)
	
	V=0
	
	utils._commonBodySetup(b,V,Vector3(0,0,0),material,pos=GridN1.state.pos,dynamic=False,fixed=True)
	b.aspherical=False # mass and inertia are lumped into the GridNodes
	b.mask=mask
	return b


#-------------------------------------------------------------------------------------

################
### ENGINES  ###
################

O.engines=[
	ForceResetter(),
	InsertionSortCollider([
		Bo1_Wall_Aabb(),
		Bo1_PFacet_Aabb(),
	],sortThenCollide=True),
	InteractionLoop(
	[
        Ig2_GridNode_GridNode_GridNodeGeom6D(),
		Ig2_Wall_PFacet_ScGeom(),Ig2_Wall_Sphere_ScGeom()
	],
	[
        Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True),
		Ip2_FrictMat_FrictMat_FrictPhys()],
	[
        Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
		Law2_ScGeom_FrictPhys_CundallStrack(),
		Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
		Law2_GridCoGridCoGeom_FrictPhys_CundallStrack()
	]),
    #GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.5,label='ts'), 
	NewtonIntegrator(gravity=(0,-9.81,0),damping=0.1,label='newton')
]

#disabled GlobalStiffnessTimeStepper due to problems, adapt dt if necessary!
O.dt=1e-6

#################
### MATERIAL  ###
#################

O.materials.append(CohFrictMat(young=1e8,poisson=0.3,density=2650,frictionAngle=radians(20),normalCohesion=1e100,shearCohesion=1e100,momentRotationLaw=True,label='gridNodeMat'))
O.materials.append(FrictMat(young=1e8,poisson=0.3,density=2650,frictionAngle=radians(20),label='pFacetMat'))

###################
### IMPORT MESH ###
###################
radius=0.5
wire=False
fixed=False

z=-1.6
color=[0,0,1]

#--- original gtsPFacet construction --- WORKING
#nodesIds0,cylIds0,pfIds0 =  gtsPFacet('box.gts',shift=(0.,0.,0.),scale=2.,radius=radius,wire=wire,fixed=fixed,materialNodes='gridNodeMat',material='pFacetMat',color=color)

#--- try to make a clump out of it --- NOT WORKING
ClumpedGtsPFacet('box.gts',shift=(0,0,0),scale=2.,radius=radius,wire=wire,fixed=fixed,materialNodes='gridNodeMat',material='pFacetMat',color=color)

#--- mass and inertia should be updated?
O.bodies.updateClumpProperties(discretization=10)

#####################
#####   Wall      ###
#####################

O.bodies.append(utils.wall(position=z,sense=0,
axis=1,color=Vector3(1,0,0),material='pFacetMat'))

##########
## VIEW ##
##########

qt.Controller()
qtv = qt.View()
qtr = qt.Renderer()
qtr.light2=True
qtr.lightPos=Vector3(1200,1500,500)
qtr.bgColor=[1,1,1]
qtv.ortho=True


O.saveTmp()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.