yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10139
[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