b2material_stress_3d_mesh_deform.H

//------------------------------------------------------------------------
// b2material_stress_3d_mesh_deform.H --
//
//     Material used for 3D Mesh Deformation (Aero-Elasticity).
//
// written by Thomas Ludwig
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//
// (c) 2009-2012,2016 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
//          Linder Höhe, 51147 Köln
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef B2MATERIAL_STRESS_3D_MESH_DEFORM_H_
#define B2MATERIAL_STRESS_3D_MESH_DEFORM_H_
 
#include "elements/properties/b2solid_mechanics.H"
#include "utils/b2tensor_calculus.H"
 
namespace b2000 {
 
class MaterialStress3DMeshDeform : public SolidMechanics3D<double> {
public:
    void get_stress(
          Model* model, const bool linear, const EquilibriumSolution equilibrium_solution,
          const double time, const double delta_time, GradientContainer* gradient_container,
          SolverHints* solver_hints, const Element* element, const double el_coordinates[3],
          const int layer_id,
          const b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation,
          const double bg_coordinates[3], const double covariant_base[3][3],
          const bool use_covariant_base, const double volume,
          const double displacement_gradient[3][3], const double strain[6],
          const double strain_dot[6], const double velocity[3], const double acceleration[3],
          double stress[6], double d_stress_d_strain[21], double d_stress_d_strain_dot[21],
          double d_stress_d_time[6], double inertia_force[3], double& density) override;
 
    inline void s_eq_c_mul_e(double s[6], const double c[21], const double e[6]) {
        const int ind[6][6] = {{0, 1, 2, 3, 4, 5},     {1, 6, 7, 8, 9, 10},
                               {2, 7, 11, 12, 13, 14}, {3, 8, 12, 15, 16, 17},
                               {4, 9, 13, 16, 18, 19}, {5, 19, 14, 17, 19, 20}};
 
        for (int i = 0; i < 6; ++i) {
            double value = 0.0;
 
            const int* iind = ind[i];
            for (int j = 0; j < 6; ++j) { value += c[iind[j]] * e[j]; }
            s[i] = value;
        }
    }
 
    // Transform the 2nd Piola-Kirchhoff stress to the Cauchy
    // stress.  The stresses are in the Voight notation [sxx, syy,
    // szz, sxy, syz, sxz].
    static inline void piola_2_stress_to_cauchy_stress(
          const double F[3][3],  // Deformation gradient, column-major.
          const double S[6],     // 2nd Piola-Kirchhoff stress (Voight not.).
          double T[6]            // Cauchy stress (Voight notation).
    ) {
        const double f00 = F[0][0];
        const double f01 = F[1][0];
        const double f02 = F[2][0];
        const double f10 = F[0][1];
        const double f11 = F[1][1];
        const double f12 = F[2][1];
        const double f20 = F[0][2];
        const double f21 = F[1][2];
        const double f22 = F[2][2];
 
        // Get the determinant of the deformation gradient tensor..
        const double det = determinant_3_3(F);
        const double det_inv = 1.0 / det;
 
        T[0] = f00 * f00 * S[0] + f10 * f10 * S[1] + f20 * f20 * S[2] + 2 * f00 * f10 * S[3]
               + 2 * f10 * f20 * S[4] + 2 * f20 * f00 * S[5];
 
        T[1] = f01 * f01 * S[0] + f11 * f11 * S[1] + f21 * f21 * S[2] + 2 * f01 * f11 * S[3]
               + 2 * f11 * f21 * S[4] + 2 * f21 * f01 * S[5];
 
        T[2] = f02 * f02 * S[0] + f12 * f12 * S[1] + f22 * f22 * S[2] + 2 * f02 * f12 * S[3]
               + 2 * f12 * f22 * S[4] + 2 * f22 * f02 * S[5];
 
        T[3] = f00 * f01 * S[0] + f10 * f11 * S[1] + f20 * f21 * S[2]
               + (f00 * f11 + f10 * f01) * S[3] + (f10 * f21 + f20 * f11) * S[4]
               + (f20 * f01 + f00 * f21) * S[5];
 
        T[4] = f01 * f02 * S[0] + f11 * f12 * S[1] + f21 * f22 * S[2]
               + (f01 * f12 + f11 * f02) * S[3] + (f11 * f22 + f21 * f12) * S[4]
               + (f21 * f02 + f01 * f22) * S[5];
 
        T[5] = f02 * f00 * S[0] + f12 * f10 * S[1] + f22 * f20 * S[2]
               + (f02 * f10 + f12 * f00) * S[3] + (f12 * f20 + f22 * f10) * S[4]
               + (f22 * f00 + f02 * f20) * S[5];
 
        for (int i = 0; i < 6; ++i) { T[i] *= det_inv; }
    }
 
public:
    double mat_e;
    double mat_p;
    double mat_density;
};
 
}  // namespace b2000
 
#endif  // B2_MATERIAL_STRESS_3D_MESH_DEFORM_H_