MatUmatImpl.hpp

Go to the documentation of this file.

/**
 * @file MatUmatImpl.hpp
 * @brief Shared UMAT implementation details for built-in and user-defined UMATs.
 */
 
#ifndef MAT_UMAT_IMPL_HPP
#define MAT_UMAT_IMPL_HPP
 
#include <MoFEM.hpp>
 
#include "MatOps.hpp"
#include "MatUmat.hpp"
#include <sstream>
#include <UmatInterface.hpp>
 
namespace MatOps {
 
using namespace MoFEM;
 
inline int umatTensorIndex(const int ii, const int jj) { return 3 * ii + jj; }
 
enum StrainType { SMALL, HENCKY, FINITE };
 
template <int DIM, int MODEL_TYPE>
struct SmallStrainHandler : public StrainHandler<DIM, MODEL_TYPE> {
 
  SmallStrainHandler(){};
 
  MoFEMErrorCode runUmat(MatUmatImpl<DIM, MODEL_TYPE> &matUmatImpl,
                         EntityHandle entity, int gg) override {
    MoFEMFunctionBegin;
    if (!matUmatImpl.umatInterfacePtr) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "UMAT interface is not initialised");
    }
    auto &umat = *matUmatImpl.umatInterfacePtr;
    umat.props = matUmatImpl.currentMaterialParameters;
    umat.noel = static_cast<int>(entity);
    umat.npt = gg + 1;
    CHKERR umat.setTimeData(matUmatImpl.currentTime,
                            matUmatImpl.currentTimeIncrement);
    CHKERR matUmatImpl.loadPreviousDeformationGradient(entity, gg);
    CHKERR matUmatImpl.loadPreviousStress(entity, gg);
    CHKERR matUmatImpl.loadPreviousStateVariables(entity, gg);
    CHKERR umat.setDeformationGradient(
        *matUmatImpl.matOpsDataPtr->getActiveDataPtr("F"),
        matUmatImpl.useDeformationGradient);
    CHKERR umat.setStrainIncrementFromDeformationGradient();
    CHKERR matUmatImpl.loadCoordinates(entity, gg);
    CHKERR matUmatImpl.loadPredef(entity, gg);
 
    CHKERR umat.callUmat();
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode
  convertStress(MatUmatImpl<DIM, MODEL_TYPE> &matUmatImpl,
                boost::shared_ptr<MatOpsData> mat_ops_data_ptr) override {
    MoFEMFunctionBegin;
    auto &umat = *matUmatImpl.umatInterfacePtr;
    auto set_stress_from_umat_cauchy = [&]() {
      MoFEMFunctionBegin;
 
      auto t_P = getFTensor2FromPtr<DIM, DIM>(
          mat_ops_data_ptr->getDependentDataPtr("P")->data().data());
 
      MatrixDouble sigma;
      CHKERR umat.getCauchyStressTensor(sigma);
 
      for (int i = 0; i != DIM; ++i)
        for (int Jidx = 0; Jidx != DIM; ++Jidx)
          t_P(i, Jidx) = sigma(i, Jidx);
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_stress_from_umat_cauchy();
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode
  convertTangent(MatUmatImpl<DIM, MODEL_TYPE> &matUmatImpl,
                 boost::shared_ptr<MatOpsData> mat_ops_data_ptr) override {
    MoFEMFunctionBegin;
    auto &umat = *matUmatImpl.umatInterfacePtr;
    auto set_tangent_from_umat_cauchy = [&]() {
      MoFEMFunctionBegin;
 
      auto t_dP = getFTensor4FromPtr<DIM, DIM, DIM, DIM>(
          mat_ops_data_ptr->getDependentDerivativesDataPtr("P_dF")
              ->data()
              .data());
 
      MatrixDouble d_sigma_d_eps;
      CHKERR umat.getCauchyTangentTensor(d_sigma_d_eps);
 
      for (int i = 0; i != DIM; ++i)
        for (int Jidx = 0; Jidx != DIM; ++Jidx)
          for (int m = 0; m != DIM; ++m)
            for (int N = 0; N != DIM; ++N)
              t_dP(i, Jidx, m, N) = d_sigma_d_eps(umatTensorIndex(i, Jidx),
                                                  umatTensorIndex(m, N));
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_tangent_from_umat_cauchy();
    MoFEMFunctionReturn(0);
  }
};
 
template <int DIM, int MODEL_TYPE>
struct HenckyStrainHandler : public StrainHandler<DIM, MODEL_TYPE> {
 
  HenckyStrainHandler(){};
 
  MoFEMErrorCode runUmat(MatUmatImpl<DIM, MODEL_TYPE> &matUmatImpl,
                         EntityHandle entity, int gg) override {
    MoFEMFunctionBegin;
    (void)matUmatImpl;
    (void)entity;
    (void)gg;
    MoFEMFunctionReturn(MOFEM_NOT_IMPLEMENTED);
  }
 
  MoFEMErrorCode
  convertStress(MatUmatImpl<DIM, MODEL_TYPE> &matUmatImpl,
                boost::shared_ptr<MatOpsData> mat_ops_data_ptr) override {
    MoFEMFunctionBegin;
    (void)matUmatImpl;
    (void)mat_ops_data_ptr;
    MoFEMFunctionReturn(MOFEM_NOT_IMPLEMENTED);
  }
 
  MoFEMErrorCode
  convertTangent(MatUmatImpl<DIM, MODEL_TYPE> &matUmatImpl,
                 boost::shared_ptr<MatOpsData> mat_ops_data_ptr) override {
    MoFEMFunctionBegin;
    (void)matUmatImpl;
    (void)mat_ops_data_ptr;
    MoFEMFunctionReturn(MOFEM_NOT_IMPLEMENTED);
  }
};
 
template <int DIM, int MODEL_TYPE>
struct FiniteStrainHandler : public StrainHandler<DIM, MODEL_TYPE> {
 
  FiniteStrainHandler(){};
 
  MoFEMErrorCode runUmat(MatUmatImpl<DIM, MODEL_TYPE> &matUmatImpl,
                         EntityHandle entity, int gg) override {
    MoFEMFunctionBegin;
    (void)matUmatImpl;
    (void)entity;
    (void)gg;
    MoFEMFunctionReturn(MOFEM_NOT_IMPLEMENTED);
  }
 
  MoFEMErrorCode
  convertStress(MatUmatImpl<DIM, MODEL_TYPE> &matUmatImpl,
                boost::shared_ptr<MatOpsData> mat_ops_data_ptr) override {
    MoFEMFunctionBegin;
    (void)matUmatImpl;
    (void)mat_ops_data_ptr;
    MoFEMFunctionReturn(MOFEM_NOT_IMPLEMENTED);
  }
 
  MoFEMErrorCode
  convertTangent(MatUmatImpl<DIM, MODEL_TYPE> &matUmatImpl,
                 boost::shared_ptr<MatOpsData> mat_ops_data_ptr) override {
    MoFEMFunctionBegin;
    (void)matUmatImpl;
    (void)mat_ops_data_ptr;
    MoFEMFunctionReturn(MOFEM_NOT_IMPLEMENTED);
  }
};
 
template <int DIM, int MODEL_TYPE>
struct MatUmatImpl : public PhysicalEquations {
  using PhysicalEquations::PhysicalEquations;
  using A = PhysicalEquations;
  using UmatInterfaceType = MoFEM::UmatInterface<MODEL_TYPE>;
  static constexpr auto DFGRD0_TAG_NAME = "UMAT_DFGRD0";
  static constexpr auto STRESS_TAG_NAME = "UMAT_STRESS";
  static constexpr auto STATEV_TAG_NAME = "UMAT_STATEV";
 
  static std::string getStateVariableTagName(const int index) {
    return (boost::format("%s_%d") % STATEV_TAG_NAME % index).str();
  }
 
  inline static bool useDeformationGradient =
      false; // by default use common data pass gradient of displacement
 
  ForcesAndSourcesCore::UserDataOperator *
  createOp(boost::shared_ptr<PhysicalEquations> physical_ptr, bool eval_stress,
           bool eval_tangent, bool update) override;
 
  MoFEMErrorCode getOptions(MoFEM::Interface *m_field_ptr = nullptr) override {
    MoFEMFunctionBegin;
 
    mField = m_field_ptr;
 
    MOFEM_LOG_CHANNEL("WORLD");
 
    char umat_library[PETSC_MAX_PATH_LEN] = "./umat.so";
    char umat_name[80] = "UMAT";
    PetscInt nstatev = 0;
    PetscInt nprops = 2;
    const char *list_strain_type[] = {"small", "hencky", "finite"};
    PetscInt choice_strain = StrainType::SMALL;
    PetscOptionsBegin(PETSC_COMM_WORLD, "umat_", "", "none");
    CHKERR PetscOptionsString("-library", "Path to UMAT shared library", "",
                              umat_library, umat_library, sizeof(umat_library),
                              PETSC_NULLPTR);
    CHKERR PetscOptionsString("-name", "UMAT material name", "", umat_name,
                              umat_name, sizeof(umat_name), PETSC_NULLPTR);
    CHKERR PetscOptionsInt("-nstatev", "Number of UMAT state variables", "",
                           nstatev, &nstatev, PETSC_NULLPTR);
    CHKERR PetscOptionsInt("-nprops", "Number of UMAT material properties", "",
                           nprops, &nprops, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-young_modulus", "Default first property", "",
                              youngModulus, &youngModulus, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-poisson_ratio", "Default second property", "",
                              poissonRatio, &poissonRatio, PETSC_NULLPTR);
    CHKERR PetscOptionsEList(
        "-strain_type", "Strain type", "", list_strain_type, 3,
        list_strain_type[choice_strain], &choice_strain, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    strainType = static_cast<StrainType>(choice_strain);
 
    switch (strainType) {
    case SMALL:
      strainHandler = std::make_unique<SmallStrainHandler<DIM, MODEL_TYPE>>();
      break;
    case HENCKY:
      strainHandler = std::make_unique<HenckyStrainHandler<DIM, MODEL_TYPE>>();
      break;
    case FINITE:
      strainHandler = std::make_unique<FiniteStrainHandler<DIM, MODEL_TYPE>>();
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unsupported UMAT strain type");
    }
 
    if (nprops <= 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "UMAT requires a positive number of properties");
    }
 
    if (nstatev < 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "UMAT requires a non-negative number of state variables");
    }
 
    numberStateVariables = static_cast<int>(nstatev);
    stateVar.resize(numberStateVariables, 0.0);
    numberMaterialProperties = static_cast<int>(nprops);
    materialName = umat_name;
    libraryPath = umat_library;
 
    defaultMaterialParameters.assign(numberMaterialProperties, 0.0);
    defaultMaterialParameters[0] = youngModulus;
    if (numberMaterialProperties > 1)
      defaultMaterialParameters[1] = poissonRatio;
    currentMaterialParameters = defaultMaterialParameters;
 
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "Default Umat parameters:")
        << " umat_library = " << libraryPath << " name = " << materialName
        << " umat_nstatev = " << numberStateVariables
        << " umat_nprops = " << numberMaterialProperties;
 
    std::string block_name = "MAT_UMAT";
    for (auto &m :
         m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
             std::regex((boost::format("%s(.*)") % block_name).str()))) {
      std::vector<double> block_data;
      CHKERR m->getAttributes(block_data);
      if ((int)block_data.size() < numberMaterialProperties) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "MAT_UMAT block requires at least %d attributes", nprops);
      }
 
      auto get_block_ents = [&]() {
        Range ents;
        CHK_MOAB_THROW(m_field_ptr->get_moab().get_entities_by_handle(
                           m->meshset, ents, true),
                       "can not get block entities");
        return ents;
      };
 
      A::paramVecByRange.push_back(
          {get_block_ents(),
           std::vector<double>(block_data.begin(),
                               block_data.begin() + numberMaterialProperties)});
 
      std::ostringstream properties_stream;
      properties_stream << "[";
      for (int ii = 0; ii != numberMaterialProperties; ++ii) {
        if (ii)
          properties_stream << ", ";
        properties_stream << block_data[ii];
      }
      properties_stream << "]";
 
      MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "MatBlock for Umat")
          << *m << " properties = " << properties_stream.str();
    }
 
    umatInterfacePtr = std::make_unique<UmatInterfaceType>(
        numberStateVariables, numberMaterialProperties);
    CHKERR umatInterfacePtr->initialise(libraryPath, materialName);
 
    if (m_field_ptr) {
      CHKERR A::matOpsDataPtr->bindStateTag(*m_field_ptr, DFGRD0_TAG_NAME,
                                            DIM * DIM);
      CHKERR A::matOpsDataPtr->bindStateTag(*m_field_ptr, STRESS_TAG_NAME,
                                            UmatInterfaceType::NTENS);
      if (numberStateVariables > 0) {
        for (int ii = 0; ii != numberStateVariables; ++ii) {
          CHKERR A::matOpsDataPtr->bindStateTag(*m_field_ptr,
                                                getStateVariableTagName(ii), 1);
        }
      }
      CHKERR bindAuxiliaryStateTags(*m_field_ptr);
      CHKERR A::matOpsDataPtr->setupStateData();
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode setParams(FEMethod *fe_ptr, int gg) override {
    MoFEMFunctionBeginHot;
    (void)gg;
    currentTimeStep = fe_ptr->ts_step;
    currentTime = fe_ptr->ts_t;
    currentTimeIncrement = fe_ptr->ts_dt;
    const auto ent = fe_ptr->getFEEntityHandle();
    currentMaterialParameters = defaultMaterialParameters;
    for (auto &[range, param_vec] : A::paramVecByRange) {
      if (range.find(ent) != range.end()) {
        currentMaterialParameters = param_vec;
        break;
      }
    }
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode recordTape() override {
    MoFEMFunctionBegin;
 
    A::matOpsDataPtr->insertCommonData("grad", MatrixDouble());
    A::matOpsDataPtr->insertCommonData("P", MatrixDouble());
    A::matOpsDataPtr->insertCommonData("P_dF", MatrixDouble());
 
    A::matOpsDataPtr->insertActiveData("F", MatrixDouble());
    A::matOpsDataPtr->insertDependentData("P", MatrixDouble());
    A::matOpsDataPtr->insertDependentDerivativesData("P_dF", MatrixDouble());
 
    A::matOpsDataPtr->getActiveDataPtr("F")->resize(DIM, DIM, false);
    A::matOpsDataPtr->getDependentDataPtr("P")->resize(DIM, DIM, false);
    A::matOpsDataPtr->getDependentDerivativesDataPtr("P_dF")->resize(
        DIM * DIM, DIM * DIM, false);
 
    CHKERR recordAuxiliaryMatOpsData();
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode evaluateVariable(int, EntityHandle entity, int gg) override {
    MoFEMFunctionBegin;
    CHKERR strainHandler->runUmat(*this, entity, gg);
    CHKERR strainHandler->convertStress(*this, A::matOpsDataPtr);
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode evaluateDerivatives(int, EntityHandle entity,
                                     int gg) override {
    MoFEMFunctionBegin;
    CHKERR strainHandler->runUmat(*this, entity, gg);
    CHKERR strainHandler->convertTangent(*this, A::matOpsDataPtr);
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode updateState(int, EntityHandle entity, int gg) override {
    MoFEMFunctionBegin;
    CHKERR strainHandler->runUmat(*this, entity, gg);
    CHKERR storeCurrentDeformationGradient(entity, gg);
    CHKERR storeCurrentStress(entity, gg);
    CHKERR storeCurrentStateVariables(entity, gg);
    MoFEMFunctionReturn(0);
  }
 
protected:
  virtual MoFEMErrorCode bindAuxiliaryStateTags(MoFEM::Interface &m_field) {
    MoFEMFunctionBeginHot;
    (void)m_field;
    MoFEMFunctionReturnHot(0);
  }
 
  virtual MoFEMErrorCode recordAuxiliaryMatOpsData() {
    MoFEMFunctionBeginHot;
    MoFEMFunctionReturnHot(0);
  }
 
  // TODO: Populate UMAT COORDS from integration-point coordinates when needed.
  // The default keeps COORDS zero; coordinate-dependent UMATs should override it.
  virtual MoFEMErrorCode loadCoordinates(EntityHandle entity, int gg) {
    MoFEMFunctionBeginHot;
    (void)entity;
    (void)gg;
    for (int dd = 0; dd != 3; ++dd)
      umatInterfacePtr->coords[dd] = 0.0;
    MoFEMFunctionReturnHot(0);
  }
 
  virtual MoFEMErrorCode loadPredef(EntityHandle entity, int gg) {
    MoFEMFunctionBeginHot;
    (void)entity;
    (void)gg;
    if (umatInterfacePtr->predef.empty()) {
      umatInterfacePtr->predef.push_back(0.0);
    }
    if (umatInterfacePtr->dpred.empty()) {
      umatInterfacePtr->dpred.push_back(0.0);
    }
    std::fill(umatInterfacePtr->predef.begin(), umatInterfacePtr->predef.end(),
              0.0);
    std::fill(umatInterfacePtr->dpred.begin(), umatInterfacePtr->dpred.end(),
              0.0);
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode loadPreviousDeformationGradient(EntityHandle entity, int gg) {
    MoFEMFunctionBegin;
 
    CHKERR umatInterfacePtr->setIdentityDeformationGradient();
 
    auto state_ptr =
        A::matOpsDataPtr->getStateDataPtr(DFGRD0_TAG_NAME, entity, gg);
    bool has_stored_gradient = false;
    for (int row = 0; row != DIM; ++row) {
      for (int col = 0; col != DIM; ++col) {
        const double value = (*state_ptr)(row * DIM + col, 0);
        umatInterfacePtr->dfgrd0[row + 3 * col] = value;
        has_stored_gradient = has_stored_gradient || std::abs(value) > 0;
      }
    }
 
    if (!has_stored_gradient) {
      CHKERR umatInterfacePtr->setIdentityDeformationGradient();
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode storeCurrentDeformationGradient(EntityHandle entity, int gg) {
    MoFEMFunctionBegin;
 
    auto state_ptr =
        A::matOpsDataPtr->getStateDataPtr(DFGRD0_TAG_NAME, entity, gg);
    for (int row = 0; row != DIM; ++row) {
      for (int col = 0; col != DIM; ++col) {
        (*state_ptr)(row *DIM + col, 0) =
            umatInterfacePtr->dfgrd1[row + 3 * col];
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode loadPreviousStress(EntityHandle entity, int gg) {
    MoFEMFunctionBegin;
 
    auto state_ptr =
        A::matOpsDataPtr->getStateDataPtr(STRESS_TAG_NAME, entity, gg);
    for (int ii = 0; ii != UmatInterfaceType::NTENS; ++ii) {
      umatInterfacePtr->stress[ii] = (*state_ptr)(ii, 0);
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode loadPreviousStateVariables(EntityHandle entity, int gg) {
    MoFEMFunctionBegin;
 
    umatInterfacePtr->statev.assign(numberStateVariables, -1);
    if (!numberStateVariables) {
      MoFEMFunctionReturnHot(0);
    }
 
    for (int ii = 0; ii != numberStateVariables; ++ii) {
      auto state_ptr = A::matOpsDataPtr->getStateDataPtr(
          getStateVariableTagName(ii), entity, gg);
      umatInterfacePtr->statev[ii] = (*state_ptr)(0, 0);
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode storeCurrentStress(EntityHandle entity, int gg) {
    MoFEMFunctionBegin;
 
    auto state_ptr =
        A::matOpsDataPtr->getStateDataPtr(STRESS_TAG_NAME, entity, gg);
    for (int ii = 0; ii != UmatInterfaceType::NTENS; ++ii) {
      (*state_ptr)(ii, 0) = umatInterfacePtr->stress[ii];
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode storeCurrentStateVariables(EntityHandle entity, int gg) {
    MoFEMFunctionBegin;
 
    if (!numberStateVariables) {
      MoFEMFunctionReturnHot(0);
    }
 
    for (int ii = 0; ii != numberStateVariables; ++ii) {
      auto state_ptr = A::matOpsDataPtr->getStateDataPtr(
          getStateVariableTagName(ii), entity, gg);
      (*state_ptr)(0, 0) = umatInterfacePtr->statev[ii];
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEM::Interface *mField = nullptr;
  std::unique_ptr<UmatInterfaceType> umatInterfacePtr;
  std::string libraryPath = "./umat.so";
  std::string materialName = "UMAT";
  StrainType strainType = StrainType::SMALL;
  PetscInt currentTimeStep = 0;
  PetscReal currentTime = 0;
  PetscReal currentTimeIncrement = 0;
  double youngModulus = 1.0;
  double poissonRatio = 0.3;
  int numberStateVariables = 1;
  int numberMaterialProperties = 2;
  std::vector<double> defaultMaterialParameters = {youngModulus, poissonRatio};
  std::vector<double> currentMaterialParameters = defaultMaterialParameters;
  std::vector<double> stateVar = std::vector<double>(numberStateVariables, 0.0);
  friend struct SmallStrainHandler<DIM, MODEL_TYPE>;
  friend struct HenckyStrainHandler<DIM, MODEL_TYPE>;
  friend struct FiniteStrainHandler<DIM, MODEL_TYPE>;
  std::unique_ptr<StrainHandler<DIM, MODEL_TYPE>> strainHandler;
};
 
} // namespace MatOps
 
#endif // MAT_UMAT_IMPL_HPP