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;

  // for contact
  double cn_cont;
  double spring_stiffness;
  double gravity_y;
  double rollers_stop_time;
  int closest_roller_id;

  // for rotation
  Tensor2<double, 3, 3> Omega;
  Tensor4<double, 3, 3, 3, 3> Is; // symmetric
  Tensor4<double, 3, 3, 3, 3> II; // unit

  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), gravity_y(0),
        rollers_stop_time(INFINITY), closest_roller_id(0) {}
};

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  Number<0> N0;
static  Number<1> N1;
static  Number<2> N2;

struct BlockParamData : public CacheParams {

  MoFEMErrorCode getOptionsFromCommandLine() {
    MoFEMFunctionBeginHot;

    // CHKERR PetscOptionsInsertString(NULL, default_options);

    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;
      for (int dd : {0, 1, 2})
        t1_omega(dd) = rot_vec[dd];
      // 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("-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 = 6;
    // if (with_contact) {
    //   for (int i = 0; i < MAX_NB_OF_ROLLERS; ++i) {
    //     int nb_coord = 2;
    //     int nb_disp = 2;
    //     double center_coords[2] = {0, -1};
    //     double disp[2] = {0, 0};
    //     PetscBool flg = PETSC_FALSE;
    //     PetscBool rflg;
    //     string param = "-roller" + to_string(i);
    //     CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(),
    //     center_coords,
    //                                     &nb_coord, &rflg);
    //     if (rflg) {
    //       param = "-move" + to_string(i);
    //       CHKERR PetscOptionsGetRealArray(NULL, NULL, param.c_str(), disp,
    //                                       &nb_disp, &flg);
    //       double radi;
    //       param = "-radius" + to_string(i);
    //       CHKERR PetscOptionsScalar(param.c_str(), "set roller radius", "",
    //       radi,
    //                                 &radi, &flg);
    //       if (nb_coord < 2 || !flg) {
    //         SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
    //                 "provide an appropriate number of coords for roller
    //                 center " "and corresponding radius");
    //       }
    //       VectorDouble center;
    //       VectorDouble displ;
    //       center.data().assign(center_coords, center_coords + 2);
    //       displ.data().assign(disp, disp + 2);
    //       rollerData.push_back(RollerData(center, displ, radi));
    //     }
    //   }

    //   if (rollerData.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);

    MoFEMFunctionReturnHot(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;
// static std::map<int, BlockParamData> cache_blocks; //for use in future

struct DataFromMove {
  size_t comp;
  double values;
  VectorDouble scaled_values;
  Range ents;
  DataFromMove(size_t c, double val, Range ents_bc)
      : comp(c), values(val), scaled_values(1), ents(ents_bc) {
    scaled_values(0) = val;
  }
  DataFromMove() = delete;
};

struct RollerData {
  const VectorDouble cCoords;
  const VectorDouble rollerDisp;
  VectorDouble rollerDispScaled;
  const double rAdius;
  double gAp0;
  Index<'i', 2> i;
  Tensor1<double, 2> t_norm;

  RollerData(VectorDouble c_coords, VectorDouble roller_disp, double radius)
      : cCoords(c_coords), rollerDisp(roller_disp), rAdius(radius) {
    rollerDispScaled = rollerDisp;
    rollerDispScaled.clear();
  }

  inline Tensor1<double, 2> getRollerCenter() {
    VectorDouble pos = cCoords + rollerDispScaled;
    return Tensor1<double, 2>(pos(0), pos(1));
  };

  template <typename T1, typename T2>
  inline Tensor1<double, 2> &getNormal(Tensor1<T1, 3> &t_coords,
                                       Tensor1<T2, 2> &t_disp) {
    auto t_center = getRollerCenter();
    Tensor1<double, 2> t_d;
    t_d(i) = t_coords(i) - t_center(i) + t_disp(i);
    double dist = sqrt(t_d(i) * t_d(i));
    this->t_norm(i) = t_d(i) / dist;

    // calculate the gap0 (from reference)
    t_d(i) = t_d(i) - t_disp(i);
    dist = sqrt(t_d(i) * t_d(i));
    double gap0 = -t_norm(i) * (t_d(i) - t_norm(i) * rAdius);
    this->gAp0 = gap0;
    return t_norm;
  };

  inline double getGap0() {
    // auto norm = getNormal(t_coords);
    return gAp0;
  };
};

template <typename T1, typename T2>
int getClosestRollerID(Tensor1<T1, 3> &t_coords, Tensor1<T2, 2> &t_disp,
                       vector<RollerData> &rollerData) {

  std::vector<double> dists;
  Index<'i', 2> i;

  for (auto &roller : rollerData) {
    // auto t_center = roller.getRollerCenter();
    // Tensor1<double, 2> t_d;
    // t_d(i) = t_coords(i) - t_center(i) + t_disp(i);
    auto t_normal = roller.getNormal(t_coords, t_disp);
    double dist = t_disp(i) * t_normal(i) - roller.getGap0();
    // Tensor1<double, 2> t_normal;
    // t_normal(i) = t_d(i) / dist;
    // dist = -t_normal(i) * (t_d(i) - t_normal(i) * roller.rAdius);
    // dist = abs(dist);
    dists.push_back(dist);
  }
  auto it = std::min_element(dists.begin(), dists.end());
  auto nb = distance(dists.begin(), it);

  return nb;
};

/**
 * \brief Not used at this stage. Could be used to do some calculations,
 * before assembly of local elements.
 */
struct FePrePostProcess : public FEMethod {

  vector<RollerData> &rollerData;
  vector<shared_ptr<DataFromMove>> bcData;
  ;
  boost::ptr_vector<MethodForForceScaling> methodsOp;

  FePrePostProcess(vector<RollerData> &roller_data,
                   vector<shared_ptr<DataFromMove>> bc_data)
      : rollerData(roller_data), bcData(bc_data) {}

  MoFEMErrorCode preProcess() {
    //
    MoFEMFunctionBegin;

    CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
    // CHKERR MatSetOption(ts_B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);

    switch (ts_ctx) {
    case CTX_TSSETIJACOBIAN: {
      snes_ctx = CTX_SNESSETJACOBIAN;
      snes_B = ts_B;
      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->scaled_values(0) = bdata->values;
      CHKERR MethodForForceScaling::applyScale(this, methodsOp,
                                               bdata->scaled_values);
    }

    // if (time <= rollers_stop_time) {
    //   if (with_contact)
    //     for (auto &roller : rollerData) {
    //       roller.rollerDispScaled = roller.rollerDisp;
    //       CHKERR MethodForForceScaling::applyScale(this, methodsOp,
    //                                                roller.rollerDispScaled);
    //     }
    // }
    MoFEMFunctionReturn(0);
  }

  MoFEMErrorCode postProcess() { return 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 EraseLHSRows : public FEMethod {
  EraseLHSRows(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 = [&](Range &boundary_ents, int comp) {
      MoFEMFunctionBegin;
      // if (boundary_ents.empty())
      //   MoFEMFunctionReturnHot(0);
      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);
    };
    if (eraseFlag) {

      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_FLUSH_ASSEMBLY);
      // CHKERR MatAssemblyEnd(ts_B, MAT_FLUSH_ASSEMBLY);

      CHKERR MatAssemblyBegin(ts_B, MAT_FINAL_ASSEMBLY);
      CHKERR MatAssemblyEnd(ts_B, MAT_FINAL_ASSEMBLY);
      CHKERR pre_proc(entX, 0);
      CHKERR pre_proc(entY, 1);
      CHKERR pre_proc(entZ, 2);
      // CHKERR MatAssemblyBegin(ts_B, MAT_FLUSH_ASSEMBLY);
      // CHKERR MatAssemblyEnd(ts_B, MAT_FLUSH_ASSEMBLY);
      // CHKERR MatAssemblyBegin(ts_B, MAT_FINAL_ASSEMBLY);
      // CHKERR MatAssemblyEnd(ts_B, MAT_FINAL_ASSEMBLY);
      eraseFlag = false;
    }
    MoFEMFunctionReturn(0);
  }

  MoFEMErrorCode postProcess() { return 0; }

private:
  Range entX;
  Range entY;
  Range entZ;
  bool eraseFlag;
  MoFEM::Interface &mField;
  std::string problemName;
};

struct EraseRHSRows : public FEMethod {
  EraseRHSRows(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 postProcess() { return 0; }

  MoFEMErrorCode preProcess() {
    eraseFlag = true;
    return 0;
  }

  MoFEMErrorCode operator()() {
    MoFEMFunctionBegin;

    auto pre_proc = [&](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) {
      // CHKERR VecGhostUpdateBegin(ts_F, ADD_VALUES, SCATTER_REVERSE);
      // CHKERR VecGhostUpdateEnd(ts_F, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecAssemblyBegin(ts_F);
      CHKERR VecAssemblyEnd(ts_F);
      CHKERR pre_proc(entX, 0);
      CHKERR pre_proc(entY, 1);
      CHKERR pre_proc(entZ, 2);

      // CHKERR VecGhostUpdateBegin(d0, INSERT_VALUES, SCATTER_FORWARD);
      // CHKERR VecGhostUpdateEnd(d0, INSERT_VALUES, SCATTER_FORWARD);
      // CHKERR DMoFEMPreProcessFiniteElements(sub_dm_rho,
      // density_bc_ptr.get()); CHKERR VecGhostUpdateBegin(d0, ADD_VALUES,
      // SCATTER_REVERSE); CHKERR VecGhostUpdateEnd(d0, ADD_VALUES,
      // SCATTER_REVERSE);

      CHKERR VecAssemblyBegin(ts_F);
      CHKERR VecAssemblyEnd(ts_F);
      CHKERR VecGhostUpdateBegin(ts_F, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(ts_F, ADD_VALUES, SCATTER_REVERSE);

      // CHKERR VecAssemblyBegin(ts_F);
      // CHKERR VecAssemblyEnd(ts_F);
      // CHKERR VecGhostUpdateBegin(ts_F, ADD_VALUES, SCATTER_FORWARD);
      // CHKERR VecGhostUpdateEnd(ts_F, ADD_VALUES, SCATTER_FORWARD);

      eraseFlag = false;
    }
    MoFEMFunctionReturn(0);
  }

private:
  Range entX;
  Range entY;
  Range entZ;
  bool eraseFlag;
  MoFEM::Interface &mField;
  std::string problemName;
};

// }; // namespace OpContactTools