fabricTensor(): re introducing all kind of interactions in the loop, with the...

fabricTensor(): re introducing all kind of interactions in the loop, with the new possibility defining a cutoff. See http://www.mail-archive.com/yade-dev@lists.launchpad.net/msg11982.html

pkg/dem/Shop.hpp

@@ -132,9 +132,9 @@ class Shop{
 		//! (contribution of normal and shear forces)
 		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, Real thresholdForce=NaN);
-		static py::tuple fabricTensor(bool splitTensor=false, Real thresholdForce=NaN);
+		//! Function to compute fabric tensor
+		static void fabricTensor(Real& Fmean, Matrix3r& fabric, Matrix3r& fabricStrong, Matrix3r& fabricWeak, Real cutoff=0.0, bool splitTensor=false, Real thresholdForce=NaN);
+		static py::tuple fabricTensor(Real cutoff=0.0, bool splitTensor=false, Real thresholdForce=NaN);
 		
 		//! Function to set translational and rotational velocities of all bodies to zero
 		static void calm(const shared_ptr<Scene>& rb=shared_ptr<Scene>(), int mask=-1);

pkg/dem/Shop_02.cpp

@@ -225,7 +225,7 @@ py::tuple Shop::normalShearStressTensors(bool compressionPositive, bool splitNor
 	
 	// *** 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;
-	fabricTensor(Fmean,f,fs,fw,false,NaN);
+	fabricTensor(Fmean,f,fs,fw); // 0,false,NaN for cutoff, split and thresholdForce as default arguments
 	Matrix3r sigNStrong(Matrix3r::Zero()), sigNWeak(Matrix3r::Zero());
 	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
 		if(!I->isReal()) continue;
@@ -257,55 +257,53 @@ py::tuple Shop::normalShearStressTensors(bool compressionPositive, bool splitNor
 
 /* Return the fabric tensor as according to [Satake1982]. */
 /* as side-effect, set Gl1_NormShear::strongWeakThresholdForce */
-void Shop::fabricTensor(Real& Fmean, Matrix3r& fabric, Matrix3r& fabricStrong, Matrix3r& fabricWeak, bool splitTensor, Real thresholdForce){
+void Shop::fabricTensor(Real& Fmean, Matrix3r& fabric, Matrix3r& fabricStrong, Matrix3r& fabricWeak, Real cutoff, bool splitTensor, Real thresholdForce){
 	Scene* scene=Omega::instance().getScene().get();
 
 	// *** Fabric tensor ***/
 	fabric=Matrix3r::Zero(); 
 	int count=0; // number of interactions
+	py::tuple aabb = Shop::aabbExtrema(cutoff);
+	Vector3r bbMin = py::extract<Vector3r>(aabb[0]), bbMax = py::extract<Vector3r>(aabb[1]);
+	Vector3r cp;
+	
+	Fmean=0; // initialize average contact force for split = 1 fabric measurements
+	// interactions loop to compute the fabric tensor returned when split = 0, and also measures average force for subsequent computations for split = 1:
 	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
 		if(!I->isReal()) continue;
-		if( !dynamic_cast<Sphere*>(Body::byId(I->getId1(),scene)->shape.get()) || !dynamic_cast<Sphere*>(Body::byId(I->getId2(),scene)->shape.get()) ) continue; // test intended to disregard boundary interactions (in non-periodic simulations)
 		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
+		cp = geom->contactPoint;
+		if( !(cp[0]>=bbMin[0] && cp[0]<=bbMax[0] && cp[1]>=bbMin[1] && cp[1]<=bbMax[1] && cp[2]>=bbMin[2] && cp[2]<=bbMax[2]) ) continue; // possible to use isInBB() from _utils.cpp ? (NB: would exclude the contact points exactly along the BB)
 		const Vector3r& n=geom->normal;
 		for(int i=0; i<3; i++) for(int j=i; j<3; j++){
 			fabric(i,j)+=n[i]*n[j];
 		}
+		NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->phys.get());
+		Real f=-phys->normalForce.dot(n);  // will be < 0 in compression
+		Fmean+=f;
 		count++;
 	}
+	Fmean/=count;
 	// fill terms under the diagonal
 	fabric(1,0)=fabric(0,1); fabric(2,0)=fabric(0,2); fabric(2,1)=fabric(1,2);
 	fabric/=count;
 	
-	// *** Average contact force ***/
-	// calculate average contact force
-	Fmean=0; // initialize
-	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
-		if(!I->isReal()) continue;
-		if( !dynamic_cast<Sphere*>(Body::byId(I->getId1(),scene)->shape.get()) || !dynamic_cast<Sphere*>(Body::byId(I->getId2(),scene)->shape.get()) ) continue; // test intended to disregard boundary interactions (in non-periodic simulations)
-		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
-		NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->phys.get());
-		const Vector3r& n=geom->normal;
-		Real f=-phys->normalForce.dot(n);  // will be < 0 in compression
-		Fmean+=f;
-	}
-	Fmean/=count; 
-
 	#ifdef YADE_OPENGL
 		Gl1_NormPhys::maxWeakFn=Fmean;
 	#endif
 	
 	// *** Weak and strong fabric tensors ***/
 	// evaluate two different parts of the fabric tensor 
-	// make distinction between strong and weak network of contact forces
+	// making distinction between strong and weak network of contact forces
 	fabricStrong=Matrix3r::Zero(); 
 	fabricWeak=Matrix3r::Zero(); 
 	int nStrong(0), nWeak(0); // number of strong and weak contacts respectively
 	if (!splitTensor & !std::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){
 		if(!I->isReal()) continue;
-		if( !dynamic_cast<Sphere*>(Body::byId(I->getId1(),scene)->shape.get()) || !dynamic_cast<Sphere*>(Body::byId(I->getId2(),scene)->shape.get()) ) continue; // test intended to disregard boundary interactions (in non-periodic simulations)
 		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
+		cp = geom->contactPoint;
+		if( !(cp[0]>=bbMin[0] && cp[0]<=bbMax[0] && cp[1]>=bbMin[1] && cp[1]<=bbMax[1] && cp[2]>=bbMin[2] && cp[2]<=bbMax[2]) ) continue; // see above for idea about isInBB() from _utils.cpp
 		NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->phys.get());
 		const Vector3r& n=geom->normal;
 		Real  f=-phys->normalForce.dot(n);
@@ -339,9 +337,9 @@ void Shop::fabricTensor(Real& Fmean, Matrix3r& fabric, Matrix3r& fabricStrong, M
 	}
 }
 
-py::tuple Shop::fabricTensor(bool splitTensor, Real thresholdForce){
+py::tuple Shop::fabricTensor(Real cutoff,bool splitTensor, Real thresholdForce){
 	Real Fmean; Matrix3r fabric, fabricStrong, fabricWeak;
-	fabricTensor(Fmean,fabric,fabricStrong,fabricWeak,splitTensor,thresholdForce);
+	fabricTensor(Fmean,fabric,fabricStrong,fabricWeak,cutoff,splitTensor,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);}

py/_utils.cpp

@@ -345,7 +345,7 @@ Real Shop__getPorosity(Real volume){ return Shop::getPorosity(Omega::instance().
 Real Shop__getVoxelPorosity(int resolution, Vector3r start,Vector3r end){ return Shop::getVoxelPorosity(Omega::instance().getScene(),resolution,start,end); }
 
 //Matrix3r Shop__stressTensorOfPeriodicCell(bool smallStrains=false){return Shop::stressTensorOfPeriodicCell(smallStrains);}
-py::tuple Shop__fabricTensor(bool splitTensor, Real thresholdForce){return Shop::fabricTensor(splitTensor,thresholdForce);}
+py::tuple Shop__fabricTensor(Real cutoff, bool splitTensor, Real thresholdForce){return Shop::fabricTensor(cutoff, splitTensor,thresholdForce);}
 py::tuple Shop__normalShearStressTensors(bool compressionPositive, bool splitNormalTensor, Real thresholdForce){return Shop::normalShearStressTensors(compressionPositive,splitNormalTensor,thresholdForce);}
 
 py::list Shop__getStressLWForEachBody(){return Shop::getStressLWForEachBody();}
@@ -493,7 +493,7 @@ BOOST_PYTHON_MODULE(_utils){
 	py::def("getViscoelasticFromSpheresInteraction",getViscoelasticFromSpheresInteraction,(py::arg("tc"),py::arg("en"),py::arg("es")),"Attention! The function is deprecated! Compute viscoelastic interaction parameters from analytical solution of a pair spheres collision problem:\n\n.. math:: 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=\\frac{2}{7}\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_t)^2\\right) \\\\ c_t=-\\frac{2}{7}\\frac{m}{t_c}\\ln e_t \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::getStress,(py::args("volume")=0),"Deprecated, use utils.getStress instead |ydeprecated|");
 	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: if true the function returns 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("thresholdForce")=NaN),"Computes the fabric tensor $F_{ij}=\\frac{1}{n_c}\\sum_c n_i n_j$ [Satake1982]_, for all sphere-sphere interactions $c$.\n\n:param bool splitTensor: split the fabric tensor into two parts related to the strong (greatest compressive normal forces) and weak contact forces respectively.\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. Use negative signed values for compressive states. To note that this value could be set to zero if one wanted to make distinction between compressive and tensile forces.");
+	py::def("fabricTensor",Shop__fabricTensor,(py::args("cutoff")=0.0,py::args("splitTensor")=false,py::args("thresholdForce")=NaN),"Computes the fabric tensor $F_{ij}=\\frac{1}{n_c}\\sum_c n_i n_j$ [Satake1982]_, for all interactions $c$.\n\n:param Real cutoff: intended to disregard boundary effects: to define in [0;1] to focus on the interactions located in the centered inner (1-cutoff)^3*$V$ part of the spherical packing $V$.\n\n:param bool splitTensor: split the fabric tensor into two parts related to the strong (greatest compressive normal forces) and weak contact forces respectively.\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. Use negative signed values for compressive states. 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,"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_{ik}\\sigma_{kj}=x_{i,k}\\sigma_{kj}=(x_{i}\\sigma_{kj})_{,k}-x_{i}\\sigma_{kj,k}$.\n\nAt equilibrium, the divergence of stress is null: $\\sigma_{kj,k}=\\vec{0}$. Consequently, after divergence theorem: $\\frac{1}{V}\\int_V \\sigma_{ij}dV = \\frac{1}{V}\\int_V (x_{i}\\sigma_{kj})_{,k}dV = \\frac{1}{V}\\int_{\\partial V}x_i\\sigma_{kj}n_kdS = \\frac{1}{V}\\sum_bx_i^bf_j^b$.\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\nThis last relation may not be very useful if porosity is not homogeneous. If it happens, one can define the equivalent bulk stress a the particles scale by assigning a volume to each particle. This volume can be obtained from :yref:`TesselationWrapper` (see e.g. [Catalano2014a]_)");
 	py::def("getStress",Shop::getStress,(py::args("volume")=0),"Compute and return Love-Weber stress tensor:\n\n $\\sigma_{ij}=\\frac{1}{V}\\sum_b f_i^b l_j^b$, where the sum is over all interactions, with $f$ the contact force and $l$ the branch vector (joining centers of the bodies). Stress is negativ for repulsive contact forces, i.e. compression. $V$ can be passed to the function. If it is not, it will be equal to the volume of the cell in periodic cases, or to the one deduced from utils.aabbDim() in non-periodic cases.");
 	py::def("getStressProfile",Shop::getStressProfile,(py::args("volume"),py::args("nCell"),py::args("dz"),py::args("zRef"),py::args("vPartAverageX"),py::args("vPartAverageY"),py::args("vPartAverageZ")),"Compute and return the stress tensor depth profile, including the contribution from Love-Weber stress tensor and the dynamic stress tensor taking into account the effect of particles inertia. For each defined cell z, the stress tensor reads: \n\n $\\sigma_{ij}^z= \\frac{1}{V}\\sum_c f_i^c l_j^{c,z} - \\frac{1}{V}\\sum_p m^p u'^p_i u'^p_j$,\n\n where the first sum is made over the contacts which are contained or cross the cell z, f^c is the contact force from particle 1 to particle 2, and l^{c,z} is the part of the branch vector from particle 2 to particle 1, contained in the cell. The second sum is made over the particles, and u'^p is the velocity fluctuations of the particle p with respect to the spatial averaged particle velocity at this point (given as input parameters). The expression of the stress tensor is the same as the one given in getStress plus the inertial contribution. Apart from that, the main difference with getStress stands in the fact that it gives a depth profile of stress tensor, i.e. from the reference horizontal plane at elevation zRef (input parameter) until the plane of elevation zRef+nCell*dz (input parameters), it is computing the stress tensor for each cell of height dz. For the love-Weber stress contribution, the branch vector taken into account in the calculations is only the part of the branch vector contained in the cell considered.\n To validate the formulation, it has been checked that activating only the Love-Weber stress tensor, and suming all the contributions at the different altitude, we recover the same stress tensor as when using getStress. For my own use, I have troubles with strong overlap between fixed object, so that I made a condition to exclude the contribution to the stress tensor of the fixed objects, this can be desactivated easily if needed (and should be desactivated for the comparison with getStress)." );

py/_utils.hpp

@@ -126,7 +126,7 @@ Real Shop__getPorosity(Real volume=-1);
 Real Shop__getVoxelPorosity(int resolution=200, Vector3r start=Vector3r(0,0,0),Vector3r end=Vector3r(0,0,0));
 
 //Matrix3r Shop__stressTensorOfPeriodicCell(bool smallStrains=false){return Shop::stressTensorOfPeriodicCell(smallStrains);}
-py::tuple Shop__fabricTensor(bool splitTensor=false, Real thresholdForce=NaN);
+py::tuple Shop__fabricTensor(Real cutoff=0.0,bool splitTensor=false, Real thresholdForce=NaN);
 py::tuple Shop__normalShearStressTensors(bool compressionPositive=false, bool splitNormalTensor=false, Real thresholdForce=NaN);
 
 py::list Shop__getStressLWForEachBody();