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 -mat_mumps_icntl_20 0");
 
    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);
  }
};
 
 
// }; // namespace OpContactTools