[WIP] Coupled solver for thermoelasticity

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -28,6 +28,7 @@
 #pragma once
 
 #include "CFEASolverBase.hpp"
+#include "CHeatSolver.hpp"
 
 /*!
  * \class CFEASolver
@@ -100,6 +101,22 @@ class CFEASolver : public CFEASolverBase {
   bool initial_calc = true;    /*!< \brief Becomes false after first call to Preprocessing. */
   bool body_forces = false;    /*!< \brief Whether any body force is active. */
 
+  /*!
+ * \brief Pointer to the heat solver for coupled simulations.
+ *
+ * This member stores a pointer to the heat solver, which handles
+ * the solution of the heat equation in weakly coupled simulations.
+ * It is initialized during the preprocessing step if the configuration
+ * enables the weak coupling of heat and elasticity solvers. This solver
+ * provides temperature information to the finite element elasticity solver
+ * and contributes to the coupled residuals.
+ *
+ * The `heat_solver` pointer remains `nullptr` when the heat solver is not enabled
+ * in the configuration. Memory for the heat solver is dynamically allocated
+ * during initialization and released in the destructor to avoid memory leaks.
+ */
+  CSolver* heat_solver = nullptr;
+
   /*!
    * \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.

SU2_CFD/src/iteration/CFEAIteration.cpp

@@ -48,6 +48,13 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom
   CIntegration* feaIntegration = integration[val_iZone][val_iInst][FEA_SOL];
   CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL];
 
+  /*--- Add heat solver integration step ---*/
+    if (config[val_iZone]->GetWeakly_Coupled_Heat()) {
+    config[val_iZone]->SetGlobalParam(MAIN_SOLVER::HEAT_EQUATION, RUNTIME_HEAT_SYS);
+    integration[val_iZone][val_iInst][HEAT_SOL]->SingleGrid_Iteration(
+        geometry, solver, numerics, config, RUNTIME_HEAT_SYS, val_iZone, val_iInst);
+    }
+
   /*--- FEA equations ---*/
   config[val_iZone]->SetGlobalParam(MAIN_SOLVER::FEM_ELASTICITY, RUNTIME_FEA_SYS);
 

SU2_CFD/src/solvers/CFEASolver.cpp

@@ -30,6 +30,7 @@
 #include "../../include/numerics/elasticity/CFEAElasticity.hpp"
 #include "../../../Common/include/toolboxes/printing_toolbox.hpp"
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
+#include "../../include/solvers/CHeatSolver.hpp"
 #include <algorithm>
 
 using namespace GeometryToolbox;
@@ -563,6 +564,10 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container,
   const bool disc_adj_fem = (config->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_FEM);
   const bool topology_mode = config->GetTopology_Optimization();
 
+  if (config->GetWeakly_Coupled_Heat()) {
+        heat_solver = solver_container[HEAT_SOL];
+    }
+
   /*
    * For topology optimization we apply a filter on the design density field to avoid
    * numerical issues (checkerboards), ensure mesh independence, and impose a length scale.
@@ -687,6 +692,10 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics,
             su2double val_Sol = nodes->GetSolution(indexNode[iNode],iDim) + val_Coord;
             element->SetRef_Coord(iNode, iDim, val_Coord);
             element->SetCurr_Coord(iNode, iDim, val_Sol);
+          }  
+          if (heat_solver) {
+                        auto nodal_temperatures = dynamic_cast<CSolver *>(heat_solver)->GetNodalTemperature();
+                        element->SetTemperature(iNode, nodal_temperatures[indexNode[iNode]]);
           }
         }
 

SU2_CFD/src/solvers/CSolverFactory.cpp

@@ -173,6 +173,9 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon
       break;
     case MAIN_SOLVER::FEM_ELASTICITY:
       solver[FEA_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel);
+      if (config->GetWeakly_Coupled_Heat()) {
+        solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel);
+      }
       break;
     case MAIN_SOLVER::DISC_ADJ_FEM:
       solver[FEA_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::FEA, solver, geometry, config, iMGLevel);