[WIP] RBF mesh deformation

.gitignore

@@ -103,3 +103,10 @@ su2preconfig.timestamp
 
 # Clangd server files
 .cache
+
+
+
+ninja-win.zip
+ninja.exe
+.gitignore
+Docs/html

Common/include/CConfig.hpp

@@ -189,6 +189,7 @@ class CConfig {
   nMarker_ActDiskBemOutlet_Axis,  /*!< \brief Number of actuator disk BEM outlet markers passed to MARKER_ACTDISK_BEM_AXIS. */
   nMarker_Deform_Mesh_Sym_Plane,  /*!< \brief Number of markers with symmetric deformation */
   nMarker_Deform_Mesh,            /*!< \brief Number of deformable markers at the boundary. */
+  nMarker_Deform_Mesh_Internal,   /*!< \brief Number of internal markers allowed to freely deform. */
   nMarker_Fluid_Load,             /*!< \brief Number of markers in which the flow load is computed/employed. */
   nMarker_Fluid_InterfaceBound,   /*!< \brief Number of fluid interface markers. */
   nMarker_CHTInterface,           /*!< \brief Number of conjugate heat transfer interface markers. */
@@ -238,6 +239,7 @@ class CConfig {
   *Marker_NearFieldBound,         /*!< \brief Near Field boundaries markers. */
   *Marker_Deform_Mesh,            /*!< \brief Deformable markers at the boundary. */
   *Marker_Deform_Mesh_Sym_Plane,  /*!< \brief Marker with symmetric deformation. */
+  *Marker_Deform_Mesh_Internal,   /*!< \brief Internal marker allowed to freely deform. */
   *Marker_Fluid_Load,             /*!< \brief Markers in which the flow load is computed/employed. */
   *Marker_Fluid_InterfaceBound,   /*!< \brief Fluid interface markers. */
   *Marker_CHTInterface,           /*!< \brief Conjugate heat transfer interface markers. */
@@ -652,6 +654,10 @@ class CConfig {
   su2double Deform_Tol_Factor;       /*!< \brief Factor to multiply smallest volume for deform tolerance (0.001 default) */
   su2double Deform_Coeff;            /*!< \brief Deform coeffienct */
   su2double Deform_Limit;            /*!< \brief Deform limit */
+  DEFORM_KIND Deform_Kind;           /*!< \brief Type of mesh deformation */
+  bool RBF_DataReduction;            /*!< \brief Determines use of data reduction methods for RBF mesh deformation. */
+  su2double RBF_GreedyTolerance;      /*!< \brief Tolerance used in the greedy data reduction for RBF mesh deformation. */
+  su2double RBF_GreedyCorrectionFactor;   /*!< \brief Correction factor used in the greedy algorithm for RBF mesh deformation. */
   unsigned short FFD_Continuity;     /*!< \brief Surface continuity at the intersection with the FFD */
   unsigned short FFD_CoordSystem;    /*!< \brief Define the coordinates system */
   su2double Deform_ElasticityMod,    /*!< \brief Young's modulus for volume deformation stiffness model */
@@ -746,6 +752,7 @@ class CConfig {
   *Marker_All_Moving,                /*!< \brief Global index for moving surfaces using the grid information. */
   *Marker_All_Deform_Mesh,           /*!< \brief Global index for deformable markers at the boundary. */
   *Marker_All_Deform_Mesh_Sym_Plane, /*!< \brief Global index for markers with symmetric deformations. */
+  *Marker_All_Deform_Mesh_Internal, /*!< \brief Global index for internal markers with free deformation. */
   *Marker_All_Fluid_Load,            /*!< \brief Global index for markers in which the flow load is computed/employed. */
   *Marker_All_PyCustom,              /*!< \brief Global index for Python customizable surfaces using the grid information. */
   *Marker_All_Designing,             /*!< \brief Global index for moving using the grid information. */
@@ -762,6 +769,7 @@ class CConfig {
   *Marker_CfgFile_Moving,             /*!< \brief Global index for moving surfaces using the config information. */
   *Marker_CfgFile_Deform_Mesh,        /*!< \brief Global index for deformable markers at the boundary. */
   *Marker_CfgFile_Deform_Mesh_Sym_Plane, /*!< \brief Global index for markers with symmetric deformations. */
+  *Marker_CfgFile_Deform_Mesh_Internal, /*!< \brief Global index for internal markers with free deformation. */
   *Marker_CfgFile_Fluid_Load,         /*!< \brief Global index for markers in which the flow load is computed/employed. */
   *Marker_CfgFile_PyCustom,           /*!< \brief Global index for Python customizable surfaces using the config information. */
   *Marker_CfgFile_DV,                 /*!< \brief Global index for design variable markers using the config information. */
@@ -3487,6 +3495,14 @@ class CConfig {
    */
   void SetMarker_All_Deform_Mesh_Sym_Plane(unsigned short val_marker, unsigned short val_deform) { Marker_All_Deform_Mesh_Sym_Plane[val_marker] = val_deform; }
 
+  /*!
+   * \brief Set if a marker <i>val_marker</i> allows deformation at the boundary.
+   * \param[in] val_marker - Index of the marker in which we are interested.
+   * \param[in] val_interface - 0 or 1 depending if the the marker is or not a DEFORM_MESH_SYM_PLANE marker.
+   */
+  void SetMarker_All_Deform_Mesh_Internal(unsigned short val_marker, unsigned short val_deform) { Marker_All_Deform_Mesh_Internal[val_marker] = val_deform; }
+
+
   /*!
    * \brief Set if a in marker <i>val_marker</i> the flow load will be computed/employed.
    * \param[in] val_marker - Index of the marker in which we are interested.
@@ -3637,6 +3653,13 @@ class CConfig {
    * \return 0 or 1 depending if the marker belongs to the DEFORM_MESH_SYM_PLANE subset.
    */
   unsigned short GetMarker_All_Deform_Mesh_Sym_Plane(unsigned short val_marker) const { return Marker_All_Deform_Mesh_Sym_Plane[val_marker]; }
+
+  /*!
+   * \brief Get whether marker <i>val_marker</i> is a DEFORM_MESH_SYM_PLANE marker
+   * \param[in] val_marker - 0 or 1 depending if the the marker belongs to the DEFORM_MESH_SYM_PLANE subset.
+   * \return 0 or 1 depending if the marker belongs to the DEFORM_MESH_SYM_PLANE subset.
+   */
+  unsigned short GetMarker_All_Deform_Mesh_Internal(unsigned short val_marker) const { return Marker_All_Deform_Mesh_Internal[val_marker]; }
 
   /*!
    * \brief Get whether marker <i>val_marker</i> is a Fluid_Load marker
@@ -4375,6 +4398,30 @@ class CConfig {
    */
   bool GetFFD_Symmetry_Plane(void) const { return FFD_Symmetry_Plane; }
 
+  /*!
+   * \brief Get the type of mesh deformation method.
+   * \return type of mesh deformation.
+   */
+  DEFORM_KIND GetDeform_Kind() const { return Deform_Kind; }
+
+  /*!
+   * \brief Determines use of data reduction methods for RBF mesh deformation.
+   * \return <code>TRUE</code> means that data reduction is used.
+   */
+  bool GetRBF_DataReduction(void) const { return RBF_DataReduction; }
+
+  /*!
+   * \brief Determines use of data reduction methods for RBF mesh deformation.
+   * \return <code>TRUE</code> means that data reduction is used.
+   */
+  su2double GetRBF_DataRedTolerance(void) const { return RBF_GreedyTolerance; }
+
+  /*!
+   * \brief Determines use of data reduction methods for RBF mesh deformation.
+   * \return <code>TRUE</code> means that data reduction is used.
+   */
+  su2double GetRBF_DataRedCorrectionFactor(void) const { return RBF_GreedyCorrectionFactor; }
+
   /*!
    * \brief Get the kind of SU2 software component.
    * \return Kind of the SU2 software component.
@@ -6359,6 +6406,12 @@ class CConfig {
    */
   unsigned short GetMarker_CfgFile_Deform_Mesh_Sym_Plane(const string& val_marker) const;
 
+  /*!
+   * \brief Get the DEFORM_MESH_INTERNAL information from the config definition for the marker <i>val_marker</i>.
+   * \return DEFORM_MESH_INTERNAL information of the boundary in the config information for the marker <i>val_marker</i>.
+   */
+  unsigned short GetMarker_CfgFile_Deform_Mesh_Internal(const string& val_marker) const;
+
   /*!
    * \brief Get the Fluid_Load information from the config definition for the marker <i>val_marker</i>.
    * \return Fluid_Load information of the boundary in the config information for the marker <i>val_marker</i>.
@@ -6714,6 +6767,12 @@ class CConfig {
    */
   unsigned short GetMarker_Deform_Mesh_Sym_Plane(const string& val_marker) const;
 
+  /*!
+   * \brief Get the internal index for a DEFORM_MESH_SYM_PLANE boundary <i>val_marker</i>.
+   * \return Internal index for a DEFORM_MESH_SYM_PLANE boundary <i>val_marker</i>.
+   */
+  unsigned short GetMarker_Deform_Mesh_Internal(const string& val_marker) const;
+
   /*!
    * \brief Get a bool for whether the marker is deformed. <i>val_marker</i>.
    * \param[in] val_marker - Name of the marker to test.

Common/include/grid_movement/CLinearElasticity.hpp

@@ -0,0 +1,246 @@
+/*!
+ * \file CLinearElasticity.hpp
+ * \brief Headers of the CLinearElasticity class.
+ * \author F. Palacios, A. Bueno, T. Economon, S. Padron.
+ * \version 8.0.1 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2023, 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 "CVolumetricMovement.hpp"
+#include "../linear_algebra/CSysMatrix.hpp"
+#include "../linear_algebra/CSysVector.hpp"
+#include "../linear_algebra/CSysSolve.hpp"
+
+/*!
+ * \class CLinearElasticity
+ * \brief Class for moving the volumetric numerical grid using the linear elasticity analogy.
+ * \author F. Palacios, A. Bueno, T. Economon, S. Padron.
+ */
+
+class CLinearElasticity final: public CVolumetricMovement{
+ protected: 
+  unsigned short nVar; /*!< \brief Number of variables. */
+
+  unsigned long nPoint;       /*!< \brief Number of points. */
+  unsigned long nPointDomain; /*!< \brief Number of points in the domain. */
+
+  unsigned long nIterMesh; /*!< \brief Number of iterations in the mesh update. +*/
+
+
+ #ifndef CODI_FORWARD_TYPE
+  CSysMatrix<su2mixedfloat> StiffMatrix; /*!< \brief Stiffness matrix of the elasticity problem. */
+  CSysSolve<su2mixedfloat> System;       /*!< \brief Linear solver/smoother. */
+#else
+  CSysMatrix<su2double> StiffMatrix;
+  CSysSolve<su2double> System;
+#endif
+  CSysVector<su2double> LinSysSol;
+  CSysVector<su2double> LinSysRes;
+
+ public:
+  /*!
+  * \brief Constructor of the class.
+  */
+  CLinearElasticity(CGeometry* geometry, CConfig* config);
+
+  /*!
+   * \brief Destructor of the class.
+   */
+  ~CLinearElasticity() override;
+
+  /*!
+   * \brief Grid deformation using the spring analogy method.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] UpdateGeo - Update geometry.
+   * \param[in] Derivative - Compute the derivative (disabled by default). Does not actually deform the grid if enabled.
+   */
+  void SetVolume_Deformation(CGeometry* geometry, CConfig* config, bool UpdateGeo, bool Derivative,
+                                              bool ForwardProjectionDerivative);
+
+  /*!
+   * \brief Update the value of the coordinates after the grid movement.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void UpdateGridCoord(CGeometry* geometry, CConfig* config);
+
+  /*!
+   * \brief Update the derivatives of the coordinates after the grid movement.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void UpdateGridCoord_Derivatives(CGeometry* geometry, CConfig* config, bool ForwardProjectionDerivative);
+
+  /*!
+   * \brief Compute the minimum distance to the nearest solid surface.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void ComputeSolid_Wall_Distance(CGeometry* geometry, CConfig* config, su2double& MinDistance,
+                                su2double& MaxDistance) const;
+
+  /*!
+   * \brief Compute the stiffness matrix for grid deformation using spring analogy.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \return Value of the length of the smallest edge of the grid.
+   */                              
+  su2double SetFEAMethodContributions_Elem(CGeometry* geometry, CConfig* config);
+
+  /*!
+   * \brief Build the stiffness matrix for a 3-D hexahedron element. The result will be placed in StiffMatrix_Elem.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] StiffMatrix_Elem - Element stiffness matrix to be filled.
+   * \param[in] CoordCorners - Index value for Node 1 of the current hexahedron.
+   * \param[in] PointCorners - Index values for element corners
+   * \param[in] nNodes - Number of nodes defining the element.
+   * \param[in] scale
+   */
+  void SetFEA_StiffMatrix2D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem,
+                                            unsigned long PointCorners[8], su2double CoordCorners[8][3],
+                                            unsigned short nNodes, su2double ElemVolume, su2double ElemDistance) ;
+
+  /*!
+   * \brief Build the stiffness matrix for a 3-D hexahedron element. The result will be placed in StiffMatrix_Elem.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] StiffMatrix_Elem - Element stiffness matrix to be filled.
+   * \param[in] CoordCorners - Index value for Node 1 of the current hexahedron.
+   * \param[in] PointCorners - Index values for element corners
+   * \param[in] nNodes - Number of nodes defining the element.
+   * \param[in] scale
+   */
+  void SetFEA_StiffMatrix3D(CGeometry* geometry, CConfig* config, su2double** StiffMatrix_Elem,
+                                              unsigned long PointCorners[8], su2double CoordCorners[8][3],
+                                              unsigned short nNodes, su2double ElemVolume, su2double ElemDistance);
+
+  /*!
+   * \brief Add the stiffness matrix for a 2-D triangular element to the global stiffness matrix for the entire mesh
+   * (node-based). \param[in] geometry - Geometrical definition of the problem. \param[in] StiffMatrix_Elem - Element
+   * stiffness matrix to be filled. \param[in] PointCorners - Index values for element corners \param[in] nNodes -
+   * Number of nodes defining the element.
+   */                                
+  void AddFEA_StiffMatrix(CGeometry* geometry, su2double** StiffMatrix_Elem,
+                                            unsigned long PointCorners[8], unsigned short nNodes);                
+
+    /*!
+   * \brief Shape functions and derivative of the shape functions
+   * \param[in] Xi - Local coordinates.
+   * \param[in] Eta - Local coordinates.
+   * \param[in] Zeta - Local coordinates.
+   * \param[in] CoordCorners - Coordiantes of the corners.
+   * \param[in] DShapeFunction - Shape function information
+   */
+  su2double ShapeFunc_Hexa(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3],
+                          su2double DShapeFunction[8][4]);
+
+  /*!
+   * \brief Shape functions and derivative of the shape functions
+   * \param[in] Xi - Local coordinates.
+   * \param[in] Eta - Local coordinates.
+   * \param[in] Zeta - Local coordinates.
+   * \param[in] CoordCorners - Coordiantes of the corners.
+   * \param[in] DShapeFunction - Shape function information
+   */
+  su2double ShapeFunc_Tetra(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3],
+                            su2double DShapeFunction[8][4]);
+
+  /*!
+   * \brief Shape functions and derivative of the shape functions
+   * \param[in] Xi - Local coordinates.
+   * \param[in] Eta - Local coordinates.
+   * \param[in] Zeta - Local coordinates.
+   * \param[in] CoordCorners - Coordiantes of the corners.
+   * \param[in] DShapeFunction - Shape function information
+   */
+  su2double ShapeFunc_Pyram(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3],
+                            su2double DShapeFunction[8][4]);
+
+  /*!
+   * \brief Shape functions and derivative of the shape functions
+   * \param[in] Xi - Local coordinates.
+   * \param[in] Eta - Local coordinates.
+   * \param[in] Zeta - Local coordinates.
+   * \param[in] CoordCorners - Coordiantes of the corners.
+   * \param[in] DShapeFunction - Shape function information
+   */
+  su2double ShapeFunc_Prism(su2double Xi, su2double Eta, su2double Zeta, su2double CoordCorners[8][3],
+                            su2double DShapeFunction[8][4]);
+
+  /*!
+   * \brief Shape functions and derivative of the shape functions
+   * \param[in] Xi - Local coordinates.
+   * \param[in] Eta - Local coordinates.
+   * \param[in] CoordCorners - Coordiantes of the corners.
+   * \param[in] DShapeFunction - Shape function information
+   */
+  su2double ShapeFunc_Triangle(su2double Xi, su2double Eta, su2double CoordCorners[8][3],
+                              su2double DShapeFunction[8][4]);
+
+  /*!
+   * \brief Shape functions and derivative of the shape functions
+   * \param[in] Xi - Local coordinates.
+   * \param[in] Eta - Local coordinates.
+   * \param[in] CoordCorners - Coordiantes of the corners.
+   * \param[in] DShapeFunction - Shape function information
+   */
+  su2double ShapeFunc_Quadrilateral(su2double Xi, su2double Eta, su2double CoordCorners[8][3],
+                                    su2double DShapeFunction[8][4]);
+
+  /*!
+  * \brief Check the domain points vertex that are going to be moved.
+  * \param[in] geometry - Geometrical definition of the problem.
+  * \param[in] config - Definition of the particular problem.
+  */
+  void SetDomainDisplacements(CGeometry* geometry, CConfig* config);
+
+  /*!
+   * \brief Check the boundary vertex that are going to be moved.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetBoundaryDisplacements(CGeometry* geometry, CConfig* config);
+
+  /*!
+   * \brief Set the derivatives of the boundary nodes.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetBoundaryDerivatives(CGeometry* geometry, CConfig* config, bool ForwardProjectionDerivative);
+
+    /*!
+   * \brief Store the number of iterations when moving the mesh.
+   * \param[in] val_nIterMesh - Number of iterations.
+   */
+  inline void Set_nIterMesh(unsigned long val_nIterMesh) { nIterMesh = val_nIterMesh; }
+
+  /*!
+   * \brief Retrieve the number of iterations when moving the mesh.
+   * \param[out] Number of iterations.
+   */
+  inline unsigned long Get_nIterMesh() const { return nIterMesh; }
+};
+
+