← Back to team overview

yade-dev team mailing list archive

  1. yade-dev team
  2. Mailing list archive
  3. Message #06636

[Branch ~yade-dev/yade/trunk] Rev 2627: 1. Fix an error in Law2_L3Geom_FrictPhys_ElPerfPl (forces were applied at the step the interactio...

  Thread Previous • Date Previous • • Thread Next 
------------------------------------------------------------
revno: 2627
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Thu 2010-12-23 15:01:33 +0100
message:
  1. Fix an error in Law2_L3Geom_FrictPhys_ElPerfPl (forces were applied at the step the interaction broke); results are now identical to ScGeom/CundallStrack (in the normal sense at least)
  2. Add the possibility to specify y-data as functions returning y-data names themselves (see scripts/tests/energy.py)
  3. Merge 2 consecutive loops in NewtonIntegrator for clumps
  4. Rename plot.legendPosition, plot.legendPositionSecondary to plot.legendLoc (tuple of 2 values), to make it consistent with matplotlib terminology (pylab.legend(...,loc=...))
modified:
  pkg/dem/Clump.cpp
  pkg/dem/Clump.hpp
  pkg/dem/L3Geom.cpp
  pkg/dem/NewtonIntegrator.cpp
  py/_utils.cpp
  py/plot.py
  py/utils.py
  scripts/test/energy.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 'pkg/dem/Clump.cpp'
--- pkg/dem/Clump.cpp	2010-12-01 11:09:53 +0000
+++ pkg/dem/Clump.cpp	2010-12-23 14:01:33 +0000
@@ -41,25 +41,25 @@
 void Clump::moveMembers(const shared_ptr<Body>& clumpBody, Scene* scene){
 	const shared_ptr<Clump>& clump=YADE_PTR_CAST<Clump>(clumpBody->shape);
 	const shared_ptr<State>& clumpState=clumpBody->state;
-
 	FOREACH(MemberMap::value_type& B, clump->members){
-		// B->first is Body::id_t, B->second is local Se3r of that body in the clump
-		State* subState=Body::byId(B.first,scene)->state.get();
-		Vector3r& subPos(B.second.position); Quaternionr& subOri(B.second.orientation);
-		//LOG_TRACE("Old #"<<I->first<<"position: "<<subRBP->se3.position);
-
-		// position update
-		subState->pos=clumpState->pos+clumpState->ori*subPos;
-		subState->ori=clumpState->ori*subOri;
-		// velocity update
-		subState->vel=clumpState->vel+clumpState->angVel.cross(subState->pos-clumpState->pos);
-		subState->angVel=clumpState->angVel;
-
-		//LOG_TRACE("New #"<<I->first<<"position: "<<subRBP->se3.position);
-		//LOG_TRACE("Clump #"<<getId()<<" moved #"<<I->first<<".");
+		// B.first is Body::id_t, B.second is local Se3r of that body in the clump
+		Clump::moveMember(clumpState, Body::byId(B.first,scene)->state, B.second);
 	}
 }
 
+void Clump::moveMember(const shared_ptr<State>& clumpState, const shared_ptr<State>& subState, const Se3r& relSe3){
+	const Vector3r& subPos(relSe3.position); const Quaternionr& subOri(relSe3.orientation);
+	//LOG_TRACE("Old #"<<I->first<<"position: "<<subRBP->se3.position);
+	// position update
+	subState->pos=clumpState->pos+clumpState->ori*subPos;
+	subState->ori=clumpState->ori*subOri;
+	// velocity update
+	subState->vel=clumpState->vel+clumpState->angVel.cross(subState->pos-clumpState->pos);
+	subState->angVel=clumpState->angVel;
+	//LOG_TRACE("New #"<<I->first<<"position: "<<subRBP->se3.position);
+	//LOG_TRACE("Clump #"<<getId()<<" moved #"<<I->first<<".");
+}
+
 /*! Clump's se3 will be updated (origin at centroid and axes coincident with principal inertia axes) and subSe3 modified in such a way that members positions in world coordinates will not change.
 
 	TODO: numerical integration of inertia based on regular space sampling with relative tolerance WRT minimum sphere radius

=== modified file 'pkg/dem/Clump.hpp'
--- pkg/dem/Clump.hpp	2010-12-01 11:09:53 +0000
+++ pkg/dem/Clump.hpp	2010-12-23 14:01:33 +0000
@@ -60,7 +60,9 @@
 		static void del(const shared_ptr<Body>& clump, const shared_ptr<Body>& subBody);
 		//! Recalculate physical properties of Clump.
 		static void updateProperties(const shared_ptr<Body>& clump, bool intersecting);
-		//! Calculate positions and orientations of members based on my own Se3.
+		//! Calculate positions and orientations of members based on relative Se3
+		static void moveMember(const shared_ptr<State>& clumpState, const shared_ptr<State>& subState, const Se3r& relSe3);
+		//! Loop around moveMember
 		static void moveMembers(const shared_ptr<Body>& clump, Scene* scene);
 		//! update member positions after clump being moved by mouse (in case simulation is paused and engines will not do that).
 		void userForcedDisplacementRedrawHook(){ throw runtime_error("Clump::userForcedDisplacementRedrawHook not yet implemented (with Clump as subclass of Shape).");}

=== modified file 'pkg/dem/L3Geom.cpp'
--- pkg/dem/L3Geom.cpp	2010-12-21 15:31:40 +0000
+++ pkg/dem/L3Geom.cpp	2010-12-23 14:01:33 +0000
@@ -179,7 +179,7 @@
 	// wall interacting from both sides: normal depends on sphere's position
 	assert(sense==-1 || sense==0 || sense==1);
 	if(sense==0) normal[ax]=dist>0?1.:-1.;
-	else normal[ax]=sense==1?1.:-1;
+	else normal[ax]=(sense==1?1.:-1);
 	Real uN=normal[ax]*dist-radius; // takes in account sense, radius and distance
 
 	if(!I->geom){
@@ -196,7 +196,7 @@
 		return true;
 	}
 	L3Geom& g=I->geom->cast<L3Geom>();
-	if(g.normal!=normal) throw std::logic_error(("Ig2_Wall_Sphere_L3Geom_Inc: normal changed in Wall+Sphere ##"+lexical_cast<string>(I->getId1())+"+"+lexical_cast<string>(I->getId2())+" (Wall.sense≠0?)").c_str());
+	if(g.normal!=normal) throw std::logic_error(("Ig2_Wall_Sphere_L3Geom_Inc: normal changed in Wall+Sphere ##"+lexical_cast<string>(I->getId1())+"+"+lexical_cast<string>(I->getId2())+" (with Wall.sense=0, a particle might cross the Wall plane, if Δt is too high)").c_str());
 	Vector3r normTwistVec=normal*scene->dt*.5*normal.dot(state1.angVel+state2.angVel);
 
 	/* TODO: below is a lot of code common with Sphere+Sphere; identify, factor out in a separate func */
@@ -230,7 +230,7 @@
 	Vector3r& localF(geom->F);
 	localF=geom->relU().cwise()*Vector3r(phys->kn,phys->ks,phys->ks);
 	// break if necessary
-	if(localF[0]>0 && !noBreak){ scene->interactions->requestErase(I->getId1(),I->getId2()); }
+	if(localF[0]>0 && !noBreak){ scene->interactions->requestErase(I->getId1(),I->getId2()); return; }
 
 	if(!noSlip){
 		// plastic slip, if necessary; non-zero elastic limit only for compression

=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp	2010-12-22 10:55:23 +0000
+++ pkg/dem/NewtonIntegrator.cpp	2010-12-23 14:01:33 +0000
@@ -170,10 +170,10 @@
 				// rotation
 				if(aspherical) leapfrogAsphericalRotate(state,id,dt,M);
 				else leapfrogSphericalRotate(state,id,dt);
-				Clump::moveMembers(b,scene);
-				// save max velocity for all members of the clump
-				FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
+				// move individual members of the clump, save maxima velocity (for collider stride)
+				FOREACH(Clump::MemberMap::value_type mm, b->shape->cast<Clump>().members){
 					const shared_ptr<Body>& member=Body::byId(mm.first,scene); assert(member->isClumpMember());
+					Clump::moveMember(b->state,member->state,mm.second);
 					saveMaximaVelocity(mm.first,member->state.get());
 				}
 			}

=== modified file 'py/_utils.cpp'
--- py/_utils.cpp	2010-12-21 12:19:50 +0000
+++ py/_utils.cpp	2010-12-23 14:01:33 +0000
@@ -436,7 +436,7 @@
 
 	py::def("PWaveTimeStep",PWaveTimeStep,"Get timestep accoring to the velocity of P-Wave propagation; computed from sphere radii, rigidities and masses.");
 	py::def("getSpheresVolume",Shop__getSpheresVolume,"Compute the total volume of spheres in the simulation (might crash for now if dynamic bodies are not spheres)");
-	py::def("porosity",Shop__getPorosity,(py::arg("volume")=-1),"Compute packing porosity $\\frac{V-V_s}{V}$ where $V$ is overall volume and $V_s$ is volume of spheres.\\:param float volume: overall volume which must be specified for aperiodic simulations. For periodic simulations, current volume of the :yref:`Cell` is used.");
+	py::def("porosity",Shop__getPorosity,(py::arg("volume")=-1),"Compute packing porosity $\\frac{V-V_s}{V}$ where $V$ is overall volume and $V_s$ is volume of spheres.\n\n:param float volume: overall volume which must be specified for aperiodic simulations. For periodic simulations, current volume of the :yref:`Cell` is used.\n");
 	py::def("aabbExtrema",aabbExtrema,(py::arg("cutoff")=0.0,py::arg("centers")=false),"Return coordinates of box enclosing all bodies\n\n:param bool centers: do not take sphere radii in account, only their centroids\n:param float∈〈0…1〉 cutoff: relative dimension by which the box will be cut away at its boundaries.\n\n\n:return: (lower corner, upper corner) as (Vector3,Vector3)\n\n");
 	py::def("ptInAABB",isInBB,"Return True/False whether the point p is within box given by its min and max corners");
 	py::def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(py::args("axis","distFactor"),"Return list of ids for spheres (only) that are on extremal ends of the specimen along given axis; distFactor multiplies their radius so that sphere that do not touch the boundary coordinate can also be returned."));
@@ -451,7 +451,7 @@
 	py::def("elasticEnergy",elasticEnergyInAABB);
 	py::def("inscribedCircleCenter",inscribedCircleCenter,(py::arg("v1"),py::arg("v2"),py::arg("v3")),"Return center of inscribed circle for triangle given by its vertices *v1*, *v2*, *v3*.");
 	py::def("unbalancedForce",&Shop__unbalancedForce,(py::args("useMaxForce")=false),"Compute the ratio of mean (or maximum, if *useMaxForce*) summary force on bodies and maximum force magnitude on interactions. For perfectly static equilibrium, summary force on all bodies is zero (since forces from interactions cancel out and induce no acceleration of particles); this ratio will tend to zero as simulation stabilizes, though zero is never reached because of finite precision computation. Sufficiently small value can be e.g. 1e-2 or smaller, depending on how much equilibrium it should be.");
-	py::def("kineticEnergy",Shop__kineticEnergy,(py::args("findMaxId")=false),"Compute overall kinetic energy of the simulation as\n\n.. math:: \\sum\\frac{1}{2}\\left(m_i\\vec{v}_i^2+\\vec{\\omega}(\\mat{I}\\vec{\\omega}^T)\\right).\n\n.. warning::\n\n\tNo transformation of inertia tensor (in local frame) $\\mat{I}$ is done, although it is multiplied by angular velocity $\\vec{\\omega}$ (in global frame); the value will not be accurate for aspherical particles.\n");
+	py::def("kineticEnergy",Shop__kineticEnergy,(py::args("findMaxId")=false),"Compute overall kinetic energy of the simulation as\n\n.. math:: \\sum\\frac{1}{2}\\left(m_i\\vec{v}_i^2+\\vec{\\omega}(\\mat{I}\\vec{\\omega}^T)\\right).\n\nFor :yref:`aspherical<Body.aspherical>` bodies, the inertia tensor $\\mat{I}$ is transformed to global frame, before multiplied by $\\vec{\\omega}$, therefore the value should be accurate.\n");
 	py::def("sumForces",sumForces,(py::arg("ids"),py::arg("direction")),"Return summary force on bodies with given *ids*, projected on the *direction* vector.");
 	py::def("sumTorques",sumTorques,(py::arg("ids"),py::arg("axis"),py::arg("axisPt")),"Sum forces and torques on bodies given in *ids* with respect to axis specified by a point *axisPt* and its direction *axis*.");
 	py::def("sumFacetNormalForces",sumFacetNormalForces,(py::arg("ids"),py::arg("axis")=-1),"Sum force magnitudes on given bodies (must have :yref:`shape<Body.shape>` of the :yref:`Facet` type), considering only part of forces perpendicular to each :yref:`facet's<Facet>` face; if *axis* has positive value, then the specified axis (0=x, 1=y, 2=z) will be used instead of facet's normals.");
@@ -467,7 +467,7 @@
 	py::def("wireNone",wireNone,"Set :yref:`Shape::wire` on all bodies to False, rendering them as solids.");
 	py::def("wireNoSpheres",wireNoSpheres,"Set :yref:`Shape::wire` to True on non-spherical bodies (:yref:`Facets<Facet>`, :yref:`Walls<Wall>`).");
 	py::def("flipCell",&Shop::flipCell,(py::arg("flip")=Matrix3r(Matrix3r::Zero())),"Flip periodic cell so that angles between $R^3$ axes and transformed axes are as small as possible. This function relies on the fact that periodic cell defines by repetition or its corners regular grid of points in $R^3$; however, all cells generating identical grid are equivalent and can be flipped one over another. This necessiatates adjustment of :yref:`Interaction.cellDist` for interactions that cross boundary and didn't before (or vice versa), and re-initialization of collider. The *flip* argument can be used to specify desired flip: integers, each column for one axis; if zero matrix, best fit (minimizing the angles) is computed automatically.\n\nIn c++, this function is accessible as ``Shop::flipCell``.\n\n.. warning::\n\t This function is currently broken and should not be used.");
-	py::def("getViscoelasticFromSpheresInteraction",getViscoelasticFromSpheresInteraction,(py::arg("tc"),py::arg("en"),py::arg("es")),"Get viscoelastic interaction parameters from analytical solution of a pair spheres collision problem. \n\n:Parameters:\n\t`m` : float\n\t\tsphere mass\n\t`tc` : float\n\t\tcollision time\n\t`en` : float\n\t\tnormal restitution coefficient\n\t`es` : float\n\t\ttangential restitution coefficient.\n\n\n:return: \n\tdict with keys:\n\n\tkn : float\n\t\tnormal elastic coefficient computed as:\n\n.. math:: k_n=\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_n)^2\\right)\n\n\tcn : float\n\t\tnormal viscous coefficient computed as:\n\n.. math:: c_n=-\\frac{2m}{t_c}\\ln e_n\n\n\n\tkt : float\n\t\ttangential elastic coefficient computed as:\n\n.. math:: k_t=\\frac27\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_t)^2\\right)\n\n\tct : float\n\t\ttangential viscous coefficient computed as:\n\n.. math:: c_t=-\\frac27\\frac{m}{t_c}\\ln e_t.\n\n\nFor details see [Pournin2001]_. ");
+	py::def("getViscoelasticFromSpheresInteraction",getViscoelasticFromSpheresInteraction,(py::arg("tc"),py::arg("en"),py::arg("es")),"Compute viscoelastic interaction parameters from analytical solution of a pair spheres collision problem:\n\n\n.. math::\n\t:nowrap:\n\n\n\t\\begin{align*}k_n&=\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_n)^2\\right)\\\\  c_n&=-\\frac{2m}{t_c}\\ln e_n \\\\k_t&=\\frac27\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_t)^2\\right)  \\\\ c_t=-\\frac27\\frac{m}{t_c}\\ln e_t \\end{align*}\n\n\nwhere $k_n$, $c_n$ are normal elastic and viscous coefficients and $k_t$, $c_t$ shear elastic and viscous coefficients. For details see [Pournin2001]_.\n\n:param float m: sphere mass $m$\n:param float tc: collision time $t_c$\n:param float en: normal restitution coefficient $e_n$\n:param float es: tangential restitution coefficient $e_s$\n:return: dictionary with keys ``kn`` (the value of $k_n$), ``cn`` ($c_n$), ``kt`` ($k_t$), ``ct`` ($c_t$).");
 	py::def("stressTensorOfPeriodicCell",Shop__stressTensorOfPeriodicCell,(py::args("smallStrains")=false),"Compute overall (macroscopic) stress of periodic cell using equation published in [Kuhl2001]_:\n\n.. math:: \\vec{\\sigma}=\\frac{1}{V}\\sum_cl^c[\\vec{N}^cf_N^c+\\vec{T}^{cT}\\cdot\\vec{f}^c_T],\n\nwhere $V$ is volume of the cell, $l^c$ length of interaction $c$, $f^c_N$ normal force and $\\vec{f}^c_T$ shear force. Sumed are values over all interactions $c$. $\\vec{N}^c$ and $\\vec{T}^{cT}$ are projection tensors (see the original publication for more details):\n\n.. math:: \\vec{N}=\\vec{n}\\otimes\\vec{n}\\rightarrow N_{ij}=n_in_j\n\n.. math:: \\vec{T}^T=\\vec{I}_{sym}\\cdot\\vec{n}-\\vec{n}\\otimes\\vec{n}\\otimes\\vec{n}\\rightarrow T^T_{ijk}=\\frac{1}{2}(\\delta_{ik}\\delta_{jl}+\\delta_{il}\\delta_{jk})n_l-n_in_jn_k\n\n.. math:: \\vec{T}^T\\cdot\\vec{f}_T\\equiv T^T_{ijk}f_k=(\\delta_{ik}n_j/2+\\delta_{jk}n_i/2-n_in_jn_k)f_k=n_jf_i/2+n_if_j/2-n_in_jn_kf_k,\n\nwhere $n$ is unit vector oriented along the interaction (:yref:`normal<GenericSpheresContact::normal>`) and $\\delta$ is Kronecker's delta. As $\\vec{n}$ and $\\vec{f}_T$ are perpendicular (therfore $n_if_i=0$) we can write\n\n.. math:: \\sigma_{ij}=\\frac{1}{V}\\sum l[n_in_jf_N+n_jf^T_i/2+n_if^T_j/2]\n\n:param bool smallStrains: if false (large strains), real values of volume and interaction lengths are computed. If true, only :yref:`refLength<Dem3DofGeom::refLength>` of interactions and initial volume are computed (can save some time).\n\n:return: macroscopic stress tensor as Matrix3");
 	py::def("normalShearStressTensors",Shop__normalShearStressTensors,(py::args("compressionPositive")=false),"Compute overall stress tensor of the periodic cell decomposed in 2 parts, one contributed by normal forces, the other by shear forces. The formulation can be found in [Thornton2000]_, eq. (3):\n\n.. math:: \\tens{\\sigma}_{ij}=\\frac{2}{V}\\sum R N \\vec{n}_i \\vec{n}_j+\\frac{2}{V}\\sum R T \\vec{n}_i\\vec{t}_j\n\nwhere $V$ is the cell volume, $R$ is \"contact radius\" (in our implementation, current distance between particle centroids), $\\vec{n}$ is the normal vector, $\\vec{t}$ is a vector perpendicular to $\\vec{n}$, $N$ and $T$ are norms of normal and shear forces.");
 	py::def("maxOverlapRatio",maxOverlapRatio,"Return maximum overlap ration in interactions (with :yref:`ScGeom`) of two :yref:`spheres<Sphere>`. The ratio is computed as $\\frac{u_N}{2(r_1 r_2)/r_1+r_2}$, where $u_N$ is the current overlap distance and $r_1$, $r_2$ are radii of the two spheres in contact.");

=== modified file 'py/plot.py'
--- py/plot.py	2010-12-21 20:37:08 +0000
+++ py/plot.py	2010-12-23 14:01:33 +0000
@@ -12,7 +12,7 @@
 # safe to import even if Tk will not be used
 import mtTkinter as Tkinter
 
-import matplotlib,os,time,math
+import matplotlib,os,time,math,itertools
 
 # running in batch
 #
@@ -43,10 +43,8 @@
 xylabels={}
 "Dictionary of 2-tuples specifying (xlabel,ylabel) for respective plots; if either of them is None, the default auto-generated title is used."
 
-legendPosition='upper left'
-"Placement of the legend on the plot"
-legendPositionSecondary='upper right'
-"Placement of the scondary legend on the plot"
+legendLoc=('upper left','upper right')
+"Location of the y1 and y2 legends on the plot, if y2 is active."
 
 live=True if yade.runtime.hasDisplay else False
 "Enable/disable live plot updating. Disabled by default for now, since it has a few rough edges."
@@ -142,13 +140,13 @@
 		except IndexError: x,y=0,0
 		# this could be written in a nicer way, very likely
 		pt=numpy.ndarray((2,),buffer=numpy.array([x,y]))
-		self.scatter.set_offsets(pt)
+		if self.scatter: self.scatter.set_offsets(pt)
 
 currLineRefs=[]
 liveTimeStamp=0 # timestamp when live update was started, so that the old thread knows to stop if that changes
 nan=float('nan')
 
-def createPlots(subPlots=False):
+def createPlots(subPlots=True):
 	global currLineRefs
 	figs=set([l.line.get_axes().get_figure() for l in currLineRefs]) # get all current figures
 	for f in figs: pylab.close(f) # close those
@@ -171,27 +169,55 @@
 				y1=False; continue
 			if y1: plots_p_y1.append(d)
 			else: plots_p_y2.append(d)
-			if d[0] not in data.keys(): missing.add(d[0])
+			if d[0] not in data.keys() and not callable(d[0]): missing.add(d[0])
 		if missing:
 			if len(data.keys())==0 or len(data[data.keys()[0]])==0: # no data at all yet, do not add garbage NaNs
 				for m in missing: data[m]=[]
 			else:
 				print 'Missing columns in plot.data, adding NaN: ',','.join(list(missing))
 				addDataColumns(missing)
-		# create y1 lines
-		for d in plots_p_y1:
-			line,=pylab.plot(data[pStrip],data[d[0]],d[1])
-			# use (0,0) if there are no data yet
-			scatterPt=[0,0] if len(data[pStrip])==0 else (data[pStrip][current],data[d[0]][current])
-			# if current value is NaN, use zero instead
-			scatter=pylab.scatter(scatterPt[0] if not math.isnan(scatterPt[0]) else 0,scatterPt[1] if not math.isnan(scatterPt[1]) else 0,color=line.get_color())
-			currLineRefs.append(LineRef(line,scatter,data[pStrip],data[d[0]]))
-		# create the legend
-		l=pylab.legend([xlateLabel(_p[0]) for _p in plots_p_y1],loc=(legendPosition if len(plots_p_y2)>0 else 'best'))
-		if hasattr(l,'draggable'): l.draggable(True)
-		pylab.ylabel((', '.join([xlateLabel(_p[0]) for _p in plots_p_y1])) if p not in xylabels or not xylabels[p][1] else xylabels[p][1])
-		pylab.xlabel(xlateLabel(pStrip) if (p not in xylabels or not xylabels[p][0]) else xylabels[p][0])
-		if scientific: pylab.ticklabel_format(style='sci',scilimits=(0,0),axis='both')
+		def createLines(pStrip,ySpecs,isY1=True,y2Exists=False):
+			'''Create data lines from specifications; this code is common for y1 and y2 axes;
+			it handles y-data specified as callables, which might create additional lines when updated with liveUpdate.
+			'''
+			# save the original specifications; they will be smuggled into the axes object
+			# the live updated will run yNameFuncs to see if there are new lines to be added
+			# and will add them if necessary
+			yNameFuncs=set([d[0] for d in ySpecs if callable(d[0])])
+			yNames=set()
+			ySpecs2=[]
+			for ys in ySpecs:
+				if not callable(ys[0]): ySpecs2.append(ys)
+				# ys[0]() must return list of strings, which are added to ySpecs2; line specifier is synthesized by tuplifyYAxis and cannot be specified by the user
+				else: ySpecs2+=[(ret,ys[1]) for ret in ys[0]()]
+			if len(ySpecs2)==0:
+				print 'yade.plot: creating fake plot, since there are no y-data yet'
+				line,=pylab.plot([nan],[nan])
+				currLineRefs.append(LineRef(line,None,[nan],[nan]))
+			for d in ySpecs2:
+				yNames.add(d)
+				line,=pylab.plot(data[pStrip],data[d[0]],d[1],label=xlateLabel(d[0]))
+				# use (0,0) if there are no data yet
+				scatterPt=[0,0] if len(data[pStrip])==0 else (data[pStrip][current],data[d[0]][current])
+				# if current value is NaN, use zero instead
+				scatter=pylab.scatter(scatterPt[0] if not math.isnan(scatterPt[0]) else 0,scatterPt[1] if not math.isnan(scatterPt[1]) else 0,color=line.get_color())
+				currLineRefs.append(LineRef(line,scatter,data[pStrip],data[d[0]]))
+			labelLoc=(legendLoc[0 if isY1 else 1] if y2Exists>0 else 'best')
+			l=pylab.legend(loc=labelLoc)
+			if hasattr(l,'draggable'): l.draggable(True)
+			if isY1:
+				pylab.ylabel((', '.join([xlateLabel(_p[0]) for _p in ySpecs2])) if p not in xylabels or not xylabels[p][1] else xylabels[p][1])
+				pylab.xlabel(xlateLabel(pStrip) if (p not in xylabels or not xylabels[p][0]) else xylabels[p][0])
+				## should be done for y2 as well, but in that case the 10^.. label goes to y1 axis (bug in matplotlib, present in versions .99--1.0.5, and possibly beyond)
+				if scientific: pylab.ticklabel_format(style='sci',scilimits=(0,0),axis='both')
+			else:
+				pylab.rcParams['lines.color']=origLinesColor
+				pylab.ylabel((', '.join([xlateLabel(_p[0]) for _p in ySpecs2])) if (p not in xylabels or len(xylabels[p])<3 or not xylabels[p][2]) else xylabels[p][2])
+			# if there are callable ySpecs, save them inside the axes object, so that the live updater can use those
+			if yNameFuncs:
+				axes=line.get_axes()
+				axes.yadeYNames,axes.yadeYFuncs,axes.yadeXName,axes.yadeLabelLoc=yNames,yNameFuncs,pStrip,labelLoc # prepend yade to avoid clashes
+		createLines(pStrip,plots_p_y1,isY1=True,y2Exists=len(plots_p_y2)>0)
 		if axesWd>0:
 			pylab.axhline(linewidth=axesWd,color='k')
 			pylab.axvline(linewidth=axesWd,color='k')
@@ -199,21 +225,8 @@
 		if len(plots_p_y2)>0:
 			# try to move in the color palette a little further (magenta is 5th): r,g,b,c,m,y,k
 			origLinesColor=pylab.rcParams['lines.color']; pylab.rcParams['lines.color']='m'
-			# create the y2 axis
-			pylab.twinx()
-			for d in plots_p_y2:
-				line,=pylab.plot(data[pStrip],data[d[0]],d[1])
-				scatterPt=[0,0] if len(data[pStrip])==0 else (data[pStrip][current],data[d[0]][current])
-				scatter=pylab.scatter(scatterPt[0] if not math.isnan(scatterPt[0]) else 0,scatterPt[1] if not math.isnan(scatterPt[1]) else 0,color=line.get_color())
-				currLineRefs.append(LineRef(line,scatter,data[pStrip],data[d[0]]))
-			# legend
-			l=pylab.legend([xlateLabel(_p[0]) for _p in plots_p_y2],loc=legendPositionSecondary)
-			if hasattr(l,'draggable'): l.draggable(True)
-			pylab.rcParams['lines.color']=origLinesColor
-			pylab.ylabel((', '.join([xlateLabel(_p[0]) for _p in plots_p_y2])) if p not in xylabels or len(xylabels[p])<3 or not xylabels[p][2] else xylabels[p][2])
-			## should be repeated for y2 axis, but in that case the 10^.. label goes to y1 axis (bug in matplotlib, it seems)
-			# if scientific: pylab.ticklabel_format(style='sci',scilimits=(0,0),axis='both')
-
+			pylab.twinx() # create the y2 axis
+			createLines(pStrip,plots_p_y2,isY1=False,y2Exists=True)
 		if 'title' in O.tags.keys(): pylab.title(O.tags['title'])
 
 
@@ -228,6 +241,26 @@
 			l.update()
 			figs.add(l.line.get_figure())
 			axes.add(l.line.get_axes())
+		# find callables in y specifiers, create new lines if necessary
+		for ax in axes:
+			if not hasattr(ax,'yadeYFuncs') or not ax.yadeYFuncs: continue # not defined of empty
+			yy=set();
+			for f in ax.yadeYFuncs: yy.update(f())
+			#print 'callables y names:',yy
+			news=yy-ax.yadeYNames
+			if not news: continue
+			for new in news:
+				print 'yade.plot: creating new line for',new
+				if not new in data.keys(): data[new]=[] # create data entry if necessary
+				line,=ax.plot(data[ax.yadeXName],data[new],label=xlateLabel(new)) # no line specifier
+				scatterPt=(0 if len(data[ax.yadeXName])==0 and not math.isnan(data[ax.yadeXName]) else data[ax.yadeXName][current]),(0 if len(data[new])==0 and not math.isnan(data[new][current]) else data[new][current])
+				scatter=ax.scatter(scatterPt[0],scatterPt[1],color=line.get_color())
+				currLineRefs.append(LineRef(line,scatter,data[ax.yadeXName],data[new]))
+				ax.yadeYNames.add(new)
+				ax.set_ylabel(ax.get_ylabel()+(', ' if ax.get_ylabel() else '')+xlateLabel(new))
+			# it is possible that the legend has not yet been created
+			l=ax.legend(loc=ax.yadeLabelLoc)
+			if hasattr(l,'draggable'): l.draggable(True)
 		if autozoom:
 			for ax in axes:
 				try:

=== modified file 'py/utils.py'
--- py/utils.py	2010-12-21 12:19:50 +0000
+++ py/utils.py	2010-12-23 14:01:33 +0000
@@ -37,17 +37,15 @@
 
 		>>> utils.loadVars('mark')
 
-	and they will be defined in the yade.params.\ *mark* module
-
-	m*==True, variables will be loaded immediately after saving. That effectively makes *\*\*kw* available in builtin namespace.
+	and they will be defined in the yade.params.\ *mark* module. The *loadNow* parameter calls :yref:`yade.utils.loadVars` after saving automatically.
 	"""
 	import cPickle
 	Omega().tags['pickledPythonVariablesDictionary'+mark]=cPickle.dumps(kw)
 	if loadNow: loadVars(mark)
 
 def loadVars(mark=None):
-	"""Load variables from saveVars, which are saved inside the simulation.
-	If mark==None, all save variables are loaded. Otherwise only those with
+	"""Load variables from :yref:`yade.utils.saveVars`, which are saved inside the simulation.
+	If ``mark==None``, all save variables are loaded. Otherwise only those with
 	the mark passed."""
 	import cPickle, types, sys, warnings
 	def loadOne(d,mark=None):
@@ -98,12 +96,6 @@
 	"""
 	return [e for e in Omega().engines if e.__class__.__name__==name][0]
 
-def downCast(obj,newClassName):
-	"""Cast given object to class deriving from the same yade root class and copy all parameters from given object.
-	Obj should be up in the inheritance tree, otherwise some attributes may not be defined in the new class."""
-	return obj.__class__(newClassName,dict([ (key,obj[key]) for key in obj.keys() ]))
-
-
 def defaultMaterial():
 	"""Return default material, when creating bodies with :yref:`yade.utils.sphere` and friends, material is unspecified and there is no shared material defined yet. By default, this function returns::
 
@@ -263,14 +255,13 @@
 	else: pos2=position
 	b.state.pos=b.state.refPos=pos2
 	b.dynamic=False
-	b.aspherical=True # never used, since the wall is not dynamic
+	b.aspherical=False # wall never moves dynamically
 	b.mask=mask
 	return b
 
 def facet(vertices,dynamic=False,wire=True,color=None,highlight=False,noBound=False,material=-1,mask=1):
 	"""Create facet with given parameters.
 
-	:Parameters:
 	:param [Vector3,Vector3,Vector3] vertices: coordinates of vertices in the global coordinate system.
 	:param bool wire: if ``True``, facets are shown as skeleton; otherwise facets are filled
 	:param bool noBound: set :yref:`Body.bounded`
@@ -284,8 +275,7 @@
 	_commonBodySetup(b,0,Vector3(0,0,0),material,noBound=noBound)
 	b.state.pos=b.state.refPos=center
 	b.dynamic=dynamic
-	#b.aspherical=True
-	b.aspherical=False # FIXME: is it reasonably for a facets?
+	b.aspherical=False # mass and inertia are 0 anyway; fell free to change to ``True`` if needed
 	b.mask=mask
 	return b
 
@@ -474,7 +464,9 @@
 	.. math:: Z_m=\frac{2C-N_1}{N-N_0-N_1}
 
 	:param cutoff: cut some relative part of the sample's bounding box away.
-	:param skipFree: 	"""
+	:param skipFree: see above.
+	
+"""
 	if cutoff==0 and not skipFree: return 2*O.interactions.countReal()*1./len(O.bodies)
 	else:
 		nums,counts=bodyNumInteractionsHistogram(aabbExtrema(cutoff))
@@ -582,15 +574,10 @@
 
 #. Find the bodies that are on the negative/positive boundary, to which the straining condition should be applied.
 
-:parameters:
-	`filename`:
-		if given, spheres will be loaded from this file (ASCII format); if not, current simulation will be used.
-	`areaSection`:
-		number of section that will be used to estimate cross-section
-	`axis`:
-		if given, force strained axis, rather than computing it from predominant length
-
-:return: dictionary with keys 'negIds', 'posIds', 'axis', 'area'.
+:param filename: if given, spheres will be loaded from this file (ASCII format); if not, current simulation will be used.
+:param float areaSection: number of section that will be used to estimate cross-section
+:param ∈{0,1,2} axis: if given, force strained axis, rather than computing it from predominant length
+:return: dictionary with keys ``negIds``, ``posIds``, ``axis``, ``area``.
 
 .. warning::
 	The function :yref:`yade.utils.approxSectionArea` uses convex hull algorithm to find the area, but the implementation is reported to be *buggy* (bot works in some cases). Always check this number, or fix the convex hull algorithm (it is documented in the source, see :ysrc:`py/_utils.cpp`).
@@ -729,7 +716,7 @@
 	'Block the simulation if running inside a batch. Typically used at the end of script so that it does not finish prematurely in batch mode (the execution would be ended in such a case).'
 	if runningInBatch(): O.wait()
 
-def readParamsFromTable(tableFileLine=None,noTableOk=False,unknownOk=False,**kw):
+def readParamsFromTable(tableFileLine=None,noTableOk=True,unknownOk=False,**kw):
 	"""
 	Read parameters from a file and assign them to __builtin__ variables.
 
@@ -742,29 +729,21 @@
 		val2  val2  … # 2nd
 		…
 
-	Assigned tags:
+	Assigned tags (the ``description`` column is synthesized if absent,see :yref:`yade.utils.TableParamReader`); 
 
-	* *description* column is assigned to Omega().tags['description']; this column is synthesized if absent (see :yref:`yade.utils.TableParamReader`);
-	* ``Omega().tags['params']="name1=val1,name2=val2,…"``
-	* ``Omega().tags['defaultParams']="unassignedName1=defaultValue1,…"``
-	* ``Omega().tags['d.id']=O.tags['id']+'.'+O.tags['description']``
-	* ``Omega().tags['id.d']=O.tags['description']+'.'+O.tags['id']``
+		O.tags['description']=…                                      # assigns the description column; might be synthesized
+		O.tags['params']="name1=val1,name2=val2,…"                   # all explicitly assigned parameters
+		O.tags['defaultParams']="unassignedName1=defaultValue1,…"    # parameters that were left at their defaults
+		O.tags['d.id']=O.tags['id']+'.'+O.tags['description']
+		O.tags['id.d']=O.tags['description']+'.'+O.tags['id']
 
 	All parameters (default as well as settable) are saved using :yref:`yade.utils.saveVars`\ ``('table')``.
 
-	:parameters:
-		`tableFile`:
-			text file (with one value per blank-separated columns)
-		`tableLine`:
-			number of line where to get the values from.
-		`noTableOk`: bool
-			do not raise exception if the file cannot be open; use default values
-		`unknownOk`: bool
-			do not raise exception if unknown column name is found in the file; assign it as well
-
-	:return:
-		number of assigned parameters.
-
+	:param tableFile: text file (with one value per blank-separated columns)
+	:param int tableLine: number of line where to get the values from
+	:param bool noTableOk: if False, raise exception if the file cannot be open; use default values otherwise
+	:param bool unknownOk: do not raise exception if unknown column name is found in the file, and assign it as well
+	:return: number of assigned parameters
 	"""
 	tagsParams=[]
 	# dictParams is what eventually ends up in yade.params.table (default+specified values)

=== modified file 'scripts/test/energy.py'
--- scripts/test/energy.py	2010-12-21 12:19:50 +0000
+++ scripts/test/energy.py	2010-12-23 14:01:33 +0000
@@ -1,23 +1,17 @@
 from yade import pack,plot
 
-utils.readParamsFromTable(useL3Geom=True,nonviscDamp=0,frictAngle=0,noTableOk=True)
+utils.readParamsFromTable(useL3Geom=True,nonviscDamp=0,frictAngle=0,useClumps=False,noTableOk=True)
 from yade.params import table
 
-sp=pack.SpherePack();
-# bunch of balls, with an infinite plane just underneath
 if 1:
-	if 1: sp.makeCloud((0,0,0),(1,1,1),.05,.5);
+	sp=pack.SpherePack();
+	# bunch of balls, with an infinite plane just underneath
+	if not table.useClumps: sp.makeCloud((0,0,0),(1,1,1),.05,.5);
 	# use clumps of 2 spheres instead, to have rotation without friction 
 	else: sp.makeClumpCloud((0,0,0),(1,1,1),[pack.SpherePack([((0,0,0),.05),((0,0,.08),.02)])],periodic=False)
 	sp.toSimulation()
-	#if useL3Geom:
-	#O.bodies.append(pack.regularHexa(pack.inAlignedBox((-1,-1,-1),(2,2,0)),.5,-.2,dynamic=False))
-	O.bodies.append(utils.wall(position=0,axis=2))
-else:
-	#O.bodies.append(utils.sphere((0,0,1),.2))
-	O.bodies.append(utils.sphere((0,0,0.1035818200000008),.2))
-	O.bodies[0].state.vel,O.bodies[0].state.accel=(0,0,-4.188869999999981),(0,0,-9.8100000000000005)
-	O.bodies.append(utils.sphere((0.30000000000000004,0.19282032302755092,-0.5),.5,dynamic=False))
+else: O.bodies.append(utils.sphere((0,0,2),radius=.5)) # one single bouncing ball
+O.bodies.append(utils.wall(position=0,axis=2,sense=1))
 
 O.engines=[
 	ForceResetter(),
@@ -27,7 +21,7 @@
 	NewtonIntegrator(damping=table.nonviscDamp,kinSplit=True),
 	PyRunner(iterPeriod=1,command='addPlotData()'),
 ]
-O.dt=.2*utils.PWaveTimeStep()
+O.dt=.1*utils.PWaveTimeStep()
 
 def addPlotData():
 	Ek,maxId=utils.kineticEnergy(findMaxId=True)
@@ -36,25 +30,10 @@
 # turn on energy tracking
 O.trackEnergy=True
 O.saveTmp()
-# run a bit to have all energy categories in O.energy.keys().
-
-# The number of steps when all energy contributions are already non-zero
-# is only empirical; you can always force replot by redefining plot.plots
-# as below, closing plots and running plot.plot() again
-#
-# (unfortunately even if we were changing plot.plots periodically,
-# plots would not pick up changes in plot.plots during live plotting)
-#O.run(427,True)  
-#print O.bodies[0].state.pos,O.bodies[0].state.vel,O.bodies[0].state.accel
-O.run(5,True)
-#print list(set(['total',]+O.energy.keys()+['gravWork','kinRot','kinTrans','plastDissip']))+[None,'maxId']
-plot.plots={'i':list(set(['total',]+O.energy.keys()+['gravWork','kinRot','kinTrans','plastDissip','elastPotential']))} #+[None,'maxId']}
-#plot.plots={'i':['total',]+O.energy.keys()+[None,'maxId']}
-
-#O.run()
-
-from yade import timing
-O.timingEnabled=True
-O.run(20000,True)
-timing.stats()
+# the callable should return list of strings, plots will be updated automatically
+plot.plots={'i':[O.energy.keys,None,'total']}
+
+#from yade import timing
+#O.timingEnabled=True
+#timing.stats()
 plot.plot(subPlots=True)
  Thread Previous • Date Previous • • Thread Next