yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2433
committer: Jan Stránský <honzik@honzik-laptop>
branch nick: yade
timestamp: Sun 2010-09-12 10:44:00 +0200
message:
  1. add function yade.utils.stressTensorOfPeriodicCell
  2. removed BoundDispatcher from pack.randomDensePack
modified:
  pkg/dem/meta/Shop.cpp
  pkg/dem/meta/Shop.hpp
  py/_utils.cpp
  py/pack/pack.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/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp	2010-09-10 11:41:37 +0000
+++ pkg/dem/meta/Shop.cpp	2010-09-12 08:44:00 +0000
@@ -1203,3 +1203,41 @@
 		bodyStates[id2].shearStress+=shearStress;
 	}
 }
+
+Matrix3r Shop::stressTensorOfPeriodicCell(bool smallStrains){
+	Scene* scene=Omega::instance().getScene().get();
+	if (!scene->isPeriodic){ throw runtime_error("Can't compute stress of periodic cell in aperiodic simulation."); }
+	Real volume;
+	if (smallStrains){volume = scene->cell->refSize[0]*scene->cell->refSize[1]*scene->cell->refSize[2];}
+	else volume = scene->cell->trsf.determinant()*scene->cell->refSize[0]*scene->cell->refSize[1]*scene->cell->refSize[2];
+	Matrix3r stress = Matrix3r::Zero();
+	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+		if(!I->isReal()) continue;
+		Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());
+		NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->interactionPhysics.get());
+		Real l;
+		if (smallStrains){l = geom->refLength;}
+		else l=(geom->se31.position-geom->se32.position).norm();
+		Vector3r& n=geom->normal;
+		Vector3r& fT=phys->shearForce;
+		Real fN=phys->normalForce.dot(n);
+		Real fT0=fT[0]; Real fT1=fT[1]; Real fT2=fT[2];
+		Real n0=n[0]; Real n1=n[1]; Real n2=n[2];
+		
+		Real s00 = n0*n0*fN + fT0*n0;
+		Real s01 = n0*n1*fN + .5*(fT0*n1 + fT1*n0);
+		Real s02 = n0*n2*fN + .5*(fT0*n2 + fT2*n0);
+		Real s11 = n1*n1*fN + fT1*n1;
+		Real s12 = n1*n2*fN + .5*(fT1*n2 + fT2*n1);
+		Real s22 = n2*n2*fN + fT2*n2;
+		
+		stress(0,0) += s00;
+		stress(0,1) += stress(1,0) = s01;
+		stress(0,2) += stress(2,0) = s02;
+		stress(1,1) += s11;
+		stress(1,2) += stress(2,1) = s12;
+		stress(2,2) += s22;
+	}
+	stress/=volume;
+	return stress;
+}

=== modified file 'pkg/dem/meta/Shop.hpp'
--- pkg/dem/meta/Shop.hpp	2010-09-09 09:51:23 +0000
+++ pkg/dem/meta/Shop.hpp	2010-09-12 08:44:00 +0000
@@ -102,6 +102,7 @@
 		static Real kineticEnergy(Scene* _rb=NULL, Body::id_t* maxId=NULL);
 		static Vector3r totalForceInVolume(Real& avgIsoStiffness, Scene *_rb=NULL);
 
+
 		//! create transientInteraction between 2 bodies, using existing Dispatcher in Omega
 		static shared_ptr<Interaction> createExplicitInteraction(Body::id_t id1, Body::id_t id2, bool force);
 
@@ -129,4 +130,6 @@
 		//! Function of getting stresses for each body
 		static void getStressForEachBody(vector<Shop::bodyState>&);
 
+		//! Function to compute overall ("macroscopic") stress of periodic cell
+		static Matrix3r stressTensorOfPeriodicCell(bool smallStrains=true);
 };

=== modified file 'py/_utils.cpp'
--- py/_utils.cpp	2010-09-09 09:51:23 +0000
+++ py/_utils.cpp	2010-09-12 08:44:00 +0000
@@ -407,6 +407,8 @@
 
 Real Shop__getPorosity(Real volume=-1){ return Shop::getPorosity(Omega::instance().getScene(),volume); }
 
+Matrix3r Shop__stressTensorOfPeriodicCell(bool smallStrains=false){return Shop::stressTensorOfPeriodicCell(smallStrains);}
+
 BOOST_PYTHON_MODULE(_utils){
 	// http://numpy.scipy.org/numpydoc/numpy-13.html mentions this must be done in module init, otherwise we will crash
 	import_array();
@@ -448,4 +450,5 @@
 	py::def("wireNoSpheres",wireNoSpheres,"Set :yref:`Shape::wire` to True on non-spherical bodies (:yref:`Facets<Facet>`, :yref:`Walls<Wall>`).");
 	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("m"),py::arg("tc"),py::arg("en"),py::arg("es")),"Get viscoelastic interaction parameters from analytical solution of a pair spheres collision problem. \n\n:Parameters:\n\t`m` : float\n\t\tsphere mass\n\t`tc` : float\n\t\tcollision time\n\t`en` : float\n\t\tnormal restitution coefficient\n\t`es` : float\n\t\ttangential restitution coefficient.\n\n\n:return: \n\tdict with keys:\n\n\tkn : float\n\t\tnormal elastic coefficient computed as:\n\n.. math:: k_n=\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_n)^2\\right)\n\n\tcn : float\n\t\tnormal viscous coefficient computed as:\n\n.. math:: c_n=-\\frac{2m}{t_c}\\ln e_n\n\n\n\tkt : float\n\t\ttangential elastic coefficient computed as:\n\n.. math:: k_t=\\frac27\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_t)^2\\right)\n\n\tct : float\n\t\ttangential viscous coefficient computed as:\n\n.. math:: c_t=-\\frac27\\frac{m}{t_c}\\ln e_t.\n\n\nFor details see [Pournin2001]_. ");
+	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");
 }

=== modified file 'py/pack/pack.py'
--- py/pack/pack.py	2010-09-09 14:02:49 +0000
+++ py/pack/pack.py	2010-09-12 08:44:00 +0000
@@ -410,7 +410,7 @@
 		#print x1,y1,z1,radius,rRelFuzz
 		O.materials.append(FrictMat(young=3e10,density=2400))
 		num=sp.makeCloud(Vector3().Zero,O.cell.refSize,radius,rRelFuzz,spheresInCell,True)
-		O.engines=[ForceResetter(),BoundDispatcher([Bo1_Sphere_Aabb()]),InsertionSortCollider(nBins=5,sweepLength=.05*radius),InteractionDispatchers([Ig2_Sphere_Sphere_Dem3DofGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_Basic()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5,keepProportions=True),NewtonIntegrator(damping=.6)]
+		O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb()],nBins=5,sweepLength=.05*radius),InteractionDispatchers([Ig2_Sphere_Sphere_Dem3DofGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_Basic()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5,keepProportions=True),NewtonIntegrator(damping=.6)]
 		O.materials.append(FrictMat(young=30e9,frictionAngle=.5,poisson=.3,density=1e3))
 		for s in sp: O.bodies.append(utils.sphere(s[0],s[1]))
 		O.dt=utils.PWaveTimeStep()