diff --git a/src/aero/fsi/FSIturbine.C b/src/aero/fsi/FSIturbine.C index 775880c345..6f0ec984af 100644 --- 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 index faaa34ba31..38f646dbac 100644 --- 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 #include #include @@ -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 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