From c27dbdef3005bdf4626694b7100034e94cb3530e Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Thu, 31 Oct 2024 15:49:41 +0100
Subject: [PATCH] Fix growth vector calculation (do more vector reduction to
 avoid strange results from inverting badly conditioned matrix)

---
 libsrc/meshing/boundarylayer.cpp | 25 ++++++++++++++-----------
 1 file changed, 14 insertions(+), 11 deletions(-)

diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index 5e32b034..8bdc44aa 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -63,22 +63,25 @@ Vec<3> CalcGrowthVector (FlatArray<Vec<3>> ns)
         for (auto j : Range(3))
           mat(i, j) = ns[i][j];
 
-      if (fabs(mat.Det()) > 1e-6)
+      if (fabs(mat.Det()) > 1e-2)
         {
           DenseMatrix mat(3, 3);
           for (auto i : Range(3))
             for (auto j : Range(3))
               mat(i, j) = ns[i] * ns[j];
-          Vector rhs(3);
-          rhs = 1.;
-          Vector res(3);
-          DenseMatrix inv(3, ns.Size());
-          CalcInverse(mat, inv);
-          inv.Mult(rhs, res);
-          Vec<3> growth = 0.;
-          for (auto i : Range(ns))
-            growth += res[i] * ns[i];
-          return growth;
+          if (fabs(mat.Det()) > 1e-2)
+            {
+              Vector rhs(3);
+              rhs = 1.;
+              Vector res(3);
+              DenseMatrix inv(3, ns.Size());
+              CalcInverse(mat, inv);
+              inv.Mult(rhs, res);
+              Vec<3> growth = 0.;
+              for (auto i : Range(ns))
+                growth += res[i] * ns[i];
+              return growth;
+            }
         }
     }
   auto [maxpos1, maxpos2] = FindCloseVectors(ns);