Various suggestions from code review.

 - Adding a common CFeaSolverBase class for finite elements solvers.
 - making functions private and variables constant in CGradientSmoothingSolver.
 - reworked the storage of boundary vertices in CSobolevSmoothingVariable.

Author: thomasdick

Common/include/linear_algebra/CSysMatrix.hpp

@@ -498,6 +498,42 @@ class CSysMatrix {
     AddBlock(block_i, block_j, val_block, OtherType(-1));
   }
 
+  /*!
+   * \brief Set the value of a scaled block in the sparse matrix.
+   * \note This is an templated overload for C2Dcontainer specialization su2matrix.
+   *       It assumes that MatrixType supports a member type Scalar and access operator[][].
+   *       If the template param Overwrite is false we add to the block (bij += alpha*b).
+   * \param[in] block_i - Row index.
+   * \param[in] block_j - Column index.
+   * \param[in] val_block - Block to set to A(i, j).
+   * \param[in] alpha - Scale factor.
+   */
+  template <class MatrixType, bool Overwrite = true>
+  inline void SetBlock(unsigned long block_i, unsigned long block_j, MatrixType& val_block,
+                       typename MatrixType::Scalar alpha = 1.0) {
+    auto mat_ij = GetBlock(block_i, block_j);
+    if (!mat_ij) return;
+    for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+      for (auto jVar = 0ul; jVar < nEqn; ++jVar) {
+        *mat_ij = (Overwrite ? ScalarType(0) : *mat_ij) + PassiveAssign(alpha * val_block[iVar][jVar]);
+        ++mat_ij;
+      }
+    }
+  }
+
+  /*!
+   * \brief Adds a scaled block to the sparse matrix (see SetBlock).
+   * \param[in] block_i - Row index.
+   * \param[in] block_j - Column index.
+   * \param[in] val_block - Block to add to A(i, j).
+   * \param[in] alpha - Scale factor.
+   */
+  template <class MatrixType>
+  inline void AddBlock(unsigned long block_i, unsigned long block_j, MatrixType& val_block,
+                       typename MatrixType::Scalar alpha = 1.0) {
+    SetBlock<MatrixType, false>(block_i, block_j, val_block, alpha);
+  }
+
   /*!
    * \brief Update 4 blocks ii, ij, ji, jj (add to i* sub from j*).
    * \note This method assumes an FVM-type sparse pattern.

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -99,6 +99,17 @@ class CFEASolver : public CFEASolverBase {
   bool topol_filter_applied;   /*!< \brief True if density filtering has been performed. */
   bool initial_calc = true;    /*!< \brief Becomes false after first call to Preprocessing. */
 
+  /*!
+   * \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

@@ -47,17 +47,6 @@ class CFEASolverBase : public CSolver {
 
   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.

SU2_CFD/include/solvers/CGradientSmoothingSolver.hpp

@@ -29,7 +29,8 @@
 
 #include "../../../Common/include/linear_algebra/CMatrixVectorProduct.hpp"
 #include "../../../Common/include/toolboxes/CSquareMatrixCM.hpp"
-#include "../../../SU2_CFD/include/solvers/CFEASolverBase.hpp"
+#include "../variables/CSobolevSmoothingVariable.hpp"
+#include "CFEASolverBase.hpp"
 
 /*! \class CGradientSmoothingSolver
  *  \brief Main class for defining a gradient smoothing.
@@ -65,6 +66,17 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
 
   CSquareMatrixCM hessian;                   /*!< \brief The approximated Hessian with respect to the design variables. */
 
+  /*!
+   * \brief The highest level in the variable hierarchy all derived solvers can safely use,
+   * CVariable is the common denominator between the FEA and Mesh deformationd variables.
+   */
+  CSobolevSmoothingVariable* 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.
@@ -129,12 +141,12 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
   /*!
    * \brief Set Dirichlet boundary conditions
    */
-  void BC_Dirichlet(CGeometry* geometry, const CConfig* config, unsigned int val_marker);
+  void BC_Dirichlet(const CGeometry* geometry, const CConfig* config, unsigned int val_marker);
 
   /*!
    * \brief Set Dirichlet boundary conditions for the surface solver
    */
-  void BC_Surface_Dirichlet(CGeometry* geometry, const CConfig* config, unsigned int val_marker);
+  void BC_Surface_Dirichlet(const CGeometry* geometry, const CConfig* config, unsigned int val_marker);
 
   /*!
    * \brief Call the linear systems solver
@@ -174,7 +186,7 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
   /*!
    * \brief write the DV gradient into a file
    */
-  void OutputDVGradient(CConfig *config, string out_file="delta_p.txt");
+  void OutputDVGradient(const CConfig* config, string out_file = "delta_p.txt");
 
   /*!
    * \brief Record a tape containing the parameter Jacobian.
@@ -211,18 +223,6 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
                        CSysVector<su2double> &registeredCoord,
                        CConfig *config);
 
-  /*!
-   * \brief Extract and set the geometrical sensitivity.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*) override;
-
-  /*!
-   * \brief Write the solution of the linear solver into the sensitivities of the nodes
-   */
-  void WriteSensitivity(CGeometry* geometry, const CConfig* config, unsigned long val_marker = 0);
-
   /*!
    * \brief Write the content of sensitivity in the nodes to the sensitivity in the geometry
    * \param[in] geometry - Geometrical definition of the problem.
@@ -235,6 +235,12 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
    */
   void ReadSensFromGeometry(const CGeometry *geometry) override;
 
+ private:
+  /*!
+   * \brief Write the solution of the linear solver into the sensitivities of the nodes
+   */
+  void WriteSensitivity(CGeometry* geometry, const CConfig* config, unsigned long val_marker = 0);
+
   /*!
    * \brief Copy sensitivities from a vector into the geometry
    */
@@ -257,13 +263,6 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
     return geometry->nodes->GetCoord(indexNode, iDim);
   }
 
-  /*!
-   * \brief Extract the Coordinates of the element from geometry
-   */
-  su2activematrix GetElementCoordinates(CGeometry *geometry,
-                                        std::vector<unsigned long>& indexNode,
-                                        int EL_KIND = 0);
-
   /*!
    * \brief Extra entries to eliminate in the linear system
    */
@@ -289,12 +288,17 @@ class CGradientSmoothingSolver final : public CFEASolverBase {
     }
   }
 
-private:
-
+  /*!
+   * \brief Set the current working dimension, if the seperate dimension option is set.
+   * \param[in] iDim - the dimension we are currently working in.
+   */
   inline void SetCurrentDim(unsigned int iDim) {
     curDim=iDim;
   }
 
+  /*!
+   * \brief Return the current working dimension.
+   */
   inline const unsigned int& GetCurrentDim() const {
     return curDim;
   }

SU2_CFD/include/variables/CSobolevSmoothingVariable.hpp

@@ -27,13 +27,14 @@
 
 #pragma once
 
-#include "../../../SU2_CFD/include/variables/CVariable.hpp"
+#include "../../../Common/include/containers/CVertexMap.hpp"
+#include "CVariable.hpp"
 
 class CSobolevSmoothingVariable final : public CVariable {
  public:
   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. */
+  CVertexMap<unsigned> BoundaryVertexMap; /*!< \brief Stores if a point belongs to the boundary of a boundary. */
 
   /*!
    * \brief Constructor of the class.
@@ -46,7 +47,7 @@ class CSobolevSmoothingVariable final : public CVariable {
   /*!
    * \brief Destructor of the class.
    */
-  ~CSobolevSmoothingVariable();
+  ~CSobolevSmoothingVariable() override = default;
 
   /*!
    * \brief Set the sensitivity at the node
@@ -65,10 +66,15 @@ class CSobolevSmoothingVariable final : public CVariable {
   /*!
    * \brief Mark a point as boundary of a boundary
    */
-  void MarkAsBoundaryPoint(unsigned long iPoint) override;
+  void MarkAsBoundaryPoint(unsigned long iPoint);
 
   /*!
    * \brief return wether a point is a boundary of a boundary
    */
-  bool GetIsBoundaryPoint(unsigned long iPoint) const override;
+  bool GetIsBoundaryPoint(unsigned long iPoint) const;
+
+  /*!
+   * \brief Allocate member variables for points marked as vertex (via "MarkAsBoundaryPoint").
+   */
+  void AllocateBoundaryVariables();
 };

SU2_CFD/include/variables/CVariable.hpp

@@ -2242,14 +2242,4 @@ class CVariable {
   virtual su2double GetSourceTerm_DispAdjoint(unsigned long iPoint, unsigned long iDim) const { return 0.0; }
   virtual su2double GetSourceTerm_VelAdjoint(unsigned long iPoint, unsigned long iDim) const { return 0.0; }
 
-  /*!
-   * \brief Mark a point as boundary of a design surface boundary
-   */
-  virtual void MarkAsBoundaryPoint(unsigned long iPoint) {}
-
-  /*!
-   * \brief return wether a point is a boundary of a design surface boundary
-   */
-  virtual bool GetIsBoundaryPoint(unsigned long iPoint) const { return false; }
-
 };

SU2_CFD/src/solvers/CFEASolver.cpp

@@ -229,6 +229,8 @@ CFEASolver::~CFEASolver(void) {
     for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++)
       omp_destroy_lock(&UpdateLocks[iPoint]);
   }
+
+  delete nodes;
 }
 
 void CFEASolver::HybridParallelInitialization(CGeometry* geometry) {

SU2_CFD/src/solvers/CFEASolverBase.cpp

@@ -69,7 +69,6 @@ CFEASolverBase::~CFEASolverBase() {
     delete [] element_container;
   }
 
-  delete nodes;
 }
 
 void CFEASolverBase::CommunicateExtraEliminationVertices(const CGeometry* geometry, vector<unsigned long>& myPoints) {