Merge pull request #1372 from su2code/improve_mg

Cleanup and document the multigrid agglomeration strategy

Common/include/geometry/CMultiGridGeometry.hpp

@@ -36,23 +36,7 @@
  * \author F. Palacios
  */
 class CMultiGridGeometry final : public CGeometry {
-
-public:
-  /*--- This is to suppress Woverloaded-virtual, omitting it has no negative impact. ---*/
-  using CGeometry::SetVertex;
-  using CGeometry::SetControlVolume;
-  using CGeometry::SetBoundControlVolume;
-  using CGeometry::SetPoint_Connectivity;
-
-  /*!
-   * \brief Constructor of the class.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] iMesh - Level of the multigrid.
-   * \param[in] iZone - Current zone in the mesh.
-   */
-  CMultiGridGeometry(CGeometry **geometry, CConfig *config_container, unsigned short iMesh);
-
+private:
   /*!
    * \brief Determine if a CVPoint van be agglomerated, if it have the same marker point as the seed.
    * \param[in] CVPoint - Control volume to be agglomerated.
@@ -61,101 +45,159 @@ class CMultiGridGeometry final : public CGeometry {
    * \param[in] config - Definition of the particular problem.
    * \return <code>TRUE</code> or <code>FALSE</code> depending if the control volume can be agglomerated.
    */
-  bool SetBoundAgglomeration(unsigned long CVPoint, short marker_seed, CGeometry *fine_grid, CConfig *config);
+  bool SetBoundAgglomeration(unsigned long CVPoint, short marker_seed, const CGeometry *fine_grid,
+                             const CConfig *config) const;
 
   /*!
    * \brief Determine if a can be agglomerated using geometrical criteria.
    * \param[in] iPoint - Seed point.
    * \param[in] fine_grid - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  bool GeometricalCheck(unsigned long iPoint, CGeometry *fine_grid, CConfig *config);
+  bool GeometricalCheck(unsigned long iPoint, const CGeometry *fine_grid, const CConfig *config) const;
 
   /*!
    * \brief Determine if a CVPoint van be agglomerated, if it have the same marker point as the seed.
-   * \param[in] Suitable_Indirect_Neighbors - List of Indirect Neighbours that can be agglomerated.
+   * \param[out] Suitable_Indirect_Neighbors - List of Indirect Neighbours that can be agglomerated.
    * \param[in] iPoint - Seed point.
    * \param[in] Index_CoarseCV - Index of agglomerated point.
    * \param[in] fine_grid - Geometrical definition of the problem.
    */
-  void SetSuitableNeighbors(vector<unsigned long> *Suitable_Indirect_Neighbors, unsigned long iPoint,
-                            unsigned long Index_CoarseCV, CGeometry *fine_grid);
+  void SetSuitableNeighbors(vector<unsigned long>& Suitable_Indirect_Neighbors, unsigned long iPoint,
+                            unsigned long Index_CoarseCV, const CGeometry *fine_grid) const;
 
   /*!
-   * \brief Set boundary vertex.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
+   * \brief Set a representative wall value of the agglomerated control volumes on a particular boundary marker.
+   * \param[in] fine_grid - Geometrical definition of the problem.
+   * \param[in] val_marker - Index of the boundary marker.
+   * \param[in] wall_quantity - Object with methods Get(iVertex_fine) and Set(iVertex_coarse, val).
    */
-  void SetVertex(CGeometry *geometry, CConfig *config) override;
+  template <class T>
+  void SetMultiGridWallQuantity(const CGeometry *fine_grid, unsigned short val_marker, T& wall_quantity) {
+
+    for (auto iVertex = 0ul; iVertex < nVertex[val_marker]; iVertex++) {
+      const auto Point_Coarse = vertex[val_marker][iVertex]->GetNode();
+
+      if (!nodes->GetDomain(Point_Coarse)) continue;
+
+      su2double Area_Parent = 0.0;
+
+      /*--- Compute area parent by taking into account only volumes that are on the marker. ---*/
+      for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
+        const auto Point_Fine = nodes->GetChildren_CV(Point_Coarse, iChildren);
+        const auto isVertex = fine_grid->nodes->GetDomain(Point_Fine) &&
+                             (fine_grid->nodes->GetVertex(Point_Fine, val_marker) != -1);
+        if (isVertex) {
+          Area_Parent += fine_grid->nodes->GetVolume(Point_Fine);
+        }
+      }
+
+      su2double Quantity_Coarse = 0.0;
+
+      /*--- Loop again to average coarser value. ---*/
+      for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(Point_Coarse); iChildren++) {
+        const auto Point_Fine = nodes->GetChildren_CV(Point_Coarse, iChildren);
+        const auto isVertex = fine_grid->nodes->GetDomain(Point_Fine) &&
+                             (fine_grid->nodes->GetVertex(Point_Fine, val_marker) != -1);
+        if (isVertex) {
+          const auto Vertex_Fine = fine_grid->nodes->GetVertex(Point_Fine, val_marker);
+          const auto Area_Children = fine_grid->nodes->GetVolume(Point_Fine);
+          Quantity_Coarse += wall_quantity.Get(Vertex_Fine) * Area_Children / Area_Parent;
+        }
+      }
+
+      /*--- Set the value at the coarse level. ---*/
+      wall_quantity.Set(iVertex, Quantity_Coarse);
+    }
+
+  }
+
+public:
+  /*--- This is to suppress Woverloaded-virtual, omitting it has no negative impact. ---*/
+  using CGeometry::SetVertex;
+  using CGeometry::SetControlVolume;
+  using CGeometry::SetBoundControlVolume;
+  using CGeometry::SetPoint_Connectivity;
 
   /*!
-   * \brief Set points which surround a point.
-   * \param[in] geometry - Geometrical definition of the problem.
+   * \brief Constructor of the class.
+   * \param[in] fine_grid - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Level of the multigrid.
    */
-  void SetPoint_Connectivity(CGeometry *geometry) override;
+  CMultiGridGeometry(CGeometry *fine_grid, CConfig *config, unsigned short iMesh);
 
   /*!
-   * \brief Set the edge structure of the agglomerated control volume.
+   * \brief Set boundary vertex.
+   * \param[in] fine_grid - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] action - Allocate or not the new elements.
    */
-  void SetControlVolume(CConfig *config, CGeometry *geometry, unsigned short action) override;
+  void SetVertex(const CGeometry *fine_grid, const CConfig *config) override;
 
   /*!
-   * \brief Mach the near field boundary condition.
-   * \param[in] config - Definition of the particular problem.
+   * \brief Set points which surround a point.
+   * \param[in] fine_grid - Geometrical definition of the child grid.
    */
-  void MatchActuator_Disk(CConfig *config) override;
+  void SetPoint_Connectivity(const CGeometry *fine_grid) override;
 
   /*!
-   * \brief Mach the periodic boundary conditions.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] val_periodic - Index of the first periodic face in a pair.
+   * \brief Set the edge structure of the agglomerated control volume.
+   * \param[in] fine_grid - Geometrical definition of the problem.
+   * \param[in] action - Allocate or not the new elements.
    */
-  void MatchPeriodic(CConfig *config, unsigned short val_periodic) override;
+  void SetControlVolume(const CGeometry *fine_grid, unsigned short action) override;
 
   /*!
    * \brief Set boundary vertex structure of the agglomerated control volume.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] fine_grid - Geometrical definition of the problem.
    * \param[in] action - Allocate or not the new elements.
    */
-  void SetBoundControlVolume(CConfig *config, CGeometry *geometry, unsigned short action) override;
+  void SetBoundControlVolume(const CGeometry *fine_grid, unsigned short action) override;
 
   /*!
    * \brief Set a representative coordinates of the agglomerated control volume.
-   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] fine_grid - Geometrical definition of the problem.
    */
-  void SetCoord(CGeometry *geometry) override;
+  void SetCoord(const CGeometry *fine_grid) override;
 
   /*!
-   * \brief Set a representative wall normal heat flux of the agglomerated control volume on a particular boundary marker.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] val_marker - Index of the boundary marker.
+   * \brief Set the grid velocity at each node in the coarse mesh level based
+   *        on a restriction from a finer mesh.
+   * \param[in] fine_grid - Geometry container for the finer mesh level.
    */
-  void SetMultiGridWallHeatFlux(CGeometry *geometry, unsigned short val_marker) override;
+  void SetRestricted_GridVelocity(const CGeometry *fine_grid) override;
 
   /*!
-   * \brief Set a representative wall temperature of the agglomerated control volume on a particular boundary marker.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] val_marker - Index of the boundary marker.
+   * \brief Find and store the closest neighbor to a vertex.
+   * \param[in] config - Definition of the particular problem.
    */
-  void SetMultiGridWallTemperature(CGeometry *geometry, unsigned short val_marker) override;
+  void FindNormal_Neighbor(const CConfig *config) override;
 
   /*!
-   * \brief Set the grid velocity at each node in the coarse mesh level based
-   *        on a restriction from a finer mesh.
-   * \param[in] fine_mesh - Geometry container for the finer mesh level.
+   * \brief Mach the near field boundary condition.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetRestricted_GridVelocity(CGeometry *fine_mesh, const CConfig *config) override;
+  void MatchActuator_Disk(const CConfig *config) override;
 
   /*!
-   * \brief Find and store the closest neighbor to a vertex.
+   * \brief Mach the periodic boundary conditions.
    * \param[in] config - Definition of the particular problem.
+   * \param[in] val_periodic - Index of the first periodic face in a pair.
    */
-  void FindNormal_Neighbor(CConfig *config) override;
+  void MatchPeriodic(const CConfig *config, unsigned short val_periodic) override;
 
-};
+  /*!
+   * \brief Set a representative wall normal heat flux of the agglomerated control volume on a particular boundary marker.
+   * \param[in] fine_grid - Geometrical definition of the problem.
+   * \param[in] val_marker - Index of the boundary marker.
+   */
+  void SetMultiGridWallHeatFlux(const CGeometry *fine_grid, unsigned short val_marker) override;
 
+  /*!
+   * \brief Set a representative wall temperature of the agglomerated control volume on a particular boundary marker.
+   * \param[in] fine_grid - Geometrical definition of the problem.
+   * \param[in] val_marker - Index of the boundary marker.
+   */
+  void SetMultiGridWallTemperature(const CGeometry *fine_grid, unsigned short val_marker) override;
+
+};

Common/include/geometry/CPhysicalGeometry.hpp

@@ -400,7 +400,7 @@ class CPhysicalGeometry final : public CGeometry {
    * \brief Set boundary vertex.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetVertex(CConfig *config) override;
+  void SetVertex(const CConfig *config) override;
 
   /*!
    * \brief Set number of span wise level for turbomachinery computation.
@@ -449,14 +449,14 @@ class CPhysicalGeometry final : public CGeometry {
    * \brief Mach the near field boundary condition.
    * \param[in] config - Definition of the particular problem.
    */
-  void MatchActuator_Disk(CConfig *config) override;
+  void MatchActuator_Disk(const CConfig *config) override;
 
   /*!
    * \brief Mach the periodic boundary conditions.
    * \param[in] config - Definition of the particular problem.
    * \param[in] val_periodic - Index of the first periodic face in a pair.
    */
-  void MatchPeriodic(CConfig *config, unsigned short val_periodic) override;
+  void MatchPeriodic(const CConfig *config, unsigned short val_periodic) override;
 
   /*!
    * \brief Set boundary vertex structure of the control volume.
@@ -580,7 +580,7 @@ class CPhysicalGeometry final : public CGeometry {
    * \brief Find and store the closest neighbor to a vertex.
    * \param[in] config - Definition of the particular problem.
    */
-  void FindNormal_Neighbor(CConfig *config) override;
+  void FindNormal_Neighbor(const CConfig *config) override;
 
   /*!
    * \brief Read the sensitivity from an input file.

Common/src/fem/geometry_structure_fem_part.cpp

@@ -583,7 +583,6 @@ void CPhysicalGeometry::Read_SU2_Format_Parallel_FEM(CConfig        *config,
 
   /*--- Allocate the memory for the coordinates to be stored on this rank. ---*/
   nPoint     = nodeIDsElemLoc.size();
-  nPointNode = nPoint;
   nodes = new CPoint(nPoint, nDim);
 
   /*--- Open the grid file again and go to the position where
@@ -1146,7 +1145,6 @@ void CPhysicalGeometry::Read_CGNS_Format_Parallel_FEM(CConfig        *config,
 
   /*--- Allocate the memory for the coordinates to be stored on this rank. ---*/
   nPoint     = nodeIDsElemLoc.size();
-  nPointNode = nPoint;
   nodes = new CPoint(nPoint, nDim);
 
   /*--- Store the global ID's of the nodes in such a way that they can
@@ -1283,7 +1281,6 @@ void CPhysicalGeometry::Read_CGNS_Format_Parallel_FEM(CConfig        *config,
   /*--- Sequential mode. Create the data for the points. The global
         number of points equals the local number of points. ---*/
   nPoint     = Global_nPoint;
-  nPointNode = nPoint;
   nodes = new CPoint(nPoint, nDim);
 
   for(unsigned long i=0; i<nPoint; ++i) {

Common/src/geometry/CGeometry.cpp

@@ -1487,7 +1487,7 @@ void CGeometry::TestGeometry(void) const {
 }
 
 bool CGeometry::SegmentIntersectsPlane(const su2double *Segment_P0, const su2double *Segment_P1, su2double Variable_P0, su2double Variable_P1,
-                                                           const su2double *Plane_P0, const su2double *Plane_Normal, su2double *Intersection, su2double &Variable_Interp) {
+                                       const su2double *Plane_P0, const su2double *Plane_Normal, su2double *Intersection, su2double &Variable_Interp) {
   su2double u[3], v[3], Denominator, Numerator, Aux, ModU;
   su2double epsilon = 1E-6; // An epsilon is added to eliminate, as much as possible, the posibility of a line that intersects a point
   unsigned short iDim;
@@ -2430,8 +2430,8 @@ void CGeometry::UpdateGeometry(CGeometry **geometry_container, CConfig *config)
   for (unsigned short iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) {
     /*--- Update the control volume structures ---*/
 
-    geometry_container[iMesh]->SetControlVolume(config,geometry_container[iMesh-1], UPDATE);
-    geometry_container[iMesh]->SetBoundControlVolume(config,geometry_container[iMesh-1], UPDATE);
+    geometry_container[iMesh]->SetControlVolume(geometry_container[iMesh-1], UPDATE);
+    geometry_container[iMesh]->SetBoundControlVolume(geometry_container[iMesh-1], UPDATE);
     geometry_container[iMesh]->SetCoord(geometry_container[iMesh-1]);
 
   }