Yade

O.cell.volume=0 during flipCell

Asked by Sacha Duverger

Hello,

I am trying to use flipCell to avoid a RuntimeError "Body larger than half of the cell size encountered" during some periodic simulations (see my previous question about the documentation [1]). I am having a more technical issue now, it appears that in some simulations the new O.cell.hSize computed in flipCell has a determinant equal to 0.

A bit of context:
The periodic DEM simulations I am running are actually RVEs for a larger scale method (the Material Point Method), I have thus no direct control on the velocity gradient imposed to the DEM periodic cells. During a MPMxDEM granular column collapse, some RVEs get deformed a lot in unpredictable directions so I’d like to flip them whenever the cell is too much deformed. I tried to generalise the flipping condition used in [2] so it takes into account all directions, see the MWE below.

A MWE with a problematic RVE’s hSize:
```
O.periodic = True

O.cell.hSize = Matrix3(0.001767398585873738732,-3.402651098340004036e-22,-2.455726386595863492e-23, 2.840638172828420518e-23,0.004990738513508447341,-0.0009967331209356602214, -1.550207193852494163e-22,-0.001442690301487606743,0.00144229021430840711)

cond_x = abs(O.cell.hSize[0,1])>abs(O.cell.hSize[0,0]) or abs(O.cell.hSize[0,2])>abs(O.cell.hSize[0,0])
cond_y = abs(O.cell.hSize[1,0])>abs(O.cell.hSize[1,1]) or abs(O.cell.hSize[1,2])>abs(O.cell.hSize[1,1])
cond_z = abs(O.cell.hSize[2,0])>abs(O.cell.hSize[2,2]) or abs(O.cell.hSize[2,1])>abs(O.cell.hSize[2,2])

print("\nInitial hSize:", O.cell.hSize, "\nFlipping condition: ", cond_x or cond_y or cond_z, “\n”)

try: flipCell()
except RuntimeError as err:
    print("\nNew hSize: ", O.cell.hSize, "\n")
    raise(err)
```

Which outputs with yadedaily (git revision 406fb7a):
```
Welcome to Yade 20220206-6304~406fb7a~focal1
Using python version: 3.8.10 (default, Mar 15 2022, 12:22:08)
[GCC 9.4.0]
Warning: no X rendering available (see https://bbs.archlinux.org/viewtopic.php?id=13189)
TCP python prompt on localhost:9001, auth cookie `duacks'
XMLRPC info provider on http://localhost:21001
/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py:935: UserWarning: Attempting to work in a virtualenv. If you encounter problems, please install IPython inside the virtualenv.
  warn("Attempting to work in a virtualenv. If you encounter problems, please "
Running script mwe.py

Initial hSize: Matrix3(0.001767398585873738732,-3.402651098340004036e-22,-2.455726386595863492e-23, 2.840638172828420518e-23,0.004990738513508447341,-0.0009967331209356602214, -1.550207193852494163e-22,-0.001442690301487606743,0.00144229021430840711)
Flipping condition: True

New hSize: Matrix3(0.001767398585873738732,0.001767398585873738732,0.001767398585873738732, 0.003994005392572786903,0.002000539150701466894,0.003994005392572786903, -4.000871791996330812e-07,0.002884180341437615237,-4.000871791996329224e-07)

Traceback (most recent call last):
  File "/usr/bin/yadedaily", line 343, in runScript
    execfile(script,globals())
  File "/usr/lib/python3/dist-packages/past/builtins/misc.py", line 87, in execfile
    exec_(code, myglobals, mylocals)
  File "mwe.py", line 13, in <module>
    raise(err)
  File "mwe.py", line 10, in <module>
    try: flipCell()
RuntimeError: Cell is degenerate (zero volume)
```

This error is raised because the determinant of O.cell.hSize is zero [3], which is tested because Cell::integrateAndUpdate is called in cell::postLoad [4], which is called right after changing cell::hSize in flipCell [5].

Should flipCell be able to find a flip matrix for any hSize ? Is there a problem the way I use it ?

Thanks in advance,

Sacha

[1]: https://answers.launchpad.net/yade/+question/701764
[2]: https://gitlab.com/yade-dev/trunk/-/blob/277c4f0c4a94dd1df423e8d499e56a457d356ea0/examples/PeriodicBoundaries/cellFlipping.py
[3]: https://gitlab.com/yade-dev/trunk/-/blob/277c4f0c4a94dd1df423e8d499e56a457d356ea0/core/Cell.cpp#L19
[4]: https://gitlab.com/yade-dev/trunk/-/blob/277c4f0c4a94dd1df423e8d499e56a457d356ea0/core/Cell.hpp#L178
[5]: https://gitlab.com/yade-dev/trunk/-/blob/277c4f0c4a94dd1df423e8d499e56a457d356ea0/preprocessing/dem/Shop_01.cpp#L84

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
Last query:
Last reply:
Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#1

In principle all cells can be flipped, yes. What you get suggests that the automatic determination of the flip matrix fails.
Not sure why right now. You could pass one flip matrix explicitly.
Is cell.velGrad always symmetric or arbitrary (do you have rotation)?

Thanks for posting example numbers, it might help investigating.

Bruno

Revision history for this message
Sacha Duverger (schmxprr) said : 		#2

Thank you for your answer, yes cell.velGrad is always symmetric in my simulations.

The automatic determination of the flip matrix was very convenient for me since I don't really understand yet how it is computed. From what I understand of [6]:
 - all diagonal coefficients a_ii are always zero
 - non diagonal coefficient a_ij are a rounded ratio between the scalar product of the i^th and j^th base cell vectors and the squared norm of the i^th base cell vector

So I guess a_ij is the cosine of a rotation angle, but not quite the angle between the i^th and j^th base cell vectors because they don't necessarily have the same norm. Also there is this round operation: I don't understand why it is performed, and why in two steps (floor + int).

Is there a simpler expression of the flip matrix when cell.hSize has always been deformed with a symmetric cell.velGrad ?

Thanks,

Sacha

[6]: https://gitlab.com/yade-dev/trunk/-/blob/277c4f0c4a94dd1df423e8d499e56a457d356ea0/preprocessing/dem/Shop_01.cpp#L66-72

Revision history for this message
Sacha Duverger (schmxprr) said (last edit ): 		#3

I tried something that seems to work: I changed the way flipCell computes hSize, it uses sums between the base cell vectors to flip the cell towards the most favorable neighbour grid point, see patch [7]. So my version of flipCell doesn't use a flip matrix anymore, and it re-computes interaction.cellDist the same way it is done in the Dispatcher.
With the MWE's hSize in my original question, this version of flipCell effectively find the flip that makes the second base cell vector almost axis aligned.

My RVEs are not excessively sheared anymore so the MPMxDEM simulation goes a little bit further, however they are now excessively stretched: one of the base cell vector has a norm more than 6 times larger than the norm of the 2 others, so I still get "Body larger than half of the cell size encountered".
I'am afraid flipCell can't do anything for that, right ?

Cheers,

Sacha

[7]: patch on preprocessing/dem/Shop_01.cpp and preprocessing/dem/Shop.hpp with respect to git revision 406fb7afb768bea4bd01def605fdac82f58122df

```
diff --git a/preprocessing/dem/Shop.hpp b/preprocessing/dem/Shop.hpp
old mode 100644
new mode 100755
index 6380ed7e5..a4624c589
--- a/preprocessing/dem/Shop.hpp
+++ b/preprocessing/dem/Shop.hpp
@@ -11,6 +11,8 @@
 #include <lib/base/AliasNamespaces.hpp>
 #include <boost/function.hpp>

+#include <cmath>
+
 namespace yade { // Cannot have #include directive inside.

 class Scene;
diff --git a/preprocessing/dem/Shop_01.cpp b/preprocessing/dem/Shop_01.cpp
old mode 100644
new mode 100755
index dc95f389b..3635e9263
--- a/preprocessing/dem/Shop_01.cpp
+++ b/preprocessing/dem/Shop_01.cpp
@@ -59,39 +59,49 @@ Matrix3r Shop::flipCell(const 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) = -int(math::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();
- }
- } else {
- if ((_flip + Matrix3r::Identity()).determinant() != 1) LOG_WARN("Flipping cell needs det(Id+flip)=1, check your input.");
- flip = _flip.cast<int>();
+ Matrix3r& new_hSize = cell->hSize;
+ Matrix3r flip;
+ int j,k;
+ double phi_1, phi_2a, phi_2b, theta;
+ Vector3r vect_a, vect_b;
+ shared_ptr<Body> b1, b2;
+ std::vector<std::tuple<int, int>> planes = {{0,1}, {0,2}, {1,2}, {1,0}, {2,0}, {2,1}};
+ for (auto plane : planes) {
+ k = std::get<0>(plane);
+ j = std::get<1>(plane);
+
+ // Get the angle between `the projection of the ith base cell vector on the current plane` and `the ith axis`
+ phi_1 = acos(cell->hSize(k,k) / pow((pow(cell->hSize(k,k), 2)+pow(cell->hSize(j,k), 2)), .5) );
+
+ // Compute the angles if we flip to neighboor grid points
+ // If we flip left
+ vect_a = cell->hSize.col(k) - cell->hSize.col(j);
+ phi_2a = acos(vect_a(k) / pow((pow(vect_a(k), 2)+pow(vect_a(j), 2)), .5));
+
+ // If we flip right
+ vect_b = cell->hSize.col(k) + cell->hSize.col(j);
+ phi_2b = acos(vect_b(k) / pow((pow(vect_b(k), 2)+pow(vect_b(j), 2)), .5));
+
+ // Keep the best angle (the smallest)
+ theta = std::min({phi_1, phi_2a, phi_2b});
+
+ // Flip in the best direction (if there is a grid point which gives a lower angle)
+ if (theta==phi_2a) new_hSize.col(k) = vect_a;
+ else if (theta==phi_2b) new_hSize.col(k) = vect_b;
  }
- cell->hSize += cell->hSize * flip.cast<Real>();
+
+ flip = _flip; // Just so -Werror=unused-parameter doesn't bother me for now (to be improved)
+ cell->hSize = new_hSize;
+
  cell->postLoad(*cell);

- // adjust Interaction::cellDist for interactions;
- // adjunct matrix of (Id + flip) is the inverse since det=1, below is the transposed co-factor matrix of (Id+flip).
- // note that Matrix3::adjoint is not the adjunct, hence the in-place adjunct below
- Matrix3i invFlip;
- invFlip << 1 - flip(2, 1) * flip(1, 2), flip(2, 1) * flip(0, 2) - flip(0, 1), flip(0, 1) * flip(1, 2) - flip(0, 2),
- flip(1, 2) * flip(2, 0) - flip(1, 0), 1 - flip(0, 2) * flip(2, 0), flip(0, 2) * flip(1, 0) - flip(1, 2), flip(1, 0) * flip(2, 1) - flip(2, 0),
- flip(2, 0) * flip(0, 1) - flip(2, 1), 1 - flip(1, 0) * flip(0, 1);
- for (const auto& i : *scene->interactions)
- i->cellDist = invFlip * i->cellDist;
+ // Recompute Interaction::cellDist for interactions;
+ for (const auto& I : *scene->interactions) {
+ b1 = Body::byId(I->getId1(), scene);
+ b2 = Body::byId(I->getId2(), scene);
+ for (int i = 0; i < 3; i++)
+ I->cellDist[i] = -(int)((b2->state->pos[i] - b1->state->pos[i]) / cell->getSize()[i] + .5);
+ }

  // force reinitialization of the collider
  bool colliderFound = false;

```

Revision history for this message
Sacha Duverger (schmxprr) said : 		#4

An update:

I worked around the excessive stretch problem by changing the initial aspect ratio of my RVEs, since I know in which direction they will be stretched.

As I am pretty confident with my version of flipCell, I pushed my branch on YADE's repo [8] (it includes a new example that illustrate my changes [9]).

On this branch, the usage of flipCell in [2] is still valid. However, the construction of the input flip matrix is different (or at least I think so, since I don't quite understand how it should have been constructed before).

Cheers,

Sacha

[8]: https://gitlab.com/yade-dev/trunk/-/tree/change/cellFlip
[9]: https://gitlab.com/yade-dev/trunk/-/blob/change/cellFlip/examples/PeriodicBoundaries/cellFlip_possibilities.py

Revision history for this message
Jérôme Duriez (jduriez) said : 		#5

Note that Sacha's new branch "cellFlip" (thank you Sacha !) might be a bit complicated to compare with current master, at the moment. We will simplify things and try to propose a proper MR at some point.

Revision history for this message
Bruno Chareyre (bruno-chareyre) said (last edit ): 		#6

Hi,

I'll check the suggested MR whenever I have a chance.

The case of stretching can be overcomed in a few special cases. It needs to define an initial cell size not aligned with the principal strain axis. It has limitations and it will probably not work for non-monotonic multi-directional deformation.

See more here: https://www.sciencedirect.com/science/article/abs/pii/030193229290074Q

B

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#7

See more here: https://www.sciencedirect.com/science/article/abs/pii/030193229290074Q

Revision history for this message
Best Bruno Chareyre (bruno-chareyre) said : 		#8

At the moment it is a bit convolved to read your changes but it seems I was able to break your algorithm.
A suspicious feature of yours is that it tries to align the period on the reference frame, which isn't helping in general.
Try a rigid body rotation (velGrad = [[0,1,0], [-1,0,0], [0,0,0]] ) and you will see the problem:
- older version leads to zero volume (wrong)
- new version does flip without "zero volume" error, but it shouldn't: it turns the cube into a more elongated thing, the opposite of what we want; it also breaks my contact model (most likely a problem in updating cellDist).
My impression is that there must be is a much simpler fix to the issue.
I'll need to cleam my test script before showing (right now it's a mess).
Let's continue the discussion on gitlab?

B

Revision history for this message
Sacha Duverger (schmxprr) said : 		#9

Ah, right, if there is rotations the grid points that give lower angles between the cell's base vectors and the global axes are not necessarily the ones that make the cell's base vectors orthogonal between each others.
I guess my version of flipCell could work if it computed the angles with respect to the rotated global axes (rotated with the non-symmetric part of O.cell. trsf).

Anyway, my version of flipCell performs many tests, so a fixed version of the original would certainly be better.

I created and pushed a new branch (change/flipCell) from master where I cherry-picked all my changes on flipCell [10], so I deleted my previous branch (change/cellFlip).

Thank you for your help,

Sacha

[10]: https://gitlab.com/yade-dev/trunk/-/compare/master...change%2FflipCell?from_project_id=10133144

Revision history for this message
Sacha Duverger (schmxprr) said : 		#10

Thanks Bruno Chareyre, that solved my question.