yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2380
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Thu 2010-07-22 15:52:50 +0200
message:
  1. add noShow to uytils.plotDirections 
  2. Some fixes for ScGeom+Cpm (still problematic)
modified:
  examples/concrete/dem3dof-scgeom.table
  examples/concrete/interaction-histogram.py
  pkg/dem/meta/ConcretePM.cpp
  py/utils.py
  scripts/test/cpm-dem3dof-scgeom.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 'examples/concrete/dem3dof-scgeom.table'
--- examples/concrete/dem3dof-scgeom.table	2010-07-18 13:20:57 +0000
+++ examples/concrete/dem3dof-scgeom.table	2010-07-22 13:52:50 +0000
@@ -10,6 +10,10 @@
 # but are less efficient in multi-threading
 # (that might be due to the implementation of Law2_ScGeom_CpmPhys_Cpm, though)
 #
+# Run the batch using
+#
+# $ yade-trunk-multi --gnuplot=d3d-scg.gnuplot dem3dof-scgeom.table uniax.py
+#
 scGeom	dtSafety
 False	.8
 False .2

=== modified file 'examples/concrete/interaction-histogram.py'
--- examples/concrete/interaction-histogram.py	2010-06-03 09:28:04 +0000
+++ examples/concrete/interaction-histogram.py	2010-07-22 13:52:50 +0000
@@ -9,7 +9,7 @@
 import os.path
 
 # run uniax.py to get this file
-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)
 
@@ -25,7 +25,7 @@
 	angles.append(angle)
 	forces.append(force)
 # easier: plain histogram
-pylab.hist(angles,weights=forces,bins=20)
+#pylab.hist(angles,weights=forces,bins=20)
 
 
 # polar histogram
@@ -41,6 +41,6 @@
 
 # predefined function
 pylab.figure()
-utils.plotDirections()
+utils.plotDirections(noShow=True).savefig('/tmp/a.pdf')
 
 pylab.show()

=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp	2010-07-21 08:00:51 +0000
+++ pkg/dem/meta/ConcretePM.cpp	2010-07-22 13:52:50 +0000
@@ -239,33 +239,35 @@
 
 CREATE_LOGGER(Law2_ScGeom_CpmPhys_Cpm);
 void Law2_ScGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I){
-	ScGeom* contGeom=static_cast<ScGeom*>(_geom.get());
+	ScGeom* geom=static_cast<ScGeom*>(_geom.get());
 	CpmPhys* BC=static_cast<CpmPhys*>(_phys.get());
-	// FIXME: better way to get current positions? Needed for
-	// (1) getting intial equilibrium distance (not stored in ScGeom and distFactor not accessible from here)
-	// (2) applying contact forces back on particles
-	const Vector3r& pos1(Body::byId(I->getId1(),scene)->state->pos); const Vector3r& pos2(Body::byId(I->getId2(),scene)->state->pos);
-	Real refLength;
 	// just the first time
 	if(I->isFresh(scene)){
 		// done with real sphere radii
-		Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
+		Real minRad=(geom->refR1<=0?geom->refR2:(geom->refR2<=0?geom->refR1:min(geom->refR1,geom->refR2)));
 		BC->crossSection=Mathr::PI*pow(minRad,2);
+		// FIXME: a better way to get current positions? Needed for getting intial equilibrium distance (not stored in ScGeom and distFactor not accessible from here)
+		const Vector3r& pos1(Body::byId(I->getId1(),scene)->state->pos); const Vector3r& pos2(Body::byId(I->getId2(),scene)->state->pos);
 		// scale sphere's radii to effective radii (intial equilibrium)
-		Real refLength=(pos2-pos1).norm(); Real distCurr=contGeom->radius1+contGeom->radius2;
-		contGeom->radius1*=refLength/distCurr; contGeom->radius2*=refLength/distCurr;
-		contGeom->penetrationDepth=0;
+		Real refLength=(pos2-pos1).norm(); Real distCurr=geom->radius1+geom->radius2;
+		geom->radius1*=refLength/distCurr; geom->radius2*=refLength/distCurr;
+		geom->penetrationDepth=0;
+		geom->contactPoint=pos1+(geom->radius1/refLength)*(pos2-pos1);
 		BC->kn=BC->crossSection*BC->E/refLength;
 		BC->ks=BC->crossSection*BC->G/refLength;
 	}
+	// compute particle positions from contact point and penetrationDepth
+	Vector3r pos1=geom->contactPoint-(geom->radius1-.5*geom->penetrationDepth)*geom->normal;
+	Vector3r pos2=geom->contactPoint+(geom->radius2-.5*geom->penetrationDepth)*geom->normal;
+	Real refLength=geom->radius1+geom->radius2;
 	// shorthands
 	Real& epsN(BC->epsN);
 	Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(Law2_ScGeom_CpmPhys_Cpm::omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); Real& omega(BC->omega); Real& sigmaN(BC->sigmaN);  Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
 	const bool& isCohesive(BC->isCohesive);
-	epsN=contGeom->penetrationDepth/refLength;
+	epsN=-geom->penetrationDepth/refLength;
 	
-	epsT=contGeom->rotate(epsT);
-	epsT+=contGeom->shearIncrement()/refLength; 
+	epsT=geom->rotate(epsT);
+	epsT+=geom->shearIncrement()/refLength; 
 
 	// simplified public model
 	epsN+=BC->isoPrestress/E;
@@ -297,10 +299,10 @@
 		return;
 	}
 
-	Fn=sigmaN*crossSection; BC->normalForce=Fn*contGeom->normal;
+	Fn=sigmaN*crossSection; BC->normalForce=Fn*geom->normal;
 	Fs=sigmaT*crossSection; BC->shearForce=Fs;
 
-	applyForceAtContactPoint(BC->normalForce+BC->shearForce,contGeom->contactPoint,I->getId1(),pos1,I->getId2(),pos2);
+	applyForceAtContactPoint(BC->normalForce+BC->shearForce,geom->contactPoint,I->getId1(),pos1,I->getId2(),pos2);
 }
 
 

=== modified file 'py/utils.py'
--- py/utils.py	2010-07-18 19:21:08 +0000
+++ py/utils.py	2010-07-22 13:52:50 +0000
@@ -509,9 +509,12 @@
 	pylab.ylabel('Body count')
 	pylab.show()
 
-def plotDirections(aabb=(),mask=0,bins=20,numHist=True):
+def plotDirections(aabb=(),mask=0,bins=20,numHist=True,noShow=False):
 	"""Plot 3 histograms for distribution of interaction directions, in yz,xz and xy planes and
-	(optional but default) histogram of number of interactions per body."""
+	(optional but default) histogram of number of interactions per body.
+	
+	:returns: If *noShow* is ``False``, displays the figure and returns nothing. If *noShow*, the figure object is returned without being displayed (works the same way as :yref:`yade.plot.plot`).
+	"""
 	import pylab,math
 	from yade import utils
 	for axis in [0,1,2]:
@@ -530,7 +533,8 @@
 		pylab.xlabel('Interactions per body (avg. %g)'%avg)
 		pylab.axvline(x=avg,linewidth=3,color='r')
 		pylab.ylabel('Body count')
-	pylab.show()
+	if noShow: return pylab.gcf()
+	else: pylab.show()
 
 
 

=== modified file 'scripts/test/cpm-dem3dof-scgeom.py'
--- scripts/test/cpm-dem3dof-scgeom.py	2010-07-18 11:59:57 +0000
+++ scripts/test/cpm-dem3dof-scgeom.py	2010-07-22 13:52:50 +0000
@@ -4,7 +4,7 @@
 # one is handled by Law2_Dem3DofGeom_CpmPhys_Cpm and the other by Law2_ScGeom_CpmPhys_Cpm
 
 # move the second sphere tangentially or rotate it, pick [0] or [1]
-mode=['mov','rot'][0]
+mode=['mov','rot'][1]
 # number of steps to do (some influence on the incremental computation)
 nSteps=100