b2material_stress_3d_mixed.H

//------------------------------------------------------------------------
// b2material_stress_3d_mixed --
//
//     Base class for materials used by 3D displacement/pressure
//     elements.
//
// written by Thomas Ludwig
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//
// (c) 2009-2015,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_MIXED_H_
#define B2MATERIAL_STRESS_3D_MIXED_H_
 
#include "elements/properties/b2solid_mechanics.H"
#include "utils/b2tensor_calculus.H"
 
namespace b2000 {
 
class MaterialStress3DMixed : virtual public SolidMechanics {
public:
    // Return the number of layers.
    virtual int number_of_layer() const { return 1; }
 
    // Return the laminate geometric position relative to the
    // mid-surface of the element.
    virtual void laminate(
          b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation,
          double& laminate_begin, double& laminate_end) const {
        laminate_begin = -1;
        laminate_end = 1;
    }
 
    // Return the layer geometric position relative to the
    // mid-surface of the element.
    virtual void layer(
          b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation,
          const int layer_id, double& layer_begin, double& layer_end) const {
        layer_begin = -1;
        layer_end = 1;
    }
 
    virtual double get_density(const int layer_id) const {
        UnimplementedError() << THROW;
        return 0.0;
    }
 
    virtual 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,
          b2linalg::Vector<double, b2linalg::Vdense_constref> nodes_interpolation,
          const double bg_coordinates[3], const double covariant_base[3][3], const double volume,
          const double deformation_gradient[3][3], const double velocity[3],
          const double acceleration[3], const double interpolated_pressure, double stress[6],
          double& pressure, double CUU[21], double CUP[6], double& CPP, double inertia_force[3],
          double& density) = 0;
 
protected:
    inline void eval_linear_strain(const double F[3][3], double E[6]) const {
        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];
 
        E[0] = f00 - 1.0;
        E[1] = f11 - 1.0;
        E[2] = f22 - 1.0;
        E[3] = f01 + f10;
        E[4] = f12 + f21;
        E[5] = f02 + f20;
    }
 
    inline void eval_green_lagrange_strain(const double F[3][3], double E[6]) const {
        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];
 
        E[0] = 0.5 * (f00 * f00 + f10 * f10 + f20 * f20 - 1.0);
        E[1] = 0.5 * (f01 * f01 + f11 * f11 + f21 * f21 - 1.0);
        E[2] = 0.5 * (f02 * f02 + f12 * f12 + f22 * f22 - 1.0);
        E[3] = f00 * f01 + f10 * f11 + f20 * f21;
        E[4] = f01 * f02 + f11 * f12 + f21 * f22;
        E[5] = f00 * f02 + f10 * f12 + f20 * f22;
    }
 
    inline double delta(const int a, const int b) { return (a == b ? 1.0 : 0.0); }
 
    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] + f01 * f01 * S[1] + f02 * f02 * S[2] + 2 * f00 * f01 * S[3]
               + 2 * f01 * f02 * S[4] + 2 * f02 * f00 * S[5];
 
        T[1] = f10 * f10 * S[0] + f11 * f11 * S[1] + f12 * f12 * S[2] + 2 * f10 * f11 * S[3]
               + 2 * f11 * f12 * S[4] + 2 * f12 * f10 * S[5];
 
        T[2] = f20 * f20 * S[0] + f21 * f21 * S[1] + f22 * f22 * S[2] + 2 * f20 * f21 * S[3]
               + 2 * f21 * f22 * S[4] + 2 * f22 * f20 * S[5];
 
        T[3] = f00 * f10 * S[0] + f01 * f11 * S[1] + f02 * f12 * S[2]
               + (f00 * f11 + f01 * f10) * S[3] + (f01 * f12 + f02 * f11) * S[4]
               + (f02 * f10 + f00 * f12) * S[5];
 
        T[4] = f10 * f20 * S[0] + f11 * f21 * S[1] + f12 * f22 * S[2]
               + (f10 * f21 + f11 * f20) * S[3] + (f11 * f22 + f12 * f21) * S[4]
               + (f12 * f20 + f10 * f22) * S[5];
 
        T[5] = f20 * f00 * S[0] + f21 * f01 * S[1] + f22 * f02 * S[2]
               + (f20 * f01 + f21 * f00) * S[3] + (f21 * f02 + f22 * f01) * S[4]
               + (f22 * f00 + f20 * f02) * S[5];
 
        for (int i = 0; i < 6; ++i) { T[i] *= det_inv; }
    }
};
 
}  // namespace b2000
 
#endif  // B2_MATERIAL_STRESS_3D_MIXED_H_