UmatInterface.hpp

Go to the documentation of this file.

/** \file UmatInterface.hpp
 * \brief UmatInterface
 *
 * This file defines the UmatInterface base class, which provides
 * an interface for integrating Umat material models with MoFEM.
 *
 * \ingroup umat
 */
 
#ifndef UMAT_INTERFACE_HPP
#define UMAT_INTERFACE_HPP
 
#include <dlfcn.h>
#include <MatOps.hpp>
 
// Abaqus UMAT function pointer type
// Last arg is the hidden gfortran character length for cmname
using umat_f = void (*)(
    double *stress, double *statev, double *ddsdde, double *sse, double *spd,
    double *scd, double *rpl, double *ddsddt, double *drplde, double *drpldt,
    double *stran, double *dstran, double *time, double *dtime, double *temp,
    double *dtemp, double *predef, double *dpred, char *cmname, int *ndi,
    int *nshr, int *ntens, int *nstatv, double *props, int *nprops,
    double *coords, double *drot, double *pnewdt, double *celent,
    double *dfgrd0, double *dfgrd1, int *noel, int *npt, int *layer, int *kspt,
    int *kstep, int *kinc, std::size_t cmname_len);
 
namespace MoFEM {
 
template <int MODEL_TYPE> struct UmatInterface {
 
  static constexpr int DIM = (MODEL_TYPE == MatOps::MODEL_3D ||
                              MODEL_TYPE == MatOps::MODEL_AXISYMMETRIC) ? 3 : 2;
 
  static constexpr int NTENS = (MODEL_TYPE == MatOps::MODEL_3D) ? 6
                               : (MODEL_TYPE == MatOps::MODEL_2D_PLANE_STRESS)
                                   ? 3
                                   : 4;
 
  static constexpr int NTANGENT = NTENS * NTENS;
 
  UmatInterface(int nstatv = 0, int nprops = 2);
 
  ~UmatInterface();
 
  MoFEMErrorCode initialise(const std::string &library_path = "./umat.so",
                            const std::string &material_name = "UMAT");
 
  MoFEMErrorCode setIdentityDeformationGradient();
  MoFEMErrorCode setIdentityRotationIncrement();
  MoFEMErrorCode setProperties(const std::vector<double> &properties);
  MoFEMErrorCode storeStateVar(int n_state_var,
                               const std::vector<double> &state_var_values);
  MoFEMErrorCode setTimeData(double current_time, double time_increment);
  MoFEMErrorCode setStepData(int step_number, int increment_number);
 
  MoFEMErrorCode setDeformationGradient(const MatrixDouble &grad,
                                        bool use_deformation_gradient);
 
  MoFEMErrorCode setStrainIncrementFromDeformationGradient();
 
  MoFEMErrorCode callUmat();
 
  MoFEMErrorCode getDeformationGradient(MatrixDouble &F) const;
  MoFEMErrorCode getCauchyStressTensor(MatrixDouble &sigma) const;
  MoFEMErrorCode getCauchyTangentTensor(MatrixDouble &cauchy_tangent) const;
 
  boost::shared_ptr<MatOps::MatOpsData> matOpsDataPtr;
 
  umat_f umat = nullptr;
  void *handle = nullptr;
  char filename[256] = {};
 
  static constexpr int NDI = (MODEL_TYPE == MatOps::MODEL_3D) ? 3
                             : (MODEL_TYPE == MatOps::MODEL_2D_PLANE_STRESS)
                                 ? 2
                                 : 3;
  static constexpr int NSHR = (MODEL_TYPE == MatOps::MODEL_3D) ? 3 : 1;
 
  int ndi = NDI, nshr = NSHR, ntens = NTENS, nstatv, nprops;
  int noel = 1, npt = 1, layer = 1, kspt = 1, kstep = 1, kinc = 1;
 
  std::vector<double> props;
  std::vector<double> statev;
 
  // Voigt stress/strain buffers
  std::array<double, NTENS> stress{}; // Cauchy stress at end of step
  std::array<double, NTENS> stran{};  // total strain at start of increment
  std::array<double, NTENS> dstran{}; // strain increment
  // Material tangent (column-major for Fortran: ddsdde[i + NTENS*j] = D(i,j))
  std::array<double, NTANGENT> ddsdde{};
  // Deformation gradients (finite-strain path, column-major DIM×DIM)
  std::array<double, 9> dfgrd0{}; // F at start of increment
  std::array<double, 9> dfgrd1{}; // F at end   of increment
  std::array<double, 9> drot{};
 
  double ddsddt[6];
  double drplde[6];
  double coords[3] = {};
  char cmname[80];
  double time[2] = {};
  double dtime = 0.0;
  double temp = 0.0, dtemp = 0.0;
  std::vector<double> predef{0.0};
  std::vector<double> dpred{0.0};
  double pnewdt = 1.0, celent = 0.0;
  double sse = 0.0, spd = 0.0, scd = 0.0, rpl = 0.0, drpldt = 0.0;
};
 
} // namespace MoFEM
 
#endif // UMAT_INTERFACE_HPP