Merge pull request #2776 from su2code/pedro/cleanup_elasticity

Cleanup elasticity numerics

Author: pcarruscag

Common/include/geometry/elements/CElement.hpp

@@ -249,7 +249,8 @@ class CElement {
    * \param[in] nodeB - index of Node b.
    * \param[in] val_Kab - value of the matrix K.
    */
-  inline void Add_Kab(unsigned short nodeA, unsigned short nodeB, su2double** val_Kab) {
+  template <typename Matrix>
+  inline void Add_Kab(unsigned short nodeA, unsigned short nodeB, Matrix& val_Kab) {
     for (unsigned short iDim = 0; iDim < nDim; iDim++)
       for (unsigned short jDim = 0; jDim < nDim; jDim++) Kab[nodeA](nodeB, iDim * nDim + jDim) += val_Kab[iDim][jDim];
   }
@@ -259,7 +260,8 @@ class CElement {
    * transpose) \param[in] nodeA - index of Node a. \param[in] nodeB - index of Node b. \param[in] val_Kab - value of
    * the matrix K.
    */
-  inline void Add_Kab_T(unsigned short nodeA, unsigned short nodeB, su2double** val_Kab) {
+  template <typename Matrix>
+  inline void Add_Kab_T(unsigned short nodeA, unsigned short nodeB, Matrix& val_Kab) {
     for (unsigned short iDim = 0; iDim < nDim; iDim++)
       for (unsigned short jDim = 0; jDim < nDim; jDim++) Kab[nodeA](nodeB, iDim * nDim + jDim) += val_Kab[jDim][iDim];
   }

Common/include/linear_algebra/CPastixWrapper.hpp

@@ -57,10 +57,10 @@ class CPastixWrapper {
   pastix_int_t nCols;            /*!< \brief Local number of columns. */
   vector<pastix_int_t> colptr;   /*!< \brief Equiv. to our "row_ptr". */
   vector<pastix_int_t> rowidx;   /*!< \brief Equiv. to our "col_ind". */
-  vector<passivedouble> values;  /*!< \brief Equiv. to our "matrix". */
+  vector<su2mixedfloat> values;  /*!< \brief Equiv. to our "matrix". */
   vector<pastix_int_t> loc2glb;  /*!< \brief Global index of the columns held by this rank. */
   vector<pastix_int_t> perm;     /*!< \brief Ordering computed by PaStiX. */
-  vector<passivedouble> workvec; /*!< \brief RHS vector which then becomes the solution. */
+  vector<su2mixedfloat> workvec; /*!< \brief RHS vector which then becomes the solution. */
 
   pastix_int_t iparm[IPARM_SIZE];  /*!< \brief Integer parameters for PaStiX. */
   passivedouble dparm[DPARM_SIZE]; /*!< \brief Floating point parameters for PaStiX. */

Common/src/linear_algebra/CPastixWrapper.cpp

@@ -181,7 +181,7 @@ void CPastixWrapper<ScalarType>::Initialize(CGeometry* geometry, const CConfig*
 
   spmInitDist(&spm, SU2_MPI::GetComm());
   spm.mtxtype = SpmGeneral;  // Despite being symmetric, we store the entire matrix.
-  spm.flttype = SpmDouble;
+  spm.flttype = std::is_same_v<su2mixedfloat, double> ? SpmDouble : SpmFloat;
   spm.fmttype = SpmCSC;
   spm.layout = SpmColMajor;
   spm.baseval = 1;

SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp

@@ -27,6 +27,7 @@
 
 #pragma once
 
+#include <memory>
 #include "../CNumerics.hpp"
 #include "../../../../Common/include/geometry/elements/CElement.hpp"
 
@@ -61,23 +62,17 @@ class CFEAElasticity : public CNumerics {
   su2double Kappa     = 0.0;              /*!< \brief Aux. variable, Compressibility constant. */
   su2double ThermalStressTerm = 0.0;      /*!< \brief Aux. variable, Relationship between stress and delta T. */
 
-  su2double *E_i      = nullptr;          /*!< \brief Young's modulus of elasticity. */
-  su2double *Nu_i     = nullptr;          /*!< \brief Poisson's ratio. */
-  su2double *Rho_s_i  = nullptr;          /*!< \brief Structural density. */
-  su2double *Rho_s_DL_i = nullptr;        /*!< \brief Structural density (for dead loads). */
-  su2double *Alpha_i  = nullptr;          /*!< \brief Thermal expansion coefficient. */
+  std::unique_ptr<su2double[]> E_i;        /*!< \brief Young's modulus of elasticity. */
+  std::unique_ptr<su2double[]> Nu_i;       /*!< \brief Poisson's ratio. */
+  std::unique_ptr<su2double[]> Rho_s_i;    /*!< \brief Structural density. */
+  std::unique_ptr<su2double[]> Rho_s_DL_i; /*!< \brief Structural density (for dead loads). */
+  std::unique_ptr<su2double[]> Alpha_i;    /*!< \brief Thermal expansion coefficient. */
 
   su2double ReferenceTemperature = 0.0;   /*!< \brief Reference temperature for thermal expansion. */
 
-  su2double **Ba_Mat = nullptr;           /*!< \brief Matrix B for node a - Auxiliary. */
-  su2double **Bb_Mat = nullptr;           /*!< \brief Matrix B for node b - Auxiliary. */
-  su2double *Ni_Vec  = nullptr;           /*!< \brief Vector of shape functions - Auxiliary. */
-  su2double **D_Mat  = nullptr;           /*!< \brief Constitutive matrix - Auxiliary. */
-  su2double **KAux_ab = nullptr;          /*!< \brief Node ab stiffness matrix - Auxiliary. */
-  su2double **GradNi_Ref_Mat = nullptr;   /*!< \brief Gradients of Ni - Auxiliary. */
-  su2double **GradNi_Curr_Mat = nullptr;  /*!< \brief Gradients of Ni - Auxiliary. */
+  su2double D_Mat[DIM_STRAIN_3D][DIM_STRAIN_3D];  /*!< \brief Constitutive matrix - Auxiliary. */
 
-  su2double *DV_Val = nullptr;            /*!< \brief For optimization cases, value of the design variables. */
+  std::unique_ptr<su2double[]> DV_Val;    /*!< \brief For optimization cases, value of the design variables. */
   unsigned short n_DV = 0;                /*!< \brief For optimization cases, number of design variables. */
 
   bool plane_stress = false;              /*!< \brief Checks if we are solving a plane stress case. */
@@ -97,11 +92,6 @@ class CFEAElasticity : public CNumerics {
    */
   CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CFEAElasticity(void) override;
-
   /*!
    * \brief Set elasticity modulus and Poisson ratio.
    * \param[in] iVal - Index of the property.
@@ -245,4 +235,25 @@ class CFEAElasticity : public CNumerics {
     return static_cast<passivedouble>(iVar == jVar);
   }
 
+  template <typename Mat1, typename Mat2>
+  void FillBMat(unsigned short iNode, const Mat1& GradNi_Mat, Mat2& B_Mat) const {
+    if (nDim == 2) {
+      B_Mat[0][0] = GradNi_Mat[iNode][0];
+      B_Mat[1][1] = GradNi_Mat[iNode][1];
+      B_Mat[2][0] = GradNi_Mat[iNode][1];
+      B_Mat[2][1] = GradNi_Mat[iNode][0];
+    }
+    else {
+      B_Mat[0][0] = GradNi_Mat[iNode][0];
+      B_Mat[1][1] = GradNi_Mat[iNode][1];
+      B_Mat[2][2] = GradNi_Mat[iNode][2];
+      B_Mat[3][0] = GradNi_Mat[iNode][1];
+      B_Mat[3][1] = GradNi_Mat[iNode][0];
+      B_Mat[4][0] = GradNi_Mat[iNode][2];
+      B_Mat[4][2] = GradNi_Mat[iNode][0];
+      B_Mat[5][1] = GradNi_Mat[iNode][2];
+      B_Mat[5][2] = GradNi_Mat[iNode][1];
+    }
+  }
+
 };

SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp

@@ -55,11 +55,6 @@ class CFEALinearElasticity : public CFEAElasticity {
    */
   CFEALinearElasticity(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CFEALinearElasticity(void) override = default;
-
   /*!
    * \brief Build the tangent stiffness matrix of an element.
    * \param[in,out] element_container - Element whose tangent matrix is being built.
@@ -110,11 +105,6 @@ class CFEAMeshElasticity final : public CFEALinearElasticity {
    */
   CFEAMeshElasticity(unsigned short val_nDim, unsigned short val_nVar, unsigned long val_nElem, const CConfig *config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CFEAMeshElasticity(void) override = default;
-
   /*!
    * \brief Set the element-based local Young's modulus in mesh problems
    * \param[in] iElem - Element index.

SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp

@@ -43,39 +43,29 @@ class CFEANonlinearElasticity : public CFEAElasticity {
 
 protected:
 
-  su2double **F_Mat;             /*!< \brief Deformation gradient. */
-  su2double **b_Mat;             /*!< \brief Left Cauchy-Green Tensor. */
-  su2double **currentCoord;      /*!< \brief Current coordinates. */
-  su2double **Stress_Tensor;     /*!< \brief Cauchy stress tensor */
-
-  su2double **FmT_Mat;           /*!< \brief Deformation gradient inverse and transpose. */
-
-  su2double **KAux_P_ab;         /*!< \brief Auxiliar matrix for the pressure term */
-  su2double *KAux_t_a;           /*!< \brief Auxiliar matrix for the pressure term */
+  su2double F_Mat[MAXNDIM][MAXNDIM];   /*!< \brief Deformation gradient. */
+  su2double b_Mat[MAXNDIM][MAXNDIM];   /*!< \brief Left Cauchy-Green Tensor. */
+  su2double Stress_Tensor[MAXNDIM][MAXNDIM];
 
   su2double J_F;                 /*!< \brief Jacobian of the transformation (determinant of F) */
 
   su2double f33;                 /*!< \brief Plane stress term for non-linear 2D plane stress analysis */
 
   bool nearly_incompressible;    /*!< \brief Boolean to consider nearly_incompressible effects */
 
-  su2double **F_Mat_Iso;         /*!< \brief Isocoric component of the deformation gradient. */
-  su2double **b_Mat_Iso;         /*!< \brief Isocoric component of the left Cauchy-Green tensor. */
-
-  su2double C10, D1;             /*!< \brief C10 = Mu/2. D1 = Kappa/2. */
-  su2double J_F_Iso;             /*!< \brief J_F_Iso: det(F)^-1/3. */
+  su2double b_Mat_Iso[MAXNDIM][MAXNDIM]; /*!< \brief Isocoric component of the left Cauchy-Green tensor. */
 
   su2double cijkl[3][3][3][3];   /*!< \brief Constitutive tensor i,j,k,l (defined only for incompressibility - near inc.). */
 
   bool maxwell_stress;           /*!< \brief Consider the effects of the dielectric loads */
 
-  su2double *EField_Ref_Unit,    /*!< \brief Electric Field, unitary, in the reference configuration. */
-  *EField_Ref_Mod;               /*!< \brief Electric Field, modulus, in the reference configuration. */
-  su2double *EField_Curr_Unit;   /*!< \brief Auxiliary vector for the unitary Electric Field in the current configuration. */
+  std::unique_ptr<su2double[]> EField_Ref_Unit;    /*!< \brief Electric Field, unitary, in the reference configuration. */
+  std::unique_ptr<su2double[]> EField_Ref_Mod;     /*!< \brief Electric Field, modulus, in the reference configuration. */
+  std::unique_ptr<su2double[]> EField_Curr_Unit;   /*!< \brief Auxiliary vector for the unitary Electric Field in the current configuration. */
   unsigned short nElectric_Field,
   nDim_Electric_Field;
 
-  su2double *ke_DE_i;            /*!< \brief Electric Constant for Dielectric Elastomers. */
+  std::unique_ptr<su2double[]> ke_DE_i;  /*!< \brief Electric Constant for Dielectric Elastomers. */
 
   su2double ke_DE;               /*!< \brief Electric Constant for Dielectric Elastomers. */
   su2double EFieldMod_Ref;       /*!< \brief Modulus of the electric field in the reference configuration. */
@@ -94,11 +84,6 @@ class CFEANonlinearElasticity : public CFEAElasticity {
    */
   CFEANonlinearElasticity(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CFEANonlinearElasticity(void) override;
-
   /*!
    * \brief Set element electric field.
    * \param[in] i_DV - Index of the variable.
@@ -161,9 +146,48 @@ class CFEANonlinearElasticity : public CFEAElasticity {
   void SetElectric_Properties(const CElement *element_container, const CConfig *config);
 
   /*!
-   * \brief TODO: Describe what this does.
+   * \brief Computes b_Mat.
+   */
+  void ComputeLeftCauchyGreenTensor() {
+    for (unsigned short iVar = 0; iVar < MAXNDIM; iVar++) {
+      for (unsigned short jVar = 0; jVar < MAXNDIM; jVar++) {
+        b_Mat[iVar][jVar] = 0;
+        for (unsigned short kVar = 0; kVar < MAXNDIM; kVar++) {
+          b_Mat[iVar][jVar] += F_Mat[iVar][kVar]*F_Mat[jVar][kVar];
+        }
+      }
+    }
+  }
+
+  /*!
+   * \brief Computes the determinant of the deformation gradient.
+   */
+  void ComputeJ_F() {
+    J_F = F_Mat[0][0]*F_Mat[1][1]*F_Mat[2][2]+
+          F_Mat[0][1]*F_Mat[1][2]*F_Mat[2][0]+
+          F_Mat[0][2]*F_Mat[1][0]*F_Mat[2][1]-
+          F_Mat[0][2]*F_Mat[1][1]*F_Mat[2][0]-
+          F_Mat[1][2]*F_Mat[2][1]*F_Mat[0][0]-
+          F_Mat[2][2]*F_Mat[0][1]*F_Mat[1][0];
+  }
+
+  /*!
+   * \brief Computes the deformation gradient transpose inverse.
    */
-  void Compute_FmT_Mat(void);
+  template <typename Mat>
+  void Compute_FmT_Mat(const Mat& F_Mat, const su2double& J_F, Mat& FmT_Mat) const {
+    FmT_Mat[0][0] = (F_Mat[1][1]*F_Mat[2][2] - F_Mat[1][2]*F_Mat[2][1]) / J_F;
+    FmT_Mat[0][1] = (F_Mat[1][2]*F_Mat[2][0] - F_Mat[2][2]*F_Mat[1][0]) / J_F;
+    FmT_Mat[0][2] = (F_Mat[1][0]*F_Mat[2][1] - F_Mat[1][1]*F_Mat[2][0]) / J_F;
+
+    FmT_Mat[1][0] = (F_Mat[0][2]*F_Mat[2][1] - F_Mat[0][1]*F_Mat[2][2]) / J_F;
+    FmT_Mat[1][1] = (F_Mat[0][0]*F_Mat[2][2] - F_Mat[2][0]*F_Mat[0][2]) / J_F;
+    FmT_Mat[1][2] = (F_Mat[0][1]*F_Mat[2][1] - F_Mat[0][0]*F_Mat[2][0]) / J_F;
+
+    FmT_Mat[2][0] = (F_Mat[0][1]*F_Mat[1][2] - F_Mat[0][2]*F_Mat[1][1]) / J_F;
+    FmT_Mat[2][1] = (F_Mat[0][2]*F_Mat[1][0] - F_Mat[0][0]*F_Mat[1][2]) / J_F;
+    FmT_Mat[2][2] = (F_Mat[0][0]*F_Mat[1][1] - F_Mat[0][1]*F_Mat[1][0]) / J_F;
+  }
 
   /*!
    * \brief TODO: Describe what this does.

SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp

@@ -48,11 +48,6 @@ class CFEM_NeoHookean_Comp final : public CFEANonlinearElasticity {
    */
   CFEM_NeoHookean_Comp(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CFEM_NeoHookean_Comp(void) override = default;
-
 private:
   /*!
    * \brief Compute the plane stress term.
@@ -100,11 +95,6 @@ class CFEM_Knowles_NearInc final : public CFEANonlinearElasticity {
    */
   CFEM_Knowles_NearInc(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CFEM_Knowles_NearInc(void) override = default;
-
 private:
   /*!
    * \brief Compute the plane stress term.
@@ -150,11 +140,6 @@ class CFEM_DielectricElastomer final : public CFEANonlinearElasticity {
    */
   CFEM_DielectricElastomer(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CFEM_DielectricElastomer(void) override = default;
-
 private:
   /*!
    * \brief Compute the plane stress term.
@@ -202,11 +187,6 @@ class CFEM_IdealDE final : public CFEANonlinearElasticity {
    */
   CFEM_IdealDE(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CFEM_IdealDE(void) override = default;
-
 private:
   /*!
    * \brief Compute the plane stress term.

SU2_CFD/include/output/filewriter/CTecplotBinaryFileWriter.hpp

@@ -27,9 +27,11 @@
 
 #pragma once
 
+#include <assert.h>
+#include <cstdint>
+
 #include "CFileWriter.hpp"
 
-#include <assert.h>
 
 class CTecplotBinaryFileWriter final: public CFileWriter{
 

SU2_CFD/src/iteration/CFEAIteration.cpp

@@ -123,11 +123,13 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom
         config[val_iZone]->SetInnerIter(CurIter);
         break;
       }
-      /*--- Linear elasticity without thermal effects only needs one iteration. ---*/
+      /*--- Linear elasticity without thermal effects and double precision only needs one iteration. ---*/
+#ifndef USE_MIXED_PRECISION
       if (linear && !heat) {
         output->SetConvergence(true);
         break;
       }
+#endif
       /*--- Normal stopping criteria. ---*/
       if (StopCalc && IntIter > 0) break;
     }