← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-pkg/yade/git-trunk] Rev 3716: adding examples with tetrehadron modeled by new Polyhedron class

 

------------------------------------------------------------
revno: 3716
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Tue 2013-10-15 19:11:37 +0200
message:
  adding examples with tetrehadron modeled by new Polyhedron class
added:
  examples/tetra/oneTetraPoly.py
  examples/tetra/twoTetrasPoly.py
modified:
  examples/tetra/oneTetra.py
  py/utils.py


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'examples/tetra/oneTetra.py'
--- examples/tetra/oneTetra.py	2013-10-04 12:03:13 +0000
+++ examples/tetra/oneTetra.py	2013-10-15 17:11:37 +0000
@@ -8,7 +8,7 @@
 v3 = (0,1,0)
 v4 = (0,0,1)
 
-t1 = tetra((v1,v2,v3,v4),color=(0,1,0))
+t1 = tetraOld((v1,v2,v3,v4),color=(0,1,0))
 O.bodies.append((t1))
 
 def changeVelocity():

=== added file 'examples/tetra/oneTetraPoly.py'
--- examples/tetra/oneTetraPoly.py	1970-01-01 00:00:00 +0000
+++ examples/tetra/oneTetraPoly.py	2013-10-15 17:11:37 +0000
@@ -0,0 +1,43 @@
+################################################################################
+# 
+# Script to test tetra gl functions and prescribed motion
+#
+################################################################################
+v1 = (0,0,0)
+v2 = (1,0,0)
+v3 = (0,1,0)
+v4 = (0,0,1)
+
+t1 = tetraPoly((v1,v2,v3,v4),color=(0,1,0))
+O.bodies.append((t1))
+
+def changeVelocity():
+	if O.iter == 50000:
+		t1.state.vel = (-1,0,0)
+	if O.iter == 100000:
+		t1.state.vel = (0,1,-1)
+	if O.iter == 150000:
+		t1.state.vel = (0,0,0)
+		t1.state.angMom = (0,0,10)
+	if O.iter == 200000:
+		O.pause()
+
+O.engines = [
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Polyhedra_Aabb()]),
+	InteractionLoop([],[],[]),
+	NewtonIntegrator(),
+	PyRunner(iterPeriod=1,command="changeVelocity()"),
+]
+O.step()
+
+
+try:
+	from yade import qt
+	qt.View()
+except:
+	pass
+
+O.dt = 1e-5
+t1.state.vel = (1,0,0)
+O.run()

=== added file 'examples/tetra/twoTetrasPoly.py'
--- examples/tetra/twoTetrasPoly.py	1970-01-01 00:00:00 +0000
+++ examples/tetra/twoTetrasPoly.py	2013-10-15 17:11:37 +0000
@@ -0,0 +1,116 @@
+################################################################################
+#
+# Python script to test tetra-tetra contact detection for different possible
+# contact modes. Some tests are run twice to test the symmetry of the law.
+#
+# It runs several tests, making pause before each one. If run with GUI, you can
+# adjust the viewer
+#
+# During the test, momentum, angular momentum and kinetic energy is tracked in
+# plot.data. You can use plot.plot() to see the results (in sime modes the
+# energy is not conserved for instace).
+#
+################################################################################
+from yade import qt,plot
+
+mat = PolyhedraMat()
+mat.density = 2600 #kg/m^3 
+mat.Ks = 2000000
+mat.Kn = 1E9 #Pa
+mat.frictionAngle = 0.5 #rad
+O.materials.append(mat)
+
+O.dt = 4e-5
+O.engines = [
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Polyhedra_Aabb()]),
+	InteractionLoop(
+		[Ig2_Polyhedra_Polyhedra_PolyhedraGeom()],
+		[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
+		[PolyhedraVolumetricLaw()]
+	),
+	NewtonIntegrator(damping=0),
+	PyRunner(iterPeriod=500,command="addPlotData()"),
+	PyRunner(iterPeriod=1,command="runExamples()"),
+]
+
+def addPlotData():
+	p = utils.momentum()
+	l = utils.angularMomentum()
+	plot.addData(
+		i = O.iter,
+		e = utils.kineticEnergy(),
+		px = p[0],
+		py = p[1],
+		pz = p[2],
+		lx = l[0],
+		ly = l[1],
+		lz = l[2],
+	)
+
+def prepareExample(vertices1,vertices2,vel1=(0,0,0),vel2=(0,0,0),amom1=(0,0,0),amom2=(0,0,0),label=''):
+	O.interactions.clear()
+	O.bodies.clear()
+	t1 = tetraPoly(vertices1,color=(0,1,0),wire=False)
+	t2 = tetraPoly(vertices2,color=(0,1,1),wire=False)
+	O.bodies.append((t1,t2))
+	t1.state.vel = vel1
+	t2.state.vel = vel2
+	t1.state.angMom = amom1
+	t2.state.angMom = amom2
+	if label: print label
+	O.pause()
+	O.step()
+
+def runExamples():
+	dt = 20000
+	if O.iter == 1:
+		vertices1 = ((0,0,0),(2,0,0),(0,2,0),(0,0,2))
+		vertices2 = ((1,1,1),(3,1,1),(1,3,1),(1,1,3))
+		prepareExample(vertices1,vertices2,vel2=(-1,-1,-1),label='\ntesting vertex-triangle contact...\n')
+	if O.iter == 1*dt:
+		plot.data = {}
+		vertices1 = ((1,1,1),(3,1,1),(1,3,1),(1,1,3))
+		vertices2 = ((0,0,0),(2,0,0),(0,2,0),(0,0,2))
+		prepareExample(vertices1,vertices2,vel1=(-1,-1,-1),label='\ntesting vertex-triangle contact 2...\n')
+	elif O.iter == 2*dt:
+		vertices1 = ((0,0,0),(0,0,1.1),(1,0,1),(0,1,.9))
+		vertices2 = ((0,.5,1.4),(-.5,.5,.6),(-1.6,0,1),(-1.6,1.1,1))
+		prepareExample(vertices1,vertices2,vel2=(1,0,0),amom2=(0,10,0),label='\ntesting edge-edge contact\n')
+	elif O.iter == 3*dt:
+		vertices1 = ((0,0,0),(0,0,1.1),(1,0,1),(0,1,.9))
+		vertices2 = ((-.5,.5,.6),(0,.5,1.4),(-1.6,1.1,1),(-1.6,0,1))
+		prepareExample(vertices1,vertices2,vel2=(1,0,0),amom2=(0,10,0),label='\ntesting edge-edge contact\n')
+	elif O.iter == 4*dt:
+		vertices1 = ((.1,-.4,-.3),(-.3,-.4,2),(3,-.2,2),(-.1,3,2))
+		vertices2 = ((.5,1.5,2.3),(1.5,.5,2.3),(2,2,3),(0,0,3))
+		prepareExample(vertices1,vertices2,vel2=(0,0,-1),amom2=(0,0,0),label='\ntesting edge-triangle contact\n')
+	elif O.iter == 5*dt:
+		vertices1 = ((.1,-.4,-.3),(-.3,-.4,2),(3,-.2,2),(-.1,3,2))
+		vertices2 = ((-.5,2.5,2.3),(2.5,-.5,2.3),(2,2,3),(0,0,3))
+		prepareExample(vertices1,vertices2,vel2=(0,0,-1),amom2=(0,0,0),label='\ntesting edge-triangle contact 2\n')
+	elif O.iter == 6*dt:
+		vertices1 = ((1,0,0),(0,1,0),(0,0,1),(1,1,1))
+		vertices2 = ((.5,.5,1.2),(0,.1,2),(2,0,2),(0,1,2))
+		prepareExample(vertices1,vertices2,vel2=(0,0,-1),label='\ntesting vertex-edge contact\n')
+	elif O.iter == 7*dt:
+		vertices1 = ((.5,.5,1.2),(0,.1,2),(2,0,2),(0,1,2))
+		vertices2 = ((1,0,0),(0,1,0),(0,0,1),(1,1,1))
+		prepareExample(vertices1,vertices2,vel1=(0,0,-1),label='\ntesting vertex-edge contact 2\n')
+	elif O.iter == 8*dt:
+		vertices1 = ((0,0,0),(1,0,0),(0,1,0),(0,0,1))
+		vertices2 = ((0,0,1.2),(.9,.8,2),(-.7,.61,2.3),(0,-.7,2.1))
+		prepareExample(vertices1,vertices2,vel2=(0,0,-1),label='\ntesting vertex-edge contact\n')
+	elif O.iter == 9*dt:
+		vertices1 = ((0,0,1.2),(.9,.8,2),(-.7,.61,2.3),(0,-.7,2.1))
+		vertices2 = ((0,0,0),(1,0,0),(0,1,0),(0,0,1))
+		prepareExample(vertices1,vertices2,vel1=(0,0,-1),label='\ntesting vertex-edge contact 2\n')
+	elif O.iter == 10*dt:
+		O.pause()
+
+viewer = qt.View()
+plot.plots = {'i':'e','i ':('px','py','pz'),'i  ':('lx','ly','lz')}
+
+O.step()
+O.step()
+O.step()

=== modified file 'py/utils.py'
--- py/utils.py	2013-10-07 09:10:44 +0000
+++ py/utils.py	2013-10-15 17:11:37 +0000
@@ -348,6 +348,24 @@
 	b.chain=chain
 	return b
 
+def tetraPoly(vertices,strictCheck=True,dynamic=True,fixed=False,wire=True,color=None,highlight=False,noBound=False,material=-1,mask=1,chain=-1):
+	"""Create tetrahedron (actually simple Polyhedra) with given parameters.
+
+	:param [Vector3,Vector3,Vector3,Vector3] vertices: coordinates of vertices in the global coordinate system.
+
+	See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""
+	b=Body()
+	b.shape = Polyhedra(v=vertices,color=color if color else randomColor(),wire=wire,highlight=highlight)
+	volume = b.shape.GetVolume()
+	inertia = b.shape.GetInertia()
+	center = b.shape.GetCentroid()
+	_commonBodySetup(b,volume,inertia,material,noBound=noBound,pos=center,fixed=fixed)
+	b.aspherical=False # mass and inertia are 0 anyway; fell free to change to ``True`` if needed
+	b.state.ori = b.shape.GetOri()
+	b.mask=mask
+	b.chain=chain
+	return b
+
 def tetra(vertices,strictCheck=True,dynamic=True,fixed=False,wire=True,color=None,highlight=False,noBound=False,material=-1,mask=1,chain=-1):
 	"""Create tetrahedron with given parameters.
 
@@ -379,6 +397,8 @@
 
 
 
+
+
 #def setNewVerticesOfFacet(b,vertices):
 #	center = inscribedCircleCenter(vertices[0],vertices[1],vertices[2])
 #	vertices = Vector3(vertices[0])-center,Vector3(vertices[1])-center,Vector3(vertices[2])-center