Apply more changes from code review.

Create a common template CFEASolverBase for solvers using the internal finite elements.

Author: thomasdick

Common/include/CConfig.hpp

@@ -446,7 +446,7 @@ class CConfig {
   bool SmoothSepDim;                /*!< \brief Flag for enabling separated calculation for every dimension. */
   bool SmoothOnSurface;             /*!< \brief Flag for assembling the system only on the surface. */
   bool SmoothDirichletSurfaceBound; /*!< \brief Flag for using zero Dirichlet boundary in the surface case. */
-  unsigned short SmoothNumMode;     /*!< \brief temporary flag for some debuging stuff */
+  ENUM_SOBOLEV_MODUS SmoothNumMode; /*!< \brief The mode in which the Sobolev smoothing solver is applied. */
 
   unsigned short  Kind_Grad_Linear_Solver,  /*!< Numerical method to smooth the gradient */
   Kind_Grad_Linear_Solver_Prec;             /*!< \brief Preconditioner of the linear solver. */
@@ -9616,7 +9616,7 @@ class CConfig {
    * \brief The modus of operation for the Sobolev solver
    * \return returns on what level we operate
    */
-  unsigned short GetSobMode(void) const { return SmoothNumMode; }
+  ENUM_SOBOLEV_MODUS GetSobMode(void) const { return SmoothNumMode; }
 
   /*!
    * \brief Get the name of the file with the hessian of the objective function.

Common/include/option_structure.hpp

@@ -2367,18 +2367,18 @@ enum class LINEAR_SOLVER_MODE {
 /*!
  * \brief mode of operation for the sobolev smoothing solver.
  */
-enum ENUM_SOBOLEV_MODUS {
-  NO_MODUS,             /*!< \brief Default option if none is choosen. */
+enum class ENUM_SOBOLEV_MODUS {
+  NONE,                 /*!< \brief Default option if none is choosen. */
   PARAM_LEVEL_COMPLETE, /*!< \brief Operate on parameter level. */
   MESH_LEVEL,           /*!< \brief Operate on mesh level. */
-  DEBUG,                /*!< \brief Special flag for debugging. */
-  ONLY_GRAD,            /*!< \brief Flag for OneShot to only compute the original gradient. */
+  ONLY_GRAD,            /*!< \brief Flag to only compute the original gradient. */
 };
 static const MapType<std::string, ENUM_SOBOLEV_MODUS> Sobolev_Modus_Map = {
-    MakePair("NONE", ENUM_SOBOLEV_MODUS::NO_MODUS)
-        MakePair("PARAM_LEVEL_COMPLETE", ENUM_SOBOLEV_MODUS::PARAM_LEVEL_COMPLETE)
-            MakePair("MESH_LEVEL", ENUM_SOBOLEV_MODUS::MESH_LEVEL) MakePair("DEBUG", ENUM_SOBOLEV_MODUS::DEBUG)
-                MakePair("ONLY_GRADIENT", ENUM_SOBOLEV_MODUS::ONLY_GRAD)};
+  MakePair("NONE",                 ENUM_SOBOLEV_MODUS::NONE)
+  MakePair("PARAM_LEVEL_COMPLETE", ENUM_SOBOLEV_MODUS::PARAM_LEVEL_COMPLETE)
+  MakePair("MESH_LEVEL",           ENUM_SOBOLEV_MODUS::MESH_LEVEL)
+  MakePair("ONLY_GRADIENT",        ENUM_SOBOLEV_MODUS::ONLY_GRAD)
+};
 
 #undef MakePair
 /* END_CONFIG_ENUMS */

Common/src/CConfig.cpp

@@ -2012,7 +2012,7 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Switch to activate zero Dirichlet boundary for surface mode */
   addBoolOption("DIRICHLET_SURFACE_BOUNDARY", SmoothDirichletSurfaceBound, false);
   /* DESCRIPTION: Switch to activate the debbuging modus */
-  addEnumOption("SOBOLEV_MODE", SmoothNumMode, Sobolev_Modus_Map, ENUM_SOBOLEV_MODUS::NO_MODUS);
+  addEnumOption("SOBOLEV_MODE", SmoothNumMode, Sobolev_Modus_Map, ENUM_SOBOLEV_MODUS::NONE);
   /*!\brief HESS_OBJFUNC_FILENAME
    *  \n DESCRIPTION: Output filename for the Sobolev Hessian approximation.  \ingroup Config*/
   addStringOption("HESS_OBJFUNC_FILENAME", ObjFunc_Hess_FileName, string("of_hess.dat"));

Common/src/geometry/elements/CTETRA1.cpp

@@ -26,7 +26,7 @@
  */
 
 #include "../../../include/geometry/elements/CElement.hpp"
-
+#include "../../../include/toolboxes/geometry_toolbox.hpp"
 
 CTETRA1::CTETRA1() : CElementWithKnownSizes<NGAUSS,NNODE,NDIM>() {
 
@@ -85,11 +85,8 @@ su2double CTETRA1::ComputeVolume(const FrameType mode) const {
     r3[iDim] = Coord[3][iDim] - Coord[0][iDim];
   }
 
-  CrossProduct[0] = (r1[1]*r2[2] - r1[2]*r2[1])*r3[0];
-  CrossProduct[1] = (r1[2]*r2[0] - r1[0]*r2[2])*r3[1];
-  CrossProduct[2] = (r1[0]*r2[1] - r1[1]*r2[0])*r3[2];
-
-  Volume = fabs(CrossProduct[0] + CrossProduct[1] + CrossProduct[2])/6.0;
+  GeometryToolbox::CrossProduct(r1, r2, CrossProduct);
+  Volume = fabs(GeometryToolbox::DotProduct(3, CrossProduct, r3))/6.0;
 
   return Volume;
 

Common/src/geometry/elements/CTETRA4.cpp

@@ -26,7 +26,7 @@
  */
 
 #include "../../../include/geometry/elements/CElement.hpp"
-
+#include "../../../include/toolboxes/geometry_toolbox.hpp"
 
 CTETRA4::CTETRA4() : CElementWithKnownSizes<NGAUSS,NNODE,NDIM>() {
 
@@ -90,11 +90,8 @@ su2double CTETRA4::ComputeVolume(const FrameType mode) const {
     r3[iDim] = Coord[3][iDim] - Coord[0][iDim];
   }
 
-  CrossProduct[0] = (r1[1]*r2[2] - r1[2]*r2[1])*r3[0];
-  CrossProduct[1] = (r1[2]*r2[0] - r1[0]*r2[2])*r3[1];
-  CrossProduct[2] = (r1[0]*r2[1] - r1[1]*r2[0])*r3[2];
-
-  Volume = fabs(CrossProduct[0] + CrossProduct[1] + CrossProduct[2])/6.0;
+  GeometryToolbox::CrossProduct(r1, r2, CrossProduct);
+  Volume = fabs(GeometryToolbox::DotProduct(3, CrossProduct, r3))/6.0;
 
   return Volume;
 

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -27,23 +27,15 @@
 
 #pragma once
 
-#include "CSolver.hpp"
-#include "../../../Common/include/geometry/elements/CElement.hpp"
-#include "../../../Common/include/parallelization/omp_structure.hpp"
+#include "CFEASolverBase.hpp"
 
 /*!
  * \class CFEASolver
  * \brief Main class for defining a FEM solver for elastic structural problems.
  * \author R. Sanchez.
  */
-class CFEASolver : public CSolver {
+class CFEASolver : public CFEASolverBase {
 protected:
-  enum : size_t {MAXNNODE_2D = 4};
-  enum : size_t {MAXNNODE_3D = 8};
-  enum : size_t {MAXNVAR = 3};
-  enum : size_t {MAXNDIM = 3};
-  enum : size_t {OMP_MIN_SIZE = 32};
-  enum : size_t {OMP_MAX_SIZE = 512};
 
   unsigned long omp_chunk_size;     /*!< \brief Chunk size used in light point loops. */
 
@@ -91,7 +83,6 @@ class CFEASolver : public CSolver {
   CSysMatrix<su2double> MassMatrix;
 #endif
 
-  CElement*** element_container = nullptr;  /*!< \brief Vector which the define the finite element structure for each problem. */
   CProperty** element_properties = nullptr; /*!< \brief Vector which stores the properties of each element */
 
 #ifdef HAVE_OMP
@@ -108,22 +99,6 @@ class CFEASolver : public CSolver {
   bool topol_filter_applied;   /*!< \brief True if density filtering has been performed. */
   bool initial_calc = true;    /*!< \brief Becomes false after first call to Preprocessing. */
 
-  unsigned long nElement;      /*!< \brief Number of elements. */
-
-  /*--- Extra vertices for row/column elimination, see Set_VertexEliminationSchedule. ---*/
-  vector<unsigned long> ExtraVerticesToEliminate;
-
-  /*!
-   * \brief The highest level in the variable hierarchy this solver can safely use,
-   * CVariable is the common denominator between the FEA and Mesh deformation variables.
-   */
-  CVariable* nodes = nullptr;
-
-  /*!
-   * \brief Return nodes to allow CSolver::base_nodes to be set.
-   */
-  inline CVariable* GetBaseClassPointerToNodes() override { return nodes; }
-
   /*!
    * \brief Get the element container index and number of nodes of a given VTK type.
    * \param[in] VTK_Type - Type of element.

SU2_CFD/include/solvers/CFEASolverBase.hpp

@@ -0,0 +1,89 @@
+/*!
+ * \file CFEASolver.hpp
+ * \brief Base class template for all FEA solvers using the SU2 internal finite elements.
+ * \author T. Dick
+ * \version 7.2.1 "Blackbird"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2021, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "CSolver.hpp"
+#include "../../../Common/include/geometry/elements/CElement.hpp"
+#include "../../../Common/include/parallelization/omp_structure.hpp"
+
+class CFEASolverBase : public CSolver {
+ public:
+  enum : size_t {MAXNNODE_2D = 4};
+  enum : size_t {MAXNNODE_3D = 8};
+  enum : size_t {MAXNDIM = 3};
+  enum : size_t {MAXNVAR = 3};
+  enum : size_t {OMP_MIN_SIZE = 32};
+  enum : size_t {OMP_MAX_SIZE = 512};
+
+ protected:
+
+  CElement*** element_container = nullptr;  /*!< \brief Vector which the define the finite element structure for each problem. */
+  unsigned long nElement;                   /*!< \brief Number of elements. */
+
+  vector<unsigned long> ExtraVerticesToEliminate; /*!< Extra vertices for row/column elimination, see CommunicateExtraEliminationVertices. */
+
+  /*!
+   * \brief The highest level in the variable hierarchy all derived solvers can safely use,
+   * CVariable is the common denominator between the FEA, Mesh deformation, and smoothing variables.
+   */
+  CVariable* nodes = nullptr;
+
+  /*!
+   * \brief Return nodes to allow CSolver::base_nodes to be set.
+   */
+  inline CVariable* GetBaseClassPointerToNodes() override { return nodes; }
+
+ public:
+  /*!
+   * \brief Constructor of the class.
+   * \param[in] mesh_deform_mode - mode of operation for the linear solver.
+   */
+  CFEASolverBase(LINEAR_SOLVER_MODE mesh_deform_mode = LINEAR_SOLVER_MODE::STANDARD);
+
+  /*!
+   * \overload
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] mesh_deform_mode - mode of operation for the linear solver.
+   */
+  CFEASolverBase(CGeometry *geometry, CConfig *config, LINEAR_SOLVER_MODE mesh_deform_mode = LINEAR_SOLVER_MODE::STANDARD);
+
+  /*!
+   * \brief Destructor of the class.
+   */
+  ~CFEASolverBase();
+
+  /*!
+   * \brief Communicate extra vertices for elimination in the linear system.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] myPoints - List of global point indeces to eliminate.
+   * \param[in] val_coord - Location (x, y, z) of the max residual point.
+   */
+  void CommunicateExtraEliminationVertices(const CGeometry* geometry, vector<unsigned long>& myPoints);
+
+};

SU2_CFD/include/solvers/CGradientSmoothingSolver.hpp

@@ -27,21 +27,17 @@
 
 #pragma once
 
-#include "../../../Common/include/geometry/elements/CElement.hpp"
 #include "../../../Common/include/linear_algebra/CMatrixVectorProduct.hpp"
 #include "../../../Common/include/toolboxes/CSquareMatrixCM.hpp"
-#include "../../../SU2_CFD/include/solvers/CSolver.hpp"
+#include "../../../SU2_CFD/include/solvers/CFEASolverBase.hpp"
 
 /*! \class CGradientSmoothingSolver
  *  \brief Main class for defining a gradient smoothing.
  *  \author T. Dick.
  *  \date March 25, 2019.
  */
-class CGradientSmoothingSolver : public CSolver {
+class CGradientSmoothingSolver final : public CFEASolverBase {
 public:
-  enum : size_t {MAXNNODE_2D = 4};
-  enum : size_t {MAXNNODE_3D = 8};
-  enum : size_t {MAXNDIM = 3};
 
 /** Introduction of a new alias for the data type to allow compilation with forward mode.
  *
@@ -57,32 +53,18 @@ class CGradientSmoothingSolver : public CSolver {
   using MatrixType = C2DContainer<unsigned long, su2double, StorageType::RowMajor, 64, DynamicSize, DynamicSize>;
 
  private:
-  unsigned long nElement;
-  CElement*** element_container  = nullptr;  /*!< \brief Container which stores the element information. */
+  unsigned int curDim;                       /*!< \brief If we separate dimensions this tells us in what dimension we currently are. */
 
-  unsigned int dir;                          /*!< \brief If we separate dimensions this tells us in what dimension we currently are. */
-
-  CSysVector<su2double> auxVec;              /*!< \brief Auxiliar vectors for output and debugging */
   CSysVector<su2double> activeCoord;         /*!< \brief Auxiliar vector to keep the indeces of geometry->vertex->Coord */
 
   CSysVector<su2matvecscalar> helperVecIn;   /*!< \brief Helper vectors for projection and matrix vector product (must be su2mixedfloat) */
   CSysVector<su2matvecscalar> helperVecOut;  /*!< \brief Helper vectors for projection and matrix vector product (must be su2mixedfloat) */
-  CSysVector<su2matvecscalar> helperVecAux;     /*!< \brief Helper vectors for matrix vector product if working on surface (smaller dim) */
-
-  CVariable* nodes = nullptr;                /*!< \brief The highest level in the variable hierarchy this solver can safely use. */
+  CSysVector<su2matvecscalar> helperVecAux;  /*!< \brief Helper vectors for matrix vector product if working on surface (smaller dim) */
 
   std::vector<su2double> deltaP;             /*!< \brief The smoothed gradient with respect to the design variables. */
 
   CSquareMatrixCM hessian;                   /*!< \brief The approximated Hessian with respect to the design variables. */
 
-  /*--- Extra vertices for row/column elimination, see Set_VertexEliminationSchedule. ---*/
-  vector<unsigned long> ExtraVerticesToEliminate;
-
-  /*!
-   * \brief Return nodes to allow CSolver::base_nodes to be set.
-   */
-  inline CVariable* GetBaseClassPointerToNodes() override { return nodes; }
-
  public:
   /*!
    * \brief Constructor of the class.
@@ -244,26 +226,24 @@ class CGradientSmoothingSolver : public CSolver {
   /*!
    * \brief Write the content of sensitivity in the nodes to the sensitivity in the geometry
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
    */
-  void WriteSens2Geometry(CGeometry* geometry, const CConfig* config);
+  void WriteSensToGeometry(CGeometry* geometry) const override;
 
   /*!
    * \brief Read the sensitivity for the geometry into the nodes
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
    */
-  void ReadSens2Geometry(CGeometry* geometry, const CConfig* config);
+  void ReadSensFromGeometry(const CGeometry *geometry) override;
 
   /*!
    * \brief Copy sensitivities from a vector into the geometry
    */
-  void WriteVector2Geometry(CGeometry* geometry, const CConfig* config, CSysVector<su2matvecscalar>& vector);
+  void WriteVectorToGeometry(CGeometry* geometry, const CSysVector<su2matvecscalar>& vector) const;
 
   /*!
    * \brief Copy sensitivities from the geometry into a vector
    */
-  void ReadVector2Geometry(CGeometry* geometry, const CConfig* config, CSysVector<su2matvecscalar>& vector);
+  void ReadVectorToGeometry(const CGeometry* geometry, CSysVector<su2matvecscalar>& vector);
 
   /*!
    * \brief Get the value of the reference coordinate to set on the element structure.
@@ -309,4 +289,14 @@ class CGradientSmoothingSolver : public CSolver {
     }
   }
 
+private:
+
+  inline void SetCurrentDim(unsigned int iDim) {
+    curDim=iDim;
+  }
+
+  inline const unsigned int& GetCurrentDim() const {
+    return curDim;
+  }
+
 };

SU2_CFD/include/solvers/CSolver.hpp

@@ -38,6 +38,8 @@
 #include <set>
 #include <stdlib.h>
 #include <stdio.h>
+#include <unordered_set>
+#include <unordered_map>
 
 #include "../fluid/CFluidModel.hpp"
 #include "../task_definition.hpp"
@@ -4200,16 +4202,15 @@ class CSolver {
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
    */
-  virtual void ReadSens2Geometry(CGeometry* geometry, const CConfig* config) {}
+  virtual void ReadSensFromGeometry(const CGeometry* geometry) {}
 
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  virtual void WriteSens2Geometry(CGeometry* geometry, const CConfig* config) {}
+  virtual void WriteSensToGeometry(CGeometry* geometry) const {}
 
   /*!
    * \brief Routine that sets the flag controlling implicit treatment for periodic BCs.

SU2_CFD/include/variables/CSobolevSmoothingVariable.hpp

@@ -34,7 +34,6 @@ class CSobolevSmoothingVariable final : public CVariable {
   MatrixType Sensitivity;  /*!< Vector holding the derivative of target functional with respect to the coordinates at this node */
 
   bool* boundary_vertex;  /*!< \brief Stores if a point belongs to the boundary of a boundary. */
-  unsigned long nBoundPoints;
 
   /*!
    * \brief Constructor of the class.
@@ -72,9 +71,4 @@ class CSobolevSmoothingVariable final : public CVariable {
    * \brief return wether a point is a boundary of a boundary
    */
   bool GetIsBoundaryPoint(unsigned long iPoint) const override;
-
-  /*!
-   * \brief return the number of marked points
-   */
-  unsigned int GetNBoundPoints() const override;
 };