yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12981
[Branch ~yade-pkg/yade/git-trunk] Rev 3998: fix the flipCell() function
------------------------------------------------------------
revno: 3998
committer: bchareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
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).");
