b2material_solid_mechanics_umat.H

//------------------------------------------------------------------------
// b2material_solid_mechanics_umat.H --
//
//     UMAT material api wrapper.
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//
// (c) 2009,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.
//------------------------------------------------------------------------
 
// TODO: set DROT, apply the rot to the strain, stress and ddsdde
// Thus this wrapper is only valid for small displacement.
 
#include "elements/properties/b2fortran_element_property.H"
#include "elements/properties/b2solid_mechanics.H"
#include "io/b2000_db/b2000_db_v3/b2fortran_element_property.H"
#include "ip/b2ip_constants.H"
#include "model/b2solution.H"
 
#ifndef B2_MATERIAL_SOLID_MECHANICS_UMAT_H_
#define B2_MATERIAL_SOLID_MECHANICS_UMAT_H_
 
namespace b2000 { namespace element {
 
typedef void(umat_subroutine_f)(
      // Output in all situations (at the and of the increment).
      double* STRESS,  // (NTENS), stress at the end of the increment
      double* STATEV,  // (NSTATV), solution-dependent state variables
      double* DDSDDE,  // (NTENS,NTENS), Jacobian matrix of the constitutive model
      double* SSE,     // elastic strain energy
      double* SPD,     // plastic dissipation
      double* SCD,     // “creep” dissipation
      // Output, only in a fully coupled temperature-displacement analysis
      double* RPL,  // volumetric heat generation per unit time
      double*
            DDSDDT,  // (NTENS), variation of the stress increments with respect to the temperature
      double* DRPLDE,  // (NTENS), variation of RPL with respect to the strain increments
      double* DRPLDT,  // variation of RPL with respect to the temperature
      // Input (at the beginning of the increment)
      const double* STRAN,   // (NTENS), total strains at the beginning of the increment
      const double* DSTRAN,  // (NTENS), strains increment
      const double TIME[2],  // step time and total time
      const double* DTIME,   // time increment
      const double* TEMP,    // temperature
      const double* DTEMP,   // temperature increment
      const double* PREDEF,  // array of interpolated values of predefined field variables at this
                             // point at the start of the increment
      const double* DPRED,   // array of increments of predefined field variables
      const char* CMNAME,    // Name given on the *MATERIAL option
      const int* NDI,        // number of direct stress components at this point
      const int* NSHR,       // number of engineering shear stress components at this point
      const int* NTENS,      // size of the stress or strain component array
      const int* NSTATV,    // number of solution-dependent state variables that are associated with
                            // this material type
      const double* PROPS,  // array of material constants
      const int* NPROPS,    // number of material constants
      const double COORDS[3],     // array containing the coordinates of this point
      const double DROT[3][3],    // rotation increment matrix
      double* PNEWDT,             // ratio of suggested new time increment
      const double* CELENT,       // characteristic element length
      const double DFGRD0[3][3],  // deformation gradient at the beginning of the increment
      const double DFGRD1[3][3],  // deformation gradient at the end of the increment
      const int* NOEL,            // element number
      const int* NPT,             // integration point number
      const int* LAYER,           // layer number
      const int* KSPT,            // section point number
      const int* KSTEP,           // step number
      const int* KINC,            // increment number
      const int CNAME_len);
 
template <umat_subroutine_f UMAT_SUBROUTINE, const char MATERIAL_TYPE_NAME[], int NSTATV>
class MaterialSolidMechanics3DUMAT_copy;
 
template <umat_subroutine_f UMAT_SUBROUTINE, const char MATERIAL_TYPE_NAME[], int NSTATV>
class MaterialSolidMechanics3DUMAT : virtual public b2000::b2dbv3::FortranElementProperty,
                                     public b2000::SolidMechanics3D {
public:
    LinearType linear(const int layer_id = -1) const { return nonlinear; }
 
    bool path_dependent(const int layer_id = -1) const { return true; }
 
    SolidMechanics3D* copy(int layer_) {
        return new MaterialSolidMechanics3DUMAT_copy<UMAT_SUBROUTINE, MATERIAL_TYPE_NAME, NSTATV>(
              *this);
    }
 
    bool isotropic(const int layer_id = -1) const { return false; }
 
    void get_dynamic_nonlinear_stress(
          Model* model, const Element* element,
          b2000::b2linalg::Vector<double, b2000::b2linalg::Vdense_constref> nodes_interpolation,
          const double element_coordinate[3], const int layer_id,
          const double covariant_natural_base[3][3], const double displacement_gradient[3][3],
          const double green_lagrange_strain[6], const double green_lagrange_strain_dot[6],
          const double velocity[3], const double acceleration[3], const double time,
          const double delta_time, const double area,
          const EquilibriumSolution equilibrium_solution, double piola_kirchhoff_stress[6],
          double d_piola_kirchhoff_stress_d_strain[21],
          double d_piola_kirchhoff_stress_d_strain_dot[21],
          double d_piola_kirchhoff_stress_d_time[6], double inertial_force[3], double& density,
          GradientContainer* gradient_container, SolverHints* solver_hints,
          const double covariant_base[3][3], const double contravariant_base[3][3]) {
        get_static_nonlinear_stress(
              model, element, nodes_interpolation, element_coordinate, layer_id,
              covariant_natural_base, displacement_gradient, green_lagrange_strain, time,
              delta_time, area, equilibrium_solution, piola_kirchhoff_stress,
              d_piola_kirchhoff_stress_d_strain, d_piola_kirchhoff_stress_d_time,
              gradient_container, solver_hints, covariant_base, contravariant_base);
    }
 
    void get_static_nonlinear_stress(
          Model* model, const Element* element, Vector<double, Vdense_constref> nodes_interpolation,
          const double element_coordinate[3], const int layer_id,
          const double covariant_natural_base[3][3], const double displacement_gradient[3][3],
          const double green_lagrange_strain[6], const double time, const double delta_time,
          const double area, const EquilibriumSolution equilibrium_solution,
          double piola_kirchhoff_stress[6], double d_piola_kirchhoff_stress_d_strain[21],
          double d_piola_kirchhoff_stress_d_time[6], GradientContainer* gradient_container,
          SolverHints* solver_hints, const double covariant_base[3][3] = 0,
          const double contravariant_base[3][3] = 0) {
        Exception() << "A copy of the object must be used." << THROW;
    }
 
    typedef ObjectTypeComplete<
          MaterialSolidMechanics3DUMAT<UMAT_SUBROUTINE, MATERIAL_TYPE_NAME, NSTATV>,
          FortranElementProperty::type_t>
          type_t;
    static type_t type;
};
 
template <umat_subroutine_f UMAT_SUBROUTINE, const char MATERIAL_TYPE_NAME[], int NSTATV>
typename MaterialSolidMechanics3DUMAT<UMAT_SUBROUTINE, MATERIAL_TYPE_NAME, NSTATV>::type_t
      MaterialSolidMechanics3DUMAT<UMAT_SUBROUTINE, MATERIAL_TYPE_NAME, NSTATV>::type(
            MATERIAL_TYPE_NAME, "", StringList(), element::module, &b2000::ElementProperty::type);
 
template <umat_subroutine_f UMAT_SUBROUTINE, const char MATERIAL_TYPE_NAME[], int NSTATV>
class MaterialSolidMechanics3DUMAT_copy : public b2000::SolidMechanics3D {
public:
    typedef MaterialSolidMechanics3DUMAT<UMAT_SUBROUTINE, MATERIAL_TYPE_NAME, NSTATV> parent_type;
 
    MaterialSolidMechanics3DUMAT_copy(parent_type& parent_)
        : parent(parent_), conv_time(0), conv_sse(0), conv_spd(0), conv_scd(0) {
        std::fill(conv_deformation_gradient[0], conv_deformation_gradient[3], 0);
        std::fill(conv_strain_ea, conv_strain_ea + 6, 0);
        std::fill(conv_stress_ea, conv_stress_ea + 6, 0);
        std::fill(conv_statv, conv_statv + NSTATV, 0);
    }
 
    LinearType linear(const int layer_id = -1) const { return nonlinear; }
 
    bool path_dependent(const int layer_id = -1) const { return true; }
 
    SolidMechanics3D* copy(int layer_) { return new MaterialSolidMechanics3DUMAT_copy(*this); }
 
    bool isotropic(const int layer_id = -1) const { return false; }
 
    void get_dynamic_nonlinear_stress(
          Model* model, const Element* element, Vector<double, Vdense_constref> nodes_interpolation,
          const double element_coordinate[3], const int layer_id,
          const double covariant_natural_base[3][3], const double displacement_gradient[3][3],
          const double green_lagrange_strain[6], const double green_lagrange_strain_dot[6],
          const double velocity[3], const double acceleration[3], const double time,
          const double delta_time, const double area,
          const EquilibriumSolution equilibrium_solution, double piola_kirchhoff_stress[6],
          double d_piola_kirchhoff_stress_d_strain[21],
          double d_piola_kirchhoff_stress_d_strain_dot[21],
          double d_piola_kirchhoff_stress_d_time[6], double inertial_force[3], double& density,
          GradientContainer* gradient_container, SolverHints* solver_hints,
          const double covariant_base[3][3], const double contravariant_base[3][3]) {
        get_static_nonlinear_stress(
              model, element, nodes_interpolation, element_coordinate, layer_id,
              covariant_natural_base, displacement_gradient, green_lagrange_strain, time,
              delta_time, area, equilibrium_solution, piola_kirchhoff_stress,
              d_piola_kirchhoff_stress_d_strain, d_piola_kirchhoff_stress_d_time,
              gradient_container, solver_hints, covariant_base, contravariant_base);
    }
 
    void get_static_nonlinear_stress(
          Model* model, const Element* element, Vector<double, Vdense_constref> nodes_interpolation,
          const double element_coordinate[3], const int layer_id,
          const double covariant_natural_base[3][3], const double displacement_gradient[3][3],
          const double green_lagrange_strain[6], const double time, const double delta_time,
          const double area, const EquilibriumSolution equilibrium_solution,
          double piola_kirchhoff_stress[6], double d_piola_kirchhoff_stress_d_strain[21],
          double d_piola_kirchhoff_stress_d_time[6], GradientContainer* gradient_container,
          SolverHints* solver_hints, const double covariant_base[3][3] = 0,
          const double contravariant_base[3][3] = 0) {
        // Obtain info on the material constant and orientation of the material.
        bool material_ref_identity;
        double material_ref[3][3];
        const double* props = parent.get_material_and_material_referential(
              element, layer_id, covariant_natural_base, material_ref, material_ref_identity, true);
        const int nprops = int(props[B2MATPOSSTART]);
 
        // Compute the Linear/Green-Lagrange strain in the branch-global
        // orthogonal system.
        double strain_bg[6];
        if (green_lagrange_strain == 0) {
            std::fill_n(strain_bg, 6, 0.0);
        } else if (covariant_base == 0) {
            std::copy(green_lagrange_strain, green_lagrange_strain + 6, strain_bg);
        } else {
            transform_strain_from_base_A_to_I(covariant_base, green_lagrange_strain, strain_bg)
                :
 
                  // Compute the Linear/Green-Lagrange strain in the element-aligned
                  // orthogonal system.
                  double strain_ea[6];
        }
        if (material_ref_identity) {
            std::copy(strain_bg, strain_bg + 6, strain_ea);
        } else {
            transform_strain_from_I_to_base_B(material_ref, strain_bg, strain_ea);
        }
 
        // Compute the strain increment in the element-aligned
        // orthogonal system.
        double dstrain[6];
        for (int i = 0; i != 6; ++i) { dstrain[i] = strain_ea[i] - conv_strain_ea[i]; }
 
        // Get the integration point coordinates.
        double coords[3] = {};
        std::pair<int, Node* const*> nodes = element->get_nodes();
        for (int i = 0; i != nodes.first; ++i) {
            const double* c = nodes.second[i]->get_coor().second;
            for (int j = 0; j != 3; ++j) { coords[j] += nodes_interpolation[i] * c[j]; }
        }
 
        // Compute the deformation gradient
        double deformation_gradient[3][3];
        for (int i = 0; i != 3; ++i) {
            for (int j = 0; j != 3; ++j) {
                deformation_gradient[j][i] = displacement_gradient[j][i];
            }
            deformation_gradient[i][i] += 1;
        }
 
        double dtime = time - conv_time;
 
        double pnewdt = delta_time;
 
        double statv[NSTATV];
        std::copy(conv_statv, conv_statv + NSTATV, statv);
 
        double sse = conv_sse;
        double spd = conv_spd;
        double scd = conv_scd;
 
        double stress_ea[6];
        std::copy(conv_stress_ea, conv_stress_ea + 6, stress_ea);
 
        double drot[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
        const double celent = area;
 
        const int noel = element->get_id();
        const int ntp = 0;
        const int layer = layer_id + 1;
 
        const int val3 = 3;
        const int val6 = 6;
        const int nstatv = NSTATV;
        const int val_unknown = 0;
        const double time_all[2] = {conv_time, conv_time};
 
        double ddsdde[6][6];
 
        // Call the Fortran umat subroutine.
        UMAT_SUBROUTINE(
              stress_ea, statv, ddsdde[0], &sse, &spd, &scd, 0, 0, 0,
              0,  // no coupled temperature-displacement analysis
              conv_strain_ea, dstrain, time_all, &dtime, 0,
              0,     // no coupled temperature-displacement analysis
              0, 0,  // no array of interpolated values of predefined field
              MATERIAL_TYPE_NAME, &val3, &val3, &val6, &nstatv, props + B2MATPOSSTART + 1, &nprops,
              coords, drot, &pnewdt, &celent, conv_deformation_gradient, deformation_gradient,
              &noel, &val_unknown, &layer, &val_unknown, &val_unknown, &val_unknown,
              strlen(MATERIAL_TYPE_NAME));
 
        if (equilibrium_solution) {
            std::copy(
                  deformation_gradient[0], deformation_gradient[3], conv_deformation_gradient[0]);
            std::copy(strain_ea, strain_ea + 6, conv_strain_ea);
            std::copy(stress_ea, stress_ea + 6, conv_stress_ea);
            conv_time = time;
            std::copy(statv, statv + NSTATV, conv_statv);
        }
 
        double G[3][3];
        if (covariant_base != 0) {
            copy_3_3(covariant_base, G);
        } else {
            set_identity_3_3(G);
        }
        double G_d[3][3];
        get_base_opposite_variance(G, G_d);
 
        //
        // Transform the Jacobian matrix of the constitutive model
        // either to the branch-global orthogonal system or to
        // system required by the element and store it in the
        // 'd_stress_d_strain' array.
        //
        if (d_piola_kirchhoff_stress_d_strain != NULL) {
            double d_stress_d_strain_ea[21];
            // Convert the ddsdde rectangular 6x6 matrix to the
            // d_stress_d_strain_ea lower packed matrix.
            // The symmetric part of the matrix is calculated by taking one
            // half the sum of the ddsdde matrix and its transpose.
            {
                double* ptr = d_stress_d_strain_ea;
                for (int j = 0; j != 6; ++j) {
                    for (int i = j; i != 6; ++i, ++ptr) {
                        *ptr = 0.5 * (ddsdde[j][i] + ddsdde[i][j]);
                    }
                }
            }
 
            // Transform the constitutive tensor from the element
            // covariant basis to the material reference frame.
            transform_constitutive_tensor_from_base_A_to_base_B(
                  G, material_ref, d_stress_d_strain_ea, d_piola_kirchhoff_stress_d_strain);
        }
 
        // Compute the Linear/2nd-Piola-Kirchhoff stress in the
        // branch-global orthogonal system.
        double stress_bg[6];
        transform_stress_from_base_A_to_I(material_ref, stress_ea, stress_bg);
 
        // Store the Linear/2nd-Piola-Kirchhoff stress in the 'stress'
        // array, either in the branch-global orthogonal system or in the
        // contravariant base given by the element.
        if (piola_kirchhoff_stress != 0) {
            transform_stress_from_I_to_base_B(G_d, stress_bg, piola_kirchhoff_stress);
        }
 
        if (gradient_container == 0) { return; }
 
        // Remove the x2 in the shear components of the
        // Linear/Green-Lagrange strain in the branch-global orthogonal
        // system.
        double strain_bg_storage[6];
        std::copy(strain_bg, strain_bg + 6, strain_bg_storage);
        for (int i = 3; i != 6; ++i) { strain_bg_storage[i] *= 0.5; }
 
        // Transform the Linear/2nd-Piola-Kirchhoff stress to the
        // Linear/Cauchy stress in branch-global orthogonal system.
        double stress_bg_storage[6];
        if (!nonlinear) {
            std::copy(stress_bg, stress_bg + 6, stress_bg_storage);
        } else {
            // Transform the 2nd-Piola-Kirchhoff stress in the
            // element-aligned orthogonal system to the Cauchy stress in
            // the branch-global orthogonal system.
 
            // 1. Compute the deformation gradient from the displacement
            // gradient.
            double deformation_gradient[3][3];
            if (displacement_gradient != NULL) {
                if (contravariant_base) {
                    inner_product_3_3_NN(
                          displacement_gradient, contravariant_base, deformation_gradient);
                } else {
                    std::copy(
                          displacement_gradient[0], displacement_gradient[3],
                          deformation_gradient[0]);
                }
            } else {
                std::fill(deformation_gradient[0], deformation_gradient[3], 0);
            }
            deformation_gradient[0][0] += 1;
            deformation_gradient[1][1] += 1;
            deformation_gradient[2][2] += 1;
 
            // 2. Compute the inverse of the determinant of the
            // deformation gradient.
            const double det_inv = 1.0 / determinant_3_3(deformation_gradient);
 
            // 3. Compute X*S*X^T.
            if (material_ref_identity) {
                transform_stress_from_base_A_ov_to_I(
                      deformation_gradient, stress_ea, stress_bg_storage);
            } else {
                double deformation_gradient_aligned_ortho[3][3];
                inner_product_3_3_NN(
                      deformation_gradient, material_ref, deformation_gradient_aligned_ortho);
                transform_stress_from_base_A_ov_to_I(
                      deformation_gradient_aligned_ortho, stress_ea, stress_bg_storage);
            }
 
            // 4. Divide by the determinant of the deformation gradient.
            for (int i = 0; i != 6; ++i) { stress_bg_storage[i] *= det_inv; }
        }
 
        //
        // Storage in the database.
        //
        static const GradientContainer::FieldDescription strain_descr = {
              "STRAIN",
              "Exx.Eyy.Ezz.Exy.Eyz.Exz",
              (nonlinear ? "Green Lagrange strain" : "Linear strain"),
              3,
              6,
              2,
              3,
              true,
              typeid(double)};
 
        static const GradientContainer::FieldDescription stress_descr = {
              "STRESS",
              "Sxx.Syy.Szz.Sxy.Syz.Sxz",
              (nonlinear ? "Cauchy stress" : "Linear stress"),
              3,
              6,
              2,
              3,
              true,
              typeid(double)};
 
        // Fill strain and stress into the GradientContainer object.
        if (gradient_container->compute_field_value(strain_descr.name)) {
            gradient_container->set_field_value(strain_descr, strain_bg_storage);
        }
        if (gradient_container->compute_field_value(stress_descr.name)) {
            gradient_container->set_field_value(stress_descr, stress_bg_storage);
        }
    }
 
private:
    parent_type& parent;
    double conv_deformation_gradient[3][3];
    double conv_strain_ea[6];
    double conv_stress_ea[6];
    double conv_time;
    double conv_statv[NSTATV];
    double conv_sse;
    double conv_spd;
    double conv_scd;
};
 
}}  // namespace b2000::element
 
#define DEFINE_NEW_MATERIAL_SOLID_MECHANICS_3D_UMAT(                                          \
      subroutine_name, SUBROUTINE_NAME, material_name, NSTATV)                                \
    extern "C" b2000::element::umat_subroutine_f FC_GLOBAL(subroutine_name, SUBROUTINE_NAME); \
    namespace {                                                                               \
    char mname[] = material_name;                                                             \
    }                                                                                         \
    template class b2000::element::MaterialSolidMechanics3DUMAT<                              \
          FC_GLOBAL(subroutine_name, SUBROUTINE_NAME), mname, NSTATV>;
 
#endif  // B2_MATERIAL_SOLID_MECHANICS_UMAT_H_