BasicFeTools.hpp

Go to the documentation of this file.

/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
// #include <BasicFiniteElements.hpp>
// namespace ContactPlasticity {
#pragma once
 
struct CacheParams {
 
  double young_modulus;
  double young_modulus_inv;
  double poisson;
  double density;
  // for plastic
  double sigmaY;
  double cn_pl;
  // hardening
  double H;
  double C1_k;
  double Q_inf;
  double b_iso;
  double scale_constraint;
 
  Ddg<double, 3, 3> tD;
  Ddg<double, 3, 3> scaling_mat; // for flow rule
 
  // for contact
  double cn_cont;
  double rollers_stop_time;
 
  // for rotation
  static double gravity_y;
  static double spring_stiffness;
 
  static int closest_roller_id;
  static Tensor2<double, 3, 3> Omega;
  static Tensor4<double, 3, 3, 3, 3> Is; // symmetric
 
  static boost::ptr_vector<RigidBodyData> rollerDataVec;
  static RigidBodyData *closest_roller;
 
  CacheParams()
      : young_modulus(10), young_modulus_inv(1), poisson(0.), density(1),
        sigmaY(1e9), cn_pl(1), H(0), C1_k(0), Q_inf(0), b_iso(0),
        cn_cont(1. / young_modulus), /*spring_stiffness(0),
                                        closest_roller_id(0), gravity_y(0),*/
        rollers_stop_time(INFINITY), scale_constraint(1) {}
 
  virtual MoFEMErrorCode scaleParameters() = 0;
 
  template <typename T1, typename T2>
  int getClosestRollerID(Tensor1<T1, 3> &t_coords, Tensor1<T2, 3> &t_disp) {
 
    std::vector<double> gaps;
    Index<'i', 3> i;
    for (auto &roller : rollerDataVec) {
      auto t_normal = roller.getNormal(t_coords, t_disp);
      double dist = roller.getGap(t_coords);
      gaps.push_back(dist);
    }
 
    auto it = std::min_element(gaps.begin(), gaps.end());
    auto nb = distance(gaps.begin(), it);
 
    return nb;
  };
};
 
static Index<'i', 3> i;
static Index<'j', 3> j;
static Index<'k', 3> k;
static Index<'l', 3> l;
static Index<'m', 3> m;
static Index<'n', 3> n;
static Index<'o', 3> o;
static Index<'p', 3> p;
 
static Index<'I', 3> I;
static Index<'J', 3> J;
static Index<'K', 3> K;
static Index<'L', 3> L;
static Index<'M', 3> M;
static Index<'N', 3> N;
 
static Number<0> N0;
static Number<1> N1;
static Number<2> N2;
 
struct BlockParamData : public CacheParams {
 
  MoFEMErrorCode getOptionsFromCommandLine() {
    MoFEMFunctionBegin;
 
    // CHKERR PetscOptionsInsertString(NULL, default_options);
    CHKERR PetscOptionsInsertString(
        NULL,
        "-mat_mumps_icntl_14 800 -mat_mumps_icntl_24 1 -mat_mumps_icntl_13 1 ");
 
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
 
    CHKERR PetscOptionsScalar("-gravity_y", "Gravity load in Y direction", "",
                              gravity_y, &gravity_y, PETSC_NULL);
 
    // for elastic
    CHKERR PetscOptionsScalar("-young", "Young's modulus", "", young_modulus,
                              &young_modulus, PETSC_NULL);
    CHKERR PetscOptionsScalar("-poisson", "Poisson ratio", "", poisson,
                              &poisson, PETSC_NULL);
    CHKERR PetscOptionsScalar("-density", "Density", "", density, &density,
                              PETSC_NULL);
 
    double rot_vec[3] = {0, 0, 1};
    int nb_ent = 3;
    PetscBool is_rotating;
    // get rotation velocity vector
    CHKERR PetscOptionsGetRealArray(NULL, NULL, "-rot_omega", rot_vec, &nb_ent,
                                    &is_rotating);
    if (nb_ent < 3 && is_rotating) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "provide an appropriate number of entries for omega vector (3)");
    }
    {
      Tensor1<double, 3> t1_omega(rot_vec[0], rot_vec[1], rot_vec[2]);
      Omega(i, k) = levi_civita<double>(i, j, k) * t1_omega(j);
    }
 
    if (poisson >= 0.5)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "poisson ratio value not supported");
 
    CHKERR PetscOptionsScalar("-sigmaY", "SigmaY - limit stress", "", sigmaY,
                              &sigmaY, PETSC_NULL);
    CHKERR PetscOptionsScalar("-H", "Hardening modulus", "", H, &H, PETSC_NULL);
    CHKERR PetscOptionsScalar("-cn_pl",
                              "cn_pl parameter for active set constrains", "",
                              cn_pl, &cn_pl, PETSC_NULL);
 
    CHKERR PetscOptionsScalar("-scale_constraint", "constraint scale", "",
                              scale_constraint, &scale_constraint, PETSC_NULL);
 
    CHKERR PetscOptionsScalar("-C1_k", "1st Backstress", "", C1_k, &C1_k,
                              PETSC_NULL);
    CHKERR PetscOptionsScalar("-Q_inf", "Saturation stress", "", Q_inf, &Q_inf,
                              PETSC_NULL);
    CHKERR PetscOptionsScalar("-b_iso", "Isotropic exponent", "", b_iso, &b_iso,
                              PETSC_NULL);
    // for contact
    CHKERR PetscOptionsScalar("-spring",
                              "Springs stiffness on the boundary (contact)", "",
                              spring_stiffness, &spring_stiffness, PETSC_NULL);
    cn_cont = 1. / young_modulus;
 
    CHKERR PetscOptionsScalar("-cn_cont",
                              "cn_cont parameter for active set constrains", "",
                              cn_cont, &cn_cont, PETSC_NULL);
    CHKERR PetscOptionsScalar(
        "-rollers_stop", "Time at which rollers will stop moving", "",
        rollers_stop_time, &rollers_stop_time, PETSC_NULL);
 
    constexpr int MAX_NB_OF_ROLLERS = 10;
    // if (with_contact) {
 
    const char *body_list[] = {"plane",        "sphere", "cylinder", "torus",
                               "superquadric", "stl",    "nurbs"};
    for (int ii = 0; ii < MAX_NB_OF_ROLLERS; ++ii) {
      int nb_coord = 3;
      int nb_disp = 3;
      int nb_dirs = 3;
      VectorDouble center_coords({0, 0, 0});
      VectorDouble disp({0, 0, 0});
      VectorDouble direction({0, 0, 1});
      PetscBool flg = PETSC_FALSE;
      PetscBool rflg;
      string param = "-body" + to_string(ii);
      // char mesh_file_name[255];
      // CHKERR PetscOptionsString(param.c_str(), "file name", "", "mesh.vtk",
      //                           mesh_file_name, 255, &flg);
      int body_type_id = LASTBODY;
      CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, param.c_str(), body_list,
                                  LASTBODY, &body_type_id, &flg);
 
      param = "-coords" + to_string(ii);
      CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(),
                                      &center_coords(0), &nb_coord, &rflg);
      param = "-move" + to_string(ii);
      CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(), &disp(0),
                                      &nb_disp, &rflg);
      if (flg) {
 
        switch (body_type_id) {
        case PLANE:
          rollerDataVec.push_back(new PlaneRigidBody(center_coords, disp, ii));
          break;
        case SPHERE:
          rollerDataVec.push_back(new SphereRigidBody(center_coords, disp, ii));
          break;
        case CYLINDER:
          rollerDataVec.push_back(
              new CylinderRigidBody(center_coords, disp, ii));
          break;
        case TORUS:
          rollerDataVec.push_back(new TorusRigidBody(center_coords, disp, ii));
          break;
        case SUPERQUADRIC:
          rollerDataVec.push_back(
              new SuperquadricRigidBody(center_coords, disp, ii));
          break;
        case STL:
          rollerDataVec.push_back(new STLRigidBody(center_coords, disp, ii));
          break;
        case NURBS:
          // code to be executed if
          break;
        default:
          //  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
          //     "this type of rigid body is not supported, check spelling "
          //     "(plane, sphere, cylinder, torus, superquadric, stl, nurbs)");
          break;
        }
        CHKERR rollerDataVec.back().getBodyOptions();
      }
    }
 
    // if (rollerDataVec.empty()) {
    //   SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
    //           "provide at least one roller data for contact (-roller "
    //           "x,y    -radius r and opional -move x,y) ");
    // }
    // }
 
    ierr = PetscOptionsEnd();
    CHKERRQ(ierr);
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode scaleParameters() {
    MoFEMFunctionBeginHot;
 
    young_modulus_inv = 1. / young_modulus;
    young_modulus = 1;
    sigmaY *= young_modulus_inv;
    H *= young_modulus_inv;
    Q_inf *= young_modulus_inv;
    C1_k *= young_modulus_inv;
 
    MoFEMFunctionReturnHot(0);
  }
};
 
extern BlockParamData *cache;
extern std::map<int, BlockParamData> mat_blocks; // for use in future
 
struct DataFromMove {
  size_t comp;
  double values;
  VectorDouble scaledValues;
  Range ents;
  DataFromMove(size_t c, double val, Range ents_bc)
      : comp(c), values(val), scaledValues(1), ents(ents_bc) {
    scaledValues(0) = val;
  }
  DataFromMove() = delete;
};
 
struct MpSnesCtx : public SnesCtx {
  SmartPetscObj<Vec> mDiagVec;
  MpSnesCtx(MoFEM::Interface &m_field, const std::string &problem_name,
            SmartPetscObj<Vec> m_diag_vec)
      : SnesCtx(m_field, problem_name), mDiagVec(m_diag_vec) {}
};
 
/**
 * \brief Not used at this stage. Could be used to do some calculations,
 * before assembly of local elements.
 */
struct FePrePostProcess : public FEMethod {
 
  MoFEM::Interface &mField;
  boost::ptr_vector<RigidBodyData> &rollerDataVec;
  boost::ptr_vector<DataFromMove> &bcData;
  boost::ptr_vector<MethodForForceScaling> methodsOp;
  bool scaleMat;
  bool runOnceFlag;
  FePrePostProcess(MoFEM::Interface &m_field,
                   boost::ptr_vector<RigidBodyData> &roller_data,
                   boost::ptr_vector<DataFromMove> &bc_data, bool scale_mat)
      : mField(m_field), rollerDataVec(roller_data), bcData(bc_data),
        scaleMat(scale_mat) {}
 
  MoFEMErrorCode preProcess() {
    //
    MoFEMFunctionBegin;
    runOnceFlag = true;
    CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
 
 
    switch (ts_ctx) {
    case CTX_TSSETIJACOBIAN: {
      snes_ctx = CTX_SNESSETJACOBIAN;
      snes_B = ts_B;
      // CHKERR MatSetOption(ts_B, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
      // CHKERR MatSetOption(ts_B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
      break;
    }
    case CTX_TSSETIFUNCTION: {
      snes_ctx = CTX_SNESSETFUNCTION;
      snes_f = ts_F;
      break;
    }
    default:
      break;
    }
 
    auto time = ts_t;
 
    // forceDataScaled = forceData;
    // pressureDataScaled = pressureData;
 
    for (auto &bdata : bcData) {
      bdata.scaledValues(0) = bdata.values;
      CHKERR MethodForForceScaling::applyScale(this, methodsOp,
                                               bdata.scaledValues);
    }
 
    if (time <= (*cache).rollers_stop_time) {
      for (auto &roller : (*cache).rollerDataVec) {
        roller.BodyDispScaled = roller.rollerDisp;
        CHKERR MethodForForceScaling::applyScale(this, methodsOp,
                                                 roller.BodyDispScaled);
      }
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode postProcess() {
 
    MoFEMFunctionBegin;
    if (runOnceFlag) {
      runOnceFlag = false;
    } else
      MoFEMFunctionReturnHot(0);
 
    if (!scaleMat)
      MoFEMFunctionReturnHot(0);
 
    Simple *simple = mField.getInterface<Simple>();
    const Problem *problem_ptr;
    auto dm = simple->getDM();
    CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
 
    SnesCtx *snes_ctx;
    CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
    MpSnesCtx *mps_ctx = (MpSnesCtx *)snes_ctx;
 
    auto diag_m_vec = mps_ctx->mDiagVec;
 
    vector<string> fields;
    const Field_multiIndex *fields_ptr;
    CHKERR mField.get_fields(&fields_ptr);
    for (auto const &field : (*fields_ptr)) {
      const std::string field_name = field->getName();
      fields.push_back(field_name);
    }
 
    // if (fields.size() < 2)
    //   MoFEMFunctionReturnHot(0);
 
    auto calculate_diagonal_scale = [&]() {
      MoFEMFunctionBegin;
 
      const double *f_array;
      CHKERR VecGetArrayRead(ts_F, &f_array);
      auto &rows_index =
          problem_ptr->getNumeredRowDofsPtr()->get<Unique_mi_tag>();
      const int nb_local = problem_ptr->getNbLocalDofsRow();
 
      size_t nb_fields = fields.size();
      vector<double> lnorms(nb_fields, 0);
      vector<double> norms(nb_fields, 0);
 
      for (size_t f = 0; f != nb_fields; ++f) {
        const auto bit_number = mField.get_field_bit_number(fields[f]);
        const auto lo_uid = FieldEntity::getLoBitNumberUId(bit_number);
        const auto hi_uid = FieldEntity::getHiBitNumberUId(bit_number);
        for (auto lo = rows_index.lower_bound(lo_uid);
             lo != rows_index.upper_bound(hi_uid); ++lo) {
          const int local_idx = (*lo)->getPetscLocalDofIdx();
          if (local_idx >= 0 && local_idx < nb_local) {
            const double val = f_array[local_idx];
            lnorms[f] += val * val;
          }
        }
      }
 
      CHKERR VecRestoreArrayRead(ts_F, &f_array);
      CHKERR MPIU_Allreduce(lnorms.data(), norms.data(), nb_fields, MPIU_REAL,
                            MPIU_SUM, PetscObjectComm((PetscObject)dm));
 
      for (auto &v : norms)
        v = sqrt(v);
 
      // for (size_t f = 0; f != nb_fields; ++f)
      //   MOFEM_LOG_C("WORLD", Sev::verbose,
      //               "Scaling diagonal for field[ %s ] by %9.8e",
      //               fields[f].c_str(), norms[f]);
 
      CHKERR VecSet(diag_m_vec, 1.);
 
      double *diag_m_array;
      CHKERR VecGetArray(diag_m_vec, &diag_m_array);
 
      for (size_t f = 0; f != nb_fields; ++f) {
        const auto bit_number = mField.get_field_bit_number(fields[f]);
        const auto lo_uid = FieldEntity::getLoBitNumberUId(bit_number);
        const auto hi_uid = FieldEntity::getHiBitNumberUId(bit_number);
        for (auto lo = rows_index.lower_bound(lo_uid);
             lo != rows_index.upper_bound(hi_uid); ++lo) {
          const int local_idx = (*lo)->getPetscLocalDofIdx();
          if (local_idx >= 0 && local_idx < nb_local)
            diag_m_array[local_idx] = 1. / norms[f];
        }
      }
 
      CHKERR VecRestoreArray(diag_m_vec, &diag_m_array);
 
      switch (ts_ctx) {
      case CTX_TSSETIJACOBIAN: {
        CHKERR MatAssemblyBegin(ts_B, MAT_FINAL_ASSEMBLY);
        CHKERR MatAssemblyEnd(ts_B, MAT_FINAL_ASSEMBLY);
        CHKERR MatDiagonalScale(ts_B, diag_m_vec, PETSC_NULL);
        break;
      }
      case CTX_TSSETIFUNCTION: {
        CHKERR VecPointwiseMult(ts_F, ts_F, diag_m_vec);
      }
      default:
        break;
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR calculate_diagonal_scale();
 
    MoFEMFunctionReturn(0);
  }
};
 
static MoFEMErrorCode custom_snes_rhs(SNES snes, Vec x, Vec f, void *ctx) {
  MpSnesCtx *snes_ctx = (MpSnesCtx *)ctx;
  MoFEMFunctionBegin;
  CHKERR SnesRhs(snes, x, f, ctx);
  CHKERR VecPointwiseMult(f, f, snes_ctx->mDiagVec);
  MoFEMFunctionReturn(0);
}
 
static MoFEMErrorCode custom_snes_mat(SNES snes, Vec x, Mat A, Mat B,
                                      void *ctx) {
  MpSnesCtx *snes_ctx = (MpSnesCtx *)ctx;
  MoFEMFunctionBegin;
  CHKERR SnesMat(snes, x, A, B, ctx);
  CHKERR MatDiagonalScale(B, snes_ctx->mDiagVec, PETSC_NULL);
 
  MoFEMFunctionReturn(0);
};
 
struct LoadScale : public MethodForForceScaling {
  double &scale;
  LoadScale(double &my_scale) : scale(my_scale){};
  MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &nf) {
    MoFEMFunctionBegin;
    nf *= scale;
    MoFEMFunctionReturn(0);
  }
};
 
struct EraseRows : public FEMethod {
  EraseRows(MoFEM::Interface &m_field, std::string problem_name,
            Range &essential_boundary_entsX, Range &essential_boundary_entsY,
            Range &essential_boundary_entsZ)
      : mField(m_field), entX(essential_boundary_entsX),
        entY(essential_boundary_entsY), entZ(essential_boundary_entsZ),
        problemName(problem_name) {}
 
  MoFEMErrorCode preProcess() {
    eraseFlag = true;
    return 0;
  }
 
  MoFEMErrorCode operator()() {
    MoFEMFunctionBegin;
 
    auto pre_proc_lhs = [&](Range &boundary_ents, int comp) {
      MoFEMFunctionBegin;
 
      IS bc_is;
 
      CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
          problemName, ROW, "U", comp, comp, &bc_is, &boundary_ents);
 
      // CHKERR ISView(bc_is, PETSC_VIEWER_STDOUT_SELF);
 
      CHKERR MatZeroRowsIS(ts_B, bc_is, 0.0, PETSC_NULL, PETSC_NULL);
 
      CHKERR ISDestroy(&bc_is);
      MoFEMFunctionReturn(0);
    };
 
    auto pre_proc_rhs = [&](Range &boundary_ents, int comp) {
      MoFEMFunctionBegin;
 
      IS bc_is;
 
      CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
          problemName, ROW, "U", comp, comp, &bc_is, &boundary_ents);
 
      // CHKERR ISView(bc_is, PETSC_VIEWER_STDOUT_SELF);
      // CHKERR VecView(ts_F,PETSC_VIEWER_STDOUT_WORLD);
 
      CHKERR VecISSet(ts_F, bc_is, 0.0);
 
      CHKERR ISDestroy(&bc_is);
      MoFEMFunctionReturn(0);
    };
 
    if (eraseFlag) {
 
      switch (ts_ctx) {
      case CTX_TSSETIJACOBIAN: {
        CHKERR MatSetOption(ts_B, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
        CHKERR MatSetOption(ts_B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
 
        CHKERR MatAssemblyBegin(ts_B, MAT_FINAL_ASSEMBLY);
        CHKERR MatAssemblyEnd(ts_B, MAT_FINAL_ASSEMBLY);
        CHKERR pre_proc_lhs(entX, 0);
        CHKERR pre_proc_lhs(entY, 1);
        CHKERR pre_proc_lhs(entZ, 2);
 
        break;
      }
      case CTX_TSSETIFUNCTION: {
        CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
        CHKERR VecAssemblyBegin(ts_F);
        CHKERR VecAssemblyEnd(ts_F);
        CHKERR pre_proc_rhs(entX, 0);
        CHKERR pre_proc_rhs(entY, 1);
        CHKERR pre_proc_rhs(entZ, 2);
 
        CHKERR VecAssemblyBegin(ts_F);
        CHKERR VecAssemblyEnd(ts_F);
 
        CHKERR VecGhostUpdateBegin(ts_F, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecGhostUpdateEnd(ts_F, ADD_VALUES, SCATTER_REVERSE);
 
        break;
      }
      default:
        break;
      }
 
      eraseFlag = false;
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode postProcess() { return 0; }
 
private:
  Range entX;
  Range entY;
  Range entZ;
  bool eraseFlag;
  MoFEM::Interface &mField;
  std::string problemName;
};
 
// }; // namespace OpContactTools