Merge branch 'develop' into supersonic_profile

Common/include/geometry/CGeometry.hpp

@@ -240,8 +240,10 @@ class CGeometry {
   unsigned long* nVertex{nullptr};     /*!< \brief Number of vertex for each marker. */
   unsigned long* nElem_Bound{nullptr}; /*!< \brief Number of elements of the boundary. */
   string* Tag_to_Marker{nullptr};      /*!< \brief Names of boundary markers. */
-  vector<bool>
-      bound_is_straight; /*!< \brief Bool if boundary-marker is straight(2D)/plane(3D) for each local marker. */
+
+  /*!< \brief Corrected normals on nodes with shared symmetry markers. */
+  vector<std::unordered_map<unsigned long, std::array<su2double, MAXNDIM>>> symmetryNormals;
+
   vector<su2double> SurfaceAreaCfgFile; /*!< \brief Total Surface area for all markers. */
 
   /*--- Partitioning-specific variables ---*/
@@ -819,6 +821,12 @@ class CGeometry {
    */
   inline virtual void SetBoundControlVolume(const CConfig* config, unsigned short action) {}
 
+  /*!
+   * \brief Computes modified normals at intersecting symmetry planes.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void ComputeModifiedSymmetryNormals(const CConfig* config);
+
   /*!
    * \brief A virtual member.
    * \param[in] config_filename - Name of the file where the tecplot information is going to be stored.
@@ -936,9 +944,10 @@ class CGeometry {
   /*!
    * \brief A virtual member.
    * \param[in] fine_grid - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
    * \param[in] action - Allocate or not the new elements.
    */
-  inline virtual void SetBoundControlVolume(const CGeometry* fine_grid, unsigned short action) {}
+  inline virtual void SetBoundControlVolume(const CGeometry* fine_grid, const CConfig* config, unsigned short action) {}
 
   /*!
    * \brief A virtual member.
@@ -1005,15 +1014,6 @@ class CGeometry {
    */
   su2double GetSurfaceArea(const CConfig* config, unsigned short val_marker) const;
 
-  /*!
-   * \brief Check if a boundary is straight(2D) / plane(3D) for EULER_WALL and SYMMETRY_PLANE
-   *        only and store the information in bound_is_straight. For all other boundary types
-   *        this will return false and could therfore be wrong. Used ultimately for BC_Slip_Wall.
-   * \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);
-
   /*!
    * \brief Find and store all vertices on a sharp corner in the geometry.
    * \param[in] config - Definition of the particular problem.

Common/include/geometry/CMultiGridGeometry.hpp

@@ -104,9 +104,10 @@ class CMultiGridGeometry final : public CGeometry {
   /*!
    * \brief Set boundary vertex structure of the agglomerated control volume.
    * \param[in] fine_grid - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
    * \param[in] action - Allocate or not the new elements.
    */
-  void SetBoundControlVolume(const CGeometry* fine_grid, unsigned short action) override;
+  void SetBoundControlVolume(const CGeometry* fine_grid, const CConfig* config, unsigned short action) override;
 
   /*!
    * \brief Set a representative coordinates of the agglomerated control volume.

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

@@ -314,10 +314,8 @@ class CPoint {
    * \return Index of the vertex.
    */
   inline long GetVertex(unsigned long iPoint, unsigned long iMarker) const {
-    if (Boundary(iPoint))
-      return Vertex[iPoint][iMarker];
-    else
-      return -1;
+    if (Boundary(iPoint)) return Vertex[iPoint][iMarker];
+    return -1;
   }
 
   /*!
@@ -369,7 +367,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

@@ -2338,7 +2338,7 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config)
     /*--- Update the control volume structures ---*/
 
     geometry_container[iMesh]->SetControlVolume(geometry_container[iMesh - 1], UPDATE);
-    geometry_container[iMesh]->SetBoundControlVolume(geometry_container[iMesh - 1], UPDATE);
+    geometry_container[iMesh]->SetBoundControlVolume(geometry_container[iMesh - 1], config, UPDATE);
     geometry_container[iMesh]->SetCoord(geometry_container[iMesh - 1]);
   }
 
@@ -2453,123 +2453,71 @@ su2double CGeometry::GetSurfaceArea(const CConfig* config, unsigned short val_ma
   return 0.0;
 }
 
-void CGeometry::ComputeSurf_Straightness(CConfig* config, bool print_on_screen) {
-  bool RefUnitNormal_defined;
-  unsigned short iDim, iMarker, iMarker_Global, nMarker_Global = config->GetnMarker_CfgFile();
-  unsigned long iVertex;
-  constexpr passivedouble epsilon = 1.0e-6;
-  su2double Area;
-  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
-        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. ---*/
-  bound_is_straight.resize(nMarker);
-  fill(bound_is_straight.begin(), bound_is_straight.end(), true);
-
-  /*--- Loop over all local markers ---*/
-  for (iMarker = 0; iMarker < nMarker; iMarker++) {
-    Local_TagBound = config->GetMarker_All_TagBound(iMarker);
-
-    /*--- Marker has to be Symmetry or Euler. Additionally marker can't be a
-          moving surface and Grid Movement Elasticity is forbidden as well. All
-          other GridMovements are rigid. ---*/
-    if ((config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE ||
-         config->GetMarker_All_KindBC(iMarker) == EULER_WALL) &&
-        !config->GetMarker_Moving_Bool(Local_TagBound) && !config->GetMarker_Deform_Mesh_Bool(Local_TagBound)) {
-      /*--- Loop over all global markers, and find the local-global pair via
-            matching unique string tags. ---*/
-      for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
-        Global_TagBound = config->GetMarker_CfgFile_TagBound(iMarker_Global);
-        if (Local_TagBound == Global_TagBound) {
-          RefUnitNormal_defined = false;
-          iVertex = 0;
-
-          while (bound_is_straight[iMarker] && iVertex < nVertex[iMarker]) {
-            vertex[iMarker][iVertex]->GetNormal(Normal.data());
-            UnitNormal = Normal;
-
-            /*--- Compute unit normal. ---*/
-            Area = 0.0;
-            for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim] * Normal[iDim];
-            Area = sqrt(Area);
-
-            /*--- Negate for outward convention. ---*/
-            for (iDim = 0; iDim < nDim; iDim++) UnitNormal[iDim] /= -Area;
-
-            /*--- Check if unit normal is within tolerance of the Reference unit normal.
-                  Reference unit normal = first unit normal found. ---*/
-            if (RefUnitNormal_defined) {
-              for (iDim = 0; iDim < nDim; iDim++) {
-                if (abs(RefUnitNormal[iDim] - UnitNormal[iDim]) > epsilon) {
-                  bound_is_straight[iMarker] = false;
-                  break;
-                }
-              }
-            } else {
-              RefUnitNormal = UnitNormal;  // deep copy of values
-              RefUnitNormal_defined = true;
-            }
+void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) {
+  /* 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. */
+  symmetryNormals.resize(nMarker);
+  std::vector<unsigned short> symMarkers;
 
-            iVertex++;
-          }  // while iVertex
-        }    // if Local == Global
-      }      // for iMarker_Global
-    } else {
-      /*--- Enforce default value: false ---*/
-      bound_is_straight[iMarker] = false;
-    }  // if sym or euler ...
-  }    // for iMarker
-
-  /*--- Communicate results and print on screen. ---*/
-  if (print_on_screen) {
-    /*--- Additional vector which can later be MPI::Allreduce(d) to pring the results
-          on screen as nMarker (local) can vary across ranks. Default 'true' as it can
-          happen that a local rank does not contain an element of each surface marker.  ---*/
-    vector<bool> bound_is_straight_Global(nMarker_Global, true);
-    /*--- Match local with global tag bound and fill a Global Marker vector. ---*/
-    for (iMarker = 0; iMarker < nMarker; iMarker++) {
-      Local_TagBound = config->GetMarker_All_TagBound(iMarker);
-      for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
-        Global_TagBound = config->GetMarker_CfgFile_TagBound(iMarker_Global);
-
-        if (Local_TagBound == Global_TagBound) bound_is_straight_Global[iMarker_Global] = bound_is_straight[iMarker];
-
-      }  // for iMarker_Global
-    }    // for iMarker
-
-    vector<int> Buff_Send_isStraight(nMarker_Global), Buff_Recv_isStraight(nMarker_Global);
-
-    /*--- Cast to int as std::vector<boolean> can be a special construct. MPI handling using <int>
-          is more straight-forward. ---*/
-    for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++)
-      Buff_Send_isStraight[iMarker_Global] = static_cast<int>(bound_is_straight_Global[iMarker_Global]);
-
-    /*--- Product of type <int>(bool) is equivalnt to a 'logical and' ---*/
-    SU2_MPI::Allreduce(Buff_Send_isStraight.data(), Buff_Recv_isStraight.data(), nMarker_Global, MPI_INT, MPI_PROD,
-                       SU2_MPI::GetComm());
-
-    /*--- Print results on screen. ---*/
-    if (rank == MASTER_NODE) {
-      for (iMarker_Global = 0; iMarker_Global < nMarker_Global; iMarker_Global++) {
-        if (config->GetMarker_CfgFile_KindBC(config->GetMarker_CfgFile_TagBound(iMarker_Global)) == SYMMETRY_PLANE ||
-            config->GetMarker_CfgFile_KindBC(config->GetMarker_CfgFile_TagBound(iMarker_Global)) == EULER_WALL) {
-          cout << "Boundary marker " << config->GetMarker_CfgFile_TagBound(iMarker_Global) << " is";
-          if (!static_cast<bool>(Buff_Recv_isStraight[iMarker_Global])) cout << " NOT";
-          if (nDim == 2) cout << " a single straight." << endl;
-          if (nDim == 3) cout << " a single plane." << endl;
-        }  // if sym or euler
-      }    // for iMarker_Global
-    }      // if rank==MASTER
-  }        // if print_on_scren
+  for (auto iMarker = 0u; iMarker < nMarker; ++iMarker) {
+    if ((config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE) ||
+        (config->GetMarker_All_KindBC(iMarker) == EULER_WALL)) {
+      symMarkers.push_back(iMarker);
+    }
+  }
+
+  /*--- Loop over all markers and find nodes on symmetry planes that are shared with other symmetries. ---*/
+  /*--- The first symmetry does not need a corrected normal vector, hence start at 1. ---*/
+  for (size_t i = 1; i < symMarkers.size(); ++i) {
+    const auto iMarker = symMarkers[i];
+
+    for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
+      const auto iPoint = vertex[iMarker][iVertex]->GetNode();
+
+      /*--- Halo points do not need to be considered. ---*/
+      if (!nodes->GetDomain(iPoint)) continue;
+
+      /*--- Get the vertex normal on the current symmetry. ---*/
+      std::array<su2double, MAXNDIM> iNormal = {};
+      vertex[iMarker][iVertex]->GetNormal(iNormal.data());
+
+      /*--- Loop over previous symmetries and if this point shares them, make this normal orthogonal to them. ---*/
+      bool isShared = false;
+
+      for (size_t j = 0; j < i; ++j) {
+        const auto jMarker = symMarkers[j];
+        const auto jVertex = nodes->GetVertex(iPoint, jMarker);
+        if (jVertex < 0) continue;
+
+        isShared = true;
+
+        std::array<su2double, MAXNDIM> jNormal = {};
+        const auto it = symmetryNormals[jMarker].find(jVertex);
+
+        if (it != symmetryNormals[jMarker].end()) {
+          jNormal = it->second;
+        } else {
+          vertex[jMarker][jVertex]->GetNormal(jNormal.data());
+          const su2double area = GeometryToolbox::Norm(nDim, jNormal.data());
+          for (auto iDim = 0u; iDim < nDim; iDim++) jNormal[iDim] /= area;
+        }
+
+        const auto proj = GeometryToolbox::DotProduct(nDim, jNormal.data(), iNormal.data());
+        for (auto iDim = 0u; iDim < nDim; iDim++) iNormal[iDim] -= proj * jNormal[iDim];
+      }
+
+      if (!isShared) continue;
+
+      /*--- 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(nDim, iNormal.data());
+      if (area > EPS) {
+        for (auto iDim = 0u; iDim < nDim; iDim++) iNormal[iDim] /= area;
+        symmetryNormals[iMarker][iVertex] = iNormal;
+      }
+    }
+  }
 }
 
 void CGeometry::ComputeSurf_Curvature(CConfig* config) {