yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04941
[Branch ~yade-dev/yade/trunk] Rev 2300: 1. Bugs in Cell (!!! no yet fixed: see scripts/test/peri8.py !!!)
------------------------------------------------------------
revno: 2300
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Tue 2010-06-29 16:32:03 +0200
message:
1. Bugs in Cell (!!! no yet fixed: see scripts/test/peri8.py !!!)
2. Add optional devirtualized functors (experimental)
3. Adjust Peri3dController for large strains and test (modulo the cell bug)
4. doc fixes
added:
scripts/test/peri8.py
modified:
SConstruct
core/Cell.hpp
core/Interaction.hpp
debian/control-template
doc/references.bib
doc/sphinx/conf.py
doc/sphinx/yadeSphinx.py
examples/concrete/uniax-post.py
examples/concrete/uniax.py
gui/qt3/GLViewer.cpp
pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp
pkg/common/Engine/Dispatcher/InteractionGeometryFunctor.hpp
pkg/common/Engine/GlobalEngine/FieldApplier.hpp
pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp
pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp
pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.hpp
pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp
pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.hpp
pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp
pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp
pkg/dem/PreProcessor/TriaxialTest.cpp
pkg/dem/PreProcessor/TriaxialTest.hpp
pkg/dem/meta/ConcretePM.cpp
pkg/dem/meta/ConcretePM.hpp
py/_extraDocs.py
py/pack/pack.py
py/plot.py
scripts/test/peri3d.py
yadeSCons.py
--
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 'SConstruct'
--- SConstruct 2010-06-08 21:22:07 +0000
+++ SConstruct 2010-06-29 14:32:03 +0000
@@ -139,11 +139,12 @@
BoolVariable('optimize','Turn on heavy optimizations',defOptions['optimize']),
ListVariable('exclude','Yade components that will not be built','none',names=['gui','extra','common','dem','lattice','snow']),
EnumVariable('PGO','Whether to "gen"erate or "use" Profile-Guided Optimization','',['','gen','use'],{'no':'','0':'','false':''},1),
- ListVariable('features','Optional features that are turned on','log4cxx,opengl,gts,openmp,vtk',names=['opengl','log4cxx','cgal','openmp','gts','vtk','python','gl2ps','never_use_this_one']),
+ ListVariable('features','Optional features that are turned on','log4cxx,opengl,gts,openmp,vtk',names=['opengl','log4cxx','cgal','openmp','gts','vtk','python','gl2ps','devirt-functors','never_use_this_one']),
('jobs','Number of jobs to run at the same time (same as -j, but saved)',2,None,int),
#('extraModules', 'Extra directories with their own SConscript files (must be in-tree) (whitespace separated)',None,None,Split),
('buildPrefix','Where to create build-[version][variant] directory for intermediary files','..'),
EnumVariable('linkStrategy','How to link plugins together',defOptions['linkStrategy'],['per-class','per-pkg[broken]','monolithic','static[broken]']),
+ ('hotPlugins','Files (without the .cpp extension) that will be compiled separately even in the monolithic build (use for those that you modify frequently); comma-separated.',''),
('chunkSize','Maximum files to compile in one translation unit when building plugins. (unlimited if <= 0)',7,None,int),
('version','Yade version (if not specified, guess will be attempted)',None),
('realVersion','Revision (usually bzr revision); guessed automatically unless specified',None),
@@ -577,7 +578,7 @@
import yadeSCons
allPlugs=yadeSCons.scanAllPlugins(None,feats=env['features'])
-buildPlugs=yadeSCons.getWantedPlugins(allPlugs,env['exclude'],env['features'],env['linkStrategy'])
+buildPlugs=yadeSCons.getWantedPlugins(allPlugs,env['exclude'],env['features'],env['linkStrategy'],env['hotPlugins'].split(','))
def linkPlugins(plugins):
"""Given list of plugins we need to link to, return list of real libraries that we should link to."""
ret=set()
=== modified file 'core/Cell.hpp'
--- core/Cell.hpp 2010-05-26 08:18:36 +0000
+++ core/Cell.hpp 2010-06-29 14:32:03 +0000
@@ -69,11 +69,11 @@
//! Apply shear on point.
Vector3r shearPt(const Vector3r& pt){ return _shearTrsf*pt; }
/*! Wrap point to inside the periodic cell; don't compute number of periods wrapped */
- Vector3r wrapPt(const Vector3r pt)const{
+ Vector3r wrapPt(const Vector3r& pt)const{
Vector3r ret; for(int i=0;i<3;i++) ret[i]=wrapNum(pt[i],_size[i]); return ret;
}
/*! Wrap point to inside the periodic cell; period will contain by how many cells it was wrapped. */
- Vector3r wrapPt(const Vector3r pt, Vector3i& period)const{
+ Vector3r wrapPt(const Vector3r& pt, Vector3i& period)const{
Vector3r ret; for(int i=0; i<3; i++){ ret[i]=wrapNum(pt[i],_size[i],period[i]); } return ret;
}
/*! Wrap number to interval 0…sz */
@@ -91,8 +91,11 @@
void setTrsf(const Matrix3r& m){ trsf=m; integrateAndUpdate(0); }
void postProcessAttributes(bool deserializing){ if(deserializing) integrateAndUpdate(0); }
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY
- (Cell,Serializable,"Parameters of periodic boundary conditions. Only applies if O.isPeriodic==True.",
+
+ Vector3r wrapShearedPt_py(const Vector3r& pt){ return wrapShearedPt(pt);}
+ Vector3r wrapPt_py(const Vector3r& pt){ return wrapPt(pt);}
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Cell,Serializable,"Parameters of periodic boundary conditions. Only applies if O.isPeriodic==True.",
((Vector3r,refSize,Vector3r(1,1,1),"[will be overridden below]"))
((Matrix3r,trsf,Matrix3r::Identity(),"[will be overridden below]"))
((Matrix3r,velGrad,Matrix3r::Zero(),"Velocity gradient of the transformation; used in NewtonIntegrator."))
@@ -104,6 +107,13 @@
/* accessors that ensure cache coherence */
.add_property("refSize",&Cell::getRefSize,&Cell::setRefSize,"Reference size of the cell.")
.add_property("trsf",&Cell::getTrsf,&Cell::setTrsf,"Transformation matrix of the cell.")
+ // debugging only
+ .def("wrap",&Cell::wrapShearedPt_py,"Transform an arbitrary point into a point in the reference cell")
+ .def("unshearPt",&Cell::unshearPt,"Apply inverse shear on the point (removes skew+rot of the cell)")
+ .def("shearPt",&Cell::shearPt,"Apply shear (cell skew+rot) on the point")
+ .def("wrapPt",&Cell::wrapPt_py,"Wrap point inside the reference cell, assuming the cell has no skew+rot.")
+ .def_readonly("shearTrsf",&Cell::_shearTrsf,"Current skew+rot transformation (no resize)")
+ .def_readonly("unshearTrsf",&Cell::_shearTrsf,"Inverse of the current skew+rot transformation (no resize)")
);
};
REGISTER_SERIALIZABLE(Cell);
=== modified file 'core/Interaction.hpp'
--- core/Interaction.hpp 2010-05-04 13:56:05 +0000
+++ core/Interaction.hpp 2010-06-29 14:32:03 +0000
@@ -7,7 +7,6 @@
#include"InteractionPhysics.hpp"
-
/////////////////////////////////
// FIXME - this is in wrong file!
//#include<boost/strong_typedef.hpp>
@@ -49,6 +48,11 @@
// Whether geometry dispatcher exists at all; this is different from !geom, since that can mean we haven't populated the cache yet.
// Therefore, geomExists must be initialized to true first (done in Interaction::reset() called from ctor).
bool geomExists;
+ #ifdef YADE_DEVIRT_FUNCTORS
+ // is a InteractionGeometryFunctor::StaticFuncPtr, but we would have to #include a file from pkg-common here
+ // cast at those few places instead, for now
+ void* geomPtr;
+ #endif
// shared_ptr's are initialized to NULLs automagically
shared_ptr<InteractionGeometryFunctor> geom;
shared_ptr<InteractionPhysicsFunctor> phys;
=== modified file 'debian/control-template'
--- debian/control-template 2010-06-16 16:39:17 +0000
+++ debian/control-template 2010-06-29 14:32:03 +0000
@@ -7,7 +7,6 @@
Package: yade@_VERSION@
Architecture: any
-Provides: yade@SNAPSHOT@
Depends: ${shlibs:Depends}, ${misc:Depends}, python-numpy, ipython, python-matplotlib, python-tk
Description: Platform for dynamical modeling.
Yet Another Dynamic Engine. etc.
@@ -16,7 +15,6 @@
Package: yade@_VERSION@-dbg
Architecture: any
-Provides: yade@SNAPSHOT@-dbg
Depends: ${shlibs:Depends}, ${misc:Depends}, python-numpy, ipython, gdb, python-matplotlib, python-tk
Description: Platform for dynamical modeling.
Yet Another Dynamic Engine. etc.
@@ -25,10 +23,30 @@
Package: yade@_VERSION@-doc
Architecture: all
-Provides: yade@SNAPSHOT@-doc
Recommends: yade@_VERSION@ | yade@_VERSION@-dbg
Depends:
Description: Platform for dynamical modeling.
Yet Another Dynamic Engine. etc.
.
This package contains examples, test scripts and documentation.
+
+Package: yade@SNAPSHOT@
+Architecture: all
+Depends: yade@_VERSION@
+Description: Platform for dynamical modeling.
+ This package depends on the latest
+ @'snapshot' if SNAPSHOT else 'stable release'@ of Yade.
+
+Package: yade@SNAPSHOT@-dbg
+Architecture: all
+Depends: yade@_VERSION@-dbg
+Description: Platform for dynamical modeling.
+ This package depends on the latest debug build of
+ @'snapshot' if SNAPSHOT else 'stable release'@ of Yade.
+
+Package: yade@SNAPSHOT@-doc
+Architecture: all
+Depends: yade@_VERSION@-doc
+Description: Platform for dynamical modeling.
+ This package depends on the latest documentation of
+ @'snapshot' if SNAPSHOT else 'stable release'@ of Yade.
=== modified file 'doc/references.bib'
--- doc/references.bib 2010-06-12 20:39:00 +0000
+++ doc/references.bib 2010-06-29 14:32:03 +0000
@@ -13,6 +13,16 @@
Thanks.
***************************
+@article{Lu1998,
+ author = {Ya Yan Lu},
+ title = {Computing the Logarithm of a Symmetric Positive Definite Matrix},
+ journal = {Appl. Numer. Math},
+ year = {1998},
+ volume = {26},
+ pages = {483--496},
+ doi={10.1016/S0168-9274(97)00103-7},
+ url={http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.759&rep=rep1&type=pdf}
+}
@inproceedings{Alonso2004,
title = {Micro-mechanical investigation of the granular ratcheting},
=== modified file 'doc/sphinx/conf.py'
--- doc/sphinx/conf.py 2010-06-18 13:02:01 +0000
+++ doc/sphinx/conf.py 2010-06-29 14:32:03 +0000
@@ -445,7 +445,7 @@
'''
-pngmath_latex_preamble=r'\usepackage[active]{preview}\usepackage{amsmath}\usepackage{amssymb}'+my_latex_preamble
+pngmath_latex_preamble=r'\usepackage[active]{preview}'+my_latex_preamble
pngmath_use_preview=True
=== modified file 'doc/sphinx/yadeSphinx.py'
--- doc/sphinx/yadeSphinx.py 2010-06-13 21:25:46 +0000
+++ doc/sphinx/yadeSphinx.py 2010-06-29 14:32:03 +0000
@@ -33,7 +33,7 @@
#
# don't forget to put the module in index.rst as well!
#
-mods={'eudoxos':['_eudoxos'],'log':[],'post2d':[],'pack':['_packSpheres','_packPredicates','_packObb','_packSpherePadder'],'plot':[],'timing':[],'utils':['_utils'],'ymport':[],'qt':['_qt'],'linterpolation':[]}
+mods={'export':[],'eudoxos':['_eudoxos'],'log':[],'post2d':[],'pack':['_packSpheres','_packPredicates','_packObb','_packSpherePadder'],'plot':[],'timing':[],'utils':['_utils'],'ymport':[],'qt':['_qt'],'linterpolation':[]}
#
# generate documentation, in alphabetical order
mm=mods.keys(); mm.sort()
@@ -60,6 +60,7 @@
docClasses.add(klass)
return ".. autoclass:: %s\n\t:members:\n\t:undoc-members:\n\n"%(klass) # \t:inherited-members:\n # \n\t:show-inheritance:
def autodocDerived(klass,doSections=True):
+ global docClasses
ret=''
if doSections: ret+='%s\n'%klass+'^'*100+'\n\n'
#ret+='.. inheritance-diagram:: %s\n\n'%klass
@@ -68,22 +69,29 @@
childs=list(yade.system.childClasses(klass));
childs.sort();
for k in childs:
- ret+=autodocClass(k)
+ if not k in docClasses: ret+=autodocClass(k)
return ret
def inheritanceDiagram(klass):
- "Generate simple inheritance graph for given classname *klass* using dot (http://www.graphviz.org) syntax."
+ """Generate simple inheritance graph for given classname *klass* using dot (http://www.graphviz.org) syntax.
+ Classes in the docClasses set (already documented) are framed differently and the inheritance tree stops there.
+ """
+ global docClasses
def mkUrl(c):
global writer
# useless, doesn't work in LaTeX anyway...
return ('#yade.wrapper.%s'%c if writer=='html' else '%%yade.wrapper#yade.wrapper.%s'%c) # HTML/LaTeX
- def mkNode(c): return '\t\t"%s" [shape="box",fontsize=8,style="setlinewidth(0.5)",height=0.2,URL="yade.wrapper.html#yade.wrapper.%s"];\n'%(c,c)
+ def mkNode(c,style='solid',fillcolor=None): return '\t\t"%s" [shape="box",fontsize=8,style="setlinewidth(0.5),%s",%sheight=0.2,URL="yade.wrapper.html#yade.wrapper.%s"];\n'%(c,style,'fillcolor=%s,'%fillcolor if fillcolor else '',c)
ret=".. graphviz::\n\n\tdigraph %s {\n\t\trankdir=RL;\n\t\tmargin=.2;\n"%klass+mkNode(klass)
childs=yade.system.childClasses(klass)
if len(childs)==0: return ''
for c in childs:
try:
base=eval(c).__bases__[0].__name__
- ret+=mkNode(c)
+ if base!=klass and base in docClasses:
+ continue # skip classes deriving from classes that are already documented
+ if c not in docClasses: ret+=mkNode(c)
+ else: # classes of which childs are documented elsewhere are marked specially
+ ret+=mkNode(c,style='filled,dashed',fillcolor='grey')
ret+='\t\t"%s" -> "%s" [arrowsize=0.5,style="setlinewidth(0.5)"];\n'%(c,base)
except NameError:
pass
@@ -91,11 +99,15 @@
return ret+'\t}\n\n'
-def sect(title,text,tops):
- return title+'\n'+100*'-'+'\n\n'+text+'\n\n'+'\n\n'.join([autodocDerived(top,doSections=(len(tops)>1)) for top in tops])+'\n'
+def sect(title,text,tops,reverse=False):
+ subsects=[autodocDerived(top,doSections=(len(tops)>1)) for top in tops]
+ if reverse: subsects.reverse()
+ return title+'\n'+100*'-'+'\n\n'+text+'\n\n'+'\n\n'.join(subsects)+'\n'
def genWrapperRst():
+ global docClasses
+ docClasses=set() # reset globals
wrapper=file('yade.wrapper.rst','w')
wrapper.write(""".. _yade.wrapper
Class reference (yade.wrapper module)
@@ -111,7 +123,7 @@
"""+
sect('Bodies','',['Body','Shape','State','Material','Bound'])+
sect('Interactions','',['Interaction','InteractionGeometry','InteractionPhysics'])+
- sect('Global engines','',['GlobalEngine'])+
+ sect('Global engines','',['FieldApplier','Collider','BoundaryController','GlobalEngine'],reverse=True)+
sect('Partial engines','',['PartialEngine'])+
sect('Bounding volume creation','',['BoundFunctor','BoundDispatcher'])+
sect('Interaction Geometry creation','',['InteractionGeometryFunctor','InteractionGeometryDispatcher'])+
=== modified file 'examples/concrete/uniax-post.py'
--- examples/concrete/uniax-post.py 2010-05-29 20:52:10 +0000
+++ examples/concrete/uniax-post.py 2010-06-29 14:32:03 +0000
@@ -6,7 +6,7 @@
from yade import post2d
import pylab # the matlab-like interface of matplotlib
-loadFile='/tmp/uniax-tension.xml.bz2'
+loadFile='/tmp/uniax-tension.yade.gz'
if not os.path.exists(loadFile): raise RuntimeError("Run uniax.py first so that %s is created"%loadFile)
O.load(loadFile)
=== modified file 'examples/concrete/uniax.py'
--- examples/concrete/uniax.py 2010-05-29 20:52:10 +0000
+++ examples/concrete/uniax.py 2010-06-29 14:32:03 +0000
@@ -53,7 +53,7 @@
dtSafety=.8,
damping=0.4,
strainRateTension=.05,
- strainRateCompression=1,
+ strainRateCompression=.5,
setSpeeds=True,
# 1=tension, 2=compression (ANDed; 3=both)
doModes=3,
@@ -97,7 +97,7 @@
NewtonIntegrator(damping=damping,label='damper'),
CpmStateUpdater(realPeriod=1),
UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
- PeriodicPythonRunner(virtPeriod=3e-5/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
+ PeriodicPythonRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
PeriodicPythonRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
]
#O.miscParams=[Gl1_CpmPhys(dmgLabel=False,colorStrain=False,epsNLabel=False,epsT=False,epsTAxes=False,normal=False,contactLine=True)]
@@ -117,7 +117,7 @@
global mode
print "init"
if O.iter>0:
- O.pause(); O.wait();
+ O.wait();
O.loadTmp('initial')
print "Reversing plot data"; plot.reverseData()
strainer.strainRate=abs(strainRateTension) if mode=='tension' else -abs(strainRateCompression)
@@ -147,9 +147,9 @@
if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if isoPrestress==0 else 5e-2):
if mode=='tension' and doModes & 2: # only if compression is enabled
mode='compression'
- O.save('/tmp/uniax-tension.xml.bz2')
- print "Saved /tmp/uniax-tension.xml.bz2 (for use with interaction-histogram.py and uniax-post.py)"
- print "Damaged, switching to compression... "
+ O.save('/tmp/uniax-tension.yade.gz')
+ print "Saved /tmp/uniax-tension.yade.gz (for use with interaction-histogram.py and uniax-post.py)"
+ print "Damaged, switching to compression... "; O.pause()
# important! initTest must be launched in a separate thread;
# otherwise O.load would wait for the iteration to finish,
# but it would wait for initTest to return and deadlock would result
@@ -162,8 +162,8 @@
title=O.tags['description'] if 'description' in O.tags.keys() else O.tags['params']
print 'gnuplot',plot.saveGnuplot(O.tags['id'],title=title)
print 'Bye.'
- # O.pause()
- sys.exit(0)
+ O.pause()
+ #sys.exit(0)
def addPlotData():
yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
=== modified file 'gui/qt3/GLViewer.cpp'
--- gui/qt3/GLViewer.cpp 2010-06-08 22:25:00 +0000
+++ gui/qt3/GLViewer.cpp 2010-06-29 14:32:03 +0000
@@ -87,6 +87,7 @@
//xyPlaneConstraint->setRotationConstraint(qglviewer::AxisPlaneConstraint::FORBIDDEN,qglviewer::Vec(0,0,1));
manipulatedFrame()->setConstraint(NULL);
+ setKeyDescription(Qt::Key_A,"Toggle visibility of global axes.");
setKeyDescription(Qt::Key_C,"Set scene center so that all bodies are visible; if a body is selected, center around this body.");
setKeyDescription(Qt::Key_C & Qt::ALT,"Set scene center to median body position");
setKeyDescription(Qt::Key_D,"Toggle time display mask");
@@ -228,6 +229,7 @@
if(false){}
/* special keys: Escape and Space */
//else if(e->key()==Qt::Key_F9 || e->key()==Qt::Key_F10 || e->key()==Qt::Key_F11 || e->key()==Qt::Key_F12){ YadeQtMainWindow::self->closeView(this); }
+ else if(e->key()==Qt::Key_A){ toggleAxisIsDrawn(); return; }
else if(e->key()==Qt::Key_Escape){
if(!isManipulating()){ /* CRASH, deleting this: YadeQtMainWindow::self->closeView(this); */ return; }
else { resetManipulation(); displayMessage("Manipulating scene."); }
=== modified file 'pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp'
--- pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2010-06-06 10:18:07 +0000
+++ pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2010-06-29 14:32:03 +0000
@@ -59,7 +59,7 @@
for(long i=0; i<size; i++){
const shared_ptr<Interaction>& I=(*scene->interactions)[i];
#else
- FOREACH(shared_ptr<Interaction> I, *scene->interactions){
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
#endif
if(removeUnseenIntrs && !I->isReal() && I->iterLastSeen<scene->currentIteration) {
eraseAfterLoop(I->getId1(),I->getId2());
@@ -80,13 +80,19 @@
// InteractionGeometryDispatcher
if(!I->functorCache.geom || !I->functorCache.phys){
I->functorCache.geom=geomDispatcher->getFunctor2D(b1_->shape,b2_->shape,swap);
- // returns NULL ptr if no functor exists; remember that and shortcut
- if(!I->functorCache.geom) { I->functorCache.geomExists=false; continue; }
+ #ifdef YADE_DEVIRT_FUNCTORS
+ if(I->functorCache.geom){ I->functorCache.geomPtr=I->functorCache.geom->getStaticFuncPtr(); /* cerr<<"["<<I->functorCache.geomPtr<<"]"; */ }
+ else
+ #else
+ if(!I->functorCache.geom)
+ #endif
+ // returns NULL ptr if no functor exists; remember that and shortcut
+ {I->functorCache.geomExists=false; continue; }
}
// arguments for the geom functor are in the reverse order (dispatcher would normally call goReverse).
// we don't remember the fact that is reverse, so we swap bodies within the interaction
// and can call go in all cases
- if(swap){I->swapOrder();}
+ if(swap){I->swapOrder(); }
// body pointers must be updated, in case we swapped
const shared_ptr<Body>& b1=Body::byId(I->getId1(),scene);
const shared_ptr<Body>& b2=Body::byId(I->getId2(),scene);
@@ -94,21 +100,22 @@
assert(I->functorCache.geom);
bool wasReal=I->isReal();
bool geomCreated;
- if(!scene->isPeriodic) geomCreated=I->functorCache.geom->go(b1->shape,b2->shape, *b1->state, *b2->state, Vector3r::Zero(), /*force*/false, I);
- else{ // handle periodicity
+ if(!scene->isPeriodic){
+ #ifdef YADE_DEVIRT_FUNCTORS
+ geomCreated=(*((InteractionGeometryFunctor::StaticFuncPtr)I->functorCache.geomPtr))(I->functorCache.geom.get(),b1->shape,b2->shape, *b1->state, *b2->state, Vector3r::Zero(), /*force*/false, I);
+ #else
+ geomCreated=I->functorCache.geom->go(b1->shape,b2->shape, *b1->state, *b2->state, Vector3r::Zero(), /*force*/false, I);
+ #endif
+ } else { // handle periodicity
Vector3r shift2(I->cellDist[0]*cellSize[0],I->cellDist[1]*cellSize[1],I->cellDist[2]*cellSize[2]);
-
// in sheared cell, apply shear on the mutual position as well
shift2=scene->cell->shearPt(shift2);
- /* suggested change to avoid one matrix multiplication (with Hsize updated in cell ofc), I'll make the change cleanly if ok.
- Same sorts of simplifications are possible in many places. Just putting one example here.
- Hsize will contain colums with transformed base vectors
- Matrix3r Hsize(scene->cell->refSize[0],scene->cell->refSize[1],scene->cell->refSize[2]); Hsize=scene->cell->trsf*Hsize;
- Vector3r shift3((Real) I->cellDist[0]*Hsize.GetColumn(0)+(Real)I->cellDist[1]*Hsize.GetColumn(1)+(Real)I->cellDist[2]*Hsize.GetColumn(2));
- if ((Omega::instance().getCurrentIteration() % 100 == 0)) LOG_DEBUG(shift2 << " vs " << shift3);
- */
-
- geomCreated=I->functorCache.geom->go(b1->shape,b2->shape,*b1->state,*b2->state,shift2,/*force*/false,I);
+ #ifdef YADE_DEVIRT_FUNCTORS
+ // cast back from void* first
+ geomCreated=(*((InteractionGeometryFunctor::StaticFuncPtr)I->functorCache.geomPtr))(I->functorCache.geom.get(),b1->shape,b2->shape,*b1->state,*b2->state,shift2,/*force*/false,I);
+ #else
+ geomCreated=I->functorCache.geom->go(b1->shape,b2->shape,*b1->state,*b2->state,shift2,/*force*/false,I);
+ #endif
}
if(!geomCreated){
if(wasReal) LOG_WARN("InteractionGeometryFunctor returned false on existing interaction!");
=== modified file 'pkg/common/Engine/Dispatcher/InteractionGeometryFunctor.hpp'
--- pkg/common/Engine/Dispatcher/InteractionGeometryFunctor.hpp 2010-02-09 20:22:04 +0000
+++ pkg/common/Engine/Dispatcher/InteractionGeometryFunctor.hpp 2010-06-29 14:32:03 +0000
@@ -46,6 +46,12 @@
>
{
public: virtual ~InteractionGeometryFunctor();
+ #ifdef YADE_DEVIRT_FUNCTORS
+ // type of the pointer to devirtualized functor (static method taking the functor instance as the first argument)
+ typedef bool(*StaticFuncPtr)(InteractionGeometryFunctor*, const shared_ptr<Shape>&, const shared_ptr<Shape>&, const State&, const State&, const Vector3r&, const bool&, const shared_ptr<Interaction>&);
+ // return devirtualized functor (static method); must be overridden in derived classes
+ virtual void* getStaticFuncPtr(){ throw runtime_error(("InteractionGeometryFunctor::getStaticFuncPtr() not overridden in class "+getClassName()+".").c_str()); }
+ #endif
YADE_CLASS_BASE_DOC(InteractionGeometryFunctor,Functor,"Functor for creating/updating :yref:`Interaction::interactionGeometry` objects.");
};
REGISTER_SERIALIZABLE(InteractionGeometryFunctor);
=== modified file 'pkg/common/Engine/GlobalEngine/FieldApplier.hpp'
--- pkg/common/Engine/GlobalEngine/FieldApplier.hpp 2010-04-14 08:09:02 +0000
+++ pkg/common/Engine/GlobalEngine/FieldApplier.hpp 2010-06-29 14:32:03 +0000
@@ -2,7 +2,7 @@
#include<yade/core/GlobalEngine.hpp>
class FieldApplier: public GlobalEngine{
virtual void action();
- YADE_CLASS_BASE_DOC(FieldApplier,GlobalEngine,"Base for engines controlling boundary conditions of simulations. Not to be used directly.");
+ YADE_CLASS_BASE_DOC(FieldApplier,GlobalEngine,"Base for engines applying force files on particles. Not to be used directly.");
};
REGISTER_SERIALIZABLE(FieldApplier);
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-06-21 17:31:31 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-06-29 14:32:03 +0000
@@ -200,7 +200,7 @@
Vector3r ScGeom::relRotVector() const{
Quaternionr relOri12=ori1.Conjugate()*ori2;
Quaternionr oriDiff=initRelOri12.Conjugate()*relOri12;
- Vector3r axis; Real angle;
+ Vector3r axis; Real angle;
oriDiff.ToAxisAngle(axis,angle);
if(angle>Mathr::PI)angle-=Mathr::TWO_PI;
return angle*axis;
=== modified file 'pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-05-04 13:56:05 +0000
+++ pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-06-29 14:32:03 +0000
@@ -17,6 +17,10 @@
+#ifdef YADE_DEVIRT_FUNCTORS
+bool Ig2_Box_Sphere_ScGeom::go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){ throw runtime_error("Do not call Ig2_Box_Sphere_ScGeom::go, use getStaticFunctorPtr and call that function instead."); }
+bool Ig2_Box_Sphere_ScGeom::goStatic(InteractionGeometryFunctor* self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
+#else
bool Ig2_Box_Sphere_ScGeom::go(
const shared_ptr<Shape>& cm1,
const shared_ptr<Shape>& cm2,
@@ -26,6 +30,7 @@
const bool& force,
const shared_ptr<Interaction>& c)
{
+#endif
const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;
bool inside=true;
=== modified file 'pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.hpp'
--- pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.hpp 2010-02-09 16:50:30 +0000
+++ pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.hpp 2010-06-29 14:32:03 +0000
@@ -30,6 +30,12 @@
const Vector3r& shift2,
const bool& force,
const shared_ptr<Interaction>& c);
+
+ #ifdef YADE_DEVIRT_FUNCTORS
+ void* getStaticFuncPtr(){ return (void*)&Ig2_Box_Sphere_ScGeom::goStatic; }
+ static bool goStatic(InteractionGeometryFunctor* self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& se32, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ #endif
+
YADE_CLASS_BASE_DOC(Ig2_Box_Sphere_ScGeom,InteractionGeometryFunctor,"Create an interaction geometry :yref:`ScGeom` from :yref:`Box` and :yref:`Sphere`")
FUNCTOR2D(Box,Sphere);
DEFINE_FUNCTOR_ORDER_2D(Box,Sphere);
=== modified file 'pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-05-24 19:56:22 +0000
+++ pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-06-29 14:32:03 +0000
@@ -10,11 +10,18 @@
#include<yade/lib-base/Math.hpp>
#include<yade/core/Omega.hpp>
+#ifdef YADE_DEVIRT_FUNCTORS
+bool Ig2_Sphere_Sphere_ScGeom::go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){ throw runtime_error("Do not call Ig2_Sphere_Sphere_ScGeom::go, use getStaticFunctorPtr and call that function instead."); }
+bool Ig2_Sphere_Sphere_ScGeom::goStatic(InteractionGeometryFunctor* _self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
+ const Ig2_Sphere_Sphere_ScGeom* self=static_cast<Ig2_Sphere_Sphere_ScGeom*>(_self);
+ const Real& interactionDetectionFactor=self->interactionDetectionFactor;
+#else
bool Ig2_Sphere_Sphere_ScGeom::go( const shared_ptr<Shape>& cm1,
const shared_ptr<Shape>& cm2,
const State& state1, const State& state2, const Vector3r& shift2, const bool& force,
const shared_ptr<Interaction>& c)
{
+#endif
const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;
Sphere *s1=static_cast<Sphere*>(cm1.get()), *s2=static_cast<Sphere*>(cm2.get());
@@ -44,6 +51,7 @@
return false;
}
+
bool Ig2_Sphere_Sphere_ScGeom::goReverse( const shared_ptr<Shape>& cm1,
const shared_ptr<Shape>& cm2,
const State& state1,
=== modified file 'pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.hpp'
--- pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.hpp 2010-02-18 21:08:54 +0000
+++ pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.hpp 2010-06-29 14:32:03 +0000
@@ -12,8 +12,12 @@
class Ig2_Sphere_Sphere_ScGeom: public InteractionGeometryFunctor{
public:
- virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& se32, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
- virtual bool goReverse( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& se32, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ #ifdef YADE_DEVIRT_FUNCTORS
+ void* getStaticFuncPtr(){ return (void*)&Ig2_Sphere_Sphere_ScGeom::goStatic; }
+ static bool goStatic(InteractionGeometryFunctor* self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& se32, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ #endif
YADE_CLASS_BASE_DOC_ATTRS(Ig2_Sphere_Sphere_ScGeom,InteractionGeometryFunctor,"Create/update a :yref:`ScGeom` instance representing intersection of two :yref:`Spheres<Sphere>`.",
((Real,interactionDetectionFactor,1,"Enlarge both radii by this factor (if >1), to permit creation of distant interactions.\n\nInteractionGeometry will be computed when interactionDetectionFactor*(rad1+rad2) > distance.\n\n.. note::\n\t This parameter is functionally coupled with :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor`, which will create larger bounding boxes and should be of the same value.\n\n.. warning::\n\tFunctionally equal class :yref:`Ig2_Sphere_Sphere_Dem3DofGeom` (which creates :yref:`Dem3DofGeom` rather than :yref:`ScGeom`) calls this parameter :yref:`distFactor<Ig2_Sphere_Sphere_Dem3DofGeom::distFactor>`, but its semantics is *different* in some aspects."))
);
=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2010-06-08 22:25:00 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2010-06-29 14:32:03 +0000
@@ -258,16 +258,23 @@
Matrix3r rot,nonrot; //nonrot=skew+normal deformation
const Matrix3r& trsf=scene->cell->trsf;
Eigen::SVD<Matrix3r>(trsf).computeUnitaryPositive(&rot,&nonrot);
- // FIXME: should be rot*(log(nonrot)), but matrix logarithm is not yet implemented in eigen
- strain=rot*(nonrot-Matrix3r::Identity());
+ // compute matrix logarithm (see documentation)
+ Eigen::SelfAdjointEigenSolver<Matrix3r> eigDec(nonrot);
+ strain=rot*eigDec.eigenvectors()*eigDec.eigenvalues().cwise().log().asDiagonal()*eigDec.eigenvectors().transpose();
+ #if 0
+ // small strains only
+ strain=rot*(nonrot-Matrix3r::Identity());
+ #endif
LOG_TRACE("Updated strain value\n"<<strain);
/* stress and stiffness
*/
K=Matrix6r::Zero();
stress=Matrix3r::Zero();
const Real volume=scene->cell->trsf.determinant()*scene->cell->refSize[0]*scene->cell->refSize[1]*scene->cell->refSize[2];
+ int nIntr=0;
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
+ nIntr++;
Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());
NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->interactionPhysics.get());
// not clear whether this should be the reference or the current distance
@@ -299,7 +306,7 @@
stress/=volume;
K/=volume;
for(int p=0; p<6; p++)for(int q=p+1;q<6;q++) K(q,p)=K(p,q);
- LOG_TRACE("Updated stress\n"<<stress);
+ LOG_TRACE("Updated stress (from "<<nIntr<<" interactions)\n"<<stress);
LOG_TRACE("Updated stiffness tensor\n"<<K);
}
@@ -308,6 +315,7 @@
update();
typedef Eigen::Matrix<Real,Eigen::Dynamic,Eigen::Dynamic> MatrixXr;
typedef Eigen::Matrix<Real,Eigen::Dynamic,1> VectorXr;
+ const Real& dt=scene->dt;
/*
sigma = K * eps
@@ -336,12 +344,12 @@
for(int j=0; j<6; j++){
// prescribed stress, i.e. un-prescribed strain
if(stressMask&(1<<j)){
- sigU[jU]=goal(mapI[j],mapJ[j]);
+ sigU[jU]=(1/dt)*(goal(mapI[j],mapJ[j])-stress(mapI[j],mapJ[j]));
int iU=0;
for(int i=0;i<6;i++){ if((stressMask&(1<<i))) Kuu(iU++,jU)=K(i,j);}
jU++;
} else {
- epsP[jP]=(j<3?1.:2.)*goal(mapI[j],mapJ[j]); /* multiply tensor shear by 2, see http://en.wikiversity.org/wiki/Introduction_to_Elasticity/Constitutive_relations */
+ epsP[jP]=(1/dt)*(j<3?1.:2.)*(goal(mapI[j],mapJ[j])-strain(mapI[j],mapJ[j])); /* multiply tensor shear by 2, see http://en.wikiversity.org/wiki/Introduction_to_Elasticity/Constitutive_relations */
int iP=0;
for(int i=0;i<6;i++){ if((stressMask&(1<<i))) Kup(iP++,jP)=K(i,j);}
jP++;
@@ -354,9 +362,9 @@
// FIXME: find a better way for this than determinant (expensive for larger matrices)
if(Kuu.rows()*Kuu.cols()>0 && Kuu.determinant()<1e-20){ Kuu+=MatrixXr(Kuu.cols(),Kuu.rows()).setIdentity(); LOG_TRACE("Kuu after sanitization\n"<<Kuu); }
epsU=Kuu.inverse()*(sigU-Kup*epsP);
- // assemble strain tensor (as matrix), from prescribed values epsP and computed values epsU
+ // assemble strain rate tensor (as matrix), from prescribed values epsP and computed values epsU
jU=0; jP=0;
- Matrix3r eps;
+ Matrix3r eps; // rate!
for(int j=0; j<6; j++){
if(stressMask&(1<<j)) eps(mapI[j],mapJ[j])=epsU[jU++];
else eps(mapI[j],mapJ[j])=(j<3?1.:.5)*epsP[jP++]; /* multiply shear components back by 1/2 when converting from Voigt vector back to tensor */
@@ -364,7 +372,7 @@
eps(2,1)=eps(1,2); eps(0,2)=eps(2,0); eps(1,0)=eps(0,1);
Matrix3r& velGrad=scene->cell->velGrad;
// rate of (goal strain - current strain)
- velGrad=(eps-strain)/scene->dt;
+ velGrad=eps;
Real mx=max(abs(velGrad.minCoeff()),abs(velGrad.maxCoeff()));
if(mx>abs(maxStrainRate)) velGrad*=abs(maxStrainRate)/mx;
LOG_TRACE("epsU=\n"<<epsU<<"\neps=\n"<<"\nabs(maxCoeff)="<<mx<<"\nvelGrad=\n"<<velGrad);
=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp 2010-06-07 17:09:08 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp 2010-06-29 14:32:03 +0000
@@ -41,6 +41,7 @@
See scripts/test/periodic-triax.py for a simple example.
*/
+
class PeriTriaxController: public BoundaryController{
public:
virtual void action();
@@ -85,7 +86,7 @@
YADE_CLASS_BASE_DOC_ATTRS(Peri3dController,BoundaryController,"Experimental controller of full strain/stress tensors on periodic cell. Detailed documentation is in py/_extraDocs.py.",
((Matrix3r,strain,Matrix3r::Zero(),"Current deformation tensor |yupdate|"))
((Matrix3r,stress,Matrix3r::Zero(),"Current stress tensor |yupdate|"))
- ((Matrix3r,goal,Matrix3r::Zero(),"Goal state."))
+ ((Matrix3r,goal,Matrix3r::Zero(),"Goal state; only the upper triangular matrix is considered; each component is either prescribed stress or strain, depending on :yref:`stressMask<Peri3dController.stressMask>`."))
((int,stressMask,((void)"all strains",0),"mask determining whether components of :yref:`goal<Peri3dController.goal>` are strain (0) or stress (1). The order is 00,11,22,12,02,01 from the least significant bit. (e.g. 0b000011 is stress 00 and stress 11)."))
((Real,maxStrainRate,1,"Maximum absolute value of strain rate (both normal and shear components of :yref:`Cell.velGrad`)"))
// not yet used
=== modified file 'pkg/dem/PreProcessor/TriaxialTest.cpp'
--- pkg/dem/PreProcessor/TriaxialTest.cpp 2010-06-11 22:22:13 +0000
+++ pkg/dem/PreProcessor/TriaxialTest.cpp 2010-06-29 14:32:03 +0000
@@ -112,8 +112,7 @@
if(thickness<0) thickness=radiusMean;
if(facetWalls || wallWalls) thickness=0;
- if(boxWalls)
- {
+
// bottom box
Vector3r center = Vector3r(
(lowerCorner[0]+upperCorner[0])/2,
@@ -205,7 +204,6 @@
triaxialcompressionEngine->wall_front_id = body->getId();
//triaxialStateRecorder->wall_front_id = body->getId();
}
- }
size_t imax=sphere_pack.pack.size();
for(size_t i=0; i<imax; i++){
const SpherePack::Sph& sp(sphere_pack.pack[i]);
@@ -345,8 +343,8 @@
triaxialcompressionEngine->Key = Key;
triaxialcompressionEngine->noFiles=noFiles;
triaxialcompressionEngine->frictionAngleDegree = sphereFrictionDeg;
- triaxialcompressionEngine->fixedPorosity = fixedPorosity;
- triaxialcompressionEngine->fixedPoroCompaction = fixedPoroCompaction;
+ triaxialcompressionEngine->fixedPoroCompaction = false;
+ triaxialcompressionEngine->fixedPorosity=1;
// recording global stress
if(recordIntervalIter>0 && !noFiles){
triaxialStateRecorder = shared_ptr<TriaxialStateRecorder>(new TriaxialStateRecorder);
=== modified file 'pkg/dem/PreProcessor/TriaxialTest.hpp'
--- pkg/dem/PreProcessor/TriaxialTest.hpp 2010-05-21 13:46:33 +0000
+++ pkg/dem/PreProcessor/TriaxialTest.hpp 2010-06-29 14:32:03 +0000
@@ -84,16 +84,13 @@
((string,WallStressRecordFile,"./WallStresses"+Key,""))
((bool,internalCompaction,false,"flag for choosing between moving boundaries or increasing particles sizes during the compaction stage."))
((bool,biaxial2dTest,false,"FIXME : what is that?"))
- ((bool,fixedPoroCompaction,false,"flag to choose an isotropic compaction until a fixed porosity choosing a same translation speed for the six walls"))
((bool,autoCompressionActivation,true,"Do we just want to generate a stable packing under isotropic pressure (false) or do we want the triaxial loading to start automatically right after compaction stage (true)?"))
((bool,autoUnload,true,"auto adjust the isotropic stress state from :yref:`TriaxialTest::sigmaIsoCompaction` to :yref:`TriaxialTest::sigmaLateralConfinement` if they have different values. See docs for :yref:`TriaxialCompressionEngine::autoUnload`"))
((bool,autoStopSimulation,false,"freeze the simulation when conditions are reached (don't activate this if you want to be able to run/stop from Qt GUI)"))
((bool,noFiles,false,"Do not create any files during run (.xml, .spheres, wall stress records)"))
((bool,facetWalls,false,"Use facets for boundaries (not tested)"))
((bool,wallWalls,false,"Use walls for boundaries (not tested)"))
- ((bool,boxWalls,true,"Use boxes for boundaries (recommended)."))
- ((Real,fixedPorosity,1,"FIXME : what is that?"))
((Real,thickness,0.001,"thickness of boundaries. It is arbitrary and should have no effect"))
((Real,maxMultiplier,1.01,"max multiplier of diameters during internal compaction (initial fast increase)"))
((Real,finalMaxMultiplier,1.001,"max multiplier of diameters during internal compaction (secondary precise adjustment)"))
=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp 2010-06-08 22:25:00 +0000
+++ pkg/dem/meta/ConcretePM.cpp 2010-06-29 14:32:03 +0000
@@ -250,31 +250,30 @@
bool Gl1_CpmPhys::epsT=false;
bool Gl1_CpmPhys::epsTAxes=false;
bool Gl1_CpmPhys::normal=false;
- bool Gl1_CpmPhys::colorStrain=false;
+ Real Gl1_CpmPhys::colorStrainRatio=-1;
void Gl1_CpmPhys::go(const shared_ptr<InteractionPhysics>& ip, const shared_ptr<Interaction>& i, const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool wireFrame){
const shared_ptr<CpmPhys>& BC=static_pointer_cast<CpmPhys>(ip);
const shared_ptr<Dem3DofGeom>& geom=YADE_PTR_CAST<Dem3DofGeom>(i->interactionGeometry);
+ // FIXME: get the scene for periodicity; ugly!
+ Scene* scene=Omega::instance().getScene().get();
//Vector3r lineColor(BC->omega,1-BC->omega,0.0); /* damaged links red, undamaged green */
Vector3r lineColor=Shop::scalarOnColorScale(1.-BC->relResidualStrength);
- if(colorStrain) lineColor=Vector3r(
- min((Real)1.,max((Real)0.,-BC->epsTrans/BC->epsCrackOnset)),
- min((Real)1.,max((Real)0.,BC->epsTrans/BC->epsCrackOnset)),
- min((Real)1.,max((Real)0.,abs(BC->epsTrans)/BC->epsCrackOnset-1)));
+ if(colorStrainRatio>0) lineColor=Shop::scalarOnColorScale(BC->epsN/(BC->epsCrackOnset*colorStrainRatio));
// FIXME: should be computed by the renderer; for now, use the real values
- const Se3r& dispSe31=b1->state->se3;
- const Se3r& dispSe32=b2->state->se3;
+ Vector3r dispPt1=geom->se31.position, dispPt2=geom->se32.position;
+ if(scene->isPeriodic){ dispPt1=scene->cell->wrapShearedPt(dispPt1); dispPt2=dispPt1+(geom->se32.position-geom->se31.position); }
- if(contactLine) GLUtils::GLDrawLine(dispSe31.position,dispSe32.position,lineColor);
- if(dmgLabel){ GLUtils::GLDrawNum(BC->omega,0.5*(dispSe32.position+dispSe32.position),lineColor); }
- else if(epsNLabel){ GLUtils::GLDrawNum(BC->epsN,0.5*(dispSe31.position+dispSe32.position),lineColor); }
+ if(contactLine) GLUtils::GLDrawLine(dispPt1,dispPt2,lineColor);
+ if(dmgLabel){ GLUtils::GLDrawNum(BC->omega,0.5*(dispPt1+dispPt2),lineColor); }
+ else if(epsNLabel){ GLUtils::GLDrawNum(BC->epsN,0.5*(dispPt1+dispPt2),lineColor); }
if(BC->omega>0 && dmgPlane){
Real halfSize=sqrt(1-BC->relResidualStrength)*.5*.705*sqrt(BC->crossSection);
- Vector3r midPt=.5*Vector3r(dispSe31.position+dispSe32.position);
+ Vector3r midPt=.5*Vector3r(dispPt1+dispPt2);
glDisable(GL_CULL_FACE);
glPushMatrix();
glTranslatev(midPt);
@@ -293,7 +292,8 @@
glPopMatrix();
}
- const Vector3r& cp=static_pointer_cast<Dem3DofGeom>(i->interactionGeometry)->contactPoint;
+ Vector3r cp=static_pointer_cast<Dem3DofGeom>(i->interactionGeometry)->contactPoint;
+ if(scene->isPeriodic){cp=scene->cell->wrapShearedPt(cp);}
if(epsT){
Real maxShear=(BC->undamagedCohesion-BC->sigmaN*BC->tanFrictionAngle)/BC->G;
Real relShear=BC->epsT.norm()/maxShear;
=== modified file 'pkg/dem/meta/ConcretePM.hpp'
--- pkg/dem/meta/ConcretePM.hpp 2010-04-25 13:18:11 +0000
+++ pkg/dem/meta/ConcretePM.hpp 2010-06-29 14:32:03 +0000
@@ -213,7 +213,7 @@
((bool,epsT,false,"Show shear strain "))
((bool,epsTAxes,false,"Show axes of shear plane "))
((bool,normal,false,"Show contact normal"))
- ((bool,colorStrain,false,"Render elements using transverse strain (relative to max elastic strain) [no longer used]"))
+ ((Real,colorStrainRatio,-1,"If positive, set the interaction (wire) color based on $\\eps_N$ normalized by $\\eps_0$ × `colorStrainRatio` ($\\eps_0$=:yref:`CpmPhys.epsCrackOnset`). Otherwise, color based on the residual strength."))
((bool,epsNLabel,false,"Numerically show normal strain"))
);
};
=== modified file 'py/_extraDocs.py'
--- py/_extraDocs.py 2010-06-03 21:01:34 +0000
+++ py/_extraDocs.py 2010-06-29 14:32:03 +0000
@@ -72,35 +72,37 @@
K_{ijkl}=\frac{1}{V}\sum_{c}{l}^2\left[K_N n_i n_j n_k n_l+K_T\left(\frac{n_j n_k \delta_{il}+n_j n_l\delta_{ik}+n_i n_k\delta_{jl}+n_i n_l\delta_{jk}}{4}-n_i n_j n_kn_l\right)\right].
-Global strain tensor $\tens{\eps}$ is symmetric and has 6 degrees of freedom, whereas the current transformation matrix $\mathcal{T}$ of the cell (:yref:`Cell.trsf`) has 9 independent components -- it has 3 additional for arbitrary rotation, which induces no strain. We perform `polar decomposition <http://en.wikipedia.org/wiki/Polar_decomposition#Matrix_polar_decomposition>`_ $\mathcal{T}=\mat{U}\mat{P}$, where $\mat{U}$ is unitary (representing rotation) and $\mat{P}$ is deformation in the proper sense; after converting $\mat{P}$ to strain, it is rotated back by $\mat{U}$ so that it is in global coordinates. To have `true strain <http://en.wikipedia.org/wiki/True_strain#True_strain>`_, the strain matrix should be obtained via
+Global strain tensor $\tens{\eps}$ is symmetric and has 6 degrees of freedom, whereas the current transformation matrix $\mathcal{T}$ of the cell (:yref:`Cell.trsf`) has 9 independent components -- it has 3 additional for arbitrary rotation, which induces no strain. We perform `polar decomposition <http://en.wikipedia.org/wiki/Polar_decomposition#Matrix_polar_decomposition>`_ $\mathcal{T}=\mat{U}\mat{P}$, where $\mat{U}$ is unitary (representing rotation) and $\mat{P}$ is deformation in the proper sense; after converting $\mat{P}$ to strain, it is rotated back by $\mat{U}$ so that it is in global coordinates. To have `true strain <http://en.wikipedia.org/wiki/True_strain#True_strain>`_, the strain matrix is obtained via
.. math:: \mat{\eps}=\mat{U}\log \mat{P}
-but as matrix logarithm is not yet implemented in Eigen, we approximate this (which works only for "small strain" scenarios) as
-
-.. math:: \mat{\eps}=\mat{U}(\mat{P}-\mat{I}_3)
-
-where $\mat{I}_3$ is 3×3 identity matrix.
-
-Once current values of $\tens{\sigma}$, $\tens{K}$ and $\tens{\eps}$ are known, the control algorithm comes. These tensors are converted to 6-vectors using the `Voigt notation <http://en.wikiversity.org/wiki/Introduction_to_Elasticity/Constitutive_relations#Voigt_notation>`_. Since for each component $i\in\{0\dots5\}$ either $\vec{\sig}_i$ of $\vec{\eps}_i$ is prescribed (determined by :yref:`stressMask<Peri3dController.stressMask>`), the equation $\vec{\sig}=\mat{K}\vec{\eps}$ (dealing with large strains, we should use the rate formulation; supposing $\mat{K}$ is quasi-constant in time (linearization), we have $\vec{\dot\sig}=\mat{K}\vec{\dot\eps}$, where all rate terms differ from their counterparts by the same factor $\Delta t$. For that reason, it is only at the end that the rate form is assumed in the proper sense) is decomposed in 2 parts:
-
-* $p$ denotes the part where strain component $\vec{\eps}_i$ is *prescribed*;
-* $u$ denotes the part where strain component $\vec{\eps}_i$ is *unprescribed* (and $\vec{\sig}_i$ is prescribed).
+computed (see `wikipedia <http://en.wikipedia.org/wiki/Logarithm_of_a_matrix#Calculating_the_logarithm_of_a_diagonalizable_matrix>`_ and [Lu1998]_) by spectral decomposition
+
+.. math::
+ :nowrap:
+
+ \begin{align*}
+ \mat{\eps}'&=\mat{V}^{-1}\mat{\eps}\mat{V} \\
+ \log(\mat{\eps})&=\mat{V}\log(\mat{\eps}')\mat{V}^{-1}
+ \end{align*}
+
+where $\mat{\eps}'$ is diagonal and $\log(\mat{\eps}')$ is a diagonal matrix with $\log(\mat{\eps}')_{ii}=\log(\mat{\eps}'_{ii})$. $\mat{V}$ is orthonormal, therefore $\mat{V}^{-1}=\mat{V}^T$.
+
+Once current values of $\tens{\sigma}$, $\tens{K}$ and $\tens{\eps}$ are known, the control algorithm comes. These tensors are converted to 6-vectors using the `Voigt notation <http://en.wikiversity.org/wiki/Introduction_to_Elasticity/Constitutive_relations#Voigt_notation>`_. Since for each component $i\in\{0\dots5\}$ either $\vec{\sig}_i$ of $\vec{\eps}_i$ is prescribed (determined by :yref:`stressMask<Peri3dController.stressMask>`), the equation $\vec{\sig}=\mat{K}\vec{\eps}$; assuming large strains, we take the rate form $\vec{\dot\sig}=\mat{K}\vec{\dot\eps}$ is decomposed in 2 parts:
+
+* $p$ denotes the part where strain component $\vec{\dot\eps}_i$ is *prescribed* (via $\frac{\vec{\eps}^g_i-\vec{\eps}_i}{\Dt}$, where $\vec{\eps}^g_i$ is the corresponding :yref:`goal<Peri3dController.goal>` strain component value);
+* $u$ denotes the part where strain component $\vec{\dot\eps}_i$ is *unprescribed* (and $\vec{\sig}_i$ is prescribed, again via $\frac{\vec{\sig}^g_i-\vec{\sig}_i}{\Dt}$).
The equation then is re-arranged as
.. math::
- \begin{Bmatrix}\vec{\sig}_u \\ \vec{\sig}_p\end{Bmatrix}=\begin{bmatrix}\mat{K}_{uu} & \mat{K}_{up} \\ \mat{K}_{pu} & \mat{K}_{pp}\end{bmatrix}\begin{Bmatrix}\vec{\eps}_u \\ \vec{\eps}_p\end{Bmatrix}
+ \begin{Bmatrix}\vec{\dot\sig}_u \\ \vec{\dot\sig}_p\end{Bmatrix}=\begin{bmatrix}\mat{K}_{uu} & \mat{K}_{up} \\ \mat{K}_{pu} & \mat{K}_{pp}\end{bmatrix}\begin{Bmatrix}\vec{\dot\eps}_u \\ \vec{\dot\eps}_p\end{Bmatrix}
Since we can only apply strain (via the :yref:`velocity gradient<Cell.velGrad>`), we must evaluate
-.. math :: \vec{\eps}_u=\mat{K}_{uu}^{-1}(\vec{\sig}_u-\mat{K}_{up}\vec{\eps}_p)
-
-The goal strain tensor $\tens{\eps}^g$ is the assembled from $\vec{\eps}_u$ and $\vec{\eps}_p$ and is used to prescribe the current :yref:`velocity gradient <Cell.velGrad>`
-
-.. math:: \nabla\vec{v}=\frac{\tens{\eps}^g-\tens{\eps}}{\Delta t}.
-
-If magnitude of any component of $\nabla\vec{v}$ exceeds :yref:`maxStrainRate<Peri3dController.maxStrainRate>`, the whole tensor is scaled accordingly.
+.. math :: \vec{\dot\eps}_u=\mat{K}_{uu}^{-1}(\vec{\dot\sig}_u-\mat{K}_{up}\vec{\dot\eps}_p)
+
+The :yref:`velocity gradient <Cell.velGrad>` is then assembled from $\vec{\eps}_u$ and $\vec{\eps}_p$. If magnitude of any component of $\nabla\vec{v}$ exceeds :yref:`maxStrainRate<Peri3dController.maxStrainRate>`, the whole tensor is scaled accordingly.
'''
=== modified file 'py/pack/pack.py'
--- py/pack/pack.py 2010-06-17 21:13:59 +0000
+++ py/pack/pack.py 2010-06-29 14:32:03 +0000
@@ -266,9 +266,10 @@
rDev*=scale; X*=scale; Y*=scale; Z*=scale
memoDbgMsg("Considering packing (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(R,.5*rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp))))
if not isPeri and wantPeri: memoDbgMsg("REJECT: is not periodic, which is requested."); continue
+ if wantPeri and (X/x1>0.9 or X/x1<0.6): memoDbgMsg("REJECT: initSize differs too much from scaled packing size."); continue
if (rRelFuzz==0 and rDev!=0) or (rRelFuzz!=0 and rDev==0) or (rRelFuzz!=0 and abs((rDev-rRelFuzz)/rRelFuzz)>1e-2): memoDbgMsg("REJECT: radius fuzz differs too much (%g, %g desired)"%(rDev,rRelFuzz)); continue # radius fuzz differs too much
if isPeri and wantPeri:
- if spheresInCell>NN: memoDbgMsg("REJECT: Number of spheres in the packing too small"); continue
+ if spheresInCell>NN and spheresInCell>0: memoDbgMsg("REJECT: Number of spheres in the packing too small"); continue
if abs((y1/x1)/(Y/X)-1)>0.3 or abs((z1/x1)/(Z/X)-1)>0.3: memoDbgMsg("REJECT: proportions (y/x=%g, z/x=%g) differ too much from what is desired (%g, %g)."%(Y/X,Z/X,y1/x1,z1/x1)); continue
else:
if (X<fullDim[0] or Y<fullDim[1] or Z<fullDim[2]): memoDbgMsg("REJECT: not large enough"); continue # not large enough
@@ -420,7 +421,7 @@
O.run(); O.wait()
ret=SpherePack()
ret.fromSimulation()
- _memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri=True,fullDim=Vector3(0,0,0)) # fullDim unused
+ _memoizePacking(memoizeDb,ret,radius,rRelFuzz,wantPeri=True,fullDim=Vector3(0,0,0)) # fullDim unused
O.switchScene()
return ret
=== modified file 'py/plot.py'
--- py/plot.py 2010-05-06 10:11:20 +0000
+++ py/plot.py 2010-06-29 14:32:03 +0000
@@ -4,7 +4,21 @@
Module containing utility functions for plotting inside yade. See :ysrc:`scripts/simple-scene-plot.py` or :ysrc:`examples/concrete/uniax.py` for example of usage.
"""
-import matplotlib
+import matplotlib,os
+# running in batch
+#
+# If GtkAgg is the default, X must be working, which is not the case
+# with batches (DISPLAY is unset in such case) and importing pylab fails then.
+#
+# Agg does not require the GUI part and works withou any DISPLAY active
+# just fine.
+#
+# see http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg04320.html
+# and https://lists.launchpad.net/yade-users/msg03289.html
+#
+if os.environ.has_key('PARAM_TABLE'):
+ matplotlib.use('Agg')
+
#matplotlib.use('TkAgg')
#matplotlib.use('GTKCairo')
#matplotlib.use('QtAgg')
=== modified file 'scripts/test/peri3d.py'
--- scripts/test/peri3d.py 2010-05-26 17:49:15 +0000
+++ scripts/test/peri3d.py 2010-06-29 14:32:03 +0000
@@ -1,4 +1,5 @@
from yade import pack,log
+O.materials.append(CpmMat(neverDamage=False,young=1e7,frictionAngle=.5,G_over_E=.2,sigmaT=3e6,epsCrackOnset=1e-2,relDuctility=10))
sp=pack.randomPeriPack(.5,.5,Vector3(10,10,10),memoizeDb='/tmp/packDb.sqlite')
O.periodic=True
O.cell.refSize=sp.cellSize
@@ -7,23 +8,33 @@
log.setLevel('Peri3dController',log.TRACE)
O.engines=[
ForceResetter(),
- BoundDispatcher([Bo1_Sphere_Aabb()]),
+ BoundDispatcher([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='bo1s')]),
InsertionSortCollider(),
InteractionDispatchers(
- [Ig2_Sphere_Sphere_Dem3DofGeom()],
- [Ip2_FrictMat_FrictMat_FrictPhys()],
- [Law2_Dem3DofGeom_FrictPhys_Basic()]
+ [Ig2_Sphere_Sphere_Dem3DofGeom(distFactor=1.5,label='ig2ss')],
+ #[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_Basic()]
+ [Ip2_CpmMat_CpmMat_CpmPhys()],[Law2_Dem3DofGeom_CpmPhys_Cpm()]
+
),
- NewtonIntegrator(homotheticCellResize=1,damping=.4),
- #Peri3dController(goal=Matrix3(-.1,0,0, 0,-.1,0, 0,0,-.1),stressMask=0,maxStrainRate=1e-1,label='ctrl'),
- #Peri3dController(goal=Matrix3(0,.2,0, 0,0,0, 0,0,0),stressMask=0,maxStrainRate=1,label='ctrl'),
- Peri3dController(goal=Matrix3(-1e3,.1,0, 0,.1,0, 0,0,0),stressMask=0b0001,maxStrainRate=1,label='ctrl'),
+ NewtonIntegrator(homotheticCellResize=1,damping=.8),
+ #Peri3dController(goal=Matrix3(-.1,0,0, 0,-.1,0, 0,0,-.1),stressMask=0,maxStrainRate=1e-1,label='p3d'),
+ #Peri3dController(goal=Matrix3(0,.2,0, 0,0,0, 0,0,0),stressMask=0,maxStrainRate=1,label='p3d'),
+ Peri3dController(goal=Matrix3(.2,0,0, 0,0,0, 0,0,0),stressMask=0b000110,maxStrainRate=.1,label='p3d'),
PeriodicPythonRunner(iterPeriod=10,command='addData()')
]
O.step()
+bo1s.aabbEnlargeFactor=ig2ss.distFactor=-1
+zRot=-pi/4.
+O.cell.trsf=Matrix3(cos(zRot),-sin(zRot),0,sin(zRot),cos(zRot),0,0,0,1)
+O.saveTmp()
from yade import plot
def addData():
- plot.addData(sxx=ctrl.stress[0,0],syy=ctrl.stress[1,1],szz=ctrl.stress[2,2],exx=ctrl.strain[0,0],eyy=ctrl.strain[1,1],ezz=ctrl.strain[2,2],t=O.time)
+ plot.addData(sxx=p3d.stress[0,0],syy=p3d.stress[1,1],szz=p3d.stress[2,2],exx=p3d.strain[0,0],eyy=p3d.strain[1,1],ezz=p3d.strain[2,2],t=O.time)
-plot.plots={'exx':('sxx'),}
+plot.plots={'t':('sxx',None,'eyy','ezz')}#,'exx':('sxx'),}
+from yade import qt
+rrr=qt.Renderer()
+rrr.intrPhys,rrr.shape=True,False
+Gl1_CpmPhys.dmgLabel,Gl1_CpmPhys.colorStrainRatio=False,10
+qt.View()
=== added file 'scripts/test/peri8.py'
--- scripts/test/peri8.py 1970-01-01 00:00:00 +0000
+++ scripts/test/peri8.py 2010-06-29 14:32:03 +0000
@@ -0,0 +1,9 @@
+a=1.; r=a/4.
+O.bodies.append(
+ [utils.sphere(c,r) for c in [(r,r,r),(3*r,r,r),(3*r,3*r,r),(r,3*r,r),(r,r,3*r),(3*r,r,3*r),(3*r,3*r,3*r),(r,3*r,3*r)]]
+)
+O.periodic=True
+O.cell.refSize=(a,a,a)
+zRot=-pi/4.
+O.cell.trsf=Matrix3(cos(zRot),-sin(zRot),0,sin(zRot),cos(zRot),0,0,0,1)
+p7=O.bodies[7].state.pos
=== modified file 'yadeSCons.py'
--- yadeSCons.py 2010-04-13 21:09:57 +0000
+++ yadeSCons.py 2010-06-29 14:32:03 +0000
@@ -98,7 +98,7 @@
if cacheFile: plugInfo.close()
return pp
-def getWantedPlugins(plugInfo,excludes,features,linkStrategy):
+def getWantedPlugins(plugInfo,excludes,features,linkStrategy,hotPlugins):
"""Use pluginInfo (generated by scanAllPlugins) and return only plugins that we should build,
based on excludes and available features.
@@ -112,15 +112,17 @@
if not plug.feats<=feats:
continue # plugin needs more feature than we have
ret[plug.name]=plug
- for p in plugInfo.values(): p.obj=getPluginObj(p,linkStrategy)
+ for p in plugInfo.values(): p.obj=getPluginObj(p,linkStrategy,hotPlugins)
for p in plugInfo.values(): p.libs=getPluginLibs(p,plugInfo)
return ret
-def getPluginObj(plug,linkStrategy):
+def getPluginObj(plug,linkStrategy,hotPlugins):
"""Return name of library this plugin will be compiled into, based on current linkStrategy."""
if linkStrategy=='per-class': return plug.name
elif linkStrategy=='per-pkg': return plug.module
- elif linkStrategy=='monolithic': return 'plugins'
+ elif linkStrategy=='monolithic':
+ if plug.name in hotPlugins: return plug.name
+ else: return 'plugins'
elif linkStrategy=='static': return 'plugins'
def getPluginLibs(p,plugInfo):
Follow ups
