From 22de863da999a45f7fd3d16c270dedf5337a2665 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 1 Aug 2022 13:12:44 -0600 Subject: [PATCH] Fixing compilation errors and finishing model classes --- src/GRANULAR/contact.cpp | 209 ++++++++++------ src/GRANULAR/contact.h | 61 +++-- src/GRANULAR/contact_damping_models.cpp | 40 ++-- src/GRANULAR/contact_damping_models.h | 40 ++-- src/GRANULAR/contact_heat_models.cpp | 184 ++------------ src/GRANULAR/contact_heat_models.h | 83 ++----- src/GRANULAR/contact_normal_models.cpp | 184 ++++++++------ src/GRANULAR/contact_normal_models.h | 68 +++--- src/GRANULAR/contact_rolling_models.cpp | 236 +++++------------- src/GRANULAR/contact_rolling_models.h | 79 ++---- src/GRANULAR/contact_sub_models.cpp | 55 +---- src/GRANULAR/contact_sub_models.h | 42 ++-- src/GRANULAR/contact_tangential_models.cpp | 130 +++++----- src/GRANULAR/contact_tangential_models.h | 92 ++++--- src/GRANULAR/contact_twisting_models.cpp | 265 ++++++--------------- src/GRANULAR/contact_twisting_models.h | 81 ++----- src/GRANULAR/pair_granular.cpp | 169 +++++++------ src/GRANULAR/pair_granular.h | 6 +- src/contact.h | 129 +++++----- src/fix_neigh_history.cpp | 6 +- src/pair.h | 2 +- 21 files changed, 903 insertions(+), 1258 deletions(-) diff --git a/src/GRANULAR/contact.cpp b/src/GRANULAR/contact.cpp index 12a9476585..4358f6df1c 100644 --- a/src/GRANULAR/contact.cpp +++ b/src/GRANULAR/contact.cpp @@ -18,99 +18,140 @@ */ #include "contact.h" +#include "pointers.h" #include "math_const.h" #include "math_extra.h" -#include "pointers.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 -using namespace LAMMPS_NS; using namespace MathExtra; -using namespace MathConst; +using namespace LAMMPS_NS::MathConst; +using namespace LAMMPS_NS::Contact; -namespace Contact { - -enum {NORMAL, TANGENTIAL, DAMPING, ROLLING, TWISTING, HEAT} - -ContactModel::ContactModel() +ContactModel::ContactModel() : + Pointers(lmp) { limit_damping = 0; beyond_contact = 0; + nondefault_history_transfer = 0; cutoff_type = 0.0; - normal_model = nullptr; - tangential_model = nullptr; - damping_model = nullptr; - rolling_model = nullptr; - twisting_model = nullptr; + nmodels = 5; + reset_contact(); - sub_models = {nullptr}; + for (int i = 0; i < nmodels; i++) sub_models[i] = nullptr; } /* ---------------------------------------------------------------------- */ -void ContactModel::init_model(char *model_name, int model_type) +ContactModel::~ContactModel() +{ + delete[] transfer_history_factor; +} + +/* ---------------------------------------------------------------------- */ + +void ContactModel::init_model(std::string model_name, ModelType model_type) { if (model_type == NORMAL) { - if (strcmp(model_name, "hooke") == 0) normal_model = new NormalHooke(); - else if (strcmp(model_name, "hertz") == 0) normal_model = new NormalHertz(); - else if (strcmp(model_name, "hertz/material") == 0) normal_model = new NormalHertzMaterial(); - else if (strcmp(model_name, "dmt") == 0) normal_model = new NormalDMT(); - else if (strcmp(model_name, "jkr") == 0) normal_model = new NormalJKR(); + if (model_name == "hooke") normal_model = new NormalHooke(); + else if (model_name == "hertz") normal_model = new NormalHertz(); + else if (model_name == "hertz/material") normal_model = new NormalHertzMaterial(); + else if (model_name == "dmt") normal_model = new NormalDMT(); + else if (model_name == "jkr") normal_model = new NormalJKR(); else error->all(FLERR, "Normal model name not recognized"); - sub_models[model_type] = &normal_model; + sub_models[model_type] = normal_model; } else if (model_type == TANGENTIAL) { - if (strcmp(model_name, "linear_nohistory") == 0) tangential_model = new TangentialLinearNoHistory(); - else if (strcmp(model_name, "linear_history") == 0) tangential_model = new TangentialLinearHistory(); - else if (strcmp(model_name, "mindlin") == 0) tangential_model = new TangentialMindlin(); - else if (strcmp(model_name, "mindlin/force") == 0) tangential_model = new TangentialMindlinForce(); - else if (strcmp(model_name, "mindlin_rescale") == 0) tangential_model = new TangentialMindlinRescale(); - else if (strcmp(model_name, "mindlin_rescale/force") == 0) tangential_model = new TangentialMindlinRescaleForce(); + if (model_name =="linear_nohistory") tangential_model = new TangentialLinearNoHistory(); + else if (model_name =="linear_history") tangential_model = new TangentialLinearHistory(); + else if (model_name =="mindlin") tangential_model = new TangentialMindlin(); + else if (model_name =="mindlin/force") tangential_model = new TangentialMindlinForce(); + else if (model_name =="mindlin_rescale") tangential_model = new TangentialMindlinRescale(); + else if (model_name =="mindlin_rescale/force") tangential_model = new TangentialMindlinRescaleForce(); else error->all(FLERR, "Tangential model name not recognized"); - sub_models[model_type] = &tangential_model; + sub_models[model_type] = tangential_model; } else if (model_type == DAMPING) { - if (strcmp(model_name, "velocity") == 0) damping_model = new DampingVelocity(); - else if (strcmp(model_name, "mass_velocity") == 0) damping_model = new DampingMassVelocity(); - else if (strcmp(model_name, "viscoelastic") == 0) damping_model = new DampingViscoelastic(); - else if (strcmp(model_name, "tsuji") == 0) damping_model = new DampingTsuji(); + if (model_name =="velocity") damping_model = new DampingVelocity(); + else if (model_name =="mass_velocity") damping_model = new DampingMassVelocity(); + else if (model_name =="viscoelastic") damping_model = new DampingViscoelastic(); + else if (model_name =="tsuji") damping_model = new DampingTsuji(); else error->all(FLERR, "Damping model name not recognized"); - sub_models[model_type] = &damping_model; + sub_models[model_type] = damping_model; } else if (model_type == ROLLING) { - if (strcmp(model_name, "none") == 0) delete rolling_model; - else if (strcmp(model_name, "sds") == 0) rolling_model = new RollingSDS(); + if (model_name =="none") delete rolling_model; + else if (model_name =="sds") rolling_model = new RollingSDS(); else error->all(FLERR, "Rolling model name not recognized"); - sub_models[model_type] = &rolling_model; + sub_models[model_type] = rolling_model; } else if (model_type == TWISTING) { - if (strcmp(model_name, "none") == 0) delete twisting_model; - else if (strcmp(model_name, "sds") == 0) twisting_model = new TwistingSDS(); - else if (strcmp(model_name, "marshall") == 0) twisting_model = new TwistingMarshall(); + if (model_name =="none") delete twisting_model; + else if (model_name =="sds") twisting_model = new TwistingSDS(); + else if (model_name =="marshall") twisting_model = new TwistingMarshall(); else error->all(FLERR, "Twisting model name not recognized"); - sub_models[model_type] = &twisting_model; + sub_models[model_type] = twisting_model; } else if (model_type == HEAT) { - if (strcmp(model_name, "none") == 0) delete heat_model; - else if (strcmp(model_name, "area") == 0) heat_model = new HeatArea(); + if (model_name =="none") delete heat_model; + else if (model_name =="area") heat_model = new HeatArea(); else error->all(FLERR, "Heat model name not recognized"); - sub_models[model_type] = &heat_model; + sub_models[model_type] = heat_model; } sub_models[model_type]->name.assign(model_name); - sub_models[model_type]->contact = *this; + sub_models[model_type]->contact = this; sub_models[model_type]->allocate_coeffs(); } /* ---------------------------------------------------------------------- */ +void ContactModel::init() +{ + int i, j, size_cumulative; + for (i = 0; i < nmodels; i++) { + if (sub_models[i]->nondefault_history_transfer) + nondefault_history_transfer = 1; + if (sub_models[i]->beyond_contact) + beyond_contact = 1; + } + + transfer_history_factor = new double(size_history); + for (i = 0; i < size_history; i++) transfer_history_factor[i] = -1.0; + + if (nondefault_history_transfer) { + for (i = 0; i < size_history; i++) { + // Find which model controls this history value + size_cumulative = 0; + for (j = 0; j < nmodels; 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 + if (sub_models[j]->nondefault_history_transfer) { + transfer_history_factor[i] = sub_models[j]->transfer_history_factor[i - size_cumulative]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + void ContactModel::mix_coeffs(ContactModel *c1, ContactModel *c2) { - for (int i = 0; i < 5; i++) - sum_model[i]->mix_coeffs(c1->submodel[i], c2->submodel[i]); + for (int i = 0; i < nmodels; i++) + sub_models[i]->mix_coeffs(c1->sub_models[i], c2->sub_models[i]); limit_damping = MAX(c1->limit_damping, c2->limit_damping); cutoff_type = MAX(c1->cutoff_type, c2->cutoff_type); @@ -120,11 +161,17 @@ void ContactModel::mix_coeffs(ContactModel *c1, ContactModel *c2) void ContactModel::write_restart(FILE *fp) { - int num_char = -1; + int num_char, num_coeffs; - for (int i = 0; i < 5; i++) { + for (int i = 0; i < nmodels; i++) { + num_char = -1; if (sub_models[i]) { - sub_models[i]->write_restart(fp); + 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(int), num_coeffs, fp); } else { fwrite(&num_char, sizeof(int), 1, fp); } @@ -135,21 +182,34 @@ void ContactModel::write_restart(FILE *fp) void ContactModel::read_restart(FILE *fp) { - int num_char; + int num_char, num_coeff; - for (int i = 0; i < 5; i++) { + for (int i = 0; i < nmodels; 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); + utils::sfread(FLERR, &num_char, sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&num_char, 1, MPI_INT, 0, world); if (num_char != -1) { - std::string model_name(num_char); + 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, world); + 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(const_cast(model_name.data(), i)); - sub_models[i]->read_restart(); + 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(int), num_coeff, fp, nullptr, error); + } + MPI_Bcast(sub_models[i]->coeffs, num_coeff, MPI_DOUBLE, 0, world); + + sub_models[i]->coeffs_to_local(); } } } @@ -158,7 +218,7 @@ void ContactModel::read_restart(FILE *fp) void ContactModel::reset_contact() { - radi = radj = Fntot = Fncrit = magtortwist = 0.0; + radi = radj = Fntot = magtortwist = 0.0; xi = xj = vi = vj = omegai = omegaj = nullptr; forces = torquesi = torquesj = history = nullptr; @@ -177,10 +237,7 @@ bool ContactModel::check_contact() radsum = radi + radj; Reff = radi * radj / radsum; - touch = false; - if (normal_model.beyond_contact) normal_model.touch(touch); - else touch = (rsq < radsum*radsum); - + touch = normal_model->touch(); return touch; } @@ -220,7 +277,7 @@ void ContactModel::prep_contact() sub3(vt, temp, vtr); vrel = len3(vtr); - if (roll_model || twist_model) + if (rolling_model || twisting_model) sub3(omegai, omegaj, relrot); if (rolling_model) { @@ -250,20 +307,24 @@ void ContactModel::calculate_forces() } //********************************************** - // normal forces + // calculate forces //********************************************** double Fne, Fdamp; - Fne = normal_model.calculate_forces(); - Fdamp = damping_model.calculate_forces(); + area = normal_model->calculate_area(); + Fne = normal_model->calculate_forces(); + normal_model->set_knfac(); // Needed for damping + + Fdamp = damping_model->calculate_forces(); + + normal_model->set_fncrit(); // Needed for tangential, rolling, twisting Fntot = Fne + Fdamp; if (limit_damping && Fntot < 0.0) Fntot = 0.0; - normal_model.set_fncrit(); - tangential_model.calculate_forces(); - if (rolling_model) rolling_model.calculate_forces(); - if (twisting_model) twisting_model.calculate_forces(); + tangential_model->calculate_forces(); + if (rolling_model) rolling_model->calculate_forces(); + if (twisting_model) twisting_model->calculate_forces(); //********************************************** // sum contributions @@ -301,7 +362,7 @@ void ContactModel::calculate_forces() double ContactModel::calculate_heat() { - return heat_model.calculate_heat(); + return heat_model->calculate_heat(); } /* ---------------------------------------------------------------------- @@ -311,7 +372,5 @@ double ContactModel::calculate_heat() double ContactModel::pulloff_distance(double radi, double radj) { - return normal_model.pulloff_distance(radi, radj); -} - + return normal_model->pulloff_distance(radi, radj); } diff --git a/src/GRANULAR/contact.h b/src/GRANULAR/contact.h index f9f18a1eb4..96758771b4 100644 --- a/src/GRANULAR/contact.h +++ b/src/GRANULAR/contact.h @@ -14,32 +14,44 @@ #ifndef LMP_CONTACT_H #define LMP_CONTACT_H -#include "pointers.h" -#include "normal_contact_models.h" -#include "tangential_contact_models.h" -#include "damping_contact_models.h" -#include "rolling_contact_models.h" -#include "twisting_contact_models.h" -#include "heat_models.h" - -using namespace LAMMPS_NS; +#include "pointers.h" // IWYU pragma: export +namespace LAMMPS_NS { namespace Contact { +enum ModelType { + NORMAL = 0, + DAMPING = 1, + TANGENTIAL = 2, + ROLLING = 3, + TWISTING = 4, + HEAT = 5 +}; + #define EPSILON 1e-10 +// forward declaration +class NormalModel; +class DampingModel; +class TangentialModel; +class RollingModel; +class TwistingModel; +class HeatModel; +class SubModel; + class ContactModel : protected Pointers { -public: + public: ContactModel(); ~ContactModel(); - int init(); + void init(); bool check_contact(); void reset_contact(); + void prep_contact(); void calculate_forces(); double calculate_heat(); double pulloff_distance(double, double); - void init_model(char*, int); + void init_model(std::string, ModelType); void mix_coeffs(ContactModel*, ContactModel*); @@ -47,21 +59,25 @@ public: void read_restart(FILE *); // Sub models - NormalModel normal_model; - DampingModel damping_model; - TangentialModel tangential_model; - RollingModel rolling_model; - TwistingModel twisting_model; - HeatModel heat_model; + NormalModel *normal_model; + DampingModel *damping_model; + TangentialModel *tangential_model; + RollingModel *rolling_model; + TwistingModel *twisting_model; + HeatModel *heat_model; SubModel *sub_models[6]; // Need to resize if we add more model flavors // Extra options - int beyond_contact, limit_damping; + int beyond_contact, limit_damping, history_update; double cutoff_type; + // History variables + int size_history, nondefault_history_transfer; + double *transfer_history_factor; + double *history; + // Contact properties/output double *forces, *torquesi, *torquesj; - double *history; double radi, radj, meff, dt, Ti, Tj, area; double Fntot, magtortwist; @@ -69,15 +85,18 @@ public: double *xi, *xj, *vi, *vj, *omegai, *omegaj; double fs[3], fr[3], ft[3]; -private: 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 prep_flag, check_flag; + int nmodels; }; } // namespace Contact +} // namespace LAMMPS_NS #endif diff --git a/src/GRANULAR/contact_damping_models.cpp b/src/GRANULAR/contact_damping_models.cpp index 1e233234b2..7a79c6946c 100644 --- a/src/GRANULAR/contact_damping_models.cpp +++ b/src/GRANULAR/contact_damping_models.cpp @@ -11,33 +11,22 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "damping_contact_models.h" -#include "math_const.h" +#include "contact_damping_models.h" +#include "contact_normal_models.h" #include "contact.h" - -#include #include "math_special.h" -using namespace MathConst; +using namespace LAMMPS_NS; +using namespace Contact; using namespace MathSpecial; -namespace Contact{ - - /* ---------------------------------------------------------------------- Default damping model ------------------------------------------------------------------------- */ -void DampingModel::DampingModel() -{ - num_coeffs = 0; -} - -/* ---------------------------------------------------------------------- */ - void DampingModel::coeffs_to_local() { - damp = contact.normal_model.damp; + damp = contact->normal_model->damp; } /* ---------------------------------------------------------------------- @@ -46,34 +35,34 @@ void DampingModel::coeffs_to_local() double DampingVelocity::calculate_forces() { - return -damp * contact.vnnr; + return -damp * contact->vnnr; } /* ---------------------------------------------------------------------- Mass velocity damping ------------------------------------------------------------------------- */ -double MassVelocity::calculate_forces() +double DampingMassVelocity::calculate_forces() { - return -damp * contact.meff * contact.vnnr; + return -damp * contact->meff * contact->vnnr; } /* ---------------------------------------------------------------------- Default, viscoelastic damping ------------------------------------------------------------------------- */ -double ViscoElastic::calculate_forces() +double DampingViscoelastic::calculate_forces() { - return -damp * contact.meff * contact.area * contact.vnnr; + return -damp * contact->meff * contact->area * contact->vnnr; } /* ---------------------------------------------------------------------- Tsuji damping ------------------------------------------------------------------------- */ -void Tsuji::coeffs_to_local() +void DampingTsuji::coeffs_to_local() { - double tmp = contact.normal_model.damp; + 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); @@ -81,8 +70,7 @@ void Tsuji::coeffs_to_local() /* ---------------------------------------------------------------------- */ -double Tsuji::calculate_forces() +double DampingTsuji::calculate_forces() { - return -damp_scaled * sqrt(contact.meff * contact.normal_model.knfac) * contact.vnnr; + return -damp * sqrt(contact->meff * contact->normal_model->knfac) * contact->vnnr; } - diff --git a/src/GRANULAR/contact_damping_models.h b/src/GRANULAR/contact_damping_models.h index 0fc4fd2c17..440ab36369 100644 --- a/src/GRANULAR/contact_damping_models.h +++ b/src/GRANULAR/contact_damping_models.h @@ -11,60 +11,54 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef DAMPING_CONTACT_MODELS_H_ -#define DAMPING_CONTACT_MODELS_H_ +#ifndef CONTACT_DAMPING_MODELS_H_ +#define CONTACT_DAMPING_MODELS_H_ -#include "contact.h"; -#include "sub_model.h" +#include "contact_sub_models.h" +namespace LAMMPS_NS { namespace Contact { -class DampingModel:SubModel -{ +class DampingModel : public SubModel { public: - DampingModel(); - virtual ~DampingModel() {}; - virtual double calculate_forces() = 0; + DampingModel() {}; + ~DampingModel() {}; virtual void coeffs_to_local(); - virtual void mix_coeffs(NormalModel*, NormalModel*); + virtual void mix_coeffs(DampingModel*, DampingModel*) {}; + virtual double calculate_forces() = 0; double damp; }; /* ---------------------------------------------------------------------- */ -class DampingVelocity:DampingModel -{ +class DampingVelocity: public DampingModel { public: double calculate_forces(); }; /* ---------------------------------------------------------------------- */ -class DampingMassVelocity:DampingModel -{ +class DampingMassVelocity: public DampingModel { public: double calculate_forces(); }; /* ---------------------------------------------------------------------- */ -class DampingViscoElastic:DampingModel -{ +class DampingViscoelastic: public DampingModel { public: double calculate_forces(); }; /* ---------------------------------------------------------------------- */ -class DampingTsuji:DampingModel -{ +class DampingTsuji: public DampingModel { public: - double calculate_forces(); void coeffs_to_local(); + double calculate_forces(); }; +} // namespace Contact +} // namespace LAMMPS_NS -} - -#endif /*DAMPING_CONTACT_MODELS_H_ */ - +#endif /*CONTACT_DAMPING_MODELS_H_ */ diff --git a/src/GRANULAR/contact_heat_models.cpp b/src/GRANULAR/contact_heat_models.cpp index e6bafe4b3c..29134b3e73 100644 --- a/src/GRANULAR/contact_heat_models.cpp +++ b/src/GRANULAR/contact_heat_models.cpp @@ -11,175 +11,39 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "normal_contact_models.h" -#include "math_const.h" +#include "contact_heat_models.h" #include "contact.h" -#include +using namespace LAMMPS_NS; +using namespace Contact; -using namespace MathConst; +/* ---------------------------------------------------------------------- + Area-based heat conduction +------------------------------------------------------------------------- */ -namespace Contact{ - -#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 behaviors where needed -// ************************ -void NormalModel::set_fncrit(){ - contact->Fncrit = fabs(contact->Fntot); +HeatArea::HeatArea() +{ + num_coeffs = 1; } -void NormalModel::mix_coeffs(NormalModel* imodel, NormalModel* jmodel){ - for (int i = 0; i < num_coeffs; i++){ - coeffs[i] = sqrt(imodel->coeffs[i]*jmodel->coeffs[i]); - } +/* ---------------------------------------------------------------------- */ + +void HeatArea::coeffs_to_local() +{ + conductivity = coeffs[0]; } -//----------------------------------------- +/* ---------------------------------------------------------------------- */ -//****************** -// Hooke -//****************** -void Hooke::Hooke(ContactModel &c){ - contact = c; - num_coeffs = 2; - allocate_coeffs(); +void HeatArea::mix_coeffs(HeatModel* imodel, HeatModel* jmodel) +{ + coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]); + coeffs_to_local(); } -void Hooke::coeffs_to_local(){ - k_norm = coeffs[0]; - damp = coeffs[1]; +/* ---------------------------------------------------------------------- */ + +double HeatArea::calculate_heat() +{ + return conductivity * contact->area * (contact->Ti - contact->Tj); } - -double Hooke::calculate_forces(){ - contact.area = sqrt(contact.dR); - contact.knfac = k_norm * contact.area; - return contact.knfac * contact.delta; -} - - -//****************** -// Hertz -//****************** -void Hertz::Hertz(ContactModel &c, int mat_flag){ - contact = c; - material_prop_flag = mat_flag; - if (material_prop_flag){ - num_coeffs = 3; - } - else{ - num_coeffs = 2; - } - allocate_coeffs(); -} - -void Hertz::coeffs_to_local(){ - if (material_prop_flag){ - Emod = coeffs[0]; - poiss = coeffs[1]; - k_norm = 4/3*Emod; - } - else{ - k_norm = coeffs[0]; - } -} - -double Hertz::calculate_forces(){ - contact.area = sqrt(contact.dR); - contact.knfac = contact.k_norm * contact.area; - return contact.knfac * contact.delta; -} - -//****************** -// DMT -//****************** -void DMT::DMT(ContactModel &c){ - contact = c; - material_prop_flag = 1; - num_coeffs = 4; - allocate_coeffs(); -} - -double DMT::calculate_forces(){ - contact.area = sqrt(contact.dR); - contact.knfac = contact.k_norm * contact.area; - double Fne = contact.knfac * contact.delta; - F_pulloff = 4 * MathConst::MY_PI * cohesion * contact.Reff; - Fne -= F_pulloff; - return Fne; -} - -void DMT::set_fncrit(){ - contact.Fncrit = fabs(contact.Fne + 2* F_pulloff); -} - -//****************** -// JKR -//****************** -void JKR::JKR(ContactModel &c){ - contact = c; - material_prop_flag = 1; - beyond_contact = 1; - num_coeffs = 4; - allocate_coeffs(); -} - -bool JKR::touch(int touch){ - double Escaled, R2, delta_pulloff, dist_pulloff; - bool touchflag; - - Escaled = k_norm * THREEQUARTERS; - if (touch) { - R2 = contact.Reff * contact.Reff; - a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled)); - delta_pulloff = a * a / contact.Reff - 2 * sqrt(MY_PI * cohesion * a / Escaled); - dist_pulloff = contact.radsum - delta_pulloff; - touchflag = (contact.rsq < dist_pulloff * dist_pulloff); - } else { - touchflag = (rcontact.sq < contact.radsum * contact.radsum); - } - return touchflag; -} - -double JKR::calculate_forces(){ - double Escaled, R2, dR2, t0, t1, t2, t3, t4, t5, t6; - double sqrt1, sqrt2, sqrt3, a2, F_pulloff, Fne; - - Escaled = k_norm * THREEQUARTERS; - - R2 = Reff * Reff; - dR2 = dR * dR; - t0 = cohesion * cohesion * R2 * R2 * Escaled; - t1 = PI27SQ*t0; - t2 = 8 * dR * dR2 * Escaled * Escaled * Escaled; - t3 = 4 * dR2 * Escaled; - - // in case sqrt(0) < 0 due to precision issues - sqrt1 = MAX(0, t0 * (t1 + 2 * t2)); - t4 = cbrt(t1 + t2 + THREEROOT3 * MY_PI * sqrt(sqrt1)); - t5 = t3 / t4 + t4 / Escaled; - sqrt2 = MAX(0, 2 * dR + t5); - t6 = sqrt(sqrt2); - sqrt3 = MAX(0, 4 * dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Escaled * t6)); - a = INVROOT6 * (t6 + sqrt(sqrt3)); - a2 = a * a; - Fne = Escaled * a * a2 / Reff - MY_2PI * a2 * sqrt(4 * cohesion * Escaled / (MY_PI * a)); - F_pulloff = 3 * MY_PI * cohesion * Reff; - - knfac = Escaled * a; - return Fne; -} - -void JKR::set_fncrit(){ - contact.Fncrit = fabs(contact.Fne + 2 * F_pulloff); -} - - - diff --git a/src/GRANULAR/contact_heat_models.h b/src/GRANULAR/contact_heat_models.h index 1a2f1ddb7d..469a53cb13 100644 --- a/src/GRANULAR/contact_heat_models.h +++ b/src/GRANULAR/contact_heat_models.h @@ -11,73 +11,36 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef NORMAL_CONTACT_MODELS_H_ -#define NORMAL_CONTACT_MODELS_H_ +#ifndef CONTACT_HEAT_MODELS_H_ +#define CONTACT_HEAT_MODELS_H_ -#include "contact.h"; -#include "sub_model.h" +#include "contact_sub_models.h" +namespace LAMMPS_NS { namespace Contact { -class NormalModel:SubModel{ -public: - NormalModel(){}; - virtual ~NormalModel(){}; - virtual bool check_contact(); - virtual void prep_contact(); - virtual void set_fncrit(); - virtual double calculate_forces() = 0; - virtual void coeffs_to_local(); - virtual void mix_coeffs(NormalModel*, NormalModel*); //When mixing is needed - -private: - int beyond_contact = 0; - int allow_limit_damping = 1; - +class HeatModel : public SubModel { + public: + HeatModel() {}; + ~HeatModel() {}; + virtual void coeffs_to_local() {}; + virtual void mix_coeffs(HeatModel*, HeatModel*) {}; + virtual double calculate_heat() = 0; }; -class Hooke:NormalModel{ -public: - Hooke(ContactModel &c); - ~Hooke(){}; +/* ---------------------------------------------------------------------- */ + +class HeatArea: public HeatModel { + public: + HeatArea(); void coeffs_to_local(); - double calculate_forces(); -private: - double k_norm, damp; + void mix_coeffs(HeatModel*, HeatModel*); + double calculate_heat(); + private: + double conductivity; }; -class Hertz:NormalModel{ -public: - Hertz(ContactModel&, int); - ~Hertz(){}; - void coeffs_to_local(); - double calculate_forces(); -private: - double k_norm, damp, Emod, poiss; -}; - -class DMT:NormalModel{ -public: - DMT(ContactModel &c); - ~DMT(){}; - void coeffs_to_local(); - void coeffs_to_local(NormalModel*, NormalModel*); - double calculate_forces(); -private: - double k_norm, damp, Emod, poiss, coh; -}; - -class JKR:NormalModel{ -public: - JKR(ContactModel &c); - ~JKR(){}; - void coeffs_to_local(); - double calculate_forces(); -private: - double k_norm, damp, Emod, poiss, coh; - -}; -} - -#endif /*NORMAL_CONTACT_MODELS_H_ */ +} // namespace Contact +} // namespace LAMMPS_NS +#endif /*CONTACT_HEAT_MODELS_H_ */ diff --git a/src/GRANULAR/contact_normal_models.cpp b/src/GRANULAR/contact_normal_models.cpp index 9b1d85d123..c11fddb1d7 100644 --- a/src/GRANULAR/contact_normal_models.cpp +++ b/src/GRANULAR/contact_normal_models.cpp @@ -12,15 +12,13 @@ ------------------------------------------------------------------------- */ #include "contact_normal_models.h" -#include "math_const.h" #include "contact.h" +#include "math_const.h" -#include - +using namespace LAMMPS_NS; +using namespace Contact; using namespace MathConst; -namespace Contact{ - #define PI27SQ 266.47931882941264802866 // 27*PI**2 #define THREEROOT3 5.19615242270663202362 // 3*sqrt(3) #define SIXROOT6 14.69693845669906728801 // 6*sqrt(6) @@ -40,23 +38,40 @@ NormalModel::NormalModel() /* ---------------------------------------------------------------------- */ -void NormalModel::set_fncrit() +bool NormalModel::touch() { - Fncrit = fabs(contact.Fntot); + bool touchflag = (contact->rsq < contact->radsum * contact->radsum); + return touchflag; +} + +/* ---------------------------------------------------------------------- + called outside of compute(), do not assume geometry defined in contact +------------------------------------------------------------------------- */ + +double NormalModel::pulloff_distance(double radi, double radj) +{ + return radi + radj; } /* ---------------------------------------------------------------------- */ -void NormalModel::pulloff_distance(double radi, double radj) +double NormalModel::calculate_area() { - return radi + radj; + return sqrt(contact->dR); +} + +/* ---------------------------------------------------------------------- */ + +void NormalModel::set_fncrit() +{ + Fncrit = fabs(contact->Fntot); } /* ---------------------------------------------------------------------- Hookean normal force ------------------------------------------------------------------------- */ -void NormalHooke::NormalHooke() +NormalHooke::NormalHooke() { num_coeffs = 2; allocate_coeffs(); @@ -66,7 +81,7 @@ void NormalHooke::NormalHooke() void NormalHooke::coeffs_to_local() { - k_norm = coeffs[0]; + k = coeffs[0]; damp = coeffs[1]; } @@ -83,17 +98,22 @@ void NormalHooke::mix_coeffs(NormalModel* imodel, NormalModel* jmodel) double NormalHooke::calculate_forces() { - contact.area = sqrt(contact.dR); - knfac = k_norm * contact.area; - Fne = knfac * contact.delta; + Fne = knfac * contact->delta; return Fne; } +/* ---------------------------------------------------------------------- */ + +void NormalHooke::set_knfac() +{ + knfac = k * contact->area; +} + /* ---------------------------------------------------------------------- Hertzian normal force ------------------------------------------------------------------------- */ -void NormalHertz::NormalHertz() +NormalHertz::NormalHertz() { num_coeffs = 2; } @@ -102,7 +122,7 @@ void NormalHertz::NormalHertz() void NormalHertz::coeffs_to_local() { - k_norm = coeffs[0]; + k = coeffs[0]; damp = coeffs[1]; } @@ -119,17 +139,23 @@ void NormalHertz::mix_coeffs(NormalModel* imodel, NormalModel* jmodel) double NormalHertz::calculate_forces() { - contact.area = sqrt(contact.dR); - knfac = contact.k_norm * contact.area; - Fne = knfac * contact.delta; + contact->area = sqrt(contact->dR); + Fne = knfac * contact->delta; return Fne; } +/* ---------------------------------------------------------------------- */ + +void NormalHertz::set_knfac() +{ + knfac = k * contact->area; +} + /* ---------------------------------------------------------------------- Hertzian normal force with material properties ------------------------------------------------------------------------- */ -void NormalHertzMaterial::NormalHertzMaterial() +NormalHertzMaterial::NormalHertzMaterial() { material_properties = 1; num_coeffs = 3; @@ -142,7 +168,7 @@ void NormalHertzMaterial::coeffs_to_local() Emod = coeffs[0]; damp = coeffs[1]; poiss = coeffs[2]; - k_norm = 4 / 3 * Emod; + k = 4 / 3 * Emod; } /* ---------------------------------------------------------------------- */ @@ -159,7 +185,7 @@ void NormalHertzMaterial::mix_coeffs(NormalModel* imodel, NormalModel* jmodel) DMT normal force ------------------------------------------------------------------------- */ -void NormalDMT::NormalDMT(ContactModel &c) +NormalDMT::NormalDMT() { allow_limit_damping = 0; material_properties = 1; @@ -174,7 +200,7 @@ void NormalDMT::coeffs_to_local() damp = coeffs[1]; poiss = coeffs[2]; cohesion = coeffs[3]; - k_norm = 4 / 3 * Emod; + k = 4 / 3 * Emod; } /* ---------------------------------------------------------------------- */ @@ -192,30 +218,35 @@ void NormalDMT::mix_coeffs(NormalModel* imodel, NormalModel* jmodel) double NormalDMT::calculate_forces() { - contact.area = sqrt(contact.dR); - knfac = k_norm * contact.area; - Fne = knfac * contact.delta; - F_pulloff = 4 * MathConst::MY_PI * cohesion * contact.Reff; + Fne = knfac * contact->delta; + F_pulloff = 4 * 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* F_pulloff); + Fncrit = fabs(Fne + 2 * F_pulloff); } /* ---------------------------------------------------------------------- JKR normal force ------------------------------------------------------------------------- */ -void NormalJKR::NormalJKR(ContactModel &c) +NormalJKR::NormalJKR() { allow_limit_damping = 0; material_properties = 1; - contact.beyond_contact = beyond_contact = 1; + beyond_contact = 1; num_coeffs = 4; } @@ -227,8 +258,8 @@ void NormalJKR::coeffs_to_local() damp = coeffs[1]; poiss = coeffs[2]; cohesion = coeffs[3]; - k_norm = 4/3*Emod; - Escaled = k_norm * THREEQUARTERS; + k = 4/3*Emod; + Escaled = k * THREEQUARTERS; } /* ---------------------------------------------------------------------- */ @@ -244,70 +275,83 @@ void NormalJKR::mix_coeffs(NormalModel* imodel, NormalModel* jmodel) /* ---------------------------------------------------------------------- */ -bool NormalJKR::touch(int touch) +bool NormalJKR::touch() { - double R2, delta_pulloff, dist_pulloff; + double area_at_pulloff, R2, delta_pulloff, dist_pulloff; bool touchflag; - if (touch) { - R2 = contact.Reff * contact.Reff; - a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled)); - delta_pulloff = a * a / contact.Reff - 2 * sqrt(MY_PI * cohesion * a / Escaled); - dist_pulloff = contact.radsum - delta_pulloff; - touchflag = (contact.rsq < dist_pulloff * dist_pulloff); - } else { - touchflag = (rcontact.sq < contact.radsum * contact.radsum); - } + R2 = contact->Reff * contact->Reff; + area_at_pulloff = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled)); + delta_pulloff = area_at_pulloff * area_at_pulloff / contact->Reff - 2 * sqrt(MY_PI * cohesion * area_at_pulloff /Escaled); + dist_pulloff = contact->radsum - delta_pulloff; + touchflag = (contact->rsq < dist_pulloff * dist_pulloff); + 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 * MY_PI * cohesion * Reff_tmp * Reff_tmp / (4 * Escaled)); + return area_at_pulloff * area_at_pulloff / Reff_tmp - 2 * sqrt(MY_PI * cohesion * area_at_pulloff / Escaled); +} + /* ---------------------------------------------------------------------- */ -double NormalJKR::calculate_forces() +double NormalJKR::calculate_area() { double R2, dR2, t0, t1, t2, t3, t4, t5, t6; - double sqrt1, sqrt2, sqrt3, a2; + double sqrt1, sqrt2, sqrt3; - R2 = Reff * Reff; - dR2 = dR * dR; + R2 = contact->Reff * contact->Reff; + dR2 = contact->dR * contact->dR; t0 = cohesion * cohesion * R2 * R2 * Escaled; - t1 = PI27SQ*t0; - t2 = 8 * dR * dR2 * Escaled * Escaled * Escaled; + t1 = PI27SQ * t0; + t2 = 8 * contact->dR * dR2 * Escaled * Escaled * Escaled; t3 = 4 * dR2 * Escaled; // in case sqrt(0) < 0 due to precision issues sqrt1 = MAX(0, t0 * (t1 + 2 * t2)); t4 = cbrt(t1 + t2 + THREEROOT3 * MY_PI * sqrt(sqrt1)); t5 = t3 / t4 + t4 / Escaled; - sqrt2 = MAX(0, 2 * dR + t5); + sqrt2 = MAX(0, 2 * contact->dR + t5); t6 = sqrt(sqrt2); - sqrt3 = MAX(0, 4 * dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Escaled * t6)); - contact.area = INVROOT6 * (t6 + sqrt(sqrt3)); - a2 = contact.area * contact.area; - Fne = Escaled * contact.area * a2 / Reff - MY_2PI * a2 * sqrt(4 * cohesion * Escaled / (MY_PI * contact.area)); - F_pulloff = 3 * MY_PI * cohesion * Reff; + sqrt3 = MAX(0, 4 * contact->dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Escaled * t6)); + + return INVROOT6 * (t6 + sqrt(sqrt3)); +} + +/* ---------------------------------------------------------------------- */ + +double NormalJKR::calculate_forces() +{ + double a2; + + a2 = contact->area * contact->area; + Fne = Escaled * contact->area * a2 / contact->Reff - MY_2PI * a2 * sqrt(4 * cohesion * Escaled / (MY_PI * contact->area)); + F_pulloff = 3 * MY_PI * cohesion * contact->Reff; - knfac = Escaled * contact.area; return Fne; } /* ---------------------------------------------------------------------- */ +void NormalJKR::set_knfac() +{ + knfac = Escaled * contact->area; +} + +/* ---------------------------------------------------------------------- */ + void NormalJKR::set_fncrit() { Fncrit = fabs(Fne + 2 * F_pulloff); } - -/* ---------------------------------------------------------------------- */ - -void NormalJKR::pulloff_distance(double radi, double radj) -{ - double a_tmp, Reff_tmp; - - Reff_tmp = radi * radj / (radi + radj); - if (Reff_tmp <= 0) return 0; - - a_tmp = cbrt(9 * MY_PI * cohesion * Reff_tmp * Reff_tmp / (4 * Ecaled)); - return a_tmp * a_tmp / Reff_tmp - 2 * sqrt(MY_PI * cohesion * a_tmp / Ecaled); -} - diff --git a/src/GRANULAR/contact_normal_models.h b/src/GRANULAR/contact_normal_models.h index e91d5af1c9..74e2800eaa 100644 --- a/src/GRANULAR/contact_normal_models.h +++ b/src/GRANULAR/contact_normal_models.h @@ -11,102 +11,108 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef NORMAL_CONTACT_MODELS_H_ -#define NORMAL_CONTACT_MODELS_H_ +#ifndef CONTACT_NORMAL_MODELS_H_ +#define CONTACT_NORMAL_MODELS_H_ -#include "contact.h"; -#include "sub_model.h" +#include "contact_sub_models.h" +namespace LAMMPS_NS { namespace Contact { -class NormalModel:SubModel{ +class NormalModel : public SubModel { public: NormalModel(); - virtual ~NormalModel() {}; - virtual bool check_contact(); - virtual void set_fncrit(); + ~NormalModel() {}; + virtual void coeffs_to_local() {}; + virtual void mix_coeffs(NormalModel*, NormalModel*) {}; + virtual bool touch(); + virtual double pulloff_distance(double, double); + virtual double calculate_area(); virtual double calculate_forces() = 0; - virtual void coeffs_to_local(); - virtual void mix_coeffs(NormalModel*, NormalModel*); - void pulloff_distance(double, double); + virtual void set_knfac() = 0; + virtual void set_fncrit(); double damp; // Vestigial argument needed by damping + double Emod, poiss; double Fncrit, Fne, knfac; int material_properties; }; /* ---------------------------------------------------------------------- */ -class NormalHooke:NormalModel -{ +class NormalHooke: public NormalModel { public: NormalHooke(); ~NormalHooke() {}; void coeffs_to_local(); void mix_coeffs(NormalModel*, NormalModel*); double calculate_forces(); - double Emod, poiss; + void set_knfac(); + private: + double k; }; /* ---------------------------------------------------------------------- */ -class NormalHertz:NormalModel -{ +class NormalHertz: public NormalModel { public: NormalHertz(); ~NormalHertz() {}; void coeffs_to_local(); void mix_coeffs(NormalModel*, NormalModel*); double calculate_forces(); + void set_knfac(); private: - double k_norm; + double k; }; /* ---------------------------------------------------------------------- */ -class NormalHertzMaterial:NormalHertz -{ +class NormalHertzMaterial: public NormalHertz { public: NormalHertzMaterial(); ~NormalHertzMaterial() {}; void coeffs_to_local(); void mix_coeffs(NormalModel*, NormalModel*); private: - double k_norm; + double k; }; /* ---------------------------------------------------------------------- */ -class NormalDMT:NormalModel -{ +class NormalDMT: public NormalModel { public: NormalDMT(); ~NormalDMT() {}; void coeffs_to_local(); void mix_coeffs(NormalModel*, NormalModel*); - void set_fncrit(); double calculate_forces(); + void set_knfac(); + void set_fncrit(); private: - double k_norm, cohesion; + double k, cohesion; double F_pulloff; }; /* ---------------------------------------------------------------------- */ -class NormalJKR:NormalModel -{ +class NormalJKR: public NormalModel { public: NormalJKR(); ~NormalJKR() {}; void coeffs_to_local(); void mix_coeffs(NormalModel*, NormalModel*); - void set_fncrit(); + bool touch(); + double pulloff_distance(double, double); + double calculate_area(); double calculate_forces(); - void pulloff_distance(double, double); + void set_knfac(); + void set_fncrit(); private: - double k_norm, cohesion; + double k, cohesion; double Escaled, F_pulloff; }; -} -#endif /*NORMAL_CONTACT_MODELS_H_ */ +} // namespace Contact +} // namespace LAMMPS_NS +#endif /*CONTACT_NORMAL_MODELS_H_ */ diff --git a/src/GRANULAR/contact_rolling_models.cpp b/src/GRANULAR/contact_rolling_models.cpp index 5521e10aee..0c256b0fbf 100644 --- a/src/GRANULAR/contact_rolling_models.cpp +++ b/src/GRANULAR/contact_rolling_models.cpp @@ -11,202 +11,71 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "normal_contact_models.h" +#include "contact_normal_models.h" +#include "contact_rolling_models.h" #include "math_const.h" +#include "math_extra.h" #include "contact.h" -#include - +using namespace LAMMPS_NS; +using namespace Contact; using namespace MathConst; +using namespace MathExtra; -namespace Contact{ +/* ---------------------------------------------------------------------- + SDS rolling friction model +------------------------------------------------------------------------- */ -#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 behaviors where needed -// ************************ -void NormalModel::set_fncrit(){ - contact->Fncrit = fabs(contact->Fntot); +RollingSDS::RollingSDS() +{ + num_coeffs = 3; + size_history = 3; } -void NormalModel::mix_coeffs(NormalModel* imodel, NormalModel* jmodel){ - for (int i = 0; i < num_coeffs; i++){ - coeffs[i] = sqrt(imodel->coeffs[i]*jmodel->coeffs[i]); - } +/* ---------------------------------------------------------------------- */ + +void RollingSDS::coeffs_to_local() +{ + k = coeffs[0]; + mu = coeffs[1]; + gamma = coeffs[2]; } -//----------------------------------------- +/* ---------------------------------------------------------------------- */ -//****************** -// Hooke -//****************** -void Hooke::Hooke(ContactModel &c){ - contact = c; - num_coeffs = 2; - allocate_coeffs(); +void RollingSDS::mix_coeffs(RollingModel* imodel, RollingModel* jmodel) +{ + coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]); + coeffs[1] = mix_geom(imodel->coeffs[1], jmodel->coeffs[1]); + coeffs[2] = mix_geom(imodel->coeffs[2], jmodel->coeffs[2]); + coeffs_to_local(); } -void Hooke::coeffs_to_local(){ - k_norm = coeffs[0]; - damp = coeffs[1]; -} +/* ---------------------------------------------------------------------- */ -double Hooke::calculate_forces(){ - contact.area = sqrt(contact.dR); - contact.knfac = k_norm * contact.area; - return contact.knfac * contact.delta; -} - - -//****************** -// Hertz -//****************** -void Hertz::Hertz(ContactModel &c, int mat_flag){ - contact = c; - material_prop_flag = mat_flag; - if (material_prop_flag){ - num_coeffs = 3; - } - else{ - num_coeffs = 2; - } - allocate_coeffs(); -} - -void Hertz::coeffs_to_local(){ - if (material_prop_flag){ - Emod = coeffs[0]; - poiss = coeffs[1]; - k_norm = 4/3*Emod; - } - else{ - k_norm = coeffs[0]; - } -} - -double Hertz::calculate_forces(){ - contact.area = sqrt(contact.dR); - contact.knfac = contact.k_norm * contact.area; - return contact.knfac * contact.delta; -} - -//****************** -// DMT -//****************** -void DMT::DMT(ContactModel &c){ - contact = c; - material_prop_flag = 1; - num_coeffs = 4; - allocate_coeffs(); -} - -double DMT::calculate_forces(){ - contact.area = sqrt(contact.dR); - contact.knfac = contact.k_norm * contact.area; - double Fne = contact.knfac * contact.delta; - F_pulloff = 4 * MathConst::MY_PI * cohesion * contact.Reff; - Fne -= F_pulloff; - return Fne; -} - -void DMT::set_fncrit(){ - contact.Fncrit = fabs(contact.Fne + 2* F_pulloff); -} - -//****************** -// JKR -//****************** -void JKR::JKR(ContactModel &c){ - contact = c; - material_prop_flag = 1; - beyond_contact = 1; - num_coeffs = 4; - allocate_coeffs(); -} - -bool JKR::touch(int touch){ - double Escaled, R2, delta_pulloff, dist_pulloff; - bool touchflag; - - Escaled = k_norm * THREEQUARTERS; - if (touch) { - R2 = contact.Reff * contact.Reff; - a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled)); - delta_pulloff = a * a / contact.Reff - 2 * sqrt(MY_PI * cohesion * a / Escaled); - dist_pulloff = contact.radsum - delta_pulloff; - touchflag = (contact.rsq < dist_pulloff * dist_pulloff); - } else { - touchflag = (rcontact.sq < contact.radsum * contact.radsum); - } - return touchflag; -} - -double JKR::calculate_forces(){ - double Escaled, R2, dR2, t0, t1, t2, t3, t4, t5, t6; - double sqrt1, sqrt2, sqrt3, a2, F_pulloff, Fne; - - Escaled = k_norm * THREEQUARTERS; - - R2 = Reff * Reff; - dR2 = dR * dR; - t0 = cohesion * cohesion * R2 * R2 * Escaled; - t1 = PI27SQ*t0; - t2 = 8 * dR * dR2 * Escaled * Escaled * Escaled; - t3 = 4 * dR2 * Escaled; - - // in case sqrt(0) < 0 due to precision issues - sqrt1 = MAX(0, t0 * (t1 + 2 * t2)); - t4 = cbrt(t1 + t2 + THREEROOT3 * MY_PI * sqrt(sqrt1)); - t5 = t3 / t4 + t4 / Escaled; - sqrt2 = MAX(0, 2 * dR + t5); - t6 = sqrt(sqrt2); - sqrt3 = MAX(0, 4 * dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Escaled * t6)); - a = INVROOT6 * (t6 + sqrt(sqrt3)); - a2 = a * a; - Fne = Escaled * a * a2 / Reff - MY_2PI * a2 * sqrt(4 * cohesion * Escaled / (MY_PI * a)); - F_pulloff = 3 * MY_PI * cohesion * Reff; - - knfac = Escaled * a; - return Fne; -} - -void JKR::set_fncrit(){ - contact.Fncrit = fabs(contact.Fne + 2 * F_pulloff); -} - - - - -void ContactModel::rolling(double *history) +double RollingSDS::calculate_forces() { int rhist0, rhist1, rhist2, frameupdate; - double rolldotn, rollmag, prjmag, magfr, hist_temp[3], temp_dbl, temp_array[3]; + double Frcrit, rolldotn, rollmag, prjmag, magfr, hist_temp[3], temp_dbl, temp_array[3]; - rhist0 = roll_history_index; + rhist0 = history_index; rhist1 = rhist0 + 1; rhist2 = rhist1 + 1; - Frcrit = mu_roll * Fncrit; + Frcrit = mu * contact->normal_model->Fncrit; - if (history_update) { - hist_temp[0] = history[rhist0]; - hist_temp[1] = history[rhist1]; - hist_temp[2] = history[rhist2]; - rolldotn = dot3(hist_temp, nx); + 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_roll > EPSILON*Frcrit; + frameupdate = fabs(rolldotn) * k > EPSILON * Frcrit; if (frameupdate) { // rotate into tangential plane rollmag = len3(hist_temp); // projection temp_dbl = -rolldotn; - scale3(temp_dbl, nx, temp_array); + scale3(temp_dbl, contact->nx, temp_array); sub3(hist_temp, temp_array, hist_temp); // also rescale to preserve magnitude @@ -215,35 +84,36 @@ void ContactModel::rolling(double *history) else temp_dbl = 0; scale3(temp_dbl, hist_temp); } - scale3(dt, vrl, temp_array); + scale3(contact->dt, contact->vrl, temp_array); add3(hist_temp, temp_array, hist_temp); } - scaleadd3(k_roll, hist_temp, gamma_roll, vrl, fr); - negate3(fr); + scaleadd3(k, hist_temp, gamma, contact->vrl, contact->fr); + negate3(contact->fr); // rescale frictional displacements and forces if needed - magfr = len3(fr); + magfr = len3(contact->fr); if (magfr > Frcrit) { rollmag = len3(hist_temp); if (rollmag != 0.0) { - temp_dbl = -Frcrit / (magfr * k_roll); - scale3(temp_dbl, fr, temp_array); + temp_dbl = -Frcrit / (magfr * k); + scale3(temp_dbl, contact->fr, temp_array); add3(hist_temp, temp_array, hist_temp); - temp_dbl = -gamma_roll/k_roll; - scale3(temp_dbl, vrl, temp_array); + temp_dbl = -gamma/k; + scale3(temp_dbl, contact->vrl, temp_array); add3(hist_temp, temp_array, hist_temp); temp_dbl = Frcrit / magfr; - scale3(temp_dbl, fr); + scale3(temp_dbl, contact->fr); } else { - zero3(fr); + zero3(contact->fr); } } - history[rhist0] = hist_temp[0]; - history[rhist1] = hist_temp[1]; - history[rhist2] = hist_temp[2]; + contact->history[rhist0] = hist_temp[0]; + contact->history[rhist1] = hist_temp[1]; + contact->history[rhist2] = hist_temp[2]; + return 0; } diff --git a/src/GRANULAR/contact_rolling_models.h b/src/GRANULAR/contact_rolling_models.h index 1a2f1ddb7d..1cd23e1d5d 100644 --- a/src/GRANULAR/contact_rolling_models.h +++ b/src/GRANULAR/contact_rolling_models.h @@ -11,73 +11,36 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef NORMAL_CONTACT_MODELS_H_ -#define NORMAL_CONTACT_MODELS_H_ +#ifndef CONTACT_ROLLING_MODELS_H_ +#define CONTACT_ROLLING_MODELS_H_ -#include "contact.h"; -#include "sub_model.h" +#include "contact_sub_models.h" +namespace LAMMPS_NS { namespace Contact { -class NormalModel:SubModel{ -public: - NormalModel(){}; - virtual ~NormalModel(){}; - virtual bool check_contact(); - virtual void prep_contact(); - virtual void set_fncrit(); +class RollingModel : public SubModel { + public: + RollingModel() {}; + ~RollingModel() {}; + virtual void coeffs_to_local() {}; + virtual void mix_coeffs(RollingModel*, RollingModel*) {}; virtual double calculate_forces() = 0; - virtual void coeffs_to_local(); - virtual void mix_coeffs(NormalModel*, NormalModel*); //When mixing is needed - -private: - int beyond_contact = 0; - int allow_limit_damping = 1; - }; -class Hooke:NormalModel{ -public: - Hooke(ContactModel &c); - ~Hooke(){}; +/* ---------------------------------------------------------------------- */ + +class RollingSDS: public RollingModel { + public: + RollingSDS(); void coeffs_to_local(); + void mix_coeffs(RollingModel*, RollingModel*); double calculate_forces(); -private: - double k_norm, damp; + private: + double k, mu, gamma; }; -class Hertz:NormalModel{ -public: - Hertz(ContactModel&, int); - ~Hertz(){}; - void coeffs_to_local(); - double calculate_forces(); -private: - double k_norm, damp, Emod, poiss; -}; - -class DMT:NormalModel{ -public: - DMT(ContactModel &c); - ~DMT(){}; - void coeffs_to_local(); - void coeffs_to_local(NormalModel*, NormalModel*); - double calculate_forces(); -private: - double k_norm, damp, Emod, poiss, coh; -}; - -class JKR:NormalModel{ -public: - JKR(ContactModel &c); - ~JKR(){}; - void coeffs_to_local(); - double calculate_forces(); -private: - double k_norm, damp, Emod, poiss, coh; - -}; -} - -#endif /*NORMAL_CONTACT_MODELS_H_ */ +} // namespace Contact +} // namespace LAMMPS_NS +#endif /*CONTACT_ROLLING_MODELS_H_ */ diff --git a/src/GRANULAR/contact_sub_models.cpp b/src/GRANULAR/contact_sub_models.cpp index 99fcd638b4..2809204de7 100644 --- a/src/GRANULAR/contact_sub_models.cpp +++ b/src/GRANULAR/contact_sub_models.cpp @@ -16,23 +16,25 @@ Multiple models can be defined and used to calculate forces and torques based on contact geometry */ -#include "contact_sub_model.h" -#include "pointers.h" + +#include "contact_sub_models.h" #include "utils.h" -#include "error.h" -#include "comm.h" using namespace LAMMPS_NS; +using namespace Contact; -namespace Contact{ - -SubModel::SubModel() +SubModel::SubModel() : + 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; } /* ---------------------------------------------------------------------- */ @@ -40,6 +42,7 @@ SubModel::SubModel() SubModel::~SubModel() { if (allocated) delete [] coeffs; + delete [] transfer_history_factor; } /* ---------------------------------------------------------------------- */ @@ -62,39 +65,6 @@ void SubModel::parse_coeffs(char **arg, int iarg) coeffs_to_local(); } -/* ---------------------------------------------------------------------- */ - -void SubModel::write_restart(FILE *fp) -{ - fwrite(&model_name.length(),sizeof(int),1,fp); - fwrite(model_name.data(),sizeof(char),model_name.length(),fp); - fwrite(&num_coeffs,sizeof(int),1,fp); - fwrite(coeffs,sizeof(int),num_coeffs,fp); -} - -/* ---------------------------------------------------------------------- */ - -void SubModel::read_restart(FILE *fp, int num_char) -{ - if (comm->me == 0) { - utils::sfread(FLERR,&num_coeffs,sizeof(int),1,fp,nullptr,error); - } - MPI_BCast(const_cast(model_name.data()), num_char, MPI_CHAR, world); - allocate_coeffs(); -} - -/* ---------------------------------------------------------------------- */ - -void SubModel::read_restart(FILE *fp) -{ - int num_char; - if (me == 0) { - utils::sfread(FLERR,&num_char,sizeof(int),1,fp,nullptr,error); - } - MPI_BCast(&num_char, 1, MPI_INT, 0, world); - read_restart(fp, num_char); -} - /* ---------------------------------------------------------------------- mixing of Young's modulus (E) ------------------------------------------------------------------------- */ @@ -115,7 +85,7 @@ 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) + double factor2 = 2 * (2 - pois2) * (1 + pois2) / E2; return 1 / (factor1 + factor2); } @@ -127,6 +97,3 @@ double SubModel::mix_geom(double val1, double val2) { return sqrt(val1 * val2); } - - -} \ No newline at end of file diff --git a/src/GRANULAR/contact_sub_models.h b/src/GRANULAR/contact_sub_models.h index f659ef1504..7eca7f4152 100644 --- a/src/GRANULAR/contact_sub_models.h +++ b/src/GRANULAR/contact_sub_models.h @@ -11,44 +11,48 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef GRANULAR_SUB_MODEL_H_ -#define GRANULAR_SUB_MODEL_H_ +#ifndef CONTACT_SUB_MODEL_H_ +#define CONTACT_SUB_MODEL_H_ #include "contact.h" +#include "pointers.h" -using namespace LAMMPS_NS; +namespace LAMMPS_NS { +namespace Contact { -namespace Contact{ - -class SubModel : Pointers{ +class SubModel : protected Pointers { + public: SubModel(); virtual ~SubModel(); -public: + int num_coeffs; double *coeffs; - virtual double calculate_forces() = 0; void read_restart(); - virtual void parse_coeffs(char **, int); - void mix_coeff(SubModel*, SubModel*); - void write_restart(FILE*); - void read_restart(FILE*); - void read_restart(FILE*, int); - virtual void coeffs_to_local(); + void parse_coeffs(char **, int); + virtual void mix_coeffs(SubModel*, SubModel*) {}; + virtual void coeffs_to_local() {}; void allocate_coeffs(); std::string name; -private: - ContactModel &contact; - int allocated; + 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_geom(double, double); }; -} +} // namespace Contact +} // namespace LAMMPS_NS -#endif /* GRANULAR_SUB_MODEL_H_ */ +#endif /* CONTACT_SUB_MODEL_H_ */ diff --git a/src/GRANULAR/contact_tangential_models.cpp b/src/GRANULAR/contact_tangential_models.cpp index 9f092aaf29..45ed993eeb 100644 --- a/src/GRANULAR/contact_tangential_models.cpp +++ b/src/GRANULAR/contact_tangential_models.cpp @@ -11,15 +11,18 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "tangential_contact_models.h" -#include "math_const.h" +#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" -#include - +using namespace LAMMPS_NS; +using namespace Contact; using namespace MathConst; - -namespace Contact{ +using namespace MathExtra; /* ---------------------------------------------------------------------- Linear model with no history @@ -35,14 +38,15 @@ TangentialLinearNoHistory::TangentialLinearNoHistory() void TangentialLinearNoHistory::coeffs_to_local() { + k = 0.0; // No tangential stiffness with no history xt = coeffs[0]; mu = coeffs[1]; - damp = xt * contact.damping_model.damp; + damp = xt * contact->damping_model->damp; } /* ---------------------------------------------------------------------- */ -void TangentialLinearNoHistory::mix_coeffs(NormalModel* imodel, NormalModel* jmodel) +void TangentialLinearNoHistory::mix_coeffs(TangentialModel* imodel, TangentialModel* jmodel) { coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]); coeffs[1] = mix_geom(imodel->coeffs[1], jmodel->coeffs[1]); @@ -56,13 +60,13 @@ void TangentialLinearNoHistory::calculate_forces() double Fscrit, fsmag, Ft; // classic pair gran/hooke (no history) - Fscrit = mu * contact.normal_model.Fncrit - fsmag = damp * contact.vrel; - if (contact.vrel != 0.0) Ft = MIN(Fscrit, fsmag) / contact.vrel; + Fscrit = mu * contact->normal_model->Fncrit; + fsmag = damp * contact->vrel; + if (contact->vrel != 0.0) Ft = MIN(Fscrit, fsmag) / contact->vrel; else Ft = 0.0; Ft = -Ft; - scale3(Ft, contact.vtr, contact.fs); + scale3(Ft, contact->vtr, contact->fs); } /* ---------------------------------------------------------------------- @@ -79,15 +83,15 @@ TangentialLinearHistory::TangentialLinearHistory() void TangentialLinearHistory::coeffs_to_local() { - kt = coeffs[0]; + k = coeffs[0]; xt = coeffs[1]; mu = coeffs[2]; - damp = xt * contact.damping_model.damp; + damp = xt * contact->damping_model->damp; } /* ---------------------------------------------------------------------- */ -void TangentialLinearHistory::mix_coeffs(NormalModel* imodel, NormalModel* jmodel) +void TangentialLinearHistory::mix_coeffs(TangentialModel* imodel, TangentialModel* jmodel) { coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]); coeffs[1] = mix_geom(imodel->coeffs[1], jmodel->coeffs[1]); @@ -102,20 +106,20 @@ void TangentialLinearHistory::calculate_forces() double Fscrit, magfs, rsht, shrmag, prjmag, temp_dbl, temp_array[3]; int frame_update = 0; - Fscrit = contact.normal_model.Fncrit * mu; + Fscrit = contact->normal_model->Fncrit * mu; - double *history = & contact.history[history_index]; + 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) * kt > EPSILON * Fscrit; + 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, history); + scale3(rsht, contact->nx, history); // also rescale to preserve magnitude prjmag = len3(history); if (prjmag > 0) temp_dbl = shrmag / prjmag; @@ -126,28 +130,28 @@ void TangentialLinearHistory::calculate_forces() // update history // tangential force // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 - temp_dbl = kt * contact.dt; - scale3(temp_dbl, contact.vtr, temp_array); + temp_dbl = k * contact->dt; + scale3(temp_dbl, contact->vtr, temp_array); sub3(history, temp_array, history); } // tangential forces = history + tangential velocity damping temp_dbl = -damp; - scale3(temp_dbl, contact.vtr, contact.fs); + scale3(temp_dbl, contact->vtr, contact->fs); // rescale frictional displacements and forces if needed - magfs = len3(contact.fs); + magfs = len3(contact->fs); if (magfs > Fscrit) { shrmag = len3(history); if (shrmag != 0.0) { temp_dbl = Fscrit / magfs; - scale3(temp_dbl, contact.fs, history); - scale3(damp, contact.vtr, temp_array); + scale3(temp_dbl, contact->fs, history); + scale3(damp, contact->vtr, temp_array); add3(history, temp_array, history); temp_dbl = Fscrit / magfs; - scale3(temp_dbl, contact.fs); + scale3(temp_dbl, contact->fs); } else { - zero3(contact.fs); + zero3(contact->fs); } } } @@ -168,20 +172,22 @@ TangentialMindlin::TangentialMindlin() void TangentialMindlin::coeffs_to_local() { - kt = coeffs[0]; + k = coeffs[0]; xt = coeffs[1]; mu = coeffs[2]; - if (ke == -1) - kt = 8.0 * mix_stiffness(contact.normal_model.Emod, contact.normal_model.Emod, - contact.normal_model.poiss, contact.normal_model.poiss); + 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"); + k = 8.0 * mix_stiffnessE(contact->normal_model->Emod, contact->normal_model->Emod, contact->normal_model->poiss, contact->normal_model->poiss); + } - damp = xt * contact.damping_model.damp; + damp = xt * contact->damping_model->damp; } /* ---------------------------------------------------------------------- */ -void TangentialMindlin::mix_coeffs(NormalModel* imodel, NormalModel* jmodel) +void TangentialMindlin::mix_coeffs(TangentialModel* imodel, TangentialModel* jmodel) { if (imodel->coeffs[0] == -1 || imodel->coeffs[0] == -1) coeffs[0] = -1; else coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]); @@ -198,22 +204,22 @@ void TangentialMindlin::calculate_forces() double temp_array[3]; int frame_update = 0; - double *history = & contact.history[history_index]; - Fscrit = contact.normal_model.Fncrit * mu; + double *history = & contact->history[history_index]; + Fscrit = contact->normal_model->Fncrit * mu; - k_scaled = k_tang * contact.area; + k_scaled = k * contact->area; if (mindlin_rescale) { // on unloading, rescale the shear displacements/force - if (contact.area < history[3]) { - temp_dbl = contact.area / history[3]; + if (contact->area < history[3]) { + temp_dbl = contact->area / history[3]; scale3(temp_dbl, 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 (contact->history_update) { + rsht = dot3(history, contact->nx); if (mindlin_force) frame_update = fabs(rsht) > EPSILON * Fscrit; else @@ -222,7 +228,7 @@ void TangentialMindlin::calculate_forces() if (frame_update) { shrmag = len3(history); // projection - scale3(rsht, contact.nx, history); + scale3(rsht, contact->nx, history); // also rescale to preserve magnitude prjmag = len3(history); if (prjmag > 0) temp_dbl = shrmag / prjmag; @@ -234,42 +240,42 @@ void TangentialMindlin::calculate_forces() if (mindlin_force) { // tangential force // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 - temp_dbl = -k_scaled * contact.dt; - scale3(temp_dbl, contact.vtr, temp_array); + temp_dbl = -k_scaled * contact->dt; + scale3(temp_dbl, contact->vtr, temp_array); } else { - scale3(contact.dt, contact.vtr, temp_array); + scale3(contact->dt, contact->vtr, temp_array); } add3(history, temp_array, history); - if (mindlin_rescale) history[3] = contact.area; + if (mindlin_rescale) history[3] = contact->area; } // tangential forces = history + tangential velocity damping temp_dbl = -damp; - scale3(temp_dbl, contact.vtr, contact.fs); + scale3(temp_dbl, contact->vtr, contact->fs); if (! mindlin_force) { scale3(k_scaled, history, temp_array); - add3(contact.fs, temp_array, contact.fs); + add3(contact->fs, temp_array, contact->fs); } // rescale frictional displacements and forces if needed - magfs = len3(contact.fs); + magfs = len3(contact->fs); if (magfs > Fscrit) { shrmag = len3(history); if (shrmag != 0.0) { temp_dbl = Fscrit / magfs; - scale3(temp_dbl, contact.fs, history); - scale3(damp, contact.vtr, temp_array); + scale3(temp_dbl, contact->fs, history); + scale3(damp, contact->vtr, temp_array); add3(history, temp_array, history); if (! mindlin_force) { temp_dbl = -1.0 / k_scaled; scale3(temp_dbl, history); } temp_dbl = Fscrit / magfs; - scale3(temp_dbl, contact.fs); + scale3(temp_dbl, contact->fs); } else { - zero3(contact.fs); + zero3(contact->fs); } } } @@ -278,7 +284,7 @@ void TangentialMindlin::calculate_forces() Mindlin force model ------------------------------------------------------------------------- */ -void TangentialMindlinForce::TangentialMindlinForce() +TangentialMindlinForce::TangentialMindlinForce() { num_coeffs = 3; size_history = 3; @@ -290,22 +296,32 @@ void TangentialMindlinForce::TangentialMindlinForce() Mindlin rescale model ------------------------------------------------------------------------- */ -void TangentialMindlinRescale::TangentialMindlinForce() +TangentialMindlinRescale::TangentialMindlinRescale() { 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 ------------------------------------------------------------------------- */ -void TangentialMindlinRescaleForce::TangentialMindlinForce() +TangentialMindlinRescaleForce::TangentialMindlinRescaleForce() { 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/GRANULAR/contact_tangential_models.h b/src/GRANULAR/contact_tangential_models.h index c0a5a727f0..f5d83a57f7 100644 --- a/src/GRANULAR/contact_tangential_models.h +++ b/src/GRANULAR/contact_tangential_models.h @@ -11,92 +11,84 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef TANGENTIAL_CONTACT_MODELS_H_ -#define TANGENTIAL_CONTACT_MODELS_H_ +#ifndef CONTACT_TANGENTIAL_MODELS_H_ +#define CONTACT_TANGENTIAL_MODELS_H_ -#include "contact.h"; -#include "sub_model.h" +#include "contact_sub_models.h" +namespace LAMMPS_NS { namespace Contact { -class TangentialModel:SubModel -{ -public: +class TangentialModel : public SubModel { + public: TangentialModel() {}; virtual ~TangentialModel() {}; - virtual double calculate_forces() = 0; - virtual void coeffs_to_local(); - virtual void mix_coeffs(TangentialModel*, TangentialModel*); + virtual void coeffs_to_local() {}; + virtual void mix_coeffs(TangentialModel*, TangentialModel*) {}; + virtual void calculate_forces() = 0; int rescale_flag; -private: - int beyond_contact; - int allow_limit_damping; + double k, damp, mu; // Used by Marshall twisting model }; /* ---------------------------------------------------------------------- */ -class TangentialLinearNoHistory:TangentialModel -{ -public: - TangentialLinearNoHistory() - virtual void coeffs_to_local(); - virtual void mix_coeffs(TangentialModel*, TangentialModel*); - double calculate_forces(); -private: - double xt, damp, mu; +class TangentialLinearNoHistory: public TangentialModel { + public: + TangentialLinearNoHistory(); + void coeffs_to_local(); + void mix_coeffs(TangentialModel*, TangentialModel*); + void calculate_forces(); + private: + double xt; }; /* ---------------------------------------------------------------------- */ -class TangentialLinearHistory:TangentialModel -{ -public: - TangentialLinearHistory() - virtual void coeffs_to_local(); - virtual void mix_coeffs(TangentialModel*, TangentialModel*); - double calculate_forces(); -private: - double kt, xt, damp, mu; +class TangentialLinearHistory: public TangentialModel { + public: + TangentialLinearHistory(); + void coeffs_to_local(); + void mix_coeffs(TangentialModel*, TangentialModel*); + void calculate_forces(); + private: + double xt; }; /* ---------------------------------------------------------------------- */ -class TangentialMindlin:TangentialModel -{ -public: +class TangentialMindlin: public TangentialModel { + public: TangentialMindlin(); void coeffs_to_local(); - void coeffs_to_local(TangentialModel*, TangentialModel*); - double calculate_forces(); -private: - double k_norm, damp, Emod, poiss, coh; + void mix_coeffs(TangentialModel*, TangentialModel*); + void calculate_forces(); + protected: int mindlin_rescale, mindlin_force; + double xt; }; /* ---------------------------------------------------------------------- */ -class TangentialMindlinForce:TangentialMindlin -{ -public: +class TangentialMindlinForce: public TangentialMindlin { + public: TangentialMindlinForce(); }; /* ---------------------------------------------------------------------- */ -class TangentialMindlinRescale:TangentialMindlin -{ -public: +class TangentialMindlinRescale: public TangentialMindlin { + public: TangentialMindlinRescale(); }; /* ---------------------------------------------------------------------- */ -class TangentialMindlinRescaleForce:TangentialMindlinRescale -{ -public: - TangentialMindlinForceRescale(); +class TangentialMindlinRescaleForce: public TangentialMindlin { + public: + TangentialMindlinRescaleForce(); }; -} -#endif /*TANGENTIAL_CONTACT_MODELS_H_ */ +} // namespace Contact +} // namespace LAMMPS_NS +#endif /*CONTACT_TANGENTIAL_MODELS_H_ */ diff --git a/src/GRANULAR/contact_twisting_models.cpp b/src/GRANULAR/contact_twisting_models.cpp index 0e9e9115ac..04d89a5b88 100644 --- a/src/GRANULAR/contact_twisting_models.cpp +++ b/src/GRANULAR/contact_twisting_models.cpp @@ -11,206 +11,99 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "normal_contact_models.h" -#include "math_const.h" +#include "contact_normal_models.h" +#include "contact_tangential_models.h" +#include "contact_twisting_models.h" #include "contact.h" +#include "math_const.h" -#include - +using namespace LAMMPS_NS; +using namespace Contact; using namespace MathConst; -namespace Contact{ +/* ---------------------------------------------------------------------- + Marshall twisting model +------------------------------------------------------------------------- */ -#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 behaviors where needed -// ************************ -void NormalModel::set_fncrit(){ - contact->Fncrit = fabs(contact->Fntot); -} - -void NormalModel::mix_coeffs(NormalModel* imodel, NormalModel* jmodel){ - for (int i = 0; i < num_coeffs; i++){ - coeffs[i] = sqrt(imodel->coeffs[i]*jmodel->coeffs[i]); - } -} - -//----------------------------------------- - -//****************** -// Hooke -//****************** -void Hooke::Hooke(ContactModel &c){ - contact = c; - num_coeffs = 2; - allocate_coeffs(); -} - -void Hooke::coeffs_to_local(){ - k_norm = coeffs[0]; - damp = coeffs[1]; -} - -double Hooke::calculate_forces(){ - contact.area = sqrt(contact.dR); - contact.knfac = k_norm * contact.area; - return contact.knfac * contact.delta; -} - - -//****************** -// Hertz -//****************** -void Hertz::Hertz(ContactModel &c, int mat_flag){ - contact = c; - material_prop_flag = mat_flag; - if (material_prop_flag){ - num_coeffs = 3; - } - else{ - num_coeffs = 2; - } - allocate_coeffs(); -} - -void Hertz::coeffs_to_local(){ - if (material_prop_flag){ - Emod = coeffs[0]; - poiss = coeffs[1]; - k_norm = 4/3*Emod; - } - else{ - k_norm = coeffs[0]; - } -} - -double Hertz::calculate_forces(){ - contact.area = sqrt(contact.dR); - contact.knfac = contact.k_norm * contact.area; - return contact.knfac * contact.delta; -} - -//****************** -// DMT -//****************** -void DMT::DMT(ContactModel &c){ - contact = c; - material_prop_flag = 1; - num_coeffs = 4; - allocate_coeffs(); -} - -double DMT::calculate_forces(){ - contact.area = sqrt(contact.dR); - contact.knfac = contact.k_norm * contact.area; - double Fne = contact.knfac * contact.delta; - F_pulloff = 4 * MathConst::MY_PI * cohesion * contact.Reff; - Fne -= F_pulloff; - return Fne; -} - -void DMT::set_fncrit(){ - contact.Fncrit = fabs(contact.Fne + 2* F_pulloff); -} - -//****************** -// JKR -//****************** -void JKR::JKR(ContactModel &c){ - contact = c; - material_prop_flag = 1; - beyond_contact = 1; - num_coeffs = 4; - allocate_coeffs(); -} - -bool JKR::touch(int touch){ - double Escaled, R2, delta_pulloff, dist_pulloff; - bool touchflag; - - Escaled = k_norm * THREEQUARTERS; - if (touch) { - R2 = contact.Reff * contact.Reff; - a = cbrt(9.0 * MY_PI * cohesion * R2 / (4 * Escaled)); - delta_pulloff = a * a / contact.Reff - 2 * sqrt(MY_PI * cohesion * a / Escaled); - dist_pulloff = contact.radsum - delta_pulloff; - touchflag = (contact.rsq < dist_pulloff * dist_pulloff); - } else { - touchflag = (rcontact.sq < contact.radsum * contact.radsum); - } - return touchflag; -} - -double JKR::calculate_forces(){ - double Escaled, R2, dR2, t0, t1, t2, t3, t4, t5, t6; - double sqrt1, sqrt2, sqrt3, a2, F_pulloff, Fne; - - Escaled = k_norm * THREEQUARTERS; - - R2 = Reff * Reff; - dR2 = dR * dR; - t0 = cohesion * cohesion * R2 * R2 * Escaled; - t1 = PI27SQ*t0; - t2 = 8 * dR * dR2 * Escaled * Escaled * Escaled; - t3 = 4 * dR2 * Escaled; - - // in case sqrt(0) < 0 due to precision issues - sqrt1 = MAX(0, t0 * (t1 + 2 * t2)); - t4 = cbrt(t1 + t2 + THREEROOT3 * MY_PI * sqrt(sqrt1)); - t5 = t3 / t4 + t4 / Escaled; - sqrt2 = MAX(0, 2 * dR + t5); - t6 = sqrt(sqrt2); - sqrt3 = MAX(0, 4 * dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Escaled * t6)); - a = INVROOT6 * (t6 + sqrt(sqrt3)); - a2 = a * a; - Fne = Escaled * a * a2 / Reff - MY_2PI * a2 * sqrt(4 * cohesion * Escaled / (MY_PI * a)); - F_pulloff = 3 * MY_PI * cohesion * Reff; - - knfac = Escaled * a; - return Fne; -} - -void JKR::set_fncrit(){ - contact.Fncrit = fabs(contact.Fne + 2 * F_pulloff); -} - - - - -/* ---------------------------------------------------------------------- */ - -void ContactModel::twisting_marshall(double *history) +TwistingMarshall::TwistingMarshall() { - // Overwrite twist coefficients with derived values - k_twist = 0.5 * k_tang * a * a; // eq 32 of Marshall paper - gamma_twist = 0.5 * gamma_tang * a * a; - mu_twist = TWOTHIRDS * a * mu_tang; - - twisting_SDS(history); + num_coeffs = 0; + size_history = 3; } /* ---------------------------------------------------------------------- */ -void ContactModel::twisting_SDS(double *history) +double TwistingMarshall::calculate_forces() { double signtwist, Mtcrit; - if (history_update) { - history[twist_history_index] += magtwist * dt; + // Calculate twist coefficients from tangential model & contact geometry + // eq 32 of Marshall paper + double k = 0.5 * contact->tangential_model->k * contact->area * contact->area; + double damp = 0.5 * contact->tangential_model->damp * contact->area * contact->area; + double mu = TWOTHIRDS * contact->area * contact->tangential_model->mu; + + if (contact->history_update) { + contact->history[history_index] += contact->magtwist * contact->dt; } - magtortwist = -k_twist * history[twist_history_index] - gamma_twist*magtwist; // M_t torque (eq 30) - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit = mu_twist * Fncrit; // critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit) { - history[twist_history_index] = (Mtcrit * signtwist - gamma_twist * magtwist) / k_twist; - magtortwist = -Mtcrit * signtwist; // eq 34 + // 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() +{ + num_coeffs = 3; + size_history = 3; +} + +/* ---------------------------------------------------------------------- */ + +void TwistingSDS::coeffs_to_local() +{ + k = coeffs[0]; + mu = coeffs[1]; + damp = coeffs[2]; +} + +/* ---------------------------------------------------------------------- */ + +void TwistingSDS::mix_coeffs(TwistingModel* imodel, TwistingModel* jmodel) +{ + coeffs[0] = mix_geom(imodel->coeffs[0], jmodel->coeffs[0]); + coeffs[1] = mix_geom(imodel->coeffs[1], jmodel->coeffs[1]); + coeffs[2] = mix_geom(imodel->coeffs[2], jmodel->coeffs[2]); + coeffs_to_local(); +} + +/* ---------------------------------------------------------------------- */ + +double 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/GRANULAR/contact_twisting_models.h b/src/GRANULAR/contact_twisting_models.h index 1a2f1ddb7d..e689c8ed35 100644 --- a/src/GRANULAR/contact_twisting_models.h +++ b/src/GRANULAR/contact_twisting_models.h @@ -11,73 +11,44 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef NORMAL_CONTACT_MODELS_H_ -#define NORMAL_CONTACT_MODELS_H_ +#ifndef CONTACT_TWISTING_MODELS_H_ +#define CONTACT_TWISTING_MODELS_H_ -#include "contact.h"; -#include "sub_model.h" +#include "contact_sub_models.h" +namespace LAMMPS_NS { namespace Contact { -class NormalModel:SubModel{ -public: - NormalModel(){}; - virtual ~NormalModel(){}; - virtual bool check_contact(); - virtual void prep_contact(); - virtual void set_fncrit(); +class TwistingModel : public SubModel { + public: + TwistingModel() {}; + virtual ~TwistingModel() {}; + virtual void coeffs_to_local() {}; + virtual void mix_coeffs(TwistingModel*, TwistingModel*) {}; virtual double calculate_forces() = 0; - virtual void coeffs_to_local(); - virtual void mix_coeffs(NormalModel*, NormalModel*); //When mixing is needed - -private: - int beyond_contact = 0; - int allow_limit_damping = 1; - }; -class Hooke:NormalModel{ -public: - Hooke(ContactModel &c); - ~Hooke(){}; - void coeffs_to_local(); +/* ---------------------------------------------------------------------- */ + +class TwistingMarshall: public TwistingModel { + public: + TwistingMarshall(); double calculate_forces(); -private: - double k_norm, damp; }; -class Hertz:NormalModel{ -public: - Hertz(ContactModel&, int); - ~Hertz(){}; +/* ---------------------------------------------------------------------- */ + +class TwistingSDS: public TwistingModel { + public: + TwistingSDS(); void coeffs_to_local(); + void mix_coeffs(TwistingModel*, TwistingModel*); double calculate_forces(); -private: - double k_norm, damp, Emod, poiss; + private: + double k, mu, damp; }; -class DMT:NormalModel{ -public: - DMT(ContactModel &c); - ~DMT(){}; - void coeffs_to_local(); - void coeffs_to_local(NormalModel*, NormalModel*); - double calculate_forces(); -private: - double k_norm, damp, Emod, poiss, coh; -}; - -class JKR:NormalModel{ -public: - JKR(ContactModel &c); - ~JKR(){}; - void coeffs_to_local(); - double calculate_forces(); -private: - double k_norm, damp, Emod, poiss, coh; - -}; -} - -#endif /*NORMAL_CONTACT_MODELS_H_ */ +} // namespace Contact +} // namespace LAMMPS_NS +#endif /*CONTACT_TWISTING_MODELS_H_ */ diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 5816850ef2..e401f8cb84 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -23,14 +23,19 @@ #include "atom.h" #include "comm.h" #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 "error.h" #include "fix.h" #include "fix_dummy.h" #include "fix_neigh_history.h" #include "force.h" -#include "math_const.h" #include "math_extra.h" -#include "math_special.h" #include "memory.h" #include "modify.h" #include "neigh_list.h" @@ -38,14 +43,12 @@ #include "neighbor.h" #include "update.h" -#include #include #include using namespace LAMMPS_NS; -using namespace MathConst; -using namespace MathSpecial; using namespace Contact; +using namespace MathExtra; /* ---------------------------------------------------------------------- */ @@ -69,8 +72,6 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) maxrad_dynamic = nullptr; maxrad_frozen = nullptr; - history_transfer_factors = nullptr; - dt = update->dt; // set comm size needed by this Pair if used with fix rigid @@ -97,7 +98,6 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) PairGranular::~PairGranular() { delete[] svector; - delete[] history_transfer_factors; if (!fix_history) modify->delete_fix("NEIGH_HISTORY_GRANULAR_DUMMY"); else modify->delete_fix("NEIGH_HISTORY_GRANULAR"); @@ -121,7 +121,7 @@ void PairGranular::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double factor_lj,mi,mj,meff,delx,dely,delz; - double forces[3], torquesi[3], torquesj[3]; + double *forces, *torquesi, *torquesj; int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; @@ -259,17 +259,17 @@ void PairGranular::compute(int eflag, int vflag) torquesj = models[itype][jtype]->torquesj; // apply forces & torques - MathExtra::scale3(factor_lj, forces); - MathExtra::add3(f[i], forces, f[i]); + scale3(factor_lj, forces); + add3(f[i], forces, f[i]); - MathExtra::scale3(factor_lj, torquesi); - MathExtra::add3(torque[i], torquesi, torque[i]); + scale3(factor_lj, torquesi); + add3(torque[i], torquesi, torque[i]); if (heat_flag) heatflux[i] += dq; if (force->newton_pair || j < nlocal) { - MathExtra::sub3(f[j], forces, f[j]); - MathExtra::scale3(factor_lj, torquesj); - MathExtra::add3(torque[j], torquesj, torque[j]); + sub3(f[j], forces, f[j]); + scale3(factor_lj, torquesj); + add3(torque[j], torquesj, torque[j]); if (heat_flag) heatflux[j] -= dq; } @@ -301,8 +301,7 @@ void PairGranular::allocate() memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(cutoff_type,n+1,n+1,"pair:cutoff_type"); - - models = (ContactModel **) memory->srealloc(fix, n+1,n+1, "pair:contact_models"); + memory->create(models,n+1,n+1,"pair:contact_models"); onerad_dynamic = new double[n+1]; onerad_frozen = new double[n+1]; @@ -349,26 +348,26 @@ void PairGranular::coeff(int narg, char **arg) //Parse mandatory normal and tangential specifications int iarg = 2; - vec_models.back().init(arg[iarg], Contact::NORMAL); - ncoeff = vec_models.back().normal_model.num_coeffs + vec_models.back().init_model(std::string(arg[iarg]), NORMAL); + ncoeff = vec_models.back().normal_model->num_coeffs; iarg += 1; if (iarg + ncoeff >= narg) error->all(FLERR,"Illegal pair_coeff command" "Insufficient arguments provided for normal model."); - vec_models.back().normal_model.parse_coeffs(arg, iarg); + vec_models.back().normal_model->parse_coeffs(arg, iarg); iarg += ncoeff; if (strcmp(arg[iarg], "tangential") == 0) { if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_coeff command, must specify " "tangential model after tangential keyword"); - vec_models.back().init(arg[iarg], Contact::TANGENTIAL); - ncoeff = vec_models.back().tangential_model.num_coeffs; + vec_models.back().init_model(std::string(arg[iarg]), TANGENTIAL); + ncoeff = vec_models.back().tangential_model->num_coeffs; iarg += 1; if (iarg + ncoeff >= narg) error->all(FLERR, "Illegal pair_coeff command" "Insufficient arguments provided for tangential model."); - vec_models.back().tangential_model.parse_coeffs(arg, iarg); + vec_models.back().tangential_model->parse_coeffs(arg, iarg); iarg += ncoeff; } else{ error->all(FLERR, "Illegal pair_coeff command, 'tangential' keyword expected"); @@ -379,49 +378,49 @@ void PairGranular::coeff(int narg, char **arg) if (strcmp(arg[iarg], "damping") == 0) { if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - vec_models.back().init(arg[iarg], Contact::DAMPING); - ncoeff = vec_models.back().damping_model.num_coeffs; + vec_models.back().init_model(std::string(arg[iarg]), DAMPING); + ncoeff = vec_models.back().damping_model->num_coeffs; iarg += 1; if (iarg + ncoeff >= narg) error->all(FLERR, "Illegal pair_coeff command" "Insufficient arguments provided for damping model."); - vec_models.back().damping_model.parse_coeffs(arg, iarg); + vec_models.back().damping_model->parse_coeffs(arg, iarg); iarg += ncoeff; } else if (strcmp(arg[iarg], "rolling") == 0) { if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - vec_models.back().init(arg[iarg], Contact::ROLLING); - ncoeff = vec_models.back().rolling_model.num_coeffs; + vec_models.back().init_model(std::string(arg[iarg]), ROLLING); + ncoeff = vec_models.back().rolling_model->num_coeffs; iarg += 1; if (iarg + ncoeff >= narg) error->all(FLERR, "Illegal pair_coeff command" "Insufficient arguments provided for rolling model."); - vec_models.back().rolling_model.parse_coeffs(arg, iarg); + vec_models.back().rolling_model->parse_coeffs(arg, iarg); iarg += ncoeff; } else if (strcmp(arg[iarg], "twisting") == 0) { if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - vec_models.back().init(arg[iarg], Contact::TWISTING); - ncoeff = vec_models.back().twisting_model.num_coeffs; + vec_models.back().init_model(std::string(arg[iarg]), TWISTING); + ncoeff = vec_models.back().twisting_model->num_coeffs; iarg += 1; if (iarg + ncoeff >= narg) error->all(FLERR, "Illegal pair_coeff command" "Insufficient arguments provided for twisting model."); - vec_models.back().twisting_model.parse_coeffs(arg, iarg); + vec_models.back().twisting_model->parse_coeffs(arg, iarg); iarg += ncoeff; } else if (strcmp(arg[iarg], "heat") == 0) { if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - vec_models.back().init(arg[iarg], Contact::HEAT); - ncoeff = vec_models.back().heat_model.num_coeffs; + vec_models.back().init_model(std::string(arg[iarg]), HEAT); + ncoeff = vec_models.back().heat_model->num_coeffs; iarg += 1; if (iarg + ncoeff >= narg) error->all(FLERR, "Illegal pair_coeff command" "Insufficient arguments provided for heat model."); - vec_models.back().heat_model.parse_coeffs(arg, iarg); + vec_models.back().heat_model->parse_coeffs(arg, iarg); iarg += ncoeff; heat_flag = 1; @@ -440,8 +439,8 @@ void PairGranular::coeff(int narg, char **arg) // Define default damping model if unspecified, takes no args if (!vec_models.back().damping_model) { - vec_models.back().init("viscoelastic", Contact::DAMPING); - vec_models.back().damping_model.parse_coeffs(arg, 0); + vec_models.back().init_model("viscoelastic", DAMPING); + vec_models.back().damping_model->parse_coeffs(arg, 0); } if (vec_models.back().limit_damping && !vec_models.back().normal_model->allow_limit_damping) @@ -475,38 +474,41 @@ void PairGranular::init_style() if (comm->ghost_velocity == 0) error->all(FLERR,"Pair granular requires ghost atoms store velocity"); - // allocate history + // allocate history and initialize models int size_normal_history = 0; + int size_damping_history = 0; int size_tangential_history = 0; int size_rolling_history = 0; int size_twisting_history = 0; for (int i = 1; i <= atom->ntypes; i++) { for (int j = i; j <= atom->ntypes; j++) { - if (models[i][j]->normal_model.history_flag || - models[i][j]->tangential_model.history_flag || - models[i][j]->rolling_model.history_flag || - models[i][j]->twisting_model.history_flag) use_history = 1; - - if (models[i][j]->nondefault_history_transfer) // ??? + if (models[i][j]->normal_model->size_history != 0 || + models[i][j]->damping_model->size_history != 0 || + models[i][j]->tangential_model->size_history != 0 || + models[i][j]->rolling_model->size_history != 0 || + models[i][j]->twisting_model->size_history != 0) use_history = 1; if (models[i][j]->normal_model->size_history > size_normal_history) - size_normal_history = models[i][j]->normal_model.size_history; + size_normal_history = models[i][j]->damping_model->size_history; + if (models[i][j]->damping_model->size_history > size_damping_history) + size_damping_history = models[i][j]->normal_model->size_history; if (models[i][j]->tangential_model->size_history > size_tangential_history) - size_tangential_history = models[i][j]->tangential_model.size_history; + size_tangential_history = models[i][j]->tangential_model->size_history; if (models[i][j]->rolling_model->size_history > size_rolling_history) - size_rolling_history = models[i][j]->rolling_model.size_history; + size_rolling_history = models[i][j]->rolling_model->size_history; if (models[i][j]->twisting_model->size_history > size_twisting_history) - size_twisting_history = models[i][j]->twisting_model.size_history; + size_twisting_history = models[i][j]->twisting_model->size_history; } } - size_history = size_normal_history + size_tangential_history + - size_rolling_history + size_twisting_history; + size_history = size_normal_history + size_damping_history + + size_tangential_history + size_rolling_history + size_twisting_history; - tangential_history_index = size_normal_history; - roll_history_index = size_normal_history + size_tangential_history; - twist_history_index = size_normal_history + size_tangential_history + size_rolling_history; + damping_history_index = size_normal_history; + tangential_history_index = size_normal_history + damping_history_index; + roll_history_index = size_normal_history + damping_history_index + size_tangential_history; + twist_history_index = size_normal_history + damping_history_index + size_tangential_history + size_rolling_history; if (use_history) neighbor->add_request(this, NeighConst::REQ_SIZE|NeighConst::REQ_HISTORY); else neighbor->add_request(this, NeighConst::REQ_SIZE); @@ -601,11 +603,11 @@ double PairGranular::init_one(int i, int j) double cutoff = 0.0; if (setflag[i][j] == 0) { - if ((models[i][i]->normal_model.name != models[j][j]->normal_model.name) || - (models[i][i]->damping_model.name != models[j][j]->damping_model.name) || - (models[i][i]->tangential_model.name != models[j][j]->tangential_model.name) || - (models[i][i]->rolling_model.name != models[j][j]->rolling_model.name) || - (models[i][i]->twisting_model.name != models[j][j]->twisting_model.name)) { + if ((models[i][i]->normal_model->name != models[j][j]->normal_model->name) || + (models[i][i]->damping_model->name != models[j][j]->damping_model->name) || + (models[i][i]->tangential_model->name != models[j][j]->tangential_model->name) || + (models[i][i]->rolling_model->name != models[j][j]->rolling_model->name) || + (models[i][i]->twisting_model->name != models[j][j]->twisting_model->name)) { error->all(FLERR,"Granular pair style functional forms are different, " "cannot mix coefficients for types {} and {}. \n" "This combination must be set explicitly via a " @@ -614,12 +616,12 @@ double PairGranular::init_one(int i, int j) vec_models.push_back(ContactModel()); models[i][j] = models[j][i] = & vec_models.back(); - vec_models.back().init(models[i][i]->normal_model.name, Contact::NORMAL); - vec_models.back().init(models[i][i]->tangential_model.name, Contact::TANGENTIAL); - vec_models.back().init(models[i][i]->damping_model.name, Contact::DAMPING); - vec_models.back().init(models[i][i]->rolling_model.name, Contact::ROLLING); - vec_models.back().init(models[i][i]->twisting_model.name, Contact::TWISTING); - vec_models.back().init(models[i][i]->heat_model.name, Contact::HEAT); + vec_models.back().init_model(models[i][i]->normal_model->name, NORMAL); + vec_models.back().init_model(models[i][i]->tangential_model->name, TANGENTIAL); + vec_models.back().init_model(models[i][i]->damping_model->name, DAMPING); + vec_models.back().init_model(models[i][i]->rolling_model->name, ROLLING); + vec_models.back().init_model(models[i][i]->twisting_model->name, TWISTING); + vec_models.back().init_model(models[i][i]->heat_model->name, HEAT); vec_models.back().mix_coeffs(models[i][i], models[j][j]); } @@ -675,10 +677,16 @@ double PairGranular::init_one(int i, int j) // Copy global options models[i][j]->dt = models[j][i]->dt = dt; - models[i][j]->normal_model.history_index = models[j][i]->normal_model.history_index = normal_history_index; - models[i][j]->tangential_model.history_index = models[j][i]->tangential_model.history_index = tangential_history_index; - models[i][j]->rolling_model.history_index = models[j][i]->rolling_model.history_index = rolling_history_index; - models[i][j]->twisting_model.history_index = models[j][i]->twisting_model.history_index = twisting_history_index; + models[i][j]->normal_model->history_index = models[j][i]->normal_model->history_index = normal_history_index; + models[i][j]->tangential_model->history_index = models[j][i]->tangential_model->history_index = tangential_history_index; + models[i][j]->rolling_model->history_index = models[j][i]->rolling_model->history_index = roll_history_index; + models[i][j]->twisting_model->history_index = models[j][i]->twisting_model->history_index = twist_history_index; + + models[i][j]->size_history = models[j][i]->size_history = size_history; + models[i][j]->init(); // Calculates cumulative properties of sub models + models[j][i]->init(); + + if (models[i][j]->nondefault_history_transfer) nondefault_history_transfer = 1; return cutoff; } @@ -718,7 +726,7 @@ void PairGranular::read_restart(FILE *fp) if (setflag[i][j]) { vec_models.push_back(ContactModel()); models[i][j] = & vec_models.back(); - models[i][j]->read_restart(); + models[i][j]->read_restart(fp); } } } @@ -729,6 +737,12 @@ void PairGranular::read_restart(FILE *fp) void PairGranular::reset_dt() { dt = update->dt; + + for (int i = 1; i <= atom->ntypes; i++) { + for (int j = i; j <= atom->ntypes; j++) { + models[i][j]->dt = dt; + } + } } /* ---------------------------------------------------------------------- */ @@ -819,7 +833,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, history = &allhistory[size_history*neighprev]; } - double forces[3], torquesi[3], torquesj[3]; + double *forces, *torquesi, *torquesj; models[itype][jtype]->calculate_forces(); forces = models[itype][jtype]->forces; torquesi = models[itype][jtype]->torquesi; @@ -893,11 +907,16 @@ double PairGranular::memory_usage() only needed if any history entries i-j are not just negative of j-i entries ------------------------------------------------------------------------- */ -void PairGranular::transfer_history(double* source, double* target, int source_type, int target_type) +void PairGranular::transfer_history(double* source, double* target, int itype, int jtype) { - models[itype][jtype]->transfer_history(); - for (int i = 0; i < size_history; i++){ - target[i] = history_transfer_factors[i]*source[i]; + if (models[itype][jtype]->nondefault_history_transfer) { + for (int i = 0; i < size_history; i++) { + target[i] = models[itype][jtype]->transfer_history_factor[i] * source [i]; + } + } else { + for (int i = 0; i < size_history; i++) { + target[i] = -source[i]; + } } } diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index e018f9fe7c..1032cf6b7c 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -22,7 +22,7 @@ PairStyle(granular,PairGranular); #include "contact.h" #include "pair.h" -#include "vector.h" +#include namespace LAMMPS_NS { @@ -65,11 +65,10 @@ class PairGranular : public Pair { int nmax; // allocated size of mass_rigid void allocate(); - void transfer_history(double *, double *) override; + void transfer_history(double *, double *, int, int) override; private: int size_history; - int *history_transfer_factors; int heat_flag; // contact models @@ -81,6 +80,7 @@ class PairGranular : public Pair { // indices of history entries int normal_history_index; + int damping_history_index; int tangential_history_index; int roll_history_index; int twist_history_index; diff --git a/src/contact.h b/src/contact.h index 167f88e9de..96758771b4 100644 --- a/src/contact.h +++ b/src/contact.h @@ -14,78 +14,89 @@ #ifndef LMP_CONTACT_H #define LMP_CONTACT_H +#include "pointers.h" // IWYU pragma: export + +namespace LAMMPS_NS { namespace Contact { - enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; - enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; - enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, - TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, - TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; - enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; - enum {ROLL_NONE, ROLL_SDS}; +enum ModelType { + NORMAL = 0, + DAMPING = 1, + TANGENTIAL = 2, + ROLLING = 3, + TWISTING = 4, + HEAT = 5 +}; - #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 +#define EPSILON 1e-10 - #define EPSILON 1e-10 +// forward declaration +class NormalModel; +class DampingModel; +class TangentialModel; +class RollingModel; +class TwistingModel; +class HeatModel; +class SubModel; - class ContactModel { - public: - ContactModel(); - void reset_contact(); - bool check_contact(); - void prep_contact(); - void calculate_forces(double *, double *, double *, double *); - double calculate_heat(); - double pulloff_distance(double, double); +class ContactModel : protected Pointers { + public: + ContactModel(); + ~ContactModel(); + void init(); + bool check_contact(); + void reset_contact(); + void prep_contact(); + void calculate_forces(); + double calculate_heat(); + double pulloff_distance(double, double); - int normal_model, damping_model, tangential_model; - int roll_model, twist_model; - int limit_damping; - double cutoff_type; - double Emod, poisson; // variables used in defining mixed interactions - double k_norm, gamma_norm, cohesion; // normal_coeffs - double k_tang, gamma_tang, mu_tang; // tangential_coeffs - wutang? - double k_roll, gamma_roll, mu_roll; // roll_coeffs - double k_twist, gamma_twist, mu_twist; // twist_coeffs - double conductivity; + void init_model(std::string, ModelType); - double radi, radj, meff, dt, Ti, Tj; - double *xi, *xj, *vi, *vj, *omegai, *omegaj; - int history_update, roll_history_index, twist_history_index; + void mix_coeffs(ContactModel*, ContactModel*); - double fs[3], fr[3], ft[3], magtortwist; + void write_restart(FILE *); + void read_restart(FILE *); - private: - double a, knfac, Fntot, Fncrit, Fscrit, Frcrit, damp_normal_prefactor; - 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; + // Sub models + NormalModel *normal_model; + DampingModel *damping_model; + TangentialModel *tangential_model; + RollingModel *rolling_model; + TwistingModel *twisting_model; + HeatModel *heat_model; + SubModel *sub_models[6]; // Need to resize if we add more model flavors - int prep_flag, check_flag; - int mindlin_rescale, mindlin_force; + // Extra options + int beyond_contact, limit_damping, history_update; + double cutoff_type; - bool touch_JKR(int); - double normal_JKR(); - double normal_DMT(); - double normal_Hertz(); - double normal_Hooke(); - double normal_damping(); - void tangential_no_history(); - void tangential_history(double *); - void tangential_mindlin(double *); - void rolling(double *); - void twisting_marshall(double *); - void twisting_SDS(double *); + // History variables + int size_history, nondefault_history_transfer; + double *transfer_history_factor; + double *history; - }; + // Contact properties/output + double *forces, *torquesi, *torquesj; + + double radi, radj, meff, dt, Ti, Tj, area; + double Fntot, magtortwist; + + 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 prep_flag, check_flag; + int nmodels; +}; } // namespace Contact +} // namespace LAMMPS_NS #endif diff --git a/src/fix_neigh_history.cpp b/src/fix_neigh_history.cpp index 37d1b37946..975d93ff83 100644 --- a/src/fix_neigh_history.cpp +++ b/src/fix_neigh_history.cpp @@ -337,6 +337,7 @@ void FixNeighHistory::pre_exchange_newton() int *ilist,*jlist,*numneigh,**firstneigh; int *allflags; double *allvalues,*onevalues,*jvalues; + int *type = atom->type; // NOTE: all operations until very end are with // nlocal_neigh <= current nlocal and nall_neigh @@ -430,7 +431,7 @@ void FixNeighHistory::pre_exchange_newton() partner[j][m] = tag[i]; jvalues = &valuepartner[j][dnum*m]; if (pair->nondefault_history_transfer) - pair->transfer_history(onevalues,jvalues); + pair->transfer_history(onevalues,jvalues,type[i],type[j]); else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; } } @@ -469,6 +470,7 @@ void FixNeighHistory::pre_exchange_no_newton() int *ilist,*jlist,*numneigh,**firstneigh; int *allflags; double *allvalues,*onevalues,*jvalues; + int *type = atom->type; // NOTE: all operations until very end are with nlocal_neigh <= current nlocal // because previous neigh list was built with nlocal_neigh @@ -544,7 +546,7 @@ void FixNeighHistory::pre_exchange_no_newton() partner[j][m] = tag[i]; jvalues = &valuepartner[j][dnum*m]; if (pair->nondefault_history_transfer) - pair->transfer_history(onevalues, jvalues); + pair->transfer_history(onevalues, jvalues,type[i],type[j]); else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; } } diff --git a/src/pair.h b/src/pair.h index 4763a225d2..9e72f03d23 100644 --- a/src/pair.h +++ b/src/pair.h @@ -211,7 +211,7 @@ class Pair : protected Pointers { virtual void min_xf_pointers(int, double **, double **) {} virtual void min_xf_get(int) {} virtual void min_x_set(int) {} - virtual void transfer_history(double *, double *) {} + virtual void transfer_history(double *, double *, int, int) {} virtual double atom2cut(int) { return 0.0; } virtual double radii2cut(double, double) { return 0.0; }