yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-03-31 22:23:18 +0200 (Tue, 31 Mar 2009)
New Revision: 1741

Added:
   trunk/scripts/test/interpolating-force.py
Modified:
   trunk/extra/Brefcom.cpp
   trunk/extra/Brefcom.hpp
   trunk/gui/py/pyAttrUtils.hpp
Log:
1. A few fixes in Brefcom.
2. Make python wrapper correctly handle long long attributes.
3. Added test script for InterpolatingDirectedForceEngine


Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp	2009-03-31 18:56:32 UTC (rev 1740)
+++ trunk/extra/Brefcom.cpp	2009-03-31 20:23:18 UTC (rev 1741)
@@ -134,11 +134,17 @@
 
 CREATE_LOGGER(ef2_Spheres_Brefcom_BrefcomLaw);
 
+long BrefcomContact::cummBetaIter=0, BrefcomContact::cummBetaCount=0;
+
 Real BrefcomContact::solveBeta(const Real c, const Real N){
+	#ifdef YADE_DEBUG
+		cummBetaCount++;
+	#endif
 	const int maxIter=20;
 	const Real maxError=1e-12;
 	Real f, ret=0.;
 	for(int i=0; i<maxIter; i++){
+		cummBetaIter++;
 		Real aux=c*exp(N*ret)+exp(ret);
 		f=log(aux);
 		if(fabs(f)<maxError) return ret;
@@ -152,7 +158,7 @@
 Real BrefcomContact::computeDmgOverstress(Real dt){
 	if(dmgStrain>=epsN*omega){ // unloading, no viscous stress
 		dmgStrain=epsN*omega;
-		LOG_DEBUG("Unloading, no viscous overstress");
+		LOG_DEBUG("Elastic/unloading, no viscous overstress");
 		return 0.;
 	}
 	Real c=epsCrackOnset*(1-omega)*pow(dmgTau/dt,dmgRateExp)*pow(epsN*omega-dmgStrain,dmgRateExp-1.);

Modified: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp	2009-03-31 18:56:32 UTC (rev 1740)
+++ trunk/extra/Brefcom.hpp	2009-03-31 20:23:18 UTC (rev 1741)
@@ -89,6 +89,8 @@
 		bool neverDamage;
 		/*! cummulative plastic strain measure (scalar) on this contact */
 		Real epsPlSum;
+		//! debugging, to see convergence rate
+		static long cummBetaIter, cummBetaCount;
 
 		/*! auxiliary variable for visualization, recalculated in BrefcomLaw at every iteration */
 		// FIXME: Fn and Fs are stored as Vector3r normalForce, shearForce in NormalShearInteraction 
@@ -122,6 +124,9 @@
 			(plRateExp)
 			(transStrainCoeff)
 
+			(cummBetaIter)
+			(cummBetaCount)
+
 			(kappaD)
 			(neverDamage)
 			(epsT)

Modified: trunk/gui/py/pyAttrUtils.hpp
===================================================================
--- trunk/gui/py/pyAttrUtils.hpp	2009-03-31 18:56:32 UTC (rev 1740)
+++ trunk/gui/py/pyAttrUtils.hpp	2009-03-31 20:23:18 UTC (rev 1741)
@@ -113,7 +113,7 @@
 					if      (any_cast<string*>(&instance)) { desc.type=AttrAccess::STRING; goto found; }
 					else if (any_cast<bool*>(&instance)) { desc.type=AttrAccess::BOOL; goto found; }
 					else if (any_cast<Real*>(&instance) || any_cast<long double*>(&instance) || any_cast<double*>(&instance) || any_cast<float*>(&instance)) { desc.type=AttrAccess::FLOAT; goto found;}
-					else if (any_cast<int*>(&instance) || any_cast<unsigned int*>(&instance) || any_cast<long*>(&instance) || any_cast<unsigned long*>(&instance)) {desc.type=AttrAccess::INTEGER; goto found; }
+					else if (any_cast<int*>(&instance) || any_cast<unsigned int*>(&instance) || any_cast<long*>(&instance) || any_cast<unsigned long*>(&instance) || any_cast<long long*>(&instance) || any_cast<unsigned long long*>(&instance)) {desc.type=AttrAccess::INTEGER; goto found; }
 					else if (any_cast<vector<string>*>(&instance)) { desc.type=AttrAccess::SEQ_STRING; goto found; }
 				#if 0
 					else if (any_cast<vector<Vector3r>*>(&instance)) { cerr<<"WWWWWWWWWWWWW"<<endl;}

Added: trunk/scripts/test/interpolating-force.py
===================================================================
--- trunk/scripts/test/interpolating-force.py	2009-03-31 18:56:32 UTC (rev 1740)
+++ trunk/scripts/test/interpolating-force.py	2009-03-31 20:23:18 UTC (rev 1741)
@@ -0,0 +1,69 @@
+# 2009 Václav Šmilauer <eudoxos@xxxxxxxx>
+# encoding: utf-8
+
+# script showing how to use InterpolatingDirectedForceEngine,
+# simply pushing one free sphere agains a fixed one with varying force
+#
+# The force evolution is sine wave, but it could really be any data
+
+from numpy import arange
+
+nPulses=2 # run for total of 2 pulses
+freq=10. # 5 pulses per second
+times=arange(0,1/freq,.01/freq) # generate 100 points equally distributed over the period (can be much more)
+maxMag=1e10 # maximum magnitude of applied force
+magnitudes=[.5*maxMag*(sin(t*(freq*2*pi))+1) for t in times] # generate points on sine wave over 1 period, but shifted up to be ∈(0,2)
+
+O.engines=[
+	StandAloneEngine('PhysicalActionContainerReseter'),
+	MetaEngine('BoundingVolumeMetaEngine',[EngineUnit('InteractingSphere2AABB'),]),
+	StandAloneEngine('PersistentSAPCollider'),
+	InteractionDispatchers(
+		[EngineUnit('InteractingSphere2InteractingSphere4SpheresContactGeometry')],
+		[EngineUnit('SimpleElasticRelationships')],
+		[EngineUnit('ef2_Spheres_Elastic_ElasticLaw'),]
+	),
+	# subscribedBodies: what bodies is force applied to
+	# direction: direction of the force (normalized automatically), constant
+	# magnitudes: series of force magnitude
+	# times: time points at which magnitudes are defined
+	# wrap: continue from t0 once t_last is reached
+	# label: automatically defines python variable of that name pointing to this engine
+	DeusExMachina('InterpolatingDirectedForceEngine',{'subscribedBodies':[1],'direction':[0,0,-1],'magnitudes':magnitudes,'times':times,'wrap':True,'label':'forcer'}),
+	# without damping, the interaction never stabilizes and oscillates wildly… try it
+	DeusExMachina('NewtonsDampedLaw',{'damping':0.01}),
+	# collect some data to plot periodically (every 50 steps)
+	StandAloneEngine('PeriodicPythonRunner',{'iterPeriod':50,'command':'myAddPlotData()','label':'plotDataCollector'})
+]
+
+O.bodies.append([
+	utils.sphere([0,0,0],1,dynamic=False,color=[1,0,0]),
+	utils.sphere([0,0,2],1,color=[0,1,0])
+])
+
+# elastic timestep
+O.dt=utils.PWaveTimeStep()
+
+# callback for plotDataCollector
+import yade.plot as yp
+def myAddPlotData():
+	yp.addData({'t':O.time,'F_applied':forcer['force'],'supportReaction':O.actions.f(0)[2]})
+
+O.saveTmp('init')
+
+# try open 3d view, if not running in pure console mode
+try:
+	import yade.qt
+	yade.qt.View()
+except ImportError: pass
+
+# run so many steps such that prescribed number of pulses is done
+O.run(int((nPulses/freq)/O.dt),True)
+
+# plot the time-series of force magnitude
+import pylab
+pylab.plot(times,magnitudes,label='Force magnitude over 1 pulse'); pylab.legend(('Force magnitude',)); pylab.xlabel('t'); pylab.grid(True)
+# plot some recorded data
+yp.plots={'t':('F_applied','supportReaction')}
+yp.plot()
+