Commit eec1c999 authored by Jan Stransky's avatar Jan Stransky

adding examples with tetrehadron modeled by new Polyhedron class

parent edb5e0d1
......@@ -8,7 +8,7 @@ v2 = (1,0,0)
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():
......
################################################################################
#
# 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()
################################################################################
#
# 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()
......@@ -348,6 +348,24 @@ def facet(vertices,dynamic=None,fixed=True,wire=True,color=None,highlight=False,
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 tetra(vertices,strictCheck=True,dynamic=True,fixed=False,wire=True,color=Non
#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
......
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