yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2589
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Mon 2010-12-06 12:54:38 +0100
message:
  - PeriTriax : remove the constant-distance approximation in stress definition, using wrapped coordinates instead of radii.
  - Get back standalone energy tracing in ECL. 
modified:
  pkg/dem/ElasticContactLaw.cpp
  pkg/dem/ElasticContactLaw.hpp
  pkg/dem/PeriIsoCompressor.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 'pkg/dem/ElasticContactLaw.cpp'
--- pkg/dem/ElasticContactLaw.cpp	2010-11-22 17:58:56 +0000
+++ pkg/dem/ElasticContactLaw.cpp	2010-12-06 11:54:38 +0000
@@ -14,7 +14,7 @@
 
 YADE_PLUGIN((Law2_ScGeom_FrictPhys_CundallStrack)(Law2_Dem3DofGeom_FrictPhys_CundallStrack)(ElasticContactLaw));
 
-#if 0
+#if 1
 Real Law2_ScGeom_FrictPhys_CundallStrack::getPlasticDissipation() {return (Real) plasticDissipation;}
 void Law2_ScGeom_FrictPhys_CundallStrack::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;}
 Real Law2_ScGeom_FrictPhys_CundallStrack::elasticEnergy()
@@ -39,7 +39,7 @@
 		if(!I->isReal()) continue;
 		#ifdef YADE_DEBUG
 			// these checks would be redundant in the functor (LawDispatcher does that already)
-			if(!dynamic_cast<ScGeom*>(I->geom.get()) || !dynamic_cast<FrictPhys*>(I->phys.get())) continue;	
+			if(!dynamic_cast<ScGeom*>(I->geom.get()) || !dynamic_cast<FrictPhys*>(I->phys.get())) continue;
 		#endif
 			functor->go(I->geom, I->phys, I.get());
 	}
@@ -68,7 +68,7 @@
 	shearForce -= phys->ks*shearDisp;
 	Real maxFs = phys->normalForce.squaredNorm()*std::pow(phys->tangensOfFrictionAngle,2);
 
-	if (likely(!scene->trackEnergy)){//Update force but don't compute energy terms (see below))
+	if (likely(!scene->trackEnergy  && !traceEnergy)){//Update force but don't compute energy terms (see below))
 		// PFC3d SlipModel, is using friction angle. CoulombCriterion
 		if( shearForce.squaredNorm() > maxFs ){
 			Real ratio = sqrt(maxFs) / shearForce.norm();
@@ -81,7 +81,8 @@
 			//define the plastic work input and increment the total plastic energy dissipated
 			shearForce *= ratio;
 			Real dissip=((1/phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
-			if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,/*reset*/false);
+			if (traceEnergy) plasticDissipation += dissip;
+			else if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,/*reset*/false);
 		}
 		// compute elastic energy as well
 		scene->energy->add(0.5*(phys->normalForce.squaredNorm()/phys->kn+phys->shearForce.squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,/*reset at every timestep*/true);

=== modified file 'pkg/dem/ElasticContactLaw.hpp'
--- pkg/dem/ElasticContactLaw.hpp	2010-11-12 18:44:15 +0000
+++ pkg/dem/ElasticContactLaw.hpp	2010-12-06 11:54:38 +0000
@@ -16,20 +16,20 @@
 
 class Law2_ScGeom_FrictPhys_CundallStrack: public LawFunctor{
 	public:
-		//OpenMPAccumulator<Real> plasticDissipation;
+		OpenMPAccumulator<Real> plasticDissipation;
 		virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
-		//Real elasticEnergy ();
-		//Real getPlasticDissipation();
-		//void initPlasticDissipation(Real initVal=0);
+		Real elasticEnergy ();
+		Real getPlasticDissipation();
+		void initPlasticDissipation(Real initVal=0);
 		YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_CundallStrack,LawFunctor,"Law for linear compression, and Mohr-Coulomb plasticity surface without cohesion.\nThis law implements the classical linear elastic-plastic law from [CundallStrack1979]_ (see also [Pfc3dManual30]_). The normal force is (with the convention of positive tensile forces) $F_n=\\min(k_n u_n, 0)$. The shear force is $F_s=k_s u_s$, the plasticity condition defines the maximum value of the shear force : $F_s^{\\max}=F_n\\tan(\\phi)$, with $\\phi$ the friction angle.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_CundallStrack`, which uses :yref:`Dem3DofGeom` (sphere-box interactions are not implemented for the latest).\n\n.. note::\n This law is well tested in the context of triaxial simulation, and has been used for a number of published results (see e.g. [Scholtes2009b]_ and other papers from the same authors). It is generalised by :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`, which adds cohesion and moments at contact.",
 		((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
-		((bool,traceEnergy,false,Attr::hidden,"Define the total energy dissipated in plastic slips at all contacts."))
+		((bool,traceEnergy,false,,"Define the total energy dissipated in plastic slips at all contacts. This will trace only plastic energy in this law, see O.trackEnergy for a more complete energies tracing"))
 		((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
 		((int,elastPotentialIx,-1,(Attr::hidden|Attr::noSave),"Index for elastic potential energy (with O.trackEnergy)"))
 		,,
-		//.def("elasticEnergy",&Law2_ScGeom_FrictPhys_CundallStrack::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
-		//.def("plasticDissipation",&Law2_ScGeom_FrictPhys_CundallStrack::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_CundallStrack::traceEnergy` is true.")
-		//.def("initPlasticDissipation",&Law2_ScGeom_FrictPhys_CundallStrack::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).")
+		.def("elasticEnergy",&Law2_ScGeom_FrictPhys_CundallStrack::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
+		.def("plasticDissipation",&Law2_ScGeom_FrictPhys_CundallStrack::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_CundallStrack::traceEnergy` is true.")
+		.def("initPlasticDissipation",&Law2_ScGeom_FrictPhys_CundallStrack::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).")
 	);
 	FUNCTOR2D(ScGeom,FrictPhys);
 	DECLARE_LOGGER;

=== modified file 'pkg/dem/PeriIsoCompressor.cpp'
--- pkg/dem/PeriIsoCompressor.cpp	2010-11-13 21:11:39 +0000
+++ pkg/dem/PeriIsoCompressor.cpp	2010-12-06 11:54:38 +0000
@@ -1,5 +1,5 @@
 
-// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx> 
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
 
 #include<yade/pkg/dem/PeriIsoCompressor.hpp>
 #include<yade/pkg/dem/Shop.hpp>
@@ -29,7 +29,7 @@
 			if(!b || !b->bound) continue;
 			for(int i=0; i<3; i++) maxSpan=max(maxSpan,b->bound->max[i]-b->bound->min[i]);
 		}
-		
+
 	}
 	if(maxDisplPerStep<0) maxDisplPerStep=1e-2*charLen; // this should be tuned somehow…
 	const long& step=scene->iter;
@@ -99,9 +99,9 @@
 	//"Natural" strain, correct for large deformations, only used for comparison with goals
 	for (int i=0;i<3;i++) strain[i]=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 
+	// NOTE : needs refSize, could be generalized to arbitrary initial shapes using trsf*refHsize
 	// → initial cell size is always box, and will be. The cell repeats  periodically, initial shape doesn't (shouldn't, at least) dinfluence interactions at all/
 	// → this is one more place where Hsize would make things shorter : volume=Hsize.Determinant; The full code relies on the fact that initial Hsize is a box. You didn't modify the equation btw, result is the same as in r1936, which was (with spaces...)
 	//trsf*Matrix3r ( scene->cell->refSize[0],0,0, 0,scene->cell->refSize[1],0,0,0,scene->cell->refSize[2] ) ).Determinant()
@@ -122,8 +122,9 @@
 		//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=(Body::byId(I->getId1())->state->pos-Body::byId(I->getId2())->state->pos);
-		Vector3r branch= gsc->normal* ( gsc->refR1+gsc->refR2 );
+
+		Vector3r branch=(-Body::byId(I->getId1())->state->pos + Body::byId(I->getId2())->state->pos + scene->cell->Hsize*I->cellDist.cast<Real>());
+// 		Vector3r branch= gsc->normal* ( gsc->refR1+gsc->refR2 );
 		#if 0
 			// remove this block later
 			// tensorial product f*branch (hand-write the tensor product to prevent matrix instanciation inside the loop by makeTensorProduct)
@@ -195,18 +196,18 @@
 				strain_rate+=dampFactor*scene->dt* ( goal[axis]-stress[axis] ) /mass;
 				//if ((scene->iter%5000)==0){cerr << axis<<": stress="<<stress[axis]<<", goal="<<goal[axis]<<", velGrad="<<strain_rate<<endl;}
 				LOG_TRACE ( axis<<": stress="<<stress[axis]<<", goal="<<goal[axis]<<", velGrad="<<strain_rate );}
-			
+
 		} else {    // control strain, see "true strain" definition here http://en.wikipedia.org/wiki/Finite_strain_theory
 			///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.
 			strain_rate = (exp ( goal[axis]-strain[axis] ) -1)/scene->dt;
 			LOG_TRACE ( axis<<": strain="<<strain[axis]<<", goal="<<goal[axis]<<", cellGrow="<<strain_rate*scene->dt);
-		}		
+		}
 		// steady evolution with fluctuations; see TriaxialStressController
 		if (!dynCell) strain_rate=(1-growDamping)*strain_rate+.8*prevGrow[axis];
 		// limit maximum strain rate
 		if (abs(strain_rate)>maxStrainRate[axis]) strain_rate = Mathr::Sign(strain_rate)*maxStrainRate[axis];
 		// do not shrink below minimum cell size (periodic collider condition), although it is suboptimal WRT resulting stress
-		
+
 		//if ((scene->iter%5000)==0){cerr<< axis <<": velGrad="<<strain_rate<<", maxCellsize"<<-(cellSize[axis]-2.1*maxBodySpan[axis])/scene->dt<<endl;}
 		strain_rate=max(strain_rate,-(cellSize[axis]-2.1*maxBodySpan[axis])/scene->dt);
 		//if ((scene->iter%5000)==0){cerr <<"velGrad="<<strain_rate<<endl<<endl;}
@@ -230,7 +231,7 @@
 		}
 	}
 	// update stress and strain
-	if (!dynCell) for ( int axis=0; axis<3; axis++ ){		
+	if (!dynCell) for ( int axis=0; axis<3; axis++ ){
 		// take in account something like poisson's effect here…
 		//Real bogusPoisson=0.25; int ax1=(axis+1)%3,ax2=(axis+2)%3;
 		//don't modify stress if dynCell, testing only stiff[axis]>0 would not allow switching the control mode in simulations,
@@ -243,7 +244,7 @@
 	externalWork+=dW;
 	if(scene->trackEnergy) scene->energy->add(dW,"velGradWork",velGradWorkIx,/*non-incremental*/false);
 	prevGrow = strainRate;
-	
+
 	if(allOk){
 		if(doUpdate || currUnbalanced<0){
 			currUnbalanced=Shop::unbalancedForce(/*useMaxForce=*/false,scene);
@@ -302,7 +303,7 @@
 		// PATH_OP_OP(0,j,k) = path[0]->operator[](j).operator(k) is k-th element of j-th Vector2r of xxPath
 		#define PATH_OP_OP(pi,op1i,op2i) paths[pi]->operator[](op1i).operator()(op2i)
 
-		for (int i=0; i<6; i++) { 
+		for (int i=0; i<6; i++) {
 			for (int j=1; j<pathSizes[i]; j++) {
 				// check if the user defined time axis is monothonically increasing
 				{ if ( PATH_OP_OP(i,j-1,0) >= PATH_OP_OP(i,j,0) ) {
@@ -321,12 +322,12 @@
 
 		// set weather the simulation is "stress based" (all stress components prescribed or all prescribed strains equal zero)
 		if (lenPe == 0) { stressBasedSimulation = true; }
-		else { 
+		else {
 			stressBasedSimulation = true;
 			for (int i=0; i<lenPe; i++) { stressBasedSimulation = stressBasedSimulation && PATH_OP_OP(pe(i),1,1)<1e9; }
 		}
 	}
-	
+
 	// increase the pathCounter by one if we cross to the next part of path
 	for (int i=0; i<6; i++) {
 		if (progress >= PATH_OP_OP(i,pathsCounter[i],0)) { pathsCounter[i]++; }
@@ -412,7 +413,7 @@
 			}
 		}
 	}
-	
+
 	// correction coefficient ix strainRate.maxabs() > maxStrainRate
 	Real srCorr = (strainRate.cwise().abs().maxCoeff() > maxStrainRate)? (maxStrainRate/strainRate.cwise().abs().maxCoeff()):1.;
 	strainRate *= srCorr;
@@ -421,12 +422,12 @@
 	const Matrix3r& trsf=scene->cell->trsf;
 	// compute rotational and nonrotational (strain in local coordinates) part of trsf
 	Matrix_computeUnitaryPositive(trsf,&rot,&nonrot);
-	 
+
 	// prescribed velocity gradient (strain tensor rate) in global coordinates
-	epsilonRate = voigt_toSymmTensor(strainRate,/*strain=*/true); 
+	epsilonRate = voigt_toSymmTensor(strainRate,/*strain=*/true);
 	/* transformation of prescribed strain rate (computed by predictor) into local cell coordinates,
 	   multiplying by time to obtain strain increment and adding it to nonrot (current strain in local coordinates)*/
-	nonrot += rot.transpose()*(epsilonRate*dt)*rot; 
+	nonrot += rot.transpose()*(epsilonRate*dt)*rot;
 	Matrix3r& velGrad=scene->cell->velGrad;
 	// compute new trsf as rot*nonrot, substract actual trsf (= trsf increment), divide by dt (=trsf rate = velGrad
 	//trsf = rot*nonrot;
@@ -508,7 +509,7 @@
 	typedef Eigen::Matrix<Real,Eigen::Dynamic,Eigen::Dynamic> MatrixXr;
 	typedef Eigen::Matrix<Real,Eigen::Dynamic,1> VectorXr;
 	const Real& dt=scene->dt;
-	/ * 
+	/ *
 	sigma = K * eps
 
 	decompose the stiffness matrix depending on what is prescribed