UmatInterface.cpp

Go to the documentation of this file.

/** \file UmatInterface.cpp
 * \brief UmatInterface
 *
 * UmatInterface
 *
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include "UmatInterface.hpp"
 
namespace UmatOps {
 
inline int tensorIndex(const int ii, const int jj) { return 3 * ii + jj; }
 
inline MoFEMErrorCode resizeTensor2(MatrixDouble &mat) {
  MoFEMFunctionBeginHot;
  mat.resize(3, 3, false);
  mat.clear();
  MoFEMFunctionReturnHot(0);
}
 
inline MoFEMErrorCode resizeTensor4(MatrixDouble &mat) {
  MoFEMFunctionBeginHot;
  mat.resize(9, 9, false);
  mat.clear();
  MoFEMFunctionReturnHot(0);
}
 
}; // namespace UmatOps
// MatUmat can do the conversions depending of the strain measure used.
 
using namespace UmatOps;
 
template <int MODEL_TYPE>
MoFEM::UmatInterface<MODEL_TYPE>::UmatInterface(int nstatv, int nprops)
    : nstatv(nstatv), nprops(nprops), props(nprops, 0.0), statev(nstatv, 0.0) {}
 
template <int MODEL_TYPE> MoFEM::UmatInterface<MODEL_TYPE>::~UmatInterface() {
  if (handle) {
    dlclose(handle);
    handle = nullptr;
    umat = nullptr;
  }
}
 
template <int MODEL_TYPE>
MoFEMErrorCode
MoFEM::UmatInterface<MODEL_TYPE>::initialise(const std::string &library_path,
                                             const std::string &material_name) {
  MoFEMFunctionBegin;
 
  if (handle) {
    dlclose(handle);
    handle = nullptr;
    umat = nullptr;
  }
 
  std::memset(filename, 0, sizeof(filename));
  std::strncpy(filename, library_path.c_str(), sizeof(filename) - 1);
 
  std::memset(cmname, ' ', sizeof(cmname));
  const auto material_name_size = material_name.size() < sizeof(cmname)
                                      ? material_name.size()
                                      : sizeof(cmname);
  std::memcpy(cmname, material_name.data(), material_name_size);
 
  handle = dlopen(filename, RTLD_NOW | RTLD_LOCAL);
  if (!handle) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Failed to open UMAT library %s: %s", filename, dlerror());
  }
 
  umat = (umat_f)dlsym(handle, "umat_");
  if (!umat) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Failed to resolve symbol umat_ in %s: %s", filename, dlerror());
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode MoFEM::UmatInterface<MODEL_TYPE>::callUmat() {
  MoFEMFunctionBegin;
 
  if (!umat) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "UMAT function pointer is not initialised");
  }
 
  umat(stress.data(), statev.data(), ddsdde.data(), &sse, &spd, &scd, &rpl,
       ddsddt, drplde, &drpldt, stran.data(), dstran.data(), time, &dtime,
       &temp, &dtemp, predef.data(), dpred.data(), cmname, &ndi, &nshr, &ntens,
       &nstatv, props.data(), &nprops, coords, drot.data(), &pnewdt, &celent,
       dfgrd0.data(), dfgrd1.data(), &noel, &npt, &layer, &kspt, &kstep, &kinc,
       sizeof(cmname));
 
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode
MoFEM::UmatInterface<MODEL_TYPE>::setIdentityDeformationGradient() {
  MoFEMFunctionBegin;
 
  dfgrd0.fill(0);
  for (int dd = 0; dd != 3; ++dd)
    dfgrd0[dd + 3 * dd] = 1.0;
 
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode
MoFEM::UmatInterface<MODEL_TYPE>::setIdentityRotationIncrement() {
  MoFEMFunctionBegin;
 
  drot.fill(0);
  for (int dd = 0; dd != 3; ++dd)
    drot[dd + 3 * dd] = 1.0;
 
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode MoFEM::UmatInterface<MODEL_TYPE>::setDeformationGradient(
    const MatrixDouble &grad, bool use_deformation_gradient) {
  MoFEMFunctionBegin;
 
  dfgrd1.fill(0);
  for (int dd = 0; dd != 3; ++dd)
    dfgrd1[dd + 3 * dd] = 1.0;
 
  for (int row = 0; row != DIM; ++row) {
    for (int col = 0; col != DIM; ++col) {
      const double grad_value = grad(row, col);
      dfgrd1[row + 3 * col] = use_deformation_gradient
                                  ? grad_value
                                  : grad_value + (row == col ? 1.0 : 0.0);
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode
MoFEM::UmatInterface<MODEL_TYPE>::setTimeData(double current_time,
                                              double time_increment) {
  MoFEMFunctionBegin;
 
  time[1] = current_time;
  dtime = time_increment;
  time[0] = current_time - time_increment;
 
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode
MoFEM::UmatInterface<MODEL_TYPE>::setStepData(int step_number,
                                              int increment_number) {
  MoFEMFunctionBegin;
 
  kstep = step_number;
  kinc = increment_number;
 
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode
MoFEM::UmatInterface<MODEL_TYPE>::setStrainIncrementFromDeformationGradient() {
  MoFEMFunctionBegin;
 
  auto get_small_strain_voigt = [](const std::array<double, 9> &dfgrd,
                                   auto &strain) {
    std::array<double, 9> grad = dfgrd;
    for (int dd = 0; dd != DIM; ++dd)
      grad[dd + 3 * dd] -= 1.0;
 
    if constexpr (MODEL_TYPE == MatOps::MODEL_3D) {
      strain[0] = grad[0];
      strain[1] = grad[4];
      strain[2] = grad[8];
      strain[3] = grad[1] + grad[3];
      strain[4] = grad[2] + grad[6];
      strain[5] = grad[5] + grad[7];
    } else if constexpr (MODEL_TYPE == MatOps::MODEL_2D_PLANE_STRAIN) {
      strain[0] = grad[0];
      strain[1] = grad[4];
      strain[2] = 0.0;
      strain[3] = grad[1] + grad[3];
    } else if constexpr (MODEL_TYPE == MatOps::MODEL_2D_PLANE_STRESS) {
      strain[0] = grad[0];
      strain[1] = grad[4];
      strain[2] = grad[1] + grad[3];
    } else if constexpr (MODEL_TYPE == MatOps::MODEL_AXISYMMETRIC) {
      strain[0] = grad[0];
      strain[1] = grad[4];
      strain[2] = grad[8];
      strain[3] = grad[1] + grad[3];
    } else {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "UMAT strain increment is not implemented for this model type");
    }
  };
 
  std::array<double, NTENS> reference_strain{};
  std::array<double, NTENS> current_strain{};
  get_small_strain_voigt(dfgrd0, reference_strain);
  get_small_strain_voigt(dfgrd1, current_strain);
 
  for (int ii = 0; ii != NTENS; ++ii) {
    stran[ii] = reference_strain[ii];
    dstran[ii] = current_strain[ii] - reference_strain[ii];
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode MoFEM::UmatInterface<MODEL_TYPE>::storeStateVar(
    int n_statev, const std::vector<double> &state_var_values) {
  MoFEMFunctionBegin;
 
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode MoFEM::UmatInterface<MODEL_TYPE>::getDeformationGradient(
    MatrixDouble &F) const {
  MoFEMFunctionBegin;
  CHKERR resizeTensor2(F);
  for (int row = 0; row != DIM; ++row)
    for (int col = 0; col != DIM; ++col)
      F(row, col) = dfgrd1[row + 3 * col];
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode MoFEM::UmatInterface<MODEL_TYPE>::getCauchyStressTensor(
    MatrixDouble &sigma) const {
  MoFEMFunctionBegin;
  CHKERR resizeTensor2(sigma);
  if constexpr (MODEL_TYPE == MatOps::MODEL_3D) {
    sigma(0, 0) = stress[0];
    sigma(1, 1) = stress[1];
    sigma(2, 2) = stress[2];
    sigma(0, 1) = sigma(1, 0) = stress[3];
    sigma(0, 2) = sigma(2, 0) = stress[4];
    sigma(1, 2) = sigma(2, 1) = stress[5];
  } else if constexpr (MODEL_TYPE == MatOps::MODEL_2D_PLANE_STRAIN ||
                       MODEL_TYPE == MatOps::MODEL_AXISYMMETRIC) {
    sigma(0, 0) = stress[0];
    sigma(1, 1) = stress[1];
    sigma(2, 2) = stress[2];
    sigma(0, 1) = sigma(1, 0) = stress[3];
  } else if constexpr (MODEL_TYPE == MatOps::MODEL_2D_PLANE_STRESS) {
    sigma(0, 0) = stress[0];
    sigma(1, 1) = stress[1];
    sigma(0, 1) = sigma(1, 0) = stress[2];
  } else {
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "UMAT stress conversion is not implemented for this model type");
  }
  MoFEMFunctionReturn(0);
}
 
template <int MODEL_TYPE>
MoFEMErrorCode MoFEM::UmatInterface<MODEL_TYPE>::getCauchyTangentTensor(
    MatrixDouble &cauchy_tangent) const {
  MoFEMFunctionBegin;
  CHKERR resizeTensor4(cauchy_tangent);
  if constexpr (MODEL_TYPE == MatOps::MODEL_3D) {
    static const std::array<std::pair<int, int>, 6> pairs = {
        std::pair<int, int>{0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}};
    for (int row = 0; row != 6; ++row) {
      const auto [i, j] = pairs[row];
      for (int col = 0; col != 6; ++col) {
        const auto [k, l] = pairs[col];
        const double value = ddsdde[row + 6 * col];
        cauchy_tangent(tensorIndex(i, j), tensorIndex(k, l)) = value;
        cauchy_tangent(tensorIndex(j, i), tensorIndex(k, l)) = value;
        cauchy_tangent(tensorIndex(i, j), tensorIndex(l, k)) = value;
        cauchy_tangent(tensorIndex(j, i), tensorIndex(l, k)) = value;
      }
    }
  } else if constexpr (MODEL_TYPE == MatOps::MODEL_2D_PLANE_STRAIN ||
                       MODEL_TYPE == MatOps::MODEL_AXISYMMETRIC) {
    static const std::array<std::pair<int, int>, 4> pairs = {
        std::pair<int, int>{0, 0}, {1, 1}, {2, 2}, {0, 1}};
    for (int row = 0; row != 4; ++row) {
      const auto [i, j] = pairs[row];
      for (int col = 0; col != 4; ++col) {
        const auto [k, l] = pairs[col];
        const double value = ddsdde[row + 4 * col];
        cauchy_tangent(tensorIndex(i, j), tensorIndex(k, l)) = value;
        cauchy_tangent(tensorIndex(j, i), tensorIndex(k, l)) = value;
        cauchy_tangent(tensorIndex(i, j), tensorIndex(l, k)) = value;
        cauchy_tangent(tensorIndex(j, i), tensorIndex(l, k)) = value;
      }
    }
  } else if constexpr (MODEL_TYPE == MatOps::MODEL_2D_PLANE_STRESS) {
    static const std::array<std::pair<int, int>, 3> pairs = {
        std::pair<int, int>{0, 0}, {1, 1}, {0, 1}};
    for (int row = 0; row != 3; ++row) {
      const auto [i, j] = pairs[row];
      for (int col = 0; col != 3; ++col) {
        const auto [k, l] = pairs[col];
        const double value = ddsdde[row + 3 * col];
        cauchy_tangent(tensorIndex(i, j), tensorIndex(k, l)) = value;
        cauchy_tangent(tensorIndex(j, i), tensorIndex(k, l)) = value;
        cauchy_tangent(tensorIndex(i, j), tensorIndex(l, k)) = value;
        cauchy_tangent(tensorIndex(j, i), tensorIndex(l, k)) = value;
      }
    }
  } else {
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "UMAT tangent conversion is not implemented for this model type");
  }
  MoFEMFunctionReturn(0);
}
 
template class MoFEM::UmatInterface<MatOps::MODEL_3D>;
template class MoFEM::UmatInterface<MatOps::MODEL_2D_PLANE_STRAIN>;
template class MoFEM::UmatInterface<MatOps::MODEL_2D_PLANE_STRESS>;
template class MoFEM::UmatInterface<MatOps::MODEL_AXISYMMETRIC>;