yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2680
committer: Sergei D. sj2001@xxxxxxxxx
branch nick: yade
timestamp: Wed 2011-01-26 15:39:07 +0300
message:
  1. Rename NozzleFactory to SpheresFactory; add derived CircularFactory and QuadroFactory. Example is scripts/test/spheresFactory.py
  2. small fixes
removed:
  pkg/common/SpheresFactory.hpp
renamed:
  pkg/dem/NozzleFactory.cpp => pkg/dem/SpheresFactory.cpp
  pkg/dem/NozzleFactory.hpp => pkg/dem/SpheresFactory.hpp
  scripts/test/shots.py => scripts/test/spheresFactory.py
modified:
  py/geom.py
  scripts/test/clump-hopper-viscoelastic.py
  pkg/dem/SpheresFactory.cpp
  pkg/dem/SpheresFactory.hpp
  scripts/test/spheresFactory.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

=== removed file 'pkg/common/SpheresFactory.hpp'
--- pkg/common/SpheresFactory.hpp	2011-01-21 08:14:28 +0000
+++ pkg/common/SpheresFactory.hpp	1970-01-01 00:00:00 +0000
@@ -1,121 +0,0 @@
-/*************************************************************************
-*  Copyright (C) 2009 by Sergei Dorofeenko				 				 *
-*  sega@xxxxxxxxxxxxxxxx                                                 *
-*                                                                        *
-*  This program is free software; it is licensed under the terms of the  *
-*  GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-#pragma once
-
-#include <yade/pkg/common/PeriodicEngines.hpp>
-#include <yade/pkg/common/Dispatching.hpp>
-#include <yade/pkg/common/Collider.hpp>
-#include <yade/core/Scene.hpp>
-#include <vector>
-#include <string>
-
-using namespace std;
-
-/// @brief Produces spheres over the course of a simulation. 
-class SpheresFactory : public PeriodicEngine {
-public:
-
-	SpheresFactory();
-	virtual ~SpheresFactory();
-
-	/// @brief Create one sphere per call.
-	virtual void action();
-
-	/// @brief The geometry of the section on which spheres will be placed. 
-	vector<Body::id_t> factoryFacets; 
-
-    /// @brief Factory section may be a surface or volume (convex). 
-    /// By default it is a surface. To make its a volume set volumeSection=true
-    bool volumeSection;
-
-	/// @brief Max attempts to place sphere.
-	/// If placing the sphere in certain random position would cause an overlap with any other physical body in the model, SpheresFactory will try to find another position. Default 20 attempts allow.
-	int maxAttempts; 
-
-	/// @brief Mean radius of spheres.
-	Real radius; 
-
-	/// @brief Half size of a radii distribution interval.
-	/// New sphere will have random radius within the range radius±radiusRange.
-	Real radiusRange;
-
-	/// @brief Mean velocity of spheres.
-	Vector3r velocity;
-
-	/// @brief Half size of a velocities distribution interval.
-	/// New sphere will have random velocity within the range velocity±velocityRange.
-	Vector3r velocityRange;
-	
-	/// @brief Mean angularVelocity of spheres.
-	Vector3r angularVelocity;
-
-	/// @brief Half size of a angularVelocity distribution interval.
-	/// New sphere will have random angularVelocity within the range angularVelocity±angularVelocityRange.
-	Vector3r angularVelocityRange;
-
-	/// @brief Young modulus.
-	Real young;
-	/// @brief Poisson ratio.
-	Real poisson;
-	/// @brief Density of material.
-	Real density;
-	/// @brief Friction angle (radians).
-	Real frictionAngle;
-	/// @brief Color.
-	Vector3r color;
-
-	/// @brief python's function for a spheres creation
-	string pySpheresCreator;
-
-private:
-	/// @brief Pointer to Collider.
-	/// It is necessary in order to probe the bounding volume for new sphere.
-	Collider* bI;
-	
-	/// @brief Pointer to IGeomDispatcher.
-	/// It is necessary in order to detect a real overlap with other bodies.
-	IGeomDispatcher* iGME;
-
-	bool first_run;
-
-    Vector3r generatePositionOnSurface();
-    Vector3r generatePositionInVolume();
-
-	void createSphere(shared_ptr<Body>& body, const Vector3r& position, Real r);
-
-	typedef	boost::variate_generator<boost::minstd_rand,boost::uniform_int<> > RandomInt;
-	shared_ptr<RandomInt> randomFacet;
-
-	static boost::variate_generator<boost::mt19937,boost::uniform_real<> > 
-		randomUnit;
-	static boost::variate_generator<boost::mt19937,boost::uniform_real<> >
-		randomSymmetricUnit;
-
-	DECLARE_LOGGER;
-
-	REGISTER_ATTRIBUTES(PeriodicEngine,
-			(factoryFacets)
-			(volumeSection)
-			(maxAttempts)
-			(radius)
-			(radiusRange)
-			(velocity)
-			(velocityRange)
-			(angularVelocity)
-			(angularVelocityRange)
-			(young)
-			(poisson)
-			(density)
-			(frictionAngle)
-			(pySpheresCreator)
-			(color))
-	REGISTER_CLASS_AND_BASE(SpheresFactory, GlobalEngine);
-};
-REGISTER_SERIALIZABLE(SpheresFactory);
-
-

=== renamed file 'pkg/dem/NozzleFactory.cpp' => 'pkg/dem/SpheresFactory.cpp'
--- pkg/dem/NozzleFactory.cpp	2010-12-20 07:19:47 +0000
+++ pkg/dem/SpheresFactory.cpp	2011-01-26 12:39:07 +0000
@@ -1,13 +1,16 @@
 
-#include<yade/pkg/dem/NozzleFactory.hpp>
+#include<yade/pkg/dem/SpheresFactory.hpp>
 #include<yade/pkg/common/Sphere.hpp>
 #include <boost/random/linear_congruential.hpp>
 #include <boost/random/uniform_real.hpp>
 #include <boost/random/variate_generator.hpp>
 
 
-YADE_PLUGIN((NozzleFactory)(DragForceApplier));
-CREATE_LOGGER(NozzleFactory);
+//YADE_PLUGIN((SpheresFactory)(DragForceApplier));
+YADE_PLUGIN((SpheresFactory)(CircularFactory)(QuadroFactory)(DragForceApplier));
+CREATE_LOGGER(SpheresFactory);
+CREATE_LOGGER(CircularFactory);
+CREATE_LOGGER(QuadroFactory);
 
 void DragForceApplier::action(){
 	FOREACH(const shared_ptr<Body>& b, *scene->bodies){
@@ -19,19 +22,24 @@
 	}
 }
 
-void NozzleFactory::action(){
-	// initialize random number generator with time seed
-	static boost::minstd_rand randGen(TimingInfo::getNow(/* get the number even if timing is disabled globally */ true));
-	static boost::variate_generator<boost::minstd_rand&, boost::uniform_real<Real> > randomUnit(randGen, boost::uniform_real<Real>(0,1));
+// initialize random number generator with time seed
+static boost::minstd_rand randGen(TimingInfo::getNow(/* get the number even if timing is disabled globally */ true));
+static boost::variate_generator<boost::minstd_rand&, boost::uniform_real<Real> > randomUnit(randGen, boost::uniform_real<Real>(0,1));
+
+void SpheresFactory::pickRandomPosition(Vector3r&,Real){
+	LOG_FATAL("Engine "<<getClassName()<<" calling virtual method SpheresFactory::pickRandomPosition(). Please submit bug report at http://bugs.launchpad.net/yade.";);
+	throw std::logic_error("SpheresFactory::pickRandomPosition() called.");
+}
+
+void SpheresFactory::action(){
 
 	if(!collider){
 		FOREACH(const shared_ptr<Engine>& e, scene->engines){ collider=dynamic_pointer_cast<Collider>(e); if(collider) break; }
-		if(!collider) throw runtime_error("NozzleFactory: No Collider instance found in engines (needed for collision detection).");
+		if(!collider) throw runtime_error("SpheresFactory: No Collider instance found in engines (needed for collision detection).");
 	}
 	goalMass+=massFlowRate*scene->dt; // totalMass that we want to attain in the current step
 
 	normal.normalize();
-	const Quaternionr q(Quaternionr().setFromTwoVectors(Vector3r::UnitZ(),normal));
 
 	LOG_TRACE("totalMass="<<totalMass<<", goalMass="<<goalMass);
 
@@ -43,10 +51,7 @@
 		// until there is no overlap, pick a random position for the new particle
 		int attempt;
 		for(attempt=0; attempt<maxAttempt; attempt++){
-			Real angle=randomUnit()*2*Mathr::PI, rr=randomUnit()*(radius-r); // random polar coordinate inside the nozzle
-			c = center+q*Vector3r(cos(angle)*rr,sin(angle)*rr,0);
-			// this version places center in a box around the nozzle center (its size 1.15=2/sqrt(3) diagonal, which is not very exact
-			//c=center+Vector3r((randomUnit()-.5)*1.15*radius,(randomUnit()-.5)*1.15*radius,(randomUnit()-.5)*1.15*radius);
+			pickRandomPosition(c,r);
 			LOG_TRACE("Center "<<c);
 			Bound b; b.min=c-Vector3r(r,r,r); b.max=c+Vector3r(r,r,r);
 			vector<Body::id_t> collidingParticles=collider->probeBoundingVolume(b);
@@ -56,7 +61,7 @@
 			#endif
 		}
 		if(attempt==maxAttempt) {
-			if (silent) {massFlowRate=0;} 
+			if (silent) {massFlowRate=0;goalMass=0;LOG_INFO("Unable to place new sphere after "<<maxAttempt<<" attempts, SpheresFactory disabled.");} 
 			else {LOG_WARN("Unable to place new sphere after "<<maxAttempt<<" attempts, giving up.");}
 			return;
 		}
@@ -67,7 +72,7 @@
 
 		// create particle
 		int mId=(materialId>=0 ? materialId : scene->materials.size()+materialId);
-		if(mId<0 || (size_t) mId>=scene->materials.size()) throw std::invalid_argument(("NozzleFactory: invalid material id "+lexical_cast<string>(materialId)).c_str());
+		if(mId<0 || (size_t) mId>=scene->materials.size()) throw std::invalid_argument(("SpheresFactory: invalid material id "+lexical_cast<string>(materialId)).c_str());
 		const shared_ptr<Material>& material=scene->materials[mId];
 		shared_ptr<Body> b(new Body);
 		shared_ptr<Sphere> sphere(new Sphere); 
@@ -92,3 +97,15 @@
 	} 
 	//std::cout<<"mass flow rate: "<<totalMass<<endl;totalMass =0.0;
 };
+
+void CircularFactory::pickRandomPosition(Vector3r& c, Real r){
+	const Quaternionr q(Quaternionr().setFromTwoVectors(Vector3r::UnitZ(),normal));
+	Real angle=randomUnit()*2*Mathr::PI, rr=randomUnit()*(radius-r); // random polar coordinate inside the nozzle
+	Real l=(randomUnit()-0.5)*length;
+	c = center+q*Vector3r(cos(angle)*rr,sin(angle)*rr,0)+normal*l;
+}
+
+void QuadroFactory::pickRandomPosition(Vector3r& c, Real r){
+	const Quaternionr q(Quaternionr().setFromTwoVectors(Vector3r::UnitZ(),normal));
+	c=center+q*Vector3r((randomUnit()-.5)*2*(extents[0]-r),(randomUnit()-.5)*2*(extents[1]-r),(randomUnit()-.5)*2*(extents[2]-r));
+}

=== renamed file 'pkg/dem/NozzleFactory.hpp' => 'pkg/dem/SpheresFactory.hpp'
--- pkg/dem/NozzleFactory.hpp	2010-12-19 09:58:59 +0000
+++ pkg/dem/SpheresFactory.hpp	2011-01-26 12:39:07 +0000
@@ -4,31 +4,59 @@
 #include<yade/pkg/common/Collider.hpp>
 
 
-class NozzleFactory: public GlobalEngine {
+class SpheresFactory: public GlobalEngine {
 	shared_ptr<Collider> collider;
+	protected:
+		// Pick random position of a sphere. Should be override in derived engine.
+		virtual void pickRandomPosition(Vector3r&/*picked position*/, Real/*sphere's radius*/);
 	public:
 		virtual void action();
 	DECLARE_LOGGER;
-	YADE_CLASS_BASE_DOC_ATTRS(NozzleFactory,GlobalEngine,"Engine for spitting spheres based on mass flow rate, particle size distribution etc. The area where spehres are generated should be circular, given by radius, center and normal. For now, axis-aligned cube-shape corresponding is used instead, centered at *center* and with size $\\frac{2\\sqrt{3}}{3}$. Initial velocity of particles is given by *vMin*, *vMax*, the *massFlowRate* determines how many particles to generate at each step. When *goalMass* is attained or positive *maxParticles* is reached, the engine does not produce particles anymore.\n\nA sample script for this engine is in :ysrc:`scripts/shots.py`.",
+	YADE_CLASS_BASE_DOC_ATTRS(SpheresFactory,GlobalEngine,"Engine for spitting spheres based on mass flow rate, particle size distribution etc. Initial velocity of particles is given by *vMin*, *vMax*, the *massFlowRate* determines how many particles to generate at each step. When *goalMass* is attained or positive *maxParticles* is reached, the engine does not produce particles anymore. Geometry of the region should be defined in a derived engine by overriden SpheresFactory::pickRandomPosition(). \n\nA sample script for this engine is in :ysrc:`scripts/spheresFactory.py`.",
 		((Real,massFlowRate,NaN,,"Mass flow rate [kg/s]"))
 		((Real,rMin,NaN,,"Minimum radius of generated spheres (uniform distribution)"))
 		((Real,rMax,NaN,,"Maximum radius of generated spheres (uniform distribution)"))
 		((Real,vMin,NaN,,"Minimum velocity norm of generated spheres (uniform distribution)"))
 		((Real,vMax,NaN,,"Maximum velocity norm of generated spheres (uniform distribution)"))
 		((Real,vAngle,NaN,,"Maximum angle by which the initial sphere velocity deviates from the nozzle normal."))
-		((Real,radius,NaN,,"Radius of the nozzle"))
+		((Vector3r,normal,Vector3r(NaN,NaN,NaN),,"Spitting direction (and orientation of the region's geometry)."))
 		((int,materialId,-1,,"Shared material id to use for newly created spheres (can be negative to count from the end)"))
-		((Vector3r,center,Vector3r(NaN,NaN,NaN),,"Center of the nozzle"))
-		((Vector3r,normal,Vector3r(NaN,NaN,NaN),,"Spitting direction, i.e. normal of the circle where spheres are generated."))
 		((Real,totalMass,0,,"Mass of spheres that was produced so far. |yupdate|"))
 		((Real,goalMass,0,,"Total mass that should be attained at the end of the current step. |yupdate|"))
 		((int,maxParticles,100,,"The number of particles at which to stop generating new ones (regardless of massFlowRate"))
 		((int,numParticles,0,,"Cummulative number of particles produces so far |yupdate|"))
 		((int,maxAttempt,5000 ,,"Maximum number of attempts to position a new sphere randomly."))
-		((bool,silent,false ,,"If true no complain about excessing maxAttempt but disable factory (by set massFlowRate=0)."))
-	);
-};
-REGISTER_SERIALIZABLE(NozzleFactory);
+		((bool,silent,false ,,"If true no complain about excessing maxAttempt but disable the factory (by set massFlowRate=0)."))
+	);
+};
+REGISTER_SERIALIZABLE(SpheresFactory);
+
+class CircularFactory: public SpheresFactory {
+	protected:
+		virtual void pickRandomPosition(Vector3r&, Real);
+	public:
+		virtual ~CircularFactory(){};
+		DECLARE_LOGGER;
+		YADE_CLASS_BASE_DOC_ATTRS(CircularFactory,SpheresFactory,"Circular geometry of the SpheresFactory region. It can be disk (given by radius and center), or cylinder (given by radius, length and center).",
+		((Real,radius,NaN,,"Radius of the region"))
+		((Real,length,0,,"Length of the cylindrical region (0 by default)"))
+		((Vector3r,center,Vector3r(NaN,NaN,NaN),,"Center of the region"))
+	);
+};
+REGISTER_SERIALIZABLE(CircularFactory);
+
+class QuadroFactory: public SpheresFactory {
+	protected:
+		virtual void pickRandomPosition(Vector3r&, Real);
+	public:
+		virtual ~QuadroFactory(){};
+		DECLARE_LOGGER;
+		YADE_CLASS_BASE_DOC_ATTRS(QuadroFactory,SpheresFactory,"Quadro geometry of the SpheresFactory region, given by extents and center",
+		((Vector3r,extents,Vector3r(NaN,NaN,NaN),,"Extents of the region"))
+		((Vector3r,center,Vector3r(NaN,NaN,NaN),,"Center of the region"))
+	);
+};
+REGISTER_SERIALIZABLE(QuadroFactory);
 
 class DragForceApplier: public GlobalEngine{
 	public: virtual void action();

=== modified file 'py/geom.py'
--- py/geom.py	2011-01-13 10:21:10 +0000
+++ py/geom.py	2011-01-26 12:39:07 +0000
@@ -261,7 +261,7 @@
 		radiusBottom = radiusTop
 	# check zero dimentions
 	if (segmentsNumber<3): raise RuntimeError("The segmentsNumber should be at least 3");
-	if (height<=0): raise RuntimeError("The height should have the positive value");
+	if (height<0): raise RuntimeError("The height should have the positive value");
 	if angleRange==None: angleRange=(0,2*math.pi)
 	if (abs(angleRange[1]-angleRange[0])>2.0*math.pi): raise RuntimeError("The |angleRange| cannot be larger 2.0*math.pi");
 	

=== modified file 'scripts/test/clump-hopper-viscoelastic.py'
--- scripts/test/clump-hopper-viscoelastic.py	2010-12-26 15:42:43 +0000
+++ scripts/test/clump-hopper-viscoelastic.py	2011-01-26 12:39:07 +0000
@@ -14,10 +14,10 @@
 frictionAngle=radians(35)# 
 density=2700
 # facets material
-params=utils.getViscoelasticFromSpheresInteraction(10e3,tc,en,es)
-facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,**params)) # **params sets kn, cn, ks, cs
+params=utils.getViscoelasticFromSpheresInteraction(tc,en,es)
+facetMat=O.materials.append(ViscElMat(density=0,frictionAngle=frictionAngle,**params)) # **params sets kn, cn, ks, cs
 # default spheres material
-dfltSpheresMat=O.materials.append(ViscElMat(density=density,frictionAngle=frictionAngle)) 
+dfltSpheresMat=O.materials.append(ViscElMat(density=density,frictionAngle=frictionAngle, **params)) 
 
 O.dt=.1*tc # time step
 
@@ -50,39 +50,27 @@
 clumpColor=(0.0, 0.5, 0.5)
 for k,l in itertools.product(arange(0,10),arange(0,10)):
 	clpId,sphId=O.bodies.appendClumped([utils.sphere(Vector3(x0t+Rs*(k*4+2),y0t+Rs*(l*4+2),i*Rs*2+zt),Rs,color=clumpColor,material=dfltSpheresMat) for i in xrange(4)])
-	for id in sphId:
-		s=O.bodies[id]
-		p=utils.getViscoelasticFromSpheresInteraction(s.state.mass,tc,en,es)
-		s.mat.kn,s.mat.cn,s.mat.ks,s.mat.cs=p['kn'],p['cn'],p['ks'],p['cs']
-	#O.bodies[clpId].state.blockedDOFs='XYZ'
-	#O.bodies[clpId].state.blockedDOFs='xy'
 
 # ... and spheres
-spheresColor=(0.4, 0.4, 0.4)
-for k,l in itertools.product(arange(0,9),arange(0,9)):
-	sphAloneId=O.bodies.append( [utils.sphere( Vector3(x0t+Rs*(k*4+4),y0t+Rs*(l*4+4),i*Rs*2.3+zt),Rs,color=spheresColor,material=dfltSpheresMat) for i in xrange(4) ] )
-	for id in sphAloneId:
-		s=O.bodies[id]
-		p=utils.getViscoelasticFromSpheresInteraction(s.state.mass,tc,en,es)
-		s.mat.kn,s.mat.cn,s.mat.ks,s.mat.cs=p['kn'],p['cn'],p['ks'],p['cs']
-		#s.state.blockedDOFs='XYZ'
-		#s.state.blockedDOFs='xy'
+#spheresColor=(0.4, 0.4, 0.4)
+#for k,l in itertools.product(arange(0,9),arange(0,9)):
+	#sphAloneId=O.bodies.append( [utils.sphere( Vector3(x0t+Rs*(k*4+4),y0t+Rs*(l*4+4),i*Rs*2.3+zt),Rs,color=spheresColor,material=dfltSpheresMat) for i in xrange(4) ] )
 
 # Create engines
 O.engines=[
 	ForceResetter(),
-	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],nBins=5,sweepLength=.1*Rs),
+	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 	InteractionLoop(
 		[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
 		[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
 		[Law2_ScGeom_ViscElPhys_Basic()],
 	),
 	GravityEngine(gravity=[0,0,-9.81]),
-	NewtonIntegrator(damping=0,exactAsphericalRot=True),
+	NewtonIntegrator(damping=0,exactAsphericalRot=False),
 	#VTKRecorder(virtPeriod=0.01,fileName='/tmp/',recorders=['spheres','velocity','facets'])
 ]
 
-renderer = qt.Renderer()
+from yade import qt
 qt.View()
 O.saveTmp()
 

=== renamed file 'scripts/test/shots.py' => 'scripts/test/spheresFactory.py'
--- scripts/test/shots.py	2011-01-12 16:22:18 +0000
+++ scripts/test/spheresFactory.py	2011-01-26 12:39:07 +0000
@@ -1,5 +1,5 @@
 
-from yade import utils,log
+from yade import geom,utils,log
 
 shotsId,steelId=O.materials.append([
 	FrictMat(young=50e9,density=6000,poisson=.2,label='shots'),
@@ -24,10 +24,11 @@
 		)],
 		[Law2_ScGeom_MindlinPhys_Mindlin(label='contactLaw')]
 	),
-	# not used now
-	#GravityEngine(gravity=(0,0,0)), DragForceApplier(density=1,...)
 	NewtonIntegrator(damping=0),
-	NozzleFactory(maxParticles=10000,radius=8e-3,center=(0,-15e-3,15e-3),rMin=0.28e-3,rMax=0.29e-3,vMin=100,vMax=100,vAngle=0,massFlowRate=100./60,normal=(0,1.5,-1),label='factory',materialId=shotsId),
+	## CircularFactory: disk if length=0 or cylinder if length>0
+	CircularFactory(maxParticles=10000,radius=8e-3,length=16e-3,center=(0,-15e-3,15e-3),rMin=0.28e-3,rMax=0.29e-3,vMin=100,vMax=100,vAngle=0,massFlowRate=100./60,normal=(0,1.5,-1),label='factory',materialId=shotsId),
+	## QuadroFactory: a line, plane or cuboid  
+	#QuadroFactory(maxParticles=10000,extents=(8e-3,8e-3,8e-3),center=(0,-15e-3,15e-3),rMin=0.28e-3,rMax=0.29e-3,vMin=100,vMax=100,vAngle=0,massFlowRate=100./60,normal=(0,1.5,-1),label='factory',materialId=shotsId),
 	DomainLimiter(lo=(-30e-3,-30e-3,0),hi=(30e-3,30e-3,60e-3),iterPeriod=200),
 	#VTKRecorder(recorders=['spheres','facets','velocity'],fileName='/tmp/nozzle-',iterPeriod=500),
 
@@ -38,7 +39,7 @@
 # we cannot use utils.PWaveTimeStep directly, since there are no spheres generated yet
 O.dt=utils.SpherePWaveTimeStep(factory.rMin,O.materials[factory.materialId].density,O.materials[factory.materialId].young)
 O.saveTmp()
-O.timingEnabled=True
+#O.timingEnabled=True
 from yade import timing
 try:
 	from yade import qt