yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2657
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Mon 2011-01-17 16:44:53 +0100
message:
  - Add a function (and python wrapper) computing the exact mean stress in each sphere.
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-01-17 08:49:26 +0000
+++ pkg/dem/Shop.cpp	2011-01-17 15:44:53 +0000
@@ -572,3 +572,26 @@
 	stress/=volume;
 	return stress;
 }
+
+void Shop::getStressLWForEachBody(vector<Matrix3r>& bStresses, bool revertSign){
+	const shared_ptr<Scene>& scene=Omega::instance().getScene();
+	bStresses.resize(scene->bodies->size());
+	for (int k=0;k<scene->bodies->size();k++) bStresses[k]=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());
+		Vector3r f=phys->normalForce+phys->shearForce;
+		if (revertSign) f=-f;//revert sign for some laws, so that every interaction type can give the same result
+		//Sum f_i*l_j/V for each contact of each particle
+		bStresses[I->getId1()]-=(3.0/(4.0*Mathr::PI*pow(geom->refR1,3)))*f*((geom->contactPoint-Body::byId(I->getId1(),scene)->state->pos).transpose());
+		bStresses[I->getId2()]+=(3.0/(4.0*Mathr::PI*pow(geom->refR2,3)))*f*((geom->contactPoint-Body::byId(I->getId2(),scene)->state->pos).transpose());
+		totStress+=f*((Body::byId(I->getId1(),scene)->state->pos - Body::byId(I->getId2(),scene)->state->pos).transpose());
+	}
+}
+
+py::tuple Shop::getStressLWForEachBody(bool revertSign){
+	vector<Matrix3r> bStresses;
+	getStressLWForEachBody(bStresses,revertSign);
+	return py::make_tuple(bStresses);
+}

=== modified file 'pkg/dem/Shop.hpp'
--- pkg/dem/Shop.hpp	2011-01-09 16:34:50 +0000
+++ pkg/dem/Shop.hpp	2011-01-17 15:44:53 +0000
@@ -103,6 +103,10 @@
 		};
 		//! Function of getting stresses for each body
 		static void getStressForEachBody(vector<Shop::bodyState>&);
+		
+		//! Define the exact average stress in each particle from contour integral ("LW" stands for Love-Weber, since this is what the contour integral gives).
+		static void getStressLWForEachBody(vector<Matrix3r>& bStresses, bool revertSign=false);
+		static py::tuple getStressLWForEachBody(bool revertSign);
 
 		//! Function to compute overall ("macroscopic") stress of periodic cell
 		static Matrix3r stressTensorOfPeriodicCell(bool smallStrains=true);

=== modified file 'py/_utils.cpp'
--- py/_utils.cpp	2010-12-23 14:01:33 +0000
+++ py/_utils.cpp	2011-01-17 15:44:53 +0000
@@ -426,6 +426,7 @@
 
 Matrix3r Shop__stressTensorOfPeriodicCell(bool smallStrains=false){return Shop::stressTensorOfPeriodicCell(smallStrains);}
 py::tuple Shop__normalShearStressTensors(bool compressionPositive=false){return Shop::normalShearStressTensors(compressionPositive);}
+py::tuple Shop__getStressLWForEachBody(bool revertSign=false){return Shop::getStressLWForEachBody(revertSign);}
 
 BOOST_PYTHON_MODULE(_utils){
 	// http://numpy.scipy.org/numpydoc/numpy-13.html mentions this must be done in module init, otherwise we will crash
@@ -470,5 +471,6 @@
 	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("bodyStressTensors",Shop__getStressLWForEachBody,(py::args("revertSign")=false),"Compute the exact mean stress tensor in each sphere from the contour integral of applied load.\n\nAfter divergence theorem, at equilibrium: $\\int_V s_{ij}dV = \\int_{S_V} x_i.s_{ij}.n_j.dS = \\sum_kx_i^k.f_j^k$. This relation applies for arbitrary shapes but the result has to be divided by the solid's volume, computed here using the radii, hence assuming spheres. The (weighted) average of per-body stresses is exactly equal to the average stress in the solid phase, i.e. $\\sigma_{ij}^{macro}/compacity$.");
 	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.");
 }