DispMap.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 __DISP_MAP_HPP__
#define __DISP_MAP_HPP__
 
namespace CellEngineering {
 
struct DispMapFe : public MoFEM::FaceElementForcesAndSourcesCore {
  DispMapFe(MoFEM::Interface &mField)
      : MoFEM::FaceElementForcesAndSourcesCore(mField) {}
  /** \brief Set integrate rule
   */
  int getRule(int order) { return 2 * order + 1; };
};
 
struct DataFromFiles {
 
  string &fIle;
 
  double sCale;
  double cUbe_size;
  vector<vector<double>> main_vector;
  double lAmbda;
 
  DataFromFiles(string &file) : fIle(file) {}
 
  MoFEMErrorCode loadFileData() {
    MoFEMFunctionBegin;
    ifstream in(fIle.c_str());
    if (!in.is_open()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "file with data points not found");
      //  cout << ">> couldn't read the file <<"<< endl;
    }
 
    typedef boost::tokenizer<boost::escaped_list_separator<char>> Tokenizer;
 
    string line;
 
    main_vector.clear();
    while (getline(in, line)) {
      vector<double> vec;
      Tokenizer tok(line);
      for (Tokenizer::iterator it(tok.begin()), end(tok.end()); it != end;
           ++it) {
        vec.push_back(atof(it->c_str()));
      }
 
      main_vector.push_back(vec);
    }
    reverse(main_vector.begin(), main_vector.end());
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode fromOptions() {
    MoFEMFunctionBegin;
    PetscOptionsBegin(PETSC_COMM_WORLD, "", "Get parameters", "none");
 
    cUbe_size = 4;
    CHKERR PetscOptionsScalar("-cube_size", "get cube size for the averaging",
                              "", 0, &cUbe_size, PETSC_NULLPTR);
 
    sCale = 1;
    CHKERR PetscOptionsScalar(
        "-scale", "scale the distance between voxels (eg. from mm to m)", "",
        sCale, &sCale, PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    MoFEMFunctionReturn(0);
  }
 
  /**
   * returns vertors of data for given index of the point
   * @param  vec_idx reference vector of indices
   * @param  vec_data  reference vector of data (colors)
   * @return         error code
   */
 
  double getDataForGivenCoordinate(const double x, const double y,
                                   vector<vector<double>> &main_vector,
                                   const double sCale, const double cUbe_size) {
 
    const int size_y = main_vector.size(); // flip
    const int size_x = main_vector.begin()->size();
    double data = 0;
 
    vector<int> vec_ix;
    vector<int> vec_iy;
    vector<double> vec_data;
    vector<double> vec_dist;
 
    int nx = ceil(cUbe_size);
    int ny = ceil(cUbe_size);
 
    int ix = ceil(x / sCale);
    int iy = ceil(y / sCale);
 
    vec_ix.resize(2 * nx);
    vec_iy.resize(2 * ny);
 
    int ii = 0;
    for (int i = 0; i < 2 * nx; i++) {
      int idx = ix - nx + i;
      if (idx >= size_x / 2 || idx < -size_x / 2)
        continue;
      vec_ix[ii++] = idx;
    }
    vec_ix.resize(ii);
    ii = 0;
    for (int i = 0; i < 2 * ny; i++) {
      int idx = iy - ny + i;
      if (idx >= size_y / 2 || idx < -size_y / 2)
        continue;
      vec_iy[ii++] = idx;
    }
    vec_iy.resize(ii);
    ii = 0;
 
    vec_data.resize(vec_iy.size() * vec_ix.size());
    vec_dist.resize(vec_iy.size() * vec_ix.size());
 
    for (vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();
         it_iy++) {
      for (vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();
           it_ix++) {
        vec_data[ii] =
            main_vector[*it_iy + size_y / 2][*it_ix + size_x / 2]; // flip
        vec_dist[ii] = ((*it_iy * sCale - y) * (*it_iy * sCale - y) +
                        (*it_ix * sCale - x) * (*it_ix * sCale - x));
 
        ii++;
      }
    }
 
    // gaussian smoothing
    vector<double> kernel;
    kernel.resize(vec_data.size());
    int i = 0;
    double sigma = 10;
    double sum = 0;
    const double m = (1 / sqrt(M_PI * 2 * sigma * sigma));
    for (int ii = 0; ii < vec_dist.size(); ii++) {
      kernel[i] = m * exp(-(vec_dist[ii]) /
                          (2 * sigma * sigma)); // distance is already squared
      sum += kernel[i++];
    }
    ii = 0;
    for (vector<double>::iterator vec_itr = vec_data.begin();
         vec_itr != vec_data.end(); vec_itr++) {
      kernel[ii] /= sum;
      data += (*vec_itr) * kernel[ii++];
    }
    return data;
  }
 
};
 
struct CommonData {
 
  VectorDouble dispAtGaussPts; ///< Values at integration Pts
  VectorDouble locF;           ///< Local element rhs vector
  MatrixDouble locA;           ///< Local element matrix
  MatrixDouble transLocA;
  Mat globA; ///< Global matrix
  Vec globF;
  double lAmbda;
 
  MoFEMErrorCode getParameters() {
 
    MoFEMFunctionBegin;
    PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
 
    this->lAmbda = 1;
    CHKERR PetscOptionsScalar("-lambda", "lambda parameter", "", 1, &lAmbda,
                              PETSC_NULLPTR);
    PetscOptionsEnd();
 
    MoFEMFunctionReturn(0);
  }
};
 
struct OpCalculateLhs
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  CommonData &commonData;
  DataFromFiles &dataFromFiles;
 
  OpCalculateLhs(const std::string &row_field, const std::string &col_field,
                 CommonData &common_data, DataFromFiles &data_from_files)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            row_field, col_field, UserDataOperator::OPROWCOL),
        commonData(common_data), dataFromFiles(data_from_files) {
    sYmm = true;
  }
 
  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_rows = row_data.getIndices().size(); // number of rows
    int nb_cols = col_data.getIndices().size(); // number of columns
    if (nb_rows == 0)
      MoFEMFunctionReturnHot(0);
    if (nb_cols == 0)
      MoFEMFunctionReturnHot(0);
 
    commonData.locA.resize(nb_rows, nb_cols, false);
    commonData.locA.clear();
    commonData.getParameters();
    int nb_gauss_pts = row_data.getN().size1();
    for (int gg = 0; gg < nb_gauss_pts; gg++) {
      double area = getArea() * getGaussPts()(2, gg);
      // if(getHoGaussPtsDetJac().size()>0) {
      //   area *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
      // }
      VectorAdaptor row_base_function = row_data.getN(gg);
      VectorAdaptor col_base_function = col_data.getN(gg);
      commonData.locA +=
          area * outer_prod(row_base_function, col_base_function);
      MatrixAdaptor row_diff_base_function = row_data.getDiffN(gg);
      MatrixAdaptor col_diff_base_function = col_data.getDiffN(gg);
 
      // const double lambda = 1;
 
      commonData.locA +=
          area * commonData.lAmbda * dataFromFiles.sCale * dataFromFiles.sCale *
          prod(row_diff_base_function, trans(col_diff_base_function));
    }
    CHKERR MatSetValues(commonData.globA, row_data.getIndices().size(),
                        &row_data.getIndices()[0], col_data.getIndices().size(),
                        &col_data.getIndices()[0], &commonData.locA(0, 0),
                        ADD_VALUES);
 
    // Assemble off symmetric side
    if (row_side != col_side || row_type != col_type) {
      commonData.transLocA.resize(commonData.locA.size2(),
                                  commonData.locA.size1(), false);
      noalias(commonData.transLocA) = trans(commonData.locA);
      CHKERR MatSetValues(
          commonData.globA, col_data.getIndices().size(),
          &col_data.getIndices()[0], row_data.getIndices().size(),
          &row_data.getIndices()[0], &commonData.transLocA(0, 0), ADD_VALUES);
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
struct OpCalulatefRhoAtGaussPts
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  CommonData &commonData;
  DataFromFiles &dataFromFiles;
 
  OpCalulatefRhoAtGaussPts(const string &row_field, CommonData &common_data,
                           DataFromFiles &data_from_files)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            row_field, row_field, UserDataOperator::OPROW),
        commonData(common_data), dataFromFiles(data_from_files) {}
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        DataForcesAndSourcesCore::EntData &row_data) {
    MoFEMFunctionBegin;
 
    if (row_type != MBVERTEX)
      MoFEMFunctionReturnHot(0);
 
    int nb_rows = row_data.getIndices().size();
    if (nb_rows == 0)
      MoFEMFunctionReturnHot(0);
 
    FTensor::Index<'i', 3> i;
 
    FTensor::Tensor1<double *, 3> t1_xyz(&getCoordsAtGaussPts()(0, 0),
                                         &getCoordsAtGaussPts()(0, 1),
                                         &getCoordsAtGaussPts()(0, 2), 3);
    const int nb_gauss_pts = row_data.getN().size1();
    commonData.dispAtGaussPts.resize(nb_gauss_pts, false);
    for (int gg = 0; gg < nb_gauss_pts; gg++) {
 
      commonData.dispAtGaussPts[gg] = dataFromFiles.getDataForGivenCoordinate(
          t1_xyz(0), t1_xyz(1), dataFromFiles.main_vector, dataFromFiles.sCale,
          dataFromFiles.cUbe_size);
      ++t1_xyz;
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
struct OpAssembleRhs
    : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
 
  CommonData &commonData;
 
  OpAssembleRhs(const string &row_field, CommonData &common_data)
      : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
            row_field, row_field, UserDataOperator::OPROW),
        commonData(common_data) {}
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        DataForcesAndSourcesCore::EntData &row_data) {
    MoFEMFunctionBegin;
 
    int nb_rows = row_data.getIndices().size();
    if (nb_rows == 0)
      MoFEMFunctionReturnHot(0);
    commonData.locF.resize(nb_rows, false);
    commonData.locF.clear();
 
    int nb_gauss_pts = row_data.getN().size1();
    for (int gg = 0; gg < nb_gauss_pts; gg++) {
 
      double area = getArea() * getGaussPts()(2, gg);
      VectorAdaptor base_function = row_data.getN(gg);
      commonData.locF += area * commonData.dispAtGaussPts[gg] * base_function;
    }
 
    CHKERR VecSetValues(commonData.globF, row_data.getIndices().size(),
                        &row_data.getIndices()[0], &commonData.locF[0],
                        ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  }
};
 
} // namespace CellEngineering
 
#endif //__DISP_MAP_HPP__