yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02178
[Branch ~yade-dev/yade/trunk] Rev 1775: Yade-multi, VTKRecorder enhancements, &c
Merge authors:
Yade <yade@pc7317>
------------------------------------------------------------
revno: 1775 [merge]
committer: Václav Šmilauer <vaclav@falx>
branch nick: trunk
timestamp: Thu 2009-10-22 09:05:16 +0200
message:
Yade-multi, VTKRecorder enhancements, &c
modified:
gui/py/PythonTCPServer.py
gui/py/PythonUI_rc.py
gui/py/yade-multi
pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp
pkg/dem/meta/ConcretePM.cpp
pkg/dem/meta/ConcretePM.hpp
py/_eudoxos.cpp
py/_utils.cpp
py/plot.py
py/utils.py
py/yadeWrapper/yadeWrapper.cpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== modified file 'gui/py/PythonTCPServer.py'
--- gui/py/PythonTCPServer.py 2008-10-14 08:48:01 +0000
+++ gui/py/PythonTCPServer.py 2009-10-21 15:22:14 +0000
@@ -3,6 +3,17 @@
from yade.wrapper import *
+class InfoSocketProvider(SocketServer.BaseRequestHandler):
+ """Class providing dictionary of important simulation information,
+ without authenticating the access"""
+ def handle(self):
+ import pickle, os
+ O=Omega()
+ ret=dict(iter=O.iter,dt=O.dt,stopAtIter=O.stopAtIter,realtime=O.realtime,time=O.time,id=O.tags['id'] if O.tags.has_key('id') else None,threads=os.environ['OMP_NUM_THREADS'] if os.environ.has_key('OMP_NUM_THREADS') else '0',numBodies=len(O.bodies),numIntrs=len(O.interactions))
+ print 'returning', ret
+ self.request.send(pickle.dumps(ret))
+
+
class PythonConsoleSocketEmulator(SocketServer.BaseRequestHandler):
def setup(self):
if not self.client_address[0].startswith('127.0.0'):
@@ -61,8 +72,8 @@
print self.client_address, 'disconnected!'
self.request.send('\nBye ' + str(self.client_address) + '\n')
-class PythonTCPServer:
- def __init__(self,minPort=9000,host='',maxPort=65536,background=True):
+class GenericTCPServer:
+ def __init__(self,handler,title,cookie=True,minPort=9000,host='',maxPort=65536,background=True):
import socket, random
self.port=-1
self.host=host
@@ -70,11 +81,14 @@
if maxPort==None: maxPort=minPort
while self.port==-1 and tryPort<=maxPort:
try:
- self.server=SocketServer.ThreadingTCPServer((host,tryPort),PythonConsoleSocketEmulator)
+ self.server=SocketServer.ThreadingTCPServer((host,tryPort),handler)
self.port=tryPort
- self.server.cookie=''.join([i for i in random.sample('yadesucks',6)])
- print "Python TCP server listening on %s:%d, auth cookie is `%s'"%(host if host else 'localhost',self.port,self.server.cookie)
- self.server.authenticated=[]
+ if cookie:
+ self.server.cookie=''.join([i for i in random.sample('yadesucks',6)])
+ self.server.authenticated=[]
+ print title,"on %s:%d, auth cookie `%s'"%(host if host else 'localhost',self.port,self.server.cookie)
+ else:
+ print title,"on %s:%d"%(host if host else 'localhost',self.port)
if background:
import thread; thread.start_new_thread(self.server.serve_forever,())
else: self.server.serve_forever()
@@ -82,7 +96,8 @@
tryPort+=1
if self.port==-1: raise RuntimeError("No free port to listen on in range %d-%d"%(minPort,maxPort))
+
if __name__=='__main__':
- p=PythonTCPServer(background=False)
+ p=GenericTCPServer(PythonConsoleSocketEmulator,'Python TCP server',background=False)
#while True: time.sleep(2)
=== modified file 'gui/py/PythonUI_rc.py'
--- gui/py/PythonUI_rc.py 2009-10-11 11:45:19 +0000
+++ gui/py/PythonUI_rc.py 2009-10-21 15:22:14 +0000
@@ -110,8 +110,9 @@
## run the TCP server
import yade.PythonTCPServer
-srv=yade.PythonTCPServer.PythonTCPServer(minPort=9000)
+srv=yade.PythonTCPServer.GenericTCPServer(handler=yade.PythonTCPServer.PythonConsoleSocketEmulator,title='TCP python prompt',cookie=True,minPort=9000)
yade.runtime.cookie=srv.server.cookie
+info=yade.PythonTCPServer.GenericTCPServer(handler=yade.PythonTCPServer.InfoSocketProvider,title='TCP info provider',cookie=False,minPort=21000)
sys.stdout.flush()
runtime.argv=[runtime.script]+runtime.argv
sys.argv=runtime.argv # could be [] as well
=== modified file 'gui/py/yade-multi'
--- gui/py/yade-multi 2009-10-17 06:36:21 +0000
+++ gui/py/yade-multi 2009-10-21 15:22:14 +0000
@@ -8,7 +8,7 @@
class JobInfo():
def __init__(self,num,id,command,log,nSlots):
self.started,self.finished,self.duration,self.exitStatus=None,None,None,None
- self.command=command; self.num=num; self.log=log; self.id=id; self.nSlots=nSlots
+ self.command=command; self.num=num; self.log=log; self.id=id; self.nSlots=nSlots; self.infoSocket=None
self.status='PENDING'
self.threadNum=None
def saveInfo(self):
@@ -23,14 +23,41 @@
finished: %s
"""%(self.id,self.exitStatus,'OK' if self.exitStatus==0 else 'FAILED',self.duration,self.command,time.asctime(time.localtime(self.started)),time.asctime(time.localtime(self.finished))));
log.close()
+ def getInfoDict(self):
+ if self.status!='RUNNING': return None
+ if not self.infoSocket:
+ for l in open(self.log,'r'):
+ if not l.startswith('TCP info provider on'): continue
+ hostPort=l[:-1].split()[4].split(':')
+ self.infoSocket=hostPort[0],int(hostPort[1])
+ if not self.infoSocket:
+ return None
+ import socket,pickle
+ s=socket.socket(socket.AF_INET,socket.SOCK_STREAM)
+ s.connect(self.infoSocket)
+ data=s.recv(2048)
+ s.close()
+ return pickle.loads(data)
+
def htmlStats(self):
ret='<tr>'
ret+='<td>%s</td>'%self.id
if self.status=='PENDING': ret+='<td bgcolor="grey">(pending)</td>'
elif self.status=='RUNNING': ret+='<td bgcolor="yellow">%s</td>'%t2hhmmss(time.time()-self.started)
- elif self.status=='DONE': ret+='<td bgcolor="%s">%s</td>'%('green' if self.exitStatus==0 else 'red',self.duration)
+ elif self.status=='DONE': ret+='<td bgcolor="%s">%s</td>'%('lime' if self.exitStatus==0 else 'red',self.duration)
+ info=self.getInfoDict()
+ if info:
+ ret+='<td>'
+ if info['stopAtIter']>0:
+ ret+='<nobr>%2.2f%% done</nobr><br/><nobr>step %d/%d</nobr>'%(info['iter']*100./info['stopAtIter'],info['iter'],info['stopAtIter'])
+ else: ret+='<nobr>step %d</nobr>'%(info['iter'])
+ if info['realtime']!=0: ret+='<br/><nobr>avg %g/sec</nobr>'%(info['iter']/info['realtime'])
+ ret+='<br/><nobr>%d bodies</nobr><br/><nobr>%d intrs</nobr>'%(info['numBodies'],info['numIntrs'])
+ ret+='</td>'
+ else:
+ ret+='<td> </td>'
+ ret+='<td>%d</td>'%self.nSlots
ret+='<td>%s</td>'%self.command
- ret+='<td>%d</td>'%self.nSlots
ret+='</tr>'
return ret
def t2hhmmss(dt): return '%02d:%02d:%02d'%(dt//3600,(dt%3600)//60,(dt%60))
@@ -38,15 +65,19 @@
def globalHtmlStats():
t0=min([j.started for j in jobs if j.started!=None])
unfinished=len([j for j in jobs if j.status!='DONE'])
+ usedSlots=sum([j.nSlots for j in jobs if j.status=='RUNNING'])
+ global maxJobs
if unfinished:
ret='<p>Running for %s, since %s.</p>'%(t2hhmmss(time.time()-t0),time.ctime(t0))
else:
failed=len([j for j in jobs if j.exitStatus!=0])
lastFinished=max([j.finished for j in jobs])
- ret='<p><span style="background-color:%s">Finished</span>, idle for %s, running time %s since %s'%('red' if failed else 'green',t2hhmmss(time.time()-lastFinished),t2hhmmss(sum([j.finished-j.started for j in jobs if j.started is not None])),time.ctime(t0))
+ ret='<p><span style="background-color:%s">Finished</span>, idle for %s, running time %s since %s.</p>'%('red' if failed else 'lime',t2hhmmss(time.time()-lastFinished),t2hhmmss(sum([j.finished-j.started for j in jobs if j.started is not None])),time.ctime(t0))
+ ret+='<p>Pid %d</p>'%(os.getpid())
+ ret+='<p>%d slots available, %d used, %d free.</p>'%(maxJobs,usedSlots,maxJobs-usedSlots)
ret+='<h3>Jobs</h3>'
nFailed=len([j for j in jobs if j.status=='DONE' and j.exitStatus!=0])
- ret+='<p><b>%d</b> total, <b>%d</b> <span style="background-color:yellow">running</span>, <b>%d</b> <span style="background-color:green">done</span>%s</p>'%(len(jobs),len([j for j in jobs if j.status=='RUNNING']), len([j for j in jobs if j.status=='DONE']),' (<b>%d <span style="background-color:red"><b>failed</b></span>)'%nFailed if nFailed>0 else '')
+ ret+='<p><b>%d</b> total, <b>%d</b> <span style="background-color:yellow">running</span>, <b>%d</b> <span style="background-color:lime">done</span>%s</p>'%(len(jobs),len([j for j in jobs if j.status=='RUNNING']), len([j for j in jobs if j.status=='DONE']),' (<b>%d <span style="background-color:red"><b>failed</b></span>)'%nFailed if nFailed>0 else '')
return ret
from BaseHTTPServer import BaseHTTPRequestHandler,HTTPServer
@@ -57,9 +88,9 @@
self.send_header('Content-type','text/html')
self.send_header('refresh','5')
self.end_headers()
- self.wfile.write('<HTML><TITLE>Yade-multi overview</TITLE><BODY>')
+ self.wfile.write('<HTML><TITLE>Yade-multi at %s overview</TITLE><BODY>'%(socket.gethostname()))
self.wfile.write(globalHtmlStats())
- self.wfile.write('<TABLE border=1><tr><th>id</th><th>status</th><th>command</th><th>slots</th></tr>')
+ self.wfile.write('<TABLE border=1><tr><th>id</th><th>status</th><th>info</th><th>slots</th><th>command</th></tr>')
for j in jobs:
self.wfile.write(j.htmlStats())
self.wfile.write('</TABLE></BODY></HTML>')
@@ -105,14 +136,22 @@
thread.start_new_thread(runJob,(j,))
break
time.sleep(.5)
-
-
-import sys,re,optparse
-def getNumCores(): return len([l for l in open('/proc/cpuinfo','r') if l.find('processor')==0])
+ sys.stdout.flush()
+
+
+import sys,re,optparse,os
+def getNumCores():
+ nCpu=len([l for l in open('/proc/cpuinfo','r') if l.find('processor')==0])
+ if os.environ.has_key("OMP_NUM_THREADS"): return min(int(os.environ['OMP_NUM_THREADS']),nCpu)
+ return nCpu
+numCores=getNumCores()
parser=optparse.OptionParser(usage='%prog [options] TABLE SIMULATION.py\n\n %prog runs yade simulation multiple times with different parameters.\n See http://yade.wikia.com/wiki/ScriptParametricStudy for details.')
-parser.add_option('-j','--jobs',dest='maxJobs',type='int',help="Maximum number of simultaneous jobs to run (default: number of cores, i.e. %d)"%getNumCores(),metavar='NUM',default=getNumCores())
-parser.add_option('--log',dest='logFormat',help='Format of log files -- must contain a % or @, which will be replaced by line number or by description column respectively (default: SIMULATION.@.log)',metavar='FORMAT')
+parser.add_option('-j','--jobs',dest='maxJobs',type='int',help="Maximum number of simultaneous threads to run (default: number of cores, further limited by OMP_NUM_THREADS if set by the environment: %d)"%numCores,metavar='NUM',default=numCores)
+parser.add_option('--job-threads',dest='defaultThreads',type='int',help="Default number of threads for one job; can be overridden by per-job OMP_NUM_THREADS. Defaults to allocate all available cores (%d) for each job."%numCores,metavar='NUM',default=numCores)
+parser.add_option('--force-threads',action='store_true',dest='forceThreads')
+parser.add_option('--log',dest='logFormat',help='Format of job log files -- must contain a % or @, which will be replaced by line number or by description column respectively (default: SIMULATION.@.log)',metavar='FORMAT')
+parser.add_option('--global-log',dest='globalLog',help='Filename where to redirect output of yade-multi itself (as opposed to --log); if not specified (default), stdout/stderr are used',metavar='FILE')
parser.add_option('-l','--lines',dest='lineList',help='Lines of TABLE to use, in the format 2,3-5,8,11-13 (default: all available lines in TABLE)',metavar='LIST')
parser.add_option('--nice',dest='nice',type='int',help='Nice value of spawned jobs (default: 10)',default=10)
parser.add_option('--executable',dest='executable',help='Name of the program to run (default: %s)'%sys.argv[0][:-6],default=sys.argv[0][:-6],metavar='FILE') ## strip the '-multi' extension
@@ -120,7 +159,11 @@
parser.add_option('--dry-run',action='store_true',dest='dryRun',help='Do not actually run (useful for getting gnuplot only, for instance)',default=False)
parser.add_option('--http-wait',action='store_true',dest='httpWait',help='Do not quit if still serving overview over http repeatedly',default=False)
opts,args=parser.parse_args()
-logFormat,lineList,maxJobs,nice,executable,gnuplotOut,dryRun,httpWait=opts.logFormat,opts.lineList,opts.maxJobs,opts.nice,opts.executable,opts.gnuplotOut,opts.dryRun,opts.httpWait
+logFormat,lineList,maxJobs,nice,executable,gnuplotOut,dryRun,httpWait,globalLog=opts.logFormat,opts.lineList,opts.maxJobs,opts.nice,opts.executable,opts.gnuplotOut,opts.dryRun,opts.httpWait,opts.globalLog
+
+if globalLog:
+ sys.stderr=open(globalLog,"w")
+ sys.stdout=sys.stderr
if len(args)!=2:
#print "Exactly two non-option arguments must be specified -- parameter table and script to be run.\n"
@@ -192,18 +235,26 @@
else: logFile=logFile.replace('@',str(l))
logFile=logFile.replace('!','')
envVars=[]
- nSlots=1
+ nSlots=opts.defaultThreads
for col,head in enumerate(headings):
if head=='!EXEC': executable=values[l][col]
if head=='!OMP_NUM_THREADS':
nSlots=int(values[l][col]); maxCpu=getNumCores()
+ continue
if head[0]=='!': envVars+=['%s=%s'%(head[1:],values[l][col])]
- if nSlots>maxJobs: logging.warning('WARNING: job #%d wants %d slots but only %d are available'%(i,nSlots,maxJobs))
- jobs.append(JobInfo(i,idStrings[l] if idStrings else '#'+str(i),'PARAM_TABLE=%s:%d %s nice -n %d %s -N PythonUI -- -n -x %s > %s 2>&1'%(table,l,' '.join(envVars),nice,executable,simul,pipes.quote(logFile)),logFile,nSlots))
+ if nSlots>maxJobs:
+ if opts.forceThreads:
+ logging.info('Forcing job #%d to use only %d slots (max available) instead of %d requested'%(i,maxJobs,nSlots))
+ nSlots=maxJobs
+ else:
+ logging.warning('WARNING: job #%d will use %d slots but only %d are available'%(i,nSlots,maxJobs))
+ jobs.append(JobInfo(i,idStrings[l] if idStrings else '#'+str(i),'PARAM_TABLE=%s:%d OMP_NUM_THREADS=%d %s nice -n %d %s -N PythonUI -- -n -x %s > %s 2>&1'%(table,l,nSlots,' '.join(envVars),nice,executable,simul,pipes.quote(logFile)),logFile,nSlots))
print "Job summary:"
for job in jobs:
print ' #%d (%s%s):'%(job.num,job.id,'' if job.nSlots==1 else '/%d'%job.nSlots),job.command
+print "Master process pid",os.getpid()
+sys.stdout.flush()
httpLastServe=0
runHttpStatsServer()
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-10-17 06:36:21 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-10-21 15:22:14 +0000
@@ -81,15 +81,24 @@
vtkSmartPointer<vtkFloatArray> intrForceN = vtkSmartPointer<vtkFloatArray>::New();
intrForceN->SetNumberOfComponents(1);
intrForceN->SetName("forceN");
+ vtkSmartPointer<vtkFloatArray> intrAbsForceT = vtkSmartPointer<vtkFloatArray>::New();
+ intrAbsForceT->SetNumberOfComponents(3);
+ intrAbsForceT->SetName("absForceT");
// extras for CPM
if(recActive[REC_CPM]) CpmStateUpdater::update(rootBody);
vtkSmartPointer<vtkFloatArray> cpmDamage = vtkSmartPointer<vtkFloatArray>::New();
cpmDamage->SetNumberOfComponents(1);
cpmDamage->SetName("cpmDamage");
- vtkSmartPointer<vtkFloatArray> cpmStress = vtkSmartPointer<vtkFloatArray>::New();
- cpmStress->SetNumberOfComponents(3);
- cpmStress->SetName("cpmStress");
+ vtkSmartPointer<vtkFloatArray> cpmSigma = vtkSmartPointer<vtkFloatArray>::New();
+ cpmSigma->SetNumberOfComponents(3);
+ cpmSigma->SetName("cpmSigma");
+ vtkSmartPointer<vtkFloatArray> cpmSigmaM = vtkSmartPointer<vtkFloatArray>::New();
+ cpmSigmaM->SetNumberOfComponents(1);
+ cpmSigmaM->SetName("cpmSigmaM");
+ vtkSmartPointer<vtkFloatArray> cpmTau = vtkSmartPointer<vtkFloatArray>::New();
+ cpmTau->SetNumberOfComponents(3);
+ cpmTau->SetName("cpmTau");
if(recActive[REC_INTR]){
// save body positions, referenced by ids by vtkLine
@@ -111,6 +120,8 @@
if(recActive[REC_CPM]){
const CpmPhys* phys = YADE_CAST<CpmPhys*>(I->interactionPhysics.get());
intrForceN->InsertNextValue(phys->Fn);
+ float fs[3]={abs(phys->shearForce[0]),abs(phys->shearForce[1]),abs(phys->shearForce[2])};
+ intrAbsForceT->InsertNextTupleValue(fs);
}
}
}
@@ -138,9 +149,13 @@
}
if (recActive[REC_CPM]) {
cpmDamage->InsertNextValue(YADE_PTR_CAST<CpmMat>(b->physicalParameters)->normDmg);
- const Vector3r& ss=YADE_PTR_CAST<CpmMat>(b->physicalParameters)->avgStress;
+ const Vector3r& ss=YADE_PTR_CAST<CpmMat>(b->physicalParameters)->sigma;
+ const Vector3r& tt=YADE_PTR_CAST<CpmMat>(b->physicalParameters)->tau;
float s[3]={ss[0],ss[1],ss[2]};
- cpmStress->InsertNextTupleValue(s);
+ float t[3]={tt[0],tt[1],tt[2]};
+ cpmSigma->InsertNextTupleValue(s);
+ cpmSigmaM->InsertNextValue((ss[0]+ss[1]+ss[2])/3.);
+ cpmTau->InsertNextTupleValue(t);
}
continue;
}
@@ -185,7 +200,9 @@
if (recActive[REC_VELOCITY]) spheresUg->GetPointData()->AddArray(spheresVelocity);
if (recActive[REC_CPM]) {
spheresUg->GetPointData()->AddArray(cpmDamage);
- spheresUg->GetPointData()->AddArray(cpmStress);
+ spheresUg->GetPointData()->AddArray(cpmSigma);
+ spheresUg->GetPointData()->AddArray(cpmSigmaM);
+ spheresUg->GetPointData()->AddArray(cpmTau);
}
vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
if(compress) writer->SetCompressor(compressor);
@@ -212,7 +229,10 @@
vtkSmartPointer<vtkUnstructuredGrid> intrUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
intrUg->SetPoints(intrBodyPos);
intrUg->SetCells(VTK_LINE, intrCells);
- if (recActive[REC_CPM]) intrUg->GetCellData()->AddArray(intrForceN);
+ if (recActive[REC_CPM]){
+ intrUg->GetCellData()->AddArray(intrForceN);
+ intrUg->GetCellData()->AddArray(intrAbsForceT);
+ }
vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
if(compress) writer->SetCompressor(compressor);
string fn=fileName+"intrs."+lexical_cast<string>(rootBody->currentIteration)+".vtu";
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-10-17 06:36:21 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-10-21 15:22:14 +0000
@@ -20,7 +20,7 @@
virtual void action(MetaBody*);
private:
- REGISTER_ATTRIBUTES(PeriodicEngine,(recorders)(fileName)(compress));
+ REGISTER_ATTRIBUTES(PeriodicEngine,(recorders)(fileName)(compress)(skipNondynamic));
REGISTER_CLASS_AND_BASE(VTKRecorder,PeriodicEngine);
DECLARE_LOGGER;
};
=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp 2009-10-17 06:36:21 +0000
+++ pkg/dem/meta/ConcretePM.cpp 2009-10-21 15:22:14 +0000
@@ -357,11 +357,13 @@
Vector3r stress=(1./phys->crossSection)*(phys->normalForce+phys->shearForce);
const Vector3r& p1(geom->se31.position); const Vector3r& cp(geom->contactPoint);
+ // force towards the body is negative, away from it is positive (compression/tension)
for(int i=0; i<3; i++){
stress[i]*=cp[i]>p1[i] ? 1. : -1.;
}
- bodyStats[id1].avgStress+=stress;
- bodyStats[id2].avgStress+=stress;
+ bodyStats[id1].sigma+=stress; bodyStats[id2].sigma+=stress;
+ bodyStats[id1].tau+=stress.Cross(cp-geom->se31.position);
+ bodyStats[id2].tau+=stress.Cross(cp-geom->se32.position);
bodyStats[id1].nLinks++; bodyStats[id2].nLinks++;
if(!phys->isCohesive) continue;
@@ -374,7 +376,8 @@
// add damaged contacts that have already been deleted
CpmMat* bpp=dynamic_cast<CpmMat*>(B->physicalParameters.get());
if(!bpp) continue;
- bpp->avgStress=bodyStats[id].avgStress; // /bodyStats[id].nLinks;
+ bpp->sigma=bodyStats[id].sigma;
+ bpp->tau=bodyStats[id].tau;
int cohLinksWhenever=bodyStats[id].nCohLinks+bpp->numBrokenCohesive;
if(cohLinksWhenever>0){
bpp->normDmg=(bodyStats[id].dmgSum+bpp->numBrokenCohesive)/cohLinksWhenever;
=== modified file 'pkg/dem/meta/ConcretePM.hpp'
--- pkg/dem/meta/ConcretePM.hpp 2009-10-16 06:33:55 +0000
+++ pkg/dem/meta/ConcretePM.hpp 2009-10-21 15:22:14 +0000
@@ -69,10 +69,10 @@
Real epsPlBroken;
//! sum of plastic strains normalized by number of contacts
Real normEpsPl;
- //! averaged stress on the particle
- Vector3r avgStress;
- CpmMat(): epsVolumetric(0.), numBrokenCohesive(0), numContacts(0), normDmg(0.), epsPlBroken(0.), normEpsPl(0.), avgStress(Vector3r::ZERO) {createIndex();};
- REGISTER_ATTRIBUTES(BodyMacroParameters, (epsVolumetric) (numBrokenCohesive) (numContacts) (normDmg) (epsPlBroken) (normEpsPl) (avgStress));
+ //! stresses on the particle
+ Vector3r sigma,tau;
+ CpmMat(): epsVolumetric(0.), numBrokenCohesive(0), numContacts(0), normDmg(0.), epsPlBroken(0.), normEpsPl(0.), sigma(Vector3r::ZERO), tau(Vector3r::ZERO) {createIndex();};
+ REGISTER_ATTRIBUTES(BodyMacroParameters, (epsVolumetric) (numBrokenCohesive) (numContacts) (normDmg) (epsPlBroken) (normEpsPl) (sigma) (tau));
REGISTER_CLASS_AND_BASE(CpmMat,BodyMacroParameters);
REGISTER_CLASS_INDEX(CpmMat,BodyMacroParameters);
};
@@ -317,7 +317,7 @@
#endif
class CpmStateUpdater: public PeriodicEngine {
- struct BodyStats{ int nCohLinks; int nLinks; Real dmgSum, epsPlSum; Vector3r avgStress; BodyStats(): nCohLinks(0), nLinks(0), dmgSum(0.), epsPlSum(0.), avgStress(Vector3r::ZERO) {} };
+ struct BodyStats{ int nCohLinks; int nLinks; Real dmgSum, epsPlSum; Vector3r sigma, tau; BodyStats(): nCohLinks(0), nLinks(0), dmgSum(0.), epsPlSum(0.), sigma(Vector3r::ZERO), tau(Vector3r::ZERO) {} };
public:
//! maximum damage over all contacts
static Real maxOmega;
=== modified file 'py/_eudoxos.cpp'
--- py/_eudoxos.cpp 2009-10-11 11:22:25 +0000
+++ py/_eudoxos.cpp 2009-10-21 15:22:14 +0000
@@ -87,6 +87,8 @@
}
BOOST_PYTHON_FUNCTION_OVERLOADS(velocityTowardsAxis_overloads,velocityTowardsAxis,3,5);
+// integrated in CpmStateUpdater
+#if 0
/* Compute σxx,σyy,σxy stresses over all spheres, in plane passing through given axis,
which will be coincident with the y axis in the 2d projection.
Not sure how much is this function useful... */
@@ -124,7 +126,7 @@
}
return ret;
}
-
+#endif
void particleConfinement(){
CpmStateUpdater::update();
}
@@ -133,7 +135,7 @@
BOOST_PYTHON_MODULE(_eudoxos){
def("velocityTowardsAxis",velocityTowardsAxis,velocityTowardsAxis_overloads(args("axisPoint","axisDirection","timeToAxis","subtractDist","perturbation")));
def("yieldSigmaTMagnitude",yieldSigmaTMagnitude);
- def("spiralSphereStresses2d",spiralSphereStresses2d,(python::arg("dH_dTheta"),python::arg("axis")=2));
+ // def("spiralSphereStresses2d",spiralSphereStresses2d,(python::arg("dH_dTheta"),python::arg("axis")=2));
def("particleConfinement",particleConfinement);
}
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2009-10-11 11:22:25 +0000
+++ py/_utils.cpp 2009-10-21 15:22:14 +0000
@@ -306,6 +306,24 @@
return inside;
}
+#if 0
+/* Compute convex hull of given points, given as python list of Vector2 */
+python::list convexHull2d(const python::list& pts){
+ size_t l=python::len(pts);
+ python::list ret;
+ std::list<Vector2r> pts2;
+ for(size_t i=0; i<l; i++){
+ pts2.push_back(python::extract<Vector2r>(pts[i]));
+ cerr<<*pts2.rbegin()<<endl;
+ }
+ ConvexHull2d ch2d(pts2);
+ vector<Vector2r> hull=ch2d();
+ FOREACH(const Vector2r& v, hull) ret.append(v);
+ cout<<pts2.size()<<" "<<hull.size()<<"@@@"<<endl;
+ return ret;
+}
+#endif
+
/* Compute area of convex hull when when taking (swept) spheres crossing the plane at coord, perpendicular to axis.
All spheres that touch the plane are projected as hexagons on their circumference to the plane.
@@ -318,7 +336,7 @@
if(axis<0 || axis>2) throw invalid_argument("Axis must be ∈ {0,1,2}");
const int ax1=(axis+1)%3, ax2=(axis+2)%3;
const Real sqrt3=sqrt(3);
- Vector2r mm(-10.,0),mx(0,0);
+ Vector2r mm,mx; int i=0;
FOREACH(const shared_ptr<Body>& b, *Omega::instance().getRootBody()->bodies){
Sphere* s=dynamic_cast<Sphere*>(b->geometricalModel.get());
if(!s) continue;
@@ -328,7 +346,7 @@
cloud.push_back(c+Vector2r(r,0)); cloud.push_back(c+Vector2r(-r,0));
cloud.push_back(c+Vector2r( r/2., sqrt3*r)); cloud.push_back(c+Vector2r( r/2.,-sqrt3*r));
cloud.push_back(c+Vector2r(-r/2., sqrt3*r)); cloud.push_back(c+Vector2r(-r/2.,-sqrt3*r));
- if(mm[0]==-10.){ mm=c, mx=c;}
+ if(i++==0){ mm=c, mx=c;}
mm=Vector2r(min(c[0]-r,mm[0]),min(c[1]-r,mm[1]));
mx=Vector2r(max(c[0]+r,mx[0]),max(c[1]+r,mx[1]));
}
@@ -397,6 +415,9 @@
def("ptInAABB",ptInAABB,"Return True/False whether the point (3-tuple) p is within box given by its min (3-tuple) and max (3-tuple) corners");
def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(args("axis","distFactor"),"Return list of ids for spheres (only) that are on extremal ends of the specimen along given axis; distFactor multiplies their radius so that sphere that do not touch the boundary coordinate can also be returned."));
def("approxSectionArea",approxSectionArea,"Compute area of convex hull when when taking (swept) spheres crossing the plane at coord, perpendicular to axis.");
+ #if 0
+ def("convexHull2d",convexHull2d,"Return 2d convex hull of list of 2d points, as list of polygon vertices.");
+ #endif
def("coordsAndDisplacements",coordsAndDisplacements,coordsAndDisplacements_overloads(args("AABB"),"Return tuple of 2 same-length lists for coordinates and displacements (coordinate minus reference coordinate) along given axis (1st arg); if the AABB=((x_min,y_min,z_min),(x_max,y_max,z_max)) box is given, only bodies within this box will be considered."));
def("setRefSe3",setRefSe3,"Set reference positions and orientation of all bodies equal to their current ones.");
def("interactionAnglesHistogram",interactionAnglesHistogram,interactionAnglesHistogram_overloads(args("axis","mask","bins","aabb")));
=== modified file 'py/plot.py'
--- py/plot.py 2009-08-03 17:53:26 +0000
+++ py/plot.py 2009-10-21 15:22:14 +0000
@@ -7,7 +7,7 @@
"""
import matplotlib
-matplotlib.use('TkAgg')
+#matplotlib.use('TkAgg')
#matplotlib.use('GTKCairo')
#matplotlib.use('QtAgg')
matplotlib.rc('axes',grid=True) # put grid in all figures
=== modified file 'py/utils.py'
--- py/utils.py 2009-09-19 19:41:52 +0000
+++ py/utils.py 2009-10-21 15:22:14 +0000
@@ -387,7 +387,7 @@
if len(env)>2: tableDesc=env[3]
o.tags['line']='l'+tableLine
# the empty '#' line to make line number 1-based
- ll=[l for l in ['#']+open(tableFile).readlines()]; values=ll[int(tableLine)].split('#')[0].split()
+ ll=[l for l in ['#']+open(tableFile).readlines()]; values=ll[int(tableLine)][:-1].split('#')[0].split()
for i in range(0,len(values)):
lineNo=int(tableLine)-1
while values[i]=='=':
=== modified file 'py/yadeWrapper/yadeWrapper.cpp'
--- py/yadeWrapper/yadeWrapper.cpp 2009-08-22 20:04:23 +0000
+++ py/yadeWrapper/yadeWrapper.cpp 2009-10-21 15:22:14 +0000
@@ -177,7 +177,7 @@
public:
pyTags(const shared_ptr<MetaBody> _mb): mb(_mb){}
const shared_ptr<MetaBody> mb;
- bool hasKey(string key){ FOREACH(string val, mb->tags){ if(algorithm::starts_with(val,key+"=")){ return true;} } return false; }
+ bool hasKey(const string& key){ FOREACH(string val, mb->tags){ if(algorithm::starts_with(val,key+"=")){ return true;} } return false; }
string getItem(string key){
FOREACH(string& val, mb->tags){
if(algorithm::starts_with(val,key+"=")){ string val1(val); algorithm::erase_head(val1,key.size()+1); algorithm::replace_all(val1,"~"," "); return val1;}
@@ -684,7 +684,8 @@
python::class_<pyTags>("TagsWrapper",python::init<pyTags&>())
.def("__getitem__",&pyTags::getItem)
.def("__setitem__",&pyTags::setItem)
- .def("keys",&pyTags::keys);
+ .def("keys",&pyTags::keys)
+ .def("has_key",&pyTags::hasKey);
python::class_<pyBodyContainer>("BodyContainer",python::init<pyBodyContainer&>())
.def("__getitem__",&pyBodyContainer::pyGetitem)
.def("__len__",&pyBodyContainer::length)
