------------------------------------------------------------ revno: 3998 committer: bchareyre <bruno.chare...@grenoble-inp.fr> timestamp: Mon 2017-02-06 18:54:41 +0100 message: fix the flipCell() function modified: pkg/dem/Shop_01.cpp py/_utils.cpp
-- lp:yade https://code.launchpad.net/~yade-pkg/yade/git-trunk Your team Yade developers is subscribed to branch lp:yade. To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'pkg/dem/Shop_01.cpp' --- pkg/dem/Shop_01.cpp 2016-05-26 16:30:41 +0000 +++ pkg/dem/Shop_01.cpp 2017-02-06 17:54:41 +0000 @@ -48,63 +48,29 @@ CREATE_LOGGER(Shop); -/*! Flip periodic cell by given number of cells. - -Still broken, some interactions are missed. Should be checked. -*/ - +/*! Flip periodic cell for shearing indefinitely.*/ Matrix3r Shop::flipCell(const Matrix3r& _flip){ - Scene* scene=Omega::instance().getScene().get(); const shared_ptr<Cell>& cell(scene->cell); const Matrix3r& trsf(cell->trsf); - Vector3r size=cell->getSize(); - Matrix3r flip; + Scene* scene=Omega::instance().getScene().get(); const shared_ptr<Cell>& cell(scene->cell); + Matrix3r& hSize = cell->hSize; + Matrix3i flip; if(_flip==Matrix3r::Zero()){ bool hasNonzero=false; for(int i=0; i<3; i++) for(int j=0; j<3; j++) { if(i==j){ flip(i,j)=0; continue; } - flip(i,j)=-floor(.5+trsf(i,j)/(size[j]/size[i])); + flip(i,j)=-floor(hSize.col(j).dot(hSize.col(i))/hSize.col(i).dot(hSize.col(i))); if(flip(i,j)!=0) hasNonzero=true; } if(!hasNonzero) {LOG_TRACE("No flip necessary."); return Matrix3r::Zero();} - LOG_DEBUG("Computed flip matrix: upper "<<flip(0,1)<<","<<flip(0,2)<<","<<flip(1,2)<<"; lower "<<flip(1,0)<<","<<flip(2,0)<<","<<flip(2,1)); } else { - flip=_flip; - } - - // current cell coords of bodies - vector<Vector3i > oldCells; oldCells.resize(scene->bodies->size()); - FOREACH(const shared_ptr<Body>& b, *scene->bodies){ - if(!b) continue; cell->wrapShearedPt(b->state->pos,oldCells[b->getId()]); - } - - // change cell trsf here - Matrix3r trsfInc; - for(int i=0; i<3; i++) for(int j=0; j<3; j++){ - if(i==j) { if(flip(i,j)!=0) LOG_WARN("Non-zero diagonal term at ["<<i<<","<<j<<"] is meaningless and will be ignored."); trsfInc(i,j)=0; continue; } - // make sure non-diagonal entries are "integers" - if(flip(i,j)!=double(int(flip(i,j)))) LOG_WARN("Flip matrix entry "<<flip(i,j)<<" at ["<<i<<","<<j<<"] not integer?! (will be rounded)"); - trsfInc(i,j)=int(flip(i,j))*size[j]/size[i]; - } - cell->trsf+=trsfInc; + flip=_flip.cast<int>(); + } + Matrix3r flipFloat = flip.cast<Real>(); + cell->hSize+=cell->hSize*flipFloat; cell->postLoad(*cell); - // new cell coords of bodies - vector<Vector3i > newCells; newCells.resize(scene->bodies->size()); - FOREACH(const shared_ptr<Body>& b, *scene->bodies){ - if(!b) continue; - cell->wrapShearedPt(b->state->pos,newCells[b->getId()]); - } - - // remove all potential interactions - scene->interactions->eraseNonReal(); - // adjust Interaction::cellDist for real interactions; - FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){ - Body::id_t id1=i->getId1(),id2=i->getId2(); - // this must be the same for both old and new interaction: cell2-cell1+cellDist - // c₂-c₁+c₁₂=k; c'₂+c₁'+c₁₂'=k (new cell coords have ') - // c₁₂'=(c₂-c₁+c₁₂)-(c₂'-c₁') - i->cellDist=(oldCells[id2]-oldCells[id1]+i->cellDist)-(newCells[id2]-newCells[id1]); - } - + // adjust Interaction::cellDist for interactions; + Matrix3r invFlip = (Matrix3r::Identity() + flipFloat).inverse(); + FOREACH(const shared_ptr<Interaction>& i, *scene->interactions) i->cellDist = invFlip*i->cellDist; // force reinitialization of the collider bool colliderFound=false; @@ -113,7 +79,7 @@ if(c){ colliderFound=true; c->invalidatePersistentData(); } } if(!colliderFound) LOG_WARN("No collider found while flipping cell; continuing simulation might give garbage results."); - return flip; + return flipFloat; } /* Apply force on contact point to 2 bodies; the force is oriented as it applies on the first body and is reversed on the second. === modified file 'py/_utils.cpp' --- py/_utils.cpp 2016-10-18 13:46:31 +0000 +++ py/_utils.cpp 2017-02-06 17:54:41 +0000 @@ -487,7 +487,7 @@ py::def("wireAll",wireAll,"Set :yref:`Shape::wire` on all bodies to True, rendering them with wireframe only."); py::def("wireNone",wireNone,"Set :yref:`Shape::wire` on all bodies to False, rendering them as solids."); 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:: This function is currently broken and should not be used."); + 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``."); 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).");
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp