Fix Floating Point Precision Issue in Lifting Surface Knot Symmetrization

pygeo/pyGeo.py

@@ -291,6 +291,8 @@ def _init_lifting_surface(
         the most commonly used options.
         """
 
+        KNOT_TOL = 1e-12
+
         if X is not None:
             Xsec = np.array(X)
         else:
@@ -394,7 +396,7 @@ def _init_lifting_surface(
                     toInsert = []
                     # Now go over the indices and see if we need to add
                     for j in range(len(indices)):
-                        if abs(baseKnots[indices[j]] - knots[j]) > 1e-12:
+                        if abs(baseKnots[indices[j]] - knots[j]) > KNOT_TOL:
                             toInsert.append(knots[j])
 
                     # Finally add the new indices and resort
@@ -411,7 +413,7 @@ def _init_lifting_surface(
             while i < len(baseKnots):
                 curKnot = baseKnots[i]
                 j = 1
-                while i + j < Nmax and abs(baseKnots[i + j] - curKnot) < 1e-12:
+                while i + j < Nmax and abs(baseKnots[i + j] - curKnot) < KNOT_TOL:
                     j += 1
                 i += j
                 newKnots.append(curKnot)
@@ -423,7 +425,12 @@ def _init_lifting_surface(
             for i in range(N):
                 if curves[i] is not None:
                     for j in range(len(newKnots)):
-                        if newKnots[j] not in curves[i].t:
+                        found = False
+                        for k_knot in curves[i].t:
+                            if abs(k_knot - newKnots[j]) < KNOT_TOL:
+                                found = True
+                                break
+                        if not found:
                             curves[i].insertKnot(newKnots[j], mult[j])
 
             # If we want a pinched tip will will zero everything here.
@@ -502,12 +509,11 @@ def _init_lifting_surface(
         knotsBot = botCurves[0].t.copy()
 
         print("Symmetrizing Knot Vectors ...")
-        eps = 1e-12
         for i in range(len(knotsTop)):
-            # Check if knotsTop[i] is not in knots_bot to within eps
+            # Check if knotsTop[i] is not in knots_bot to within KNOT_TOL
             found = False
             for j in range(len(knotsBot)):
-                if abs(knotsTop[i] - knotsBot[j]) < eps:
+                if abs(knotsTop[i] - knotsBot[j]) < KNOT_TOL:
                     found = True
 
             if not found:
@@ -516,10 +522,10 @@ def _init_lifting_surface(
                     botCurves[ii].insertKnot(knotsTop[i], 1)
 
         for i in range(len(knotsBot)):
-            # Check if knotsBot[i] is not in knotsTop to within eps
+            # Check if knotsBot[i] is not in knotsTop to within KNOT_TOL
             found = False
             for j in range(len(knotsTop)):
-                if abs(knotsBot[i] - knotsTop[j]) < eps:
+                if abs(knotsBot[i] - knotsTop[j]) < KNOT_TOL:
                     found = True
 
             if not found:

tests/reg_tests/ref/test_pyGeo.ref

@@ -1,3 +1,4 @@
 {
-    "sum of surface data": 339809.4700319107
+    "clash sum of surface data": 134260.6609166085,
+    "original sum of surface data": 339809.4700319107
 }

tests/reg_tests/test_pyGeo.py

@@ -18,34 +18,40 @@ def setUp(self):
 
     def train(self):
         with BaseRegTest(self.refFile, train=True) as handler:
-            self.regTest(handler)
+            self.regTest(handler, clash=False)
+            self.regTest(handler, clash=True)
             handler.writeRef()
 
-    def test(self):
+    def test_original(self):
         with BaseRegTest(self.refFile, train=False) as handler:
-            self.regTest(handler)
+            self.regTest(handler, clash=False)
 
-    def regTest(self, handler):
-        dirName = os.path.join(baseDir, "../../input_files")
+    def test_clash(self):
+        with BaseRegTest(self.refFile, train=False) as handler:
+            self.regTest(handler, clash=True)
+
+    def regTest(self, handler, clash=False):
+        inputDir = os.path.join(baseDir, "../../input_files")
 
-        # Airfoil file
-        airfoil_list = [dirName + "/rae2822.dat"] * 2
-        naf = len(airfoil_list)  # number of airfoils
+        if clash:
+            # This case, which uses two different airfoils with very similar point distributions, led to an error due
+            # to floating point precision issues in the clash detection algorithm. So now we test that it runs without
+            # error.
+            airfoil_list = [os.path.join(inputDir, "CRMTailRoot.dat"), os.path.join(inputDir, "CRMTailTip.dat")]
+        else:
+            airfoil_list = [os.path.join(inputDir, "rae2822.dat")] * 2
 
-        # Wing definition
-        # Airfoil leading edge positions
+        naf = len(airfoil_list)
+
+        # Wing definition (Common to both)
         x = [0.0, 7.5]
         y = [0.0, 0.0]
         z = [0.0, 14.0]
-        offset = np.zeros((naf, 2))  # x-y offset applied to airfoil position before scaling
-
-        # Airfoil rotations
+        offset = np.zeros((naf, 2))
         rot_x = [0.0, 0.0]
         rot_y = [0.0, 0.0]
         rot_z = [0.0, 0.0]
-
-        # Airfoil scaling
-        chord = [5.0, 1.5]  # chord lengths
+        chord = [5.0, 1.5]
 
         # Run pyGeo
         wing = pyGeo(
@@ -68,4 +74,9 @@ def regTest(self, handler):
         for isurf in range(wing.nSurf):
             wing.surfs[isurf].computeData()
         surf = wing.surfs[isurf].data
-        handler.root_add_val("sum of surface data", sum(surf.flatten()), tol=1e-10)
+        label = "clash" if clash else "original"
+        handler.root_add_val(f"{label} sum of surface data", sum(surf.flatten()), tol=1e-10)
+
+
+if __name__ == "__main__":
+    unittest.main()