yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2537
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Thu 2010-11-11 15:22:33 +0100
message:
  1. Fix particle highlighting bug introduced lastly
  2. Enhance the SpherePack.makeClumpCloud routine to be much more efficient (uses bounding spheres)
  3. Add Ig2_Wall_Sphere_ScGeom
  4. Enhance scripts/test/clumpPack.py
modified:
  pkg/common/OpenGLRenderer.cpp
  pkg/common/Wall.cpp
  pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp
  pkg/dem/Ig2_Facet_Sphere_ScGeom.hpp
  pkg/dem/SpherePack.cpp
  pkg/dem/SpherePack.hpp
  scripts/test/clumpPack.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/common/OpenGLRenderer.cpp'
--- pkg/common/OpenGLRenderer.cpp	2010-11-09 08:27:14 +0000
+++ pkg/common/OpenGLRenderer.cpp	2010-11-11 14:22:33 +0000
@@ -323,7 +323,7 @@
 
 		// ignored in non-selection mode, use it always
 		glPushName(b->id);
-		bool highlight=(b->id==selId || b->clumpId==selId || b->shape->highlight);
+		bool highlight=(b->id==selId || (b->clumpId>=0 && b->clumpId==selId) || b->shape->highlight);
 
 		glPushMatrix();
 			AngleAxisr aa(ori);	
@@ -332,7 +332,7 @@
 			if(highlight){
 				// set hightlight
 				// different color for body highlighted by selection and by the shape attribute
-				const Vector3r& h((selId==b->id||selId==b->clumpId) ? highlightEmission0 : highlightEmission1);
+				const Vector3r& h((selId==b->id||(b->clumpId>=0 && selId==b->clumpId)) ? highlightEmission0 : highlightEmission1);
 				glMaterialv(GL_FRONT_AND_BACK,GL_EMISSION,h);
 				glMaterialv(GL_FRONT_AND_BACK,GL_SPECULAR,h);
 				shapeDispatcher(b->shape,b->state,wire || b->shape->wire,viewInfo);

=== modified file 'pkg/common/Wall.cpp'
--- pkg/common/Wall.cpp	2010-10-13 16:23:08 +0000
+++ pkg/common/Wall.cpp	2010-11-11 14:22:33 +0000
@@ -15,7 +15,7 @@
 	Wall* wall=static_cast<Wall*>(cm.get());
 	if(!bv){ bv=shared_ptr<Bound>(new Aabb); }
 	Aabb* aabb=static_cast<Aabb*>(bv.get());
-	if(scene->isPeriodic && scene->cell->hasShear()) throw logic_error(__FILE__ "Walls not (yet?) supported in sheared cell.");
+	if(scene->isPeriodic && scene->cell->hasShear()) throw logic_error(__FILE__ "Walls not supported in sheared cell.");
 	const Real& inf=std::numeric_limits<Real>::infinity();
 	aabb->min=Vector3r(-inf,-inf,-inf); aabb->min[wall->axis]=se3.position[wall->axis];
 	aabb->max=Vector3r( inf, inf, inf); aabb->max[wall->axis]=se3.position[wall->axis];

=== modified file 'pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp'
--- pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp	2010-10-27 18:54:33 +0000
+++ pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp	2010-11-11 14:22:33 +0000
@@ -10,9 +10,12 @@
 #include<yade/pkg-dem/ScGeom.hpp>
 #include<yade/pkg-common/Sphere.hpp>
 #include<yade/pkg-common/Facet.hpp>
+#include<yade/pkg-common/Wall.hpp>
 #include<yade/core/Scene.hpp>
 #include<yade/lib-base/Math.hpp>
 
+YADE_PLUGIN((Ig2_Facet_Sphere_ScGeom)(Ig2_Wall_Sphere_ScGeom));
+
 CREATE_LOGGER(Ig2_Facet_Sphere_ScGeom);
 
 bool Ig2_Facet_Sphere_ScGeom::go(const shared_ptr<Shape>& cm1,
@@ -127,7 +130,32 @@
 	return go(cm2,cm1,state2,state1,-shift2,force,c);
 }
 
-YADE_PLUGIN((Ig2_Facet_Sphere_ScGeom));
-
-//YADE_REQUIRE_FEATURE(PHYSPAR);
+
+/********* Wall + Sphere **********/
+
+bool Ig2_Wall_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){
+	Wall* wall=static_cast<Wall*>(cm1.get());
+	const Real radius=static_cast<Sphere*>(cm2.get())->radius;
+	const int& ax(wall->axis);
+	Real dist=(state2.pos)[ax]+shift2[ax]-state1.pos[ax]; // signed "distance" between centers
+	if(!c->isReal() && abs(dist)>radius && !force) { return false; }// wall and sphere too far from each other
+
+	// contact point is sphere center projected onto the wall
+	Vector3r contPt=state2.pos+shift2; contPt[ax]=state1.pos[ax];
+	Vector3r normal(0.,0.,0.);
+	// wall interacting from both sides: normal depends on sphere's position
+	assert(wall->sense==-1 || wall->sense==0 || wall->sense==1);
+	if(wall->sense==0) normal[ax]=dist>0?1.:-1.;
+	else normal[ax]=wall->sense==1?1.:-1;
+
+	bool isNew=!c->geom;
+	if(isNew) c->geom=shared_ptr<ScGeom>(new ScGeom());
+	const shared_ptr<ScGeom>& ws=YADE_PTR_CAST<ScGeom>(c->geom);
+	ws->radius1=ws->radius2=radius; // do the same as for facet-sphere: wall's "radius" is the same as the sphere's radius
+	ws->contactPoint=contPt;
+	ws->penetrationDepth=-(abs(dist)-radius);
+	// ws->normal is assigned by precompute
+	ws->precompute(state1,state2,scene,c,normal,isNew,shift2,noRatch);
+	return true;
+}
 

=== modified file 'pkg/dem/Ig2_Facet_Sphere_ScGeom.hpp'
--- pkg/dem/Ig2_Facet_Sphere_ScGeom.hpp	2010-10-13 16:23:08 +0000
+++ pkg/dem/Ig2_Facet_Sphere_ScGeom.hpp	2010-11-11 14:22:33 +0000
@@ -39,3 +39,16 @@
 REGISTER_SERIALIZABLE(Ig2_Facet_Sphere_ScGeom);
 
 
+class Ig2_Wall_Sphere_ScGeom: public IGeomFunctor{
+	public:
+		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);
+	YADE_CLASS_BASE_DOC_ATTRS(Ig2_Wall_Sphere_ScGeom,IGeomFunctor,"Create/update a :yref:`ScGeom` instance representing intersection of :yref:`Wall` and :yref:`Sphere`.",
+		((bool,noRatch,true,,"Avoid granular ratcheting"))
+	);
+	FUNCTOR2D(Wall,Sphere);
+	DEFINE_FUNCTOR_ORDER_2D(Wall,Sphere);
+};
+REGISTER_SERIALIZABLE(Ig2_Wall_Sphere_ScGeom);
+
+
+

=== modified file 'pkg/dem/SpherePack.cpp'
--- pkg/dem/SpherePack.cpp	2010-11-09 08:27:14 +0000
+++ pkg/dem/SpherePack.cpp	2010-11-11 14:22:33 +0000
@@ -277,18 +277,23 @@
 long SpherePack::makeClumpCloud(const Vector3r& mn, const Vector3r& mx, const vector<shared_ptr<SpherePack> >& _clumps, bool periodic, int num){
 	// recenter given clumps and compute their margins
 	vector<SpherePack> clumps; /* vector<Vector3r> margins; */ Vector3r boxMargins(Vector3r::Zero()); Real maxR=0;
+	vector<Real> boundRad; // squared radii of bounding sphere for each clump
 	FOREACH(const shared_ptr<SpherePack>& c, _clumps){
 		SpherePack c2(*c); 
 		c2.translate(c2.midPt()); //recenter
 		clumps.push_back(c2);
+		Real rSq=0;
+		FOREACH(const Sph& s, c2.pack) rSq=max(rSq,s.c.squaredNorm()+pow(s.r,2));
+		boundRad.push_back(sqrt(rSq));
 		Vector3r cMn,cMx; c2.aabb(cMn,cMx); // centered at zero now, this gives margin
 		//margins.push_back(periodic?cMx:Vector3r::Zero()); 
 		//boxMargins=boxMargins.cwise().max(cMx);
 		FOREACH(const Sph& s, c2.pack) maxR=max(maxR,s.r); // keep track of maximum sphere radius
 	}
+	std::list<ClumpInfo> clumpInfos;
 	Vector3r size=mx-mn;
 	if(periodic)(cellSize=size);
-	const int maxTry=1000;
+	const int maxTry=200;
 	int nGen=0; // number of clumps generated
 	// random point coordinate generator, with non-zero margins if aperiodic
 	static boost::minstd_rand randGen(TimingInfo::getNow(true));
@@ -306,8 +311,36 @@
 			Quaternionr ori(rndUnit(),rndUnit(),rndUnit(),rndUnit()); ori.normalize();
 			// copy the packing and rotate
 			SpherePack C(clumps[clumpChoice]); C.rotateAroundOrigin(ori); C.translate(pos);
+			const Real& rad(boundRad[clumpChoice]);
+			ClumpInfo ci; // to be used later, but must be here because of goto's
 
 			// find collisions
+			// check against bounding spheres of other clumps, and only check individual spheres if there is overlap
+			if(!periodic){
+				// check overlap with box margins first
+				if((pos+rad*Vector3r::Ones()).cwise().max(mx)!=mx || (pos-rad*Vector3r::Ones()).cwise().min(mn)!=mn){ FOREACH(const Sph& s, C.pack) if((s.c+s.r*Vector3r::Ones()).cwise().max(mx)!=mx || (s.c-s.r*Vector3r::Ones()).cwise().min(mn)!=mn) goto overlap; }
+				// check overlaps with bounding spheres of other clumps
+				FOREACH(const ClumpInfo& cInfo, clumpInfos){
+					bool detailedCheck=false;
+					// check overlaps between individual spheres and bounding sphere of the other clump
+					if((pos-cInfo.center).squaredNorm()<pow(rad+cInfo.rad,2)){ FOREACH(const Sph& s, C.pack) if(pow(s.r+cInfo.rad,2)>(s.c-cInfo.center).squaredNorm()){ detailedCheck=true; break; }}
+					// check sphere-by-sphere, since bounding spheres did overlap
+					if(detailedCheck){ FOREACH(const Sph& s, C.pack) for(int id=cInfo.minId; id<=cInfo.maxId; id++) if((s.c-pack[id].c).squaredNorm()<pow(s.r+pack[id].r,2)) goto overlap;}
+				}
+			} else {
+				FOREACH(const ClumpInfo& cInfo, clumpInfos){
+					// bounding spheres overlap (in the periodic space)
+					if(periPtDistSq(pos,cInfo.center)<pow(rad+cInfo.rad,2)){
+						bool detailedCheck=false;
+						// check spheres with bounding sphere of the other clump
+						FOREACH(const Sph& s, C.pack) if(pow(s.r+cInfo.rad,2)>periPtDistSq(s.c,cInfo.center)){ detailedCheck=true; break; }
+						// check sphere-by-sphere
+						if(detailedCheck){ FOREACH(const Sph& s, C.pack) for(int id=cInfo.minId; id<=cInfo.maxId; id++) if(periPtDistSq(s.c,pack[id].c)<pow(s.r+pack[id].r,2)) goto overlap; }
+					}
+				}
+			}
+
+			#if 0
 			// crude algorithm: check all spheres against all other spheres (slow!!)
 			// use vtkPointLocator, add all existing points and check distance of r+maxRadius, then refine
 			// for periodicity, duplicate points close than boxMargins to the respective boundary
@@ -332,9 +365,14 @@
 					}
 				}
 			}
+			#endif
 
 			// add the clump, if no collisions
-			FOREACH(const Sph& s, C.pack){ pack.push_back(Sph(s.c,s.r,/*number clumps consecutively*/nGen)); }
+			/*number clumps consecutively*/ ci.clumpId=nGen; ci.center=pos; ci.rad=rad; ci.minId=pack.size(); ci.maxId=pack.size()+C.pack.size(); 
+			FOREACH(const Sph& s, C.pack){
+				pack.push_back(Sph(s.c,s.r,ci.clumpId));
+			}
+			clumpInfos.push_back(ci);
 			nGen++;
 			//cerr<<"O";
 			break; // break away from the try-loop

=== modified file 'pkg/dem/SpherePack.hpp'
--- pkg/dem/SpherePack.hpp	2010-11-09 08:27:14 +0000
+++ pkg/dem/SpherePack.hpp	2010-11-11 14:22:33 +0000
@@ -27,6 +27,13 @@
 		Real xNorm=(x-x0)/(x1-x0);
 		return (xNorm-floor(xNorm))*(x1-x0);
 	}
+	Real periPtDistSq(const Vector3r& p1, const Vector3r& p2){
+		Vector3r dr;
+		for(int ax=0; ax<3; ax++) dr[ax]=min(cellWrapRel(p1[ax],p2[ax],p2[ax]+cellSize[ax]),cellWrapRel(p2[ax],p1[ax],p1[ax]+cellSize[ax]));
+		return dr.squaredNorm();
+	}
+	struct ClumpInfo{ int clumpId; Vector3r center; Real rad; int minId, maxId; };
+
 public:
 	enum {RDIST_RMEAN, RDIST_POROSITY, RDIST_PSD};
 	struct Sph{

=== modified file 'scripts/test/clumpPack.py'
--- scripts/test/clumpPack.py	2010-11-09 08:27:14 +0000
+++ scripts/test/clumpPack.py	2010-11-11 14:22:33 +0000
@@ -3,5 +3,22 @@
 c1=pack.SpherePack([((0,0,0),.5),((.5,0,0),.5),((0,.5,0),.3)])
 c2=pack.SpherePack([((0,0,0),.5),((.7,0,0),.3),((.9,0,0),.2)])
 sp=pack.SpherePack()
-sp.makeClumpCloud((0,0,0),(10,10,10),[c1,c2],periodic=False)
+print 'Generated # of clumps:',sp.makeClumpCloud((0,0,0),(15,15,15),[c1,c2],periodic=False)
 sp.toSimulation()
+
+O.bodies.append(utils.wall(position=0,axis=2))
+
+O.engines=[
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
+	InteractionLoop(
+		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
+		[Ip2_FrictMat_FrictMat_FrictPhys()],
+		[Law2_ScGeom_FrictPhys_CundallStrack()]
+	),
+	GravityEngine(gravity=(0,0,-100)),
+	NewtonIntegrator(damping=.4)
+]
+O.dt=.7*utils.PWaveTimeStep()
+O.saveTmp()
+