[WIP] Fix cornernode normal

Common/include/geometry/CMultiGridGeometry.hpp

@@ -123,7 +123,7 @@ class CMultiGridGeometry final : public CGeometry {
   void SetRestricted_GridVelocity(const CGeometry* fine_grid) override;
 
   /*!
-   * \brief Find and store the closest neighbor to a vertex.
+   * \brief Find and store the closest, most normal, neighbor to a vertex.
    * \param[in] config - Definition of the particular problem.
    */
   void FindNormal_Neighbor(const CConfig* config) override;

Common/include/geometry/CPhysicalGeometry.hpp

@@ -560,7 +560,7 @@ class CPhysicalGeometry final : public CGeometry {
   void ComputeMeshQualityStatistics(const CConfig* config) override;
 
   /*!
-   * \brief Find and store the closest neighbor to a vertex.
+   * \brief Find and store the closest, most normal, neighbor to a vertex.
    * \param[in] config - Definition of the particular problem.
    */
   void FindNormal_Neighbor(const CConfig* config) override;

Common/include/geometry/dual_grid/CVertex.hpp

@@ -45,7 +45,8 @@ class CVertex : public CDualGrid {
   long PeriodicPoint[5] = {-1};  /*!< \brief Store the periodic point of a boundary (iProcessor, iPoint) */
   bool ActDisk_Perimeter = false;      /*!< \brief Identify nodes at the perimeter of the actuator disk */
   short Rotation_Type;                 /*!< \brief Type of rotation associated with the vertex (MPI and periodic) */
-  unsigned long Normal_Neighbor;       /*!< \brief Index of the closest neighbor. */
+  unsigned long Normal_Neighbor;       /*!< \brief Index of the closest, most normal, neighbor. */
+  unsigned long Nearest_Neighbor;      /*!< \brief Index of the closest interior neighbor. */
   su2double Basis_Function[3] = {0.0}; /*!< \brief Basis function values for interpolation across zones. */
 
  public:

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -2419,6 +2419,14 @@ class CFVMFlowSolverBase : public CSolver {
     return YPlus[val_marker][val_vertex];
   }
 
+  /*!
+   * \brief Set the y plus.
+   * \param[in] val_yplus - new value of yplus
+   */
+  inline void SetYPlus(su2double val_yplus, unsigned short val_marker, unsigned long val_vertex) final {
+    YPlus[val_marker][val_vertex] = val_yplus;
+  }
+
   /*!
    * \brief Get the u_tau .
    * \param[in] val_marker - Surface marker where the coefficient is computed.

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2498,10 +2498,15 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       FrictionVel = sqrt(fabs(WallShearStress[iMarker][iVertex]) / Density);
 
       if (!wallfunctions && (MGLevel == MESH_0 || geometry->nodes->GetDomain(iPoint))) {
-        // for CMultiGridGeometry, the normal neighbor of halo nodes in not set
+        // for CMultiGridGeometry, the normal neighbor of halo nodes is not set
         iPointNormal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
         Coord_Normal = geometry->nodes->GetCoord(iPointNormal);
         WallDistMod = GeometryToolbox::Distance(nDim, Coord, Coord_Normal);
+        //nijso: we can modifiy y+ as well now?
+
+
+
+
         YPlus[iMarker][iVertex] = WallDistMod * FrictionVel / (Viscosity / Density);
       }
 

SU2_CFD/include/solvers/CSolver.hpp

@@ -2989,6 +2989,14 @@ class CSolver {
    */
   inline virtual su2double GetYPlus(unsigned short val_marker, unsigned long val_vertex) const { return 0; }
 
+  /*!
+   * \brief A virtual member.
+   * \param[in] val_marker - Surface marker where the coefficient is computed.
+   * \param[in] val_vertex - Vertex of the marker <i>val_marker</i> where the coefficient is evaluated.
+   * \return Value of the y plus.
+   */
+  inline virtual void SetYPlus(su2double val_yplus, unsigned short val_marker, unsigned long val_vertex) {};
+
   /*!
    * \brief A virtual member.
    * \param[in] val_marker - Surface marker where the coefficient is computed.

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -302,5 +302,6 @@ class CTurbSSTSolver final : public CTurbSolver {
    */
   inline su2double GetOmega_Inf(void) const override { return Solution_Inf[1]; }
 
+  su2double GetNearest_Neighbor(CGeometry *geometry, unsigned long iPoint,unsigned short iMarker, unsigned long iVertex);
 
 };

SU2_CFD/src/drivers/CDriver.cpp

@@ -833,7 +833,7 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) {
     geometry[iMGlevel]->SetBoundControlVolume(geometry[iMGlevel-1], config, ALLOCATE);
     geometry[iMGlevel]->SetCoord(geometry[iMGlevel-1]);
 
-    /*--- Find closest neighbor to a surface point ---*/
+    /*--- Find closest, most normal, neighbor to a surface point ---*/
 
     geometry[iMGlevel]->FindNormal_Neighbor(config);
 
@@ -2517,7 +2517,7 @@ void CDriver::InitializeInterface(CConfig **config, CSolver***** solver, CGeomet
             else
               interface_type = NO_TRANSFER;
           }
-          
+
           if (interface_type != NO_TRANSFER) {
             auto nVar = 4;
             interface[donor][target] = new CConjugateHeatInterface(nVar, 0);

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -259,7 +259,9 @@
 
           /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/
 
+
           if (geometry->nodes->GetDomain(iPoint)) {
+
             const auto jPoint = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
 
             su2double shearStress = 0.0;
@@ -269,7 +271,12 @@
             shearStress = sqrt(shearStress);
 
             const su2double FrictionVelocity = sqrt(shearStress/flowNodes->GetDensity(iPoint));
-            const su2double wall_dist = geometry->nodes->GetWall_Distance(jPoint);
+            const su2double wall_dist_old = geometry->nodes->GetWall_Distance(jPoint);
+
+            const su2double wall_dist = GetNearest_Neighbor(geometry,iPoint,iMarker,iVertex);
+
+            cout << "************************* wall distance old = " << wall_dist_old << " , wall distance new = " << wall_dist << endl;
+
             const su2double Derivative = flowNodes->GetLaminarViscosity(jPoint) * pow(nodes->GetSolution(jPoint, 0), 0.673) / wall_dist;
             const su2double turbulence_index = 6.1 * Derivative / pow(FrictionVelocity, 2.346);
 
@@ -459,11 +466,12 @@
       } else { // smooth wall
 
         /*--- distance to closest neighbor ---*/
-        const auto jPoint = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
+        su2double wall_dist = GetNearest_Neighbor(geometry,iPoint,val_marker, iVertex);
 
-        su2double distance2 = GeometryToolbox::SquaredDistance(nDim,
-                                                             geometry->nodes->GetCoord(iPoint),
-                                                             geometry->nodes->GetCoord(jPoint));
+        const auto jPoint = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
+        //su2double distance2 = GeometryToolbox::SquaredDistance(nDim,
+        //                                                     geometry->nodes->GetCoord(iPoint),
+        //                                                     geometry->nodes->GetCoord(jPoint));
         /*--- Set wall values ---*/
 
         su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(jPoint);
@@ -472,7 +480,8 @@
         su2double beta_1 = constants[4];
         su2double solution[MAXNVAR];
         solution[0] = 0.0;
-        solution[1] = 60.0*laminar_viscosity/(density*beta_1*distance2);
+        solution[1] = 60.0*laminar_viscosity/(density*beta_1*wall_dist*wall_dist);
+        //solution[1] = 60.0*laminar_viscosity/(density*beta_1*distance2);
 
         /*--- Set the solution values and zero the residual ---*/
         nodes->SetSolution_Old(iPoint,solution);
@@ -1012,8 +1021,7 @@
 
   su2double Normal[MAXNDIM] = {0.0};
   geometry->vertex[iMarker][iVertex]->GetNormal(Normal);
-  return GeometryToolbox::Norm(nDim, Normal);
-}
+  return GeometryToolbox::Norm(nDim, Normal);}
 
 void CTurbSSTSolver::SetUniformInlet(const CConfig* config, unsigned short iMarker) {
   if (config->GetMarker_All_KindBC(iMarker) == INLET_FLOW) {
@@ -1024,3 +1032,59 @@
   }
 
 }
+
+/*--- We determine the interior node that is closest to the wall node. If there is no
+      interior node, we have to take the closest wall node. ---*/
+su2double  CTurbSSTSolver::GetNearest_Neighbor(CGeometry *geometry, unsigned long iPoint, unsigned short iMarker, unsigned long iVertex) {
+  su2double dist_min;
+  su2double distance = 0.0;
+  unsigned long Point_Normal;
+
+  Point_Normal = 0;
+  dist_min = 0.0;
+  for (size_t iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) {
+    size_t jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+    if (geometry->nodes->GetViscousBoundary(jPoint) == false) {
+      su2double distance[MAXNDIM] = {0.0};
+      // compute distance between ipoint and jpoint (vector)
+      GeometryToolbox::Distance(nDim,
+                                geometry->nodes->GetCoord(iPoint),
+                                geometry->nodes->GetCoord(jPoint),
+                                distance);
+      // get the edge
+      size_t iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
+      const su2double* Normal = geometry->edges->GetNormal(iEdge);
+      const su2double Area = GeometryToolbox::Norm(nDim, Normal);
+      su2double projdistance = abs(GeometryToolbox::DotProduct(nDim,distance,Normal) / Area);
+
+      // Take the interior node that is closest to the wall node.
+      //  cout << "projected distance = " << projdistance << " , area="<< Area << endl;
+
+      Point_Normal = jPoint;
+      dist_min = projdistance;
+    }
+  }
+
+
+  if (Point_Normal==0) {
+    //cout << "no point found at i = " << iPoint << endl;
+    su2double Area = 0.0;
+    su2double TwoVol = 2.0* (geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint));
+
+    for (size_t iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+      size_t iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
+      size_t jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
+
+      if (geometry->nodes->GetViscousBoundary(jPoint) == true) {
+        const su2double* Normal = geometry->edges->GetNormal(iEdge);
+        Area = GeometryToolbox::Norm(nDim, Normal);
+        dist_min += TwoVol / Area;
+        //cout << "    distmin = " << dist_min << " , vol = "<<TwoVol << " , Area=" << Area<< endl;
+      }
+
+    }
+  //cout << "distmin = " << dist_min << " , vol = "<<TwoVol << " , Area=" << Area<< endl;
+  }
+
+  return (dist_min);
+}