CFEAElasticity.hpp

Go to the documentation of this file.

/*!
 * \file CFEAElasticity.hpp
 * \brief Declaration and inlines of the base class for elasticity problems.
 * \author Ruben Sanchez
 * \version 8.5.0 "Harrier"
 *
 * SU2 Project Website: https://su2code.github.io
 *
 * The SU2 Project is maintained by the SU2 Foundation
 * (http://su2foundation.org)
 *
 * Copyright 2012-2026, 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 <memory>
#include "../CNumerics.hpp"
#include "../../../../Common/include/geometry/elements/CElement.hpp"
 
/*!
 * \class CFEAElasticity
 * \ingroup Elasticity_Equations
 * \brief Abstract class for computing the tangent matrix and the residual for structural problems.
 * \note  At the next level of abstraction (linear or not) a class must define the constitutive term.
 *        The methods we override in this class with an empty implementation are here just to better
 *        document the public interface of this class hierarchy.
 * \author R.Sanchez
 * \version 8.5.0 "Harrier"
 */
class CFEAElasticity : public CNumerics {
 
protected:
 
  enum : unsigned short {DIM_STRAIN_2D = 3,   /*!< \brief Exx, Eyy, Gxy. */
                         DIM_STRAIN_3D = 6};  /*!< \brief Exx, Eyy, Ezz, Gxy, Gxz, Gyz. */
 
  enum : unsigned short {NNODES_2D = 4,   /*!< \brief Maximum number of nodes for 2D problems. */
                         NNODES_3D = 8};  /*!< \brief Maximum number of nodes for 3D problems. */
 
  su2double E         = 1.0;              /*!< \brief Aux. variable, Young's modulus of elasticity. */
  su2double Nu        = 0.0;              /*!< \brief Aux. variable, Poisson's ratio. */
  su2double Rho_s     = 0.0;              /*!< \brief Aux. variable, Structural density. */
  su2double Rho_s_DL  = 0.0;              /*!< \brief Aux. variable, Structural density (for dead loads). */
  su2double Alpha     = 0.0;              /*!< \brief Aux. variable, thermal expansion coefficient. */
 
  su2double Mu        = 0.0;              /*!< \brief Aux. variable, Lame's coeficient. */
  su2double Lambda    = 0.0;              /*!< \brief Aux. variable, Lame's coeficient. */
  su2double Kappa     = 0.0;              /*!< \brief Aux. variable, Compressibility constant. */
  su2double ThermalStressTerm = 0.0;      /*!< \brief Aux. variable, Relationship between stress and delta T. */
 
  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 D_Mat[DIM_STRAIN_3D][DIM_STRAIN_3D];  /*!< \brief Constitutive matrix - Auxiliary. */
 
  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. */
  bool linear = false;                    /*!< \brief Checks if we are solving a linear elasticity case. */
 
public:
  /*!
   * \brief Default constructor
   */
  CFEAElasticity() = default;
 
  /*!
   * \brief Constructor of the class (overload).
   * \param[in] val_nDim - Number of dimensions of the problem.
   * \param[in] val_nVar - Number of variables of the problem.
   * \param[in] config - Definition of the particular problem.
   */
  CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config);
 
  /*!
   * \brief Set elasticity modulus and Poisson ratio.
   * \param[in] iVal - Index of the property.
   * \param[in] val_E - Value of the elasticity (Young) modulus.
   * \param[in] val_Nu - Value of the Poisson ratio.
   */
  inline void SetMaterial_Properties(unsigned short iVal, su2double val_E, su2double val_Nu) final {
    E_i[iVal] = val_E;
    Nu_i[iVal] = val_Nu;
  }
 
  /*!
   * \brief Set densities, real and for gravity loading purposes.
   * \param[in] iVal - Index of the property.
   * \param[in] val_Rho - Material density.
   * \param[in] val_Rho_DL - Density for gravity (dead) loads.
   */
  inline void SetMaterial_Density(unsigned short iVal, su2double val_Rho, su2double val_Rho_DL) final {
    Rho_s_i[iVal] = val_Rho;
    Rho_s_DL_i[iVal] = val_Rho_DL;
  }
 
  /*!
   * \brief Set element electric field.
   * \param[in] i_DV - Index of the variable.
   * \param[in] val_EField - Value of the field.
   */
  inline void Set_ElectricField(unsigned short i_DV, su2double val_EField) override { }
 
  /*!
   * \brief Set the element-based local Young's modulus in mesh problems
   * \param[in] iElem - Element index.
   * \param[in] val_E - Value of elasticity modulus.
   */
  inline void SetMeshElasticProperties(unsigned long iElem, su2double val_E) override { }
 
  /*!
   * \brief Set the value of a design variable.
   * \param[in] i_DV - Index of the variable.
   * \param[in] val_DV - Value of the variable.
   */
  inline void Set_DV_Val(unsigned short i_DV, su2double val_DV) final { DV_Val[i_DV] = val_DV; }
 
  /*!
   * \brief Get the value of a design variable.
   * \param[in] i_DV - Index of the variable.
   * \return Value of the variable.
   */
  inline su2double Get_DV_Val(unsigned short i_DV) const final { return DV_Val[i_DV]; }
 
  /*!
   * \brief Build the mass matrix of an element.
   * \param[in,out] element_container - Element whose mass matrix is being built.
   * \param[in] config - Definition of the problem.
   */
  void Compute_Mass_Matrix(CElement *element_container, const CConfig *config) final;
 
  /*!
   * \brief Compute the nodal inertial loads for an element.
   * \param[in,out] element_container - The element for which the inertial loads are computed.
   * \param[in] config - Definition of the problem.
   */
  void Compute_Body_Forces(CElement *element_container, const CConfig *config) final;
 
  /*!
   * \brief Build the tangent stiffness matrix of an element.
   * \param[in,out] element_container - Element whose tangent matrix is being built.
   * \param[in] config - Definition of the problem.
   */
  inline void Compute_Tangent_Matrix(CElement *element_container, const CConfig *config) override { };
 
  /*!
   * \brief Compute averaged nodal stresses (for post processing).
   * \param[in,out] element_container - The finite element.
   * \param[in] config - Definition of the problem.
   */
  inline su2double Compute_Averaged_NodalStress(CElement *element_container, const CConfig *config) override { return 0; };
 
  /*!
   * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz.
   * \note Szz is required in 2D whereas Sxz and Syz are assumed to be 0.
   */
  template <class T>
  static su2double VonMisesStress(unsigned short nDim, const T& stress) {
    /*--- In 2D, we only have 4 components: Sxx, Syy, Sxy, Szz. ---*/
    const auto& Sxx = stress[0];
    const auto& Syy = stress[1];
    const auto& Sxy = stress[2];
    const auto& Szz = stress[3];
    su2double Sxz = 0, Syz = 0;
    if (nDim == 3) {
      Sxz = stress[4];
      Syz = stress[5];
    }
    return sqrt(0.5 * (pow(Sxx - Syy, 2) + pow(Syy - Szz, 2) + pow(Szz - Sxx, 2) + 6 * (Sxy * Sxy + Sxz * Sxz + Syz * Syz)));
  }
 
protected:
  /*!
   * \brief Compute the constitutive matrix, must be implemented by derived classes.
   * \param[in,out] element_container - The finite element.
   * \param[in] config - Definition of the problem.
   */
  virtual void Compute_Constitutive_Matrix(CElement *element_container, const CConfig *config) = 0;
 
  /*!
   * \brief Set element material properties.
   * \param[in] element_container - Element defining the properties.
   * \param[in] config - Definition of the problem.
   */
  virtual void SetElement_Properties(const CElement *element_container, const CConfig *config);
 
  /*!
   * \brief Read design variables from file.
   * \param[in] config - Definition of the problem.
   */
  void ReadDV(const CConfig *config);
 
  /*!
   * \brief Update the Lame parameters (required in AD to account for all dependencies).
   */
  inline void Compute_Lame_Parameters(void) {
    Mu     = E / (2.0*(1.0 + Nu));
    Lambda = Nu*E/((1.0+Nu)*(1.0-2.0*Nu));
    Kappa  = Lambda + (2/3)*Mu;
    /*--- https://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php
     * The stress tensor for nonlinear problems is still 3x3 and plane stress is imposed
     * by determining the deformation that makes sigma_33 = 0, hence this denominator is
     * only changed for linear elasticity. ---*/
    const auto nu_mult = (linear && plane_stress) ? 1 : 2;
    ThermalStressTerm = -Alpha * E / (1 - nu_mult * Nu);
  }
 
  /*!
   * \brief Kronecker delta.
   * \param[in] iVar - Index i.
   * \param[in] jVar - Index j.
   * \return 1 if i=j, 0 otherwise.
   */
  inline static passivedouble deltaij(unsigned short iVar, unsigned short jVar) {
    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];
    }
  }
 
};