ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Removing old contact files, fixing capitalization in dump_custom.cpp

Browse Source
This commit is contained in:
jtclemm
2022-10-28 16:38:32 -06:00
parent 1eb048ef93
commit f66179f336
17 changed files with 3 additions and 2481 deletions
Download Patch File Download Diff File Expand all files Collapse all files

479
src/contact.cpp
View File

@ -1,479 +0,0 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
This class contains a series of tools for DEM contacts
Multiple models can be defined and used to calculate forces
and torques based on contact geometry
Contributing authors:
Dan Bolintineanu (SNL), Joel Clemmer (SNL)
----------------------------------------------------------------------- */
#include "contact.h"
#include "contact_sub_models.h"
#include "contact_normal_models.h"
#include "contact_tangential_models.h"
#include "contact_damping_models.h"
#include "contact_rolling_models.h"
#include "contact_twisting_models.h"
#include "contact_heat_models.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_extra.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathExtra;
ContactModel::ContactModel(LAMMPS *lmp) : Pointers(lmp)
{
limit_damping = 0;
beyond_contact = 0;
nondefault_history_transfer = 0;
classic_model = 0;
contact_type = PAIR;
normal_model = nullptr;
damping_model = nullptr;
tangential_model = nullptr;
rolling_model = nullptr;
twisting_model = nullptr;
heat_model = nullptr;
for (int i = 0; i < NSUBMODELS; i++) sub_models[i] = nullptr;
transfer_history_factor = nullptr;
}
/* ---------------------------------------------------------------------- */
ContactModel::~ContactModel()
{
delete [] transfer_history_factor;
delete normal_model;
delete damping_model;
delete tangential_model;
delete rolling_model;
delete twisting_model;
delete heat_model;
}
/* ---------------------------------------------------------------------- */
void ContactModel::init_model(std::string model_name, ModelType model_type)
{
if (model_type == NORMAL) {
delete normal_model;
if (model_name == "none") normal_model = new NormalNone(lmp);
else if (model_name == "hooke") normal_model = new NormalHooke(lmp);
else if (model_name == "hertz") normal_model = new NormalHertz(lmp);
else if (model_name == "hertz/material") normal_model = new NormalHertzMaterial(lmp);
else if (model_name == "dmt") normal_model = new NormalDMT(lmp);
else if (model_name == "jkr") normal_model = new NormalJKR(lmp);
else error->all(FLERR, "Normal model name {} not recognized", model_name);
sub_models[model_type] = normal_model;
} else if (model_type == TANGENTIAL) {
delete tangential_model;
if (model_name == "none") tangential_model = new TangentialNone(lmp);
else if (model_name == "linear_nohistory") tangential_model = new TangentialLinearNoHistory(lmp);
else if (model_name == "linear_history") tangential_model = new TangentialLinearHistory(lmp);
else if (model_name == "linear_history_classic") tangential_model = new TangentialLinearHistoryClassic(lmp);
else if (model_name == "mindlin") tangential_model = new TangentialMindlin(lmp);
else if (model_name == "mindlin/force") tangential_model = new TangentialMindlinForce(lmp);
else if (model_name == "mindlin_rescale") tangential_model = new TangentialMindlinRescale(lmp);
else if (model_name == "mindlin_rescale/force") tangential_model = new TangentialMindlinRescaleForce(lmp);
else error->all(FLERR, "Tangential model name {} not recognized", model_name);
sub_models[model_type] = tangential_model;
} else if (model_type == DAMPING) {
delete damping_model;
if (model_name == "none") damping_model = new DampingNone(lmp);
else if (model_name == "velocity") damping_model = new DampingVelocity(lmp);
else if (model_name == "mass_velocity") damping_model = new DampingMassVelocity(lmp);
else if (model_name == "viscoelastic") damping_model = new DampingViscoelastic(lmp);
else if (model_name == "tsuji") damping_model = new DampingTsuji(lmp);
else error->all(FLERR, "Damping model name {} not recognized", model_name);
sub_models[model_type] = damping_model;
} else if (model_type == ROLLING) {
delete rolling_model;
rolling_defined = 1;
if (model_name == "none") {
rolling_model = new RollingNone(lmp);
rolling_defined = 0;
} else if (model_name == "sds") rolling_model = new RollingSDS(lmp);
else error->all(FLERR, "Rolling model name {} not recognized", model_name);
sub_models[model_type] = rolling_model;
} else if (model_type == TWISTING) {
delete twisting_model;
twisting_defined = 1;
if (model_name == "none") {
twisting_model = new TwistingNone(lmp);
twisting_defined = 0;
} else if (model_name == "sds") twisting_model = new TwistingSDS(lmp);
else if (model_name == "marshall") twisting_model = new TwistingMarshall(lmp);
else error->all(FLERR, "Twisting model name {} not recognized", model_name);
sub_models[model_type] = twisting_model;
} else if (model_type == HEAT) {
delete heat_model;
heat_defined = 1;
if (model_name == "none") {
heat_model = new HeatNone(lmp);
heat_defined = 0;
} else if (model_name == "area") heat_model = new HeatArea(lmp);
else error->all(FLERR, "Heat model name not {} recognized", model_name);
sub_models[model_type] = heat_model;
} else {
error->all(FLERR, "Illegal model type {}", model_type);
}
sub_models[model_type]->name.assign(model_name);
sub_models[model_type]->contact = this;
sub_models[model_type]->allocate_coeffs();
}
/* ---------------------------------------------------------------------- */
int ContactModel::init_classic_model(char **arg, int iarg, int narg)
{
double kn, kt, gamman, gammat, xmu;
classic_model = 1;
if (iarg + 6 >= narg)
error->all(FLERR,"Insufficient arguments provided for classic gran model command");
kn = utils::numeric(FLERR,arg[iarg + 1],false,lmp);
if (strcmp(arg[iarg + 2],"NULL") == 0) kt = kn * 2.0 / 7.0;
else kt = utils::numeric(FLERR,arg[iarg + 2],false,lmp);
gamman = utils::numeric(FLERR,arg[iarg + 3],false,lmp);
if (strcmp(arg[iarg + 4],"NULL") == 0) gammat = 0.5 * gamman;
else gammat = utils::numeric(FLERR,arg[iarg + 4],false,lmp);
xmu = utils::numeric(FLERR,arg[iarg + 5],false,lmp);
int dampflag = utils::inumeric(FLERR,arg[iarg + 6],false,lmp);
if (dampflag == 0) gammat = 0.0;
if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 ||
xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1)
error->all(FLERR,"Illegal classic gran model command");
if (strcmp(arg[iarg],"hooke") == 0) {
init_model("hooke", NORMAL);
init_model("linear_nohistory", TANGENTIAL);
init_model("mass_velocity", DAMPING);
} else if (strcmp(arg[iarg],"hooke/history") == 0) {
init_model("hooke", NORMAL);
init_model("linear_history_classic", TANGENTIAL);
init_model("mass_velocity", DAMPING);
} else if (strcmp(arg[iarg],"hertz/history") == 0) {
// convert Kn and Kt from pressure units to force/distance^2 if Hertzian
kn /= force->nktv2p;
kt /= force->nktv2p;
init_model("hertz", NORMAL);
init_model("mindlin", TANGENTIAL);
init_model("viscoelastic", DAMPING);
} else error->all(FLERR,"Invalid classic gran model");
// ensure additional models are undefined
init_model("none", ROLLING);
init_model("none", TWISTING);
init_model("none", HEAT);
// manually parse coeffs
normal_model->coeffs[0] = kn;
normal_model->coeffs[1] = gamman;
tangential_model->coeffs[0] = kt;
tangential_model->coeffs[1] = gammat / gamman;
tangential_model->coeffs[2] = xmu;
normal_model->coeffs_to_local();
tangential_model->coeffs_to_local();
damping_model->coeffs_to_local();
iarg += 7;
return iarg;
}
/* ---------------------------------------------------------------------- */
void ContactModel::init()
{
int i, j;
for (i = 0; i < NSUBMODELS; i++)
if (!sub_models[i]) init_model("none", (ModelType) i);
// Must have valid normal, damping, and tangential models
if (normal_model->name == "none") error->all(FLERR, "Must specify normal contact model");
if (damping_model->name == "none") error->all(FLERR, "Must specify damping contact model");
if (tangential_model->name == "none") error->all(FLERR, "Must specify tangential contact model");
int size_cumulative;
size_history = 0;
for (i = 0; i < NSUBMODELS; i++) {
if (sub_models[i]->nondefault_history_transfer)
nondefault_history_transfer = 1;
if (sub_models[i]->beyond_contact)
beyond_contact = 1;
size_history += sub_models[i]->size_history;
}
if (nondefault_history_transfer) {
transfer_history_factor = new double[size_history];
for (i = 0; i < size_history; i++) {
// Find which model owns this history value
size_cumulative = 0;
for (j = 0; j < NSUBMODELS; j++) {
if (size_cumulative + sub_models[j]->size_history > i) break;
size_cumulative += sub_models[j]->size_history;
}
// Check if model has nondefault transfers, if so copy its array
transfer_history_factor[i] = -1;
if (j != NSUBMODELS) {
if (sub_models[j]->nondefault_history_transfer) {
transfer_history_factor[i] = sub_models[j]->transfer_history_factor[i - size_cumulative];
}
}
}
}
for (i = 0; i < NSUBMODELS; i++) sub_models[i]->init();
}
/* ---------------------------------------------------------------------- */
int ContactModel::mix_coeffs(ContactModel *c1, ContactModel *c2)
{
int i;
for (i = 0; i < NSUBMODELS; i++) {
if (c1->sub_models[i]->name != c2->sub_models[i]->name) return i;
init_model(c1->sub_models[i]->name, (ModelType) i);
sub_models[i]->mix_coeffs(c1->sub_models[i]->coeffs, c2->sub_models[i]->coeffs);
}
limit_damping = MAX(c1->limit_damping, c2->limit_damping);
return -1;
}
/* ---------------------------------------------------------------------- */
void ContactModel::write_restart(FILE *fp)
{
int num_char, num_coeffs;
for (int i = 0; i < NSUBMODELS; i++) {
num_char = sub_models[i]->name.length();
num_coeffs = sub_models[i]->num_coeffs;
fwrite(&num_char, sizeof(int), 1, fp);
fwrite(sub_models[i]->name.data(), sizeof(char), num_char, fp);
fwrite(&num_coeffs, sizeof(int), 1, fp);
fwrite(sub_models[i]->coeffs, sizeof(double), num_coeffs, fp);
}
}
/* ---------------------------------------------------------------------- */
void ContactModel::read_restart(FILE *fp)
{
int num_char, num_coeff;
for (int i = 0; i < NSUBMODELS; i++) {
if (comm->me == 0)
utils::sfread(FLERR, &num_char, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&num_char, 1, MPI_INT, 0, world);
std::string model_name (num_char, ' ');
if (comm->me == 0)
utils::sfread(FLERR, const_cast<char*>(model_name.data()), sizeof(char),num_char, fp, nullptr, error);
MPI_Bcast(const_cast<char*>(model_name.data()), num_char, MPI_CHAR, 0, world);
init_model(model_name, (ModelType) i);
if (comm->me == 0) {
utils::sfread(FLERR, &num_coeff, sizeof(int), 1, fp, nullptr, error);
if (num_coeff != sub_models[i]->num_coeffs)
error->one(FLERR, "Invalid contact model written to restart file");
}
MPI_Bcast(&num_coeff, 1, MPI_INT, 0, world);
if (comm->me == 0) {
utils::sfread(FLERR, sub_models[i]->coeffs, sizeof(double), num_coeff, fp, nullptr, error);
}
MPI_Bcast(sub_models[i]->coeffs, num_coeff, MPI_DOUBLE, 0, world);
sub_models[i]->coeffs_to_local();
}
}
/* ---------------------------------------------------------------------- */
bool ContactModel::check_contact()
{
if (contact_type == WALL) {
// Used by fix_wall_gran.cpp
// radj = radius of wall
// dx already provided
rsq = lensq3(dx);
radsum = radi;
if (radj == 0) Reff = radi;
else Reff = radi * radj / (radi + radj);
} else if (contact_type == WALLREGION) {
// Used by fix_wall_gran_region.cpp
// radj = radius of wall
// dx and r already provided
rsq = r * r;
radsum = radi;
if (radj == 0) Reff = radi;
else Reff = radi * radj / (radi + radj);
} else {
sub3(xi, xj, dx);
rsq = lensq3(dx);
radsum = radi + radj;
Reff = radi * radj / radsum;
}
touch = normal_model->touch();
return touch;
}
/* ---------------------------------------------------------------------- */
void ContactModel::prep_contact()
{
double temp[3];
// Standard geometric quantities
if (contact_type != WALLREGION) r = sqrt(rsq);
rinv = 1.0 / r;
delta = radsum - r;
dR = delta * Reff;
scale3(rinv, dx, nx);
// relative translational velocity
sub3(vi, vj, vr);
// normal component
vnnr = dot3(vr, nx); //v_R . n
scale3(vnnr, nx, vn);
// tangential component
sub3(vr, vn, vt);
// relative rotational velocity
scaleadd3(radi, omegai, radj, omegaj, wr);
// relative tangential velocities
cross3(wr, nx, temp);
sub3(vt, temp, vtr);
vrel = len3(vtr);
if (rolling_defined || twisting_defined)
sub3(omegai, omegaj, relrot);
if (rolling_defined) {
// rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
// this is different from the Marshall papers, which use the Bagi/Kuhn formulation
// for rolling velocity (see Wang et al for why the latter is wrong)
vrl[0] = Reff * (relrot[1] * nx[2] - relrot[2] * nx[1]);
vrl[1] = Reff * (relrot[2] * nx[0] - relrot[0] * nx[2]);
vrl[2] = Reff * (relrot[0] * nx[1] - relrot[1] * nx[0]);
}
if (twisting_defined) {
// omega_T (eq 29 of Marshall)
magtwist = dot3(relrot, nx);
}
}
/* ---------------------------------------------------------------------- */
void ContactModel::calculate_forces()
{
// calculate forces/torques
forces[0] = 0.0;
double Fne, Fdamp, dist_to_contact;
area = normal_model->calculate_area();
normal_model->set_knfac();
Fne = normal_model->calculate_forces();
Fdamp = damping_model->calculate_forces();
Fntot = Fne + Fdamp;
if (limit_damping && Fntot < 0.0) Fntot = 0.0;
normal_model->set_fncrit(); // Needed for tangential, rolling, twisting
tangential_model->calculate_forces();
if (rolling_defined) rolling_model->calculate_forces();
if (twisting_defined) twisting_model->calculate_forces();
// sum contributions
scale3(Fntot, nx, forces);
add3(forces, fs, forces);
//May need to rethink eventually tris..
cross3(nx, fs, torquesi);
scale3(-1, torquesi);
if (contact_type == PAIR) copy3(torquesi, torquesj);
if (!classic_model && contact_type == PAIR) {
dist_to_contact = radi - 0.5 * delta;
scale3(dist_to_contact, torquesi);
dist_to_contact = radj - 0.5 * delta;
scale3(dist_to_contact, torquesj);
} else {
dist_to_contact = radi;
scale3(dist_to_contact, torquesi);
}
double torroll[3];
if (rolling_defined) {
cross3(nx, fr, torroll);
scale3(Reff, torroll);
add3(torquesi, torroll, torquesi);
if (contact_type == PAIR) sub3(torquesj, torroll, torquesj);
}
double tortwist[3];
if (twisting_defined) {
scale3(magtortwist, nx, tortwist);
add3(torquesi, tortwist, torquesi);
if (contact_type == PAIR) sub3(torquesj, tortwist, torquesj);
}
if (heat_defined) {
dq = heat_model->calculate_heat();
}
}
/* ----------------------------------------------------------------------
compute pull-off distance (beyond contact) for a given radius and atom type
use temporary variables since this does not use a specific contact geometry
------------------------------------------------------------------------- */
double ContactModel::pulloff_distance(double radi, double radj)
{
return normal_model->pulloff_distance(radi, radj);
}

107
src/contact.h
View File

@ -1,107 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LMP_CONTACT_H
#define LMP_CONTACT_H
#include "pointers.h" // IWYU pragma: export
namespace LAMMPS_NS {
namespace Contact {
#define EPSILON 1e-10
#define NSUBMODELS 6
enum ModelType {
NORMAL = 0,
DAMPING = 1,
TANGENTIAL = 2,
ROLLING = 3,
TWISTING = 4,
HEAT = 5
}; // Relative order matters since some derive coeffs from others
enum ContactType {
PAIR = 0,
WALL = 1,
WALLREGION = 2
};
// forward declaration
class NormalModel;
class DampingModel;
class TangentialModel;
class RollingModel;
class TwistingModel;
class HeatModel;
class SubModel;
class ContactModel : protected Pointers {
public:
ContactModel(class LAMMPS *);
~ContactModel();
void init();
bool check_contact();
void prep_contact();
void calculate_forces();
double pulloff_distance(double, double);
void init_model(std::string, ModelType);
int init_classic_model(char **, int, int);
int mix_coeffs(ContactModel*, ContactModel*);
void write_restart(FILE *);
void read_restart(FILE *);
// Sub models
NormalModel *normal_model;
DampingModel *damping_model;
TangentialModel *tangential_model;
RollingModel *rolling_model;
TwistingModel *twisting_model;
HeatModel *heat_model;
SubModel *sub_models[NSUBMODELS]; // Need to resize if we add more model flavors
// Extra options
int beyond_contact, limit_damping, history_update;
ContactType contact_type;
// History variables
int size_history, nondefault_history_transfer;
double *transfer_history_factor;
double *history;
// Contact properties/output
double forces[3], torquesi[3], torquesj[3], dq;
double radi, radj, meff, dt, Ti, Tj, area;
double Fntot, magtortwist;
int i, j;
double *xi, *xj, *vi, *vj, *omegai, *omegaj;
double fs[3], fr[3], ft[3];
double dx[3], nx[3], r, rsq, rinv, Reff, radsum, delta, dR;
double vr[3], vn[3], vnnr, vt[3], wr[3], vtr[3], vrl[3], relrot[3], vrel;
double magtwist;
bool touch;
protected:
int rolling_defined, twisting_defined, heat_defined; // Used to quickly skip undefined submodels
int classic_model;
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif

113
src/contact_damping_models.cpp
View File

@ -1,113 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "contact_damping_models.h"
#include "contact_normal_models.h"
#include "contact.h"
#include "math_special.h"
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathSpecial;
/* ----------------------------------------------------------------------
Default damping model
------------------------------------------------------------------------- */
DampingModel::DampingModel(LAMMPS *lmp) : SubModel(lmp) {}
/* ---------------------------------------------------------------------- */
void DampingModel::init()
{
damp = contact->normal_model->damp;
}
/* ----------------------------------------------------------------------
No model
------------------------------------------------------------------------- */
DampingNone::DampingNone(LAMMPS *lmp) : DampingModel(lmp) {}
/* ---------------------------------------------------------------------- */
double DampingNone::calculate_forces()
{
return 0.0;
}
/* ----------------------------------------------------------------------
Velocity damping
------------------------------------------------------------------------- */
DampingVelocity::DampingVelocity(LAMMPS *lmp) : DampingModel(lmp) {}
/* ---------------------------------------------------------------------- */
double DampingVelocity::calculate_forces()
{
damp_prefactor = damp;
return -damp_prefactor * contact->vnnr;
}
/* ----------------------------------------------------------------------
Mass velocity damping
------------------------------------------------------------------------- */
DampingMassVelocity::DampingMassVelocity(LAMMPS *lmp) : DampingModel(lmp) {}
/* ---------------------------------------------------------------------- */
double DampingMassVelocity::calculate_forces()
{
damp_prefactor = damp * contact->meff;
return -damp_prefactor * contact->vnnr;
}
/* ----------------------------------------------------------------------
Default, viscoelastic damping
------------------------------------------------------------------------- */
DampingViscoelastic::DampingViscoelastic(LAMMPS *lmp) : DampingModel(lmp) {}
/* ---------------------------------------------------------------------- */
double DampingViscoelastic::calculate_forces()
{
damp_prefactor = damp * contact->meff * contact->area;
return -damp_prefactor * contact->vnnr;
}
/* ----------------------------------------------------------------------
Tsuji damping
------------------------------------------------------------------------- */
DampingTsuji::DampingTsuji(LAMMPS *lmp) : DampingModel(lmp) {}
/* ---------------------------------------------------------------------- */
void DampingTsuji::init()
{
double tmp = contact->normal_model->damp;
damp = 1.2728 - 4.2783 * tmp + 11.087 * square(tmp);
damp += -22.348 * cube(tmp)+ 27.467 * powint(tmp, 4);
damp += -18.022 * powint(tmp, 5) + 4.8218 * powint(tmp,6);
}
/* ---------------------------------------------------------------------- */
double DampingTsuji::calculate_forces()
{
damp_prefactor = damp * sqrt(contact->meff * contact->normal_model->knfac);
return -damp_prefactor * contact->vnnr;
}

79
src/contact_damping_models.h
View File

@ -1,79 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef CONTACT_DAMPING_MODELS_H_
#define CONTACT_DAMPING_MODELS_H_
#include "contact_sub_models.h"
#include "pointers.h"
namespace LAMMPS_NS {
namespace Contact {
class DampingModel : public SubModel {
public:
DampingModel(class LAMMPS *);
~DampingModel() {};
virtual void coeffs_to_local() {};
virtual void mix_coeffs(double*, double*) {};
virtual void init();
virtual double calculate_forces() = 0;
double damp, damp_prefactor;
};
/* ---------------------------------------------------------------------- */
class DampingNone : public DampingModel {
public:
DampingNone(class LAMMPS *);
void init() override {};
double calculate_forces();
};
/* ---------------------------------------------------------------------- */
class DampingVelocity : public DampingModel {
public:
DampingVelocity(class LAMMPS *);
double calculate_forces();
};
/* ---------------------------------------------------------------------- */
class DampingMassVelocity : public DampingModel {
public:
DampingMassVelocity(class LAMMPS *);
double calculate_forces();
};
/* ---------------------------------------------------------------------- */
class DampingViscoelastic : public DampingModel {
public:
DampingViscoelastic(class LAMMPS *);
double calculate_forces();
};
/* ---------------------------------------------------------------------- */
class DampingTsuji : public DampingModel {
public:
DampingTsuji(class LAMMPS *);
void init() override;
double calculate_forces();
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_DAMPING_MODELS_H_ */

63
src/contact_heat_models.cpp
View File

@ -1,63 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "contact_heat_models.h"
#include "contact.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace Contact;
/* ----------------------------------------------------------------------
Default heat conduction
------------------------------------------------------------------------- */
HeatModel::HeatModel(LAMMPS *lmp) : SubModel(lmp) {}
/* ----------------------------------------------------------------------
Area-based heat conduction
------------------------------------------------------------------------- */
HeatNone::HeatNone(LAMMPS *lmp) : HeatModel(lmp) {}
/* ---------------------------------------------------------------------- */
double HeatNone::calculate_heat()
{
return 0.0;
}
/* ----------------------------------------------------------------------
Area-based heat conduction
------------------------------------------------------------------------- */
HeatArea::HeatArea(LAMMPS *lmp) : HeatModel(lmp)
{
num_coeffs = 1;
}
/* ---------------------------------------------------------------------- */
void HeatArea::coeffs_to_local()
{
conductivity = coeffs[0];
if (conductivity < 0.0) error->all(FLERR, "Illegal area heat model");
}
/* ---------------------------------------------------------------------- */
double HeatArea::calculate_heat()
{
return conductivity * contact->area * (contact->Tj - contact->Ti);
}

53
src/contact_heat_models.h
View File

@ -1,53 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef CONTACT_HEAT_MODELS_H_
#define CONTACT_HEAT_MODELS_H_
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class HeatModel : public SubModel {
public:
HeatModel(class LAMMPS *);
~HeatModel() {};
virtual void coeffs_to_local() {};
virtual void init() {};
virtual double calculate_heat() = 0;
};
/* ---------------------------------------------------------------------- */
class HeatNone : public HeatModel {
public:
HeatNone(class LAMMPS *);
double calculate_heat();
};
/* ---------------------------------------------------------------------- */
class HeatArea : public HeatModel {
public:
HeatArea(class LAMMPS *);
void coeffs_to_local() override;
double calculate_heat();
protected:
double conductivity;
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_HEAT_MODELS_H_ */

377
src/contact_normal_models.cpp
View File

@ -1,377 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "contact_normal_models.h"
#include "contact.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathConst;
#define PI27SQ 266.47931882941264802866 // 27*PI**2
#define THREEROOT3 5.19615242270663202362 // 3*sqrt(3)
#define SIXROOT6 14.69693845669906728801 // 6*sqrt(6)
#define INVROOT6 0.40824829046386307274 // 1/sqrt(6)
#define FOURTHIRDS (4.0/3.0) // 4/3
#define ONETHIRD (1.0/3.0) // 1/3
#define THREEQUARTERS 0.75 // 3/4
/* ----------------------------------------------------------------------
Default normal model
------------------------------------------------------------------------- */
NormalModel::NormalModel(LAMMPS *lmp) : SubModel(lmp)
{
material_properties = 0;
}
/* ---------------------------------------------------------------------- */
bool NormalModel::touch()
{
bool touchflag = (contact->rsq < contact->radsum * contact->radsum);
return touchflag;
}
/* ---------------------------------------------------------------------- */
double NormalModel::pulloff_distance(double radi, double radj)
{
//called outside of compute(), do not assume correct geometry defined in contact
return radi + radj;
}
/* ---------------------------------------------------------------------- */
double NormalModel::calculate_area()
{
return sqrt(contact->dR);
}
/* ---------------------------------------------------------------------- */
void NormalModel::set_fncrit()
{
Fncrit = fabs(contact->Fntot);
}
/* ----------------------------------------------------------------------
No model
------------------------------------------------------------------------- */
NormalNone::NormalNone(LAMMPS *lmp) : NormalModel(lmp) {}
/* ---------------------------------------------------------------------- */
double NormalNone::calculate_forces()
{
return 0.0;
}
/* ----------------------------------------------------------------------
Hookean normal force
------------------------------------------------------------------------- */
NormalHooke::NormalHooke(LAMMPS *lmp) : NormalModel(lmp)
{
num_coeffs = 2;
}
/* ---------------------------------------------------------------------- */
void NormalHooke::coeffs_to_local()
{
k = coeffs[0];
damp = coeffs[1];
if (k < 0.0 || damp < 0.0) error->all(FLERR, "Illegal Hooke normal model");
}
/* ---------------------------------------------------------------------- */
double NormalHooke::calculate_forces()
{
Fne = knfac * contact->delta;
return Fne;
}
/* ---------------------------------------------------------------------- */
void NormalHooke::set_knfac()
{
knfac = k;
}
/* ----------------------------------------------------------------------
Hertzian normal force
------------------------------------------------------------------------- */
NormalHertz::NormalHertz(LAMMPS *lmp) : NormalModel(lmp)
{
num_coeffs = 2;
}
/* ---------------------------------------------------------------------- */
void NormalHertz::coeffs_to_local()
{
k = coeffs[0];
damp = coeffs[1];
if (k < 0.0 || damp < 0.0) error->all(FLERR, "Illegal Hertz normal model");
}
/* ---------------------------------------------------------------------- */
double NormalHertz::calculate_forces()
{
Fne = knfac * contact->delta;
return Fne;
}
/* ---------------------------------------------------------------------- */
void NormalHertz::set_knfac()
{
knfac = k * contact->area;
}
/* ----------------------------------------------------------------------
Hertzian normal force with material properties
------------------------------------------------------------------------- */
NormalHertzMaterial::NormalHertzMaterial(LAMMPS *lmp) : NormalHertz(lmp)
{
material_properties = 1;
num_coeffs = 3;
}
/* ---------------------------------------------------------------------- */
void NormalHertzMaterial::coeffs_to_local()
{
Emod = coeffs[0];
damp = coeffs[1];
poiss = coeffs[2];
if (contact->contact_type == PAIR) {
k = FOURTHIRDS * mix_stiffnessE(Emod, Emod, poiss, poiss);
} else {
k = FOURTHIRDS * mix_stiffnessE_wall(Emod, poiss);
}
if (Emod < 0.0 || damp < 0.0) error->all(FLERR, "Illegal Hertz material normal model");
}
/* ---------------------------------------------------------------------- */
void NormalHertzMaterial::mix_coeffs(double* icoeffs, double* jcoeffs)
{
coeffs[0] = mix_stiffnessE(icoeffs[0], jcoeffs[0],icoeffs[2], jcoeffs[2]);
coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]);
coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]);
coeffs_to_local();
}
/* ----------------------------------------------------------------------
DMT normal force
------------------------------------------------------------------------- */
NormalDMT::NormalDMT(LAMMPS *lmp) : NormalModel(lmp)
{
allow_limit_damping = 0;
material_properties = 1;
num_coeffs = 4;
}
/* ---------------------------------------------------------------------- */
void NormalDMT::coeffs_to_local()
{
Emod = coeffs[0];
damp = coeffs[1];
poiss = coeffs[2];
cohesion = coeffs[3];
if (contact->contact_type == PAIR) {
k = FOURTHIRDS * mix_stiffnessE(Emod, Emod, poiss, poiss);
} else {
k = FOURTHIRDS * mix_stiffnessE_wall(Emod, poiss);
}
if (Emod < 0.0 || damp < 0.0) error->all(FLERR, "Illegal DMT normal model");
}
/* ---------------------------------------------------------------------- */
void NormalDMT::mix_coeffs(double* icoeffs, double* jcoeffs)
{
coeffs[0] = mix_stiffnessE(icoeffs[0], jcoeffs[0],icoeffs[2], jcoeffs[2]);
coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]);
coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]);
coeffs[3] = mix_geom(icoeffs[3], jcoeffs[3]);
coeffs_to_local();
}
/* ---------------------------------------------------------------------- */
double NormalDMT::calculate_forces()
{
Fne = knfac * contact->delta;
F_pulloff = 4.0 * MathConst::MY_PI * cohesion * contact->Reff;
Fne -= F_pulloff;
return Fne;
}
/* ---------------------------------------------------------------------- */
void NormalDMT::set_knfac()
{
knfac = k * contact->area;
}
/* ---------------------------------------------------------------------- */
void NormalDMT::set_fncrit()
{
Fncrit = fabs(Fne + 2.0 * F_pulloff);
}
/* ----------------------------------------------------------------------
JKR normal force
------------------------------------------------------------------------- */
NormalJKR::NormalJKR(LAMMPS *lmp) : NormalModel(lmp)
{
allow_limit_damping = 0;
material_properties = 1;
beyond_contact = 1;
num_coeffs = 4;
}
/* ---------------------------------------------------------------------- */
void NormalJKR::coeffs_to_local()
{
Emod = coeffs[0];
damp = coeffs[1];
poiss = coeffs[2];
cohesion = coeffs[3];
if (contact->contact_type == PAIR) {
Emix = mix_stiffnessE(Emod, Emod, poiss, poiss);
} else {
Emix = mix_stiffnessE_wall(Emod, poiss);
}
k = FOURTHIRDS * Emix;
if (Emod < 0.0 || damp < 0.0) error->all(FLERR, "Illegal JKR normal model");
}
/* ---------------------------------------------------------------------- */
void NormalJKR::mix_coeffs(double* icoeffs, double* jcoeffs)
{
coeffs[0] = mix_stiffnessE(icoeffs[0], jcoeffs[0],icoeffs[2], jcoeffs[2]);
coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]);
coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]);
coeffs[3] = mix_geom(icoeffs[3], jcoeffs[3]);
coeffs_to_local();
}
/* ---------------------------------------------------------------------- */
bool NormalJKR::touch()
{
double area_at_pulloff, R2, delta_pulloff, dist_pulloff;
bool touchflag;
if (contact->touch) {
R2 = contact->Reff * contact->Reff;
area_at_pulloff = cbrt(9.0 * MY_PI * cohesion * R2 / (4.0 * Emix));
delta_pulloff = area_at_pulloff * area_at_pulloff / contact->Reff - 2.0 * sqrt(MY_PI * cohesion * area_at_pulloff / Emix);
dist_pulloff = contact->radsum - delta_pulloff;
touchflag = contact->rsq < (dist_pulloff * dist_pulloff);
} else {
touchflag = contact->rsq < (contact->radsum * contact->radsum);
}
return touchflag;
}
/* ----------------------------------------------------------------------
called outside of compute(), do not assume geometry defined in contact
------------------------------------------------------------------------- */
double NormalJKR::pulloff_distance(double radi, double radj)
{
double area_at_pulloff, Reff_tmp;
Reff_tmp = radi * radj / (radi + radj); // May not be defined
if (Reff_tmp <= 0) return 0;
area_at_pulloff = cbrt(9.0 * MY_PI * cohesion * Reff_tmp * Reff_tmp / (4.0 * Emix));
return area_at_pulloff * area_at_pulloff / Reff_tmp - 2.0 * sqrt(MY_PI * cohesion * area_at_pulloff / Emix);
}
/* ---------------------------------------------------------------------- */
double NormalJKR::calculate_area()
{
double R2, dR2, t0, t1, t2, t3, t4, t5, t6;
double sqrt1, sqrt2, sqrt3;
R2 = contact->Reff * contact->Reff;
dR2 = contact->dR * contact->dR;
t0 = cohesion * cohesion * R2 * R2 * Emix;
t1 = PI27SQ * t0;
t2 = 8.0 * contact->dR * dR2 * Emix * Emix * Emix;
t3 = 4.0 * dR2 * Emix;
// in case sqrt(0) < 0 due to precision issues
sqrt1 = MAX(0, t0 * (t1 + 2.0 * t2));
t4 = cbrt(t1 + t2 + THREEROOT3 * MY_PI * sqrt(sqrt1));
t5 = t3 / t4 + t4 / Emix;
sqrt2 = MAX(0, 2.0 * contact->dR + t5);
t6 = sqrt(sqrt2);
sqrt3 = MAX(0, 4.0 * contact->dR - t5 + SIXROOT6 * cohesion * MY_PI * R2 / (Emix * t6));
return INVROOT6 * (t6 + sqrt(sqrt3));
}
/* ---------------------------------------------------------------------- */
double NormalJKR::calculate_forces()
{
double a2;
a2 = contact->area * contact->area;
Fne = knfac * a2 / contact->Reff - MY_2PI * a2 * sqrt(4.0 * cohesion * Emix / (MY_PI * contact->area));
F_pulloff = 3.0 * MY_PI * cohesion * contact->Reff;
return Fne;
}
/* ---------------------------------------------------------------------- */
void NormalJKR::set_knfac()
{
knfac = k * contact->area;
}
/* ---------------------------------------------------------------------- */
void NormalJKR::set_fncrit()
{
Fncrit = fabs(Fne + 2.0 * F_pulloff);
}

118
src/contact_normal_models.h
View File

@ -1,118 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef CONTACT_NORMAL_MODELS_H_
#define CONTACT_NORMAL_MODELS_H_
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class NormalModel : public SubModel {
public:
NormalModel(class LAMMPS *);
~NormalModel() {};
virtual void coeffs_to_local() {};
virtual void init() {};
virtual bool touch();
virtual double pulloff_distance(double, double);
virtual double calculate_area();
virtual void set_knfac() = 0;
virtual double calculate_forces() = 0;
virtual void set_fncrit();
double damp; // Vestigial argument needed by damping
double Emod, poiss;
double Fncrit, Fne, knfac;
int material_properties;
};
/* ---------------------------------------------------------------------- */
class NormalNone : public NormalModel {
public:
NormalNone(class LAMMPS *);
void set_knfac() {};
double calculate_forces();
};
/* ---------------------------------------------------------------------- */
class NormalHooke : public NormalModel {
public:
NormalHooke(class LAMMPS *);
void coeffs_to_local() override;
void set_knfac();
double calculate_forces();
protected:
double k;
};
/* ---------------------------------------------------------------------- */
class NormalHertz : public NormalModel {
public:
NormalHertz(class LAMMPS *);
void coeffs_to_local() override;
void set_knfac();
double calculate_forces();
protected:
double k;
};
/* ---------------------------------------------------------------------- */
class NormalHertzMaterial : public NormalHertz {
public:
NormalHertzMaterial(class LAMMPS *);
void coeffs_to_local() override;
void mix_coeffs(double*, double*) override;
};
/* ---------------------------------------------------------------------- */
class NormalDMT : public NormalModel {
public:
NormalDMT(class LAMMPS *);
void coeffs_to_local() override;
void mix_coeffs(double*, double*) override;
void set_knfac();
double calculate_forces();
void set_fncrit() override;
protected:
double k, cohesion;
double F_pulloff;
};
/* ---------------------------------------------------------------------- */
class NormalJKR : public NormalModel {
public:
NormalJKR(class LAMMPS *);
void coeffs_to_local() override;
void mix_coeffs(double*, double*) override;
bool touch() override;
double pulloff_distance(double, double) override;
double calculate_area() override;
void set_knfac();
double calculate_forces();
void set_fncrit() override;
protected:
double k, cohesion;
double Emix, F_pulloff;
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_NORMAL_MODELS_H_ */

121
src/contact_rolling_models.cpp
View File

@ -1,121 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "contact_normal_models.h"
#include "contact_rolling_models.h"
#include "contact.h"
#include "error.h"
#include "math_const.h"
#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathConst;
using namespace MathExtra;
/* ----------------------------------------------------------------------
Default rolling friction model
------------------------------------------------------------------------- */
RollingModel::RollingModel(LAMMPS *lmp) : SubModel(lmp) {}
/* ----------------------------------------------------------------------
No model
------------------------------------------------------------------------- */
RollingNone::RollingNone(LAMMPS *lmp) : RollingModel(lmp) {}
/* ----------------------------------------------------------------------
SDS rolling friction model
------------------------------------------------------------------------- */
RollingSDS::RollingSDS(LAMMPS *lmp) : RollingModel(lmp)
{
num_coeffs = 3;
size_history = 3;
}
/* ---------------------------------------------------------------------- */
void RollingSDS::coeffs_to_local()
{
k = coeffs[0];
gamma = coeffs[1];
mu = coeffs[2];
if (k < 0.0 || mu < 0.0 || gamma < 0.0)
error->all(FLERR, "Illegal SDS rolling model");
}
/* ---------------------------------------------------------------------- */
void RollingSDS::calculate_forces()
{
int rhist0, rhist1, rhist2, frameupdate;
double Frcrit, rolldotn, rollmag, prjmag, magfr, hist_temp[3], scalefac, temp_array[3];
double k_inv, magfr_inv;
rhist0 = history_index;
rhist1 = rhist0 + 1;
rhist2 = rhist1 + 1;
Frcrit = mu * contact->normal_model->Fncrit;
if (contact->history_update) {
hist_temp[0] = contact->history[rhist0];
hist_temp[1] = contact->history[rhist1];
hist_temp[2] = contact->history[rhist2];
rolldotn = dot3(hist_temp, contact->nx);
frameupdate = (fabs(rolldotn) * k) > (EPSILON * Frcrit);
if (frameupdate) { // rotate into tangential plane
rollmag = len3(hist_temp);
// projection
scale3(rolldotn, contact->nx, temp_array);
sub3(hist_temp, temp_array, hist_temp);
// also rescale to preserve magnitude
prjmag = len3(hist_temp);
if (prjmag > 0) scalefac = rollmag / prjmag;
else scalefac = 0;
scale3(scalefac, hist_temp);
}
scale3(contact->dt, contact->vrl, temp_array);
add3(hist_temp, temp_array, hist_temp);
}
scaleadd3(-k, hist_temp, -gamma, contact->vrl, contact->fr);
// rescale frictional displacements and forces if needed
magfr = len3(contact->fr);
if (magfr > Frcrit) {
rollmag = len3(hist_temp);
if (rollmag != 0.0) {
k_inv = 1.0 / k;
magfr_inv = 1.0 / magfr;
scale3(-Frcrit * k_inv * magfr_inv, contact->fr, hist_temp);
scale3(-gamma * k_inv, contact->vrl, temp_array);
add3(hist_temp, temp_array, hist_temp);
scale3(Frcrit * magfr_inv, contact->fr);
} else {
zero3(contact->fr);
}
}
if (contact->history_update) {
contact->history[rhist0] = hist_temp[0];
contact->history[rhist1] = hist_temp[1];
contact->history[rhist2] = hist_temp[2];
}
}

53
src/contact_rolling_models.h
View File

@ -1,53 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef CONTACT_ROLLING_MODELS_H_
#define CONTACT_ROLLING_MODELS_H_
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class RollingModel : public SubModel {
public:
RollingModel(class LAMMPS *);
~RollingModel() {};
virtual void coeffs_to_local() {};
virtual void init() {};
virtual void calculate_forces() = 0;
};
/* ---------------------------------------------------------------------- */
class RollingNone : public RollingModel {
public:
RollingNone(class LAMMPS *);
void calculate_forces() {};
};
/* ---------------------------------------------------------------------- */
class RollingSDS : public RollingModel {
public:
RollingSDS(class LAMMPS *);
void coeffs_to_local() override;
void calculate_forces();
protected:
double k, mu, gamma;
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_ROLLING_MODELS_H_ */

140
src/contact_sub_models.cpp
View File

@ -1,140 +0,0 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------*/
/* ----------------------------------------------------------------------
This class contains a framework for normal, damping, tangential,
rolling, twisting, and heat models used to calculate forces
and torques based on contact geometry
Contributing authors:
Dan Bolintineanu (SNL), Joel Clemmer (SNL)
----------------------------------------------------------------------- */
#include "contact_sub_models.h"
#include "error.h"
#include "utils.h"
using namespace LAMMPS_NS;
using namespace Contact;
/* ----------------------------------------------------------------------
Parent class for all types of contact forces
------------------------------------------------------------------------- */
SubModel::SubModel(LAMMPS *lmp) : Pointers(lmp)
{
allocated = 0;
size_history = 0;
history_index = 0;
allow_limit_damping = 1;
beyond_contact = 0;
num_coeffs = 0;
nondefault_history_transfer = 0;
transfer_history_factor = nullptr;
}
/* ---------------------------------------------------------------------- */
SubModel::~SubModel()
{
if (allocated) delete [] coeffs;
delete [] transfer_history_factor;
}
/* ---------------------------------------------------------------------- */
void SubModel::allocate_coeffs()
{
allocated = 1;
coeffs = new double[num_coeffs];
}
/* ---------------------------------------------------------------------- */
int SubModel::parse_coeffs(char **arg, int iarg, int narg)
{
if (iarg + num_coeffs > narg)
error->all(FLERR, "Insufficient arguments provided for {} model", name);
for (int i = 0; i < num_coeffs; i++) {
// A few parameters (eg.g kt for tangential mindlin) allow null
if (strcmp(arg[iarg+i], "NULL") == 0) coeffs[i] = -1;
else coeffs[i] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
}
coeffs_to_local();
return iarg + num_coeffs;
}
/* ---------------------------------------------------------------------- */
void SubModel::mix_coeffs(double* icoeffs, double* jcoeffs)
{
for (int i = 0; i < num_coeffs; i++)
coeffs[i] = mix_geom(icoeffs[i], jcoeffs[i]);
coeffs_to_local();
}
/* ----------------------------------------------------------------------
mixing of Young's modulus (E)
------------------------------------------------------------------------- */
double SubModel::mix_stiffnessE(double E1, double E2,
double pois1, double pois2)
{
double factor1 = (1 - pois1 * pois1) / E1;
double factor2 = (1 - pois2 * pois2) / E2;
return 1 / (factor1 + factor2);
}
/* ----------------------------------------------------------------------
mixing of shear modulus (G)
------------------------------------------------------------------------ */
double SubModel::mix_stiffnessG(double E1, double E2,
double pois1, double pois2)
{
double factor1 = 2 * (2 - pois1) * (1 + pois1) / E1;
double factor2 = 2 * (2 - pois2) * (1 + pois2) / E2;
return 1 / (factor1 + factor2);
}
/* ----------------------------------------------------------------------
mixing of Young's modulus (E) for walls
------------------------------------------------------------------------- */
double SubModel::mix_stiffnessE_wall(double E, double pois)
{
double factor = 2 * (1 - pois);
return E / factor;
}
/* ----------------------------------------------------------------------
mixing of shear modulus (G) for walls
------------------------------------------------------------------------ */
double SubModel::mix_stiffnessG_wall(double E, double pois)
{
double factor = 32.0 * (2 - pois) * (1 + pois);
return E / factor;
}
/* ----------------------------------------------------------------------
mixing of everything else
------------------------------------------------------------------------- */
double SubModel::mix_geom(double val1, double val2)
{
return sqrt(val1 * val2);
}

62
src/contact_sub_models.h
View File

@ -1,62 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef CONTACT_SUB_MODEL_H_
#define CONTACT_SUB_MODEL_H_
#include "contact.h"
#include "pointers.h" // IWYU pragma: export
namespace LAMMPS_NS {
namespace Contact {
class SubModel : protected Pointers {
public:
SubModel(class LAMMPS *);
virtual ~SubModel();
int num_coeffs;
double *coeffs;
void read_restart();
int parse_coeffs(char **, int, int);
virtual void mix_coeffs(double*, double*);
virtual void coeffs_to_local() = 0;
virtual void init() = 0; // called after all other submodel coeffs defined
void allocate_coeffs();
std::string name;
int size_history;
int nondefault_history_transfer;
double *transfer_history_factor;
int history_index;
int beyond_contact;
int allow_limit_damping;
ContactModel *contact;
protected:
int allocated;
double mix_stiffnessE(double, double, double, double);
double mix_stiffnessG(double, double, double, double);
double mix_stiffnessE_wall(double, double);
double mix_stiffnessG_wall(double, double);
double mix_geom(double, double);
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif /* CONTACT_SUB_MODEL_H_ */

412
src/contact_tangential_models.cpp
View File

@ -1,412 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "contact_damping_models.h"
#include "contact_normal_models.h"
#include "contact_tangential_models.h"
#include "contact.h"
#include "error.h"
#include "math_const.h"
#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathConst;
using namespace MathExtra;
/* ----------------------------------------------------------------------
Default model
------------------------------------------------------------------------- */
TangentialModel::TangentialModel(LAMMPS *lmp) : SubModel(lmp) {}
/* ----------------------------------------------------------------------
No model
------------------------------------------------------------------------- */
TangentialNone::TangentialNone(LAMMPS *lmp) : TangentialModel(lmp) {}
/* ----------------------------------------------------------------------
Linear model with no history
------------------------------------------------------------------------- */
TangentialLinearNoHistory::TangentialLinearNoHistory(LAMMPS *lmp) : TangentialModel(lmp)
{
num_coeffs = 2;
size_history = 3;
}
/* ---------------------------------------------------------------------- */
void TangentialLinearNoHistory::coeffs_to_local()
{
k = 0.0; // No tangential stiffness with no history
xt = coeffs[0];
mu = coeffs[1];
if (k < 0.0 || xt < 0.0 || mu < 0.0)
error->all(FLERR, "Illegal linear no history tangential model");
}
/* ---------------------------------------------------------------------- */
void TangentialLinearNoHistory::calculate_forces()
{
// classic pair gran/hooke (no history)
damp = xt * contact->damping_model->damp_prefactor;
double Fscrit = mu * contact->normal_model->Fncrit;
double fsmag = damp * contact->vrel;
double Ft;
if (contact->vrel != 0.0) Ft = MIN(Fscrit, fsmag) / contact->vrel;
else Ft = 0.0;
scale3(-Ft, contact->vtr, contact->fs);
}
/* ----------------------------------------------------------------------
Linear model with history
------------------------------------------------------------------------- */
TangentialLinearHistory::TangentialLinearHistory(LAMMPS *lmp) : TangentialModel(lmp)
{
num_coeffs = 3;
size_history = 3;
}
/* ---------------------------------------------------------------------- */
void TangentialLinearHistory::coeffs_to_local()
{
k = coeffs[0];
xt = coeffs[1];
mu = coeffs[2];
if (k < 0.0 || xt < 0.0 || mu < 0.0)
error->all(FLERR, "Illegal linear tangential model");
}
/* ---------------------------------------------------------------------- */
void TangentialLinearHistory::calculate_forces()
{
// Note: this is the same as the base Mindlin calculation except k isn't scaled by area
double magfs, magfs_inv, rsht, shrmag, prjmag, temp_dbl, temp_array[3];
int frame_update = 0;
damp = xt * contact->damping_model->damp_prefactor;
double Fscrit = contact->normal_model->Fncrit * mu;
double *history = & contact->history[history_index];
// rotate and update displacements / force.
// see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
if (contact->history_update) {
rsht = dot3(history, contact->nx);
frame_update = (fabs(rsht) * k) > (EPSILON * Fscrit);
if (frame_update) {
shrmag = len3(history);
// projection
scale3(rsht, contact->nx, temp_array);
sub3(history, temp_array, history);
// also rescale to preserve magnitude
prjmag = len3(history);
if (prjmag > 0) temp_dbl = shrmag / prjmag;
else temp_dbl = 0;
scale3(temp_dbl, history);
}
// update history, tangential force
// see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
scale3(contact->dt, contact->vtr, temp_array);
add3(history, temp_array, history);
}
// tangential forces = history + tangential velocity damping
scale3(-k, history, contact->fs);
scale3(damp, contact->vtr, temp_array);
sub3(contact->fs, temp_array, contact->fs);
// rescale frictional displacements and forces if needed
magfs = len3(contact->fs);
if (magfs > Fscrit) {
shrmag = len3(history);
if (shrmag != 0.0) {
magfs_inv = 1.0 / magfs;
scale3(Fscrit * magfs_inv, contact->fs, history);
scale3(damp, contact->vtr, temp_array);
add3(history, temp_array, history);
scale3(-1.0 / k, history);
scale3(Fscrit * magfs_inv, contact->fs);
} else {
zero3(contact->fs);
}
}
}
/* ----------------------------------------------------------------------
Linear model with history from pair gran/hooke/history
------------------------------------------------------------------------- */
TangentialLinearHistoryClassic::TangentialLinearHistoryClassic(LAMMPS *lmp) : TangentialLinearHistory(lmp)
{
scale_area = 0; // Sets gran/hooke/history behavior
}
/* ---------------------------------------------------------------------- */
void TangentialLinearHistoryClassic::calculate_forces()
{
double k_scaled, magfs, magfs_inv, rsht, shrmag, prjmag, temp_dbl;
double temp_array[3];
int frame_update = 0;
k_scaled = k;
if (scale_area) k_scaled *= contact->area;
damp = xt * contact->damping_model->damp_prefactor;
double Fscrit = contact->normal_model->Fncrit * mu;
double *history = & contact->history[history_index];
// update history
if (contact->history_update) {
scale3(contact->dt, contact->vtr, temp_array);
add3(history, temp_array, history);
}
shrmag = len3(history);
// rotate shear displacements
if (contact->history_update) {
rsht = dot3(history, contact->nx);
scale3(rsht, contact->nx, temp_array);
sub3(history, temp_array, history);
}
// tangential forces = history + tangential velocity damping
scale3(-k_scaled, history, contact->fs);
scale3(damp, contact->vtr, temp_array);
sub3(contact->fs, temp_array, contact->fs);
// rescale frictional displacements and forces if needed
magfs = len3(contact->fs);
if (magfs > Fscrit) {
if (shrmag != 0.0) {
magfs_inv = 1.0 / magfs;
scale3(Fscrit * magfs_inv, contact->fs, history);
scale3(damp, contact->vtr, temp_array);
add3(history, temp_array, history);
temp_dbl = -1.0 / k_scaled;
if (scale_area) temp_dbl /= contact->area;
scale3(temp_dbl, history);
scale3(Fscrit * magfs_inv, contact->fs);
} else {
zero3(contact->fs);
}
}
}
/* ----------------------------------------------------------------------
Mindlin from pair gran/hertz/history
------------------------------------------------------------------------- */
TangentialMindlinClassic::TangentialMindlinClassic(LAMMPS *lmp) : TangentialLinearHistoryClassic(lmp)
{
scale_area = 1; // Sets gran/hertz/history behavior
}
/* ----------------------------------------------------------------------
Mindlin model
------------------------------------------------------------------------- */
TangentialMindlin::TangentialMindlin(LAMMPS *lmp) : TangentialModel(lmp)
{
num_coeffs = 3;
size_history = 3;
mindlin_force = 0;
mindlin_rescale = 0;
}
/* ---------------------------------------------------------------------- */
void TangentialMindlin::coeffs_to_local()
{
k = coeffs[0];
xt = coeffs[1];
mu = coeffs[2];
if (k == -1) {
if (!contact->normal_model->material_properties)
error->all(FLERR, "Must either specify tangential stiffness or material properties for normal model for the Mindlin tangential style");
double Emod = contact->normal_model->Emod;
double poiss = contact->normal_model->poiss;
if (contact->contact_type == PAIR) {
k = 8.0 * mix_stiffnessG(Emod, Emod, poiss, poiss);
} else {
k = 8.0 * mix_stiffnessG_wall(Emod, poiss);
}
}
if (k < 0.0 || xt < 0.0 || mu < 0.0)
error->all(FLERR, "Illegal Mindlin tangential model");
}
/* ---------------------------------------------------------------------- */
void TangentialMindlin::mix_coeffs(double* icoeffs, double* jcoeffs)
{
if (icoeffs[0] == -1 || jcoeffs[0] == -1) coeffs[0] = -1;
else coeffs[0] = mix_geom(icoeffs[0], jcoeffs[0]);
coeffs[1] = mix_geom(icoeffs[1], jcoeffs[1]);
coeffs[2] = mix_geom(icoeffs[2], jcoeffs[2]);
coeffs_to_local();
}
/* ---------------------------------------------------------------------- */
void TangentialMindlin::calculate_forces()
{
double k_scaled, magfs, magfs_inv, rsht, shrmag, prjmag, temp_dbl;
double temp_array[3];
int frame_update = 0;
damp = xt * contact->damping_model->damp_prefactor;
double *history = & contact->history[history_index];
double Fscrit = contact->normal_model->Fncrit * mu;
k_scaled = k * contact->area;
// on unloading, rescale the shear displacements/force
if (mindlin_rescale)
if (contact->area < history[3])
scale3(contact->area / history[3], history);
// rotate and update displacements / force.
// see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
if (contact->history_update) {
rsht = dot3(history, contact->nx);
if (mindlin_force) {
frame_update = fabs(rsht) > (EPSILON * Fscrit);
} else {
frame_update = (fabs(rsht) * k_scaled) > (EPSILON * Fscrit);
}
if (frame_update) {
shrmag = len3(history);
// projection
scale3(rsht, contact->nx, temp_array);
sub3(history, temp_array, history);
// also rescale to preserve magnitude
prjmag = len3(history);
if (prjmag > 0) temp_dbl = shrmag / prjmag;
else temp_dbl = 0;
scale3(temp_dbl, history);
}
// update history
if (mindlin_force) {
// tangential force
// see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
scale3(-k_scaled * contact->dt, contact->vtr, temp_array);
} else {
scale3(contact->dt, contact->vtr, temp_array);
}
add3(history, temp_array, history);
if (mindlin_rescale) history[3] = contact->area;
}
// tangential forces = history + tangential velocity damping
scale3(-damp, contact->vtr, contact->fs);
if (!mindlin_force) {
scale3(k_scaled, history, temp_array);
sub3(contact->fs, temp_array, contact->fs);
} else {
add3(contact->fs, history, contact->fs);
}
// rescale frictional displacements and forces if needed
magfs = len3(contact->fs);
if (magfs > Fscrit) {
shrmag = len3(history);
if (shrmag != 0.0) {
magfs_inv = 1.0 / magfs;
scale3(Fscrit * magfs_inv, contact->fs, history);
scale3(damp, contact->vtr, temp_array);
add3(history, temp_array, history);
if (!mindlin_force)
scale3(-1.0 / k_scaled, history);
scale3(Fscrit * magfs_inv, contact->fs);
} else {
zero3(contact->fs);
}
}
}
/* ----------------------------------------------------------------------
Mindlin force model
------------------------------------------------------------------------- */
TangentialMindlinForce::TangentialMindlinForce(LAMMPS *lmp) : TangentialMindlin(lmp)
{
num_coeffs = 3;
size_history = 3;
mindlin_force = 1;
mindlin_rescale = 0;
}
/* ----------------------------------------------------------------------
Mindlin rescale model
------------------------------------------------------------------------- */
TangentialMindlinRescale::TangentialMindlinRescale(LAMMPS *lmp) : TangentialMindlin(lmp)
{
num_coeffs = 3;
size_history = 4;
mindlin_force = 0;
mindlin_rescale = 1;
nondefault_history_transfer = 1;
transfer_history_factor = new double[size_history];
for (int i = 0; i < size_history; i++) transfer_history_factor[i] = -1.0;
transfer_history_factor[3] = +1;
}
/* ----------------------------------------------------------------------
Mindlin rescale force model
------------------------------------------------------------------------- */
TangentialMindlinRescaleForce::TangentialMindlinRescaleForce(LAMMPS *lmp) : TangentialMindlin(lmp)
{
num_coeffs = 3;
size_history = 4;
mindlin_force = 1;
mindlin_rescale = 1;
nondefault_history_transfer = 1;
transfer_history_factor = new double[size_history];
for (int i = 0; i < size_history; i++) transfer_history_factor[i] = -1.0;
transfer_history_factor[3] = +1;
}

113
src/contact_tangential_models.h
View File

@ -1,113 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef CONTACT_TANGENTIAL_MODELS_H_
#define CONTACT_TANGENTIAL_MODELS_H_
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class TangentialModel : public SubModel {
public:
TangentialModel(class LAMMPS *);
virtual ~TangentialModel() {};
virtual void coeffs_to_local() {};
virtual void init() {};
virtual void calculate_forces() = 0;
int rescale_flag;
double k, damp, mu; // Used by Marshall twisting model
protected:
double xt;
};
/* ---------------------------------------------------------------------- */
class TangentialNone : public TangentialModel {
public:
TangentialNone(class LAMMPS *);
void calculate_forces() {};
};
/* ---------------------------------------------------------------------- */
class TangentialLinearNoHistory : public TangentialModel {
public:
TangentialLinearNoHistory(class LAMMPS *);
void coeffs_to_local() override;
void calculate_forces();
};
/* ---------------------------------------------------------------------- */
class TangentialLinearHistory : public TangentialModel {
public:
TangentialLinearHistory(class LAMMPS *);
void coeffs_to_local() override;
void calculate_forces();
};
/* ---------------------------------------------------------------------- */
class TangentialLinearHistoryClassic : public TangentialLinearHistory {
public:
TangentialLinearHistoryClassic(class LAMMPS *);
void calculate_forces();
int scale_area;
};
/* ---------------------------------------------------------------------- */
class TangentialMindlinClassic : public TangentialLinearHistoryClassic {
public:
TangentialMindlinClassic(class LAMMPS *);
};
/* ---------------------------------------------------------------------- */
class TangentialMindlin : public TangentialModel {
public:
TangentialMindlin(class LAMMPS *);
void coeffs_to_local() override;
void mix_coeffs(double*, double*) override;
void calculate_forces();
protected:
int mindlin_rescale, mindlin_force;
};
/* ---------------------------------------------------------------------- */
class TangentialMindlinForce : public TangentialMindlin {
public:
TangentialMindlinForce(class LAMMPS *);
};
/* ---------------------------------------------------------------------- */
class TangentialMindlinRescale : public TangentialMindlin {
public:
TangentialMindlinRescale(class LAMMPS *);
};
/* ---------------------------------------------------------------------- */
class TangentialMindlinRescaleForce : public TangentialMindlin {
public:
TangentialMindlinRescaleForce(class LAMMPS *);
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_TANGENTIAL_MODELS_H_ */

124
src/contact_twisting_models.cpp
View File

@ -1,124 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "contact_normal_models.h"
#include "contact_tangential_models.h"
#include "contact_twisting_models.h"
#include "contact.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace Contact;
using namespace MathConst;
/* ----------------------------------------------------------------------
Default twisting model
------------------------------------------------------------------------- */
TwistingModel::TwistingModel(LAMMPS *lmp) : SubModel(lmp) {}
/* ----------------------------------------------------------------------
No model
------------------------------------------------------------------------- */
TwistingNone::TwistingNone(LAMMPS *lmp) : TwistingModel(lmp) {}
/* ----------------------------------------------------------------------
Marshall twisting model
------------------------------------------------------------------------- */
TwistingMarshall::TwistingMarshall(LAMMPS *lmp) : TwistingModel(lmp)
{
num_coeffs = 0;
size_history = 3;
}
/* ---------------------------------------------------------------------- */
void TwistingMarshall::init()
{
k_tang = contact->tangential_model->k;
mu_tang = contact->tangential_model->mu;
}
/* ---------------------------------------------------------------------- */
void TwistingMarshall::calculate_forces()
{
double signtwist, Mtcrit;
// Calculate twist coefficients from tangential model & contact geometry
// eq 32 of Marshall paper
double k = 0.5 * k_tang * contact->area * contact->area;
double damp = 0.5 * contact->tangential_model->damp * contact->area * contact->area;
double mu = TWOTHIRDS * mu_tang * contact->area;
if (contact->history_update) {
contact->history[history_index] += contact->magtwist * contact->dt;
}
// M_t torque (eq 30)
contact->magtortwist = -k * contact->history[history_index] - damp * contact->magtwist;
signtwist = (contact->magtwist > 0) - (contact->magtwist < 0);
Mtcrit = mu * contact->normal_model->Fncrit; // critical torque (eq 44)
if (fabs(contact->magtortwist) > Mtcrit) {
contact->history[history_index] = (Mtcrit * signtwist - damp * contact->magtwist) / k;
contact->magtortwist = -Mtcrit * signtwist; // eq 34
}
}
/* ----------------------------------------------------------------------
SDS twisting model
------------------------------------------------------------------------- */
TwistingSDS::TwistingSDS(LAMMPS *lmp) : TwistingModel(lmp)
{
num_coeffs = 3;
size_history = 3;
}
/* ---------------------------------------------------------------------- */
void TwistingSDS::coeffs_to_local()
{
k = coeffs[0];
damp = coeffs[1];
mu = coeffs[2];
if (k < 0.0 || mu < 0.0 || damp < 0.0)
error->all(FLERR, "Illegal SDS twisting model");
}
/* ---------------------------------------------------------------------- */
void TwistingSDS::calculate_forces()
{
double signtwist, Mtcrit;
if (contact->history_update) {
contact->history[history_index] += contact->magtwist * contact->dt;
}
// M_t torque (eq 30)
contact->magtortwist = -k * contact->history[history_index] - damp * contact->magtwist;
signtwist = (contact->magtwist > 0) - (contact->magtwist < 0);
Mtcrit = mu * contact->normal_model->Fncrit; // critical torque (eq 44)
if (fabs(contact->magtortwist) > Mtcrit) {
contact->history[history_index] = (Mtcrit * signtwist - damp * contact->magtwist) / k;
contact->magtortwist = -Mtcrit * signtwist; // eq 34
}
}

64
src/contact_twisting_models.h
View File

@ -1,64 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef CONTACT_TWISTING_MODELS_H_
#define CONTACT_TWISTING_MODELS_H_
#include "contact_sub_models.h"
namespace LAMMPS_NS {
namespace Contact {
class TwistingModel : public SubModel {
public:
TwistingModel(class LAMMPS *);
virtual ~TwistingModel() {};
virtual void coeffs_to_local() {};
virtual void init() {};
virtual void calculate_forces() = 0;
};
/* ---------------------------------------------------------------------- */
class TwistingNone : public TwistingModel {
public:
TwistingNone(class LAMMPS *);
void calculate_forces() {};
};
/* ---------------------------------------------------------------------- */
class TwistingMarshall : public TwistingModel {
public:
TwistingMarshall(class LAMMPS *);
void calculate_forces();
void init();
protected:
double k_tang, mu_tang;
};
/* ---------------------------------------------------------------------- */
class TwistingSDS : public TwistingModel {
public:
TwistingSDS(class LAMMPS *);
void coeffs_to_local() override;
void calculate_forces();
protected:
double k, mu, damp;
};
} // namespace Contact
} // namespace LAMMPS_NS
#endif /*CONTACT_TWISTING_MODELS_H_ */

6
src/dump_custom.cpp
View File

@ -41,7 +41,7 @@ enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS,
XSU,YSU,ZSU,XSUTRI,YSUTRI,ZSUTRI,
IX,IY,IZ,
VX,VY,VZ,FX,FY,FZ,
Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,heatflow,TEMPERATURE,
Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,HEATFLOW,TEMPERATURE,
OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
TQX,TQY,TQZ,
COMPUTE,FIX,VARIABLE,IVEC,DVEC,IARRAY,DARRAY};
@ -930,7 +930,7 @@ int DumpCustom::count()
for (i = 0; i < nlocal; i++) dchoose[i] = 2.0*radius[i];
ptr = dchoose;
nstride = 1;
} else if (thresh_array[ithresh] == heatflow) {
} else if (thresh_array[ithresh] == HEATFLOW) {
if (!atom->heatflow_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
@ -1867,7 +1867,7 @@ int DumpCustom::modify_param(int narg, char **arg)
else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS;
else if (strcmp(arg[1],"diameter") == 0) thresh_array[nthresh] = DIAMETER;
else if (strcmp(arg[1],"heatflow") == 0) thresh_array[nthresh] = heatflow;
else if (strcmp(arg[1],"heatflow") == 0) thresh_array[nthresh] = HEATFLOW;
else if (strcmp(arg[1],"temperature") == 0) thresh_array[nthresh] = TEMPERATURE;
else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX;
else if (strcmp(arg[1],"omegay") == 0) thresh_array[nthresh] = OMEGAY;