Relax projection for FSI mapping to account for curvature (#1223)

* Relax projection for FSI mapping to account for curvature * Style * Format again
Exawind · Oct 18, 2023 · 5b978b2 · 5b978b2
1 parent 174f2f0
commit 5b978b2
Show file tree

Hide file tree

Showing 2 changed files with 52 additions and 0 deletions.
diff --git a/src/aero/fsi/FSIturbine.C b/src/aero/fsi/FSIturbine.C
@@ -1765,6 +1765,8 @@ fsiTurbine::computeMapping()
         vs::Vector ptCoords(xyz[0], xyz[1], xyz[2]);
         bool foundProj = false;
         double nDimCoord = -1.0;
+        double minDispMapInterp = 1.0e6;
+        int minDispMap = 1e6;
         for (int i = 0; i < nPtsBlade - 1; i++) {
           vs::Vector lStart = {
             brFSIdata_.bld_ref_pos[(iStart + i) * 6],
@@ -1776,6 +1778,11 @@ fsiTurbine::computeMapping()
             brFSIdata_.bld_ref_pos[(iStart + i + 1) * 6 + 2]};
           nDimCoord = fsi::projectPt2Line(ptCoords, lStart, lEnd);
 
+          if (std::abs(nDimCoord) < minDispMapInterp) {
+            minDispMapInterp = std::abs(nDimCoord);
+            minDispMap = i;
+          }
+
           if ((nDimCoord >= 0) && (nDimCoord <= 1.0)) {
             foundProj = true;
             *dispMapInterpNode = nDimCoord;
@@ -1785,6 +1792,14 @@ fsiTurbine::computeMapping()
           }
         }
 
+        // if we are very very close to a point then we need to use it
+        // curvature issues can break the projection
+        if (!foundProj && minDispMapInterp < 0.50) {
+          *dispMapInterpNode = 0.0;
+          *dispMapNode = minDispMap;
+          foundProj = true;
+        }
+
         // If no element in the OpenFAST mesh contains this node do some sanity
         // check on the perpendicular distance between the surface mesh node and
         // the line joining the ends of the blade

diff --git a/unit_tests/aero/UnitTestPt2Line.C b/unit_tests/aero/UnitTestPt2Line.C
@@ -7,6 +7,7 @@
 // for more details.
 //
 
+#include "vs/trig_ops.h"
 #include <gtest/gtest.h>
 #include <KokkosInterface.h>
 #include <aero/aero_utils/Pt2Line.h>
@@ -94,4 +95,40 @@ test_perpProjectDist_Pt2Line()
 
 TEST(aero_utils, perpProjectDist_Pt2Line) { test_perpProjectDist_Pt2Line(); }
 
+void
+test_projectPt2Line_relative(
+  const vs::Vector& lStart,
+  const vs::Vector& lEnd,
+  const vs::Vector& pt,
+  std::function<bool(double)> checker)
+{
+  double result = call_projectPt2Line_on_device(lStart, lEnd, pt);
+  auto w1 = lEnd - lStart;
+  auto w2 = pt - lStart;
+  auto angle = utils::degrees(vs::angle(w1, w2));
+  std::cerr << "ANGLE: " << angle << std::endl;
+  EXPECT_TRUE(checker(result));
+}
+
+TEST(aero_utils, projectPt2Line_corner_cases)
+{
+  // case 1
+  {
+    vs::Vector pt(-5.82431, 4.75596, 69.9067);
+    vs::Vector lStart(-2.94954, 0.177514, 69.7348);
+    vs::Vector lEnd(-2.87053, 0.189687, 70.7317);
+    auto checker = [&](double nDimCoord) { return nDimCoord > 0.0; };
+    test_projectPt2Line_relative(lStart, lEnd, pt, checker);
+  }
+
+  // case 2
+  {
+    vs::Vector pt(-0.570796, -4.73676, 80.6121);
+    vs::Vector lStart(-2.23001, 0.20372, 80.7139);
+    vs::Vector lEnd(-2.18333, 0.199448, 81.7136);
+    auto checker = [&](double nDimCoord) { return nDimCoord > 0.0; };
+    test_projectPt2Line_relative(lStart, lEnd, pt, checker);
+  }
+}
+
 } // namespace