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

minor changes to source and doc files

Browse Source
This commit is contained in:
Steve Plimpton
2024-02-26 15:59:03 -07:00
parent 2a95c09d77
commit 4e77556610
3 changed files with 261 additions and 205 deletions
Download Patch File Download Diff File Expand all files Collapse all files

371
src/fix_indent.cpp
View File

@ -38,21 +38,6 @@ using namespace FixConst;
enum{NONE, SPHERE, CYLINDER, PLANE, CONE};
enum{INSIDE, OUTSIDE};
static bool PointInsideCone(int dir, double *center, double lo,
double hi, double rlo, double rhi, double *point);
static void DistanceExteriorPoint(int dir, double *center,
double lo, double hi, double rlo, double rhi, double &x,
double &y, double &z);
static void DistanceInteriorPoint(int dir, double *center,
double lo, double hi, double rlo, double rhi, double &x,
double &y, double &z);
static void point_on_line_segment(double *a, double *b, double *c, double *d);
static double closest(double *x, double *near, double *nearest, double dsq);
/* ---------------------------------------------------------------------- */
FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
@ -76,10 +61,11 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
if (k < 0.0) error->all(FLERR, "Illegal fix indent force constant: {}", k);
k3 = k/3.0;
// read options from end of input line
options(narg-4,&arg[4]);
// read geometry of indenter and optional args
int iarg = geometry(narg-4,&arg[4]) + 4;
options(narg-iarg,&arg[iarg]);
// setup scaling
const double xscale { scaleflag ? domain->lattice->xlattice : 1.0};
@ -93,8 +79,8 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
if (!ystr) yvalue *= yscale;
if (!zstr) zvalue *= zscale;
if (!rstr) rvalue *= xscale;
} else if (istyle == CONE) {
if (!xstr) xvalue *= xscale;
if (!ystr) yvalue *= yscale;
if (!zstr) zvalue *= zscale;
@ -121,6 +107,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
if (cdim == 0 && !pstr) pvalue *= xscale;
else if (cdim == 1 && !pstr) pvalue *= yscale;
else if (cdim == 2 && !pstr) pvalue *= zscale;
} else error->all(FLERR,"Unknown fix indent keyword: {}", istyle);
varflag = 0;
@ -195,7 +182,6 @@ void FixIndent::init()
if (!input->variable->equalstyle(pvar))
error->all(FLERR,"Variable {} for fix indent is invalid style", pstr);
}
if (rlostr) {
rlovar = input->variable->find(rlostr);
if (rlovar < 0)
@ -203,7 +189,6 @@ void FixIndent::init()
if (!input->variable->equalstyle(rlovar))
error->all(FLERR,"Variable {} for fix indent is invalid style", rlostr);
}
if (rhistr) {
rhivar = input->variable->find(rhistr);
if (rhivar < 0)
@ -211,7 +196,6 @@ void FixIndent::init()
if (!input->variable->equalstyle(rhivar))
error->all(FLERR,"Variable {} for fix indent is invalid style", rhistr);
}
if (lostr) {
lovar = input->variable->find(lostr);
if (lovar < 0)
@ -219,7 +203,6 @@ void FixIndent::init()
if (!input->variable->equalstyle(lovar))
error->all(FLERR,"Variable {} for fix indent is invalid style", lostr);
}
if (histr) {
hivar = input->variable->find(histr);
if (hivar < 0)
@ -267,6 +250,7 @@ void FixIndent::post_force(int /*vflag*/)
indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0;
// ctr = current indenter centerz
double ctr[3] {xvalue, yvalue, zvalue};
if (xstr) ctr[0] = input->variable->compute_equal(xvar);
if (ystr) ctr[1] = input->variable->compute_equal(yvar);
@ -284,6 +268,7 @@ void FixIndent::post_force(int /*vflag*/)
if (istyle == SPHERE) {
// remap indenter center into periodic box
domain->remap(ctr);
double radius { rstr ? input->variable->compute_equal(rvar) : rvalue};
@ -387,20 +372,24 @@ void FixIndent::post_force(int /*vflag*/)
double x0[3] {delx + ctr[0], dely + ctr[1], delz + ctr[2]};
r = sqrt(delx * delx + dely * dely + delz * delz);
// find if the particle is inside or outside the cone
// check if particle is inside or outside the cone
bool point_inside_cone = PointInsideCone(cdim, ctr, lo, hi, radiuslo, radiushi, x0);
if (side == INSIDE && point_inside_cone) continue;
if (side == OUTSIDE && !point_inside_cone) continue;
// find the distance between the point and the cone
if (point_inside_cone) {
DistanceInteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]);
} else {
DistanceExteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]);
}
// compute the force from the center of the cone - it is different from the approach of fix wall/region
// compute the force from the center of the cone
// this is different from how it is done in fix wall/region
dr = sqrt(x0[0] * x0[0] + x0[1] * x0[1] + x0[2] * x0[2]);
int force_sign = { point_inside_cone ? 1 : -1 };
@ -486,10 +475,10 @@ double FixIndent::compute_vector(int n)
}
/* ----------------------------------------------------------------------
parse optional parameters at end of input line
parse input args for geometry of indenter
------------------------------------------------------------------------- */
void FixIndent::options(int narg, char **arg)
int FixIndent::geometry(int narg, char **arg)
{
if (narg < 0) utils::missing_cmd_args(FLERR, "fix indent", error);
@ -499,139 +488,168 @@ void FixIndent::options(int narg, char **arg)
scaleflag = 1;
side = OUTSIDE;
// sphere
if (strcmp(arg[0],"sphere") == 0) {
if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword");
if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent sphere", error);
if (utils::strmatch(arg[1],"^v_")) {
xstr = utils::strdup(arg[1]+2);
} else xvalue = utils::numeric(FLERR,arg[1],false,lmp);
if (utils::strmatch(arg[2],"^v_")) {
ystr = utils::strdup(arg[2]+2);
} else yvalue = utils::numeric(FLERR,arg[2],false,lmp);
if (utils::strmatch(arg[3],"^v_")) {
zstr = utils::strdup(arg[3]+2);
} else zvalue = utils::numeric(FLERR,arg[3],false,lmp);
if (utils::strmatch(arg[4],"^v_")) {
rstr = utils::strdup(arg[4]+2);
} else rvalue = utils::numeric(FLERR,arg[4],false,lmp);
istyle = SPHERE;
return 5;
}
// cylinder
if (strcmp(arg[0],"cylinder") == 0) {
if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword");
if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error);
if (strcmp(arg[1],"x") == 0) {
cdim = 0;
if (utils::strmatch(arg[2],"^v_")) {
ystr = utils::strdup(arg[2]+2);
} else yvalue = utils::numeric(FLERR,arg[2],false,lmp);
if (utils::strmatch(arg[3],"^v_")) {
zstr = utils::strdup(arg[3]+2);
} else zvalue = utils::numeric(FLERR,arg[3],false,lmp);
} else if (strcmp(arg[1],"y") == 0) {
cdim = 1;
if (utils::strmatch(arg[2],"^v_")) {
xstr = utils::strdup(arg[2]+2);
} else xvalue = utils::numeric(FLERR,arg[2],false,lmp);
if (utils::strmatch(arg[3],"^v_")) {
zstr = utils::strdup(arg[3]+2);
} else zvalue = utils::numeric(FLERR,arg[3],false,lmp);
} else if (strcmp(arg[1],"z") == 0) {
cdim = 2;
if (utils::strmatch(arg[2],"^v_")) {
xstr = utils::strdup(arg[2]+2);
} else xvalue = utils::numeric(FLERR,arg[2],false,lmp);
if (utils::strmatch(arg[3],"^v_")) {
ystr = utils::strdup(arg[3]+2);
} else yvalue = utils::numeric(FLERR,arg[3],false,lmp);
} else error->all(FLERR,"Unknown fix indent cylinder argument: {}", arg[1]);
if (utils::strmatch(arg[4],"^v_")) {
rstr = utils::strdup(arg[4]+2);
} else rvalue = utils::numeric(FLERR,arg[4],false,lmp);
istyle = CYLINDER;
return 5;
}
// cone
if (strcmp(arg[0],"cone") == 0) {
if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword");
if (8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error);
if (strcmp(arg[1],"x") == 0) {
cdim = 0;
if (utils::strmatch(arg[2],"^v_")) {
ystr = utils::strdup(arg[2]+2);
} else yvalue = utils::numeric(FLERR,arg[2],false,lmp);
if (utils::strmatch(arg[3],"^v_")) {
zstr = utils::strdup(arg[3]+2);
} else zvalue = utils::numeric(FLERR,arg[3],false,lmp);
} else if (strcmp(arg[1],"y") == 0) {
cdim = 1;
if (utils::strmatch(arg[2],"^v_")) {
xstr = utils::strdup(arg[2]+2);
} else xvalue = utils::numeric(FLERR,arg[2],false,lmp);
if (utils::strmatch(arg[3],"^v_")) {
zstr = utils::strdup(arg[3]+2);
} else zvalue = utils::numeric(FLERR,arg[3],false,lmp);
} else if (strcmp(arg[1],"z") == 0) {
cdim = 2;
if (utils::strmatch(arg[2],"^v_")) {
xstr = utils::strdup(arg[2]+2);
} else xvalue = utils::numeric(FLERR,arg[2],false,lmp);
if (utils::strmatch(arg[3],"^v_")) {
ystr = utils::strdup(arg[3]+2);
} else yvalue = utils::numeric(FLERR,arg[3],false,lmp);
} else error->all(FLERR,"Unknown fix indent cone argument: {}", arg[1]);
if (utils::strmatch(arg[4],"^v_")) {
rlostr = utils::strdup(arg[4]+2);
} else rlovalue = utils::numeric(FLERR,arg[4],false,lmp);
if (utils::strmatch(arg[5],"^v_")) {
rhistr = utils::strdup(arg[5]+2);
} else rhivalue = utils::numeric(FLERR,arg[5],false,lmp);
if (utils::strmatch(arg[6],"^v_")) {
lostr = utils::strdup(arg[6]+2);
} else lovalue = utils::numeric(FLERR,arg[6],false,lmp);
if (utils::strmatch(arg[7],"^v_")) {
histr = utils::strdup(arg[7]+2);
} else hivalue = utils::numeric(FLERR,arg[7],false,lmp);
istyle = CONE;
return 8;
}
// plane
if (strcmp(arg[0],"plane") == 0) {
if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword");
if (4 > narg) utils::missing_cmd_args(FLERR, "fix indent plane", error);
if (strcmp(arg[1],"x") == 0) cdim = 0;
else if (strcmp(arg[1],"y") == 0) cdim = 1;
else if (strcmp(arg[1],"z") == 0) cdim = 2;
else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[1]);
if (utils::strmatch(arg[2],"^v_")) {
pstr = utils::strdup(arg[2]+2);
} else pvalue = utils::numeric(FLERR,arg[2],false,lmp);
if (strcmp(arg[3],"lo") == 0) planeside = -1;
else if (strcmp(arg[3],"hi") == 0) planeside = 1;
else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[3]);
istyle = PLANE;
return 4;
}
// invalid istyle arg
error->all(FLERR,"Unknown fix indent argument: {}", arg[0]);
return 0;
}
/* ----------------------------------------------------------------------
parse optional input args
------------------------------------------------------------------------- */
void FixIndent::options(int narg, char **arg)
{
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"sphere") == 0) {
if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "fix indent sphere", error);
if (utils::strmatch(arg[iarg+1],"^v_")) {
xstr = utils::strdup(arg[iarg+1]+2);
} else xvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp);
if (utils::strmatch(arg[iarg+2],"^v_")) {
ystr = utils::strdup(arg[iarg+2]+2);
} else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
zstr = utils::strdup(arg[iarg+3]+2);
} else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (utils::strmatch(arg[iarg+4],"^v_")) {
rstr = utils::strdup(arg[iarg+4]+2);
} else rvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp);
istyle = SPHERE;
iarg += 5;
} else if (strcmp(arg[iarg],"cylinder") == 0) {
if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error);
if (strcmp(arg[iarg+1],"x") == 0) {
cdim = 0;
if (utils::strmatch(arg[iarg+2],"^v_")) {
ystr = utils::strdup(arg[iarg+2]+2);
} else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
zstr = utils::strdup(arg[iarg+3]+2);
} else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
} else if (strcmp(arg[iarg+1],"y") == 0) {
cdim = 1;
if (utils::strmatch(arg[iarg+2],"^v_")) {
xstr = utils::strdup(arg[iarg+2]+2);
} else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
zstr = utils::strdup(arg[iarg+3]+2);
} else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
} else if (strcmp(arg[iarg+1],"z") == 0) {
cdim = 2;
if (utils::strmatch(arg[iarg+2],"^v_")) {
xstr = utils::strdup(arg[iarg+2]+2);
} else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
ystr = utils::strdup(arg[iarg+3]+2);
} else yvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
} else error->all(FLERR,"Unknown fix indent cylinder argument: {}", arg[iarg+1]);
if (utils::strmatch(arg[iarg+4],"^v_")) {
rstr = utils::strdup(arg[iarg+4]+2);
} else rvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp);
istyle = CYLINDER;
iarg += 5;
} else if (strcmp(arg[iarg],"plane") == 0) {
if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "fix indent plane", error);
if (strcmp(arg[iarg+1],"x") == 0) cdim = 0;
else if (strcmp(arg[iarg+1],"y") == 0) cdim = 1;
else if (strcmp(arg[iarg+1],"z") == 0) cdim = 2;
else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[iarg+1]);
if (utils::strmatch(arg[iarg+2],"^v_")) {
pstr = utils::strdup(arg[iarg+2]+2);
} else pvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (strcmp(arg[iarg+3],"lo") == 0) planeside = -1;
else if (strcmp(arg[iarg+3],"hi") == 0) planeside = 1;
else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[iarg+3]);
istyle = PLANE;
iarg += 4;
} else if (strcmp(arg[iarg],"cone") == 0) {
if (iarg+8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error);
if (strcmp(arg[iarg+1],"x") == 0) {
cdim = 0;
if (utils::strmatch(arg[iarg+2],"^v_")) {
ystr = utils::strdup(arg[iarg+2]+2);
} else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
zstr = utils::strdup(arg[iarg+3]+2);
} else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
} else if (strcmp(arg[iarg+1],"y") == 0) {
cdim = 1;
if (utils::strmatch(arg[iarg+2],"^v_")) {
xstr = utils::strdup(arg[iarg+2]+2);
} else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
zstr = utils::strdup(arg[iarg+3]+2);
} else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
} else if (strcmp(arg[iarg+1],"z") == 0) {
cdim = 2;
if (utils::strmatch(arg[iarg+2],"^v_")) {
xstr = utils::strdup(arg[iarg+2]+2);
} else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (utils::strmatch(arg[iarg+3],"^v_")) {
ystr = utils::strdup(arg[iarg+3]+2);
} else yvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp);
} else error->all(FLERR,"Unknown fix indent cone argument: {}", arg[iarg+1]);
if (utils::strmatch(arg[iarg+4],"^v_")) {
rlostr = utils::strdup(arg[iarg+4]+2);
} else rlovalue = utils::numeric(FLERR,arg[iarg+4],false,lmp);
if (utils::strmatch(arg[iarg+5],"^v_")) {
rhistr = utils::strdup(arg[iarg+5]+2);
} else rhivalue = utils::numeric(FLERR,arg[iarg+5],false,lmp);
if (utils::strmatch(arg[iarg+6],"^v_")) {
lostr = utils::strdup(arg[iarg+6]+2);
} else lovalue = utils::numeric(FLERR,arg[iarg+6],false,lmp);
if (utils::strmatch(arg[iarg+7],"^v_")) {
histr = utils::strdup(arg[iarg+7]+2);
} else hivalue = utils::numeric(FLERR,arg[iarg+7],false,lmp);
istyle = CONE;
iarg += 8;
} else if (strcmp(arg[iarg],"units") == 0) {
if (strcmp(arg[iarg],"units") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix indent units", error);
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
@ -644,16 +662,17 @@ void FixIndent::options(int narg, char **arg)
else if (strcmp(arg[iarg+1],"out") == 0) side = OUTSIDE;
else error->all(FLERR,"Unknown fix indent side argument: {}", arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Unknown fix indent argument: {}", arg[iarg]);
}
}
/* ----------------------------------------------------------------------
determines if a point is inside (true) or outside (false) of a cone
determines if a point is inside (true) or outside (false) of a cone
------------------------------------------------------------------------- */
bool PointInsideCone(int dir, double *center, double lo,
double hi, double rlo, double rhi, double *x)
bool FixIndent::PointInsideCone(int dir, double *center, double lo,
double hi, double rlo, double rhi, double *x)
{
if ((x[dir] > hi) || (x[dir] < lo)) return false;
@ -669,10 +688,12 @@ bool PointInsideCone(int dir, double *center, double lo,
}
/* ----------------------------------------------------------------------
distance between an exterior point and a cone
distance between an exterior point and a cone
------------------------------------------------------------------------- */
void DistanceExteriorPoint(int dir, double *center, double lo, double hi,
double rlo, double rhi, double &x, double &y, double &z)
void FixIndent::DistanceExteriorPoint(int dir, double *center, double lo, double hi,
double rlo, double rhi,
double &x, double &y, double &z)
{
double xp[3], nearest[3], corner1[3], corner2[3];
@ -700,17 +721,21 @@ void DistanceExteriorPoint(int dir, double *center, double lo, double hi,
corner4[dir] = hi;
// initialize distance to a big number
double distsq = 1.0e20;
// check the first triangle
point_on_line_segment(corner1, corner2, point, xp);
distsq = closest(point, xp, nearest, distsq);
// check the second triangle
point_on_line_segment(corner1, corner3, point, xp);
distsq = closest(point, xp, nearest, distsq);
// check the third triangle
point_on_line_segment(corner2, corner4, point, xp);
distsq = closest(point, xp, nearest, distsq);
@ -722,12 +747,12 @@ void DistanceExteriorPoint(int dir, double *center, double lo, double hi,
}
/* ----------------------------------------------------------------------
distance between an interior point and a cone
distance between an interior point and a cone
------------------------------------------------------------------------- */
void DistanceInteriorPoint(int dir, double *center,
double lo, double hi, double rlo, double rhi, double &x,
double &y, double &z)
void FixIndent::DistanceInteriorPoint(int dir, double *center,
double lo, double hi, double rlo, double rhi, double &x,
double &y, double &z)
{
double r, dist_disk, dist_surf;
double surflo[3], surfhi[3], xs[3];
@ -735,6 +760,7 @@ void DistanceInteriorPoint(int dir, double *center,
double point[3] {0.0, 0.0, 0.0};
// initial check with the two disks
if ( (initial_point[dir] - lo) < (hi - initial_point[dir]) ) {
dist_disk = (initial_point[dir] - lo) * (initial_point[dir] - lo);
point[dir] = initial_point[dir] - lo;
@ -744,6 +770,7 @@ void DistanceInteriorPoint(int dir, double *center,
}
// check with the points in the conical surface
double del[3] {x - center[0], y - center[1], z - center[2]};
del[dir] = 0.0;
r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]);
@ -776,10 +803,10 @@ void DistanceInteriorPoint(int dir, double *center,
}
/* ----------------------------------------------------------------------
helper function extracted from region.cpp
helper function extracted from region.cpp
------------------------------------------------------------------------- */
void point_on_line_segment(double *a, double *b, double *c, double *d)
void FixIndent::point_on_line_segment(double *a, double *b, double *c, double *d)
{
double ba[3], ca[3];
@ -802,10 +829,10 @@ void point_on_line_segment(double *a, double *b, double *c, double *d)
}
/* ----------------------------------------------------------------------
helper function extracted from region_cone.cpp
helper function extracted from region_cone.cpp
------------------------------------------------------------------------- */
double closest(double *x, double *near, double *nearest, double dsq)
double FixIndent::closest(double *x, double *near, double *nearest, double dsq)
{
double dx = x[0] - near[0];
double dy = x[1] - near[1];