use the variable index classes in more places

Author: pcarruscag

Common/src/geometry/CPhysicalGeometry.cpp

@@ -4757,36 +4757,51 @@ void CPhysicalGeometry::SetRCM_Ordering(CConfig *config) {
   vector<unsigned long> AuxQueue, Result;
   Result.reserve(nPoint);
 
-  /*--- Select the node with the lowest degree in the grid. ---*/
+  while (Result.size() < nPointDomain) {
 
-  unsigned long AddPoint = 0;
-  auto MinDegree = nodes->GetnPoint(AddPoint);
-  for (auto iPoint = 1ul; iPoint < nPointDomain; iPoint++) {
-    auto Degree = nodes->GetnPoint(iPoint);
-    if (Degree < MinDegree) { MinDegree = Degree; AddPoint = iPoint; }
-  }
+    /*--- Select the node with the lowest degree in the grid. ---*/
 
-  /*--- Add the node in the first free position. ---*/
+    auto AddPoint = nPoint;
+    auto MinDegree = std::numeric_limits<unsigned short>::max();
+    for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) {
+      auto Degree = nodes->GetnPoint(iPoint);
+      if (!inQueue[iPoint] && Degree < MinDegree) {
+        MinDegree = Degree;
+        AddPoint = iPoint;
+      }
+    }
+    if (AddPoint == nPoint) {
+      SU2_MPI::Error("RCM ordering failed", CURRENT_FUNCTION);
+    }
 
-  Result.push_back(AddPoint); inQueue[AddPoint] = true;
+    /*--- Add the node in the first free position. ---*/
 
-  /*--- Loop until reorganize all the nodes ---*/
+    Result.push_back(AddPoint);
+    inQueue[AddPoint] = true;
 
-  do {
+    /*--- Loop until reorganize all the nodes ---*/
 
-    /*--- Add to the queue all the nodes adjacent in the increasing
-     order of their degree, checking if the element is already
-     in the Queue. ---*/
+    while (!Queue.empty()) {
 
-    AuxQueue.clear();
-    for (auto iNode = 0u; iNode < nodes->GetnPoint(AddPoint); iNode++) {
-      auto AdjPoint = nodes->GetPoint(AddPoint, iNode);
-      if ((!inQueue[AdjPoint]) && (AdjPoint < nPointDomain)) {
-        AuxQueue.push_back(AdjPoint);
-      }
-    }
+      /*--- Extract the first node from the queue and add it in the first free
+       position. ---*/
 
-    if (!AuxQueue.empty()) {
+      AddPoint = Queue.front();
+      Result.push_back(AddPoint);
+      Queue.pop();
+
+      /*--- Add to the queue all the nodes adjacent in the increasing
+       order of their degree, checking if the element is already
+       in the Queue. ---*/
+
+      AuxQueue.clear();
+      for (auto iNode = 0u; iNode < nodes->GetnPoint(AddPoint); iNode++) {
+        auto AdjPoint = nodes->GetPoint(AddPoint, iNode);
+        if (!inQueue[AdjPoint] && (AdjPoint < nPointDomain)) {
+          AuxQueue.push_back(AdjPoint);
+        }
+      }
+      if (AuxQueue.empty()) continue;
 
       /*--- Sort the auxiliar queue based on the number of neighbors ---*/
 
@@ -4800,24 +4815,15 @@ void CPhysicalGeometry::SetRCM_Ordering(CConfig *config) {
         Queue.push(iPoint);
         inQueue[iPoint] = true;
       }
-
     }
-
-    /*--- Extract the first node from the queue and add it in the first free
-     position. ---*/
-
-    if (!Queue.empty()) {
-      AddPoint = Queue.front();
-      Result.push_back(AddPoint);
-      Queue.pop();
-    }
-
-  } while (!Queue.empty());
+  }
 
   /*--- Check that all the points have been added ---*/
 
   for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) {
-    if (inQueue[iPoint] == false) Result.push_back(iPoint);
+    if (!inQueue[iPoint]) {
+      SU2_MPI::Error("RCM ordering failed", CURRENT_FUNCTION);
+    }
   }
 
   reverse(Result.begin(), Result.end());

SU2_CFD/include/solvers/CEulerSolver.hpp

@@ -279,13 +279,10 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
   void SetReferenceValues(const CConfig& config) final;
 
 public:
-  /*!
-   * \brief Constructor of the class.
-   */
-  CEulerSolver() : CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMPRESSIBLE>() {}
+  CEulerSolver() = delete;
 
   /*!
-   * \overload Main constructor of this class.
+   * \brief Main constructor of this class.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    * \param[in] iMesh - Grid level.

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -44,6 +44,9 @@ class CFVMFlowSolverBase : public CSolver {
   }
 
  protected:
+  /*!< \brief Object with indices of primitive variables. */
+  const typename VariableType::template CIndices<unsigned short> prim_idx;
+
   static constexpr size_t MAXNDIM = 3; /*!< \brief Max number of space dimensions, used in some static arrays. */
   static constexpr size_t MAXNVAR = VariableType::MAXNVAR; /*!< \brief Max number of variables, for static arrays. */
 
@@ -226,9 +229,10 @@ class CFVMFlowSolverBase : public CSolver {
   inline CVariable* GetBaseClassPointerToNodes() final { return nodes; }
 
   /*!
-   * \brief Default constructor, this class is not directly instantiable.
+   * \brief Protected constructor, this class is not directly instantiable.
    */
-  CFVMFlowSolverBase() : CSolver() {}
+  CFVMFlowSolverBase(const CGeometry& geometry, const CConfig& config)
+    : CSolver(), prim_idx(geometry.GetnDim(), config.GetnSpecies()) {}
 
   /*!
    * \brief Set reference values for pressure, forces, etc., e.g. "AeroCoeffForceRef".

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2468,7 +2468,6 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
-  unsigned short T_INDEX = 0, TVE_INDEX = 0, VEL_INDEX = 0;
   su2double Viscosity = 0.0, Area, Density = 0.0, GradTemperature = 0.0, WallDistMod, FrictionVel,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
             thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor;
@@ -2477,27 +2476,19 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
   string Marker_Tag, Monitoring_Tag;
 
-  su2double Alpha = config->GetAoA() * PI_NUMBER / 180.0;
-  su2double Beta = config->GetAoS() * PI_NUMBER / 180.0;
-  su2double RefLength = config->GetRefLength();
-  su2double RefHeatFlux = config->GetHeat_Flux_Ref();
-  su2double Gas_Constant = config->GetGas_ConstantND();
+  const su2double Alpha = config->GetAoA() * PI_NUMBER / 180.0;
+  const su2double Beta = config->GetAoS() * PI_NUMBER / 180.0;
+  const su2double RefLength = config->GetRefLength();
+  const su2double RefHeatFlux = config->GetHeat_Flux_Ref();
+  const su2double Gas_Constant = config->GetGas_ConstantND();
   auto Origin = config->GetRefOriginMoment(0);
 
-  su2double Prandtl_Lam = config->GetPrandtl_Lam();
-  bool energy = config->GetEnergy_Equation();
-  bool QCR = config->GetQCR();
-  bool axisymmetric = config->GetAxisymmetric();
-  bool roughwall = (config->GetnRoughWall() > 0);
-  bool nemo = config->GetNEMOProblem();
-
-  /*--- Get the locations of the primitive variables for NEMO ---*/
-  if (nemo) {
-    unsigned short nSpecies = config->GetnSpecies();
-    T_INDEX       = nSpecies;
-    TVE_INDEX     = nSpecies+1;
-    VEL_INDEX     = nSpecies+2;
-  }
+  const su2double Prandtl_Lam = config->GetPrandtl_Lam();
+  const bool energy = config->GetEnergy_Equation();
+  const bool QCR = config->GetQCR();
+  const bool axisymmetric = config->GetAxisymmetric();
+  const bool roughwall = (config->GetnRoughWall() > 0);
+  const bool nemo = config->GetNEMOProblem();
 
   const su2double factor = 1.0 / AeroCoeffForceRef;
   const su2double factorFric = config->GetRefArea() * factor;
@@ -2559,15 +2550,9 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
       for (iDim = 0; iDim < nDim; iDim++) {
         for (jDim = 0; jDim < nDim; jDim++) {
-          if (!nemo) Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint, iDim + 1, jDim);
-          else       Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint, iDim + VEL_INDEX, jDim);
+          Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Velocity() + iDim, jDim);
         }
-
-        /// TODO: Move the temperature index logic to a function.
-
-        if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) Grad_Temp[iDim] = nodes->GetGradient_Primitive(iPoint, 0, iDim);
-
-        if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE) Grad_Temp[iDim] = nodes->GetGradient_Primitive(iPoint, nDim + 1, iDim);
+        Grad_Temp[iDim] = nodes->GetGradient_Primitive(iPoint, prim_idx.Temperature(), iDim);
       }
 
       Viscosity = nodes->GetLaminarViscosity(iPoint);
@@ -2660,8 +2645,8 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
         const auto& thermal_conductivity_ve = nodes->GetThermalConductivity_ve(iPoint);
         const auto& Grad_PrimVar            = nodes->GetGradient_Primitive(iPoint);
 
-        su2double dTn   = GeometryToolbox::DotProduct(nDim, Grad_PrimVar[T_INDEX], UnitNormal);
-        su2double dTven = GeometryToolbox::DotProduct(nDim, Grad_PrimVar[TVE_INDEX], UnitNormal);
+        su2double dTn   = GeometryToolbox::DotProduct(nDim, Grad_PrimVar[prim_idx.Temperature()], UnitNormal);
+        su2double dTven = GeometryToolbox::DotProduct(nDim, Grad_PrimVar[prim_idx.Temperature_ve()], UnitNormal);
 
         /*--- Surface energy balance: trans-rot heat flux, vib-el heat flux,
         enthalpy transport due to mass diffusion ---*/

SU2_CFD/include/solvers/CIncEulerSolver.hpp

@@ -125,13 +125,10 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
   void SetReferenceValues(const CConfig& config) final;
 
 public:
-  /*!
-   * \brief Constructor of the class.
-   */
-  CIncEulerSolver() : CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME::INCOMPRESSIBLE>() {}
+  CIncEulerSolver() = delete;
 
   /*!
-   * \overload
+   * \brief Constructor of the class.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    * \param[in] iMesh - Grid level.

SU2_CFD/include/solvers/CIncNSSolver.hpp

@@ -98,11 +98,6 @@ class CIncNSSolver final : public CIncEulerSolver {
 public:
   /*!
    * \brief Constructor of the class.
-   */
-  CIncNSSolver() : CIncEulerSolver() {}
-
-  /*!
-   * \overload
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */

SU2_CFD/include/solvers/CNEMOEulerSolver.hpp

@@ -91,17 +91,13 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, ENUM_REGI
   void SetReferenceValues(const CConfig& config) final;
 
 public:
+  CNEMOEulerSolver() = delete;
 
   /*!
-   * \brief Constructor of the class.
+   * \brief Contructor of the class.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
    */
-  CNEMOEulerSolver() : CFVMFlowSolverBase<CNEMOEulerVariable, ENUM_REGIME::COMPRESSIBLE>() {}
-
-  /*!
-     * \overload
-     * \param[in] geometry - Geometrical definition of the problem.
-     * \param[in] config - Definition of the particular problem.
-     */
   CNEMOEulerSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh, const bool navier_stokes = false);
 
   /*!

SU2_CFD/include/variables/CEulerVariable.hpp

@@ -55,6 +55,9 @@ class CEulerVariable : public CFlowVariable {
     inline IndexType EddyViscosity() const { return nDim+6; }
     inline IndexType ThermalConductivity() const { return nDim+7; }
     inline IndexType CpTotal() const { return nDim+8; }
+
+    /*--- For compatible interface with NEMO. ---*/
+    inline IndexType Temperature_ve() const { return std::numeric_limits<IndexType>::max(); }
   };
 
  protected:

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -55,6 +55,9 @@ class CIncEulerVariable : public CFlowVariable {
     inline IndexType ThermalConductivity() const { return nDim+6; }
     inline IndexType CpTotal() const { return nDim+7; }
     inline IndexType CvTotal() const { return nDim+8; }
+
+    /*--- For compatible interface with NEMO. ---*/
+    inline IndexType Temperature_ve() const { return std::numeric_limits<IndexType>::max(); }
   };
 
  protected:

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -38,7 +38,7 @@
 
 CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
                            unsigned short iMesh, const bool navier_stokes) :
-  CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMPRESSIBLE>() {
+  CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMPRESSIBLE>(*geometry, *config) {
 
   /*--- Based on the navier_stokes boolean, determine if this constructor is
    *    being called by itself, or by its derived class CNSSolver. ---*/