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 visH;
  double scale_constraint;
 
  MatrixDouble mtD; // for stiffness matrix
 
  // for contact
  double cn_cont;
  double rollers_stop_time;
 
  // for rotation
  static vector<double> gravity;
  static double spring_stiffness;
 
  static int closest_roller_id;
  static double delta;
  static Tensor1<double, 3> omega;
  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(1), young_modulus_inv(1), poisson(0.), density(1),
        sigmaY(1e9), cn_pl(1), H(0), C1_k(0), Q_inf(0), b_iso(0), visH(0),
        cn_cont(1. / young_modulus), /*spring_stiffness(0),
                                        closest_roller_id(0),*/
        rollers_stop_time(INFINITY), scale_constraint(1) {
    mtD.resize(36, 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 1200 -mat_mumps_icntl_24 1 -mat_mumps_icntl_13 1 ");
 
    char mumps_options[] = "-fieldsplit_0_mat_mumps_icntl_14 800 "
                           "-fieldsplit_0_mat_mumps_icntl_24 1 "
                           "-fieldsplit_0_mat_mumps_icntl_13 1 "
                           "-fieldsplit_1_mat_mumps_icntl_14 800 "
                           "-fieldsplit_1_mat_mumps_icntl_24 1 "
                           "-fieldsplit_1_mat_mumps_icntl_13 1";
    CHKERR PetscOptionsInsertString(NULL, mumps_options);
 
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "none");
 
    int nb_ent_grav = 3;
    CHKERR PetscOptionsGetRealArray(NULL, NULL, "-gravity", gravity.data(),
                                    &nb_ent_grav, 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);
 
    PetscBool with_contact = PETSC_FALSE;
    CHKERR PetscOptionsBool("-with_contact", "run calculations with contact",
                            "", with_contact, &with_contact, 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) = t1_omega(i);
      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);
    CHKERR PetscOptionsScalar("-visH", "viscous hardening", "", visH, &visH,
                              PETSC_NULL);
    // for contact
    CHKERR PetscOptionsScalar("-spring",
                              "Springs stiffness on the boundary (contact)", "",
                              spring_stiffness, &spring_stiffness, PETSC_NULL);
 
    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", "cone",
                               "torus", "roller", "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 CONE:
          rollerDataVec.push_back(new ConeRigidBody(center_coords, disp, ii));
          break;
        case TORUS:
          rollerDataVec.push_back(new TorusRigidBody(center_coords, disp, ii));
          break;
        case ROLLER:
          rollerDataVec.push_back(
              new RollerRigidBody(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();
 
        char pos_data_file[255];
        PetscBool ctg_flag = PETSC_FALSE;
        param = "-position_data" + to_string(ii);
        CHKERR PetscOptionsString(param.c_str(), "", "", "", pos_data_file, 255,
                                  &ctg_flag);
        if (ctg_flag) {
          rollerDataVec.back().positionDataParamName = string(param);
          rollerDataVec.back().originCoords.clear();
          rollerDataVec.back().rollerDisp.clear();
        }
      }
    }
 
    if (with_contact && rollerDataVec.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "provide at least one body data for contact (-body1 "
              "sphere    -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;
    visH *= young_modulus_inv;
    
    gravity[0] *= young_modulus_inv;
    gravity[1] *= young_modulus_inv;
    gravity[2] *= young_modulus_inv;
    // cn_cont = 1; // perhaps these should be always assigned by user
    // cn_pl = 1;
 
    MoFEMFunctionReturnHot(0);
  }
};
 
extern BlockParamData *cache;
extern std::map<int, BlockParamData> mat_blocks; // for use in future
 
struct DataFromMove {
  size_t comp;
  VectorDouble values;
  VectorDouble scaledValues;
  Range ents;
 
  DataFromMove(size_t c, VectorDouble val, Range ents_bc)
      : comp(c), values(val), scaledValues(val), ents(ents_bc) {}
  DataFromMove(VectorDouble val, Range ents_bc)
      : comp(0), values(val), scaledValues(val), ents(ents_bc) {}
  
  DataFromMove() = delete;
};
 
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);
  }
};
 
/**
 * \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::shared_ptr<DirichletDisplacementBc> dirichletBcPtr;
  boost::shared_ptr<MethodForForceScaling> methodsOpForMove;
  boost::shared_ptr<MethodForForceScaling> methodsOpForRollersDisp;
  vector<boost::shared_ptr<MethodForForceScaling>> methodsOpForRollersPosition;
 
  boost::shared_ptr<NestedMatrixContainer> nestedContainer;
 
  SmartPetscObj<Mat> SEpEp;
  SmartPetscObj<AO> aoSEpEp;
 
  bool scaleMat;
  bool runOnceFlag;
 
  struct BcEntMethod : EntityMethod {
    // MoFEM::Interface &mField;
    DataFromBc &dataFromDirichletBc;
    boost::shared_ptr<DirichletDisplacementBc> dirichletBcPtr;
    BcEntMethod(DataFromBc &data_from_dirichlet_bc,
                boost::shared_ptr<DirichletDisplacementBc> dirichlet_bc)
        : dataFromDirichletBc(data_from_dirichlet_bc),
          dirichletBcPtr(dirichlet_bc) {}
    MoFEMErrorCode preProcess() { return 0; }
    MoFEMErrorCode postProcess() { return 0; }
    MoFEMErrorCode operator()() {
      MoFEMFunctionBegin;
      EntityHandle v = entPtr->getEnt();
      auto &bc_it = dataFromDirichletBc;
      for (int i = 0; i != 3; i++) {
        CHKERR dirichletBcPtr->calculateRotationForDof(v, bc_it);
        if (bc_it.bc_flags[i]) {
          if (entPtr->getEntType() == MBVERTEX) {
            entPtr->getEntFieldData()[i] = bc_it.scaled_values[i];
          } else {
            entPtr->getEntFieldData()[i] = 0;
          }
        }
      }
 
      MoFEMFunctionReturn(0);
    }
  };
 
  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 setNestedContainer(
      boost::shared_ptr<NestedMatrixContainer> nested_container) {
    MoFEMFunctionBegin;
    nestedContainer = nested_container;
    MoFEMFunctionReturn(0);
  }
 
  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);
      // FIXME: this is probably not necessary
      // this is only for keeping constant matrices
      if (nestedContainer) {
        auto &nested_matrices = nestedContainer->nested_matrices;
        if (nested_matrices(0, 0)) {
          CHKERR MatZeroEntries(nested_matrices(0, 0));
          CHKERR MatZeroEntries(nested_matrices(0, 1));
          CHKERR MatZeroEntries(nested_matrices(1, 1));
        }
      }
      if (SEpEp)
        CHKERR MatZeroEntries(SEpEp);
      break;
    }
    case CTX_TSSETIFUNCTION: {
      snes_ctx = CTX_SNESSETFUNCTION;
      snes_f = ts_F;
      break;
    }
    default:
      break;
    }
 
    auto time = ts_t;
 
    // forceDataScaled = forceData;
    // pressureDataScaled = pressureData;
    if (methodsOpForMove)
      for (auto &bdata : bcData) {
        bdata.scaledValues = bdata.values;
        CHKERR MethodForForceScaling::applyScale(this, methodsOpForMove,
                                                 bdata.scaledValues);
      }
 
    if (time <= (*cache).rollers_stop_time) {
 
      if (methodsOpForRollersDisp)
        for (auto &roller : (*cache).rollerDataVec) {
          roller.BodyDispScaled = roller.rollerDisp;
          CHKERR MethodForForceScaling::applyScale(
              this, methodsOpForRollersDisp, roller.BodyDispScaled);
        }
 
      if (!methodsOpForRollersPosition.empty())
 
        for (int id = 0; id != (*cache).rollerDataVec.size(); ++id) {
          auto &roller = (*cache).rollerDataVec[id];
          if (!roller.positionDataParamName.empty())
            CHKERR MethodForForceScaling::applyScale(
                this, methodsOpForRollersPosition[id], roller.BodyDispScaled);
        }
    }
 
    auto fix_dofs_values = [this]() {
      MoFEMFunctionBeginHot;
      auto &methodsOp = dirichletBcPtr->methodsOp;
 
      // for fixing dirichlet bcs
      vector<DataFromBc> bcData;
      CHKERR dirichletBcPtr->getRotationBcFromBlock(bcData);
      CHKERR dirichletBcPtr->getBcDataFromSetsAndBlocks(bcData);
 
      for (auto &bc_it : bcData) {
        Range all_bc_ents;
        CHKERR MethodForForceScaling::applyScale(this, methodsOp,
                                                 bc_it.scaled_values);
        bc_it.theta = bc_it.scaled_values[0];
        BcEntMethod method(bc_it, dirichletBcPtr);
        for (auto &ents : bc_it.bc_ents)
          all_bc_ents.merge(ents);
 
        CHKERR mField.loop_entities("U", method, &all_bc_ents);
      }
 
      MoFEMFunctionReturnHot(0);
    };
 
    CHKERR fix_dofs_values();
 
    // if (snes_ctx == CTX_SNESNONE) {
    //   if (!dofsIndices.empty()) {
    //     CHKERR VecSetValues(snes_x, dofsIndices.size(),
    //     &*dofsIndices.begin(),
    //                         &*dofsValues.begin(), INSERT_VALUES);
    //   }
    //   CHKERR VecAssemblyBegin(snes_x);
    //   CHKERR VecAssemblyEnd(snes_x);
    // }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBeginHot;
    // if (ts_ctx == CTX_TSSETIJACOBIAN)
    //   CHKERR MatView(ts_B, PETSC_VIEWER_STDOUT_SELF);
    auto assemble = [](Mat a) {
      MoFEMFunctionBeginHot;
      if (a) {
        CHKERR MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY);
        CHKERR MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY);
      }
      MoFEMFunctionReturnHot(0);
    };
 
    if (ts_ctx == CTX_TSSETIJACOBIAN) {
 
      // FIXME: this is probably not necessary
      // this is only for keeping constant matrices
      if (nestedContainer) {
        auto &nested_matrices = nestedContainer->nested_matrices;
        if (nested_matrices(0, 0)) {
          CHKERR assemble(nested_matrices(0, 0));
          CHKERR assemble(nested_matrices(0, 1));
          CHKERR assemble(nested_matrices(1, 1));
        }
 
        // cout << " nested_matrices(0, 0) " << endl;
        // CHKERR MatView(nested_matrices(0, 0), PETSC_VIEWER_STDOUT_SELF);
        // cout << " nested_matrices(0, 1) " << endl;
        // CHKERR MatView(nested_matrices(0, 1), PETSC_VIEWER_STDOUT_SELF);
        // cout << " nested_matrices(1, 1) " << endl;
        // CHKERR MatView(nested_matrices(1, 1), PETSC_VIEWER_STDOUT_SELF);
      }
      if (SEpEp)
        CHKERR assemble(SEpEp);
      // if (true)
      //   CHKERR MatView(ts_B, PETSC_VIEWER_DRAW_WORLD);
      // if (SEpEp)
      //   CHKERR MatView(SEpEp, PETSC_VIEWER_DRAW_WORLD);
 
    }
 
    MoFEMFunctionReturnHot(0);
  }
};
 
// }; // namespace OpContactTools