Correction of symmetry plane implementation

.github/workflows/regression.yml

@@ -213,7 +213,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t fix_symmetry -c fix_symmetry -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
         with:

Common/include/geometry/CGeometry.hpp

@@ -66,6 +66,18 @@ extern "C" {
 
 using namespace std;
 
+struct symNormal {
+  long unsigned int index;
+  su2double normal[3];
+  symNormal(long unsigned int index) : index(index){};
+};
+
+struct findSymNormalIndex : std::unary_function<symNormal, bool> {
+  unsigned long int index;
+  findSymNormalIndex(unsigned long int index) : index(index) {}
+  bool operator()(symNormal const& m) const { return m.index == index; }
+};
+
 /*!
  * \class CGeometry
  * \brief Parent class for defining the geometry of the problem (complete geometry,
@@ -229,6 +241,7 @@ class CGeometry {
 
  public:
   /*--- Main geometric elements of the grid. ---*/
+  vector<vector<symNormal>> symmetryNormals; /*!< \brief Corrected normals on nodes with shared symmetry markers */
 
   CPrimalGrid** elem{nullptr};   /*!< \brief Element vector (primal grid information). */
   CPrimalGrid*** bound{nullptr}; /*!< \brief Boundary vector (primal grid information). */
@@ -1012,7 +1025,7 @@ class CGeometry {
    * \param[in] config - Definition of the particular problem.
    * \param[in] print_on_screen - Boolean whether to print result on screen.
    */
-  void ComputeSurf_Straightness(CConfig* config, bool print_on_screen);
+  void ComputeSurf_Straightness(const CConfig* config, bool print_on_screen);
 
   /*!
    * \brief Find and store all vertices on a sharp corner in the geometry.

Common/include/geometry/CPhysicalGeometry.hpp

@@ -208,10 +208,13 @@ class CPhysicalGeometry final : public CGeometry {
    * \brief Routine to launch non-blocking sends and recvs amongst all processors.
    * \param[in] bufSend - Buffer of data to be sent.
    * \param[in] nElemSend - Array containing the number of elements to send to other processors in cumulative storage
-   * format. \param[in] sendReq - Array of MPI send requests. \param[in] bufRecv - Buffer of data to be received.
+   * format.
+   * \param[in] sendReq - Array of MPI send requests.
+   * \param[in] bufRecv - Buffer of data to be received.
    * \param[in] nElemSend - Array containing the number of elements to receive from other processors in cumulative
-   * storage format. \param[in] sendReq - Array of MPI recv requests. \param[in] countPerElem - Pieces of data per
-   * element communicated.
+   * storage format.
+   * \param[in] sendReq - Array of MPI recv requests.
+   * \param[in] countPerElem - Pieces of data per element communicated.
    */
   void InitiateCommsAll(void* bufSend, const int* nElemSend, SU2_MPI::Request* sendReq, void* bufRecv,
                         const int* nElemRecv, SU2_MPI::Request* recvReq, unsigned short countPerElem,

Common/include/geometry/dual_grid/CPoint.hpp

@@ -369,7 +369,7 @@ class CPoint {
   inline bool GetPhysicalBoundary(unsigned long iPoint) const { return PhysicalBoundary(iPoint); }
 
   /*!
-   * \brief Set if a point belong to the boundary.
+   * \brief Set if a point belong to the solid wall boundary.
    * \param[in] iPoint - Index of the point.
    * \param[in] boundary - <code>TRUE</code> if the point belong to the physical boundary; otherwise <code>FALSE</code>.
    */

Common/include/geometry/meshreader/CCGNSMeshReaderFVM.hpp

@@ -121,10 +121,13 @@ class CCGNSMeshReaderFVM : public CMeshReaderFVM {
    * \brief Routine to launch non-blocking sends and recvs amongst all processors.
    * \param[in] bufSend - Buffer of data to be sent.
    * \param[in] nElemSend - Array containing the number of elements to send to other processors in cumulative storage
-   * format. \param[in] sendReq - Array of MPI send requests. \param[in] bufRecv - Buffer of data to be received.
+   * format.
+   * \param[in] sendReq - Array of MPI send requests.
+   * \param[in] bufRecv - Buffer of data to be received.
    * \param[in] nElemSend - Array containing the number of elements to receive from other processors in cumulative
-   * storage format. \param[in] sendReq - Array of MPI recv requests. \param[in] countPerElem - Pieces of data per
-   * element communicated.
+   * storage format.
+   * \param[in] sendReq - Array of MPI recv requests.
+   * \param[in] countPerElem - Pieces of data per element communicated.
    */
   void InitiateCommsAll(void* bufSend, const int* nElemSend, SU2_MPI::Request* sendReq, void* bufRecv,
                         const int* nElemRecv, SU2_MPI::Request* recvReq, unsigned short countPerElem,

Common/src/geometry/CGeometry.cpp

@@ -2453,7 +2453,7 @@ su2double CGeometry::GetSurfaceArea(const CConfig* config, unsigned short val_ma
   return 0.0;
 }
 
-void CGeometry::ComputeSurf_Straightness(CConfig* config, bool print_on_screen) {
+void CGeometry::ComputeSurf_Straightness(const CConfig* config, bool print_on_screen) {
   bool RefUnitNormal_defined;
   unsigned short iDim, iMarker, iMarker_Global, nMarker_Global = config->GetnMarker_CfgFile();
   unsigned long iVertex;

Common/src/geometry/CMultiGridGeometry.cpp

@@ -108,17 +108,27 @@
          This can be improved. If there is only a marker, it is a good
          candidate for agglomeration ---*/
 
-        if (counter == 1) agglomerate_seed = true;
+        if (counter == 1) {
+          agglomerate_seed = true;
 
-        /*--- If there are two markers, we will aglomerate if any of the
+          /*--- Euler walls can be curved and agglomerating them leads to difficulties ---*/
+          if (config->GetMarker_All_KindBC(marker_seed) == EULER_WALL) agglomerate_seed = false;
+        }
+        /*--- If there are two markers, we will agglomerate if any of the
          markers is SEND_RECEIVE ---*/
 
         if (counter == 2) {
           agglomerate_seed = (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) ||
                              (config->GetMarker_All_KindBC(copy_marker[1]) == SEND_RECEIVE);
+
+          /* --- Euler walls can also not be agglomerated when the point has 2 markers ---*/
+          if ((config->GetMarker_All_KindBC(copy_marker[0]) == EULER_WALL) ||
+              (config->GetMarker_All_KindBC(copy_marker[1]) == EULER_WALL)) {
+            agglomerate_seed = false;
+          }
         }
 
-        /*--- If there are more than 2 markers, the aglomeration will be discarted ---*/
+        /*--- If there are more than 2 markers, the aglomeration will be discarded ---*/
 
         if (counter > 2) agglomerate_seed = false;
 
@@ -220,7 +230,7 @@
     const auto iPoint = MGQueue_InnerCV.NextCV();
     iteration++;
 
-    /*--- If the element has not being previously agglomerated, belongs to the physical domain,
+    /*--- If the element has not been previously agglomerated, belongs to the physical domain,
      and satisfies several geometrical criteria then the seed CV is accepted for agglomeration. ---*/
 
     if ((!fine_grid->nodes->GetAgglomerate(iPoint)) && (fine_grid->nodes->GetDomain(iPoint)) &&
@@ -438,7 +448,6 @@
 
         /*--- Be careful, it is possible that a node changes the agglomeration configuration,
          the priority is always when receiving the information. ---*/
-
         fine_grid->nodes->SetParent_CV(iPoint_Fine, iPoint_Coarse);
         nodes->SetChildren_CV(iPoint_Coarse, nChildren_MPI[iPoint_Coarse], iPoint_Fine);
         nChildren_MPI[iPoint_Coarse]++;
@@ -501,7 +510,7 @@
                                                const CConfig* config) const {
   bool agglomerate_CV = false;
 
-  /*--- Basic condition, the point has not being previously agglomerated, it belongs to the domain,
+  /*--- Basic condition, the point has not been previously agglomerated, it belongs to the domain,
    and has passed some basic geometrical checks. ---*/
 
   if ((!fine_grid->nodes->GetAgglomerate(CVPoint)) && (fine_grid->nodes->GetDomain(CVPoint)) &&
@@ -511,6 +520,7 @@
     if (fine_grid->nodes->GetBoundary(CVPoint)) {
       /*--- Identify the markers of the vertex that we want to agglomerate ---*/
 
+      // count number of markers on the agglomeration candidate
       int counter = 0;
       unsigned short copy_marker[3] = {};
       for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker() && counter < 3; jMarker++) {
@@ -526,13 +536,22 @@
       /*--- Only one marker in the vertex that is going to be aglomerated ---*/
 
       if (counter == 1) {
-        /*--- We agglomerate if there is only a marker and is the same marker as the seed marker ---*/
-
+        /*--- We agglomerate if there is only one marker and it is the same marker as the seed marker ---*/
+        // note that this should be the same marker id, not just the same marker type
         if (copy_marker[0] == marker_seed) agglomerate_CV = true;
 
         /*--- If there is only one marker, but the marker is the SEND_RECEIVE ---*/
 
-        if (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) agglomerate_CV = true;
+        if (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) {
+          agglomerate_CV = true;
+        }
+
+        if ((config->GetMarker_All_KindBC(marker_seed) == SYMMETRY_PLANE) ||
+            (config->GetMarker_All_KindBC(marker_seed) == EULER_WALL)) {
+          if (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) {
+            agglomerate_CV = false;
+          }
+        }
       }
 
       /*--- If there are two markers in the vertex that is going to be aglomerated ---*/
@@ -541,20 +560,25 @@
         /*--- First we verify that the seed is a physical boundary ---*/
 
         if (config->GetMarker_All_KindBC(marker_seed) != SEND_RECEIVE) {
-          /*--- Then we check that one of the marker is equal to the seed marker, and the other is send/receive ---*/
+          /*--- Then we check that one of the markers is equal to the seed marker, and the other is send/receive ---*/
 
           if (((copy_marker[0] == marker_seed) && (config->GetMarker_All_KindBC(copy_marker[1]) == SEND_RECEIVE)) ||
-              ((config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) && (copy_marker[1] == marker_seed)))
+              ((config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) && (copy_marker[1] == marker_seed))) {
             agglomerate_CV = true;
+          }
         }
       }
-
     }
-
-    /*--- If the element belongs to the domain, it is allways aglomerated. ---*/
-
+    /*--- If the element belongs to the domain, it is always agglomerated. ---*/
     else {
       agglomerate_CV = true;
+
+      // actually, for symmetry (and possibly other cells) we only agglomerate cells that are on the marker
+      // at this point, the seed was on the boundary and the CV was not. so we check if the seed is a symmetry
+      if ((config->GetMarker_All_KindBC(marker_seed) == SYMMETRY_PLANE) ||
+          (config->GetMarker_All_KindBC(marker_seed) == EULER_WALL)) {
+        agglomerate_CV = false;
+      }
     }
   }
 
@@ -936,6 +960,89 @@
         if (Area == 0.0)
           for (iDim = 0; iDim < nDim; iDim++) NormalFace[iDim] = EPS * EPS;
       }
+
+    /* 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. */
+    static constexpr size_t MAXNSYMS = 100;
+    bool Syms[MAXNSYMS] = {false};
+    symmetryNormals.resize(nMarker);
+
+    unsigned short nSym = 0;
+    for (iMarker = 0; iMarker < nMarker; ++iMarker) {
+      /*--- create the list with all corrected normals for all markers ---*/
+      // symmetryNormals.push_back({});
+
+      /*--- Note that Syms is a sorted list, so val(Syms[i]) > val[Syms[i-1]] ---*/
+      if (fine_grid->symmetryNormals[iMarker].size() > 0) {
+        Syms[iMarker] = true;
+        nSym++;
+      }
+    }
+
+    // /*--- Loop over all markers and find nodes on symmetry planes that are shared with other symmetries. ---*/
+    for (unsigned short val_marker = 0; val_marker < nMarker; val_marker++) {
+      if (Syms[val_marker] == false) continue;
+      for (iVertex = 0; iVertex < nVertex[val_marker]; iVertex++) {
+        iCoarsePoint = vertex[val_marker][iVertex]->GetNode();
+
+        /*--- Halo points do not need to be considered. ---*/
+        if (!nodes->GetDomain(iCoarsePoint)) continue;
+
+        /*--- Get the normal of the current symmetry ---*/
+        su2double Normal[MAXNDIM] = {0.0}, UnitNormal[MAXNDIM] = {0.0};
+        vertex[val_marker][iVertex]->GetNormal(Normal);
+        const su2double Area = GeometryToolbox::Norm(nDim, Normal);
+        for (unsigned short iDim = 0; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim] / Area;
+
+        /*--- At this point we find out if the node is shared with another symmetry. ---*/
+        /*--- We then apply a Gram-Schmidt process to compute the orthonormal basis from all
+              normal vectors of the shared symmetries at the node. ---*/
+        /*--- Step 1: do we have other symmetries? ---*/
+        if (nSym > 0) {
+          /*--- Normal of the primary symmetry plane ---*/
+          su2double NormalPrim[MAXNDIM] = {0.0}, UnitNormalPrim[MAXNDIM] = {0.0};
+          /*---  Step 2: are we on a shared node? ---*/
+          for (auto iMarker = 0; iMarker < nMarker; iMarker++) {
+            if (Syms[iMarker] == false) continue;  // not on a symmetry
+            /*--- We do not want the current symmetry ---*/
+            if (val_marker != iMarker) {
+              /*--- Loop over all points on the other symmetry and check if current iPoint is on the symmetry ---*/
+              for (auto jVertex = 0ul; jVertex < nVertex[iMarker]; jVertex++) {
+                const auto jPoint = vertex[iMarker][jVertex]->GetNode();
+                if (!nodes->GetDomain(jPoint)) continue;
+                /*--- We are on a shared node.  ---*/
+                if (iCoarsePoint == jPoint) {
+                  /*--- Does the other symmetry have a lower ID? Then that is the primary symmetry ---*/
+                  if (iMarker < val_marker) {
+                    /*--- So we have to get the normal of that other marker ---*/
+                    vertex[iMarker][jVertex]->GetNormal(NormalPrim);
+                    su2double AreaPrim = GeometryToolbox::Norm(nDim, NormalPrim);
+                    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+                      UnitNormalPrim[iDim] = NormalPrim[iDim] / AreaPrim;
+                    }
+                    /*--- Correct the current normal as n2_new = n2 - (n2.n1)n1 ---*/
+                    su2double ProjNorm = 0.0;
+                    for (auto iDim = 0u; iDim < nDim; iDim++) ProjNorm += UnitNormal[iDim] * UnitNormalPrim[iDim];
+                    symNormal sn = {iCoarsePoint};
+
+                    /*--- We check if the normal of the 2 planes coincide.
+                     * We only update the normal if the normals of the symmetry planes are different. ---*/
+                    if (fabs(1.0 - ProjNorm) > EPS) {
+                      for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] -= ProjNorm * UnitNormalPrim[iDim];
+                      /*--- Make normalized vector ---*/
+                      su2double newarea = GeometryToolbox::Norm(nDim, UnitNormal);
+                      for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = UnitNormal[iDim] / newarea;
+                    }
+                    for (auto iDim = 0u; iDim < nDim; iDim++) sn.normal[iDim] = UnitNormal[iDim];
+                    symmetryNormals[val_marker].push_back(sn);
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+    }
   }
   END_SU2_OMP_SAFE_GLOBAL_ACCESS
 }