CellForces.hpp

Go to the documentation of this file.

/* Copyright (c) 2025 Authors under the MIT License.
 * See LICENSE.md for details.
 * SPDX-License-Identifier: MIT
 */
 
#ifndef __CELL_FORCES_HPP__
#define __CELL_FORCES_HPP__
 
namespace CellEngineering {
 
struct FatPrism : public MoFEM::FatPrismElementForcesAndSourcesCore {
  FatPrism(MoFEM::Interface &m_field)
      : MoFEM::FatPrismElementForcesAndSourcesCore(m_field) {}
  int getRuleTrianglesOnly(int order) { return 2 * (order - 1); }
  int getRuleThroughThickness(int order) { return 2 * order; }
};
 
/**
 * Common data structure used to pass data between operators
 */
struct CommonData {
 
  MatrixDouble forceAtGaussPts;
  VectorDouble forcePotentialAtGaussPoints;
  MatrixDouble gradForceAtGaussPtrs;
 
  boost::shared_ptr<VectorDouble> dispX;
  boost::shared_ptr<VectorDouble> dispY;
 
  CommonData()
      : dispX(boost::shared_ptr<VectorDouble>(new VectorDouble())),
        dispY(boost::shared_ptr<VectorDouble>(new VectorDouble())) {}
};
 
/**
 * Face element used to run operators
 */
struct FaceElement : public MoFEM::FaceElementForcesAndSourcesCore {
  FaceElement(MoFEM::Interface &m_field)
      : MoFEM::FaceElementForcesAndSourcesCore(m_field) {}
  int getRule(int order) { return 2 * order; };
};
 
/**
 * Operator to get measured displacements in X-direction on the face
 */
struct OpGetDispX : public MoFEM::OpCalculateScalarFieldValues {
  OpGetDispX(CommonData &common_data)
      : MoFEM::OpCalculateScalarFieldValues("DISP_X", common_data.dispX) {}
};
 
/**
 * Operator to get measured displacements in Y-direction on the face
 */
struct OpGetDispY : public MoFEM::OpCalculateScalarFieldValues {
  OpGetDispY(CommonData &common_data)
      : MoFEM::OpCalculateScalarFieldValues("DISP_Y", common_data.dispY) {}
};
 
/**
 * \brief Calculate and assemble D matrix
 */
struct OpCellPotentialD
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  Mat A;
  const double epsRho;
 
  OpCellPotentialD(Mat a, const double rho)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            "RHO", "RHO", UserDataOperator::OPROWCOL),
        A(a), epsRho(rho) {
    sYmm = false;
  }
 
  MatrixDouble C;
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data) {
 
    MoFEMFunctionBegin;
 
    int nb_row = row_data.getIndices().size();
    int nb_col = col_data.getIndices().size();
 
    if (nb_row == 0)
      MoFEMFunctionReturnHot(0);
    if (nb_col == 0)
      MoFEMFunctionReturnHot(0);
 
    const VectorInt &row_indices = row_data.getIndices();
    const VectorInt &col_indices = col_data.getIndices();
 
    const int nb_gauss_pts = row_data.getN().size1();
 
    C.resize(nb_row, nb_col, false);
    C.clear();
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      // get integration weights
      const double val = getGaussPts()(2, gg) * getArea() * epsRho;
      for (int rr = 0; rr != nb_row; rr++) {
        const double diff_nr_x = row_data.getDiffN()(gg, 2 * rr + 0);
        const double diff_nr_y = row_data.getDiffN()(gg, 2 * rr + 1);
        for (int cc = 0; cc != nb_col; cc++) {
          const double diff_nc_x = col_data.getDiffN()(gg, 2 * cc + 0);
          const double diff_nc_y = col_data.getDiffN()(gg, 2 * cc + 1);
          C(rr, cc) += val * (diff_nr_x * diff_nc_x + diff_nr_y * diff_nc_y);
        }
      }
    }
 
    CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
                        &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  }
};
 
/**
 * \brief Calculate and assemble B matrix
 
 \f[
 \mathbf{B} = \int_{S_0}
 \mathbf{N}^\textrm{T}_\textrm{u} \frac{\partial \mathbf{N}_\rho}{\partial
 \mathbf{x}} \textrm{d}S \f]
 
 */
struct OpCellPotentialB
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  Mat A;
 
  OpCellPotentialB(Mat a, const string field_name)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            "RHO", field_name, UserDataOperator::OPROWCOL),
        A(a) {
    sYmm = false;
  }
 
  MatrixDouble C;
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data) {
 
    MoFEMFunctionBegin;
 
    int nb_row = row_data.getIndices().size();
    int nb_col = col_data.getIndices().size();
 
    if (nb_row == 0)
      MoFEMFunctionReturnHot(0);
    if (nb_col == 0)
      MoFEMFunctionReturnHot(0);
 
    const VectorInt &row_indices = row_data.getIndices();
    const VectorInt &col_indices = col_data.getIndices();
 
    const int nb_gauss_pts = row_data.getN().size1();
 
    C.resize(nb_row, nb_col, false);
    C.clear();
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      // get integration weights
      const double val = getGaussPts()(2, gg) * getArea();
      for (int rr = 0; rr != nb_row; rr++) {
        const double diff_nr_x = row_data.getDiffN()(gg, 2 * rr + 0);
        const double diff_nr_y = row_data.getDiffN()(gg, 2 * rr + 1);
        for (int cc = 0; cc != nb_col / 3; cc++) {
          const double nc = col_data.getN(gg)[cc];
          C(rr, 3 * cc + 0) -= val * diff_nr_x * nc;
          C(rr, 3 * cc + 1) -= val * diff_nr_y * nc;
        }
      }
    }
 
    CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
                        &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  }
};
 
/**
 * \brief Calculate and assemble S matrix
 
 */
struct OpCellS
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  Mat A;
  const double epsU;
 
  OpCellS(Mat a, const double eps_u)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            "UPSILON", "U", UserDataOperator::OPROWCOL),
        A(a), epsU(eps_u) {
    sYmm = false;
  }
 
  MatrixDouble C;
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data) {
 
    MoFEMFunctionBegin;
 
    int nb_row = row_data.getIndices().size();
    int nb_col = col_data.getIndices().size();
 
    if (nb_row == 0)
      MoFEMFunctionReturnHot(0);
    if (nb_col == 0)
      MoFEMFunctionReturnHot(0);
 
    const VectorInt &row_indices = row_data.getIndices();
    const VectorInt &col_indices = col_data.getIndices();
 
    const int nb_gauss_pts = row_data.getN().size1();
 
    C.resize(nb_row, nb_col, false);
    C.clear();
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      // get integration weights
      const double val = getGaussPts()(2, gg) * getArea() * pow(epsU, -1);
      for (int rr = 0; rr != nb_row / 3; rr++) {
        const double nr = row_data.getN(gg)[rr];
        for (int cc = 0; cc != nb_col / 3; cc++) {
          const double nc = col_data.getN(gg)[cc];
          const double c = val * nr * nc;
          C(3 * rr + 0, 3 * cc + 0) += c;
          C(3 * rr + 1, 3 * cc + 1) += c;
        }
      }
    }
 
    CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
                        &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  }
};
 
/**
 * \brief Calculate and assemble g vector
 
 \f[
 \mathbf{g} = \int_{S_1} \mathbf{N}^\textrm{T}_\textrm{u} \overline{\mathbf{u}}
 \textrm{d}S, \f]
 
 */
struct OpCell_g
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  Vec F;
  const double epsU;
  CommonData &commonData;
 
  OpCell_g(Vec f, const double eps_u, CommonData &common_data)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            "UPSILON", UserDataOperator::OPROW),
        F(f), epsU(eps_u), commonData(common_data) {}
 
  VectorDouble nF;
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) {
 
    MoFEMFunctionBegin;
 
    int nb_row = data.getIndices().size();
    if (nb_row == 0)
      MoFEMFunctionReturnHot(0);
 
    const VectorDouble &disp_x = *commonData.dispX;
    const VectorDouble &disp_y = *commonData.dispY;
 
    const int nb_gauss_pts = data.getN().size1();
 
    nF.resize(nb_row, false);
    nF.clear();
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      const double val = getGaussPts()(2, gg) * getArea() * pow(epsU, -1);
      for (int rr = 0; rr != nb_row / 3; rr++) {
        const double nr = val * data.getN(gg)[rr];
        nF[3 * rr + 0] += nr * disp_x[gg];
        nF[3 * rr + 1] += nr * disp_y[gg];
      }
    }
    // assemble local face vector
    CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &nF[0], ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  }
};
 
/**
 * \brief Calculate and assemble Z matrix
 */
struct OpCellCurlD
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  Mat A;
  const double epsRho;
  const double epsL;
 
  OpCellCurlD(Mat a, const double eps_rho, const double eps_l)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            "RHO", "RHO", UserDataOperator::OPROWCOL),
        A(a), epsRho(eps_rho), epsL(eps_l) {
    sYmm = false;
  }
 
  MatrixDouble C;
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data) {
 
    MoFEMFunctionBegin;
 
    int nb_row = row_data.getIndices().size();
    int nb_col = col_data.getIndices().size();
 
    if (nb_row == 0)
      MoFEMFunctionReturnHot(0);
    if (nb_col == 0)
      MoFEMFunctionReturnHot(0);
 
    const VectorInt &row_indices = row_data.getIndices();
    const VectorInt &col_indices = col_data.getIndices();
 
    const int nb_gauss_pts = row_data.getN().size1();
 
    auto t_diff_base_row = row_data.getFTensor2DiffN<3, 2>();
    auto t_base_row = row_data.getFTensor1N<3>();
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 2> j;
    FTensor::Number<0> N0;
    FTensor::Number<1> N1;
 
    C.resize(nb_row, nb_col, false);
    C.clear();
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      // get integration weights
      const double val = getGaussPts()(2, gg) * getArea();
      for (int rr = 0; rr != nb_row; rr++) {
        const double curl_row = t_diff_base_row(1, 0) - t_diff_base_row(0, 1);
        auto t_diff_base_col = col_data.getFTensor2DiffN<3, 2>(gg, 0);
        auto t_base_col = col_data.getFTensor1N<3>(gg, 0);
        for (int cc = 0; cc != nb_col; cc++) {
          const double curl_col = t_diff_base_col(1, 0) - t_diff_base_col(0, 1);
          const double diag = val * epsRho * t_base_row(i) * t_base_col(i);
          const double curl =
              epsL > 0 ? val * pow(epsL, -1) * curl_row * curl_col : 0;
          C(rr, cc) += diag + curl;
          ++t_diff_base_col;
          ++t_base_col;
        }
        ++t_diff_base_row;
        ++t_base_row;
      }
    }
 
    CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
                        &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  }
};
 
/**
 * \brief Calculate and assemble Z matrix
 */
struct OpCellCurlB
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  Mat A;
 
  OpCellCurlB(Mat a, const string field_name)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            "RHO", field_name, UserDataOperator::OPROWCOL),
        A(a) {
    sYmm = false;
  }
 
  MatrixDouble C;
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        DataForcesAndSourcesCore::EntData &row_data,
                        DataForcesAndSourcesCore::EntData &col_data) {
 
    MoFEMFunctionBegin;
 
    int nb_row = row_data.getIndices().size();
    int nb_col = col_data.getIndices().size();
 
    if (nb_row == 0)
      MoFEMFunctionReturnHot(0);
    if (nb_col == 0)
      MoFEMFunctionReturnHot(0);
 
    const VectorInt &row_indices = row_data.getIndices();
    const VectorInt &col_indices = col_data.getIndices();
 
    const int nb_gauss_pts = row_data.getN().size1();
 
    auto t_base_row = row_data.getFTensor1N<3>();
 
    C.resize(nb_row, nb_col, false);
    C.clear();
    for (int gg = 0; gg != nb_gauss_pts; gg++) {
      // get integration weights
      const double val = getGaussPts()(2, gg) * getArea();
      for (int rr = 0; rr != nb_row; rr++) {
        // This is scalar field, diff yield tensor1
        auto t_base_col = col_data.getFTensor0N(gg, 0);
        for (int cc = 0; cc != nb_col / 3; cc++) {
          C(rr, 3 * cc + 0) -= val * t_base_row(0) * t_base_col;
          C(rr, 3 * cc + 1) -= val * t_base_row(1) * t_base_col;
          ++t_base_col;
        }
        ++t_base_row;
      }
    }
 
    CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
                        &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  }
};
 
/**
 * \brief Post-process tractions
 *
 * \todo Should be done with Prisms, no need to have three kinds of elements,
 * tets for solid, prisms for constrains and triangles to do post-processing. We
 * can do post-pocessing with prisms.
 *
 */
struct OpVirtualPotentialRho
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  CommonData &commonData;
 
  OpVirtualPotentialRho(std::string field_name, CommonData &common_data)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            field_name, UserDataOperator::OPROW),
        commonData(common_data) {}
 
  /**
   * Real work is done here, i.e. assemble tractions
   */
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) {
    //
    //
    MoFEMFunctionBegin;
 
    if (data.getFieldData().empty())
      MoFEMFunctionReturnHot(0);
 
    const int nb_gauss_pts = data.getN().size1();
    const int nb_dofs = data.getFieldData().size();
 
    VectorDouble *vec;
    MatrixDouble *mat;
    MatrixDouble *grad;
 
    if (rowFieldName == "RHO") {
      vec = &(commonData.forcePotentialAtGaussPoints);
      mat = &commonData.forceAtGaussPts;
      grad = &commonData.gradForceAtGaussPtrs;
    }
 
    if (type == MBVERTEX) {
      mat->resize(nb_gauss_pts, 2, false);
      mat->clear();
      vec->resize(nb_gauss_pts, false);
      vec->clear();
      grad->resize(0, 0, false);
    }
 
    // calculate forces at integration points
    for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
      for (int rr = 0; rr != nb_dofs; rr++) {
        const double diff_base_x = data.getDiffN()(gg, 2 * rr + 0);
        const double diff_base_y = data.getDiffN()(gg, 2 * rr + 1);
        const double val = data.getFieldData()[rr];
        (*mat)(gg, 0) += diff_base_x * val;
        (*mat)(gg, 1) += diff_base_y * val;
        const double base = data.getN()(gg, rr);
        (*vec)[gg] += base * val;
      }
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
/**
 * \brief Post-process tractions
 *
 * \todo Should be done with Prisms, no need to have three kinds of elements,
 * tets for solid, prisms for constrains and triangles to do post-processing. We
 * can do post-pocessing with prisms.
 *
 */
struct OpVirtualCurlRho
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  CommonData &commonData;
 
  OpVirtualCurlRho(std::string field_name, CommonData &common_data)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            field_name, UserDataOperator::OPROW),
        commonData(common_data) {}
 
  /**
   * Real work is done here, i.e. assemble tractions
   */
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) {
    //
    //
    MoFEMFunctionBegin;
 
    if (data.getFieldData().empty())
      MoFEMFunctionReturnHot(0);
 
    const int nb_gauss_pts = data.getN().size1();
    const int nb_dofs = data.getFieldData().size();
 
    VectorDouble *vec;
    MatrixDouble *mat;
    MatrixDouble *grad;
 
    if (rowFieldName == "RHO") {
      mat = &commonData.forceAtGaussPts;
      grad = &commonData.gradForceAtGaussPtrs;
      vec = &commonData.forcePotentialAtGaussPoints;
    }
 
    if (type == MBEDGE && side == 0) {
      mat->resize(nb_gauss_pts, 2, false);
      mat->clear();
      grad->resize(nb_gauss_pts, 4, false);
      grad->clear();
      vec->resize(0, false);
    }
 
    // calculate forces at integration points
    auto t_base = data.getFTensor1N<3>();
    auto t_diff_base_fun = data.getFTensor2DiffN<3, 2>();
 
    for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
      auto t_data = data.getFTensor0FieldData();
      for (int rr = 0; rr != nb_dofs; rr++) {
        const double base_x = t_base(0);
        const double base_y = t_base(1);
        (*mat)(gg, 0) += base_x * t_data;
        (*mat)(gg, 1) += base_y * t_data;
        for (int i = 0; i != 2; i++) {
          for (int j = 0; j != 2; j++) {
            (*grad)(gg, 2 * i + j) += t_diff_base_fun(i, j) * t_data;
          }
        }
        ++t_data;
        ++t_base;
        ++t_diff_base_fun;
      }
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
/**
 * \brief Shave results on mesh tags for post-processing
 */
struct PostProcTraction
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  MoFEM::Interface &mField;
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  CommonData &commonData;
 
  /**
   * Constructor
   */
  PostProcTraction(MoFEM::Interface &m_field, moab::Interface &post_proc_mesh,
                   std::vector<EntityHandle> &map_gauss_pts,
                   CommonData &common_data)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            "RHO", ForcesAndSourcesCore::UserDataOperator::OPROW),
        mField(m_field), postProcMesh(post_proc_mesh),
        mapGaussPts(map_gauss_pts), commonData(common_data) {}
 
  /**
   * \brief Here real work is done
   */
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) {
    MoFEMFunctionBegin;
 
    if (type != MBVERTEX)
      MoFEMFunctionReturnHot(0);
 
    double def_val[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
    Tag th_force, th_force_potential, th_grad_force;
    CHKERR postProcMesh.tag_get_handle("FORCE", 3, MB_TYPE_DOUBLE, th_force,
                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
    
    if (commonData.forcePotentialAtGaussPoints.size() > 0) {
      CHKERR postProcMesh.tag_get_handle("FORCE_POTENTIAL", 1, MB_TYPE_DOUBLE,
                                         th_force_potential,
                                         MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
    }
    if (commonData.gradForceAtGaussPtrs.size1() > 0) {
      CHKERR postProcMesh.tag_get_handle("FORCE_GRADIENT", 9, MB_TYPE_DOUBLE,
                                         th_grad_force,
                                         MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
    }
 
    int nb_gauss_pts = data.getN().size1();
    if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
    }
 
    for (int gg = 0; gg < nb_gauss_pts; gg++) {
      const double force[] = {commonData.forceAtGaussPts(gg, 0),
                              commonData.forceAtGaussPts(gg, 1), 0};
      CHKERR postProcMesh.tag_set_data(th_force, &mapGaussPts[gg], 1, force);
      if (commonData.forcePotentialAtGaussPoints.size() > 0) {
        CHKERR postProcMesh.tag_set_data(
            th_force_potential, &mapGaussPts[gg], 1,
            &commonData.forcePotentialAtGaussPoints[gg]);
      }
      if (commonData.gradForceAtGaussPtrs.size1() > 0) {
        const double grad_force[] = {commonData.gradForceAtGaussPtrs(gg, 0),
                                     commonData.gradForceAtGaussPtrs(gg, 1),
                                     0,
                                     commonData.gradForceAtGaussPtrs(gg, 2),
                                     commonData.gradForceAtGaussPtrs(gg, 3),
                                     0,
                                     0,
                                     0,
                                     0};
        CHKERR postProcMesh.tag_set_data(th_grad_force, &mapGaussPts[gg], 1,
                                         grad_force);
        
      }
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
} // namespace CellEngineering
 
#endif // __CELL_FORCES_HPP__