EshelbianMatCore.cpp

Go to the documentation of this file.

/**
 * @file EshelbianMatCore.cpp
 * @author your name (you@domain.com)
 * @brief
 * @version 0.1
 * @date 2026-04-14
 *
 * @copyright Copyright (c) 2026
 *
 */
 
namespace EshelbianPlasticity {
 
struct MatPhysicalEquations : public PhysicalEquations {
 
  MatPhysicalEquations(
      boost::shared_ptr<MatOps::PhysicalEquations> mat_physical_equations_ptr)
      : PhysicalEquations(),
        matPhysicalEquationsPtr(mat_physical_equations_ptr) {
    CHK_THROW_MESSAGE(getOptions(), "get options failed");
  }
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
    PetscOptionsBegin(PETSC_COMM_WORLD, "meta_", "", "none");
    alphaGradU = 0;
    CHKERR PetscOptionsScalar("-viscosity_alpha_grad_u", "viscosity", "",
                              alphaGradU, &alphaGradU, PETSC_NULLPTR);
    PetscOptionsEnd();
    MOFEM_LOG("EP", Sev::inform) << "Mat core grad alpha u = " << alphaGradU;
    MoFEMFunctionReturn(0);
  }
 
  UserDataOperator *
  returnOpJacobian(const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr) override {
    return matPhysicalEquationsPtr->createOp(matPhysicalEquationsPtr, eval_rhs,
                                             eval_lhs, false);
  }
 
  VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha) override {
    return new OpSpatialPhysical(field_name, matPhysicalEquationsPtr, alpha,
                                 alphaGradU);
  }
 
  VolUserDataOperator *returnOpSpatialPhysical_du_du(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha) {
    return new OpSpatialPhysical_du_du(row_field, col_field,
                                       matPhysicalEquationsPtr, alpha,
                                       alphaGradU);
  }
 
  struct OpSpatialPhysicalExternalStrain : public OpAssembleVolume {
 
    OpSpatialPhysicalExternalStrain(
        const std::string &field_name,
        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
        boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
        std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
 
    MoFEMErrorCode integrate(EntData &data);
 
  private:
    boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
    std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
  };
 
  VolUserDataOperator *returnOpSpatialPhysicalExternalStrain(
      const std::string &field_name,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv) {
    return new OpSpatialPhysicalExternalStrain(field_name, data_ptr,
                                               external_strain_vec_ptr, smv);
  }
 
protected:
  boost::shared_ptr<MatOps::PhysicalEquations> matPhysicalEquationsPtr;
  double alphaGradU = 0;
 
  struct OpSpatialPhysical : public OpAssembleVolume {
 
    OpSpatialPhysical(const std::string &field_name,
                      boost::weak_ptr<MatOps::PhysicalEquations> physics_ptr,
                      const double alpha, const double alpha_grad_u)
        : OpAssembleVolume(field_name, nullptr, OPROW), alphaU(alpha),
          alphaGradU(alpha_grad_u), physicsPtr(physics_ptr) {}
 
    MoFEMErrorCode integrate(EntData &data);
 
  private:
    const double alphaU;
    const double alphaGradU;
    boost::shared_ptr<MatOps::PhysicalEquations>
        physicsPtr; ///< material physical equations
  };
 
  struct OpSpatialPhysical_du_du : public OpAssembleVolume {
    OpSpatialPhysical_du_du(
        std::string row_field, std::string col_field,
        boost::shared_ptr<MatOps::PhysicalEquations> physics_ptr,
        const double alpha, const double alpha_grad_u)
        : OpAssembleVolume(row_field, col_field, nullptr, OPROWCOL, false),
          alphaU(alpha), alphaGradU(alpha_grad_u), physicsPtr(physics_ptr) {
      sYmm = false;
    }
 
    MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
  private:
    const double alphaU;
    const double alphaGradU;
    boost::shared_ptr<MatOps::PhysicalEquations>
        physicsPtr; ///< material physical equations
  };
};
 
MoFEMErrorCode
MatPhysicalEquations::OpSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  auto physics_ptr = physicsPtr;
  auto mat_ops_data = physics_ptr->matOpsDataPtr;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
 
  auto t_P = getFTensor2FromMat<3, 3, -1, DL>(
      mat_ops_data->getCommonDataPtr("PAtPts"));
  auto t_approx_P_adjoint_log_du = getFTensor1FromMat<size_symm, -1, DL>(
      mat_ops_data->getCommonDataPtr("adjointPdUAtPts"));
  auto t_diff_u = getFTensor4FromMat<3, 3, 3, 3, -1, DL>(
      mat_ops_data->getCommonDataPtr("diffStretchH1AtPts"));
  auto t_dot_log_u = getFTensor2SymmetricFromMat<3, -1, DL>(
      mat_ops_data->getCommonDataPtr("logStretchDotTensorAtPts"));
  auto t_grad_dot_log_u = getFTensor2FromMat<size_symm, 3, -1, DL>(
      mat_ops_data->getCommonDataPtr("gradLogStretchDotTensorAtPts"));
 
  FTENSOR_INDEX(3, i);
  FTENSOR_INDEX(3, j);
  FTENSOR_INDEX(3, k);
  FTENSOR_INDEX(3, l);
  FTENSOR_INDEX(3, m);
  FTENSOR_INDEX(3, n);
 
  auto get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  auto t_row_grad_fun = data.getFTensor1DiffN<3>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor2(nF);
 
    FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
    t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
 
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) =
        a * (t_approx_P_adjoint_log_du(L) - t_Ldiff_u(i, j, L) * t_P(i, j));
    t_residual(L) += (a * alphaU) * (t_dot_log_u(i, j) * t_L(i, j, L));
 
    FTensor::Tensor2<double, size_symm, 3> t_grad_residual;
    t_grad_residual(L, i) = (a * alphaGradU) * t_grad_dot_log_u(L, i);
 
    int bb = 0;
    for (; bb != nb_dofs / size_symm; ++bb) {
      t_nf(L) -= t_row_base_fun * t_residual(L);
      t_nf(L) += t_row_grad_fun(i) * t_grad_residual(L, i);
      ++t_nf;
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
    for (; bb != nb_base_functions; ++bb) {
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
 
    ++t_w;
    ++t_approx_P_adjoint_log_du;
    ++t_P;
    ++t_diff_u;
    ++t_dot_log_u;
    ++t_grad_dot_log_u;
  }
  MoFEMFunctionReturn(0);
}
 
 
MoFEMErrorCode
MatPhysicalEquations::OpSpatialPhysical_du_du::integrate(EntData &row_data,
                                                           EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  FTensor::Index<'J', size_symm> J;
  auto t_L = FTensor::SymmLTensor<double, 3>();
  auto t_diff = FTensor::DiffTensor<double>();
  constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
 
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 6>, size_symm,
                            size_symm>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
        &m(r + 0, c + 4), &m(r + 0, c + 5),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
        &m(r + 1, c + 4), &m(r + 1, c + 5),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
        &m(r + 2, c + 4), &m(r + 2, c + 5),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
        &m(r + 3, c + 4), &m(r + 3, c + 5),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
        &m(r + 4, c + 4), &m(r + 4, c + 5),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
        &m(r + 5, c + 4), &m(r + 5, c + 5)
 
    );
  };
 
  FTENSOR_INDEX(3, i);
  FTENSOR_INDEX(3, j);
  FTENSOR_INDEX(3, k);
  FTENSOR_INDEX(3, l);
  FTENSOR_INDEX(3, m);
  FTENSOR_INDEX(3, n);
 
  MatrixDouble dP;
 
  auto get_dP = [&]() {
    auto get_tensor_for_dP =
        MatrixSizeHelper<GetFTensor2FromMatType<size_symm, size_symm, -1, DL>,
                         DL>::size(dP, nb_integration_pts);
 
    auto mat_ops_data = physicsPtr->matOpsDataPtr;
 
    auto t_P = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("PAtPts"));
    auto t_P_du = getFTensor4FromMat<3, 3, 3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("PAtPts_du"));
    auto t_approx_P_adjoint_dstretch = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("adjointPdstretchAtPts"));
    auto t_u = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("stretchH1AtPts"));
    auto t_diff_u = getFTensor4FromMat<3, 3, 3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("diffStretchH1AtPts"));
    auto t_grad_h1 = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("wGradH1AtPts"));
    auto t_eigen_vals = getFTensor1FromMat<3, -1, DL>(
        mat_ops_data->getCommonDataPtr("eigenVals"));
    auto t_eigen_vecs = getFTensor2FromMat<3, 3, -1, DL>(
        mat_ops_data->getCommonDataPtr("eigenVecs"));
 
    auto ts_a = getTSa();
    auto t_kd = FTensor::Kronecker_Delta<double>();
 
    auto t_dP = get_tensor_for_dP();
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
      FTensor::Tensor2<double, 3, 3> t_h1;
 
      switch (EshelbianCore::gradApproximator) {
      case LARGE_ROT:
      case MODERATE_ROT:
        t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
        break;
      default:
        t_h1(i, j) = t_kd(i, j);
      }
 
      FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
      t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
      t_dP(L, J) =
          -t_Ldiff_u(i, j, L) * (t_P_du(i, j, k, l) * t_Ldiff_u(k, l, J));
      t_dP(L, J) += (ts_a * alphaU) *
                    (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
 
      if (EshelbianCore::stretchSelector > LINEAR &&
          EshelbianCore::gradApproximator > MODERATE_ROT) {
        FTensor::Tensor2<double, 3, 3> t_deltaP;
        t_deltaP(i, j) =
            t_approx_P_adjoint_dstretch(i, j) - t_P(i, k) * t_h1(j, k);
        FTensor::Tensor2_symmetric<double, 3> t_deltaP_sym;
        t_deltaP_sym(i, j) = (t_deltaP(i, j) || t_deltaP(j, i));
        t_deltaP_sym(i, j) /= 2.0;
        auto nb_uniq =
            getUniqNb<3>(getVectorAdaptor(&t_eigen_vals(0), 3),
                         FTensor::Number<3>{});
        auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_deltaP_sym, nb_uniq);
        t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP2(i, j, k, l) * t_L(k, l, J));
      }
 
      ++t_P;
      ++t_dP;
      ++t_P_du;
      ++t_approx_P_adjoint_dstretch;
      ++t_u;
      ++t_diff_u;
      ++t_grad_h1;
      ++t_eigen_vals;
      ++t_eigen_vecs;
    }
 
    return get_tensor_for_dP();
  };
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
  auto t_row_grad_fun = row_data.getFTensor1DiffN<3>();
 
  auto t_dP = get_dP();
  auto v = getVolume();
  auto ts_a = getTSa();
  auto t_w = getFTensor0IntegrationWeight();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    int rr = 0;
    for (; rr != row_nb_dofs / size_symm; ++rr) {
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_col_grad_fun = col_data.getFTensor1DiffN<3>(gg, 0);
      auto t_m = get_ftensor2(K, size_symm * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / size_symm; ++cc) {
        const double b = a * t_row_base_fun * t_col_base_fun;
        t_m(L, J) -= b * t_dP(L, J);
        const double c = (a * alphaGradU * ts_a) *
                         (t_row_grad_fun(i) * t_col_grad_fun(i));
        t_m(L, J) += c * t_kd_sym(L, J);
        ++t_m;
        ++t_col_base_fun;
        ++t_col_grad_fun;
      }
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
 
    ++t_w;
    ++t_dP;
  }
  MoFEMFunctionReturn(0);
}
 
MatPhysicalEquations::OpSpatialPhysicalExternalStrain::
    OpSpatialPhysicalExternalStrain(
        const std::string &field_name,
        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
        boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
        std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
    : OpAssembleVolume(field_name, data_ptr, OPROW),
      externalStrainVecPtr(external_strain_vec_ptr), scalingMethodsMap{smv} {}
 
MoFEMErrorCode MatPhysicalEquations::OpSpatialPhysicalExternalStrain::integrate(
    EntData &data) {
  MoFEMFunctionBegin;
 
 
 
  double time = OpAssembleVolume::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
  // get entity of tet
  EntityHandle fe_ent = OpAssembleVolume::getFEEntityHandle();
  // iterate over all block data
 
  for (auto &ext_strain_block : (*externalStrainVecPtr)) {
    auto block_name = "(.*)ANALYTICAL_EXTERNALSTRAIN(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(ext_strain_block.blockName, reg_name)) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Analytical external strain not implemented for Neo-Hookean "
              "material.");
    }
    // check if finite element entity is part of the EXTERNALSTRAIN block
    if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
      double scale = 1;
      if (scalingMethodsMap.find(ext_strain_block.blockName) !=
          scalingMethodsMap.end()) {
        scale *=
            scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
      } else {
        MOFEM_LOG("SELF", Sev::warning)
            << "No scaling method found for " << ext_strain_block.blockName;
      }
 
      // get ExternalStrain block data
      double external_strain_val = scale * ext_strain_block.val;
      double K = ext_strain_block.bulkModulusK;
 
      FTensor::Index<'L', size_symm> L;
      auto t_L = FTensor::SymmLTensor<double, 3>();
      constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = data.getN().size1();
      auto vol = getVolume();
      auto t_w = getFTensor0IntegrationWeight();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
      FTENSOR_INDEX(SPACE_DIM, k);
      FTENSOR_INDEX(SPACE_DIM, l);
      FTENSOR_INDEX(SPACE_DIM, m);
      FTENSOR_INDEX(SPACE_DIM, n);
 
      int nb_base_functions = data.getN().size2();
      auto t_row_base_fun = data.getFTensor0N();
 
      const double bulk_modulus = K;
      const double diag_val =
          EshelbianCore::inv_f(std::exp(external_strain_val));
      auto fun_neohookean_bulk = [](double K, double tr) { return K * tr; };
      const double sigma_J = fun_neohookean_bulk(bulk_modulus, 3 * diag_val);
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        double a = vol * t_w;
        ++t_w;
 
        FTensor::Tensor1<double, size_symm> t_residual;
        t_residual(L) = (t_L(i, j, L) * t_kd(i, j)) * sigma_J;
        t_residual(L) *= a;
 
        auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
        int bb = 0;
        for (; bb != nb_dofs / size_symm; ++bb) {
          t_nf(L) -= t_row_base_fun * t_residual(L);
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb) {
          ++t_row_base_fun;
        }
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity