yade-dev team mailing list archive

Thread
Date

[Branch ~yade-dev/yade/trunk] Rev 1932: - merge stiffness/inertia control in PeriTriaxEngine and remove periEngine

To: Yade developers <yade-dev@xxxxxxxxxxxxxxxxxxx>
From: noreply@xxxxxxxxxxxxx
Date: Thu, 31 Dec 2009 15:25:21 -0000
Reply-to: noreply@xxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx

------------------------------------------------------------
revno: 1932
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Thu 2009-12-31 16:22:20 +0100
message:
  - merge stiffness/inertia control in PeriTriaxEngine and remove periEngine
  - show usage of inertia (commented out) in py script 
modified:
  pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp
  pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp
  scripts/test/periodic-triax.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/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp	2009-12-30 15:31:04 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp	2009-12-31 15:22:20 +0000
@@ -97,248 +97,174 @@
 
 void PeriTriaxController::strainStressStiffUpdate(){
 	// update strain first
-	const Vector3r& cellSize(scene->cell->getSize());
-	for(int i=0; i<3; i++) strain[i]=scene->cell->trsf[i][i]-1.;
-	// stress and stiffness
-	Vector3r sumForce(Vector3r::ZERO), sumStiff(Vector3r::ZERO), sumLength(Vector3r::ZERO);
+	//"Natural" strain, correct for large deformations, only used for comparison with goals
+	for ( int i=0; i<3; i++ ) strain[i]=Mathr::Log ( scene->cell->trsf[i][i] );
+
+	//stress tensor and stiffness
+	
+	//Compute volume of the deformed cell
+	// NOTE : needs refSize, could be generalized to arbitrary initial shapes using trsf*refHsize 
+	Real volume = ( scene->cell->trsf*Matrix3r ( scene->cell->refSize[0],0,0, 0,scene->cell->refSize[1],0,0,0,scene->cell->refSize[2] ) ).Determinant();
+
+	//Compute sum(fi*lj) and stiffness
+	stressTensor = Matrix3r::ZERO;
+	Vector3r sumStiff ( Vector3r::ZERO );
 	int n=0;
-	FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){
-		if(!I->isReal()) continue;
-		n++;
-		NormalShearInteraction* nsi=YADE_CAST<NormalShearInteraction*>(I->interactionPhysics.get());
-		GenericSpheresContact* gsc=YADE_CAST<GenericSpheresContact*>(I->interactionGeometry.get());
-		for(int i=0; i<3; i++){
-			const Real absNormI=abs(gsc->normal[i]);
-			Real f=(reversedForces?-1.:1.)*(nsi->normalForce[i]+nsi->shearForce[i]);
-			// force is as applied to id1, and normal is from id1 towards id2;
-			// therefore if they have the same sense, it is tensile force (positive)
-			sumForce[i]+=Mathr::Sign(f*gsc->normal[i])*abs(f);
-			sumLength[i]+=absNormI*(gsc->refR1+gsc->refR2);
-			sumStiff[i]+=absNormI*nsi->kn+(1-absNormI)*nsi->ks;
+	//NOTE : This sort of loops on interactions could be removed if we had callbacks in e.g. constitutive laws
+	FOREACH ( const shared_ptr<Interaction>&I, *scene->interactions )
+	{
+		if ( !I->isReal() ) continue;
+		NormalShearInteraction* nsi=YADE_CAST<NormalShearInteraction*> ( I->interactionPhysics.get() );
+		GenericSpheresContact* gsc=YADE_CAST<GenericSpheresContact*> ( I->interactionGeometry.get() );
+		//Contact force
+		Vector3r f= ( reversedForces?-1.:1. ) * ( nsi->normalForce+nsi->shearForce );
+		//branch vector, FIXME : the first definition generalizes to non-spherical bodies but needs wrapped coords.
+		//    Vector3r branch=
+		//   (reversedForces?-1.:1.)*(Body::byId(I->getId1())->state->pos-Body::byId(I->getId2())->state->pos);
+		Vector3r branch= ( reversedForces?-1.:1. ) *gsc->normal* ( gsc->refR1+gsc->refR2 );
+		// tensorial product f*branch
+		stressTensor+=Matrix3r ( f, branch );
+		if ( !dynCell )
+		{
+			for ( int i=0; i<3; i++ ) sumStiff[i]+=abs ( gsc->normal[i] ) *nsi->kn+ ( 1-abs ( gsc->normal[i] ) ) *nsi->ks;
+			n++;
 		}
 	}
-	if(n>0){
-		// sumForce/area is sum of stresses on fictious slices; the number of slices is (average length)/length in the respective sense
-		for(int i=0; i<3; i++) stress[i]=((1./n)*sumLength[i]/cellSize[i])*sumForce[i]/(cellSize[(i+1)%3]*cellSize[(i+2)%3]);
-		stiff=(1./n)*sumStiff;
-	} else { stiff=Vector3r::ZERO; stress=Vector3r::ZERO; }
+	//Compute stressTensor=sum(fi*lj)/Volume (Love equation)
+	stressTensor /= volume;
+	for ( int axis=0; axis<3; axis++ ) stress[axis]=stressTensor[axis][axis];
+	LOG_DEBUG ( "stressTensor : "<<endl
+				<<stressTensor[0][0]<<" "<<stressTensor[0][1]<<" "<<stressTensor[0][2]<<endl
+				<<stressTensor[1][0]<<" "<<stressTensor[1][1]<<" "<<stressTensor[1][2]<<endl
+				<<stressTensor[2][0]<<" "<<stressTensor[2][1]<<" "<<stressTensor[2][2]<<endl
+				<<"unbalanced = "<<Shop::unbalancedForce ( /*useMaxForce=*/false,scene ) );
+
+	if ( n>0 ) stiff= ( 1./n ) *sumStiff;
+	else stiff=Vector3r::ZERO;
 }
 
+
 CREATE_LOGGER(PeriTriaxController);
 
 
 
-void PeriTriaxController::action(Scene*){
-	if(!scene->isPeriodic){ throw runtime_error("PeriTriaxController run on aperiodic simulation."); }
+void PeriTriaxController::action ( Scene* )
+{
+	if ( !scene->isPeriodic ) { throw runtime_error ( "PeriTriaxController run on aperiodic simulation." ); }
 	const Vector3r& cellSize=scene->cell->getSize();
-	Vector3r cellArea=Vector3r(cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1]);
+	Vector3r cellArea=Vector3r ( cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1] );
 	// initial updates
 	const Vector3r& refSize=scene->cell->refSize;
-	if(maxBodySpan[0]<=0){
-		FOREACH(const shared_ptr<Body>& b,*scene->bodies){
-			if(!b || !b->bound) continue;
-			for(int i=0; i<3; i++) maxBodySpan[i]=max(maxBodySpan[i],b->bound->max[i]-b->bound->min[i]);
+	if ( maxBodySpan[0]<=0 )
+	{
+		FOREACH ( const shared_ptr<Body>& b,*scene->bodies )
+		{
+			if ( !b || !b->bound ) continue;
+			for ( int i=0; i<3; i++ ) maxBodySpan[i]=max ( maxBodySpan[i],b->bound->max[i]-b->bound->min[i] );
 		}
 	}
 	// check current size
-	if(2.1*maxBodySpan[0]>cellSize[0] || 2.1*maxBodySpan[1]>cellSize[1] || 2.1*maxBodySpan[2]>cellSize[2]){
-		LOG_DEBUG("cellSize="<<cellSize<<", maxBodySpan="<<maxBodySpan);
-		throw runtime_error("Minimum cell size is smaller than 2.1*maxBodySpan (periodic collider requirement)");
+	if ( 2.1*maxBodySpan[0]>cellSize[0] || 2.1*maxBodySpan[1]>cellSize[1] || 2.1*maxBodySpan[2]>cellSize[2] )
+	{
+		LOG_DEBUG ( "cellSize="<<cellSize<<", maxBodySpan="<<maxBodySpan );
+		throw runtime_error ( "Minimum cell size is smaller than 2.1*maxBodySpan (periodic collider requirement)" );
 	}
-	bool doUpdate((scene->currentIteration%globUpdate)==0);
-	if(doUpdate || min(stiff[0],min(stiff[1],stiff[2]))<=0){ strainStressStiffUpdate(); }
+	bool doUpdate ( ( scene->currentIteration%globUpdate ) ==0 );
+	if ( doUpdate || min ( stiff[0],min ( stiff[1],stiff[2] ) ) <=0 || dynCell ) { strainStressStiffUpdate(); }
 
-	bool allOk=true; Vector3r cellGrow(Vector3r::ZERO);
+	bool allOk=true; Vector3r cellGrow ( Vector3r::ZERO );
 	// apply condition along each axis separately (stress or strain)
-	for(int axis=0; axis<3; axis++){
-		Real maxGrow=refSize[axis]*maxStrainRate[axis]*scene->dt;
-		if(stressMask & (1<<axis)){ // control stress
-			// stiffness K=EA; Ïƒâ‚=goal stress; Î”Ïƒ wanted stress difference to apply
-			// Î”ÎµE=(Î”l/lâ‚€)(K/A) = Î”Ïƒ=Ïƒâ‚-Ïƒ â†” Î”l=(Ïƒâ‚-Ïƒ)lâ‚€(A/K)
-			cellGrow[axis]=(goal[axis]-stress[axis])*refSize[axis]*cellArea[axis]/(stiff[axis]>0?stiff[axis]:1.);
-			LOG_TRACE(axis<<": stress="<<stress[axis]<<", goal="<<goal[axis]<<", cellGrow="<<cellGrow[axis]);
-		} else { // control strain
-			cellGrow[axis]=(goal[axis]-strain[axis])*refSize[axis];
-			LOG_TRACE(axis<<": strain="<<strain[axis]<<", goal="<<goal[axis]<<", cellGrow="<<cellGrow[axis]);
+	for ( int axis=0; axis<3; axis++ )
+	{
+		//NOTE (1): don't multiply by refSize here (see note (2))
+		Real maxGrow=/*refSize[axis]**/maxStrainRate[axis]*scene->dt;
+		if ( stressMask & ( 1<<axis ) )   // control stress
+		{
+			if ( !dynCell )
+			{
+				// stiffness K=EA; Ïƒâ‚=goal stress; Î”Ïƒ wanted stress difference to apply
+				// Î”ÎµE=(Î”l/lâ‚€)(K/A) = Î”Ïƒ=Ïƒâ‚-Ïƒ â†” Î”l=(Ïƒâ‚-Ïƒ)lâ‚€(A/K)
+				cellGrow[axis]= ( goal[axis]-stress[axis] ) */*refSize[axis]**/cellArea[axis]/ ( stiff[axis]>0?stiff[axis]:1. );
+				LOG_TRACE ( axis<<": stress="<<stress[axis]<<", goal="<<goal[axis]<<", cellGrow="<<cellGrow[axis] );
+			}
+			else  //accelerate the deformation using the density of the period, includes Cundall's damping
+			{
+				assert ( mass>0 );//user set
+				Real dampFactor = 1 - growDamping*Mathr::Sign ( ( scene->cell->velGrad[axis][axis] ) * ( goal[axis]-stress[axis] ) );
+				scene->cell->velGrad[axis][axis]+=dampFactor*scene->dt* ( goal[axis]-stress[axis] ) /mass;
+				LOG_TRACE ( axis<<": stress="<<stress[axis]<<", goal="<<goal[axis]<<", velGrad="<<scene->cell->velGrad[axis][axis] );
+			}
+		}
+		else   // control strain
+		{
+			///NOTE : everything could be generalized to 9 independant components by comparing F[i,i] vs. Matrix3r goal[i,i], but it would be simpler in that case to let the user set the prescribed loading rates velGrad[i,i] when [i,i] is not stress-controlled. This "else" would disappear.
+			cellGrow[axis]=Mathr::Exp ( goal[axis]-strain[axis] ) /**refSize[axis]*/-1;
+			LOG_TRACE ( axis<<": strain="<<strain[axis]<<", goal="<<goal[axis]<<", cellGrow="<<cellGrow[axis] );
 		}
 		// limit maximum strain rate
-		if(abs(cellGrow[axis])>maxGrow) cellGrow[axis]=Mathr::Sign(cellGrow[axis])*maxGrow;
+		if ( abs ( cellGrow[axis] ) >maxGrow ) cellGrow[axis]=Mathr::Sign ( cellGrow[axis] ) *maxGrow;
 		// do not shrink below minimum cell size (periodic collider condition), although it is suboptimal WRT resulting stress
-		cellGrow[axis]=max(cellGrow[axis],-(cellSize[axis]-2.1*maxBodySpan[axis]));
+		/** NOTE : this is one place where size can't be removed (except for span tests). Commented out for now. Do we really need to handle this very special case here? The simulation will go crazy anyway, so we can as well throw at the next step. I'm trying to remove "size" when it's not needed, partly because I have in mind independant control on 9 components of stress/velGrad and it will be easier with a the matrix formalism (we dont have 9 sizes...). If one day size is not used anywhere, we can remove it. (B.C.)*/
+		//cellGrow[axis]=max(cellGrow[axis],-(cellSize[axis]-2.1*maxBodySpan[axis]));
+
 		// steady evolution with fluctuations; see TriaxialStressController
-		cellGrow[axis]=(1-growDamping)*cellGrow[axis]+.8*prevGrow[axis];
+		cellGrow[axis]= ( 1-growDamping ) *cellGrow[axis]+.8*prevGrow[axis];
 		// crude way of predicting stress, for steps when it is not computed from intrs
-		if(doUpdate) LOG_DEBUG(axis<<": cellGrow="<<cellGrow[axis]<<", new stress="<<stress[axis]<<", new strain="<<strain[axis]);
+		if ( doUpdate ) LOG_DEBUG ( axis<<": cellGrow="<<cellGrow[axis]<<", new stress="<<stress[axis]<<", new strain="<<strain[axis] );
 		// signal if condition not satisfied
-		if(stressMask & (1<<axis)){
+		if ( stressMask & ( 1<<axis ) )
+		{
 			Real curr=stress[axis];
-			if((goal[axis]!=0 && abs((curr-goal[axis])/goal[axis])>relStressTol) || abs(curr-goal[axis])>absStressTol) allOk=false;
-		} else {
+			if ( ( goal[axis]!=0 && abs ( ( curr-goal[axis] ) /goal[axis] ) >relStressTol ) || abs ( curr-goal[axis] ) >absStressTol ) allOk=false;
+		}
+		else
+		{
 			Real curr=strain[axis];
 			// since strain is prescribed exactly, tolerances need just to accomodate rounding issues
-			if((goal[axis]!=0 && abs((curr-goal[axis])/goal[axis])>1e-6) || abs(curr-goal[axis])>1e-6) {
+			if ( ( goal[axis]!=0 && abs ( ( curr-goal[axis] ) /goal[axis] ) >1e-6 ) || abs ( curr-goal[axis] ) >1e-6 )
+			{
 				allOk=false;
-				if(doUpdate) LOG_DEBUG("Strain not OK; "<<abs(curr-goal[axis])<<">1e-6");
+				if ( doUpdate ) LOG_DEBUG ( "Strain not OK; "<<abs ( curr-goal[axis] ) <<">1e-6" );
 			}
 		}
 	}
-	assert(scene->dt>0.);
+	assert ( scene->dt>0. );
 	// update stress and strain
-	for(int axis=0; axis<3; axis++){
+	if (!dynCell) for ( int axis=0; axis<3; axis++ )
+	{
 		// either prescribe velocity gradient
-		scene->cell->velGrad[axis][axis]=cellGrow[axis]/(scene->dt*refSize[axis]);
+		//NOTE (2) : ... and don't divide here
+		scene->cell->velGrad[axis][axis]=cellGrow[axis]/ ( scene->dt/**refSize[axis]*/ );
 		// or strain increment (but NOT both)
 		// strain[axis]+=cellGrow[axis]/refSize[axis];
 
 		// take in account something like poisson's effect hereâ€¦
 		//Real bogusPoisson=0.25; int ax1=(axis+1)%3,ax2=(axis+2)%3;
-		if(stiff[axis]>0) stress[axis]+=(cellGrow[axis]/refSize[axis])*(stiff[axis]/cellArea[axis]); //-bogusPoisson*(cellGrow[ax1]/refSize[ax1])*(stiff[ax1]/cellArea[ax1])-bogusPoisson*(cellGrow[ax2]/refSize[ax2])*(stiff[ax2]/cellArea[ax2]);
+		//don't modify stress if dynCell, testing only stiff[axis]>0 would not allow switching the control mode in simulations,
+		///NOTE : this one uses size and is hard to generalize to 9 components as well, how can we predict the full stress tensor? It would need a 9x9 stifness matrix!! At the end, I wonder if the implementation of independant control on the 9 components would not be easier and cleaner in a different engine... This stiffness approach is not a piece of cake. Controlling arbitrary components using inertia and comparing Fkk vs. goal[k][k] is really a different (and easier) problem. I feel like I'm mixing two different things, with a lot of complication to preserve a few specific features that will not be generalised at the end, just for the sake of avoiding a few duplicated lines. What do you think?
+		if ( stiff[axis]>0 && !dynCell ) stress[axis]+= ( cellGrow[axis]/refSize[axis] ) * ( stiff[axis]/cellArea[axis] ); //-bogusPoisson*(cellGrow[ax1]/refSize[ax1])*(stiff[ax1]/cellArea[ax1])-bogusPoisson*(cellGrow[ax2]/refSize[ax2])*(stiff[ax2]/cellArea[ax2]);
 	}
 	// change cell size now
 	// scene->cell->refSize+=cellGrow;
 
 	strainRate=cellGrow/scene->dt;
 
-	if(allOk){
-		if(doUpdate || currUnbalanced<0) {
-			currUnbalanced=Shop::unbalancedForce(/*useMaxForce=*/false,scene);
-			LOG_DEBUG("Stress/strain="<<(stressMask&1?stress[0]:strain[0])<<","<<(stressMask&2?stress[1]:strain[1])<<","<<(stressMask&4?stress[2]:strain[2])<<", goal="<<goal<<", unbalanced="<<currUnbalanced);
+	if ( allOk )
+	{
+		if ( doUpdate || currUnbalanced<0 )
+		{
+			currUnbalanced=Shop::unbalancedForce ( /*useMaxForce=*/false,scene );
+			LOG_DEBUG ( "Stress/strain="<< ( stressMask&1?stress[0]:strain[0] ) <<","<< ( stressMask&2?stress[1]:strain[1] ) <<","<< ( stressMask&4?stress[2]:strain[2] ) <<", goal="<<goal<<", unbalanced="<<currUnbalanced );
 		}
-		if(currUnbalanced<maxUnbalanced){
+		if ( currUnbalanced<maxUnbalanced )
+		{
 			// LOG_INFO("Goal reached, packing stable.");
-			if(!doneHook.empty()){
-				LOG_DEBUG("Running doneHook: "<<doneHook);
-				PyGILState_STATE gstate; gstate=PyGILState_Ensure(); PyRun_SimpleString(doneHook.c_str()); PyGILState_Release(gstate);
-			} else { Omega::instance().stopSimulationLoop(); }
-		}
-	}
-}
-
-CREATE_LOGGER(PeriController);
-
-void PeriController::action(Scene*){
-	if(!scene->isPeriodic){ throw runtime_error("PeriController run on aperiodic simulation."); }
-	
-	const Vector3r& cellSize=scene->cell->getSize();
-	const Vector3r& refSize=scene->cell->refSize;
-	//Matrix3r Hsize (refSize[0],0,0,0,refSize[1],0,0,0,refSize[2]);
-	//Hsize =  scene->cell->trsf*Hsize;
-	//Vector3r cellArea=Vector3r(cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1]);
-	// initial updates
-	LOG_DEBUG("Entering PeriController::action");
-	if(maxBodySpan[0]<=0){
-		FOREACH(const shared_ptr<Body>& b,*scene->bodies){
-			if(!b || !b->bound) continue;
-			for(int i=0; i<3; i++) maxBodySpan[i]=max(maxBodySpan[i],b->bound->max[i]-b->bound->min[i]);
-		}
-	}
-	// check current size
-	if(2.1*maxBodySpan[0]>cellSize[0] || 2.1*maxBodySpan[1]>cellSize[1] || 2.1*maxBodySpan[2]>cellSize[2]){
-		LOG_DEBUG("cellSize="<<cellSize<<", maxBodySpan="<<maxBodySpan);
-		throw runtime_error("Minimum cell size is smaller than 2.1*maxBodySpan (periodic collider requirement)");
-	}
-	
-	stressStrainUpdate();
-	
-	bool allOk=true; Vector3r cellGrow(Vector3r::ZERO);
-	// apply condition along each axis separately (stress or strain)
-	for(int axis=0; axis<3; axis++){
-		Real maxGrow=refSize[axis]*maxStrainRate[axis]*scene->dt;
-		if(stressMask & (1<<axis)){ // control stress
-			//accelerate the deformation using the density of the period, includes Cundall's damping
-			Real mass = 0.2;
-			Real dampFactor = 1 - growDamping*Mathr::Sign ( (scene->cell->velGrad[axis][axis])*(goal[axis]-stress[axis]) );
-			scene->cell->velGrad[axis][axis]+=dampFactor*scene->dt*(goal[axis]-stress[axis])/mass;
-			LOG_TRACE(axis<<": stress="<<stress[axis]<<", goal="<<goal[axis]<<", velGrad="<<scene->cell->velGrad[axis][axis]);
- 		} else { // don't control strainRate, strainRate should always be defined by the user
-// 			cellGrow[axis]=(goal[axis]-strain[axis])*refSize[axis];
- 			LOG_TRACE(axis<<": strain="<<strain[axis]<<", goal="<<goal[axis]<<", cellGrow="<<cellGrow[axis]);
-		}
-		// limit maximum strain rate
-		//if(abs(cellGrow[axis])>maxGrow) cellGrow[axis]=Mathr::Sign(cellGrow[axis])*maxGrow;
-		// do not shrink below minimum cell size (periodic collider condition), although it is suboptimal WRT resulting stress
-		//cellGrow[axis]=max(cellGrow[axis],-(cellSize[axis]-2.1*maxBodySpan[axis]));
-		// steady evolution with fluctuations; see TriaxialStressController
-		//cellGrow[axis]=(1-growDamping)*cellGrow[axis]+.8*prevGrow[axis];
-		// crude way of predicting stress, for steps when it is not computed from intrs
-		//LOG_DEBUG(axis<<": cellGrow="<<cellGrow[axis]<<", new stress="<<stress[axis]<<", new strain="<<strain[axis]);
-		// signal if condition not satisfied
-		if(stressMask & (1<<axis)){
-			Real curr=stress[axis];
-			if((goal[axis]!=0 && abs((curr-goal[axis])/goal[axis])>relStressTol) || abs(curr-goal[axis])>absStressTol) allOk=false;
-		} else {
-			Real curr=strain[axis];
-			// since strain is prescribed exactly, tolerances need just to accomodate rounding issues
-			if((goal[axis]!=0 && abs((curr-goal[axis])/goal[axis])>1e-6) || abs(curr-goal[axis])>1e-6) {
-				allOk=false;
-				LOG_DEBUG("Strain not OK; "<<abs(curr-goal[axis])<<">1e-6");
+			if ( !doneHook.empty() )
+			{
+				LOG_DEBUG ( "Running doneHook: "<<doneHook );
+				PyGILState_STATE gstate; gstate=PyGILState_Ensure(); PyRun_SimpleString ( doneHook.c_str() ); PyGILState_Release ( gstate );
 			}
-		}
-	}
-	
-	// update stress and strain
-// 	for(int axis=0; axis<3; axis++){
-// 		// either prescribe velocity gradient
-// 		scene->cell->velGrad[axis][axis]=cellGrow[axis]/(scene->dt*refSize[axis]);
-// 	}
-	
-// 	strainRate=cellGrow/scene->dt;
-	
-	if(allOk){
-		if(currUnbalanced<0 || (scene->currentIteration%globUpdate)==0 ) {
-			currUnbalanced=Shop::unbalancedForce(/*useMaxForce=*/false,scene);
-			LOG_DEBUG("Stress/strain="<<(stressMask&1?stress[0]:strain[0])<<","<<(stressMask&2?stress[1]:strain[1])<<","<<(stressMask&4?stress[2]:strain[2])<<", goal="<<goal<<", unbalanced="<<currUnbalanced);
-		}
-		if(currUnbalanced<maxUnbalanced){
-			// LOG_INFO("Goal reached, packing stable.");
-			if(!doneHook.empty()){
-				LOG_DEBUG("Running doneHook: "<<doneHook);
-				PyGILState_STATE gstate; gstate=PyGILState_Ensure(); PyRun_SimpleString(doneHook.c_str()); PyGILState_Release(gstate);
-			} else { Omega::instance().stopSimulationLoop(); }
-		}
-	}
-}
-
-Matrix3r PeriController::stressStrainUpdate(){
-	//Compute volume
-	Matrix3r Hsize (scene->cell->refSize[0],0,0,
-			0,scene->cell->refSize[1],0,
-  			 0,0,scene->cell->refSize[2]);
-	Hsize =  scene->cell->trsf*Hsize;
-	Real volume = Hsize.Determinant();
-	//"Natural" strain, correct for large deformations, only used for comparison with goals
-	for(int i=0; i<3; i++) strain[i]=Mathr::Log(scene->cell->trsf[i][i]);
-	
-	//Compute sum(fi*lj)
-	Matrix3r stressTensor(Matrix3r::ZERO);
-	
-	FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){
-		if(!I->isReal()) continue;
-		//n++;
-		NormalShearInteraction* nsi=YADE_CAST<NormalShearInteraction*>(I->interactionPhysics.get());
-		GenericSpheresContact* gsc=YADE_CAST<GenericSpheresContact*>(I->interactionGeometry.get());
-		for(int i=0; i<3; i++){
-			//Contact force
-			Vector3r f=(reversedForces?-1.:1.)*(nsi->normalForce+nsi->shearForce);
-			
-			//branch vector, FIXME : the first definition generalizes to non-spherical bodies but needs wrapped coords.
-// 			Vector3r branch=
-//			(reversedForces?-1.:1.)*(Body::byId(I->getId1())->state->pos-Body::byId(I->getId2())->state->pos);
-			Vector3r branch=gsc->normal*(gsc->refR1+gsc->refR2);
-
-			//Cross product
-			// tensor product f*branch
-			stressTensor+=Matrix3r (f, branch);
-			
-		}
-	}
-	//Compute stress=sum(fi*lj)/Volume (Love equation)
-	stressTensor /= volume;
-	for (int axis=0; axis<3; axis++) stress[axis]=stressTensor[axis][axis];
-	LOG_DEBUG("stressTensor : "<<endl
-			<<stressTensor[0][0]<<" "<<stressTensor[0][1]<<" "<<stressTensor[0][2]<<endl
-			<<stressTensor[1][0]<<" "<<stressTensor[1][1]<<" "<<stressTensor[1][2]<<endl
-			<<stressTensor[2][0]<<" "<<stressTensor[2][1]<<" "<<stressTensor[2][2]<<endl
-			<<"unbalanced = "<<Shop::unbalancedForce(/*useMaxForce=*/false,scene));
-	return stressTensor;
+			else { Omega::instance().stopSimulationLoop(); }
+		}
+	}
 }

=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp	2009-12-30 15:31:04 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp	2009-12-31 15:22:20 +0000
@@ -47,7 +47,10 @@
 	//! For broken constitutive laws, normalForce and shearForce on interactions are in the reverse sense
 	//! see https://bugs.launchpad.net/yade/+bug/493102
 	bool reversedForces;
-	//! Desired stress or strain values (depending on stressMask)
+	//! Imposed stress can be controlled using the packing stiffness or by applying the laws of dynamic (dynCell=true)
+	//! Don't forget to assign a mass to the cell (PeriTriaxController->mass)
+	bool dynCell;
+	//! Desired stress or strain values (depending on stressMask), strains defined as strain(i)=log(Fii)
 	Vector3r goal;
 	//! mask determining strain/stress (0/1) meaning for goal components
 	int stressMask;
@@ -68,8 +71,9 @@
 	//! maximum body dimension (set automatically)
 	Vector3r maxBodySpan;
 
-
-	//! average stresses, updated at every step (only every globUpdate steps recomputed from interactions)
+	//! average stresses, updated at every step (only every globUpdate steps recomputed from interactions if !dynCell)
+	Matrix3r stressTensor;
+	//! diagonal terms of the stress tensor
 	Vector3r stress;
 	//! cell strain, updated at every step
 	Vector3r strain;
@@ -81,65 +85,14 @@
 	Real currUnbalanced;
 	//! previous cell grow
 	Vector3r prevGrow;
+	//! mass of the cell (user set)
+	Real mass;
 
 	void action(Scene*);
 	void strainStressStiffUpdate();
-	PeriTriaxController(): reversedForces(false),goal(Vector3r::ZERO),stressMask(0),maxStrainRate(Vector3r(1,1,1)),maxUnbalanced(1e-4),absStressTol(1e3),relStressTol(3e-5),growDamping(.25),globUpdate(5),maxBodySpan(Vector3r(-1,-1,-1)),stress(Vector3r::ZERO),strain(Vector3r::ZERO),strainRate(Vector3r::ZERO),stiff(Vector3r::ZERO),currUnbalanced(-1),prevGrow(Vector3r::ZERO){}
-	REGISTER_ATTRIBUTES(GlobalEngine,(reversedForces)(goal)(stressMask)(maxStrainRate)(maxUnbalanced)(absStressTol)(relStressTol)(growDamping)(globUpdate)(doneHook)(stress)(strain)(strainRate)(stiff));
+	PeriTriaxController(): reversedForces(false),dynCell(false),goal(Vector3r::ZERO),stressMask(0),maxStrainRate(Vector3r(1,1,1)),maxUnbalanced(1e-4),absStressTol(1e3),relStressTol(3e-5),growDamping(.25),globUpdate(5),maxBodySpan(Vector3r(-1,-1,-1)),stressTensor(Matrix3r::ZERO),stress(Vector3r::ZERO),strain(Vector3r::ZERO),strainRate(Vector3r::ZERO),stiff(Vector3r::ZERO),currUnbalanced(-1),prevGrow(Vector3r::ZERO),mass(0){}
+	REGISTER_ATTRIBUTES(GlobalEngine,(reversedForces)(dynCell)(goal)(stressMask)(maxStrainRate)(maxUnbalanced)(absStressTol)(relStressTol)(growDamping)(globUpdate)(doneHook)(stressTensor)(stress)(strain)(strainRate)(stiff)(mass));
 	DECLARE_LOGGER;
 	REGISTER_CLASS_AND_BASE(PeriTriaxController,GlobalEngine);
 };
 REGISTER_SERIALIZABLE(PeriTriaxController);
-
-//TODO :(1) isotropic compression with mean stress/strain (preserve aspect ration like in other engines)
-//	(2) check vs. maxStrainRate? Not sure it is needed in that case
-class PeriController: public GlobalEngine{
-	public:
-	//! For broken constitutive laws, normalForce and shearForce on interactions are in the reverse sense
-	//! see https://bugs.launchpad.net/yade/+bug/493102
-		bool reversedForces;
-	//! Desired stress or strain values (depending on stressMask)
-		Vector3r goal;
-	//! mask determining strain/stress (0/1) meaning for goal components
-		int stressMask;
-	//! Maximum strain rate of the periodic cell
-		Vector3r maxStrainRate;
-	//! maximum unbalanced force (defaults to 1e-4)
-		Real maxUnbalanced;
-	//! Absolute stress tolerance (1e3)
-		Real absStressTol;
-	//! Relative stress tolerance (3e-5)
-		Real relStressTol;
-	//! Damping of cell resizing (0=perfect control, 1=no control at all); see also TriaxialStressController::wallDamping.
-		Real growDamping;
-	//! how often to recompute average stress, stiffness and unbalaced force (defaults to 100)
-		int globUpdate;
-	//! python command to be run when the desired state is reached
-		string doneHook;
-	//! maximum body dimension (set automatically)
-		Vector3r maxBodySpan;
-
-
-	//! average stresses, updated at every step (only every globUpdate steps recomputed from interactions)
-		Vector3r stress;
-	//! cell strain, updated at every step
-		Vector3r strain;
-	//! cell strain rate, updated at every step
-		Vector3r strainRate;
-	//! average stiffness, updated at every step (only every globUpdate steps recomputed from interactions)
-		Vector3r stiff;
-	//! current unbalanced force (updated every globUpdate)
-		Real currUnbalanced;
-	//! previous cell grow
-		Vector3r prevGrow;
-
-		void action(Scene*);
-		Matrix3r stressStrainUpdate();
-		PeriController(): reversedForces(false),goal(Vector3r::ZERO),stressMask(0),maxStrainRate(Vector3r(1,1,1)),maxUnbalanced(1e-4),absStressTol(1e3),relStressTol(3e-5),growDamping(.5),globUpdate(5),maxBodySpan(Vector3r(-1,-1,-1)),stress(Vector3r::ZERO),strain(Vector3r::ZERO),strainRate(Vector3r::ZERO),stiff(Vector3r::ZERO),currUnbalanced(-1),prevGrow(Vector3r::ZERO){}
-		REGISTER_ATTRIBUTES(GlobalEngine,(reversedForces)(goal)(stressMask)(maxStrainRate)(maxUnbalanced)(absStressTol)(relStressTol)(growDamping)(globUpdate)(doneHook)(stress)(strain)(strainRate)(stiff));
-		DECLARE_LOGGER;
-		REGISTER_CLASS_AND_BASE(PeriController,GlobalEngine);
-};
-REGISTER_SERIALIZABLE(PeriController);
-
-

=== modified file 'scripts/test/periodic-triax.py'
--- scripts/test/periodic-triax.py	2009-12-25 14:46:48 +0000
+++ scripts/test/periodic-triax.py	2009-12-31 15:22:20 +0000
@@ -3,9 +3,10 @@
 "Test and demonstrate use of PeriTriaxController."
 from yade import *
 from yade import pack,log,qt
-#log.setLevel('PeriTriaxController',log.DEBUG)
+log.setLevel('PeriTriaxController',log.TRACE)
 O.periodic=True
 O.cell.refSize=Vector3(.1,.1,.1)
+#O.cell.Hsize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
 sp=pack.SpherePack()
 radius=5e-3
 num=sp.makeCloud(Vector3().ZERO,O.cell.refSize,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
@@ -21,8 +22,10 @@
 		[SimpleElasticRelationships()],
 		[Law2_Dem3Dof_Elastic_Elastic()]
 	),
-	PeriTriaxController(goal=[-1e5,-1e5,0],stressMask=3,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
-	NewtonIntegrator(damping=.6),
+	#PeriTriaxController(maxUnbalanced=0.01,relStressTol=0.02,goal=[-1e4,-1e4,0],stressMask=3,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
+	#using cell inertia
+	PeriTriaxController(dynCell=True,mass=0.2,maxUnbalanced=0.01,relStressTol=0.02,goal=[-1e4,-1e4,0],stressMask=3,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
+	NewtonIntegrator(damping=.2),
 ]
 O.dt=utils.PWaveTimeStep()
 O.run();
@@ -34,7 +37,7 @@
 	if phase==0:
 		print 'Here we are: stress',triax['stress'],'strain',triax['strain'],'stiffness',triax['stiff']
 		print 'Now Îµz will go from 0 to .2 while Ïƒx and Ïƒy will be kept the same.'
-		triax['goal']=[-1e5,-1e5,.2]
+		triax['goal']=[-1e4,-1e4,-0.2]
 		phase+=1
 	elif phase==1:
 		print 'Here we are: stress',triax['stress'],'strain',triax['strain'],'stiffness',triax['stiff']

Follow ups

Re: [Branch ~yade-dev/yade/trunk] Rev 1932: - merge stiffness/inertia control in PeriTriaxEngine and remove periEngine
From: Václav Šmilauer, 2010-01-01