EshelbianInterface.hpp

Go to the documentation of this file.

/**
 * @file EshelbianInterface.hpp
 * @author your name (you@domain.com)
 * @brief
 * @version 0.1
 * @date 2025-11-20
 *
 * @copyright Copyright (c) 2025
 *
 */
 
#include "EshelbianOperators.hpp"
#include "MoFEM.hpp"
 
#ifndef ESHELBIANINTERFACE_HPP
  #define ESHELBIANINTERFACE_HPP
 
namespace EshelbianPlasticity {
 
struct OpGetTagData : public UserDataOperator {
  OpGetTagData(boost::shared_ptr<DataAtIntegrationPts> data_ptr, Tag th,
               boost::shared_ptr<Range> ents_ptr)
      : UserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr), tagHandle(th),
        entsPtr(ents_ptr) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
    auto fe_handle = getFEEntityHandle();
    if (entsPtr->contains(fe_handle) == false) {
      MoFEMFunctionReturn(0);
    }
    auto nb_gauss_pts = getGaussPts().size2();
    auto data_ptr = dataAtPts->getTagData(tagHandle);
    auto &moab = m_field.getMoab();
    const void *tag_data[] = {nullptr};
    const void tag_sizes[] = {0};
    CHKERR moab.tag_get_by_ptr(tagHandle, &fe_handle, 1, tag_data, tag_sizes);
    if (tag_sizes[0] % nb_gauss_pts != 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Tag length not compatible with number of gauss pts");
    }
    data_ptr->resize(tag_sizes[0] / nb_gauss_pts, nb_gauss_pts, false);
    std::copy(static_cast<const double *>(tag_data[0]),
              static_cast<const double *>(tag_data[0]) + tag_sizes[0],
              data_ptr->data().begin());
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  Tag tagHandle;
  boost::shared_ptr<Range> entsPtr;
};
 
struct OpSetTagData : public UserDataOperator {
  OpSetTagData(boost::shared_ptr<DataAtIntegrationPts> data_ptr, Tag th,
               boost::shared_ptr<Range> ents_ptr)
      : UserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr), tagHandle(th),
        entsPtr(ents_ptr) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
    auto fe_handle = getFEEntityHandle();
    if (entsPtr->contains(fe_handle) == false) {
      MoFEMFunctionReturn(0);
    }
    auto data_ptr = dataAtPts->getTagData(tagHandle);
    auto &moab = m_field.getMoab();
    int tag_length = data_ptr->data().size();
    CHKERR moab.tag_set_by_ptr(tagHandle, &fe_handle, 1,
                               data_ptr->data().data(), tag_length);
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  Tag tagHandle;
  boost::shared_ptr<Range> entsPtr;
};
 
struct CrackInteface {
  virtual MatrixDouble &
  getGapAtPts(boost::shared_ptr<DataAtIntegrationPts> data_ptr, ,
              MatrixDouble &normalsAtGaussPts, int sense) = 0;
  virtual MatrixDouble &
  getInvStiffnessAtPts(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                       MatrixDouble &normalsAtGaussPts, int sense) = 0;
};
 
struct ExponentialGriffithReleaseData : public CrackInteface {
 
  MatrixDouble &getGapAtPts(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                            MatrixDouble &normalsAtGaussPts, int sense) {
 
    auto impl = [&]() {
      MoFEMFunctionBegin;
      FTENSOR_INDEX(3, i);
      FTENSOR_INDEX(3, J);
 
      auto nb_gauss_pts = data_ptr->approxPAtPts.size2();
      auto t_approx_P = getFTensor2FromMat<3, 3>(data_ptr->approxPAtPts);
      double *ptr = normalsAtGaussPts().data().data();
      FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_normal(ptr, &ptr[1],
                                                                  &ptr[2]);
      gapAtPts.resize(3, nb_gauss_pts, false);
      auto t_delta_gap = getFTensor2FromMat<3, 3>(gapAtPts);
      for (size_t gp = 0; gp < nb_gauss_pts; ++gp) {
        auto nrm2 = t_normal.l2();
        t_traction(i) = t_approx_P(i, J) * (sense * t_normal(J) / nrm2);
        t_gap(i) = t_traction(i) * (h / E);
        ++t_approx_P;
        ++t_normal;
        ++t_gap;
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(impl(), "Error in ExponentialGriffithReleaseData");
 
    return gapAtPts;
  }
 
  MatrixDouble &
  getInvStiffnessAtPts(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                       MatrixDouble &normalsAtGaussPts, int sense) {
    MoFEMFunctionBegin;
    constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
    auto nb_gauss_pts = data_ptr->approxPAtPts.size2();
    double *ptr = normalsAtGaussPts().data().data();
    FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_normal(ptr, &ptr[1],
                                                                &ptr[2]);
    invStiffnessAtPts.resize(3, 3 * nb_gauss_pts, false);
    auto t_inv_D = getFTensor2FromMat<3, 3>(invStiffnessAtPts);
    for (size_t gp = 0; gp < nb_gauss_pts; ++gp) {
      t_inv_D(i, j) = t_kd(i, j) * (h / E);
      ++t_inv_D;
    }
    MoFEMFunctionReturn(0);
  }
 
protected:
  double Gc; //< Griffith energy release rate
  double E;  //< Young modulus
  double h;  //< Characteristic length
  double ft; //< Tensile strength
 
  MatrixDouble invStiffnessAtPts;
  MatrixDouble gapAtPts;
};
 
struct OpInterfaceRhs : public FormsIntegrators<FaceUserDataOperator>::Assembly<
                            PETSC>::OpBrokenBase {
 
  using OP = typename FormsIntegrators<FaceUserDataOperator>::Assembly<
      PETSC>::OpBrokenBase;
 
  OpInterfaceRhs(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      std::string row_field, boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<CrackInterface> release_data_ptr,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, ents_ptr), dataAtPts(data_ptr),
        releaseDataPtr(release_data_ptr) {}
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        EntitiesFieldData::EntData &row_data) {
 
    MoFEMFunctionBegin;
 
    if (OP::entsPtr) {
      if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
        MoFEMFunctionReturnHot(0);
    }
 
  #ifndef NDEBUG
    if (!brokenBaseSideData) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
    }
  #endif // NDEBUG
 
    auto do_work_rhs = [this](int row_side, EntityType row_type,
                              EntitiesFieldData::EntData &row_data, int sense) {
      MoFEMFunctionBegin;
      // get number of dofs on row
      OP::nbRows = row_data.getIndices().size();
      if (!OP::nbRows)
        MoFEMFunctionReturnHot(0);
      // get number of integration points
      OP::nbIntegrationPts = OP::getGaussPts().size2();
      // get row base functions
      OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
      // resize and clear the right hand side vector
      OP::locF.resize(OP::nbRows, false);
      OP::locF.clear();
      // integrate local vector
      this->rowSense = bd.getSense();
      CHKERR this->iNtegrate(row_data);
      // assemble local vector
      CHKERR this->aSsemble(row_data);
      MoFEMFunctionReturn(0);
    };
 
    switch (OP::opType) {
    case OP::OPROW:
      SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
              "LHS not implemented for OpEvalExponentialGriffithReleaseRhs");
      break;
    case OP::OPSPACE:
      for (auto &bd : *brokenBaseSideData) {
        CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData(),
                           bd.getSense());
      }
      break;
    default:
      CHK_MOAB_THROW(MOFEM_IMPOSSIBLE_CASE,
                     (std::string("wrong op type ") +
                      OpBaseDerivativesBase::OpTypeNames[OP::opType])
                         .c_str());
    }
 
    MoFEMFunctionReturn(0);
  };
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, J);
 
    auto nb_gauss_pts = getGaussPts().size2();
    int nb_dofs = data.getIndices().size();
    OP::locF.resize(nb_dofs, false);
    OP::locF.clear();
 
    auto &gap_map = releaseDataPtr->getGapAtPts(
        dataAtPts, getNormalsAtGaussPts(), rowSense);
    t_delta_gap = getFTensor1FromMat<3>(gap_map);
 
    auto area = OP::getMeasure();
    auto t_w = OP::getFTensor0IntegrationWeight();
    auto t_normal = getFTensor1NormalsAtGaussPts();
    auto t_row_base = data.getFTensor1N<3>();
    auto sense = bd.getSense();
 
    for (auto gp = 0; gp < nb_gauss_pts; ++gp) {
      double alpha = t_w / 2.;
      int bb = 0;
      auto t_rhs = getFTensor1FromPtr<3>(OP::locF.data().data());
      for (; bb != nb_dofs / SPACE_DIM; ++bb) {
        t_rhs(i) +=
            alpha * (t_row_base(J) * (rowSense * t_normal(J))) * t_delta_gap(i);
        ++t_rhs;
        ++t_row_base;
      }
      ++t_w;
      ++t_normal;
      ++t_delta_gap;
    }
 
    MoFEMFunctionReturn(0);
  }
 
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  boost::shared_ptr<CrackInterface> releaseDataPtr;
  int rowSense;
};
 
struct OpInterfaceLhs : public FormsIntegrators<FaceUserDataOperator>::Assembly<
                          PETSC>::OpBrokenBase {
 
  using OP = typename FormsIntegrators<FaceUserDataOperator>::Assembly<
      PETSC>::OpBrokenBase;
 
  OpInterfaceLhs(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<CrackInteface> release_data_ptr,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, ents_ptr), dataAtPts(data_ptr),
        releaseDataPtr(release_data_ptr) {}
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        EntitiesFieldData::EntData &row_data) {
    MoFEMFunctionBegin;
 
    if (OP::entsPtr) {
      if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
        MoFEMFunctionReturnHot(0);
    }
 
  #ifndef NDEBUG
    if (!brokenBaseSideData) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
    }
  #endif // NDEBUG
 
    auto do_work_lhs = [this](int row_side, int col_side, EntityType row_type,
                              EntityType col_type,
                              EntitiesFieldData::EntData &row_data,
                              EntitiesFieldData::EntData &col_data,
                              int row_sense, int col_sense) {
      MoFEMFunctionBegin;
 
      auto check_if_assemble_transpose = [&] {
        if (this->sYmm) {
          if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
            return true;
          else
            return false;
        } else if (OP::assembleTranspose) {
          return true;
        }
        return false;
      };
 
      OP::rowSide = row_side;
      OP::rowType = row_type;
      OP::colSide = col_side;
      OP::colType = col_type;
      OP::nbCols = col_data.getIndices().size();
      OP::locMat.resize(OP::nbRows, OP::nbCols, false);
      OP::locMat.clear();
      rowSense = row_sense;
      colSense = col_sense;
      CHKERR this->iNtegrate(row_data, col_data);
      CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
      MoFEMFunctionReturn(0);
    };
 
    switch (OP::opType) {
    case OP::OPSPACE:
 
      for (auto &bd : *brokenBaseSideData) {
 
  #ifndef NDEBUG
        if (bd.getData().getSpace() != HDIV &&
            bd.getData().getSpace() != HCURL) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "%s",
                  (std::string("Expect Hdiv or Hcurl space but received ") +
                   FieldSpaceNames[bd.getData().getSpace()])
                      .c_str());
        }
        if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "base functions not set");
        }
  #endif
 
        // Iterate over diagonal. Both sides of interface are coupled by
        // Lagrange multiplier.
        CHKERR do_work_lhs(
 
            // side
            bd.getSide(), bd.getSide(),
 
            // type
            bd.getType(), bd.getType(),
 
            // row_data
            bd.getData(), bd.getData(),
 
            // row sense
            bd.getSense(), 
 
            // col sense
            bd.getSense()
 
        );
      }
 
      break;
    case OP::OPROW:
      SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
              "Op LHS is done by OP space");
      break;
    default:
      CHK_MOAB_THROW(MOFEM_IMPOSSIBLE_CASE,
                     (std::string("wrong op type ") +
                      OpBaseDerivativesBase::OpTypeNames[OP::opType])
                         .c_str());
    }
 
    MoFEMFunctionReturn(0);
  };
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, j);
    FTENSOR_INDEX(SPACE_DIM, J);
 
    auto nb_gauss_pts = getGaussPts().size2();
    int nb_row_dofs = row_data.getIndices().size();
    int nb_col_dofs = col_data.getIndices().size();
    OP::locMat.resize(nb_row_dofs, nb_col_dofs, false);
    OP::locMat.clear();
 
    auto &stiffens_map = releaseDataPtr->getInvStiffnessAtPts(
        dataAtPts, normalsAtGaussPts, sense);
    t_inv_D = getFTensor2FromMat<3, 3>(stiffens_map);
 
    auto t_w = OP::getFTensor0IntegrationWeight();
    auto t_row_base = row_data.getFTensor1N<3>();
    for (auto gp = 0; gp < nb_gauss_pts; ++gp) {
      double alpha = t_w / 2.;
      auto nrm2 = t_normal.l2();
      int rr = 0;
      for (; rr != nb_dofs / SPACE_DIM; ++rr) {
        auto t_row = t_row_base(J) * (rowSense * t_normal(J));
        auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
        auto t_mat = getFTensor2FromMat<3, 3, 3>(OP::locMat(rr * 3, 0));
        for (int cc = 0; cc != nb_dofs / SPACE_DIM; ++cc) {
          auto t_col = t_col_base(J) * ((colSense * t_normal(J)) / nrm2);
          t_mat(i, j) += (alpha * (t_row * t_col)) * t_inv_D(i, j);
          ++t_col_base;
          ++t_mat;
        }
        ++t_row_base;
      }
      ++t_w;
      ++t_inv_D;
    }
 
    MoFEMFunctionReturn(0);
  }
 
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  boost::shared_ptr<CrackInterface> releaseDataPtr;
  int rowSense;
  int colSense;
};
 
} // namespace EshelbianPlasticity
 
#endif // ESHELBIANINTERFACE_HPP