yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2801
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Sun 2011-04-03 18:49:58 +0200
message:
  -Update clumps forces and torques in Newton, as in https://blueprints.launchpad.net/yade/+spec/clump-forces.
  -Update Shop::unbalancedForce accordingly
  -Remove the old unbalancedForce code from TSController and make it point to the Shop version.
modified:
  core/ForceContainer.hpp
  pkg/dem/NewtonIntegrator.cpp
  pkg/dem/Shop.cpp
  pkg/dem/TriaxialStressController.cpp


--
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 'core/ForceContainer.hpp'
--- core/ForceContainer.hpp	2011-03-17 13:47:46 +0000
+++ core/ForceContainer.hpp	2011-04-03 16:49:58 +0000
@@ -66,13 +66,7 @@
 
 		inline void ensureSynced(){ if(!synced) throw runtime_error("ForceContainer not thread-synchronized; call sync() first!"); }
 
-		#if 0
-			/*! Function to allow friend classes to get force even if not synced.
-			* Dangerous! The caller must know what it is doing! (i.e. don't read after write
-			* for a particular body id. */
-			const Vector3r& getForceUnsynced (Body::id_t id){ensureSize(id); return _force[id];}
-			const Vector3r& getTorqueUnsynced(Body::id_t id){ensureSize(id); return _force[id];}
-		#endif
+		
 		// dummy function to avoid template resolution failure
 		friend class boost::serialization::access; template<class ArchiveT> void serialize(ArchiveT & ar, unsigned int version){}
 	public:
@@ -84,7 +78,6 @@
 				sizeOfThreads.push_back(0);
 			}
 		}
-
 		const Vector3r& getForce(Body::id_t id)         { ensureSynced(); return ((size_t)id<size)?_force[id]:_zero; }
 		void  addForce(Body::id_t id, const Vector3r& f){ ensureSize(id,omp_get_thread_num()); synced=false;   _forceData[omp_get_thread_num()][id]+=f;}
 		const Vector3r& getTorque(Body::id_t id)        { ensureSynced(); return ((size_t)id<size)?_torque[id]:_zero; }
@@ -93,6 +86,13 @@
 		void  addMove(Body::id_t id, const Vector3r& m) { ensureSize(id,omp_get_thread_num()); synced=false; moveRotUsed=true; _moveData[omp_get_thread_num()][id]+=m;}
 		const Vector3r& getRot(Body::id_t id)           { ensureSynced(); return ((size_t)id<size)?_rot[id]:_zero; }
 		void  addRot(Body::id_t id, const Vector3r& r)  { ensureSize(id,omp_get_thread_num()); synced=false; moveRotUsed=true; _rotData[omp_get_thread_num()][id]+=r;}
+		
+		/*! Function to allow friend classes to get force even if not synced. Used for clumps by NewtonIntegrator.
+		* Dangerous! The caller must know what it is doing! (i.e. don't read after write
+		* for a particular body id. */
+		Vector3r& getForceUnsynced (Body::id_t id){return ((size_t)id<size)?_force[id]:_zero;}
+		Vector3r& getTorqueUnsynced(Body::id_t id){return ((size_t)id<size)?_torque[id]:_zero;}
+		
 		/* To be benchmarked: sum thread data in getForce/getTorque upon request for each body individually instead of by the sync() function globally */
 		// this function is used from python so that running simulation is not slowed down by sync'ing on occasions
 		// since Vector3r writes are not atomic, it might (rarely) return wrong value, if the computation is running meanwhile

=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp	2011-01-29 22:47:18 +0000
+++ pkg/dem/NewtonIntegrator.cpp	2011-04-03 16:49:58 +0000
@@ -107,7 +107,21 @@
 			if(unlikely(!b || b->isClumpMember())) continue;
 
 			State* state=b->state.get(); const Body::id_t& id=b->getId();
-			Vector3r f=scene->forces.getForce(id), m=scene->forces.getTorque(id);
+			Vector3r f=Vector3r::Zero(), m=Vector3r::Zero();
+			// clumps forces
+			if(b->isClump()) {
+				b->shape->cast<Clump>().addForceTorqueFromMembers(state,scene,f,m);
+				#ifdef YADE_OPENMP
+				//it is safe here, since onky one thread will read/write
+				scene->forces.getTorqueUnsynced(id)+=m;
+				scene->forces.getForceUnsynced(id)+=f;
+				#else
+				scene->forces.addTorque(id,m);
+				scene->forces.addForce(id, f);
+				#endif
+			}
+			//in most cases, the initial force on clumps will zero and next line is not changing f and m, but make sure we don't miss something (e.g. user defined forces on clumps)
+			f=scene->forces.getForce(id); m=scene->forces.getTorque(id);
 			#ifdef YADE_DEBUG
 				if(isnan(f[0])||isnan(f[1])||isnan(f[2])) throw runtime_error(("NewtonIntegrator: NaN force acting on #"+lexical_cast<string>(id)+".").c_str());
 				if(isnan(m[0])||isnan(m[1])||isnan(m[2])) throw runtime_error(("NewtonIntegrator: NaN torque acting on #"+lexical_cast<string>(id)+".").c_str());
@@ -119,6 +133,8 @@
 			// in aperiodic boundaries, it is equal to absolute velocity
 			Vector3r fluctVel=isPeriodic?scene->cell->bodyFluctuationVel(b->state->pos,b->state->vel):state->vel;
 
+
+
 			// numerical damping & kinetic energy
 			if(unlikely(trackEnergy)) updateEnergy(b,state,fluctVel,f,m);
 
@@ -127,8 +143,6 @@
 
 			// for particles not totally blocked, compute accelerations; otherwise, the computations would be useless
 			if (state->blockedDOFs!=State::DOF_ALL) {
-				// forces
-				if(b->isClump()) b->shape->cast<Clump>().addForceTorqueFromMembers(state,scene,f,m);
 				// linear acceleration
 				Vector3r linAccel=computeAccel(f,state->mass,state->blockedDOFs);
 				cundallDamp2nd(dt,f,fluctVel,linAccel);

=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp	2011-03-16 21:02:40 +0000
+++ pkg/dem/Shop.cpp	2011-04-03 16:49:58 +0000
@@ -162,8 +162,8 @@
 	Real sumF=0,maxF=0,currF; int nb=0;
 	FOREACH(const shared_ptr<Body>& b, *rb->bodies){
 		if(!b || b->isClumpMember() || !b->isDynamic()) continue;
-		if(!b->isClump()){ currF=rb->forces.getForce(b->id).norm(); }
-		else { // for clump, sum forces from members on the clump itself
+		currF=rb->forces.getForce(b->id).norm();
+		if(b->isClump() && currF==0){ // this should not happen unless the function is called by an engine whose position in the loop is before Newton (with the exception of bodies which really have null force), because clumps forces are updated in Newton. Typical triaxial loops are using such ordering unfortunately (triaxEngine before Newton). So, here we make sure that they will get correct unbalance. In the future, it is better for optimality to check unbalancedF inside scripts at the end of loops, so that this "if" is never active.
 			Vector3r f(rb->forces.getForce(b->id)),m(Vector3r::Zero());
 			b->shape->cast<Clump>().addForceTorqueFromMembers(b->state.get(),rb,f,m);
 			currF=f.norm();

=== modified file 'pkg/dem/TriaxialStressController.cpp'
--- pkg/dem/TriaxialStressController.cpp	2010-12-31 14:35:21 +0000
+++ pkg/dem/TriaxialStressController.cpp	2011-04-03 16:49:58 +0000
@@ -14,6 +14,7 @@
 #include<yade/core/State.hpp>
 #include<assert.h>
 #include<yade/core/Scene.hpp>
+#include<yade/pkg/dem/Shop.hpp>
 
 CREATE_LOGGER(TriaxialStressController);
 YADE_PLUGIN((TriaxialStressController));
@@ -244,53 +245,6 @@
 /*!
     \fn TriaxialStressController::ComputeUnbalancedForce( bool maxUnbalanced)
  */
-Real TriaxialStressController::ComputeUnbalancedForce( bool maxUnbalanced)
-{
-	scene->forces.sync();
-	//compute the mean contact force
-	Real MeanForce = 0.f;
-	long nForce = 0;
-	shared_ptr<BodyContainer>& bodies = scene->bodies;
-
-	InteractionContainer::iterator ii    = scene->interactions->begin();
-	InteractionContainer::iterator iiEnd = scene->interactions->end();
-	for(  ; ii!=iiEnd ; ++ii ) {
-		if ((*ii)->isReal()) {
-			const shared_ptr<Interaction>& contact = *ii;
-			Real f = (static_cast<FrictPhys*> ((contact->phys.get()))->normalForce+static_cast<FrictPhys*>(contact->phys.get())->shearForce).squaredNorm();
-			if (f!=0) {
-			MeanForce += sqrt(f);
-			++nForce;}
-		}
-	}
-	if (nForce!=0) MeanForce /= nForce;
-
-	if (!maxUnbalanced) {
-		//compute mean Unbalanced Force
-		Real MeanUnbalanced=0;
-		long nBodies = 0;
-		BodyContainer::iterator bi    = bodies->begin();
-		BodyContainer::iterator biEnd = bodies->end();
-		Real f;
-		for(  ; bi!=biEnd ; ++bi ) {
-			if ((*bi)->isDynamic()) {
-				f=getForce(scene,(*bi)->getId()).norm();
-				MeanUnbalanced += f;
-				if (f!=0) ++nBodies;
-			}
-		}
-		if (nBodies != 0 && MeanForce != 0) MeanUnbalanced = MeanUnbalanced/nBodies/MeanForce;
-		return  MeanUnbalanced;
-	} else {
-		//compute max Unbalanced Force
-		Real MaxUnbalanced=0;
-		BodyContainer::iterator bi    = bodies->begin();
-		BodyContainer::iterator biEnd = bodies->end();
-		for(  ; bi!=biEnd ; ++bi ) if ((*bi)->isDynamic())
-				MaxUnbalanced = std::max(getForce(scene,(*bi)->getId()).norm(),MaxUnbalanced);
-		if (MeanForce != 0) MaxUnbalanced = MaxUnbalanced/MeanForce;
-		return MaxUnbalanced;
-	}
-}
+Real TriaxialStressController::ComputeUnbalancedForce( bool maxUnbalanced) {return Shop::unbalancedForce(maxUnbalanced,scene);}