yade-users team mailing list archive

[Sacha Duverger <question701782@xxxxxxxxxxxxxxxxxxxxx>]

Question #701782 on Yade changed:
https://answers.launchpad.net/yade/+question/701782

Sacha Duverger posted a new comment:
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, but it still computes it for the interaction.cellDist adjustment. 
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..65bc35d4c
--- a/preprocessing/dem/Shop_01.cpp
+++ b/preprocessing/dem/Shop_01.cpp
@@ -59,39 +59,52 @@ 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;
+	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>();
+
+	cell->hSize = new_hSize;
+
+	flip = _flip; // Just so -Werror=unused-parameter doesn't show up (to be improved)
+	flip = cell->hSize.inverse() * new_hSize -  Matrix3r::Identity();
+
 	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;
+	Matrix3r 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;
+		i->cellDist = invFlip.cast<int>() * i->cellDist;
 
 	// force reinitialization of the collider
 	bool colliderFound = false;

```

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.