yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10079
[Branch ~yade-pkg/yade/git-trunk] Rev 3696: marginal modifications
------------------------------------------------------------
revno: 3696
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Fri 2013-10-04 14:03:13 +0200
message:
marginal modifications
removed:
examples/polyhedra/
examples/polyhedra/oneTetra.py
examples/polyhedra/twoTetras.py
added:
examples/tetra/
examples/tetra/oneTetra.py
examples/tetra/twoTetras.py
modified:
pkg/dem/ConcretePM.cpp
pkg/dem/Tetra.cpp
py/export.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
=== removed directory 'examples/polyhedra'
=== removed file 'examples/polyhedra/oneTetra.py'
--- examples/polyhedra/oneTetra.py 2013-09-11 21:43:09 +0000
+++ examples/polyhedra/oneTetra.py 1970-01-01 00:00:00 +0000
@@ -1,43 +0,0 @@
-################################################################################
-#
-# 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 = tetra((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_Tetra_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()
=== removed file 'examples/polyhedra/twoTetras.py'
--- examples/polyhedra/twoTetras.py 2013-09-08 11:12:42 +0000
+++ examples/polyhedra/twoTetras.py 1970-01-01 00:00:00 +0000
@@ -1,111 +0,0 @@
-################################################################################
-#
-# 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
-
-O.materials.append(ElastMat(young=1e10))
-O.dt = 4e-5
-O.engines = [
- ForceResetter(),
- InsertionSortCollider([Bo1_Tetra_Aabb()]),
- InteractionLoop(
- [Ig2_Tetra_Tetra_TTetraSimpleGeom()],
- [Ip2_ElastMat_ElastMat_NormPhys()],
- [Law2_TTetraSimpleGeom_NormPhys_Simple()]
- ),
- NewtonIntegrator(damping=0),
- PyRunner(iterPeriod=500,command="addPlotData()"),
- PyRunner(iterPeriod=1000,command="printt()"),
- 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 = tetra(vertices1,color=(0,1,0),wire=False)
- t2 = tetra(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()
=== added directory 'examples/tetra'
=== added file 'examples/tetra/oneTetra.py'
--- examples/tetra/oneTetra.py 1970-01-01 00:00:00 +0000
+++ examples/tetra/oneTetra.py 2013-10-04 12:03:13 +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 = tetra((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_Tetra_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/twoTetras.py'
--- examples/tetra/twoTetras.py 1970-01-01 00:00:00 +0000
+++ examples/tetra/twoTetras.py 2013-10-04 12:03:13 +0000
@@ -0,0 +1,110 @@
+################################################################################
+#
+# 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
+
+O.materials.append(ElastMat(young=1e10))
+O.dt = 4e-5
+O.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Tetra_Aabb()]),
+ InteractionLoop(
+ [Ig2_Tetra_Tetra_TTetraSimpleGeom()],
+ [Ip2_ElastMat_ElastMat_NormPhys()],
+ [Law2_TTetraSimpleGeom_NormPhys_Simple()]
+ ),
+ 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 = tetra(vertices1,color=(0,1,0),wire=False)
+ t2 = tetra(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 'pkg/dem/ConcretePM.cpp'
--- pkg/dem/ConcretePM.cpp 2013-08-26 21:42:20 +0000
+++ pkg/dem/ConcretePM.cpp 2013-10-04 12:03:13 +0000
@@ -345,9 +345,9 @@
Vector3r& Fs(phys->Fs); /* for python access */
const bool& isCohesive(phys->isCohesive);
-// Vector3r& epsTPl(phys->epsTPl);
#ifdef CPM_MATERIAL_MODEL
+ Vector3r& epsTPl(phys->epsTPl);
Real& epsNPl(phys->epsNPl);
const Real& dt = scene->dt;
const Real& dmgTau(phys->dmgTau);
@@ -394,9 +394,15 @@
#endif
sigmaN -= phys->isoPrestress;
-
- NNAN(kappaD); NNAN(epsFracture); NNAN(omega);
- NNAN(sigmaN); NNANV(sigmaT); NNAN(crossSection);
+
+ NNAN(sigmaN);
+ NNANV(sigmaT);
+ NNAN(crossSection);
+ if (!neverDamage) {
+ NNAN(kappaD);
+ NNAN(epsFracture);
+ NNAN(omega);
+ }
/* handle broken contacts */
if (epsN>0. && ((isCohesive && omega>omegaThreshold) || !isCohesive)) {
=== modified file 'pkg/dem/Tetra.cpp'
--- pkg/dem/Tetra.cpp 2013-09-11 21:43:09 +0000
+++ pkg/dem/Tetra.cpp 2013-10-04 12:03:13 +0000
@@ -280,7 +280,6 @@
Vector3r& contactPoint,
Real& penetrationVolume)
{
-
for (int i=0; i<6; i++) {
const Segment& sb = sB[i];
int ni = 0;
@@ -313,6 +312,7 @@
if ( !(b2a || b2b || b2c) ) { continue; }
int l = b1a? tsMap[j][0] : b1b? tsMap[j][1] : tsMap[j][2];
int m = b2a? tsMap[k][0] : b2b? tsMap[k][1] : tsMap[k][2];
+ if (sstMap[l][m] == -1) { continue; }
const Segment& sa1 = sA[l];
const Segment& sa2 = sA[m];
const Triangle& taN = tA[sstMap[l][m]];
@@ -549,6 +549,7 @@
geom->flag = flag; \
return true;
+
if (checkVertexToTriangleCase(t1,p2,s2,n,cp,V)) {
flag = 1;
SET_GEOM_AND_RETURN_TRUE
@@ -591,6 +592,7 @@
}
#undef SET_GEOM_AND_RETURN_TRUE
+
if (interaction->geom) {
TTetraSimpleGeom* geom = static_cast<TTetraSimpleGeom*>(interaction->geom.get());
=== modified file 'py/export.py'
--- py/export.py 2013-08-14 09:03:28 +0000
+++ py/export.py 2013-10-04 12:03:13 +0000
@@ -199,7 +199,7 @@
self.intrsSnapCount = startSnap
self.baseName = baseName
- def exportSpheres(self,ids='all',what=[],comment="comment",numLabel=None):
+ def exportSpheres(self,ids='all',what=[],comment="comment",numLabel=None,useRef=False):
"""exports spheres (positions and radius) and defined properties.
:param ids: if "all", then export all spheres, otherwise only spheres from integer list
@@ -208,6 +208,7 @@
:type what: [tuple(2)]
:param string comment: comment to add to vtk file
:param int numLabel: number of file (e.g. time step), if unspecified, the last used value + 1 will be used
+ :param bool useRef: if False (default), use current position of the spheres for export, use reference position otherwise
"""
allIds = False
if ids=='all':
@@ -226,7 +227,7 @@
outFile = open(fName, 'w')
outFile.write("# vtk DataFile Version 3.0.\n%s\nASCII\n\nDATASET POLYDATA\nPOINTS %d double\n"%(comment,n))
for b in bodies:
- pos = b.state.pos
+ pos = b.state.refPos if useRef else b.state.pos
outFile.write("%g %g %g\n"%(pos[0],pos[1],pos[2]))
outFile.write("\nPOINT_DATA %d\nSCALARS radius double 1\nLOOKUP_TABLE default\n"%(n))
for b in bodies:
=== modified file 'py/utils.py'
--- py/utils.py 2013-09-11 21:43:09 +0000
+++ py/utils.py 2013-10-04 12:03:13 +0000
@@ -956,19 +956,33 @@
self.relPositions = relPositions
class UnstructuredGrid:
- """EXPERIMENTAL. TODO docs
+ """EXPERIMENTAL.
+ Class representing triangulated FEM-like unstructured grid. It is used for transfereing data from ad to YADE and external FEM program. The main purpose of this class is to store information about individual grid vertices/nodes coords (since facets stores only coordinates of vertices in local coords) and to avaluate and/or apply nodal forces from contact forces (from actual contact force and contact point the force is distributed to nodes using linear approximation).
+
+ TODO rewrite to C++
+ TODO better docs
+
+ :param dict vertices: dict of {internal vertex label:vertex}, e.g. {5:(0,0,0),22:(0,1,0),23:(1,0,0)}
+ :param dict connectivityTable: dict of {internal element label:[indices of vertices]}, e.g. {88:[5,22,23]}
"""
def __init__(self,vertices={},connectivityTable={},**kw):
self.vertices = vertices
self.connectivityTable = connectivityTable
self.forces = dict( (i,Vector3(0,0,0)) for i in vertices)
self.build(**kw)
- def setup(self,vertices,connectivityTable,toSimulation=False,**kw):
+ def setup(self,vertices,connectivityTable,toSimulation=False,bodies=None,**kw):
+ """Sets new information to receiver
+
+ :param dict vertices: see constructor for explanation
+ :param dict connectivityTable: see constructor for explanation
+ :param bool toSimulation: if new information should be inserted to Yade simulation (create new bodies or not)
+ :param [[int]]|None bodies: list of list of bodies indices to be appended as clumps (thus no contact detection is done within one body)
+ """
self.vertices = dict( (i,v) for i,v in vertices.iteritems())
self.connectivityTable = dict( (i,e) for i,e in connectivityTable.iteritems())
self.build(**kw)
if toSimulation:
- self.toSimulation()
+ self.toSimulation(bodies)
def build(self,**kw):
self.elements = {}
for i,c in self.connectivityTable.iteritems():
@@ -977,12 +991,14 @@
elif len(c) == 4:
b = tetra([self.vertices[j] for j in c],**kw)
else:
- raise RuntimeError, "TODO"
+ raise RuntimeError, "Unsupported cell shape (should be triangle or tetrahedron)"
self.elements[i] = b
def resetForces(self):
for i in self.vertices:
self.forces[i] = Vector3(0,0,0)
def getForcesOfNodes(self):
+ """Computes forces for each vertex/node. The nodal force is computed from contact force and contact point using linear approximation
+ """
self.resetForces()
for i,e in self.elements.iteritems():
ie = self.connectivityTable[i]
@@ -1013,10 +1029,16 @@
self.forces[ie[3]] += f*w3/ww
return self.forces
def setPositionsOfNodes(self,newPoss):
+ """Sets new position of nodes and also updates all elements in the simulation
+
+ :param [Vector3] newPoss: list of new positions
+ """
for i in self.vertices:
self.vertices[i] = newPoss[i]
self.updateElements()
def updateElements(self):
+ """Updates positions of all elements in the simulation
+ """
for i,c in self.connectivityTable.iteritems():
e = self.elements[i]
e.state.pos = Vector3(0,0,0)
@@ -1025,5 +1047,12 @@
e.shape.vertices = [self.vertices[j] for j in c]
elif isinstance(e.shape,Tetra):
e.shape.v = [self.vertices[j] for j in c]
- def toSimulation(self):
- O.bodies.append(self.elements.values())
+ def toSimulation(self,bodies=None):
+ """Insert all elements to Yade simulation
+ """
+ if bodies:
+ e = self.elements.values()
+ for b in bodies:
+ O.bodies.appendClumped([e[i] for i in b])
+ else:
+ O.bodies.append(self.elements.values())
