Integrate CMeshSolver for mesh deformation based on linear elasticity…

… solver

Common/include/geometry/primal_grid/CPrimalGrid.hpp

@@ -32,6 +32,7 @@
 
 #include <iostream>
 #include <vector>
+#include <limits>
 #include <cstdlib>
 
 #include "../../CConfig.hpp"

SU2_CFD/include/solvers/CMeshSolver.hpp

@@ -129,19 +129,28 @@ class CMeshSolver final : public CFEASolver {
   /*!
    * \brief Grid deformation using the linear elasticity equations.
    * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] numerics - Numerics used in the solution.
    * \param[in] config - Definition of the particular problem.
    */
   void DeformMesh(CGeometry **geometry,
                   CNumerics **numerics,
                   CConfig *config) override;
 
   /*!
-   * \brief Set the stiffness of the mesh.
+   * \brief Grid deformation using the linear elasticity equations (no multigrid).
    * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] numerics - Numerics used in the solution.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void DeformMesh(CGeometry *geometry,
+                  CNumerics **numerics,
+                  CConfig *config) override;
+  /*!
+   * \brief Set the stiffness of the mesh.
+   * \param[in] numerics - Numerics used in the solution.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetMesh_Stiffness(CGeometry **geometry,
-                         CNumerics **numerics,
+  void SetMesh_Stiffness(CNumerics **numerics,
                          CConfig *config) override;
 
   /*!

SU2_CFD/include/solvers/CSolver.hpp

@@ -4205,8 +4205,16 @@ class CSolver {
    * \param[in] config - Definition of the particular problem.
    * \param[in] referenceCoord - Determine if the mesh is deformed from the reference or from the current coordinates.
    */
-  inline virtual void SetMesh_Stiffness(CGeometry **geometry,
-                                        CNumerics **numerics,
+  inline virtual void DeformMesh(CGeometry *geometry,
+                                 CNumerics **numerics,
+                                 CConfig *config) { }
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] referenceCoord - Determine if the mesh is deformed from the reference or from the current coordinates.
+   */
+  inline virtual void SetMesh_Stiffness(CNumerics **numerics,
                                         CConfig *config) { }
 
   /*!

SU2_CFD/src/iteration/CIteration.cpp

@@ -163,7 +163,7 @@ void CIteration::SetMesh_Deformation(CGeometry** geometry, CSolver** solver, CNu
 
   /*--- Set the stiffness of each element mesh into the mesh numerics ---*/
 
-  solver[MESH_SOL]->SetMesh_Stiffness(geometry, numerics[MESH_SOL], config);
+  solver[MESH_SOL]->SetMesh_Stiffness(numerics[MESH_SOL], config);
 
   /*--- Deform the volume grid around the new boundary locations ---*/
 

SU2_CFD/src/solvers/CMeshSolver.cpp

@@ -420,7 +420,7 @@ void CMeshSolver::SetWallDistance(CGeometry *geometry, CConfig *config) {
   END_SU2_OMP_PARALLEL
 }
 
-void CMeshSolver::SetMesh_Stiffness(CGeometry **geometry, CNumerics **numerics, CConfig *config){
+void CMeshSolver::SetMesh_Stiffness(CNumerics **numerics, CConfig *config){
 
   if (stiffness_set) return;
 
@@ -542,6 +542,62 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig
 
 }
 
+void CMeshSolver::DeformMesh(CGeometry *geometry, CNumerics **numerics, CConfig *config){
+
+  if (multizone) nodes->Set_BGSSolution_k();
+
+  /*--- Capture a few MPI dependencies for AD. ---*/
+  geometry->InitiateComms(geometry, config, COORDINATES);
+  geometry->CompleteComms(geometry, config, COORDINATES);
+
+  InitiateComms(geometry, config, SOLUTION);
+  CompleteComms(geometry, config, SOLUTION);
+
+  InitiateComms(geometry, config, MESH_DISPLACEMENTS);
+  CompleteComms(geometry, config, MESH_DISPLACEMENTS);
+
+  /*--- Compute the stiffness matrix, no point recording because we clear the residual. ---*/
+
+  const bool wasActive = AD::BeginPassive();
+
+  Compute_StiffMatrix(geometry, numerics, config);
+
+  AD::EndPassive(wasActive);
+
+  /*--- Clear residual (loses AD info), we do not want an incremental solution. ---*/
+  SU2_OMP_PARALLEL {
+    LinSysRes.SetValZero();
+    if (time_domain && config->GetFSI_Simulation()) LinSysSol.SetValZero();
+  }
+  END_SU2_OMP_PARALLEL
+
+  /*--- Impose boundary conditions (all of them are ESSENTIAL BC's - displacements). ---*/
+  SetBoundaryDisplacements(geometry, config, false);
+
+  /*--- Solve the linear system. ---*/
+  Solve_System(geometry, config);
+
+  SU2_OMP_PARALLEL {
+
+  /*--- Update the grid coordinates and cell volumes using the solution
+     of the linear system (usol contains the x, y, z displacements). ---*/
+  UpdateGridCoord(geometry, config);
+
+  /*--- Update the dual grid. ---*/
+  UpdateDualGrid(geometry, config);
+
+  /*--- Check for failed deformation (negative volumes). ---*/
+  SetMinMaxVolume(geometry, config, true);
+
+  /*--- The Grid Velocity is only computed if the problem is time domain ---*/
+  if (time_domain && !config->GetFSI_Simulation())
+    ComputeGridVelocity(geometry, config);
+
+  }
+  END_SU2_OMP_PARALLEL
+
+}
+
 void CMeshSolver::UpdateGridCoord(CGeometry *geometry, CConfig *config){
 
   /*--- Update the grid coordinates using the solution of the linear system ---*/

SU2_DEF/include/drivers/CDeformationDriver.hpp

@@ -1,4 +1,4 @@
-/*!
+ /*!
  * \file CDeformationDriver.hpp
  * \brief Headers of the main subroutines for driving the mesh deformation.
  * \author T. Economon, H. Kline, R. Sanchez
@@ -32,6 +32,7 @@
 #include "../../../Common/include/grid_movement/CSurfaceMovement.hpp"
 #include "../../../Common/include/grid_movement/CVolumetricMovement.hpp"
 #include "../../../SU2_CFD/include/output/COutput.hpp"
+#include "../../../SU2_CFD/include/numerics/CNumerics.hpp"
 #include "../../../Common/include/geometry/CGeometry.hpp"
 
 /*!
@@ -56,6 +57,8 @@ class CDeformationDriver {
   CGeometry **geometry_container;             /*!< \brief Geometrical definition of the problem. */
   CSurfaceMovement **surface_movement;          /*!< \brief Surface movement classes of the problem. */
   CVolumetricMovement **grid_movement;         /*!< \brief Volume grid movement classes of the problem. */
+  CSolver **solver_container;
+  CNumerics ***numerics_container;
   COutput **output_container;                   /*!< \brief Pointer to the COutput class. */
 
 public:
@@ -76,6 +79,11 @@ class CDeformationDriver {
    */
   void Run();
 
+  /*!
+   * \brief Output the mesh.
+   */
+  void Output();
+
   /*!
    * \brief Deallocation routine
    */
@@ -102,4 +110,24 @@ class CDeformationDriver {
    */
   void Output_Preprocessing();
 
+  /*!
+   * \brief Preprocess the mesh solver container.
+   */
+  void Solver_Preprocessing();
+
+  /*!
+   * \brief Preprocess the numerics container.
+   */
+  void Numerics_Preprocessing();
+
+  /*!
+   * \brief Mesh deformation based on linear elasticity solver (CMeshSolver).
+   */
+  void Update();
+
+  /*!
+   * \brief Mesh deformation based on legacy implementation.
+   */
+  void Update_Legacy();
+
 };