yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2719
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
branch nick: yade
timestamp: Wed 2011-02-09 11:01:01 +0000
message:
  - Average stress tensor of periodic cell can also be spit into two contributions, given a threshold value (mean force by default).
modified:
  pkg/dem/Shop.cpp
  pkg/dem/Shop.hpp
  py/_utils.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/Shop.cpp'
--- pkg/dem/Shop.cpp	2011-02-08 16:37:12 +0000
+++ pkg/dem/Shop.cpp	2011-02-09 11:01:01 +0000
@@ -500,7 +500,9 @@
 The formulation follows the [Thornton2000]_ article
 "Numerical simulations of deviatoric shear deformation of granular media", eq (3) and (4)
  */
-py::tuple Shop::normalShearStressTensors(bool compressionPositive){
+py::tuple Shop::normalShearStressTensors(bool compressionPositive, bool splitNormalTensor, Real thresholdForce){
+  
+	//*** Stress tensor split into shear and normal contribution ***/
 	Scene* scene=Omega::instance().getScene().get();
 	if (!scene->isPeriodic){ throw runtime_error("Can't compute stress of periodic cell in aperiodic simulation."); }
 	Matrix3r sigN(Matrix3r::Zero()),sigT(Matrix3r::Zero());
@@ -528,17 +530,48 @@
 	// fill terms under the diagonal
 	sigN(1,0)=sigN(0,1); sigN(2,0)=sigN(0,2); sigN(2,1)=sigN(1,2);
 	sigT(1,0)=sigT(0,1); sigT(2,0)=sigT(0,2); sigT(2,1)=sigT(1,2);
-	return py::make_tuple(sigN,sigT);
+	
+	// *** Normal stress tensor split into two parts according to subnetworks of strong and weak forces (or other distinction if a threshold value for the force is assigned) ***/
+	Real Fmean(0); Matrix3r f, fs, fw;
+	/// FIXME: this does not compute Fmean as it should?
+	fabricTensor(Fmean,f,fs,fw,false,compressionPositive,NaN);
+	Matrix3r sigNStrong(Matrix3r::Zero()), sigNWeak(Matrix3r::Zero());
+	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+		if(!I->isReal()) continue;
+		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
+		NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->phys.get());
+		const Vector3r& n=geom->normal;
+		Real N=(compressionPositive?-1:1)*phys->normalForce.dot(n);
+		// Real R=(Body::byId(I->getId2(),scene)->state->pos+cellHsize*I->cellDist.cast<Real>()-Body::byId(I->getId1(),scene)->state->pos).norm();
+		Real R=.5*(geom->refR1+geom->refR2);
+		Real Fsplit=(!isnan(thresholdForce))?thresholdForce:Fmean;
+		if (compressionPositive?(N<Fsplit):(N>Fsplit)){
+			for(int i=0; i<3; i++) for(int j=i; j<3; j++){
+				sigNStrong(i,j)+=R*N*n[i]*n[j];}
+		}
+		else{
+			for(int i=0; i<3; i++) for(int j=i; j<3; j++){
+				sigNWeak(i,j)+=R*N*n[i]*n[j];}
+		}
+	}
+	sigNStrong*=2/vol; sigNWeak*=2/vol;
+	// fill terms under the diagonal
+	sigNStrong(1,0)=sigNStrong(0,1); sigNStrong(2,0)=sigNStrong(0,2); sigNStrong(2,1)=sigNStrong(1,2);
+	sigNWeak(1,0)=sigNWeak(0,1); sigNWeak(2,0)=sigNWeak(0,2); sigNWeak(2,1)=sigNWeak(1,2);
+	
+	/// tensile forces are taken as positive!
+	if (splitNormalTensor){return py::make_tuple(sigNStrong,sigNWeak);} // return strong-weak or tensile-compressive parts of the stress tensor (only normal part)
+	return py::make_tuple(sigN,sigT); // return normal and shear components
 }
 
 /* Return the fabric tensor as according to [Satake1982]. */
 /* as side-effect, set Gl1_NormShear::strongWeakThresholdForce */
-py::tuple Shop::fabricTensor(bool splitTensor, bool revertSign, Real thresholdForce){
+void Shop::fabricTensor(Real& Fmean, Matrix3r fabric, Matrix3r fabricStrong, Matrix3r fabricWeak, bool splitTensor, bool revertSign, Real thresholdForce){
 	Scene* scene=Omega::instance().getScene().get();
 	if (!scene->isPeriodic){ throw runtime_error("Can't compute fabric tensor of periodic cell in aperiodic simulation."); }
 	
 	// *** Fabric tensor ***/
-	Matrix3r fabric(Matrix3r::Zero()); 
+	fabric=Matrix3r::Zero(); 
 	int count=0; // number of interactions
 	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
 		if(!I->isReal()) continue;
@@ -555,7 +588,7 @@
 	
 	// *** Average contact force ***/
 	// calculate average contact force
-	Real Fmean(0);
+	Fmean=0; // initialize
 	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
 		if(!I->isReal()) continue;
 		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
@@ -574,7 +607,8 @@
 	// *** Weak and strong fabric tensors ***/
 	// evaluate two different parts of the fabric tensor 
 	// make distinction between strong and weak network of contact forces
-	Matrix3r fabricStrong(Matrix3r::Zero()), fabricWeak(Matrix3r::Zero()); 
+	fabricStrong=Matrix3r::Zero(); 
+	fabricWeak=Matrix3r::Zero(); 
 	int nStrong(0), nWeak(0); // number of strong and weak contacts respectively
 	if (!splitTensor & !isnan(thresholdForce)) {LOG_WARN("The bool splitTensor should be set to True if you specified a threshold value for the contact force, otherwise the function will return only the fabric tensor and not the two separate contributions.");}
 	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
@@ -611,7 +645,11 @@
 		q=nStrong*1./count; 
 		fabricTot=(1-q)*fabricWeak+q*fabricStrong;
 	}
-	
+}
+
+py::tuple Shop::fabricTensor(bool splitTensor, bool revertSign, Real thresholdForce){
+	Real Fmean; Matrix3r fabric, fabricStrong, fabricWeak;
+	fabricTensor(Fmean,fabric,fabricStrong,fabricWeak,splitTensor,revertSign,thresholdForce);
 	// returns fabric tensor or alternatively the two distinct contributions according to strong and weak subnetworks (or, if thresholdForce is specified, the distinction is made according to that value and not the mean one)
 	if (!splitTensor){return py::make_tuple(fabric);} 
 	else{return py::make_tuple(fabricStrong,fabricWeak);}

=== modified file 'pkg/dem/Shop.hpp'
--- pkg/dem/Shop.hpp	2011-02-08 16:37:12 +0000
+++ pkg/dem/Shop.hpp	2011-02-09 11:01:01 +0000
@@ -112,8 +112,9 @@
 		static Matrix3r stressTensorOfPeriodicCell(bool smallStrains=true);
 		//! Compute overall ("macroscopic") stress of periodic cell, returning 2 tensors
 		//! (contribution of normal and shear forces)
-		static py::tuple normalShearStressTensors(bool compressionPositive=false);
+		static py::tuple normalShearStressTensors(bool compressionPositive=false, bool splitNormalTensor=false, Real thresholdForce=NaN);
 		
 		//! Function to compute fabric tensor of periodic cell
+		static void fabricTensor(Real& Fmean, Matrix3r fabric, Matrix3r fabricStrong, Matrix3r fabricWeak, bool splitTensor=false, bool revertSign=false, Real thresholdForce=NaN);
 		static py::tuple fabricTensor(bool splitTensor=false, bool revertSign=false, Real thresholdForce=NaN);
 };

=== modified file 'py/_utils.cpp'
--- py/_utils.cpp	2011-02-08 16:37:12 +0000
+++ py/_utils.cpp	2011-02-09 11:01:01 +0000
@@ -430,7 +430,8 @@
 
 Matrix3r Shop__stressTensorOfPeriodicCell(bool smallStrains=false){return Shop::stressTensorOfPeriodicCell(smallStrains);}
 py::tuple Shop__fabricTensor(bool splitTensor=false, bool revertSign=false, Real thresholdForce=NaN){return Shop::fabricTensor(splitTensor,revertSign,thresholdForce);}
-py::tuple Shop__normalShearStressTensors(bool compressionPositive=false){return Shop::normalShearStressTensors(compressionPositive);}
+py::tuple Shop__normalShearStressTensors(bool compressionPositive=false, bool splitNormalTensor=false, Real thresholdForce=NaN){return Shop::normalShearStressTensors(compressionPositive,splitNormalTensor,thresholdForce);}
+
 py::tuple Shop__getStressLWForEachBody(bool revertSign=false){return Shop::getStressLWForEachBody(revertSign);}
 
 BOOST_PYTHON_MODULE(_utils){
@@ -475,7 +476,7 @@
 	py::def("flipCell",&Shop::flipCell,(py::arg("flip")=Matrix3r(Matrix3r::Zero())),"Flip periodic cell so that angles between $R^3$ axes and transformed axes are as small as possible. This function relies on the fact that periodic cell defines by repetition or its corners regular grid of points in $R^3$; however, all cells generating identical grid are equivalent and can be flipped one over another. This necessiatates adjustment of :yref:`Interaction.cellDist` for interactions that cross boundary and didn't before (or vice versa), and re-initialization of collider. The *flip* argument can be used to specify desired flip: integers, each column for one axis; if zero matrix, best fit (minimizing the angles) is computed automatically.\n\nIn c++, this function is accessible as ``Shop::flipCell``.\n\n.. warning::\n\t This function is currently broken and should not be used.");
 	py::def("getViscoelasticFromSpheresInteraction",getViscoelasticFromSpheresInteraction,(py::arg("tc"),py::arg("en"),py::arg("es")),"Compute viscoelastic interaction parameters from analytical solution of a pair spheres collision problem:\n\n\n.. math::\n\t:nowrap:\n\n\n\t\\begin{align*}k_n&=\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_n)^2\\right)\\\\  c_n&=-\\frac{2m}{t_c}\\ln e_n \\\\k_t&=\\frac27\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_t)^2\\right)  \\\\ c_t=-\\frac27\\frac{m}{t_c}\\ln e_t \\end{align*}\n\n\nwhere $k_n$, $c_n$ are normal elastic and viscous coefficients and $k_t$, $c_t$ shear elastic and viscous coefficients. For details see [Pournin2001]_.\n\n:param float m: sphere mass $m$\n:param float tc: collision time $t_c$\n:param float en: normal restitution coefficient $e_n$\n:param float es: tangential restitution coefficient $e_s$\n:return: dictionary with keys ``kn`` (the value of $k_n$), ``cn`` ($c_n$), ``kt`` ($k_t$), ``ct`` ($c_t$).");
 	py::def("stressTensorOfPeriodicCell",Shop__stressTensorOfPeriodicCell,(py::args("smallStrains")=false),"Compute overall (macroscopic) stress of periodic cell using equation published in [Kuhl2001]_:\n\n.. math:: \\vec{\\sigma}=\\frac{1}{V}\\sum_cl^c[\\vec{N}^cf_N^c+\\vec{T}^{cT}\\cdot\\vec{f}^c_T],\n\nwhere $V$ is volume of the cell, $l^c$ length of interaction $c$, $f^c_N$ normal force and $\\vec{f}^c_T$ shear force. Sumed are values over all interactions $c$. $\\vec{N}^c$ and $\\vec{T}^{cT}$ are projection tensors (see the original publication for more details):\n\n.. math:: \\vec{N}=\\vec{n}\\otimes\\vec{n}\\rightarrow N_{ij}=n_in_j\n\n.. math:: \\vec{T}^T=\\vec{I}_{sym}\\cdot\\vec{n}-\\vec{n}\\otimes\\vec{n}\\otimes\\vec{n}\\rightarrow T^T_{ijk}=\\frac{1}{2}(\\delta_{ik}\\delta_{jl}+\\delta_{il}\\delta_{jk})n_l-n_in_jn_k\n\n.. math:: \\vec{T}^T\\cdot\\vec{f}_T\\equiv T^T_{ijk}f_k=(\\delta_{ik}n_j/2+\\delta_{jk}n_i/2-n_in_jn_k)f_k=n_jf_i/2+n_if_j/2-n_in_jn_kf_k,\n\nwhere $n$ is unit vector oriented along the interaction (:yref:`normal<GenericSpheresContact::normal>`) and $\\delta$ is Kronecker's delta. As $\\vec{n}$ and $\\vec{f}_T$ are perpendicular (therfore $n_if_i=0$) we can write\n\n.. math:: \\sigma_{ij}=\\frac{1}{V}\\sum l[n_in_jf_N+n_jf^T_i/2+n_if^T_j/2]\n\n:param bool smallStrains: if false (large strains), real values of volume and interaction lengths are computed. If true, only :yref:`refLength<Dem3DofGeom::refLength>` of interactions and initial volume are computed (can save some time).\n\n:return: macroscopic stress tensor as Matrix3");
-	py::def("normalShearStressTensors",Shop__normalShearStressTensors,(py::args("compressionPositive")=false),"Compute overall stress tensor of the periodic cell decomposed in 2 parts, one contributed by normal forces, the other by shear forces. The formulation can be found in [Thornton2000]_, eq. (3):\n\n.. math:: \\tens{\\sigma}_{ij}=\\frac{2}{V}\\sum R N \\vec{n}_i \\vec{n}_j+\\frac{2}{V}\\sum R T \\vec{n}_i\\vec{t}_j\n\nwhere $V$ is the cell volume, $R$ is \"contact radius\" (in our implementation, current distance between particle centroids), $\\vec{n}$ is the normal vector, $\\vec{t}$ is a vector perpendicular to $\\vec{n}$, $N$ and $T$ are norms of normal and shear forces.");
+	py::def("normalShearStressTensors",Shop__normalShearStressTensors,(py::args("compressionPositive")=false,py::args("splitNormalTensor")=false,py::args("thresholdForce")=NaN),"Compute overall stress tensor of the periodic cell decomposed in 2 parts, one contributed by normal forces, the other by shear forces. The formulation can be found in [Thornton2000]_, eq. (3):\n\n.. math:: \\tens{\\sigma}_{ij}=\\frac{2}{V}\\sum R N \\vec{n}_i \\vec{n}_j+\\frac{2}{V}\\sum R T \\vec{n}_i\\vec{t}_j\n\nwhere $V$ is the cell volume, $R$ is \"contact radius\" (in our implementation, current distance between particle centroids), $\\vec{n}$ is the normal vector, $\\vec{t}$ is a vector perpendicular to $\\vec{n}$, $N$ and $T$ are norms of normal and shear forces.\n\n:param bool splitNormalTensor: return normal stress tensor split into two parts according to the two subnetworks of strong an weak forces.\n\n:param Real thresholdForce: threshold value according to which the normal stress tensor can be split (e.g. a zero value would make distinction between tensile and compressive forces).");
 	py::def("fabricTensor",Shop__fabricTensor,(py::args("splitTensor")=false,py::args("revertSign")=false,py::args("thresholdForce")=NaN),"Compute the fabric tensor of the periodic cell. The original paper can be found in [Satake1982]_.\n\n:param bool splitTensor: split the fabric tensor into two parts related to the strong and weak contact forces respectively.\n\n:param bool revertSign: it must be set to true if the contact law's convention takes compressive forces as positive.\n\n:param Real thresholdForce: if the fabric tensor is split into two parts, a threshold value can be specified otherwise the mean contact force is considered by default. It is worth to note that this value has a sign and the user needs to set it according to the convention adopted for the contact law. To note that this value could be set to zero if one wanted to make distinction between compressive and tensile forces.");
 	py::def("bodyStressTensors",Shop__getStressLWForEachBody,(py::args("revertSign")=false),"Compute and return a table with per-particle stress tensors. Each tensor represents the average stress in one particle, obtained from the contour integral of applied load as detailed below. This definition is considering each sphere as a continuum. It can be considered exact in the context of spheres at static equilibrium, interacting at contact points with negligible volume changes of the solid phase (this last assumption is not restricting possible deformations and volume changes at the packing scale).\n\nProof:\n\nFirst, we remark the identity:  $\\sigma_{ij}=\\delta_{ij}\\sigma_{ij}=x_{i,j}\\sigma_{ij}=(x_{i}\\sigma_{ij})_{,j}-x_{i}\\sigma_{ij,j}$.\n\nAt equilibrium, the divergence of stress is null: $\\sigma_{ij,j}=\\vec{0}$. Consequently, after divergence theorem: $\\frac{1}{V}\\int_V \\sigma_{ij}dV = \\frac{1}{V}\\int_V (x_{i}\\sigma_{ij})_{,j}dV = \\frac{1}{V}\\int_{\\partial V}x_i.\\sigma_{ij}.\\vec{n_j}.dS = \\frac{1}{V}\\sum_kx_i^k.f_j^k$.\n\nThe last equality is implicitely based on the representation of external loads as Dirac distributions whose zeros are the so-called *contact points*: 0-sized surfaces on which the *contact forces* are applied, located at $x_i$ in the deformed configuration.\n\nA weighted average of per-body stresses will give the average stress inside the solid phase. There is a simple relation between the stress inside the solid phase and the stress in an equivalent continuum in the absence of fluid pressure. For porosity $n$, the relation reads: $\\sigma_{ij}^{equ.}=(1-n)\\sigma_{ij}^{solid}$.\n\n:param bool revertSign: invert the sign of returned tensors components.");
 	py::def("maxOverlapRatio",maxOverlapRatio,"Return maximum overlap ration in interactions (with :yref:`ScGeom`) of two :yref:`spheres<Sphere>`. The ratio is computed as $\\frac{u_N}{2(r_1 r_2)/r_1+r_2}$, where $u_N$ is the current overlap distance and $r_1$, $r_2$ are radii of the two spheres in contact.");