From f66179f3364d2df1fc034aaf37f84e95d7b31712 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 28 Oct 2022 16:38:32 -0600 Subject: [PATCH] Removing old contact files, fixing capitalization in dump_custom.cpp --- src/contact.cpp | 479 ------------------------------ src/contact.h | 107 ------- src/contact_damping_models.cpp | 113 ------- src/contact_damping_models.h | 79 ----- src/contact_heat_models.cpp | 63 ---- src/contact_heat_models.h | 53 ---- src/contact_normal_models.cpp | 377 ----------------------- src/contact_normal_models.h | 118 -------- src/contact_rolling_models.cpp | 121 -------- src/contact_rolling_models.h | 53 ---- src/contact_sub_models.cpp | 140 --------- src/contact_sub_models.h | 62 ---- src/contact_tangential_models.cpp | 412 ------------------------- src/contact_tangential_models.h | 113 ------- src/contact_twisting_models.cpp | 124 -------- src/contact_twisting_models.h | 64 ---- src/dump_custom.cpp | 6 +- 17 files changed, 3 insertions(+), 2481 deletions(-) delete mode 100644 src/contact.cpp delete mode 100644 src/contact.h delete mode 100644 src/contact_damping_models.cpp delete mode 100644 src/contact_damping_models.h delete mode 100644 src/contact_heat_models.cpp delete mode 100644 src/contact_heat_models.h delete mode 100644 src/contact_normal_models.cpp delete mode 100644 src/contact_normal_models.h delete mode 100644 src/contact_rolling_models.cpp delete mode 100644 src/contact_rolling_models.h delete mode 100644 src/contact_sub_models.cpp delete mode 100644 src/contact_sub_models.h delete mode 100644 src/contact_tangential_models.cpp delete mode 100644 src/contact_tangential_models.h delete mode 100644 src/contact_twisting_models.cpp delete mode 100644 src/contact_twisting_models.h diff --git a/src/contact.cpp b/src/contact.cpp deleted file mode 100644 index ae294d7db9..0000000000 --- a/src/contact.cpp +++ /dev/null @@ -1,479 +0,0 @@ -// clang-format off -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - This class contains a series of tools for DEM contacts - Multiple models can be defined and used to calculate forces - and torques based on contact geometry - - Contributing authors: - Dan Bolintineanu (SNL), Joel Clemmer (SNL) ------------------------------------------------------------------------ */ - -#include "contact.h" -#include "contact_sub_models.h" -#include "contact_normal_models.h" -#include "contact_tangential_models.h" -#include "contact_damping_models.h" -#include "contact_rolling_models.h" -#include "contact_twisting_models.h" -#include "contact_heat_models.h" -#include "comm.h" -#include "error.h" -#include "force.h" -#include "math_extra.h" - -#include - -using namespace LAMMPS_NS; -using namespace Contact; -using namespace MathExtra; - -ContactModel::ContactModel(LAMMPS *lmp) : Pointers(lmp) -{ - limit_damping = 0; - beyond_contact = 0; - nondefault_history_transfer = 0; - classic_model = 0; - contact_type = PAIR; - - normal_model = nullptr; - damping_model = nullptr; - tangential_model = nullptr; - rolling_model = nullptr; - twisting_model = nullptr; - heat_model = nullptr; - - for (int i = 0; i < NSUBMODELS; i++) sub_models[i] = nullptr; - transfer_history_factor = nullptr; -} - -/* ---------------------------------------------------------------------- */ - -ContactModel::~ContactModel() -{ - delete [] transfer_history_factor; - - delete normal_model; - delete damping_model; - delete tangential_model; - delete rolling_model; - delete twisting_model; - delete heat_model; -} - -/* ---------------------------------------------------------------------- */ - -void ContactModel::init_model(std::string model_name, ModelType model_type) -{ - if (model_type == NORMAL) { - delete normal_model; - if (model_name == "none") normal_model = new NormalNone(lmp); - else if (model_name == "hooke") normal_model = new NormalHooke(lmp); - else if (model_name == "hertz") normal_model = new NormalHertz(lmp); - else if (model_name == "hertz/material") normal_model = new NormalHertzMaterial(lmp); - else if (model_name == "dmt") normal_model = new NormalDMT(lmp); - else if (model_name == "jkr") normal_model = new NormalJKR(lmp); - else error->all(FLERR, "Normal model name {} not recognized", model_name); - sub_models[model_type] = normal_model; - - } else if (model_type == TANGENTIAL) { - delete tangential_model; - if (model_name == "none") tangential_model = new TangentialNone(lmp); - else if (model_name == "linear_nohistory") tangential_model = new TangentialLinearNoHistory(lmp); - else if (model_name == "linear_history") tangential_model = new TangentialLinearHistory(lmp); - else if (model_name == "linear_history_classic") tangential_model = new TangentialLinearHistoryClassic(lmp); - else if (model_name == "mindlin") tangential_model = new TangentialMindlin(lmp); - else if (model_name == "mindlin/force") tangential_model = new TangentialMindlinForce(lmp); - else if (model_name == "mindlin_rescale") tangential_model = new TangentialMindlinRescale(lmp); - else if (model_name == "mindlin_rescale/force") tangential_model = new TangentialMindlinRescaleForce(lmp); - else error->all(FLERR, "Tangential model name {} not recognized", model_name); - sub_models[model_type] = tangential_model; - - } else if (model_type == DAMPING) { - delete damping_model; - if (model_name == "none") damping_model = new DampingNone(lmp); - else if (model_name == "velocity") damping_model = new DampingVelocity(lmp); - else if (model_name == "mass_velocity") damping_model = new DampingMassVelocity(lmp); - else if (model_name == "viscoelastic") damping_model = new DampingViscoelastic(lmp); - else if (model_name == "tsuji") damping_model = new DampingTsuji(lmp); - else error->all(FLERR, "Damping model name {} not recognized", model_name); - sub_models[model_type] = damping_model; - - } else if (model_type == ROLLING) { - delete rolling_model; - rolling_defined = 1; - if (model_name == "none") { - rolling_model = new RollingNone(lmp); - rolling_defined = 0; - } else if (model_name == "sds") rolling_model = new RollingSDS(lmp); - else error->all(FLERR, "Rolling model name {} not recognized", model_name); - sub_models[model_type] = rolling_model; - - } else if (model_type == TWISTING) { - delete twisting_model; - twisting_defined = 1; - if (model_name == "none") { - twisting_model = new TwistingNone(lmp); - twisting_defined = 0; - } else if (model_name == "sds") twisting_model = new TwistingSDS(lmp); - else if (model_name == "marshall") twisting_model = new TwistingMarshall(lmp); - else error->all(FLERR, "Twisting model name {} not recognized", model_name); - sub_models[model_type] = twisting_model; - - } else if (model_type == HEAT) { - delete heat_model; - heat_defined = 1; - if (model_name == "none") { - heat_model = new HeatNone(lmp); - heat_defined = 0; - } else if (model_name == "area") heat_model = new HeatArea(lmp); - else error->all(FLERR, "Heat model name not {} recognized", model_name); - sub_models[model_type] = heat_model; - } else { - error->all(FLERR, "Illegal model type {}", model_type); - } - - sub_models[model_type]->name.assign(model_name); - sub_models[model_type]->contact = this; - sub_models[model_type]->allocate_coeffs(); -} - -/* ---------------------------------------------------------------------- */ - -int ContactModel::init_classic_model(char **arg, int iarg, int narg) -{ - double kn, kt, gamman, gammat, xmu; - - classic_model = 1; - - if (iarg + 6 >= narg) - error->all(FLERR,"Insufficient arguments provided for classic gran model command"); - - kn = utils::numeric(FLERR,arg[iarg + 1],false,lmp); - if (strcmp(arg[iarg + 2],"NULL") == 0) kt = kn * 2.0 / 7.0; - else kt = utils::numeric(FLERR,arg[iarg + 2],false,lmp); - - gamman = utils::numeric(FLERR,arg[iarg + 3],false,lmp); - if (strcmp(arg[iarg + 4],"NULL") == 0) gammat = 0.5 * gamman; - else gammat = utils::numeric(FLERR,arg[iarg + 4],false,lmp); - - xmu = utils::numeric(FLERR,arg[iarg + 5],false,lmp); - int dampflag = utils::inumeric(FLERR,arg[iarg + 6],false,lmp); - if (dampflag == 0) gammat = 0.0; - - if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 || - xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1) - error->all(FLERR,"Illegal classic gran model command"); - - if (strcmp(arg[iarg],"hooke") == 0) { - init_model("hooke", NORMAL); - init_model("linear_nohistory", TANGENTIAL); - init_model("mass_velocity", DAMPING); - } else if (strcmp(arg[iarg],"hooke/history") == 0) { - init_model("hooke", NORMAL); - init_model("linear_history_classic", TANGENTIAL); - init_model("mass_velocity", DAMPING); - } else if (strcmp(arg[iarg],"hertz/history") == 0) { - // convert Kn and Kt from pressure units to force/distance^2 if Hertzian - kn /= force->nktv2p; - kt /= force->nktv2p; - init_model("hertz", NORMAL); - init_model("mindlin", TANGENTIAL); - init_model("viscoelastic", DAMPING); - } else error->all(FLERR,"Invalid classic gran model"); - - // ensure additional models are undefined - init_model("none", ROLLING); - init_model("none", TWISTING); - init_model("none", HEAT); - - // manually parse coeffs - normal_model->coeffs[0] = kn; - normal_model->coeffs[1] = gamman; - tangential_model->coeffs[0] = kt; - tangential_model->coeffs[1] = gammat / gamman; - tangential_model->coeffs[2] = xmu; - - normal_model->coeffs_to_local(); - tangential_model->coeffs_to_local(); - damping_model->coeffs_to_local(); - - iarg += 7; - return iarg; -} - -/* ---------------------------------------------------------------------- */ - -void ContactModel::init() -{ - int i, j; - for (i = 0; i < NSUBMODELS; i++) - if (!sub_models[i]) init_model("none", (ModelType) i); - - // Must have valid normal, damping, and tangential models - if (normal_model->name == "none") error->all(FLERR, "Must specify normal contact model"); - if (damping_model->name == "none") error->all(FLERR, "Must specify damping contact model"); - if (tangential_model->name == "none") error->all(FLERR, "Must specify tangential contact model"); - - int size_cumulative; - size_history = 0; - for (i = 0; i < NSUBMODELS; i++) { - if (sub_models[i]->nondefault_history_transfer) - nondefault_history_transfer = 1; - if (sub_models[i]->beyond_contact) - beyond_contact = 1; - size_history += sub_models[i]->size_history; - } - - if (nondefault_history_transfer) { - transfer_history_factor = new double[size_history]; - - for (i = 0; i < size_history; i++) { - // Find which model owns this history value - size_cumulative = 0; - for (j = 0; j < NSUBMODELS; j++) { - if (size_cumulative + sub_models[j]->size_history > i) break; - size_cumulative += sub_models[j]->size_history; - } - - // Check if model has nondefault transfers, if so copy its array - transfer_history_factor[i] = -1; - if (j != NSUBMODELS) { - if (sub_models[j]->nondefault_history_transfer) { - transfer_history_factor[i] = sub_models[j]->transfer_history_factor[i - size_cumulative]; - } - } - } - } - - for (i = 0; i < NSUBMODELS; i++) sub_models[i]->init(); -} - -/* ---------------------------------------------------------------------- */ - -int ContactModel::mix_coeffs(ContactModel *c1, ContactModel *c2) -{ - int i; - for (i = 0; i < NSUBMODELS; i++) { - if (c1->sub_models[i]->name != c2->sub_models[i]->name) return i; - - init_model(c1->sub_models[i]->name, (ModelType) i); - sub_models[i]->mix_coeffs(c1->sub_models[i]->coeffs, c2->sub_models[i]->coeffs); - } - - limit_damping = MAX(c1->limit_damping, c2->limit_damping); - - return -1; -} - -/* ---------------------------------------------------------------------- */ - -void ContactModel::write_restart(FILE *fp) -{ - int num_char, num_coeffs; - - for (int i = 0; i < NSUBMODELS; i++) { - num_char = sub_models[i]->name.length(); - num_coeffs = sub_models[i]->num_coeffs; - fwrite(&num_char, sizeof(int), 1, fp); - fwrite(sub_models[i]->name.data(), sizeof(char), num_char, fp); - fwrite(&num_coeffs, sizeof(int), 1, fp); - fwrite(sub_models[i]->coeffs, sizeof(double), num_coeffs, fp); - } -} - -/* ---------------------------------------------------------------------- */ - -void ContactModel::read_restart(FILE *fp) -{ - int num_char, num_coeff; - - for (int i = 0; i < NSUBMODELS; i++) { - if (comm->me == 0) - utils::sfread(FLERR, &num_char, sizeof(int), 1, fp, nullptr, error); - MPI_Bcast(&num_char, 1, MPI_INT, 0, world); - - std::string model_name (num_char, ' '); - if (comm->me == 0) - utils::sfread(FLERR, const_cast(model_name.data()), sizeof(char),num_char, fp, nullptr, error); - MPI_Bcast(const_cast(model_name.data()), num_char, MPI_CHAR, 0, world); - - init_model(model_name, (ModelType) i); - - if (comm->me == 0) { - utils::sfread(FLERR, &num_coeff, sizeof(int), 1, fp, nullptr, error); - if (num_coeff != sub_models[i]->num_coeffs) - error->one(FLERR, "Invalid contact model written to restart file"); - } - MPI_Bcast(&num_coeff, 1, MPI_INT, 0, world); - - if (comm->me == 0) { - utils::sfread(FLERR, sub_models[i]->coeffs, sizeof(double), num_coeff, fp, nullptr, error); - } - MPI_Bcast(sub_models[i]->coeffs, num_coeff, MPI_DOUBLE, 0, world); - sub_models[i]->coeffs_to_local(); - } -} - -/* ---------------------------------------------------------------------- */ - -bool ContactModel::check_contact() -{ - if (contact_type == WALL) { - // Used by fix_wall_gran.cpp - // radj = radius of wall - // dx already provided - rsq = lensq3(dx); - radsum = radi; - if (radj == 0) Reff = radi; - else Reff = radi * radj / (radi + radj); - } else if (contact_type == WALLREGION) { - // Used by fix_wall_gran_region.cpp - // radj = radius of wall - // dx and r already provided - rsq = r * r; - radsum = radi; - if (radj == 0) Reff = radi; - else Reff = radi * radj / (radi + radj); - } else { - sub3(xi, xj, dx); - rsq = lensq3(dx); - radsum = radi + radj; - Reff = radi * radj / radsum; - } - - touch = normal_model->touch(); - return touch; -} - -/* ---------------------------------------------------------------------- */ - -void ContactModel::prep_contact() -{ - double temp[3]; - - // Standard geometric quantities - if (contact_type != WALLREGION) r = sqrt(rsq); - rinv = 1.0 / r; - delta = radsum - r; - dR = delta * Reff; - scale3(rinv, dx, nx); - - // relative translational velocity - sub3(vi, vj, vr); - - // normal component - vnnr = dot3(vr, nx); //v_R . n - scale3(vnnr, nx, vn); - - // tangential component - sub3(vr, vn, vt); - - // relative rotational velocity - scaleadd3(radi, omegai, radj, omegaj, wr); - - // relative tangential velocities - cross3(wr, nx, temp); - sub3(vt, temp, vtr); - vrel = len3(vtr); - - if (rolling_defined || twisting_defined) - sub3(omegai, omegaj, relrot); - - if (rolling_defined) { - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // this is different from the Marshall papers, which use the Bagi/Kuhn formulation - // for rolling velocity (see Wang et al for why the latter is wrong) - vrl[0] = Reff * (relrot[1] * nx[2] - relrot[2] * nx[1]); - vrl[1] = Reff * (relrot[2] * nx[0] - relrot[0] * nx[2]); - vrl[2] = Reff * (relrot[0] * nx[1] - relrot[1] * nx[0]); - } - - if (twisting_defined) { - // omega_T (eq 29 of Marshall) - magtwist = dot3(relrot, nx); - } -} - -/* ---------------------------------------------------------------------- */ - -void ContactModel::calculate_forces() -{ - // calculate forces/torques - - forces[0] = 0.0; - double Fne, Fdamp, dist_to_contact; - area = normal_model->calculate_area(); - normal_model->set_knfac(); - Fne = normal_model->calculate_forces(); - - Fdamp = damping_model->calculate_forces(); - Fntot = Fne + Fdamp; - if (limit_damping && Fntot < 0.0) Fntot = 0.0; - - normal_model->set_fncrit(); // Needed for tangential, rolling, twisting - tangential_model->calculate_forces(); - if (rolling_defined) rolling_model->calculate_forces(); - if (twisting_defined) twisting_model->calculate_forces(); - - // sum contributions - - scale3(Fntot, nx, forces); - add3(forces, fs, forces); - - //May need to rethink eventually tris.. - cross3(nx, fs, torquesi); - scale3(-1, torquesi); - if (contact_type == PAIR) copy3(torquesi, torquesj); - - if (!classic_model && contact_type == PAIR) { - dist_to_contact = radi - 0.5 * delta; - scale3(dist_to_contact, torquesi); - dist_to_contact = radj - 0.5 * delta; - scale3(dist_to_contact, torquesj); - } else { - dist_to_contact = radi; - scale3(dist_to_contact, torquesi); - } - - double torroll[3]; - if (rolling_defined) { - cross3(nx, fr, torroll); - scale3(Reff, torroll); - add3(torquesi, torroll, torquesi); - if (contact_type == PAIR) sub3(torquesj, torroll, torquesj); - } - - double tortwist[3]; - if (twisting_defined) { - scale3(magtortwist, nx, tortwist); - add3(torquesi, tortwist, torquesi); - if (contact_type == PAIR) sub3(torquesj, tortwist, torquesj); - } - - if (heat_defined) { - dq = heat_model->calculate_heat(); - } -} - -/* ---------------------------------------------------------------------- - compute pull-off distance (beyond contact) for a given radius and atom type - use temporary variables since this does not use a specific contact geometry -------------------------------------------------------------------------- */ - -double ContactModel::pulloff_distance(double radi, double radj) -{ - return normal_model->pulloff_distance(radi, radj); -} diff --git a/src/contact.h b/src/contact.h deleted file mode 100644 index d1a13c5eb1..0000000000 --- a/src/contact.h +++ /dev/null @@ -1,107 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef LMP_CONTACT_H -#define LMP_CONTACT_H - -#include "pointers.h" // IWYU pragma: export - -namespace LAMMPS_NS { -namespace Contact { - -#define EPSILON 1e-10 -#define NSUBMODELS 6 - -enum ModelType { - NORMAL = 0, - DAMPING = 1, - TANGENTIAL = 2, - ROLLING = 3, - TWISTING = 4, - HEAT = 5 -}; // Relative order matters since some derive coeffs from others - -enum ContactType { - PAIR = 0, - WALL = 1, - WALLREGION = 2 -}; - -// forward declaration -class NormalModel; -class DampingModel; -class TangentialModel; -class RollingModel; -class TwistingModel; -class HeatModel; -class SubModel; - -class ContactModel : protected Pointers { - public: - ContactModel(class LAMMPS *); - ~ContactModel(); - void init(); - bool check_contact(); - void prep_contact(); - void calculate_forces(); - double pulloff_distance(double, double); - - void init_model(std::string, ModelType); - int init_classic_model(char **, int, int); - int mix_coeffs(ContactModel*, ContactModel*); - - void write_restart(FILE *); - void read_restart(FILE *); - - // Sub models - NormalModel *normal_model; - DampingModel *damping_model; - TangentialModel *tangential_model; - RollingModel *rolling_model; - TwistingModel *twisting_model; - HeatModel *heat_model; - SubModel *sub_models[NSUBMODELS]; // Need to resize if we add more model flavors - - // Extra options - int beyond_contact, limit_damping, history_update; - ContactType contact_type; - - // History variables - int size_history, nondefault_history_transfer; - double *transfer_history_factor; - double *history; - - // Contact properties/output - double forces[3], torquesi[3], torquesj[3], dq; - - double radi, radj, meff, dt, Ti, Tj, area; - double Fntot, magtortwist; - - int i, j; - double *xi, *xj, *vi, *vj, *omegai, *omegaj; - double fs[3], fr[3], ft[3]; - - double dx[3], nx[3], r, rsq, rinv, Reff, radsum, delta, dR; - double vr[3], vn[3], vnnr, vt[3], wr[3], vtr[3], vrl[3], relrot[3], vrel; - double magtwist; - bool touch; - - protected: - int rolling_defined, twisting_defined, heat_defined; // Used to quickly skip undefined submodels - int classic_model; -}; - -} // namespace Contact -} // namespace LAMMPS_NS - -#endif diff --git a/src/contact_damping_models.cpp b/src/contact_damping_models.cpp deleted file mode 100644 index 07b6b88649..0000000000 --- a/src/contact_damping_models.cpp +++ /dev/null @@ -1,113 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "contact_damping_models.h" -#include "contact_normal_models.h" -#include "contact.h" -#include "math_special.h" - -using namespace LAMMPS_NS; -using namespace Contact; -using namespace MathSpecial; - -/* ---------------------------------------------------------------------- - Default damping model -------------------------------------------------------------------------- */ - -DampingModel::DampingModel(LAMMPS *lmp) : SubModel(lmp) {} - -/* ---------------------------------------------------------------------- */ - -void DampingModel::init() -{ - damp = contact->normal_model->damp; -} - -/* ---------------------------------------------------------------------- - No model -------------------------------------------------------------------------- */ - -DampingNone::DampingNone(LAMMPS *lmp) : DampingModel(lmp) {} - -/* ---------------------------------------------------------------------- */ - -double DampingNone::calculate_forces() -{ - return 0.0; -} - -/* ---------------------------------------------------------------------- - Velocity damping -------------------------------------------------------------------------- */ - -DampingVelocity::DampingVelocity(LAMMPS *lmp) : DampingModel(lmp) {} - -/* ---------------------------------------------------------------------- */ - -double DampingVelocity::calculate_forces() -{ - damp_prefactor = damp; - return -damp_prefactor * contact->vnnr; -} - -/* ---------------------------------------------------------------------- - Mass velocity damping -------------------------------------------------------------------------- */ - -DampingMassVelocity::DampingMassVelocity(LAMMPS *lmp) : DampingModel(lmp) {} - -/* ---------------------------------------------------------------------- */ - -double DampingMassVelocity::calculate_forces() -{ - damp_prefactor = damp * contact->meff; - return -damp_prefactor * contact->vnnr; -} - -/* ---------------------------------------------------------------------- - Default, viscoelastic damping -------------------------------------------------------------------------- */ - -DampingViscoelastic::DampingViscoelastic(LAMMPS *lmp) : DampingModel(lmp) {} - -/* ---------------------------------------------------------------------- */ - -double DampingViscoelastic::calculate_forces() -{ - damp_prefactor = damp * contact->meff * contact->area; - return -damp_prefactor * contact->vnnr; -} - -/* ---------------------------------------------------------------------- - Tsuji damping -------------------------------------------------------------------------- */ - -DampingTsuji::DampingTsuji(LAMMPS *lmp) : DampingModel(lmp) {} - -/* ---------------------------------------------------------------------- */ - -void DampingTsuji::init() -{ - double tmp = contact->normal_model->damp; - damp = 1.2728 - 4.2783 * tmp + 11.087 * square(tmp); - damp += -22.348 * cube(tmp)+ 27.467 * powint(tmp, 4); - damp += -18.022 * powint(tmp, 5) + 4.8218 * powint(tmp,6); -} - -/* ---------------------------------------------------------------------- */ - -double DampingTsuji::calculate_forces() -{ - damp_prefactor = damp * sqrt(contact->meff * contact->normal_model->knfac); - return -damp_prefactor * contact->vnnr; -} diff --git a/src/contact_damping_models.h b/src/contact_damping_models.h deleted file mode 100644 index 222e901160..0000000000 --- a/src/contact_damping_models.h +++ /dev/null @@ -1,79 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef CONTACT_DAMPING_MODELS_H_ -#define CONTACT_DAMPING_MODELS_H_ - -#include "contact_sub_models.h" -#include "pointers.h" - -namespace LAMMPS_NS { -namespace Contact { - -class DampingModel : public SubModel { - public: - DampingModel(class LAMMPS *); - ~DampingModel() {}; - virtual void coeffs_to_local() {}; - virtual void mix_coeffs(double*, double*) {}; - virtual void init(); - virtual double calculate_forces() = 0; - double damp, damp_prefactor; -}; - -/* ---------------------------------------------------------------------- */ - -class DampingNone : public DampingModel { - public: - DampingNone(class LAMMPS *); - void init() override {}; - double calculate_forces(); -}; - -/* ---------------------------------------------------------------------- */ - -class DampingVelocity : public DampingModel { - public: - DampingVelocity(class LAMMPS *); - double calculate_forces(); -}; - -/* ---------------------------------------------------------------------- */ - -class DampingMassVelocity : public DampingModel { - public: - DampingMassVelocity(class LAMMPS *); - double calculate_forces(); -}; - -/* ---------------------------------------------------------------------- */ - -class DampingViscoelastic : public DampingModel { - public: - DampingViscoelastic(class LAMMPS *); - double calculate_forces(); -}; - -/* ---------------------------------------------------------------------- */ - -class DampingTsuji : public DampingModel { - public: - DampingTsuji(class LAMMPS *); - void init() override; - double calculate_forces(); -}; - -} // namespace Contact -} // namespace LAMMPS_NS - -#endif /*CONTACT_DAMPING_MODELS_H_ */ diff --git a/src/contact_heat_models.cpp b/src/contact_heat_models.cpp deleted file mode 100644 index 1804ece9df..0000000000 --- a/src/contact_heat_models.cpp +++ /dev/null @@ -1,63 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "contact_heat_models.h" -#include "contact.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace Contact; - -/* ---------------------------------------------------------------------- - Default heat conduction -------------------------------------------------------------------------- */ - -HeatModel::HeatModel(LAMMPS *lmp) : SubModel(lmp) {} - -/* ---------------------------------------------------------------------- - Area-based heat conduction -------------------------------------------------------------------------- */ - -HeatNone::HeatNone(LAMMPS *lmp) : HeatModel(lmp) {} - -/* ---------------------------------------------------------------------- */ - -double HeatNone::calculate_heat() -{ - return 0.0; -} - -/* ---------------------------------------------------------------------- - Area-based heat conduction -------------------------------------------------------------------------- */ - -HeatArea::HeatArea(LAMMPS *lmp) : HeatModel(lmp) -{ - num_coeffs = 1; -} - -/* ---------------------------------------------------------------------- */ - -void HeatArea::coeffs_to_local() -{ - conductivity = coeffs[0]; - - if (conductivity < 0.0) error->all(FLERR, "Illegal area heat model"); -} - -/* ---------------------------------------------------------------------- */ - -double HeatArea::calculate_heat() -{ - return conductivity * contact->area * (contact->Tj - contact->Ti); -} diff --git a/src/contact_heat_models.h b/src/contact_heat_models.h deleted file mode 100644 index 5a0d495810..0000000000 --- a/src/contact_heat_models.h +++ /dev/null @@ -1,53 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef CONTACT_HEAT_MODELS_H_ -#define CONTACT_HEAT_MODELS_H_ - -#include "contact_sub_models.h" - -namespace LAMMPS_NS { -namespace Contact { - -class HeatModel : public SubModel { - public: - HeatModel(class LAMMPS *); - ~HeatModel() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; - virtual double calculate_heat() = 0; -}; - -/* ---------------------------------------------------------------------- */ - -class HeatNone : public HeatModel { - public: - HeatNone(class LAMMPS *); - double calculate_heat(); -}; - -/* ---------------------------------------------------------------------- */ - -class HeatArea : public HeatModel { - public: - HeatArea(class LAMMPS *); - void coeffs_to_local() override; - double calculate_heat(); - protected: - double conductivity; -}; - -} // namespace Contact -} // namespace LAMMPS_NS - -#endif /*CONTACT_HEAT_MODELS_H_ */ diff --git a/src/contact_normal_models.cpp b/src/contact_normal_models.cpp deleted file mode 100644 index b377d5d027..0000000000 --- a/src/contact_normal_models.cpp +++ /dev/null @@ -1,377 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "contact_normal_models.h" -#include "contact.h" -#include "error.h" -#include "math_const.h" - -using namespace LAMMPS_NS; -using namespace Contact; -using namespace MathConst; - -#define PI27SQ 266.47931882941264802866 // 27*PI**2 -#define THREEROOT3 5.19615242270663202362 // 3*sqrt(3) -#define SIXROOT6 14.69693845669906728801 // 6*sqrt(6) -#define INVROOT6 0.40824829046386307274 // 1/sqrt(6) -#define FOURTHIRDS (4.0/3.0) // 4/3 -#define ONETHIRD (1.0/3.0) // 1/3 -#define THREEQUARTERS 0.75 // 3/4 - -/* ---------------------------------------------------------------------- - Default normal model -------------------------------------------------------------------------- */ - -NormalModel::NormalModel(LAMMPS *lmp) : SubModel(lmp) -{ - material_properties = 0; -} - -/* ---------------------------------------------------------------------- */ - -bool NormalModel::touch() -{ - bool touchflag = (contact->rsq < contact->radsum * contact->radsum); - return touchflag; -} - -/* ---------------------------------------------------------------------- */ - -double NormalModel::pulloff_distance(double radi, double radj) -{ - //called outside of compute(), do not assume correct geometry defined in contact - return radi + radj; -} - -/* ---------------------------------------------------------------------- */ - -double NormalModel::calculate_area() -{ - return sqrt(contact->dR); -} - -/* ---------------------------------------------------------------------- */ - -void NormalModel::set_fncrit() -{ - Fncrit = fabs(contact->Fntot); -} - -/* ---------------------------------------------------------------------- - No model -------------------------------------------------------------------------- */ - -NormalNone::NormalNone(LAMMPS *lmp) : NormalModel(lmp) {} - -/* ---------------------------------------------------------------------- */ - -double NormalNone::calculate_forces() -{ - return 0.0; -} - -/* ---------------------------------------------------------------------- - Hookean normal force -------------------------------------------------------------------------- */ - -NormalHooke::NormalHooke(LAMMPS *lmp) : NormalModel(lmp) -{ - num_coeffs = 2; -} - -/* ---------------------------------------------------------------------- */ - -void NormalHooke::coeffs_to_local() -{ - k = coeffs[0]; - damp = coeffs[1]; - - if (k < 0.0 || damp < 0.0) error->all(FLERR, "Illegal Hooke normal model"); -} - -/* ---------------------------------------------------------------------- */ - -double NormalHooke::calculate_forces() -{ - Fne = knfac * contact->delta; - return Fne; -} - -/* ---------------------------------------------------------------------- */ - -void NormalHooke::set_knfac() -{ - knfac = k; -} - -/* ---------------------------------------------------------------------- - Hertzian normal force -------------------------------------------------------------------------- */ - -NormalHertz::NormalHertz(LAMMPS *lmp) : NormalModel(lmp) -{ - num_coeffs = 2; -} - -/* ---------------------------------------------------------------------- */ - -void NormalHertz::coeffs_to_local() -{ - k = coeffs[0]; - damp = coeffs[1]; - - if (k < 0.0 || damp < 0.0) error->all(FLERR, "Illegal Hertz normal model"); -} - -/* ---------------------------------------------------------------------- */ - -double NormalHertz::calculate_forces() -{ - Fne = knfac * contact->delta; - return Fne; -} - -/* ---------------------------------------------------------------------- */ - -void NormalHertz::set_knfac() -{ - knfac = k * contact->area; -} - -/* ---------------------------------------------------------------------- - Hertzian normal force with material properties -------------------------------------------------------------------------- */ - -NormalHertzMaterial::NormalHertzMaterial(LAMMPS *lmp) : NormalHertz(lmp) -{ - material_properties = 1; - num_coeffs = 3; -} - -/* ---------------------------------------------------------------------- */ - -void NormalHertzMaterial::coeffs_to_local() -{ - Emod = coeffs[0]; - damp = coeffs[1]; - poiss = coeffs[2]; - if (contact->contact_type == PAIR) { - k = FOURTHIRDS * mix_stiffnessE(Emod, Emod, poiss, poiss); - } else { - k = FOURTHIRDS * mix_stiffnessE_wall(Emod, poiss); - } - - if (Emod < 0.0 || damp < 0.0) error->all(FLERR, "Illegal Hertz material normal model"); -} - -/* ---------------------------------------------------------------------- */ - -void NormalHertzMaterial::mix_coeffs(double* icoeffs, double* jcoeffs) -{ - coeffs[0] = mix_stiffnessE(icoeffs[0], jcoeffs[0],icoeffs[2], jcoeffs[2]); - coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]); - coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]); - coeffs_to_local(); -} - -/* ---------------------------------------------------------------------- - DMT normal force -------------------------------------------------------------------------- */ - -NormalDMT::NormalDMT(LAMMPS *lmp) : NormalModel(lmp) -{ - allow_limit_damping = 0; - material_properties = 1; - num_coeffs = 4; -} - -/* ---------------------------------------------------------------------- */ - -void NormalDMT::coeffs_to_local() -{ - Emod = coeffs[0]; - damp = coeffs[1]; - poiss = coeffs[2]; - cohesion = coeffs[3]; - if (contact->contact_type == PAIR) { - k = FOURTHIRDS * mix_stiffnessE(Emod, Emod, poiss, poiss); - } else { - k = FOURTHIRDS * mix_stiffnessE_wall(Emod, poiss); - } - - if (Emod < 0.0 || damp < 0.0) error->all(FLERR, "Illegal DMT normal model"); -} - -/* ---------------------------------------------------------------------- */ - -void NormalDMT::mix_coeffs(double* icoeffs, double* jcoeffs) -{ - coeffs[0] = mix_stiffnessE(icoeffs[0], jcoeffs[0],icoeffs[2], jcoeffs[2]); - coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]); - coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]); - coeffs[3] = mix_geom(icoeffs[3], jcoeffs[3]); - coeffs_to_local(); -} - -/* ---------------------------------------------------------------------- */ - -double NormalDMT::calculate_forces() -{ - Fne = knfac * contact->delta; - F_pulloff = 4.0 * MathConst::MY_PI * cohesion * contact->Reff; - Fne -= F_pulloff; - return Fne; -} - -/* ---------------------------------------------------------------------- */ - -void NormalDMT::set_knfac() -{ - knfac = k * contact->area; -} - -/* ---------------------------------------------------------------------- */ - -void NormalDMT::set_fncrit() -{ - Fncrit = fabs(Fne + 2.0 * F_pulloff); -} - -/* ---------------------------------------------------------------------- - JKR normal force -------------------------------------------------------------------------- */ - -NormalJKR::NormalJKR(LAMMPS *lmp) : NormalModel(lmp) -{ - allow_limit_damping = 0; - material_properties = 1; - beyond_contact = 1; - num_coeffs = 4; -} - -/* ---------------------------------------------------------------------- */ - -void NormalJKR::coeffs_to_local() -{ - Emod = coeffs[0]; - damp = coeffs[1]; - poiss = coeffs[2]; - cohesion = coeffs[3]; - - if (contact->contact_type == PAIR) { - Emix = mix_stiffnessE(Emod, Emod, poiss, poiss); - } else { - Emix = mix_stiffnessE_wall(Emod, poiss); - } - - k = FOURTHIRDS * Emix; - - if (Emod < 0.0 || damp < 0.0) error->all(FLERR, "Illegal JKR normal model"); -} - -/* ---------------------------------------------------------------------- */ - -void NormalJKR::mix_coeffs(double* icoeffs, double* jcoeffs) -{ - coeffs[0] = mix_stiffnessE(icoeffs[0], jcoeffs[0],icoeffs[2], jcoeffs[2]); - coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]); - coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]); - coeffs[3] = mix_geom(icoeffs[3], jcoeffs[3]); - coeffs_to_local(); -} - -/* ---------------------------------------------------------------------- */ - -bool NormalJKR::touch() -{ - double area_at_pulloff, R2, delta_pulloff, dist_pulloff; - bool touchflag; - - if (contact->touch) { - R2 = contact->Reff * contact->Reff; - area_at_pulloff = cbrt(9.0 * MY_PI * cohesion * R2 / (4.0 * Emix)); - delta_pulloff = area_at_pulloff * area_at_pulloff / contact->Reff - 2.0 * sqrt(MY_PI * cohesion * area_at_pulloff / Emix); - dist_pulloff = contact->radsum - delta_pulloff; - touchflag = contact->rsq < (dist_pulloff * dist_pulloff); - } else { - touchflag = contact->rsq < (contact->radsum * contact->radsum); - } - - return touchflag; -} - -/* ---------------------------------------------------------------------- - called outside of compute(), do not assume geometry defined in contact -------------------------------------------------------------------------- */ - -double NormalJKR::pulloff_distance(double radi, double radj) -{ - double area_at_pulloff, Reff_tmp; - - Reff_tmp = radi * radj / (radi + radj); // May not be defined - if (Reff_tmp <= 0) return 0; - - area_at_pulloff = cbrt(9.0 * MY_PI * cohesion * Reff_tmp * Reff_tmp / (4.0 * Emix)); - return area_at_pulloff * area_at_pulloff / Reff_tmp - 2.0 * sqrt(MY_PI * cohesion * area_at_pulloff / Emix); -} - -/* ---------------------------------------------------------------------- */ - -double NormalJKR::calculate_area() -{ - double R2, dR2, t0, t1, t2, t3, t4, t5, t6; - double sqrt1, sqrt2, sqrt3; - - R2 = contact->Reff * contact->Reff; - dR2 = contact->dR * contact->dR; - t0 = cohesion * cohesion * R2 * R2 * Emix; - t1 = PI27SQ * t0; - t2 = 8.0 * contact->dR * dR2 * Emix * Emix * Emix; - t3 = 4.0 * dR2 * Emix; - - // in case sqrt(0) < 0 due to precision issues - sqrt1 = MAX(0, t0 * (t1 + 2.0 * t2)); - t4 = cbrt(t1 + t2 + THREEROOT3 * MY_PI * sqrt(sqrt1)); - t5 = t3 / t4 + t4 / Emix; - sqrt2 = MAX(0, 2.0 * contact->dR + t5); - t6 = sqrt(sqrt2); - sqrt3 = MAX(0, 4.0 * contact->dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Emix * t6)); - - return INVROOT6 * (t6 + sqrt(sqrt3)); -} - -/* ---------------------------------------------------------------------- */ - -double NormalJKR::calculate_forces() -{ - double a2; - a2 = contact->area * contact->area; - Fne = knfac * a2 / contact->Reff - MY_2PI * a2 * sqrt(4.0 * cohesion * Emix / (MY_PI * contact->area)); - F_pulloff = 3.0 * MY_PI * cohesion * contact->Reff; - - return Fne; -} - -/* ---------------------------------------------------------------------- */ - -void NormalJKR::set_knfac() -{ - knfac = k * contact->area; -} - -/* ---------------------------------------------------------------------- */ - -void NormalJKR::set_fncrit() -{ - Fncrit = fabs(Fne + 2.0 * F_pulloff); -} diff --git a/src/contact_normal_models.h b/src/contact_normal_models.h deleted file mode 100644 index 8159a4c620..0000000000 --- a/src/contact_normal_models.h +++ /dev/null @@ -1,118 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef CONTACT_NORMAL_MODELS_H_ -#define CONTACT_NORMAL_MODELS_H_ - -#include "contact_sub_models.h" - -namespace LAMMPS_NS { -namespace Contact { - -class NormalModel : public SubModel { - public: - NormalModel(class LAMMPS *); - ~NormalModel() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; - virtual bool touch(); - virtual double pulloff_distance(double, double); - virtual double calculate_area(); - virtual void set_knfac() = 0; - virtual double calculate_forces() = 0; - virtual void set_fncrit(); - double damp; // Vestigial argument needed by damping - double Emod, poiss; - double Fncrit, Fne, knfac; - int material_properties; -}; - -/* ---------------------------------------------------------------------- */ - -class NormalNone : public NormalModel { - public: - NormalNone(class LAMMPS *); - void set_knfac() {}; - double calculate_forces(); -}; - -/* ---------------------------------------------------------------------- */ - -class NormalHooke : public NormalModel { - public: - NormalHooke(class LAMMPS *); - void coeffs_to_local() override; - void set_knfac(); - double calculate_forces(); - protected: - double k; -}; - -/* ---------------------------------------------------------------------- */ - -class NormalHertz : public NormalModel { - public: - NormalHertz(class LAMMPS *); - void coeffs_to_local() override; - void set_knfac(); - double calculate_forces(); - protected: - double k; -}; - -/* ---------------------------------------------------------------------- */ - -class NormalHertzMaterial : public NormalHertz { - public: - NormalHertzMaterial(class LAMMPS *); - void coeffs_to_local() override; - void mix_coeffs(double*, double*) override; -}; - -/* ---------------------------------------------------------------------- */ - -class NormalDMT : public NormalModel { - public: - NormalDMT(class LAMMPS *); - void coeffs_to_local() override; - void mix_coeffs(double*, double*) override; - void set_knfac(); - double calculate_forces(); - void set_fncrit() override; - protected: - double k, cohesion; - double F_pulloff; -}; - -/* ---------------------------------------------------------------------- */ - -class NormalJKR : public NormalModel { - public: - NormalJKR(class LAMMPS *); - void coeffs_to_local() override; - void mix_coeffs(double*, double*) override; - bool touch() override; - double pulloff_distance(double, double) override; - double calculate_area() override; - void set_knfac(); - double calculate_forces(); - void set_fncrit() override; - protected: - double k, cohesion; - double Emix, F_pulloff; -}; - -} // namespace Contact -} // namespace LAMMPS_NS - -#endif /*CONTACT_NORMAL_MODELS_H_ */ diff --git a/src/contact_rolling_models.cpp b/src/contact_rolling_models.cpp deleted file mode 100644 index 3ee80c0e48..0000000000 --- a/src/contact_rolling_models.cpp +++ /dev/null @@ -1,121 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "contact_normal_models.h" -#include "contact_rolling_models.h" -#include "contact.h" -#include "error.h" -#include "math_const.h" -#include "math_extra.h" - -using namespace LAMMPS_NS; -using namespace Contact; -using namespace MathConst; -using namespace MathExtra; - -/* ---------------------------------------------------------------------- - Default rolling friction model -------------------------------------------------------------------------- */ - -RollingModel::RollingModel(LAMMPS *lmp) : SubModel(lmp) {} - -/* ---------------------------------------------------------------------- - No model -------------------------------------------------------------------------- */ - -RollingNone::RollingNone(LAMMPS *lmp) : RollingModel(lmp) {} - -/* ---------------------------------------------------------------------- - SDS rolling friction model -------------------------------------------------------------------------- */ - -RollingSDS::RollingSDS(LAMMPS *lmp) : RollingModel(lmp) -{ - num_coeffs = 3; - size_history = 3; -} - -/* ---------------------------------------------------------------------- */ - -void RollingSDS::coeffs_to_local() -{ - k = coeffs[0]; - gamma = coeffs[1]; - mu = coeffs[2]; - - if (k < 0.0 || mu < 0.0 || gamma < 0.0) - error->all(FLERR, "Illegal SDS rolling model"); -} - -/* ---------------------------------------------------------------------- */ - -void RollingSDS::calculate_forces() -{ - int rhist0, rhist1, rhist2, frameupdate; - double Frcrit, rolldotn, rollmag, prjmag, magfr, hist_temp[3], scalefac, temp_array[3]; - double k_inv, magfr_inv; - - rhist0 = history_index; - rhist1 = rhist0 + 1; - rhist2 = rhist1 + 1; - - Frcrit = mu * contact->normal_model->Fncrit; - - if (contact->history_update) { - hist_temp[0] = contact->history[rhist0]; - hist_temp[1] = contact->history[rhist1]; - hist_temp[2] = contact->history[rhist2]; - rolldotn = dot3(hist_temp, contact->nx); - - frameupdate = (fabs(rolldotn) * k) > (EPSILON * Frcrit); - if (frameupdate) { // rotate into tangential plane - rollmag = len3(hist_temp); - // projection - scale3(rolldotn, contact->nx, temp_array); - sub3(hist_temp, temp_array, hist_temp); - - // also rescale to preserve magnitude - prjmag = len3(hist_temp); - if (prjmag > 0) scalefac = rollmag / prjmag; - else scalefac = 0; - scale3(scalefac, hist_temp); - } - scale3(contact->dt, contact->vrl, temp_array); - add3(hist_temp, temp_array, hist_temp); - } - - scaleadd3(-k, hist_temp, -gamma, contact->vrl, contact->fr); - - // rescale frictional displacements and forces if needed - magfr = len3(contact->fr); - if (magfr > Frcrit) { - rollmag = len3(hist_temp); - if (rollmag != 0.0) { - k_inv = 1.0 / k; - magfr_inv = 1.0 / magfr; - scale3(-Frcrit * k_inv * magfr_inv, contact->fr, hist_temp); - scale3(-gamma * k_inv, contact->vrl, temp_array); - add3(hist_temp, temp_array, hist_temp); - - scale3(Frcrit * magfr_inv, contact->fr); - } else { - zero3(contact->fr); - } - } - - if (contact->history_update) { - contact->history[rhist0] = hist_temp[0]; - contact->history[rhist1] = hist_temp[1]; - contact->history[rhist2] = hist_temp[2]; - } -} diff --git a/src/contact_rolling_models.h b/src/contact_rolling_models.h deleted file mode 100644 index f5af08c1b0..0000000000 --- a/src/contact_rolling_models.h +++ /dev/null @@ -1,53 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef CONTACT_ROLLING_MODELS_H_ -#define CONTACT_ROLLING_MODELS_H_ - -#include "contact_sub_models.h" - -namespace LAMMPS_NS { -namespace Contact { - -class RollingModel : public SubModel { - public: - RollingModel(class LAMMPS *); - ~RollingModel() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; - virtual void calculate_forces() = 0; -}; - -/* ---------------------------------------------------------------------- */ - -class RollingNone : public RollingModel { - public: - RollingNone(class LAMMPS *); - void calculate_forces() {}; -}; - -/* ---------------------------------------------------------------------- */ - -class RollingSDS : public RollingModel { - public: - RollingSDS(class LAMMPS *); - void coeffs_to_local() override; - void calculate_forces(); - protected: - double k, mu, gamma; -}; - -} // namespace Contact -} // namespace LAMMPS_NS - -#endif /*CONTACT_ROLLING_MODELS_H_ */ diff --git a/src/contact_sub_models.cpp b/src/contact_sub_models.cpp deleted file mode 100644 index 1c9fe6b4f5..0000000000 --- a/src/contact_sub_models.cpp +++ /dev/null @@ -1,140 +0,0 @@ -// clang-format off -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. --------------------------------------------------------------------------*/ - -/* ---------------------------------------------------------------------- - This class contains a framework for normal, damping, tangential, - rolling, twisting, and heat models used to calculate forces - and torques based on contact geometry - - Contributing authors: - Dan Bolintineanu (SNL), Joel Clemmer (SNL) ------------------------------------------------------------------------ */ - -#include "contact_sub_models.h" -#include "error.h" -#include "utils.h" - -using namespace LAMMPS_NS; -using namespace Contact; - -/* ---------------------------------------------------------------------- - Parent class for all types of contact forces -------------------------------------------------------------------------- */ - -SubModel::SubModel(LAMMPS *lmp) : Pointers(lmp) -{ - allocated = 0; - size_history = 0; - history_index = 0; - allow_limit_damping = 1; - beyond_contact = 0; - num_coeffs = 0; - - nondefault_history_transfer = 0; - transfer_history_factor = nullptr; -} - -/* ---------------------------------------------------------------------- */ - -SubModel::~SubModel() -{ - if (allocated) delete [] coeffs; - delete [] transfer_history_factor; -} - -/* ---------------------------------------------------------------------- */ - -void SubModel::allocate_coeffs() -{ - allocated = 1; - coeffs = new double[num_coeffs]; -} - -/* ---------------------------------------------------------------------- */ - -int SubModel::parse_coeffs(char **arg, int iarg, int narg) -{ - if (iarg + num_coeffs > narg) - error->all(FLERR, "Insufficient arguments provided for {} model", name); - for (int i = 0; i < num_coeffs; i++) { - // A few parameters (eg.g kt for tangential mindlin) allow null - if (strcmp(arg[iarg+i], "NULL") == 0) coeffs[i] = -1; - else coeffs[i] = utils::numeric(FLERR,arg[iarg+i],false,lmp); - } - coeffs_to_local(); - - return iarg + num_coeffs; -} - -/* ---------------------------------------------------------------------- */ - -void SubModel::mix_coeffs(double* icoeffs, double* jcoeffs) -{ - for (int i = 0; i < num_coeffs; i++) - coeffs[i] = mix_geom(icoeffs[i], jcoeffs[i]); - coeffs_to_local(); -} - -/* ---------------------------------------------------------------------- - mixing of Young's modulus (E) -------------------------------------------------------------------------- */ - -double SubModel::mix_stiffnessE(double E1, double E2, - double pois1, double pois2) -{ - double factor1 = (1 - pois1 * pois1) / E1; - double factor2 = (1 - pois2 * pois2) / E2; - return 1 / (factor1 + factor2); -} - -/* ---------------------------------------------------------------------- - mixing of shear modulus (G) ------------------------------------------------------------------------- */ - -double SubModel::mix_stiffnessG(double E1, double E2, - double pois1, double pois2) -{ - double factor1 = 2 * (2 - pois1) * (1 + pois1) / E1; - double factor2 = 2 * (2 - pois2) * (1 + pois2) / E2; - return 1 / (factor1 + factor2); -} - -/* ---------------------------------------------------------------------- - mixing of Young's modulus (E) for walls -------------------------------------------------------------------------- */ - -double SubModel::mix_stiffnessE_wall(double E, double pois) -{ - double factor = 2 * (1 - pois); - return E / factor; -} - -/* ---------------------------------------------------------------------- - mixing of shear modulus (G) for walls ------------------------------------------------------------------------- */ - -double SubModel::mix_stiffnessG_wall(double E, double pois) -{ - double factor = 32.0 * (2 - pois) * (1 + pois); - return E / factor; -} - -/* ---------------------------------------------------------------------- - mixing of everything else -------------------------------------------------------------------------- */ - -double SubModel::mix_geom(double val1, double val2) -{ - return sqrt(val1 * val2); -} diff --git a/src/contact_sub_models.h b/src/contact_sub_models.h deleted file mode 100644 index 0f577bf08e..0000000000 --- a/src/contact_sub_models.h +++ /dev/null @@ -1,62 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef CONTACT_SUB_MODEL_H_ -#define CONTACT_SUB_MODEL_H_ - -#include "contact.h" -#include "pointers.h" // IWYU pragma: export - -namespace LAMMPS_NS { -namespace Contact { - -class SubModel : protected Pointers { - public: - SubModel(class LAMMPS *); - virtual ~SubModel(); - - int num_coeffs; - double *coeffs; - void read_restart(); - int parse_coeffs(char **, int, int); - virtual void mix_coeffs(double*, double*); - virtual void coeffs_to_local() = 0; - virtual void init() = 0; // called after all other submodel coeffs defined - - void allocate_coeffs(); - std::string name; - - int size_history; - int nondefault_history_transfer; - double *transfer_history_factor; - - int history_index; - int beyond_contact; - int allow_limit_damping; - - ContactModel *contact; - - protected: - int allocated; - - double mix_stiffnessE(double, double, double, double); - double mix_stiffnessG(double, double, double, double); - double mix_stiffnessE_wall(double, double); - double mix_stiffnessG_wall(double, double); - double mix_geom(double, double); -}; - -} // namespace Contact -} // namespace LAMMPS_NS - -#endif /* CONTACT_SUB_MODEL_H_ */ diff --git a/src/contact_tangential_models.cpp b/src/contact_tangential_models.cpp deleted file mode 100644 index ca95d1d00b..0000000000 --- a/src/contact_tangential_models.cpp +++ /dev/null @@ -1,412 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "contact_damping_models.h" -#include "contact_normal_models.h" -#include "contact_tangential_models.h" -#include "contact.h" -#include "error.h" -#include "math_const.h" -#include "math_extra.h" - -using namespace LAMMPS_NS; -using namespace Contact; -using namespace MathConst; -using namespace MathExtra; - -/* ---------------------------------------------------------------------- - Default model -------------------------------------------------------------------------- */ - -TangentialModel::TangentialModel(LAMMPS *lmp) : SubModel(lmp) {} - -/* ---------------------------------------------------------------------- - No model -------------------------------------------------------------------------- */ - -TangentialNone::TangentialNone(LAMMPS *lmp) : TangentialModel(lmp) {} - -/* ---------------------------------------------------------------------- - Linear model with no history -------------------------------------------------------------------------- */ - -TangentialLinearNoHistory::TangentialLinearNoHistory(LAMMPS *lmp) : TangentialModel(lmp) -{ - num_coeffs = 2; - size_history = 3; -} - -/* ---------------------------------------------------------------------- */ - -void TangentialLinearNoHistory::coeffs_to_local() -{ - k = 0.0; // No tangential stiffness with no history - xt = coeffs[0]; - mu = coeffs[1]; - - if (k < 0.0 || xt < 0.0 || mu < 0.0) - error->all(FLERR, "Illegal linear no history tangential model"); -} - -/* ---------------------------------------------------------------------- */ - -void TangentialLinearNoHistory::calculate_forces() -{ - // classic pair gran/hooke (no history) - damp = xt * contact->damping_model->damp_prefactor; - - double Fscrit = mu * contact->normal_model->Fncrit; - double fsmag = damp * contact->vrel; - - double Ft; - if (contact->vrel != 0.0) Ft = MIN(Fscrit, fsmag) / contact->vrel; - else Ft = 0.0; - - scale3(-Ft, contact->vtr, contact->fs); -} - -/* ---------------------------------------------------------------------- - Linear model with history -------------------------------------------------------------------------- */ - -TangentialLinearHistory::TangentialLinearHistory(LAMMPS *lmp) : TangentialModel(lmp) -{ - num_coeffs = 3; - size_history = 3; -} - -/* ---------------------------------------------------------------------- */ - -void TangentialLinearHistory::coeffs_to_local() -{ - k = coeffs[0]; - xt = coeffs[1]; - mu = coeffs[2]; - - if (k < 0.0 || xt < 0.0 || mu < 0.0) - error->all(FLERR, "Illegal linear tangential model"); -} - -/* ---------------------------------------------------------------------- */ - -void TangentialLinearHistory::calculate_forces() -{ - // Note: this is the same as the base Mindlin calculation except k isn't scaled by area - double magfs, magfs_inv, rsht, shrmag, prjmag, temp_dbl, temp_array[3]; - int frame_update = 0; - - damp = xt * contact->damping_model->damp_prefactor; - - double Fscrit = contact->normal_model->Fncrit * mu; - double *history = & contact->history[history_index]; - - // rotate and update displacements / force. - // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 - if (contact->history_update) { - rsht = dot3(history, contact->nx); - frame_update = (fabs(rsht) * k) > (EPSILON * Fscrit); - - if (frame_update) { - shrmag = len3(history); - - // projection - scale3(rsht, contact->nx, temp_array); - sub3(history, temp_array, history); - - // also rescale to preserve magnitude - prjmag = len3(history); - if (prjmag > 0) temp_dbl = shrmag / prjmag; - else temp_dbl = 0; - scale3(temp_dbl, history); - } - - // update history, tangential force - // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 - scale3(contact->dt, contact->vtr, temp_array); - add3(history, temp_array, history); - } - - - // tangential forces = history + tangential velocity damping - scale3(-k, history, contact->fs); - scale3(damp, contact->vtr, temp_array); - sub3(contact->fs, temp_array, contact->fs); - - // rescale frictional displacements and forces if needed - magfs = len3(contact->fs); - if (magfs > Fscrit) { - shrmag = len3(history); - if (shrmag != 0.0) { - magfs_inv = 1.0 / magfs; - scale3(Fscrit * magfs_inv, contact->fs, history); - scale3(damp, contact->vtr, temp_array); - add3(history, temp_array, history); - scale3(-1.0 / k, history); - scale3(Fscrit * magfs_inv, contact->fs); - } else { - zero3(contact->fs); - } - } -} - -/* ---------------------------------------------------------------------- - Linear model with history from pair gran/hooke/history -------------------------------------------------------------------------- */ - -TangentialLinearHistoryClassic::TangentialLinearHistoryClassic(LAMMPS *lmp) : TangentialLinearHistory(lmp) -{ - scale_area = 0; // Sets gran/hooke/history behavior -} - -/* ---------------------------------------------------------------------- */ - -void TangentialLinearHistoryClassic::calculate_forces() -{ - double k_scaled, magfs, magfs_inv, rsht, shrmag, prjmag, temp_dbl; - double temp_array[3]; - int frame_update = 0; - - k_scaled = k; - if (scale_area) k_scaled *= contact->area; - - damp = xt * contact->damping_model->damp_prefactor; - - double Fscrit = contact->normal_model->Fncrit * mu; - double *history = & contact->history[history_index]; - - // update history - if (contact->history_update) { - scale3(contact->dt, contact->vtr, temp_array); - add3(history, temp_array, history); - } - - shrmag = len3(history); - - // rotate shear displacements - if (contact->history_update) { - rsht = dot3(history, contact->nx); - scale3(rsht, contact->nx, temp_array); - sub3(history, temp_array, history); - } - - // tangential forces = history + tangential velocity damping - scale3(-k_scaled, history, contact->fs); - scale3(damp, contact->vtr, temp_array); - sub3(contact->fs, temp_array, contact->fs); - - // rescale frictional displacements and forces if needed - magfs = len3(contact->fs); - if (magfs > Fscrit) { - if (shrmag != 0.0) { - magfs_inv = 1.0 / magfs; - scale3(Fscrit * magfs_inv, contact->fs, history); - scale3(damp, contact->vtr, temp_array); - add3(history, temp_array, history); - temp_dbl = -1.0 / k_scaled; - if (scale_area) temp_dbl /= contact->area; - scale3(temp_dbl, history); - scale3(Fscrit * magfs_inv, contact->fs); - } else { - zero3(contact->fs); - } - } -} - -/* ---------------------------------------------------------------------- - Mindlin from pair gran/hertz/history -------------------------------------------------------------------------- */ - -TangentialMindlinClassic::TangentialMindlinClassic(LAMMPS *lmp) : TangentialLinearHistoryClassic(lmp) -{ - scale_area = 1; // Sets gran/hertz/history behavior -} - -/* ---------------------------------------------------------------------- - Mindlin model -------------------------------------------------------------------------- */ - -TangentialMindlin::TangentialMindlin(LAMMPS *lmp) : TangentialModel(lmp) -{ - num_coeffs = 3; - size_history = 3; - mindlin_force = 0; - mindlin_rescale = 0; -} - -/* ---------------------------------------------------------------------- */ - -void TangentialMindlin::coeffs_to_local() -{ - k = coeffs[0]; - xt = coeffs[1]; - mu = coeffs[2]; - - if (k == -1) { - if (!contact->normal_model->material_properties) - error->all(FLERR, "Must either specify tangential stiffness or material properties for normal model for the Mindlin tangential style"); - - double Emod = contact->normal_model->Emod; - double poiss = contact->normal_model->poiss; - - if (contact->contact_type == PAIR) { - k = 8.0 * mix_stiffnessG(Emod, Emod, poiss, poiss); - } else { - k = 8.0 * mix_stiffnessG_wall(Emod, poiss); - } - } - - if (k < 0.0 || xt < 0.0 || mu < 0.0) - error->all(FLERR, "Illegal Mindlin tangential model"); -} - -/* ---------------------------------------------------------------------- */ - -void TangentialMindlin::mix_coeffs(double* icoeffs, double* jcoeffs) -{ - if (icoeffs[0] == -1 || jcoeffs[0] == -1) coeffs[0] = -1; - else coeffs[0] = mix_geom(icoeffs[0], jcoeffs[0]); - coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]); - coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]); - coeffs_to_local(); -} - -/* ---------------------------------------------------------------------- */ - -void TangentialMindlin::calculate_forces() -{ - double k_scaled, magfs, magfs_inv, rsht, shrmag, prjmag, temp_dbl; - double temp_array[3]; - int frame_update = 0; - - damp = xt * contact->damping_model->damp_prefactor; - - double *history = & contact->history[history_index]; - double Fscrit = contact->normal_model->Fncrit * mu; - - k_scaled = k * contact->area; - - // on unloading, rescale the shear displacements/force - if (mindlin_rescale) - if (contact->area < history[3]) - scale3(contact->area / history[3], history); - - // rotate and update displacements / force. - // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 - if (contact->history_update) { - rsht = dot3(history, contact->nx); - if (mindlin_force) { - frame_update = fabs(rsht) > (EPSILON * Fscrit); - } else { - frame_update = (fabs(rsht) * k_scaled) > (EPSILON * Fscrit); - } - - if (frame_update) { - shrmag = len3(history); - // projection - scale3(rsht, contact->nx, temp_array); - sub3(history, temp_array, history); - // also rescale to preserve magnitude - prjmag = len3(history); - if (prjmag > 0) temp_dbl = shrmag / prjmag; - else temp_dbl = 0; - scale3(temp_dbl, history); - } - - // update history - if (mindlin_force) { - // tangential force - // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 - scale3(-k_scaled * contact->dt, contact->vtr, temp_array); - } else { - scale3(contact->dt, contact->vtr, temp_array); - } - add3(history, temp_array, history); - - if (mindlin_rescale) history[3] = contact->area; - } - - // tangential forces = history + tangential velocity damping - scale3(-damp, contact->vtr, contact->fs); - - if (!mindlin_force) { - scale3(k_scaled, history, temp_array); - sub3(contact->fs, temp_array, contact->fs); - } else { - add3(contact->fs, history, contact->fs); - } - - // rescale frictional displacements and forces if needed - magfs = len3(contact->fs); - if (magfs > Fscrit) { - shrmag = len3(history); - if (shrmag != 0.0) { - magfs_inv = 1.0 / magfs; - scale3(Fscrit * magfs_inv, contact->fs, history); - scale3(damp, contact->vtr, temp_array); - add3(history, temp_array, history); - - if (!mindlin_force) - scale3(-1.0 / k_scaled, history); - - scale3(Fscrit * magfs_inv, contact->fs); - } else { - zero3(contact->fs); - } - } -} - -/* ---------------------------------------------------------------------- - Mindlin force model -------------------------------------------------------------------------- */ - -TangentialMindlinForce::TangentialMindlinForce(LAMMPS *lmp) : TangentialMindlin(lmp) -{ - num_coeffs = 3; - size_history = 3; - mindlin_force = 1; - mindlin_rescale = 0; -} - -/* ---------------------------------------------------------------------- - Mindlin rescale model -------------------------------------------------------------------------- */ - -TangentialMindlinRescale::TangentialMindlinRescale(LAMMPS *lmp) : TangentialMindlin(lmp) -{ - num_coeffs = 3; - size_history = 4; - mindlin_force = 0; - mindlin_rescale = 1; - - nondefault_history_transfer = 1; - transfer_history_factor = new double[size_history]; - for (int i = 0; i < size_history; i++) transfer_history_factor[i] = -1.0; - transfer_history_factor[3] = +1; -} - -/* ---------------------------------------------------------------------- - Mindlin rescale force model -------------------------------------------------------------------------- */ - -TangentialMindlinRescaleForce::TangentialMindlinRescaleForce(LAMMPS *lmp) : TangentialMindlin(lmp) -{ - num_coeffs = 3; - size_history = 4; - mindlin_force = 1; - mindlin_rescale = 1; - - nondefault_history_transfer = 1; - transfer_history_factor = new double[size_history]; - for (int i = 0; i < size_history; i++) transfer_history_factor[i] = -1.0; - transfer_history_factor[3] = +1; -} diff --git a/src/contact_tangential_models.h b/src/contact_tangential_models.h deleted file mode 100644 index a54d03c495..0000000000 --- a/src/contact_tangential_models.h +++ /dev/null @@ -1,113 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef CONTACT_TANGENTIAL_MODELS_H_ -#define CONTACT_TANGENTIAL_MODELS_H_ - -#include "contact_sub_models.h" - -namespace LAMMPS_NS { -namespace Contact { - -class TangentialModel : public SubModel { - public: - TangentialModel(class LAMMPS *); - virtual ~TangentialModel() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; - virtual void calculate_forces() = 0; - int rescale_flag; - double k, damp, mu; // Used by Marshall twisting model - protected: - double xt; -}; - -/* ---------------------------------------------------------------------- */ - -class TangentialNone : public TangentialModel { - public: - TangentialNone(class LAMMPS *); - void calculate_forces() {}; -}; - -/* ---------------------------------------------------------------------- */ - -class TangentialLinearNoHistory : public TangentialModel { - public: - TangentialLinearNoHistory(class LAMMPS *); - void coeffs_to_local() override; - void calculate_forces(); -}; - -/* ---------------------------------------------------------------------- */ - -class TangentialLinearHistory : public TangentialModel { - public: - TangentialLinearHistory(class LAMMPS *); - void coeffs_to_local() override; - void calculate_forces(); -}; - -/* ---------------------------------------------------------------------- */ - -class TangentialLinearHistoryClassic : public TangentialLinearHistory { - public: - TangentialLinearHistoryClassic(class LAMMPS *); - void calculate_forces(); - int scale_area; -}; - -/* ---------------------------------------------------------------------- */ - -class TangentialMindlinClassic : public TangentialLinearHistoryClassic { - public: - TangentialMindlinClassic(class LAMMPS *); -}; - -/* ---------------------------------------------------------------------- */ - -class TangentialMindlin : public TangentialModel { - public: - TangentialMindlin(class LAMMPS *); - void coeffs_to_local() override; - void mix_coeffs(double*, double*) override; - void calculate_forces(); - protected: - int mindlin_rescale, mindlin_force; -}; - -/* ---------------------------------------------------------------------- */ - -class TangentialMindlinForce : public TangentialMindlin { - public: - TangentialMindlinForce(class LAMMPS *); -}; - -/* ---------------------------------------------------------------------- */ - -class TangentialMindlinRescale : public TangentialMindlin { - public: - TangentialMindlinRescale(class LAMMPS *); -}; - -/* ---------------------------------------------------------------------- */ - -class TangentialMindlinRescaleForce : public TangentialMindlin { - public: - TangentialMindlinRescaleForce(class LAMMPS *); -}; - -} // namespace Contact -} // namespace LAMMPS_NS - -#endif /*CONTACT_TANGENTIAL_MODELS_H_ */ diff --git a/src/contact_twisting_models.cpp b/src/contact_twisting_models.cpp deleted file mode 100644 index 94a88161f3..0000000000 --- a/src/contact_twisting_models.cpp +++ /dev/null @@ -1,124 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "contact_normal_models.h" -#include "contact_tangential_models.h" -#include "contact_twisting_models.h" -#include "contact.h" -#include "error.h" -#include "math_const.h" - -using namespace LAMMPS_NS; -using namespace Contact; -using namespace MathConst; - -/* ---------------------------------------------------------------------- - Default twisting model -------------------------------------------------------------------------- */ - -TwistingModel::TwistingModel(LAMMPS *lmp) : SubModel(lmp) {} - -/* ---------------------------------------------------------------------- - No model -------------------------------------------------------------------------- */ - -TwistingNone::TwistingNone(LAMMPS *lmp) : TwistingModel(lmp) {} - -/* ---------------------------------------------------------------------- - Marshall twisting model -------------------------------------------------------------------------- */ - -TwistingMarshall::TwistingMarshall(LAMMPS *lmp) : TwistingModel(lmp) -{ - num_coeffs = 0; - size_history = 3; -} - -/* ---------------------------------------------------------------------- */ - - -void TwistingMarshall::init() -{ - k_tang = contact->tangential_model->k; - mu_tang = contact->tangential_model->mu; -} - -/* ---------------------------------------------------------------------- */ - -void TwistingMarshall::calculate_forces() -{ - double signtwist, Mtcrit; - - // Calculate twist coefficients from tangential model & contact geometry - // eq 32 of Marshall paper - double k = 0.5 * k_tang * contact->area * contact->area; - double damp = 0.5 * contact->tangential_model->damp * contact->area * contact->area; - double mu = TWOTHIRDS * mu_tang * contact->area; - - if (contact->history_update) { - contact->history[history_index] += contact->magtwist * contact->dt; - } - - // M_t torque (eq 30) - contact->magtortwist = -k * contact->history[history_index] - damp * contact->magtwist; - signtwist = (contact->magtwist > 0) - (contact->magtwist < 0); - Mtcrit = mu * contact->normal_model->Fncrit; // critical torque (eq 44) - - if (fabs(contact->magtortwist) > Mtcrit) { - contact->history[history_index] = (Mtcrit * signtwist - damp * contact->magtwist) / k; - contact->magtortwist = -Mtcrit * signtwist; // eq 34 - } -} - -/* ---------------------------------------------------------------------- - SDS twisting model -------------------------------------------------------------------------- */ - -TwistingSDS::TwistingSDS(LAMMPS *lmp) : TwistingModel(lmp) -{ - num_coeffs = 3; - size_history = 3; -} - -/* ---------------------------------------------------------------------- */ - -void TwistingSDS::coeffs_to_local() -{ - k = coeffs[0]; - damp = coeffs[1]; - mu = coeffs[2]; - - if (k < 0.0 || mu < 0.0 || damp < 0.0) - error->all(FLERR, "Illegal SDS twisting model"); -} - -/* ---------------------------------------------------------------------- */ - -void TwistingSDS::calculate_forces() -{ - double signtwist, Mtcrit; - - if (contact->history_update) { - contact->history[history_index] += contact->magtwist * contact->dt; - } - - // M_t torque (eq 30) - contact->magtortwist = -k * contact->history[history_index] - damp * contact->magtwist; - signtwist = (contact->magtwist > 0) - (contact->magtwist < 0); - Mtcrit = mu * contact->normal_model->Fncrit; // critical torque (eq 44) - - if (fabs(contact->magtortwist) > Mtcrit) { - contact->history[history_index] = (Mtcrit * signtwist - damp * contact->magtwist) / k; - contact->magtortwist = -Mtcrit * signtwist; // eq 34 - } -} diff --git a/src/contact_twisting_models.h b/src/contact_twisting_models.h deleted file mode 100644 index 1fa5aa2c6a..0000000000 --- a/src/contact_twisting_models.h +++ /dev/null @@ -1,64 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef CONTACT_TWISTING_MODELS_H_ -#define CONTACT_TWISTING_MODELS_H_ - -#include "contact_sub_models.h" - -namespace LAMMPS_NS { -namespace Contact { - -class TwistingModel : public SubModel { - public: - TwistingModel(class LAMMPS *); - virtual ~TwistingModel() {}; - virtual void coeffs_to_local() {}; - virtual void init() {}; - virtual void calculate_forces() = 0; -}; - -/* ---------------------------------------------------------------------- */ - -class TwistingNone : public TwistingModel { - public: - TwistingNone(class LAMMPS *); - void calculate_forces() {}; -}; - -/* ---------------------------------------------------------------------- */ - -class TwistingMarshall : public TwistingModel { - public: - TwistingMarshall(class LAMMPS *); - void calculate_forces(); - void init(); - protected: - double k_tang, mu_tang; -}; - -/* ---------------------------------------------------------------------- */ - -class TwistingSDS : public TwistingModel { - public: - TwistingSDS(class LAMMPS *); - void coeffs_to_local() override; - void calculate_forces(); - protected: - double k, mu, damp; -}; - -} // namespace Contact -} // namespace LAMMPS_NS - -#endif /*CONTACT_TWISTING_MODELS_H_ */ diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index cc4fdbe20b..f5b86b6069 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -41,7 +41,7 @@ enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS, XSU,YSU,ZSU,XSUTRI,YSUTRI,ZSUTRI, IX,IY,IZ, VX,VY,VZ,FX,FY,FZ, - Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,heatflow,TEMPERATURE, + Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,HEATFLOW,TEMPERATURE, OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, TQX,TQY,TQZ, COMPUTE,FIX,VARIABLE,IVEC,DVEC,IARRAY,DARRAY}; @@ -930,7 +930,7 @@ int DumpCustom::count() for (i = 0; i < nlocal; i++) dchoose[i] = 2.0*radius[i]; ptr = dchoose; nstride = 1; - } else if (thresh_array[ithresh] == heatflow) { + } else if (thresh_array[ithresh] == HEATFLOW) { if (!atom->heatflow_flag) error->all(FLERR, "Threshold for an atom property that isn't allocated"); @@ -1867,7 +1867,7 @@ int DumpCustom::modify_param(int narg, char **arg) else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS; else if (strcmp(arg[1],"diameter") == 0) thresh_array[nthresh] = DIAMETER; - else if (strcmp(arg[1],"heatflow") == 0) thresh_array[nthresh] = heatflow; + else if (strcmp(arg[1],"heatflow") == 0) thresh_array[nthresh] = HEATFLOW; else if (strcmp(arg[1],"temperature") == 0) thresh_array[nthresh] = TEMPERATURE; else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX; else if (strcmp(arg[1],"omegay") == 0) thresh_array[nthresh] = OMEGAY;