fix symmetry shared node (#2719)

* fix symmetry shared node

Author: bigfooted

Common/src/geometry/CGeometry.cpp

@@ -2468,6 +2468,9 @@ su2double CGeometry::GetSurfaceArea(const CConfig* config, unsigned short val_ma
 }
 
 void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) {
+  const su2double MIN_AREA = 1e-12;            // minimum area to consider a normal valid
+  const su2double PARALLEL_TOLERANCE = 0.001;  // min. angle between symm. planes to consider them non-parallel.
+
   /* Check how many symmetry planes there are and use the first (lowest ID) as the basis to orthogonalize against.
    * All nodes that are shared by multiple symmetries have to get a corrected normal. */
 
@@ -2558,6 +2561,8 @@ void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) {
       /*--- Loop over previous symmetries and if this point shares them, make this normal orthogonal to them.
        * It's ok if we normalize merged normals against themselves, we get 0 area and this becomes a no-op. ---*/
 
+      std::vector<size_t> parallelMarkers;  // Track markers with nearly-parallel normals
+
       for (size_t j = 0; j < i; ++j) {
         const auto jMarker = symMarkers[j];
         const auto jVertex = nodes->GetVertex(iPoint, jMarker);
@@ -2567,14 +2572,48 @@ void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) {
         GetNormal(jMarker, jVertex, jNormal);
 
         const su2double proj = GeometryToolbox::DotProduct(int(MAXNDIM), jNormal.data(), iNormal.data());
+        const su2double angleDiff = std::abs(1.0 - std::abs(proj));
+
+        // Check if normals are nearly parallel (within ~2.5 degrees)
+        // cos(2.5°) ≈ 0.999, so (1 - cos(2.5°)) ≈ 0.001
+        if (angleDiff < PARALLEL_TOLERANCE) {
+          // These normals are nearly parallel - average them instead of orthogonalizing
+          parallelMarkers.push_back(j);
+          for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] += jNormal[iDim];
+          continue;
+        }
+
         for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] -= proj * jNormal[iDim];
       }
 
+      /*--- If we found parallel markers, average and store the result for all involved markers ---*/
+      if (!parallelMarkers.empty()) {
+        // Normalize the averaged normal
+        const su2double avgArea = GeometryToolbox::Norm(int(MAXNDIM), iNormal.data());
+        if (avgArea > MIN_AREA) {
+          for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] /= avgArea;
+
+          // Store the averaged normal for the current marker
+          symmetryNormals[iMarker][iVertex] = iNormal;
+
+          // Also update all parallel markers with the same averaged normal
+          for (const auto j : parallelMarkers) {
+            const auto jMarker = symMarkers[j];
+            const auto jVertex = nodes->GetVertex(iPoint, jMarker);
+            if (jVertex >= 0) {
+              symmetryNormals[jMarker][jVertex] = iNormal;
+            }
+          }
+        }
+        continue;  // Skip the normal orthogonalization path below
+      }
+
       /*--- Normalize. If the norm is close to zero it means the normal is a linear combination of previous
        * normals, in this case we don't need to store the corrected normal, using the original in the gradient
        * correction will have no effect since previous corrections will remove components in this direction). ---*/
       const su2double area = GeometryToolbox::Norm(int(MAXNDIM), iNormal.data());
-      if (area > 1e-12) {
+
+      if (area > MIN_AREA) {
         for (auto iDim = 0ul; iDim < MAXNDIM; ++iDim) iNormal[iDim] /= area;
         symmetryNormals[iMarker][iVertex] = iNormal;
       }
@@ -2591,16 +2630,14 @@ void CGeometry::ComputeSurfStraightness(const CConfig* config, bool print_on_scr
   string Local_TagBound, Global_TagBound;
 
   vector<su2double> Normal(nDim), UnitNormal(nDim), RefUnitNormal(nDim);
-
   /*--- Assume now that this boundary marker is straight. As soon as one
-        AreaElement is found that is not aligend with a Reference then it is
-        certain that the boundary marker is not straight and one can stop
-        searching. Another possibility is that this process doesn't own
+        AreaElement is found that is not aligned with a Reference then it
+        is certain that the boundary marker is not straight and one can
+        stop searching. Another possibility is that this process doesn't own
         any nodes of that boundary, in that case we also have to assume the
-        boundary is straight.
-        Any boundary type other than SYMMETRY_PLANE or EULER_WALL gets
-        the value false (or see cases specified in the conditional below)
-        which could be wrong. ---*/
+        boundary is straight. Any boundary type other than SYMMETRY_PLANE or
+        EULER_WALL gets the value false (or see cases specified in the
+        conditional below) which could be wrong. ---*/
   boundIsStraight.resize(nMarker);
   fill(boundIsStraight.begin(), boundIsStraight.end(), true);
 
@@ -3902,11 +3939,13 @@ void CGeometry::ColorMGLevels(unsigned short nMGLevels, const CGeometry* const*
     for (auto step = 0u; step < iMesh; ++step) {
       auto coarseMesh = geometry[iMesh - 1 - step];
       if (step)
-        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint)
+        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint) {
           CoarseGridColor_(iPoint, step) = CoarseGridColor_(coarseMesh->nodes->GetParent_CV(iPoint), step - 1);
+        }
       else
-        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint)
+        for (auto iPoint = 0ul; iPoint < coarseMesh->GetnPoint(); ++iPoint) {
           CoarseGridColor_(iPoint, step) = color[coarseMesh->nodes->GetParent_CV(iPoint)];
+        }
     }
   }
 }

TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref

@@ -1,9 +1,9 @@
 INDEX	GRAD
-0	-4.570880175564371e-04
-1	-2.401474104649833e-04
-2	-9.134741714239637e-05
-3	-1.628103592186675e-05
-4	-1.741046982337274e-05
-5	-9.462481735830197e-05
-6	-2.452469745130044e-04
-7	-4.635498256044697e-04
+0	-4.570880019852262e-04
+1	-2.401473990776599e-04
+2	-9.134741024986813e-05
+3	-1.628103377709482e-05
+4	-1.741047267006423e-05
+5	-9.462482550045010e-05
+6	-2.452469883016078e-04
+7	-4.635498453796528e-04

TestCases/parallel_regression.py

@@ -1526,7 +1526,7 @@ def main():
     pywrapper_deformingBump.cfg_dir = "py_wrapper/deforming_bump_in_channel"
     pywrapper_deformingBump.cfg_file = "config.cfg"
     pywrapper_deformingBump.test_iter = 1
-    pywrapper_deformingBump.test_vals = [0.500000, 0.000000, -2.811520, -1.603562, -2.074259, 2.424289, 7.616891, -0.205655]
+    pywrapper_deformingBump.test_vals = [0.500000, 0.000000, -2.556309, -1.270839, -2.350590, 2.606851, 8.002480, -0.300272]
     pywrapper_deformingBump.command = TestCase.Command("mpirun -np 2", "python", "run.py")
     pywrapper_deformingBump.unsteady = True
     test_list.append(pywrapper_deformingBump)

TestCases/parallel_regression_AD.py

@@ -276,8 +276,8 @@ def main():
     discadj_fsi2.cfg_dir   = "disc_adj_fsi/Airfoil_2d"
     discadj_fsi2.cfg_file  = "config.cfg"
     discadj_fsi2.test_iter = 8
-    discadj_fsi2.test_vals         = [-4.772768, 0.916597, -3.863369, 0.295450, 3.839800]
-    discadj_fsi2.test_vals_aarch64 = [-4.773008, 0.916312, -3.863369, 0.295450, 3.841200]
+    discadj_fsi2.test_vals         = [-4.772641, 0.917601, -3.863369, 0.295450, 3.839800]
+    discadj_fsi2.test_vals_aarch64 = [-4.772641, 0.917601, -3.863369, 0.295450, 3.839800]
     discadj_fsi2.tol       = 0.00001
     test_list.append(discadj_fsi2)
 

TestCases/tutorials.py

@@ -292,7 +292,7 @@ def main():
     tutorial_nicfd_nozzle_pinn.cfg_dir   = "../Tutorials/compressible_flow/NICFD_nozzle/PhysicsInformed"
     tutorial_nicfd_nozzle_pinn.cfg_file  = "config_NICFD_PINN.cfg"
     tutorial_nicfd_nozzle_pinn.test_iter = 20
-    tutorial_nicfd_nozzle_pinn.test_vals = [-3.181747, -1.638856, -1.277037, 2.445964, -11.759632]
+    tutorial_nicfd_nozzle_pinn.test_vals = [-2.728179, -0.849337, -1.224542, 2.898995, -11.420290]
     tutorial_nicfd_nozzle_pinn.no_restart = True
     test_list.append(tutorial_nicfd_nozzle_pinn)