~ubuntu-branches/ubuntu/trusty/yade/trusty

« back to all changes in this revision

Viewing changes to py/_utils.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky, cf3f8d9
  • Date: 2013-10-30 20:56:33 UTC
  • mfrom: (20.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20131030205633-1f01r7hjce17d723
Tags: 1.05.0-2
[cf3f8d9] Pass -ftrack-macro-expansion=0 only if gcc>=4.8. (Closes: #726009)

Show diffs side-by-side

added added

removed removed

Lines of Context:
7
7
#include<yade/pkg/dem/ScGeom.hpp>
8
8
#include<yade/pkg/dem/DemXDofGeom.hpp>
9
9
#include<yade/pkg/common/Facet.hpp>
 
10
#include<yade/pkg/dem/Tetra.hpp>
10
11
#include<yade/pkg/common/Sphere.hpp>
11
12
#include<yade/pkg/common/NormShearPhys.hpp>
12
13
#include<yade/lib/computational-geometry/Hull2d.hpp>
75
76
Real PWaveTimeStep(){return Shop::PWaveTimeStep();};
76
77
Real RayleighWaveTimeStep(){return Shop::RayleighWaveTimeStep();};
77
78
 
78
 
Real elasticEnergyInAABB(py::tuple Aabb){
79
 
        Vector3r bbMin=py::extract<Vector3r>(Aabb[0])(), bbMax=py::extract<Vector3r>(Aabb[1])();
80
 
        shared_ptr<Scene> rb=Omega::instance().getScene();
81
 
        Real E=0;
82
 
        FOREACH(const shared_ptr<Interaction>&i, *rb->interactions){
83
 
                if(!i->phys) continue;
84
 
                shared_ptr<NormShearPhys> bc=dynamic_pointer_cast<NormShearPhys>(i->phys); if(!bc) continue;
85
 
                shared_ptr<Dem3DofGeom> geom=dynamic_pointer_cast<Dem3DofGeom>(i->geom); if(!geom){LOG_ERROR("NormShearPhys contact doesn't have Dem3DofGeom associated?!"); continue;}
86
 
                const shared_ptr<Body>& b1=Body::byId(i->getId1(),rb), b2=Body::byId(i->getId2(),rb);
87
 
                bool isIn1=isInBB(b1->state->pos,bbMin,bbMax), isIn2=isInBB(b2->state->pos,bbMin,bbMax);
88
 
                if(!isIn1 && !isIn2) continue;
89
 
                LOG_DEBUG("Interaction #"<<i->getId1()<<"--#"<<i->getId2());
90
 
                Real weight=1.;
91
 
                if((!isIn1 && isIn2) || (isIn1 && !isIn2)){
92
 
                        //shared_ptr<Body> bIn=isIn1?b1:b2, bOut=isIn2?b2:b1;
93
 
                        Vector3r vIn=(isIn1?b1:b2)->state->pos, vOut=(isIn2?b1:b2)->state->pos;
94
 
                        #define _WEIGHT_COMPONENT(axis) if(vOut[axis]<bbMin[axis]) weight=min(weight,abs((vOut[axis]-bbMin[axis])/(vOut[axis]-vIn[axis]))); else if(vOut[axis]>bbMax[axis]) weight=min(weight,abs((vOut[axis]-bbMax[axis])/(vOut[axis]-vIn[axis])));
95
 
                        _WEIGHT_COMPONENT(0); _WEIGHT_COMPONENT(1); _WEIGHT_COMPONENT(2);
96
 
                        assert(weight>=0 && weight<=1);
97
 
                        //cerr<<"Interaction crosses Aabb boundary, weight is "<<weight<<endl;
98
 
                        //LOG_DEBUG("Interaction crosses Aabb boundary, weight is "<<weight);
99
 
                } else {assert(isIn1 && isIn2); /* cerr<<"Interaction inside, weight is "<<weight<<endl;*/ /*LOG_DEBUG("Interaction inside, weight is "<<weight);*/}
100
 
                E+=geom->refLength*weight*(.5*bc->kn*pow(geom->strainN(),2)+.5*bc->ks*pow(geom->strainT().norm(),2));
101
 
        }
102
 
        return E;
103
 
}
 
79
// FIXME: rewrite for ScGeom or remove
 
80
// Real elasticEnergyInAABB(py::tuple Aabb){
 
81
//      Vector3r bbMin=py::extract<Vector3r>(Aabb[0])(), bbMax=py::extract<Vector3r>(Aabb[1])();
 
82
//      shared_ptr<Scene> rb=Omega::instance().getScene();
 
83
//      Real E=0;
 
84
//      FOREACH(const shared_ptr<Interaction>&i, *rb->interactions){
 
85
//              if(!i->phys) continue;
 
86
//              shared_ptr<NormShearPhys> bc=dynamic_pointer_cast<NormShearPhys>(i->phys); if(!bc) continue;
 
87
//              shared_ptr<Dem3DofGeom> geom=dynamic_pointer_cast<Dem3DofGeom>(i->geom); if(!geom){LOG_ERROR("NormShearPhys contact doesn't have Dem3DofGeom associated?!"); continue;}
 
88
//              const shared_ptr<Body>& b1=Body::byId(i->getId1(),rb), b2=Body::byId(i->getId2(),rb);
 
89
//              bool isIn1=isInBB(b1->state->pos,bbMin,bbMax), isIn2=isInBB(b2->state->pos,bbMin,bbMax);
 
90
//              if(!isIn1 && !isIn2) continue;
 
91
//              LOG_DEBUG("Interaction #"<<i->getId1()<<"--#"<<i->getId2());
 
92
//              Real weight=1.;
 
93
//              if((!isIn1 && isIn2) || (isIn1 && !isIn2)){
 
94
//                      //shared_ptr<Body> bIn=isIn1?b1:b2, bOut=isIn2?b2:b1;
 
95
//                      Vector3r vIn=(isIn1?b1:b2)->state->pos, vOut=(isIn2?b1:b2)->state->pos;
 
96
//                      #define _WEIGHT_COMPONENT(axis) if(vOut[axis]<bbMin[axis]) weight=min(weight,abs((vOut[axis]-bbMin[axis])/(vOut[axis]-vIn[axis]))); else if(vOut[axis]>bbMax[axis]) weight=min(weight,abs((vOut[axis]-bbMax[axis])/(vOut[axis]-vIn[axis])));
 
97
//                      _WEIGHT_COMPONENT(0); _WEIGHT_COMPONENT(1); _WEIGHT_COMPONENT(2);
 
98
//                      assert(weight>=0 && weight<=1);
 
99
//                      //cerr<<"Interaction crosses Aabb boundary, weight is "<<weight<<endl;
 
100
//                      //LOG_DEBUG("Interaction crosses Aabb boundary, weight is "<<weight);
 
101
//              } else {assert(isIn1 && isIn2); /* cerr<<"Interaction inside, weight is "<<weight<<endl;*/ /*LOG_DEBUG("Interaction inside, weight is "<<weight);*/}
 
102
//              E+=geom->refLength*weight*(.5*bc->kn*pow(geom->strainN(),2)+.5*bc->ks*pow(geom->strainT().norm(),2));
 
103
//      }
 
104
//      return E;
 
105
// }
104
106
 
105
107
/* return histogram ([bin1Min,bin2Min,…],[value1,value2,…]) from projections of interactions
106
108
 * to the plane perpendicular to axis∈[0,1,2]; the number of bins can be specified and they cover
353
355
        upper part" (lower and upper parts with respect to the plane's normal).
354
356
 
355
357
        (This could be easily extended to return sum of only normal forces or only of shear forces.)
 
358
        // FIXME: adapt to ScGeom or remove
356
359
*/
357
 
Vector3r forcesOnPlane(const Vector3r& planePt, const Vector3r&  normal){
358
 
        Vector3r ret(Vector3r::Zero());
359
 
        Scene* scene=Omega::instance().getScene().get();
360
 
        FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){
361
 
                if(!I->isReal()) continue;
362
 
                NormShearPhys* nsi=dynamic_cast<NormShearPhys*>(I->phys.get());
363
 
                if(!nsi) continue;
364
 
                Vector3r pos1,pos2;
365
 
                Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(I->geom.get()); // Dem3DofGeom has copy of se3 in itself, otherwise we have to look up the bodies
366
 
                if(d3dg){ pos1=d3dg->se31.position; pos2=d3dg->se32.position; }
367
 
                else{ pos1=Body::byId(I->getId1(),scene)->state->pos; pos2=Body::byId(I->getId2(),scene)->state->pos; }
368
 
                Real dot1=(pos1-planePt).dot(normal), dot2=(pos2-planePt).dot(normal);
369
 
                if(dot1*dot2>0) continue; // both (centers of) bodies are on the same side of the plane=> this interaction has to be disregarded
370
 
                // if pt1 is on the negative plane side, d3dg->normal.Dot(normal)>0, the force is well oriented;
371
 
                // otherwise, reverse its contribution. So that we return finally
372
 
                // Sum [ ( normal(plane) dot normal(interaction= from 1 to 2) ) "nsi->force" ]
373
 
                ret+=(dot1<0.?1.:-1.)*(nsi->normalForce+nsi->shearForce);
374
 
        }
375
 
        return ret;
376
 
}
 
360
// Vector3r forcesOnPlane(const Vector3r& planePt, const Vector3r&  normal){
 
361
//      Vector3r ret(Vector3r::Zero());
 
362
//      Scene* scene=Omega::instance().getScene().get();
 
363
//      FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){
 
364
//              if(!I->isReal()) continue;
 
365
//              NormShearPhys* nsi=dynamic_cast<NormShearPhys*>(I->phys.get());
 
366
//              if(!nsi) continue;
 
367
//              Vector3r pos1,pos2;
 
368
//              Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(I->geom.get()); // Dem3DofGeom has copy of se3 in itself, otherwise we have to look up the bodies
 
369
//              if(d3dg){ pos1=d3dg->se31.position; pos2=d3dg->se32.position; }
 
370
//              else{ pos1=Body::byId(I->getId1(),scene)->state->pos; pos2=Body::byId(I->getId2(),scene)->state->pos; }
 
371
//              Real dot1=(pos1-planePt).dot(normal), dot2=(pos2-planePt).dot(normal);
 
372
//              if(dot1*dot2>0) continue; // both (centers of) bodies are on the same side of the plane=> this interaction has to be disregarded
 
373
//              // if pt1 is on the negative plane side, d3dg->normal.Dot(normal)>0, the force is well oriented;
 
374
//              // otherwise, reverse its contribution. So that we return finally
 
375
//              // Sum [ ( normal(plane) dot normal(interaction= from 1 to 2) ) "nsi->force" ]
 
376
//              ret+=(dot1<0.?1.:-1.)*(nsi->normalForce+nsi->shearForce);
 
377
//      }
 
378
//      return ret;
 
379
// }
377
380
 
378
381
/* Less general than forcesOnPlane, computes force on plane perpendicular to axis, passing through coordinate coord. */
379
 
Vector3r forcesOnCoordPlane(Real coord, int axis){
380
 
        Vector3r planePt(Vector3r::Zero()); planePt[axis]=coord;
381
 
        Vector3r normal(Vector3r::Zero()); normal[axis]=1;
382
 
        return forcesOnPlane(planePt,normal);
383
 
}
 
382
// /*Vector3r forcesOnCoordPlane(Real coord, int axis){
 
383
//      Vector3r planePt(Vector3r::Zero()); planePt[axis]=coord;
 
384
//      Vector3r normal(Vector3r::Zero()); normal[axis]=1;
 
385
//      return forcesOnPlane(planePt,normal);
 
386
// }*/
384
387
 
385
388
 
386
389
py::tuple spiralProject(const Vector3r& pt, Real dH_dTheta, int axis=2, Real periodStart=std::numeric_limits<Real>::quiet_NaN(), Real theta0=0){
507
510
        py::def("setRefSe3",setRefSe3,"Set reference :yref:`positions<State::refPos>` and :yref:`orientations<State::refOri>` of all :yref:`bodies<Body>` equal to their current :yref:`positions<State::pos>` and :yref:`orientations<State::ori>`.");
508
511
        py::def("interactionAnglesHistogram",interactionAnglesHistogram,interactionAnglesHistogram_overloads(py::args("axis","mask","bins","aabb")));
509
512
        py::def("bodyNumInteractionsHistogram",bodyNumInteractionsHistogram,bodyNumInteractionsHistogram_overloads(py::args("aabb")));
510
 
        py::def("elasticEnergy",elasticEnergyInAABB);
 
513
//      py::def("elasticEnergy",elasticEnergyInAABB);
511
514
        py::def("inscribedCircleCenter",inscribedCircleCenter,(py::arg("v1"),py::arg("v2"),py::arg("v3")),"Return center of inscribed circle for triangle given by its vertices *v1*, *v2*, *v3*.");
512
515
        py::def("unbalancedForce",&Shop__unbalancedForce,(py::args("useMaxForce")=false),"Compute the ratio of mean (or maximum, if *useMaxForce*) summary force on bodies and mean force magnitude on interactions. For perfectly static equilibrium, summary force on all bodies is zero (since forces from interactions cancel out and induce no acceleration of particles); this ratio will tend to zero as simulation stabilizes, though zero is never reached because of finite precision computation. Sufficiently small value can be e.g. 1e-2 or smaller, depending on how much equilibrium it should be.");
513
516
        py::def("kineticEnergy",Shop__kineticEnergy,(py::args("findMaxId")=false),"Compute overall kinetic energy of the simulation as\n\n.. math:: \\sum\\frac{1}{2}\\left(m_i\\vec{v}_i^2+\\vec{\\omega}(\\mat{I}\\vec{\\omega}^T)\\right).\n\nFor :yref:`aspherical<Body.aspherical>` bodies, the inertia tensor $\\mat{I}$ is transformed to global frame, before multiplied by $\\vec{\\omega}$, therefore the value should be accurate.\n");
514
517
        py::def("sumForces",sumForces,(py::arg("ids"),py::arg("direction")),"Return summary force on bodies with given *ids*, projected on the *direction* vector.");
515
518
        py::def("sumTorques",sumTorques,(py::arg("ids"),py::arg("axis"),py::arg("axisPt")),"Sum forces and torques on bodies given in *ids* with respect to axis specified by a point *axisPt* and its direction *axis*.");
516
519
        py::def("sumFacetNormalForces",sumFacetNormalForces,(py::arg("ids"),py::arg("axis")=-1),"Sum force magnitudes on given bodies (must have :yref:`shape<Body.shape>` of the :yref:`Facet` type), considering only part of forces perpendicular to each :yref:`facet's<Facet>` face; if *axis* has positive value, then the specified axis (0=x, 1=y, 2=z) will be used instead of facet's normals.");
517
 
        py::def("forcesOnPlane",forcesOnPlane,(py::arg("planePt"),py::arg("normal")),"Find all interactions deriving from :yref:`NormShearPhys` that cross given plane and sum forces (both normal and shear) on them.\n\n:param Vector3 planePt: a point on the plane\n:param Vector3 normal: plane normal (will be normalized).\n");
518
 
        py::def("forcesOnCoordPlane",forcesOnCoordPlane);
 
520
//      py::def("forcesOnPlane",forcesOnPlane,(py::arg("planePt"),py::arg("normal")),"Find all interactions deriving from :yref:`NormShearPhys` that cross given plane and sum forces (both normal and shear) on them.\n\n:param Vector3 planePt: a point on the plane\n:param Vector3 normal: plane normal (will be normalized).\n");
 
521
//      py::def("forcesOnCoordPlane",forcesOnCoordPlane);
519
522
        py::def("totalForceInVolume",Shop__totalForceInVolume,"Return summed forces on all interactions and average isotropic stiffness, as tuple (Vector3,float)");
520
523
        py::def("createInteraction",Shop__createExplicitInteraction,(py::arg("id1"),py::arg("id2")),"Create interaction between given bodies by hand.\n\nCurrent engines are searched for :yref:`IGeomDispatcher` and :yref:`IPhysDispatcher` (might be both hidden in :yref:`InteractionLoop`). Geometry is created using ``force`` parameter of the :yref:`geometry dispatcher<IGeomDispatcher>`, wherefore the interaction will exist even if bodies do not spatially overlap and the functor would return ``false`` under normal circumstances. \n\n.. warning:: This function will very likely behave incorrectly for periodic simulations (though it could be extended it to handle it farily easily).");
521
524
        py::def("spiralProject",spiralProject,(py::arg("pt"),py::arg("dH_dTheta"),py::arg("axis")=2,py::arg("periodStart")=std::numeric_limits<Real>::quiet_NaN(),py::arg("theta0")=0));
526
529
        py::def("wireNone",wireNone,"Set :yref:`Shape::wire` on all bodies to False, rendering them as solids.");
527
530
        py::def("wireNoSpheres",wireNoSpheres,"Set :yref:`Shape::wire` to True on non-spherical bodies (:yref:`Facets<Facet>`, :yref:`Walls<Wall>`).");
528
531
        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:: This function is currently broken and should not be used.");
529
 
        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.. math:: k_n=\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_n)^2\\right) \\\\ c_n=-\\frac{2m}{t_c} \\\\  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$).");
 
532
        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.. 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$).");
530
533
        py::def("stressTensorOfPeriodicCell",Shop::getStress,(py::args("volume")=0),"Deprecated, use utils.getStress instead |ydeprecated|");
531
534
        //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");
532
535
        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).");
533
536
        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.");
534
537
        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}$.");
535
538
        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 l_i^b f_j^b$, where the sum is over all interactions, with $l$ the branch vector (joining centers of the bodies) and $f$ is the contact force. $V$ can be passed to the function. If it is not, it will be equal to one in non-periodic cases, or equal to the volume of the cell in periodic cases.");
 
539
        py::def("getCapillaryStress",Shop::getCapillaryStress,(py::args("volume")=0),"Compute and return Love-Weber capillary stress tensor:\n\n $\\sigma^{cap}_{ij}=\\frac{1}{V}\\sum_b l_i^b f^{cap,b}_j$, where the sum is over all interactions, with $l$ the branch vector (joining centers of the bodies) and $f^{cap}$ is the capillary force. $V$ can be passed to the function. If it is not, it will be equal to one in non-periodic cases, or equal to the volume of the cell in periodic cases. Only the CapillaryPhys interaction type is supported presently.");
536
540
        py::def("getBodyIdsContacts",Shop__getBodyIdsContacts,(py::args("bodyID")=0),"Get a list of body-ids, which contacts the given body.");
537
541
        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.");
538
542
        py::def("shiftBodies",shiftBodies,(py::arg("ids"),py::arg("shift")),"Shifts bodies listed in ids without updating their velocities.");
539
543
        py::def("calm",Shop__calm,(py::arg("mask")=-1),"Set translational and rotational velocities of all bodies to zero.");
540
544
        py::def("setNewVerticesOfFacet",setNewVerticesOfFacet,(py::arg("b"),py::arg("v1"),py::arg("v2"),py::arg("v3")),"Sets new vertices (in global coordinates) to given facet.");
541
545
        py::def("setContactFriction",Shop::setContactFriction,py::arg("angleRad"),"Modify the friction angle (in radians) inside the material classes and existing contacts. The friction for non-dynamic bodies is not modified.");
542
 
        py::def("growParticles",Shop::growParticles,(py::args("multiplier"), py::args("updateMass")=true, py::args("dynamicOnly")=true), "Change the size of spheres and sphere clumps by the multiplier. If updateMass=True, then the mass is updated. dynamicOnly=True is mandatory in many cases since in current implementation the function would crash on the non-spherical and non-dynamic bodies (e.g. facets, boxes, etc.)");
 
546
        py::def("growParticles",Shop::growParticles,(py::args("multiplier"), py::args("updateMass")=true, py::args("dynamicOnly")=true, py::args("discretization")=15, py::args("integrateInertia")=true), "Change the size of spheres and sphere clumps by the multiplier. If updateMass=True, then the mass is updated. dynamicOnly=True is mandatory in many cases since in current implementation the function would crash on the non-spherical and non-dynamic bodies (e.g. facets, boxes, etc.). For clumps the masses and inertias are adapted automatically (for details of inertia tensor integration scheme see :yref:`clump()<BodyContainer.clump>`).");
543
547
        py::def("intrsOfEachBody",intrsOfEachBody,"returns list of lists of interactions of each body");
544
548
        py::def("numIntrsOfEachBody",numIntrsOfEachBody,"returns list of number of interactions of each body");
 
549
        py::def("TetrahedronSignedVolume",static_cast<Real (*)(const vector<Vector3r>&)>(&TetrahedronSignedVolume),"TODO");
 
550
        py::def("TetrahedronVolume",static_cast<Real (*)(const vector<Vector3r>&)>(&TetrahedronVolume),"TODO");
 
551
        py::def("TetrahedronInertiaTensor",TetrahedronInertiaTensor,"TODO");
 
552
        py::def("TetrahedronCentralInertiaTensor",TetrahedronCentralInertiaTensor,"TODO");
 
553
        py::def("TetrahedronWithLocalAxesPrincipal",TetrahedronWithLocalAxesPrincipal,"TODO");
 
554
        py::def("momentum",Shop::momentum,"TODO");
 
555
        py::def("angularMomentum",Shop::angularMomentum,(py::args("origin")=Vector3r(Vector3r::Zero())),"TODO");
545
556
}