← Back to team overview

yade-dev team mailing list archive

Re: Bugs in Cell

 

Hi,

The bug is in Cell.hpp in wrapPt method. If the cell in not rotated, it is ok. But if you apply rotation, periodic shifts are performed along new (local) coordinates, however, the cell is not periodic in this direction (sorry, Vaclav, I yesterday said tho opposite :-)

the solution would be to rotate all particles together wih the cell.

Best wishes,

Jan

______________________________________________________________
> Od: noreply@xxxxxxxxxxxxx
> Komu: Yade developers <yade-dev@xxxxxxxxxxxxxxxxxxx>
> Datum: 29.06.2010 16:33
> Předmět: [Yade-dev] [Branch ~yade-dev/yade/trunk] Rev 2300: 1. Bugs in Cell(!!! no yetfixed: 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:: %snt:members:nt:undoc-members:nn"%(klass) # t:inherited-members:n # nt:show-inheritance:
> def autodocDerived(klass,doSections=True):
>+	global docClasses
> 	ret=''
> 	if doSections: ret+='%sn'%klass+'^'*100+'nn'
> 	#ret+='.. inheritance-diagram:: %snn'%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 'tt"%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 'tt"%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::nntdigraph %s {nttrankdir=RL;nttmargin=.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+='tt"%s" -> "%s" [arrowsize=0.5,style="setlinewidth(0.5)"];n'%(c,base)
> 		except NameError:
> 			pass
>@@ -91,11 +99,15 @@
> 	return ret+'t}nn'
> 
> 
>-def sect(title,text,tops):
>-	return title+'n'+100*'-'+'nn'+text+'nn'+'nn'.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*'-'+'nn'+text+'nn'+'nn'.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.nnInteractionGeometry will be computed when interactionDetectionFactor*(rad1+rad2) > distance.nn.. note::nt This parameter is functionally coupled with :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor`, which will create larger bounding boxes and should be of the same value.nn.. warning::ntFunctionally 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 valuen"<<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 stressn"<<stress);
>+	LOG_TRACE("Updated stress (from "<<nIntr<<" interactions)n"<<stress);
> 	LOG_TRACE("Updated stiffness tensorn"<<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 sanitizationn"<<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}^2left[K_N n_i n_j n_k n_l+K_Tleft(frac{n_j n_k delta_{il}+n_j n_ldelta_{ik}+n_i n_kdelta_{jl}+n_i n_ldelta_{jk}}{4}-n_i n_j n_kn_lright)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 $iin{0dots5}$ 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{dotsig}=mat{K}vec{doteps}$, 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 $iin{0dots5}$ 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{dotsig}=mat{K}vec{doteps}$ is decomposed in 2 parts:
>+
>+* $p$ denotes the part where strain component $vec{doteps}_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{doteps}_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}_pend{Bmatrix}=begin{bmatrix}mat{K}_{uu} & mat{K}_{up} \ mat{K}_{pu} & mat{K}_{pp}end{bmatrix}begin{Bmatrix}vec{eps}_u \ vec{eps}_pend{Bmatrix}
>+	begin{Bmatrix}vec{dotsig}_u \ vec{dotsig}_pend{Bmatrix}=begin{bmatrix}mat{K}_{uu} & mat{K}_{up} \ mat{K}_{pu} & mat{K}_{pp}end{bmatrix}begin{Bmatrix}vec{doteps}_u \ vec{doteps}_pend{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:: nablavec{v}=frac{tens{eps}^g-tens{eps}}{Delta t}.
>-
>-If magnitude of any component of $nablavec{v}$ exceeds :yref:`maxStrainRate<Peri3dController.maxStrainRate>`, the whole tensor is scaled accordingly.
>+.. math :: vec{doteps}_u=mat{K}_{uu}^{-1}(vec{dotsig}_u-mat{K}_{up}vec{doteps}_p)
>+
>+The :yref:`velocity gradient <Cell.velGrad>` is then assembled from $vec{eps}_u$ and $vec{eps}_p$. If magnitude of any component of $nablavec{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):
>
>
>_______________________________________________
>Mailing list: https://launchpad.net/~yade-dev
>Post to     : yade-dev@xxxxxxxxxxxxxxxxxxx
>Unsubscribe : https://launchpad.net/~yade-dev
>More help   : https://help.launchpad.net/ListHelp
>
>



Follow ups

References