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

enable and apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-04-03 21:41:09 -04:00
parent b53a47b192
commit 27127a46cc
5 changed files with 844 additions and 734 deletions
Download Patch File Download Diff File Expand all files Collapse all files

357
src/region_block.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -24,100 +23,125 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{CONSTANT,VARIABLE}; enum { CONSTANT, VARIABLE };
#define BIG 1.0e20 static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), xlostr(nullptr), xhistr(nullptr), ylostr(nullptr), yhistr(nullptr), zlostr(nullptr), zhistr(nullptr) Region(lmp, narg, arg), xlostr(nullptr), xhistr(nullptr), ylostr(nullptr), yhistr(nullptr),
zlostr(nullptr), zhistr(nullptr)
{ {
options(narg-8,&arg[8]); options(narg - 8, &arg[8]);
xlostyle = CONSTANT; xlostyle = CONSTANT;
if (strcmp(arg[2],"INF") == 0 || strcmp(arg[2],"EDGE") == 0) { if (strcmp(arg[2], "INF") == 0 || strcmp(arg[2], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[2],"INF") == 0) xlo = -BIG; if (strcmp(arg[2], "INF") == 0)
else if (domain->triclinic == 0) xlo = domain->boxlo[0]; xlo = -BIG;
else xlo = domain->boxlo_bound[0]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[2],"^v_")) { xlo = domain->boxlo[0];
xlostr = utils::strdup(arg[2]+2); else
xlo = 0.0; xlo = domain->boxlo_bound[0];
xlostyle = VARIABLE; } else if (utils::strmatch(arg[2], "^v_")) {
varshape = 1; xlostr = utils::strdup(arg[2] + 2);
} else xlo = xscale*utils::numeric(FLERR,arg[2],false,lmp); xlo = 0.0;
xlostyle = VARIABLE;
varshape = 1;
} else
xlo = xscale * utils::numeric(FLERR, arg[2], false, lmp);
xhistyle = CONSTANT; xhistyle = CONSTANT;
if (strcmp(arg[3],"INF") == 0 || strcmp(arg[3],"EDGE") == 0) { if (strcmp(arg[3], "INF") == 0 || strcmp(arg[3], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[3],"INF") == 0) xhi = BIG; if (strcmp(arg[3], "INF") == 0)
else if (domain->triclinic == 0) xhi = domain->boxhi[0]; xhi = BIG;
else xhi = domain->boxhi_bound[0]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[3],"^v_")) { xhi = domain->boxhi[0];
xhistr = utils::strdup(arg[3]+2); else
xhi = 0.0; xhi = domain->boxhi_bound[0];
xhistyle = VARIABLE; } else if (utils::strmatch(arg[3], "^v_")) {
varshape = 1; xhistr = utils::strdup(arg[3] + 2);
} else xhi = xscale*utils::numeric(FLERR,arg[3],false,lmp); xhi = 0.0;
xhistyle = VARIABLE;
varshape = 1;
} else
xhi = xscale * utils::numeric(FLERR, arg[3], false, lmp);
ylostyle = CONSTANT; ylostyle = CONSTANT;
if (strcmp(arg[4],"INF") == 0 || strcmp(arg[4],"EDGE") == 0) { if (strcmp(arg[4], "INF") == 0 || strcmp(arg[4], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[4],"INF") == 0) ylo = -BIG; if (strcmp(arg[4], "INF") == 0)
else if (domain->triclinic == 0) ylo = domain->boxlo[1]; ylo = -BIG;
else ylo = domain->boxlo_bound[1]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[4],"^v_")) { ylo = domain->boxlo[1];
ylostr = utils::strdup(arg[4]+2); else
ylo = 0.0; ylo = domain->boxlo_bound[1];
ylostyle = VARIABLE; } else if (utils::strmatch(arg[4], "^v_")) {
varshape = 1; ylostr = utils::strdup(arg[4] + 2);
} else ylo = yscale*utils::numeric(FLERR,arg[4],false,lmp); ylo = 0.0;
ylostyle = VARIABLE;
varshape = 1;
} else
ylo = yscale * utils::numeric(FLERR, arg[4], false, lmp);
yhistyle = CONSTANT; yhistyle = CONSTANT;
if (strcmp(arg[5],"INF") == 0 || strcmp(arg[5],"EDGE") == 0) { if (strcmp(arg[5], "INF") == 0 || strcmp(arg[5], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[5],"INF") == 0) yhi = BIG; if (strcmp(arg[5], "INF") == 0)
else if (domain->triclinic == 0) yhi = domain->boxhi[1]; yhi = BIG;
else yhi = domain->boxhi_bound[1]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[5],"^v_")) { yhi = domain->boxhi[1];
yhistr = utils::strdup(arg[5]+2); else
yhi = 0.0; yhi = domain->boxhi_bound[1];
yhistyle = VARIABLE; } else if (utils::strmatch(arg[5], "^v_")) {
varshape = 1; yhistr = utils::strdup(arg[5] + 2);
} else yhi = yscale*utils::numeric(FLERR,arg[5],false,lmp); yhi = 0.0;
yhistyle = VARIABLE;
varshape = 1;
} else
yhi = yscale * utils::numeric(FLERR, arg[5], false, lmp);
zlostyle = CONSTANT; zlostyle = CONSTANT;
if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { if (strcmp(arg[6], "INF") == 0 || strcmp(arg[6], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[6],"INF") == 0) zlo = -BIG; if (strcmp(arg[6], "INF") == 0)
else if (domain->triclinic == 0) zlo = domain->boxlo[2]; zlo = -BIG;
else zlo = domain->boxlo_bound[2]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[6],"^v_")) { zlo = domain->boxlo[2];
zlostr = utils::strdup(arg[6]+2); else
zlo = 0.0; zlo = domain->boxlo_bound[2];
zlostyle = VARIABLE; } else if (utils::strmatch(arg[6], "^v_")) {
varshape = 1; zlostr = utils::strdup(arg[6] + 2);
} else zlo = zscale*utils::numeric(FLERR,arg[6],false,lmp); zlo = 0.0;
zlostyle = VARIABLE;
varshape = 1;
} else
zlo = zscale * utils::numeric(FLERR, arg[6], false, lmp);
zhistyle = CONSTANT; zhistyle = CONSTANT;
if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[7],"INF") == 0) zhi = BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) zhi = domain->boxhi[2]; zhi = BIG;
else zhi = domain->boxhi_bound[2]; else if (domain->triclinic == 0)
} else if (utils::strmatch(arg[7],"^v_")) { zhi = domain->boxhi[2];
zhistr = utils::strdup(arg[7]+2); else
zhi = 0.0; zhi = domain->boxhi_bound[2];
zhistyle = VARIABLE; } else if (utils::strmatch(arg[7], "^v_")) {
varshape = 1; zhistr = utils::strdup(arg[7] + 2);
} else zhi = zscale*utils::numeric(FLERR,arg[7],false,lmp); zhi = 0.0;
zhistyle = VARIABLE;
varshape = 1;
} else
zhi = zscale * utils::numeric(FLERR, arg[7], false, lmp);
if (varshape) { if (varshape) {
variable_check(); variable_check();
@ -126,9 +150,9 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
// error check // error check
if (xlo > xhi) error->all(FLERR,"Illegal region block xlo: {} >= xhi: {}", xlo, xhi); if (xlo > xhi) error->all(FLERR, "Illegal region block xlo: {} >= xhi: {}", xlo, xhi);
if (ylo > yhi) error->all(FLERR,"Illegal region block ylo: {} >= yhi: {}", ylo, yhi); if (ylo > yhi) error->all(FLERR, "Illegal region block ylo: {} >= yhi: {}", ylo, yhi);
if (zlo > zhi) error->all(FLERR,"Illegal region block zlo: {} >= zhi: {}", zlo, zhi); if (zlo > zhi) error->all(FLERR, "Illegal region block zlo: {} >= zhi: {}", zlo, zhi);
// extent of block // extent of block
@ -140,15 +164,18 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
extent_yhi = yhi; extent_yhi = yhi;
extent_zlo = zlo; extent_zlo = zlo;
extent_zhi = zhi; extent_zhi = zhi;
} else bboxflag = 0; } else
bboxflag = 0;
// particle could be close to all 6 planes // particle could be close to all 6 planes
// particle can only touch 3 planes // particle can only touch 3 planes
cmax = 6; cmax = 6;
contact = new Contact[cmax]; contact = new Contact[cmax];
if (interior) tmax = 3; if (interior)
else tmax = 1; tmax = 3;
else
tmax = 1;
// open face data structs // open face data structs
@ -235,13 +262,13 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
RegBlock::~RegBlock() RegBlock::~RegBlock()
{ {
if (copymode) return; if (copymode) return;
delete [] xlostr; delete[] xlostr;
delete [] xhistr; delete[] xhistr;
delete [] ylostr; delete[] ylostr;
delete [] yhistr; delete[] yhistr;
delete [] zlostr; delete[] zlostr;
delete [] zhistr; delete[] zhistr;
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -259,8 +286,7 @@ void RegBlock::init()
int RegBlock::inside(double x, double y, double z) int RegBlock::inside(double x, double y, double z)
{ {
if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) return 1;
return 1;
return 0; return 0;
} }
@ -277,8 +303,7 @@ int RegBlock::surface_interior(double *x, double cutoff)
// x is exterior to block // x is exterior to block
if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || x[2] < zlo || x[2] > zhi) return 0;
x[2] < zlo || x[2] > zhi) return 0;
// x is interior to block or on its surface // x is interior to block or on its surface
@ -352,17 +377,16 @@ int RegBlock::surface_interior(double *x, double cutoff)
int RegBlock::surface_exterior(double *x, double cutoff) int RegBlock::surface_exterior(double *x, double cutoff)
{ {
double xp,yp,zp; double xp, yp, zp;
double xc,yc,zc,dist,mindist; double xc, yc, zc, dist, mindist;
// x is far enough from block that there is no contact // x is far enough from block that there is no contact
// x is interior to block // x is interior to block
if (x[0] <= xlo-cutoff || x[0] >= xhi+cutoff || if (x[0] <= xlo - cutoff || x[0] >= xhi + cutoff || x[1] <= ylo - cutoff ||
x[1] <= ylo-cutoff || x[1] >= yhi+cutoff || x[1] >= yhi + cutoff || x[2] <= zlo - cutoff || x[2] >= zhi + cutoff)
x[2] <= zlo-cutoff || x[2] >= zhi+cutoff) return 0; return 0;
if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && x[2] > zlo && x[2] < zhi) return 0;
x[2] > zlo && x[2] < zhi) return 0;
// x is exterior to block or on its surface // x is exterior to block or on its surface
// xp,yp,zp = point on surface of block that x is closest to // xp,yp,zp = point on surface of block that x is closest to
@ -370,20 +394,29 @@ int RegBlock::surface_exterior(double *x, double cutoff)
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
if (!openflag) { if (!openflag) {
if (x[0] < xlo) xp = xlo; if (x[0] < xlo)
else if (x[0] > xhi) xp = xhi; xp = xlo;
else xp = x[0]; else if (x[0] > xhi)
if (x[1] < ylo) yp = ylo; xp = xhi;
else if (x[1] > yhi) yp = yhi; else
else yp = x[1]; xp = x[0];
if (x[2] < zlo) zp = zlo; if (x[1] < ylo)
else if (x[2] > zhi) zp = zhi; yp = ylo;
else zp = x[2]; else if (x[1] > yhi)
yp = yhi;
else
yp = x[1];
if (x[2] < zlo)
zp = zlo;
else if (x[2] > zhi)
zp = zhi;
else
zp = x[2];
} else { } else {
mindist = BIG; mindist = BIG;
for (int i = 0; i < 6; i++) { for (int i = 0; i < 6; i++) {
if (open_faces[i]) continue; if (open_faces[i]) continue;
dist = find_closest_point(i,x,xc,yc,zc); dist = find_closest_point(i, x, xc, yc, zc);
if (dist < mindist) { if (dist < mindist) {
xp = xc; xp = xc;
yp = yc; yp = yc;
@ -393,7 +426,7 @@ int RegBlock::surface_exterior(double *x, double cutoff)
} }
} }
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
return 0; return 0;
@ -403,28 +436,17 @@ int RegBlock::surface_exterior(double *x, double cutoff)
change region shape via variable evaluation change region shape via variable evaluation
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void RegBlock::shape_update() // addition void RegBlock::shape_update() // addition
{ {
if (xlostyle == VARIABLE) if (xlostyle == VARIABLE) xlo = xscale * input->variable->compute_equal(xlovar);
xlo = xscale * input->variable->compute_equal(xlovar); if (xhistyle == VARIABLE) xhi = xscale * input->variable->compute_equal(xhivar);
if (ylostyle == VARIABLE) ylo = yscale * input->variable->compute_equal(ylovar);
if (xhistyle == VARIABLE) if (yhistyle == VARIABLE) yhi = yscale * input->variable->compute_equal(yhivar);
xhi = xscale * input->variable->compute_equal(xhivar); if (zlostyle == VARIABLE) zlo = zscale * input->variable->compute_equal(zlovar);
if (zhistyle == VARIABLE) zhi = zscale * input->variable->compute_equal(zhivar);
if (ylostyle == VARIABLE)
ylo = yscale * input->variable->compute_equal(ylovar);
if (yhistyle == VARIABLE)
yhi = yscale * input->variable->compute_equal(yhivar);
if (zlostyle == VARIABLE)
zlo = zscale * input->variable->compute_equal(zlovar);
if (zhistyle == VARIABLE)
zhi = zscale * input->variable->compute_equal(zhivar);
if (xlo > xhi || ylo > yhi || zlo > zhi) if (xlo > xhi || ylo > yhi || zlo > zhi)
error->one(FLERR,"Variable evaluation in region gave bad value"); error->one(FLERR, "Variable evaluation in region gave bad value");
// face[0] // face[0]
@ -489,54 +511,48 @@ void RegBlock::shape_update() // addition
error check on existence of variable error check on existence of variable
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void RegBlock::variable_check() // addition void RegBlock::variable_check() // addition
{ {
if (xlostyle == VARIABLE) { if (xlostyle == VARIABLE) {
xlovar = input->variable->find(xlostr); xlovar = input->variable->find(xlostr);
if (xlovar < 0) if (xlovar < 0) error->all(FLERR, "Variable name for region block does not exist");
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(xlovar)) if (!input->variable->equalstyle(xlovar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable for region block is invalid style");
} }
if (xhistyle == VARIABLE) { if (xhistyle == VARIABLE) {
xhivar = input->variable->find(xhistr); xhivar = input->variable->find(xhistr);
if (xhivar < 0) if (xhivar < 0) error->all(FLERR, "Variable name for region block does not exist");
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(xhivar)) if (!input->variable->equalstyle(xhivar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable for region block is invalid style");
} }
if (ylostyle == VARIABLE) { if (ylostyle == VARIABLE) {
ylovar = input->variable->find(ylostr); ylovar = input->variable->find(ylostr);
if (ylovar < 0) if (ylovar < 0) error->all(FLERR, "Variable name for region block does not exist");
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(ylovar)) if (!input->variable->equalstyle(ylovar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable for region block is invalid style");
} }
if (yhistyle == VARIABLE) { if (yhistyle == VARIABLE) {
yhivar = input->variable->find(yhistr); yhivar = input->variable->find(yhistr);
if (yhivar < 0) if (yhivar < 0) error->all(FLERR, "Variable name for region block does not exist");
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(yhivar)) if (!input->variable->equalstyle(yhivar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable for region block is invalid style");
} }
if (zlostyle == VARIABLE) { if (zlostyle == VARIABLE) {
zlovar = input->variable->find(zlostr); zlovar = input->variable->find(zlostr);
if (zlovar < 0) if (zlovar < 0) error->all(FLERR, "Variable name for region block does not exist");
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(zlovar)) if (!input->variable->equalstyle(zlovar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable for region block is invalid style");
} }
if (zhistyle == VARIABLE) { if (zhistyle == VARIABLE) {
zhivar = input->variable->find(zhistr); zhivar = input->variable->find(zhistr);
if (zhivar < 0) if (zhivar < 0) error->all(FLERR, "Variable name for region block does not exist");
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(zhivar)) if (!input->variable->equalstyle(zhivar))
error->all(FLERR,"Variable for region block is invalid style"); error->all(FLERR, "Variable for region block is invalid style");
} }
} }
@ -545,36 +561,35 @@ void RegBlock::variable_check() // addition
store closest point in xc,yc,zc store closest point in xc,yc,zc
--------------------------------------------------------------------------*/ --------------------------------------------------------------------------*/
double RegBlock::find_closest_point(int i, double *x, double RegBlock::find_closest_point(int i, double *x, double &xc, double &yc, double &zc)
double &xc, double &yc, double &zc)
{ {
double dot,d2,d2min; double dot, d2, d2min;
double xr[3],xproj[3],p[3]; double xr[3], xproj[3], p[3];
xr[0] = x[0] - corners[i][0][0]; xr[0] = x[0] - corners[i][0][0];
xr[1] = x[1] - corners[i][0][1]; xr[1] = x[1] - corners[i][0][1];
xr[2] = x[2] - corners[i][0][2]; xr[2] = x[2] - corners[i][0][2];
dot = face[i][0]*xr[0] + face[i][1]*xr[1] + face[i][2]*xr[2]; dot = face[i][0] * xr[0] + face[i][1] * xr[1] + face[i][2] * xr[2];
xproj[0] = xr[0] - dot*face[i][0]; xproj[0] = xr[0] - dot * face[i][0];
xproj[1] = xr[1] - dot*face[i][1]; xproj[1] = xr[1] - dot * face[i][1];
xproj[2] = xr[2] - dot*face[i][2]; xproj[2] = xr[2] - dot * face[i][2];
d2min = BIG; d2min = BIG;
// check if point projects inside of face // check if point projects inside of face
if (inside_face(xproj, i)) { if (inside_face(xproj, i)) {
d2 = d2min = dot*dot; d2 = d2min = dot * dot;
xc = xproj[0] + corners[i][0][0]; xc = xproj[0] + corners[i][0][0];
yc = xproj[1] + corners[i][0][1]; yc = xproj[1] + corners[i][0][1];
zc = xproj[2] + corners[i][0][2]; zc = xproj[2] + corners[i][0][2];
// check each edge // check each edge
} else { } else {
point_on_line_segment(corners[i][0],corners[i][1],x,p); point_on_line_segment(corners[i][0], corners[i][1], x, p);
d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2]-x[2])*(p[2]-x[2]); (p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) { if (d2 < d2min) {
d2min = d2; d2min = d2;
xc = p[0]; xc = p[0];
@ -582,9 +597,9 @@ double RegBlock::find_closest_point(int i, double *x,
zc = p[2]; zc = p[2];
} }
point_on_line_segment(corners[i][1],corners[i][2],x,p); point_on_line_segment(corners[i][1], corners[i][2], x, p);
d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2]-x[2])*(p[2]-x[2]); (p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) { if (d2 < d2min) {
d2min = d2; d2min = d2;
xc = p[0]; xc = p[0];
@ -592,9 +607,9 @@ double RegBlock::find_closest_point(int i, double *x,
zc = p[2]; zc = p[2];
} }
point_on_line_segment(corners[i][2],corners[i][3],x,p); point_on_line_segment(corners[i][2], corners[i][3], x, p);
d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2]-x[2])*(p[2]-x[2]); (p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) { if (d2 < d2min) {
d2min = d2; d2min = d2;
xc = p[0]; xc = p[0];
@ -602,9 +617,9 @@ double RegBlock::find_closest_point(int i, double *x,
zc = p[2]; zc = p[2];
} }
point_on_line_segment(corners[i][3],corners[i][0],x,p); point_on_line_segment(corners[i][3], corners[i][0], x, p);
d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2]-x[2])*(p[2]-x[2]); (p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) { if (d2 < d2min) {
d2min = d2; d2min = d2;
xc = p[0]; xc = p[0];
@ -623,14 +638,12 @@ double RegBlock::find_closest_point(int i, double *x,
int RegBlock::inside_face(double *xproj, int iface) int RegBlock::inside_face(double *xproj, int iface)
{ {
if (iface < 2) { if (iface < 2) {
if (xproj[1] > 0 && (xproj[1] < yhi-ylo) && if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1;
xproj[2] > 0 && (xproj[2] < zhi-zlo)) return 1;
} else if (iface < 4) { } else if (iface < 4) {
if (xproj[0] > 0 && (xproj[0] < (xhi-xlo)) && if (xproj[0] > 0 && (xproj[0] < (xhi - xlo)) && xproj[2] > 0 && (xproj[2] < (zhi - zlo)))
xproj[2] > 0 && (xproj[2] < (zhi-zlo))) return 1; return 1;
} else { } else {
if (xproj[0] > 0 && xproj[0] < (xhi-xlo) && if (xproj[0] > 0 && xproj[0] < (xhi - xlo) && xproj[1] > 0 && xproj[1] < (yhi - ylo)) return 1;
xproj[1] > 0 && xproj[1] < (yhi-ylo)) return 1;
} }
return 0; return 0;

335
src/region_cone.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -26,102 +25,120 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define BIG 1.0e20 static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo(0.0), hi(0.0)
Region(lmp, narg, arg), lo(0.0), hi(0.0)
{ {
options(narg-9,&arg[9]); options(narg - 9, &arg[9]);
// check open face settings // check open face settings
if (openflag) if (openflag)
for (int i=3; i<6; i++) for (int i = 3; i < 6; i++)
if (open_faces[i]) error->all(FLERR,"Illegal region cone open face: {}", i+1); if (open_faces[i]) error->all(FLERR, "Illegal region cone open face: {}", i + 1);
if (strcmp(arg[2],"x") != 0 && strcmp(arg[2],"y") != 0 && strcmp(arg[2],"z") != 0) if (strcmp(arg[2], "x") != 0 && strcmp(arg[2], "y") != 0 && strcmp(arg[2], "z") != 0)
error->all(FLERR,"Illegal region cone axis: {}", arg[2]); error->all(FLERR, "Illegal region cone axis: {}", arg[2]);
axis = arg[2][0]; axis = arg[2][0];
if (axis == 'x') { if (axis == 'x') {
c1 = yscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp);
c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
radiuslo = yscale*utils::numeric(FLERR,arg[5],false,lmp); radiuslo = yscale * utils::numeric(FLERR, arg[5], false, lmp);
radiushi = yscale*utils::numeric(FLERR,arg[6],false,lmp); radiushi = yscale * utils::numeric(FLERR, arg[6], false, lmp);
} else if (axis == 'y') { } else if (axis == 'y') {
c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
radiuslo = xscale*utils::numeric(FLERR,arg[5],false,lmp); radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp);
radiushi = xscale*utils::numeric(FLERR,arg[6],false,lmp); radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp);
} else if (axis == 'z') { } else if (axis == 'z') {
c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c2 = yscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp);
radiuslo = xscale*utils::numeric(FLERR,arg[5],false,lmp); radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp);
radiushi = xscale*utils::numeric(FLERR,arg[6],false,lmp); radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp);
} }
if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (axis == 'x') { if (axis == 'x') {
if (strcmp(arg[7],"INF") == 0) lo = -BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[0]; lo = -BIG;
else lo = domain->boxlo_bound[0]; else if (domain->triclinic == 0)
lo = domain->boxlo[0];
else
lo = domain->boxlo_bound[0];
} }
if (axis == 'y') { if (axis == 'y') {
if (strcmp(arg[7],"INF") == 0) lo = -BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[1]; lo = -BIG;
else lo = domain->boxlo_bound[1]; else if (domain->triclinic == 0)
lo = domain->boxlo[1];
else
lo = domain->boxlo_bound[1];
} }
if (axis == 'z') { if (axis == 'z') {
if (strcmp(arg[7],"INF") == 0) lo = -BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[2]; lo = -BIG;
else lo = domain->boxlo_bound[2]; else if (domain->triclinic == 0)
lo = domain->boxlo[2];
else
lo = domain->boxlo_bound[2];
} }
} else { } else {
if (axis == 'x') lo = xscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'x') lo = xscale * utils::numeric(FLERR, arg[7], false, lmp);
if (axis == 'y') lo = yscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'y') lo = yscale * utils::numeric(FLERR, arg[7], false, lmp);
if (axis == 'z') lo = zscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'z') lo = zscale * utils::numeric(FLERR, arg[7], false, lmp);
} }
if (strcmp(arg[8],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[8], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (axis == 'x') { if (axis == 'x') {
if (strcmp(arg[8],"INF") == 0) hi = BIG; if (strcmp(arg[8], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[0]; hi = BIG;
else hi = domain->boxhi_bound[0]; else if (domain->triclinic == 0)
hi = domain->boxhi[0];
else
hi = domain->boxhi_bound[0];
} }
if (axis == 'y') { if (axis == 'y') {
if (strcmp(arg[8],"INF") == 0) hi = BIG; if (strcmp(arg[8], "INF") == 0) hi = BIG;
if (domain->triclinic == 0) hi = domain->boxhi[1]; if (domain->triclinic == 0)
else hi = domain->boxhi_bound[1]; hi = domain->boxhi[1];
else
hi = domain->boxhi_bound[1];
} }
if (axis == 'z') { if (axis == 'z') {
if (strcmp(arg[8],"INF") == 0) hi = BIG; if (strcmp(arg[8], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[2]; hi = BIG;
else hi = domain->boxhi_bound[2]; else if (domain->triclinic == 0)
hi = domain->boxhi[2];
else
hi = domain->boxhi_bound[2];
} }
} else { } else {
if (axis == 'x') hi = xscale*utils::numeric(FLERR,arg[8],false,lmp); if (axis == 'x') hi = xscale * utils::numeric(FLERR, arg[8], false, lmp);
if (axis == 'y') hi = yscale*utils::numeric(FLERR,arg[8],false,lmp); if (axis == 'y') hi = yscale * utils::numeric(FLERR, arg[8], false, lmp);
if (axis == 'z') hi = zscale*utils::numeric(FLERR,arg[8],false,lmp); if (axis == 'z') hi = zscale * utils::numeric(FLERR, arg[8], false, lmp);
} }
// error check // error check
if (radiuslo < 0.0) error->all(FLERR,"Illegal radius in region cone command"); if (radiuslo < 0.0) error->all(FLERR, "Illegal radius in region cone command");
if (radiushi < 0.0) error->all(FLERR,"Illegal radius in region cone command"); if (radiushi < 0.0) error->all(FLERR, "Illegal radius in region cone command");
if (radiuslo == 0.0 && radiushi == 0.0) if (radiuslo == 0.0 && radiushi == 0.0)
error->all(FLERR,"Illegal radius in region cone command"); error->all(FLERR, "Illegal radius in region cone command");
if (hi == lo) error->all(FLERR,"Illegal cone length in region cone command"); if (hi == lo) error->all(FLERR, "Illegal cone length in region cone command");
// extent of cone // extent of cone
if (radiuslo > radiushi) maxradius = radiuslo; if (radiuslo > radiushi)
else maxradius = radiushi; maxradius = radiuslo;
else
maxradius = radiushi;
if (interior) { if (interior) {
bboxflag = 1; bboxflag = 1;
@ -149,22 +166,25 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) :
extent_zlo = lo; extent_zlo = lo;
extent_zhi = hi; extent_zhi = hi;
} }
} else bboxflag = 0; } else
bboxflag = 0;
// particle could be close to cone surface and 2 ends // particle could be close to cone surface and 2 ends
// particle can only touch surface and 1 end // particle can only touch surface and 1 end
cmax = 3; cmax = 3;
contact = new Contact[cmax]; contact = new Contact[cmax];
if (interior) tmax = 2; if (interior)
else tmax = 1; tmax = 2;
else
tmax = 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegCone::~RegCone() RegCone::~RegCone()
{ {
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -174,30 +194,36 @@ RegCone::~RegCone()
int RegCone::inside(double x, double y, double z) int RegCone::inside(double x, double y, double z)
{ {
double del1,del2,dist; double del1, del2, dist;
double currentradius; double currentradius;
if (axis == 'x') { if (axis == 'x') {
del1 = y - c1; del1 = y - c1;
del2 = z - c2; del2 = z - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x - lo) * (radiushi - radiuslo) / (hi - lo);
if (dist <= currentradius && x >= lo && x <= hi) return 1; if (dist <= currentradius && x >= lo && x <= hi)
else return 0; return 1;
else
return 0;
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x - c1; del1 = x - c1;
del2 = z - c2; del2 = z - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (y-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (y - lo) * (radiushi - radiuslo) / (hi - lo);
if (dist <= currentradius && y >= lo && y <= hi) return 1; if (dist <= currentradius && y >= lo && y <= hi)
else return 0; return 1;
else
return 0;
} else if (axis == 'z') { } else if (axis == 'z') {
del1 = x - c1; del1 = x - c1;
del2 = y - c2; del2 = y - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (z-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (z - lo) * (radiushi - radiuslo) / (hi - lo);
if (dist <= currentradius && z >= lo && z <= hi) return 1; if (dist <= currentradius && z >= lo && z <= hi)
else return 0; return 1;
else
return 0;
} }
return 0; return 0;
@ -213,16 +239,16 @@ int RegCone::inside(double x, double y, double z)
int RegCone::surface_interior(double *x, double cutoff) int RegCone::surface_interior(double *x, double cutoff)
{ {
double del1,del2,r,currentradius,delx,dely,delz,dist,delta; double del1, del2, r, currentradius, delx, dely, delz, dist, delta;
double surflo[3],surfhi[3],xs[3]; double surflo[3], surfhi[3], xs[3];
int n = 0; int n = 0;
if (axis == 'x') { if (axis == 'x') {
del1 = x[1] - c1; del1 = x[1] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[0] - lo) * (radiushi - radiuslo) / (hi - lo);
// x is exterior to cone // x is exterior to cone
@ -234,23 +260,22 @@ int RegCone::surface_interior(double *x, double cutoff)
if (r > 0.0 && !open_faces[2]) { if (r > 0.0 && !open_faces[2]) {
surflo[0] = lo; surflo[0] = lo;
surflo[1] = c1 + del1*radiuslo/r; surflo[1] = c1 + del1 * radiuslo / r;
surflo[2] = c2 + del2*radiuslo/r; surflo[2] = c2 + del2 * radiuslo / r;
surfhi[0] = hi; surfhi[0] = hi;
surfhi[1] = c1 + del1*radiushi/r; surfhi[1] = c1 + del1 * radiushi / r;
surfhi[2] = c2 + del2*radiushi/r; surfhi[2] = c2 + del2 * radiushi / r;
point_on_line_segment(surflo,surfhi,x,xs); point_on_line_segment(surflo, surfhi, x, xs);
delx = x[0] - xs[0]; delx = x[0] - xs[0];
dely = x[1] - xs[1]; dely = x[1] - xs[1];
delz = x[2] - xs[2]; delz = x[2] - xs[2];
dist = sqrt(delx*delx + dely*dely + delz*delz); dist = sqrt(delx * delx + dely * dely + delz * delz);
if (dist < cutoff) { if (dist < cutoff) {
contact[n].r = dist; contact[n].r = dist;
contact[n].delx = delx; contact[n].delx = delx;
contact[n].dely = dely; contact[n].dely = dely;
contact[n].delz = delz; contact[n].delz = delz;
contact[n].radius = -2.0*(radiuslo + (xs[0]-lo)* contact[n].radius = -2.0 * (radiuslo + (xs[0] - lo) * (radiushi - radiuslo) / (hi - lo));
(radiushi-radiuslo)/(hi-lo));
contact[n].iwall = 2; contact[n].iwall = 2;
n++; n++;
} }
@ -278,9 +303,8 @@ int RegCone::surface_interior(double *x, double cutoff)
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[1]-lo)* currentradius = radiuslo + (x[1] - lo) * (radiushi - radiuslo) / (hi - lo);
(radiushi-radiuslo)/(hi-lo);
// y is exterior to cone // y is exterior to cone
@ -291,25 +315,24 @@ int RegCone::surface_interior(double *x, double cutoff)
// surfhi = pt on outer circle of top end plane, same dir as y vs axis // surfhi = pt on outer circle of top end plane, same dir as y vs axis
if (r > 0.0 && !open_faces[2]) { if (r > 0.0 && !open_faces[2]) {
surflo[0] = c1 + del1*radiuslo/r; surflo[0] = c1 + del1 * radiuslo / r;
surflo[1] = lo; surflo[1] = lo;
surflo[2] = c2 + del2*radiuslo/r; surflo[2] = c2 + del2 * radiuslo / r;
surfhi[0] = c1 + del1*radiushi/r; surfhi[0] = c1 + del1 * radiushi / r;
surfhi[1] = hi; surfhi[1] = hi;
surfhi[2] = c2 + del2*radiushi/r; surfhi[2] = c2 + del2 * radiushi / r;
point_on_line_segment(surflo,surfhi,x,xs); point_on_line_segment(surflo, surfhi, x, xs);
delx = x[0] - xs[0]; delx = x[0] - xs[0];
dely = x[1] - xs[1]; dely = x[1] - xs[1];
delz = x[2] - xs[2]; delz = x[2] - xs[2];
dist = sqrt(delx*delx + dely*dely + delz*delz); dist = sqrt(delx * delx + dely * dely + delz * delz);
if (dist < cutoff) { if (dist < cutoff) {
contact[n].r = dist; contact[n].r = dist;
contact[n].delx = delx; contact[n].delx = delx;
contact[n].dely = dely; contact[n].dely = dely;
contact[n].delz = delz; contact[n].delz = delz;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].radius = -2.0*(radiuslo + (xs[1]-lo)* contact[n].radius = -2.0 * (radiuslo + (xs[1] - lo) * (radiushi - radiuslo) / (hi - lo));
(radiushi-radiuslo)/(hi-lo));
n++; n++;
} }
} }
@ -336,8 +359,8 @@ int RegCone::surface_interior(double *x, double cutoff)
} else { } else {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[1] - c2; del2 = x[1] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[2] - lo) * (radiushi - radiuslo) / (hi - lo);
// z is exterior to cone // z is exterior to cone
@ -348,31 +371,30 @@ int RegCone::surface_interior(double *x, double cutoff)
// surfhi = pt on outer circle of top end plane, same dir as z vs axis // surfhi = pt on outer circle of top end plane, same dir as z vs axis
if (r > 0.0 && !open_faces[2]) { if (r > 0.0 && !open_faces[2]) {
surflo[0] = c1 + del1*radiuslo/r; surflo[0] = c1 + del1 * radiuslo / r;
surflo[1] = c2 + del2*radiuslo/r; surflo[1] = c2 + del2 * radiuslo / r;
surflo[2] = lo; surflo[2] = lo;
surfhi[0] = c1 + del1*radiushi/r; surfhi[0] = c1 + del1 * radiushi / r;
surfhi[1] = c2 + del2*radiushi/r; surfhi[1] = c2 + del2 * radiushi / r;
surfhi[2] = hi; surfhi[2] = hi;
point_on_line_segment(surflo,surfhi,x,xs); point_on_line_segment(surflo, surfhi, x, xs);
delx = x[0] - xs[0]; delx = x[0] - xs[0];
dely = x[1] - xs[1]; dely = x[1] - xs[1];
delz = x[2] - xs[2]; delz = x[2] - xs[2];
dist = sqrt(delx*delx + dely*dely + delz*delz); dist = sqrt(delx * delx + dely * dely + delz * delz);
if (dist < cutoff) { if (dist < cutoff) {
contact[n].r = dist; contact[n].r = dist;
contact[n].delx = delx; contact[n].delx = delx;
contact[n].dely = dely; contact[n].dely = dely;
contact[n].delz = delz; contact[n].delz = delz;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].radius = -2.0*(radiuslo + (xs[2]-lo)* contact[n].radius = -2.0 * (radiuslo + (xs[2] - lo) * (radiushi - radiuslo) / (hi - lo));
(radiushi-radiuslo)/(hi-lo));
n++; n++;
} }
} }
delta = x[2] - lo; delta = x[2] - lo;
if (delta < cutoff && !open_faces[0]) { if (delta < cutoff && !open_faces[0]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delz = delta; contact[n].delz = delta;
contact[n].delx = contact[n].dely = 0.0; contact[n].delx = contact[n].dely = 0.0;
@ -381,7 +403,7 @@ int RegCone::surface_interior(double *x, double cutoff)
n++; n++;
} }
delta = hi - x[2]; delta = hi - x[2];
if (delta < cutoff && !open_faces[1]) { if (delta < cutoff && !open_faces[1]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delz = -delta; contact[n].delz = -delta;
contact[n].delx = contact[n].dely = 0.0; contact[n].delx = contact[n].dely = 0.0;
@ -402,14 +424,14 @@ int RegCone::surface_interior(double *x, double cutoff)
int RegCone::surface_exterior(double *x, double cutoff) int RegCone::surface_exterior(double *x, double cutoff)
{ {
double del1,del2,r,currentradius,distsq,distsqprev,crad; double del1, del2, r, currentradius, distsq, distsqprev, crad;
double corner1[3],corner2[3],corner3[3],corner4[3],xp[3],nearest[3]; double corner1[3], corner2[3], corner3[3], corner4[3], xp[3], nearest[3];
if (axis == 'x') { if (axis == 'x') {
del1 = x[1] - c1; del1 = x[1] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[0] - lo) * (radiushi - radiuslo) / (hi - lo);
// radius of curvature, only used for granular walls // radius of curvature, only used for granular walls
@ -418,8 +440,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
// x is far enough from cone that there is no contact // x is far enough from cone that there is no contact
// x is interior to cone // x is interior to cone
if (r >= maxradius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff) if (r >= maxradius + cutoff || x[0] <= lo - cutoff || x[0] >= hi + cutoff) return 0;
return 0;
if (r < currentradius && x[0] > lo && x[0] < hi) return 0; if (r < currentradius && x[0] > lo && x[0] < hi) return 0;
// x is exterior to cone or on its surface // x is exterior to cone or on its surface
@ -430,11 +451,11 @@ int RegCone::surface_exterior(double *x, double cutoff)
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
corner1[0] = lo; corner1[0] = lo;
corner1[1] = c1 + del1*radiuslo/r; corner1[1] = c1 + del1 * radiuslo / r;
corner1[2] = c2 + del2*radiuslo/r; corner1[2] = c2 + del2 * radiuslo / r;
corner2[0] = hi; corner2[0] = hi;
corner2[1] = c1 + del1*radiushi/r; corner2[1] = c1 + del1 * radiushi / r;
corner2[2] = c2 + del2*radiushi/r; corner2[2] = c2 + del2 * radiushi / r;
corner3[0] = lo; corner3[0] = lo;
corner3[1] = c1; corner3[1] = c1;
corner3[2] = c2; corner3[2] = c2;
@ -445,27 +466,27 @@ int RegCone::surface_exterior(double *x, double cutoff)
distsq = BIG; distsq = BIG;
if (!open_faces[2]) { if (!open_faces[2]) {
point_on_line_segment(corner1,corner2,x,xp); point_on_line_segment(corner1, corner2, x, xp);
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
crad = -2.0*(radiuslo + (nearest[0]-lo)*(radiushi-radiuslo)/(hi-lo)); crad = -2.0 * (radiuslo + (nearest[0] - lo) * (radiushi - radiuslo) / (hi - lo));
} }
if (!open_faces[0]) { if (!open_faces[0]) {
point_on_line_segment(corner1,corner3,x,xp); point_on_line_segment(corner1, corner3, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0.0; if (distsq < distsqprev) crad = 0.0;
} }
if (!open_faces[1]) { if (!open_faces[1]) {
point_on_line_segment(corner2,corner4,x,xp); point_on_line_segment(corner2, corner4, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0.0; if (distsq < distsqprev) crad = 0.0;
} }
if (distsq == BIG) return 0; if (distsq == BIG) return 0;
add_contact(0,x,nearest[0],nearest[1],nearest[2]); add_contact(0, x, nearest[0], nearest[1], nearest[2]);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
@ -474,8 +495,8 @@ int RegCone::surface_exterior(double *x, double cutoff)
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[1]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[1] - lo) * (radiushi - radiuslo) / (hi - lo);
// radius of curvature, only used for granular walls // radius of curvature, only used for granular walls
@ -484,8 +505,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
// y is far enough from cone that there is no contact // y is far enough from cone that there is no contact
// y is interior to cone // y is interior to cone
if (r >= maxradius+cutoff || if (r >= maxradius + cutoff || x[1] <= lo - cutoff || x[1] >= hi + cutoff) return 0;
x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0;
if (r < currentradius && x[1] > lo && x[1] < hi) return 0; if (r < currentradius && x[1] > lo && x[1] < hi) return 0;
// y is exterior to cone or on its surface // y is exterior to cone or on its surface
@ -495,12 +515,12 @@ int RegCone::surface_exterior(double *x, double cutoff)
// could be edge of cone // could be edge of cone
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
corner1[0] = c1 + del1*radiuslo/r; corner1[0] = c1 + del1 * radiuslo / r;
corner1[1] = lo; corner1[1] = lo;
corner1[2] = c2 + del2*radiuslo/r; corner1[2] = c2 + del2 * radiuslo / r;
corner2[0] = c1 + del1*radiushi/r; corner2[0] = c1 + del1 * radiushi / r;
corner2[1] = hi; corner2[1] = hi;
corner2[2] = c2 + del2*radiushi/r; corner2[2] = c2 + del2 * radiushi / r;
corner3[0] = c1; corner3[0] = c1;
corner3[1] = lo; corner3[1] = lo;
corner3[2] = c2; corner3[2] = c2;
@ -511,26 +531,26 @@ int RegCone::surface_exterior(double *x, double cutoff)
distsq = BIG; distsq = BIG;
if (!open_faces[2]) { if (!open_faces[2]) {
point_on_line_segment(corner1,corner2,x,xp); point_on_line_segment(corner1, corner2, x, xp);
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
crad = -2.0*(radiuslo + (nearest[1]-lo)*(radiushi-radiuslo)/(hi-lo)); crad = -2.0 * (radiuslo + (nearest[1] - lo) * (radiushi - radiuslo) / (hi - lo));
} }
if (!open_faces[0]) { if (!open_faces[0]) {
point_on_line_segment(corner1,corner3,x,xp); point_on_line_segment(corner1, corner3, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0; if (distsq < distsqprev) crad = 0;
} }
if (!open_faces[1]) { if (!open_faces[1]) {
point_on_line_segment(corner2,corner4,x,xp); point_on_line_segment(corner2, corner4, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0; if (distsq < distsqprev) crad = 0;
} }
add_contact(0,x,nearest[0],nearest[1],nearest[2]); add_contact(0, x, nearest[0], nearest[1], nearest[2]);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
@ -539,8 +559,8 @@ int RegCone::surface_exterior(double *x, double cutoff)
} else { } else {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[1] - c2; del2 = x[1] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo); currentradius = radiuslo + (x[2] - lo) * (radiushi - radiuslo) / (hi - lo);
// radius of curvature, only used for granular walls // radius of curvature, only used for granular walls
@ -549,8 +569,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
// z is far enough from cone that there is no contact // z is far enough from cone that there is no contact
// z is interior to cone // z is interior to cone
if (r >= maxradius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff) if (r >= maxradius + cutoff || x[2] <= lo - cutoff || x[2] >= hi + cutoff) return 0;
return 0;
if (r < currentradius && x[2] > lo && x[2] < hi) return 0; if (r < currentradius && x[2] > lo && x[2] < hi) return 0;
// z is exterior to cone or on its surface // z is exterior to cone or on its surface
@ -560,11 +579,11 @@ int RegCone::surface_exterior(double *x, double cutoff)
// could be edge of cone // could be edge of cone
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
corner1[0] = c1 + del1*radiuslo/r; corner1[0] = c1 + del1 * radiuslo / r;
corner1[1] = c2 + del2*radiuslo/r; corner1[1] = c2 + del2 * radiuslo / r;
corner1[2] = lo; corner1[2] = lo;
corner2[0] = c1 + del1*radiushi/r; corner2[0] = c1 + del1 * radiushi / r;
corner2[1] = c2 + del2*radiushi/r; corner2[1] = c2 + del2 * radiushi / r;
corner2[2] = hi; corner2[2] = hi;
corner3[0] = c1; corner3[0] = c1;
corner3[1] = c2; corner3[1] = c2;
@ -576,26 +595,26 @@ int RegCone::surface_exterior(double *x, double cutoff)
distsq = BIG; distsq = BIG;
if (!open_faces[2]) { if (!open_faces[2]) {
point_on_line_segment(corner1,corner2,x,xp); point_on_line_segment(corner1, corner2, x, xp);
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
crad = -2.0*(radiuslo + (nearest[2]-lo)*(radiushi-radiuslo)/(hi-lo)); crad = -2.0 * (radiuslo + (nearest[2] - lo) * (radiushi - radiuslo) / (hi - lo));
} }
if (!open_faces[0]) { if (!open_faces[0]) {
point_on_line_segment(corner1,corner3,x,xp); point_on_line_segment(corner1, corner3, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0; if (distsq < distsqprev) crad = 0;
} }
if (!open_faces[1]) { if (!open_faces[1]) {
point_on_line_segment(corner2,corner4,x,xp); point_on_line_segment(corner2, corner4, x, xp);
distsqprev = distsq; distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq); distsq = closest(x, xp, nearest, distsq);
if (distsq < distsqprev) crad = 0; if (distsq < distsqprev) crad = 0;
} }
add_contact(0,x,nearest[0],nearest[1],nearest[2]); add_contact(0, x, nearest[0], nearest[1], nearest[2]);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
@ -610,7 +629,7 @@ double RegCone::closest(double *x, double *near, double *nearest, double dsq)
double delx = x[0] - near[0]; double delx = x[0] - near[0];
double dely = x[1] - near[1]; double dely = x[1] - near[1];
double delz = x[2] - near[2]; double delz = x[2] - near[2];
double rsq = delx*delx + dely*dely + delz*delz; double rsq = delx * delx + dely * dely + delz * delz;
if (rsq >= dsq) return dsq; if (rsq >= dsq) return dsq;
nearest[0] = near[0]; nearest[0] = near[0];

428
src/region_cylinder.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -25,94 +24,97 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define BIG 1.0e20 static constexpr double BIG = 1.0e20;
enum{CONSTANT,VARIABLE};
enum { CONSTANT, VARIABLE };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) : RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), c1str(nullptr), c2str(nullptr), rstr(nullptr) Region(lmp, narg, arg), c1str(nullptr), c2str(nullptr), rstr(nullptr)
{ {
options(narg-8,&arg[8]); options(narg - 8, &arg[8]);
// check open face settings // check open face settings
if (openflag) if (openflag)
for (int i=3; i<6; i++) for (int i = 3; i < 6; i++)
if (open_faces[i]) error->all(FLERR,"Illegal region cylinder open face: {}", i+1); if (open_faces[i]) error->all(FLERR, "Illegal region cylinder open face: {}", i + 1);
if (strcmp(arg[2],"x") != 0 && strcmp(arg[2],"y") != 0 && strcmp(arg[2],"z") != 0) if (strcmp(arg[2], "x") != 0 && strcmp(arg[2], "y") != 0 && strcmp(arg[2], "z") != 0)
error->all(FLERR,"Illegal region cylinder axis: {}", arg[2]); error->all(FLERR, "Illegal region cylinder axis: {}", arg[2]);
axis = arg[2][0]; axis = arg[2][0];
if (axis == 'x') { if (axis == 'x') {
if (utils::strmatch(arg[3],"^v_")) { if (utils::strmatch(arg[3], "^v_")) {
c1str = utils::strdup(arg[3]+2); c1str = utils::strdup(arg[3] + 2);
c1 = 0.0; c1 = 0.0;
c1style = VARIABLE; c1style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c1 = yscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp);
c1style = CONSTANT; c1style = CONSTANT;
} }
if (utils::strmatch(arg[4],"^v_")) { if (utils::strmatch(arg[4], "^v_")) {
c2str = utils::strdup(arg[4]+2); c2str = utils::strdup(arg[4] + 2);
c2 = 0.0; c2 = 0.0;
c2style = VARIABLE; c2style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
c2style = CONSTANT; c2style = CONSTANT;
} }
} else if (axis == 'y') { } else if (axis == 'y') {
if (utils::strmatch(arg[3],"^v_")) { if (utils::strmatch(arg[3], "^v_")) {
c1str = utils::strdup(arg[3]+2); c1str = utils::strdup(arg[3] + 2);
c1 = 0.0; c1 = 0.0;
c1style = VARIABLE; c1style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c1style = CONSTANT; c1style = CONSTANT;
} }
if (utils::strmatch(arg[4],"^v_")) { if (utils::strmatch(arg[4], "^v_")) {
c2str = utils::strdup(arg[4]+2); c2str = utils::strdup(arg[4] + 2);
c2 = 0.0; c2 = 0.0;
c2style = VARIABLE; c2style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp);
c2style = CONSTANT; c2style = CONSTANT;
} }
} else if (axis == 'z') { } else if (axis == 'z') {
if (utils::strmatch(arg[3],"^v_")) { if (utils::strmatch(arg[3], "^v_")) {
c1str = utils::strdup(arg[3]+2); c1str = utils::strdup(arg[3] + 2);
c1 = 0.0; c1 = 0.0;
c1style = VARIABLE; c1style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp);
c1style = CONSTANT; c1style = CONSTANT;
} }
if (utils::strmatch(arg[4],"^v_")) { if (utils::strmatch(arg[4], "^v_")) {
c2str = utils::strdup(arg[4]+2); c2str = utils::strdup(arg[4] + 2);
c2 = 0.0; c2 = 0.0;
c2style = VARIABLE; c2style = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
c2 = yscale*utils::numeric(FLERR,arg[4],false,lmp); c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp);
c2style = CONSTANT; c2style = CONSTANT;
} }
} }
if (utils::strmatch(arg[5],"^v_")) { if (utils::strmatch(arg[5], "^v_")) {
rstr = utils::strdup(arg[5]+2); rstr = utils::strdup(arg[5] + 2);
radius = 0.0; radius = 0.0;
rstyle = VARIABLE; rstyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
radius = utils::numeric(FLERR,arg[5],false,lmp); radius = utils::numeric(FLERR, arg[5], false, lmp);
if (axis == 'x') radius *= yscale; if (axis == 'x')
else radius *= xscale; radius *= yscale;
else
radius *= xscale;
rstyle = CONSTANT; rstyle = CONSTANT;
} }
@ -121,57 +123,75 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
RegCylinder::shape_update(); RegCylinder::shape_update();
} }
if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { if (strcmp(arg[6], "INF") == 0 || strcmp(arg[6], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (axis == 'x') { if (axis == 'x') {
if (strcmp(arg[6],"INF") == 0) lo = -BIG; if (strcmp(arg[6], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[0]; lo = -BIG;
else lo = domain->boxlo_bound[0]; else if (domain->triclinic == 0)
lo = domain->boxlo[0];
else
lo = domain->boxlo_bound[0];
} }
if (axis == 'y') { if (axis == 'y') {
if (strcmp(arg[6],"INF") == 0) lo = -BIG; if (strcmp(arg[6], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[1]; lo = -BIG;
else lo = domain->boxlo_bound[1]; else if (domain->triclinic == 0)
lo = domain->boxlo[1];
else
lo = domain->boxlo_bound[1];
} }
if (axis == 'z') { if (axis == 'z') {
if (strcmp(arg[6],"INF") == 0) lo = -BIG; if (strcmp(arg[6], "INF") == 0)
else if (domain->triclinic == 0) lo = domain->boxlo[2]; lo = -BIG;
else lo = domain->boxlo_bound[2]; else if (domain->triclinic == 0)
lo = domain->boxlo[2];
else
lo = domain->boxlo_bound[2];
} }
} else { } else {
if (axis == 'x') lo = xscale*utils::numeric(FLERR,arg[6],false,lmp); if (axis == 'x') lo = xscale * utils::numeric(FLERR, arg[6], false, lmp);
if (axis == 'y') lo = yscale*utils::numeric(FLERR,arg[6],false,lmp); if (axis == 'y') lo = yscale * utils::numeric(FLERR, arg[6], false, lmp);
if (axis == 'z') lo = zscale*utils::numeric(FLERR,arg[6],false,lmp); if (axis == 'z') lo = zscale * utils::numeric(FLERR, arg[6], false, lmp);
} }
if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (axis == 'x') { if (axis == 'x') {
if (strcmp(arg[7],"INF") == 0) hi = BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[0]; hi = BIG;
else hi = domain->boxhi_bound[0]; else if (domain->triclinic == 0)
hi = domain->boxhi[0];
else
hi = domain->boxhi_bound[0];
} }
if (axis == 'y') { if (axis == 'y') {
if (strcmp(arg[7],"INF") == 0) hi = BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[1]; hi = BIG;
else hi = domain->boxhi_bound[1]; else if (domain->triclinic == 0)
hi = domain->boxhi[1];
else
hi = domain->boxhi_bound[1];
} }
if (axis == 'z') { if (axis == 'z') {
if (strcmp(arg[7],"INF") == 0) hi = BIG; if (strcmp(arg[7], "INF") == 0)
else if (domain->triclinic == 0) hi = domain->boxhi[2]; hi = BIG;
else hi = domain->boxhi_bound[2]; else if (domain->triclinic == 0)
hi = domain->boxhi[2];
else
hi = domain->boxhi_bound[2];
} }
} else { } else {
if (axis == 'x') hi = xscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'x') hi = xscale * utils::numeric(FLERR, arg[7], false, lmp);
if (axis == 'y') hi = yscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'y') hi = yscale * utils::numeric(FLERR, arg[7], false, lmp);
if (axis == 'z') hi = zscale*utils::numeric(FLERR,arg[7],false,lmp); if (axis == 'z') hi = zscale * utils::numeric(FLERR, arg[7], false, lmp);
} }
// error check // error check
if (radius <= 0.0) error->all(FLERR,"Illegal radius {} in region cylinder command", radius); if (radius <= 0.0) error->all(FLERR, "Illegal radius {} in region cylinder command", radius);
// extent of cylinder // extent of cylinder
// for variable radius, uses initial radius // for variable radius, uses initial radius
@ -202,25 +222,28 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
extent_zlo = lo; extent_zlo = lo;
extent_zhi = hi; extent_zhi = hi;
} }
} else bboxflag = 0; } else
bboxflag = 0;
// particle could be close to cylinder surface and 2 ends // particle could be close to cylinder surface and 2 ends
// particle can only touch surface and 1 end // particle can only touch surface and 1 end
cmax = 3; cmax = 3;
contact = new Contact[cmax]; contact = new Contact[cmax];
if (interior) tmax = 2; if (interior)
else tmax = 1; tmax = 2;
else
tmax = 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegCylinder::~RegCylinder() RegCylinder::~RegCylinder()
{ {
delete [] c1str; delete[] c1str;
delete [] c2str; delete[] c2str;
delete [] rstr; delete[] rstr;
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -238,27 +261,33 @@ void RegCylinder::init()
int RegCylinder::inside(double x, double y, double z) int RegCylinder::inside(double x, double y, double z)
{ {
double del1,del2,dist; double del1, del2, dist;
int inside; int inside;
if (axis == 'x') { if (axis == 'x') {
del1 = y - c1; del1 = y - c1;
del2 = z - c2; del2 = z - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
if (dist <= radius && x >= lo && x <= hi) inside = 1; if (dist <= radius && x >= lo && x <= hi)
else inside = 0; inside = 1;
else
inside = 0;
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x - c1; del1 = x - c1;
del2 = z - c2; del2 = z - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
if (dist <= radius && y >= lo && y <= hi) inside = 1; if (dist <= radius && y >= lo && y <= hi)
else inside = 0; inside = 1;
else
inside = 0;
} else { } else {
del1 = x - c1; del1 = x - c1;
del2 = y - c2; del2 = y - c2;
dist = sqrt(del1*del1 + del2*del2); dist = sqrt(del1 * del1 + del2 * del2);
if (dist <= radius && z >= lo && z <= hi) inside = 1; if (dist <= radius && z >= lo && z <= hi)
else inside = 0; inside = 1;
else
inside = 0;
} }
return inside; return inside;
@ -274,14 +303,14 @@ int RegCylinder::inside(double x, double y, double z)
int RegCylinder::surface_interior(double *x, double cutoff) int RegCylinder::surface_interior(double *x, double cutoff)
{ {
double del1,del2,r,delta; double del1, del2, r, delta;
int n = 0; int n = 0;
if (axis == 'x') { if (axis == 'x') {
del1 = x[1] - c1; del1 = x[1] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// x is exterior to cylinder // x is exterior to cylinder
@ -293,9 +322,9 @@ int RegCylinder::surface_interior(double *x, double cutoff)
if (delta < cutoff && r > 0.0 && !open_faces[2]) { if (delta < cutoff && r > 0.0 && !open_faces[2]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delx = 0.0; contact[n].delx = 0.0;
contact[n].dely = del1*(1.0-radius/r); contact[n].dely = del1 * (1.0 - radius / r);
contact[n].delz = del2*(1.0-radius/r); contact[n].delz = del2 * (1.0 - radius / r);
contact[n].radius = -2.0*radius; contact[n].radius = -2.0 * radius;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].varflag = 1; contact[n].varflag = 1;
n++; n++;
@ -324,7 +353,7 @@ int RegCylinder::surface_interior(double *x, double cutoff)
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// y is exterior to cylinder // y is exterior to cylinder
@ -335,10 +364,10 @@ int RegCylinder::surface_interior(double *x, double cutoff)
delta = radius - r; delta = radius - r;
if (delta < cutoff && r > 0.0 && !open_faces[2]) { if (delta < cutoff && r > 0.0 && !open_faces[2]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delx = del1*(1.0-radius/r); contact[n].delx = del1 * (1.0 - radius / r);
contact[n].dely = 0.0; contact[n].dely = 0.0;
contact[n].delz = del2*(1.0-radius/r); contact[n].delz = del2 * (1.0 - radius / r);
contact[n].radius = -2.0*radius; contact[n].radius = -2.0 * radius;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].varflag = 1; contact[n].varflag = 1;
n++; n++;
@ -367,7 +396,7 @@ int RegCylinder::surface_interior(double *x, double cutoff)
} else { } else {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[1] - c2; del2 = x[1] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// z is exterior to cylinder // z is exterior to cylinder
@ -378,10 +407,10 @@ int RegCylinder::surface_interior(double *x, double cutoff)
delta = radius - r; delta = radius - r;
if (delta < cutoff && r > 0.0 && !open_faces[2]) { if (delta < cutoff && r > 0.0 && !open_faces[2]) {
contact[n].r = delta; contact[n].r = delta;
contact[n].delx = del1*(1.0-radius/r); contact[n].delx = del1 * (1.0 - radius / r);
contact[n].dely = del2*(1.0-radius/r); contact[n].dely = del2 * (1.0 - radius / r);
contact[n].delz = 0.0; contact[n].delz = 0.0;
contact[n].radius = -2.0*radius; contact[n].radius = -2.0 * radius;
contact[n].iwall = 2; contact[n].iwall = 2;
contact[n].varflag = 1; contact[n].varflag = 1;
n++; n++;
@ -419,12 +448,12 @@ int RegCylinder::surface_interior(double *x, double cutoff)
int RegCylinder::surface_exterior(double *x, double cutoff) int RegCylinder::surface_exterior(double *x, double cutoff)
{ {
double del1,del2,r; double del1, del2, r;
double xp,yp,zp; double xp, yp, zp;
double dx, dr, dr2, d2, d2prev; double dx, dr, dr2, d2, d2prev;
// radius of curvature for granular // radius of curvature for granular
// 0 for flat surfaces (infinite case), 2*radius for curved portion // 0 for flat surfaces (infinite case), 2*radius for curved portion
double crad = 0.0; double crad = 0.0;
int varflag = 0; int varflag = 0;
@ -432,12 +461,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (axis == 'x') { if (axis == 'x') {
del1 = x[1] - c1; del1 = x[1] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// x is far enough from cylinder that there is no contact // x is far enough from cylinder that there is no contact
// x is interior to cylinder // x is interior to cylinder
if (r >= radius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff) return 0; if (r >= radius + cutoff || x[0] <= lo - cutoff || x[0] >= hi + cutoff) return 0;
if (r < radius && x[0] > lo && x[0] < hi) return 0; if (r < radius && x[0] > lo && x[0] < hi) return 0;
// x is exterior to cylinder or on its surface // x is exterior to cylinder or on its surface
@ -448,48 +477,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
d2prev = BIG; d2prev = BIG;
if (!openflag) { if (!openflag) {
if (r > radius) { if (r > radius) {
yp = c1 + del1*radius/r; yp = c1 + del1 * radius / r;
zp = c2 + del2*radius/r; zp = c2 + del2 * radius / r;
crad = 2.0*radius; crad = 2.0 * radius;
varflag = 1; varflag = 1;
} else { } else {
yp = x[1]; yp = x[1];
zp = x[2]; zp = x[2];
} }
if (x[0] < lo) xp = lo; if (x[0] < lo)
else if (x[0] > hi) xp = hi; xp = lo;
else xp = x[0]; else if (x[0] > hi)
xp = hi;
else
xp = x[0];
} else { } else {
// closest point on curved surface // closest point on curved surface
dr = r - radius; dr = r - radius;
dr2 = dr*dr; dr2 = dr * dr;
if (!open_faces[2]) { if (!open_faces[2]) {
yp = c1 + del1*radius/r; yp = c1 + del1 * radius / r;
zp = c2 + del2*radius/r; zp = c2 + del2 * radius / r;
if (x[0] < lo) { if (x[0] < lo) {
dx = lo-x[0]; dx = lo - x[0];
xp = lo; xp = lo;
} } else if (x[0] > hi) {
else if (x[0] > hi) { dx = x[0] - hi;
dx = x[0]-hi;
xp = hi; xp = hi;
} } else {
else {
dx = 0; dx = 0;
xp = x[0]; xp = x[0];
} }
d2 = d2prev = dr2 + dx*dx; d2 = d2prev = dr2 + dx * dx;
} }
// closest point on bottom cap // closest point on bottom cap
if (!open_faces[0]) { if (!open_faces[0]) {
dx = lo - x[0]; dx = lo - x[0];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
xp = lo; xp = lo;
if (r < radius) { if (r < radius) {
@ -504,8 +536,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (!open_faces[1]) { if (!open_faces[1]) {
dx = hi - x[0]; dx = hi - x[0];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
xp = hi; xp = hi;
if (r < radius) { if (r < radius) {
@ -516,7 +550,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} }
} }
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].varflag = varflag; contact[0].varflag = varflag;
contact[0].iwall = 0; contact[0].iwall = 0;
@ -526,12 +560,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} else if (axis == 'y') { } else if (axis == 'y') {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[2] - c2; del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// y is far enough from cylinder that there is no contact // y is far enough from cylinder that there is no contact
// y is interior to cylinder // y is interior to cylinder
if (r >= radius+cutoff || x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0; if (r >= radius + cutoff || x[1] <= lo - cutoff || x[1] >= hi + cutoff) return 0;
if (r < radius && x[1] > lo && x[1] < hi) return 0; if (r < radius && x[1] > lo && x[1] < hi) return 0;
// y is exterior to cylinder or on its surface // y is exterior to cylinder or on its surface
@ -542,48 +576,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
d2prev = BIG; d2prev = BIG;
if (!openflag) { if (!openflag) {
if (r > radius) { if (r > radius) {
xp = c1 + del1*radius/r; xp = c1 + del1 * radius / r;
zp = c2 + del2*radius/r; zp = c2 + del2 * radius / r;
crad = 2.0*radius; crad = 2.0 * radius;
varflag = 1; varflag = 1;
} else { } else {
xp = x[0]; xp = x[0];
zp = x[2]; zp = x[2];
} }
if (x[1] < lo) yp = lo; if (x[1] < lo)
else if (x[1] > hi) yp = hi; yp = lo;
else yp = x[1]; else if (x[1] > hi)
yp = hi;
else
yp = x[1];
} else { } else {
// closest point on curved surface // closest point on curved surface
dr = r - radius; dr = r - radius;
dr2 = dr*dr; dr2 = dr * dr;
if (!open_faces[2]) { if (!open_faces[2]) {
xp = c1 + del1*radius/r; xp = c1 + del1 * radius / r;
zp = c2 + del2*radius/r; zp = c2 + del2 * radius / r;
if (x[1] < lo) { if (x[1] < lo) {
dx = lo-x[1]; dx = lo - x[1];
yp = lo; yp = lo;
} } else if (x[1] > hi) {
else if (x[1] > hi) { dx = x[1] - hi;
dx = x[1]-hi;
yp = hi; yp = hi;
} } else {
else {
dx = 0; dx = 0;
yp = x[1]; yp = x[1];
} }
d2 = d2prev = dr2 + dx*dx; d2 = d2prev = dr2 + dx * dx;
} }
// closest point on bottom cap // closest point on bottom cap
if (!open_faces[0]) { if (!open_faces[0]) {
dx = lo - x[1]; dx = lo - x[1];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
yp = lo; yp = lo;
if (r < radius) { if (r < radius) {
@ -598,8 +635,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (!open_faces[1]) { if (!open_faces[1]) {
dx = hi - x[1]; dx = hi - x[1];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
yp = hi; yp = hi;
if (r < radius) { if (r < radius) {
@ -610,7 +649,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} }
} }
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].varflag = varflag; contact[0].varflag = varflag;
contact[0].iwall = 0; contact[0].iwall = 0;
@ -620,12 +659,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} else { } else {
del1 = x[0] - c1; del1 = x[0] - c1;
del2 = x[1] - c2; del2 = x[1] - c2;
r = sqrt(del1*del1 + del2*del2); r = sqrt(del1 * del1 + del2 * del2);
// z is far enough from cylinder that there is no contact // z is far enough from cylinder that there is no contact
// z is interior to cylinder // z is interior to cylinder
if (r >= radius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff) return 0; if (r >= radius + cutoff || x[2] <= lo - cutoff || x[2] >= hi + cutoff) return 0;
if (r < radius && x[2] > lo && x[2] < hi) return 0; if (r < radius && x[2] > lo && x[2] < hi) return 0;
// z is exterior to cylinder or on its surface // z is exterior to cylinder or on its surface
@ -636,46 +675,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
d2prev = BIG; d2prev = BIG;
if (!openflag) { if (!openflag) {
if (r > radius) { if (r > radius) {
xp = c1 + del1*radius/r; xp = c1 + del1 * radius / r;
yp = c2 + del2*radius/r; yp = c2 + del2 * radius / r;
crad = 2.0*radius; crad = 2.0 * radius;
varflag = 1; varflag = 1;
} else { } else {
xp = x[0]; xp = x[0];
yp = x[1]; yp = x[1];
} }
if (x[2] < lo) zp = lo; if (x[2] < lo)
else if (x[2] > hi) zp = hi; zp = lo;
else zp = x[2]; else if (x[2] > hi)
zp = hi;
else
zp = x[2];
} else { } else {
// closest point on curved surface // closest point on curved surface
dr = r - radius; dr = r - radius;
dr2 = dr*dr; dr2 = dr * dr;
if (!open_faces[2]) { if (!open_faces[2]) {
xp = c1 + del1*radius/r; xp = c1 + del1 * radius / r;
yp = c2 + del2*radius/r; yp = c2 + del2 * radius / r;
if (x[2] < lo) { if (x[2] < lo) {
dx = lo-x[2]; dx = lo - x[2];
zp = lo; zp = lo;
} else if (x[2] > hi) { } else if (x[2] > hi) {
dx = x[2]-hi; dx = x[2] - hi;
zp = hi; zp = hi;
} else { } else {
dx = 0; dx = 0;
zp = x[2]; zp = x[2];
} }
d2prev = dr2 + dx*dx; d2prev = dr2 + dx * dx;
} }
// closest point on bottom cap // closest point on bottom cap
if (!open_faces[0]) { if (!open_faces[0]) {
dx = lo - x[2]; dx = lo - x[2];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
zp = lo; zp = lo;
if (r < radius) { if (r < radius) {
@ -690,8 +734,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (!open_faces[1]) { if (!open_faces[1]) {
dx = hi - x[2]; dx = hi - x[2];
if (r < radius) d2 = dx*dx; if (r < radius)
else d2 = dr2 + dx*dx; d2 = dx * dx;
else
d2 = dr2 + dx * dx;
if (d2 < d2prev) { if (d2 < d2prev) {
zp = hi; zp = hi;
if (r < radius) { if (r < radius) {
@ -702,7 +748,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
} }
} }
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].radius = crad; contact[0].radius = crad;
contact[0].varflag = varflag; contact[0].varflag = varflag;
contact[0].iwall = 0; contact[0].iwall = 0;
@ -721,22 +767,21 @@ void RegCylinder::shape_update()
if (c2style == VARIABLE) c2 = input->variable->compute_equal(c2var); if (c2style == VARIABLE) c2 = input->variable->compute_equal(c2var);
if (rstyle == VARIABLE) { if (rstyle == VARIABLE) {
radius = input->variable->compute_equal(rvar); radius = input->variable->compute_equal(rvar);
if (radius < 0.0) if (radius < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
error->one(FLERR,"Variable evaluation in region gave bad value");
} }
if (axis == 'x') { if (axis == 'x') {
if (c1style == VARIABLE) c1 *= yscale; if (c1style == VARIABLE) c1 *= yscale;
if (c2style == VARIABLE) c2 *= zscale; if (c2style == VARIABLE) c2 *= zscale;
if (rstyle == VARIABLE) radius *= yscale; if (rstyle == VARIABLE) radius *= yscale;
} else if (axis == 'y') { } else if (axis == 'y') {
if (c1style == VARIABLE) c1 *= xscale; if (c1style == VARIABLE) c1 *= xscale;
if (c2style == VARIABLE) c2 *= zscale; if (c2style == VARIABLE) c2 *= zscale;
if (rstyle == VARIABLE) radius *= xscale; if (rstyle == VARIABLE) radius *= xscale;
} else { // axis == 'z' } else { // axis == 'z'
if (c1style == VARIABLE) c1 *= xscale; if (c1style == VARIABLE) c1 *= xscale;
if (c2style == VARIABLE) c2 *= yscale; if (c2style == VARIABLE) c2 *= yscale;
if (rstyle == VARIABLE) radius *= xscale; if (rstyle == VARIABLE) radius *= xscale;
} }
} }
@ -748,30 +793,26 @@ void RegCylinder::variable_check()
{ {
if (c1style == VARIABLE) { if (c1style == VARIABLE) {
c1var = input->variable->find(c1str); c1var = input->variable->find(c1str);
if (c1var < 0) if (c1var < 0) error->all(FLERR, "Variable name for region cylinder does not exist");
error->all(FLERR,"Variable name for region cylinder does not exist");
if (!input->variable->equalstyle(c1var)) if (!input->variable->equalstyle(c1var))
error->all(FLERR,"Variable for region cylinder is invalid style"); error->all(FLERR, "Variable for region cylinder is invalid style");
} }
if (c2style == VARIABLE) { if (c2style == VARIABLE) {
c2var = input->variable->find(c2str); c2var = input->variable->find(c2str);
if (c2var < 0) if (c2var < 0) error->all(FLERR, "Variable name for region cylinder does not exist");
error->all(FLERR,"Variable name for region cylinder does not exist");
if (!input->variable->equalstyle(c2var)) if (!input->variable->equalstyle(c2var))
error->all(FLERR,"Variable for region cylinder is invalid style"); error->all(FLERR, "Variable for region cylinder is invalid style");
} }
if (rstyle == VARIABLE) { if (rstyle == VARIABLE) {
rvar = input->variable->find(rstr); rvar = input->variable->find(rstr);
if (rvar < 0) if (rvar < 0) error->all(FLERR, "Variable name for region cylinder does not exist");
error->all(FLERR,"Variable name for region cylinder does not exist");
if (!input->variable->equalstyle(rvar)) if (!input->variable->equalstyle(rvar))
error->all(FLERR,"Variable for region cylinder is invalid style"); error->all(FLERR, "Variable for region cylinder is invalid style");
} }
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Set values needed to calculate velocity due to shape changes. Set values needed to calculate velocity due to shape changes.
These values do not depend on the contact, so this function is These values do not depend on the contact, so this function is
@ -795,35 +836,34 @@ void RegCylinder::set_velocity_shape()
xcenter[2] = 0; xcenter[2] = 0;
} }
forward_transform(xcenter[0], xcenter[1], xcenter[2]); forward_transform(xcenter[0], xcenter[1], xcenter[2]);
if (update->ntimestep > 0) rprev = prev[4]; if (update->ntimestep > 0)
else rprev = radius; rprev = prev[4];
else
rprev = radius;
prev[4] = radius; prev[4] = radius;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
add velocity due to shape change to wall velocity add velocity due to shape change to wall velocity
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void RegCylinder::velocity_contact_shape(double *vwall, double *xc) void RegCylinder::velocity_contact_shape(double *vwall, double *xc)
{ {
double delx, dely, delz; // Displacement of contact point in x,y,z double delx, dely, delz; // Displacement of contact point in x,y,z
if (axis == 'x') { if (axis == 'x') {
delx = 0; delx = 0;
dely = (xc[1] - xcenter[1])*(1 - rprev/radius); dely = (xc[1] - xcenter[1]) * (1 - rprev / radius);
delz = (xc[2] - xcenter[2])*(1 - rprev/radius); delz = (xc[2] - xcenter[2]) * (1 - rprev / radius);
} else if (axis == 'y') { } else if (axis == 'y') {
delx = (xc[0] - xcenter[0])*(1 - rprev/radius); delx = (xc[0] - xcenter[0]) * (1 - rprev / radius);
dely = 0; dely = 0;
delz = (xc[2] - xcenter[2])*(1 - rprev/radius); delz = (xc[2] - xcenter[2]) * (1 - rprev / radius);
} else { } else {
delx = (xc[0] - xcenter[0])*(1 - rprev/radius); delx = (xc[0] - xcenter[0]) * (1 - rprev / radius);
dely = (xc[1] - xcenter[1])*(1 - rprev/radius); dely = (xc[1] - xcenter[1]) * (1 - rprev / radius);
delz = 0; delz = 0;
} }
vwall[0] += delx/update->dt; vwall[0] += delx / update->dt;
vwall[1] += dely/update->dt; vwall[1] += dely / update->dt;
vwall[2] += delz/update->dt; vwall[2] += delz / update->dt;
} }

339
src/region_prism.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -26,105 +25,126 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define BIG 1.0e20 static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
{ {
options(narg-11,&arg[11]); options(narg - 11, &arg[11]);
if (strcmp(arg[2],"INF") == 0 || strcmp(arg[2],"EDGE") == 0) { if (strcmp(arg[2], "INF") == 0 || strcmp(arg[2], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[2],"INF") == 0) xlo = -BIG; if (strcmp(arg[2], "INF") == 0)
else xlo = domain->boxlo[0]; xlo = -BIG;
} else xlo = xscale*utils::numeric(FLERR,arg[2],false,lmp); else
xlo = domain->boxlo[0];
} else
xlo = xscale * utils::numeric(FLERR, arg[2], false, lmp);
if (strcmp(arg[3],"INF") == 0 || strcmp(arg[3],"EDGE") == 0) { if (strcmp(arg[3], "INF") == 0 || strcmp(arg[3], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[3],"INF") == 0) xhi = BIG; if (strcmp(arg[3], "INF") == 0)
else xhi = domain->boxhi[0]; xhi = BIG;
} else xhi = xscale*utils::numeric(FLERR,arg[3],false,lmp); else
xhi = domain->boxhi[0];
} else
xhi = xscale * utils::numeric(FLERR, arg[3], false, lmp);
if (strcmp(arg[4],"INF") == 0 || strcmp(arg[4],"EDGE") == 0) { if (strcmp(arg[4], "INF") == 0 || strcmp(arg[4], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[4],"INF") == 0) ylo = -BIG; if (strcmp(arg[4], "INF") == 0)
else ylo = domain->boxlo[1]; ylo = -BIG;
} else ylo = yscale*utils::numeric(FLERR,arg[4],false,lmp); else
ylo = domain->boxlo[1];
} else
ylo = yscale * utils::numeric(FLERR, arg[4], false, lmp);
if (strcmp(arg[5],"INF") == 0 || strcmp(arg[5],"EDGE") == 0) { if (strcmp(arg[5], "INF") == 0 || strcmp(arg[5], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[5],"INF") == 0) yhi = BIG; if (strcmp(arg[5], "INF") == 0)
else yhi = domain->boxhi[1]; yhi = BIG;
} else yhi = yscale*utils::numeric(FLERR,arg[5],false,lmp); else
yhi = domain->boxhi[1];
} else
yhi = yscale * utils::numeric(FLERR, arg[5], false, lmp);
if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { if (strcmp(arg[6], "INF") == 0 || strcmp(arg[6], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[6],"INF") == 0) zlo = -BIG; if (strcmp(arg[6], "INF") == 0)
else zlo = domain->boxlo[2]; zlo = -BIG;
} else zlo = zscale*utils::numeric(FLERR,arg[6],false,lmp); else
zlo = domain->boxlo[2];
} else
zlo = zscale * utils::numeric(FLERR, arg[6], false, lmp);
if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); error->all(FLERR, "Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[7],"INF") == 0) zhi = BIG; if (strcmp(arg[7], "INF") == 0)
else zhi = domain->boxhi[2]; zhi = BIG;
} else zhi = zscale*utils::numeric(FLERR,arg[7],false,lmp); else
zhi = domain->boxhi[2];
} else
zhi = zscale * utils::numeric(FLERR, arg[7], false, lmp);
xy = xscale*utils::numeric(FLERR,arg[8],false,lmp); xy = xscale * utils::numeric(FLERR, arg[8], false, lmp);
xz = xscale*utils::numeric(FLERR,arg[9],false,lmp); xz = xscale * utils::numeric(FLERR, arg[9], false, lmp);
yz = yscale*utils::numeric(FLERR,arg[10],false,lmp); yz = yscale * utils::numeric(FLERR, arg[10], false, lmp);
// error check // error check
// prism cannot be 0 thickness in any dim, else inverse blows up // prism cannot be 0 thickness in any dim, else inverse blows up
// non-zero tilt values cannot be used if either dim is INF on both ends // non-zero tilt values cannot be used if either dim is INF on both ends
if (xlo >= xhi) error->all(FLERR,"Illegal region prism xlo: {} >= xhi: {}", xlo, xhi); if (xlo >= xhi) error->all(FLERR, "Illegal region prism xlo: {} >= xhi: {}", xlo, xhi);
if (ylo >= yhi) error->all(FLERR,"Illegal region prism ylo: {} >= yhi: {}", ylo, yhi); if (ylo >= yhi) error->all(FLERR, "Illegal region prism ylo: {} >= yhi: {}", ylo, yhi);
if (zlo >= zhi) error->all(FLERR,"Illegal region prism zlo: {} >= zhi: {}", zlo ,zhi); if (zlo >= zhi) error->all(FLERR, "Illegal region prism zlo: {} >= zhi: {}", zlo, zhi);
if (xy != 0.0 && xlo == -BIG && xhi == BIG) if (xy != 0.0 && xlo == -BIG && xhi == BIG)
error->all(FLERR,"Illegal region prism non-zero xy tilt with infinite x size"); error->all(FLERR, "Illegal region prism non-zero xy tilt with infinite x size");
if (xy != 0.0 && ylo == -BIG && yhi == BIG) if (xy != 0.0 && ylo == -BIG && yhi == BIG)
error->all(FLERR,"Illegal region prism non-zero xy tilt with infinite y size"); error->all(FLERR, "Illegal region prism non-zero xy tilt with infinite y size");
if (xz != 0.0 && xlo == -BIG && xhi == BIG) if (xz != 0.0 && xlo == -BIG && xhi == BIG)
error->all(FLERR,"Illegal region prism non-zero xz tilt with infinite x size"); error->all(FLERR, "Illegal region prism non-zero xz tilt with infinite x size");
if (xz != 0.0 && zlo == -BIG && zhi == BIG) if (xz != 0.0 && zlo == -BIG && zhi == BIG)
error->all(FLERR,"Illegal region prism non-zero xz tilt with infinite z size"); error->all(FLERR, "Illegal region prism non-zero xz tilt with infinite z size");
if (yz != 0.0 && ylo == -BIG && yhi == BIG) if (yz != 0.0 && ylo == -BIG && yhi == BIG)
error->all(FLERR,"Illegal region prism non-zero yz tilt with infinite y size"); error->all(FLERR, "Illegal region prism non-zero yz tilt with infinite y size");
if (yz != 0.0 && zlo == -BIG && zhi == BIG) if (yz != 0.0 && zlo == -BIG && zhi == BIG)
error->all(FLERR,"Illegal region prism non-zero yz tilt with infinite z size"); error->all(FLERR, "Illegal region prism non-zero yz tilt with infinite z size");
// extent of prism // extent of prism
if (interior) { if (interior) {
bboxflag = 1; bboxflag = 1;
extent_xlo = MIN(xlo,xlo+xy); extent_xlo = MIN(xlo, xlo + xy);
extent_xlo = MIN(extent_xlo,extent_xlo+xz); extent_xlo = MIN(extent_xlo, extent_xlo + xz);
extent_ylo = MIN(ylo,ylo+yz); extent_ylo = MIN(ylo, ylo + yz);
extent_zlo = zlo; extent_zlo = zlo;
extent_xhi = MAX(xhi,xhi+xy); extent_xhi = MAX(xhi, xhi + xy);
extent_xhi = MAX(extent_xhi,extent_xhi+xz); extent_xhi = MAX(extent_xhi, extent_xhi + xz);
extent_yhi = MAX(yhi,yhi+yz); extent_yhi = MAX(yhi, yhi + yz);
extent_zhi = zhi; extent_zhi = zhi;
} else bboxflag = 0; } else
bboxflag = 0;
// particle could be close to all 6 planes // particle could be close to all 6 planes
// particle can only touch 3 planes // particle can only touch 3 planes
cmax = 6; cmax = 6;
contact = new Contact[cmax]; contact = new Contact[cmax];
if (interior) tmax = 3; if (interior)
else tmax = 1; tmax = 3;
else
tmax = 1;
// h = transformation matrix from tilt coords (0-1) to box coords (xyz) // h = transformation matrix from tilt coords (0-1) to box coords (xyz)
// columns of h are edge vectors of tilted box // columns of h are edge vectors of tilted box
@ -140,26 +160,26 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
h[1][2] = yz; h[1][2] = yz;
h[2][2] = zhi - zlo; h[2][2] = zhi - zlo;
hinv[0][0] = 1.0/h[0][0]; hinv[0][0] = 1.0 / h[0][0];
hinv[0][1] = -h[0][1] / (h[0][0]*h[1][1]); hinv[0][1] = -h[0][1] / (h[0][0] * h[1][1]);
hinv[0][2] = (h[0][1]*h[1][2] - h[0][2]*h[1][1]) / (h[0][0]*h[1][1]*h[2][2]); hinv[0][2] = (h[0][1] * h[1][2] - h[0][2] * h[1][1]) / (h[0][0] * h[1][1] * h[2][2]);
hinv[1][1] = 1.0/h[1][1]; hinv[1][1] = 1.0 / h[1][1];
hinv[1][2] = -h[1][2] / (h[1][1]*h[2][2]); hinv[1][2] = -h[1][2] / (h[1][1] * h[2][2]);
hinv[2][2] = 1.0/h[2][2]; hinv[2][2] = 1.0 / h[2][2];
// corners = 8 corner points of prism // corners = 8 corner points of prism
// order = x varies fastest, then y, finally z // order = x varies fastest, then y, finally z
// clo/chi = lo and hi corner pts of prism // clo/chi = lo and hi corner pts of prism
a[0] = xhi-xlo; a[0] = xhi - xlo;
a[1] = 0.0; a[1] = 0.0;
a[2] = 0.0; a[2] = 0.0;
b[0] = xy; b[0] = xy;
b[1] = yhi-ylo; b[1] = yhi - ylo;
b[2] = 0.0; b[2] = 0.0;
c[0] = xz; c[0] = xz;
c[1] = yz; c[1] = yz;
c[2] = zhi-zlo; c[2] = zhi - zlo;
clo[0] = corners[0][0] = xlo; clo[0] = corners[0][0] = xlo;
clo[1] = corners[0][1] = ylo; clo[1] = corners[0][1] = ylo;
@ -196,19 +216,18 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
// face = 6 inward-facing unit normals to prism faces // face = 6 inward-facing unit normals to prism faces
// order = xy plane, xz plane, yz plane // order = xy plane, xz plane, yz plane
MathExtra::cross3(a,b,face[0]); MathExtra::cross3(a, b, face[0]);
MathExtra::cross3(b,a,face[1]); MathExtra::cross3(b, a, face[1]);
MathExtra::cross3(c,a,face[2]); MathExtra::cross3(c, a, face[2]);
MathExtra::cross3(a,c,face[3]); MathExtra::cross3(a, c, face[3]);
MathExtra::cross3(b,c,face[4]); MathExtra::cross3(b, c, face[4]);
MathExtra::cross3(c,b,face[5]); MathExtra::cross3(c, b, face[5]);
// remap open face indices to be consistent // remap open face indices to be consistent
if (openflag) { if (openflag) {
int temp[6]; int temp[6];
for (int i = 0; i < 6; i++) for (int i = 0; i < 6; i++) temp[i] = open_faces[i];
temp[i] = open_faces[i];
open_faces[0] = temp[4]; open_faces[0] = temp[4];
open_faces[1] = temp[5]; open_faces[1] = temp[5];
open_faces[2] = temp[2]; open_faces[2] = temp[2];
@ -223,27 +242,51 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
// verts in each tri are ordered so that right-hand rule gives inward norm // verts in each tri are ordered so that right-hand rule gives inward norm
// order = xy plane, xz plane, yz plane // order = xy plane, xz plane, yz plane
tri[0][0] = 0; tri[0][1] = 1; tri[0][2] = 3; tri[0][0] = 0;
tri[1][0] = 0; tri[1][1] = 3; tri[1][2] = 2; tri[0][1] = 1;
tri[2][0] = 4; tri[2][1] = 7; tri[2][2] = 5; tri[0][2] = 3;
tri[3][0] = 4; tri[3][1] = 6; tri[3][2] = 7; tri[1][0] = 0;
tri[1][1] = 3;
tri[1][2] = 2;
tri[2][0] = 4;
tri[2][1] = 7;
tri[2][2] = 5;
tri[3][0] = 4;
tri[3][1] = 6;
tri[3][2] = 7;
tri[4][0] = 0; tri[4][1] = 4; tri[4][2] = 5; tri[4][0] = 0;
tri[5][0] = 0; tri[5][1] = 5; tri[5][2] = 1; tri[4][1] = 4;
tri[6][0] = 2; tri[6][1] = 7; tri[6][2] = 6; tri[4][2] = 5;
tri[7][0] = 2; tri[7][1] = 3; tri[7][2] = 7; tri[5][0] = 0;
tri[5][1] = 5;
tri[5][2] = 1;
tri[6][0] = 2;
tri[6][1] = 7;
tri[6][2] = 6;
tri[7][0] = 2;
tri[7][1] = 3;
tri[7][2] = 7;
tri[8][0] = 2; tri[8][1] = 6; tri[8][2] = 4; tri[8][0] = 2;
tri[9][0] = 2; tri[9][1] = 4; tri[9][2] = 0; tri[8][1] = 6;
tri[10][0] = 1; tri[10][1] = 5; tri[10][2] = 7; tri[8][2] = 4;
tri[11][0] = 1; tri[11][1] = 7; tri[11][2] = 3; tri[9][0] = 2;
tri[9][1] = 4;
tri[9][2] = 0;
tri[10][0] = 1;
tri[10][1] = 5;
tri[10][2] = 7;
tri[11][0] = 1;
tri[11][1] = 7;
tri[11][2] = 3;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegPrism::~RegPrism() RegPrism::~RegPrism()
{ {
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -258,12 +301,11 @@ RegPrism::~RegPrism()
int RegPrism::inside(double x, double y, double z) int RegPrism::inside(double x, double y, double z)
{ {
double a = hinv[0][0]*(x-xlo) + hinv[0][1]*(y-ylo) + hinv[0][2]*(z-zlo); double a = hinv[0][0] * (x - xlo) + hinv[0][1] * (y - ylo) + hinv[0][2] * (z - zlo);
double b = hinv[1][1]*(y-ylo) + hinv[1][2]*(z-zlo); double b = hinv[1][1] * (y - ylo) + hinv[1][2] * (z - zlo);
double c = hinv[2][2]*(z-zlo); double c = hinv[2][2] * (z - zlo);
if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0) if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0) return 1;
return 1;
return 0; return 0;
} }
@ -283,10 +325,12 @@ int RegPrism::surface_interior(double *x, double cutoff)
// x is exterior to prism // x is exterior to prism
for (i = 0; i < 6; i++) { for (i = 0; i < 6; i++) {
if (i % 2) corner = chi; if (i % 2)
else corner = clo; corner = chi;
dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + else
(x[2]-corner[2])*face[i][2]; corner = clo;
dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] +
(x[2] - corner[2]) * face[i][2];
if (dot < 0.0) return 0; if (dot < 0.0) return 0;
} }
@ -296,15 +340,17 @@ int RegPrism::surface_interior(double *x, double cutoff)
for (i = 0; i < 6; i++) { for (i = 0; i < 6; i++) {
if (open_faces[i]) continue; if (open_faces[i]) continue;
if (i % 2) corner = chi; if (i % 2)
else corner = clo; corner = chi;
dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + else
(x[2]-corner[2])*face[i][2]; corner = clo;
dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] +
(x[2] - corner[2]) * face[i][2];
if (dot < cutoff) { if (dot < cutoff) {
contact[n].r = dot; contact[n].r = dot;
contact[n].delx = dot*face[i][0]; contact[n].delx = dot * face[i][0];
contact[n].dely = dot*face[i][1]; contact[n].dely = dot * face[i][1];
contact[n].delz = dot*face[i][2]; contact[n].delz = dot * face[i][2];
contact[n].radius = 0; contact[n].radius = 0;
contact[n].iwall = i; contact[n].iwall = i;
n++; n++;
@ -325,25 +371,29 @@ int RegPrism::surface_exterior(double *x, double cutoff)
int i; int i;
double dot; double dot;
double *corner; double *corner;
double xp,yp,zp; double xp, yp, zp;
// x is far enough from prism that there is no contact // x is far enough from prism that there is no contact
for (i = 0; i < 6; i++) { for (i = 0; i < 6; i++) {
if (i % 2) corner = chi; if (i % 2)
else corner = clo; corner = chi;
dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + else
(x[2]-corner[2])*face[i][2]; corner = clo;
dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] +
(x[2] - corner[2]) * face[i][2];
if (dot <= -cutoff) return 0; if (dot <= -cutoff) return 0;
} }
// x is interior to prism // x is interior to prism
for (i = 0; i < 6; i++) { for (i = 0; i < 6; i++) {
if (i % 2) corner = chi; if (i % 2)
else corner = clo; corner = chi;
dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + else
(x[2]-corner[2])*face[i][2]; corner = clo;
dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] +
(x[2] - corner[2]) * face[i][2];
if (dot <= 0.0) break; if (dot <= 0.0) break;
} }
@ -354,8 +404,8 @@ int RegPrism::surface_exterior(double *x, double cutoff)
// could be edge or corner pt of prism // could be edge or corner pt of prism
// do not add contact point if r >= cutoff // do not add contact point if r >= cutoff
find_nearest(x,xp,yp,zp); find_nearest(x, xp, yp, zp);
add_contact(0,x,xp,yp,zp); add_contact(0, x, xp, yp, zp);
contact[0].radius = 0; contact[0].radius = 0;
contact[0].iwall = 0; contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1; if (contact[0].r < cutoff) return 1;
@ -369,8 +419,8 @@ int RegPrism::surface_exterior(double *x, double cutoff)
void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp) void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp)
{ {
int i,j,k,iface; int i, j, k, iface;
double xproj[3],xline[3],nearest[3]; double xproj[3], xline[3], nearest[3];
double dot; double dot;
// generate successive xnear points, one nearest to x is (xp,yp,zp) // generate successive xnear points, one nearest to x is (xp,yp,zp)
@ -384,27 +434,25 @@ void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp)
double distsq = BIG; double distsq = BIG;
for (int itri = 0; itri < 12; itri++) { for (int itri = 0; itri < 12; itri++) {
iface = itri/2; iface = itri / 2;
if (open_faces[iface]) continue; if (open_faces[iface]) continue;
i = tri[itri][0]; i = tri[itri][0];
j = tri[itri][1]; j = tri[itri][1];
k = tri[itri][2]; k = tri[itri][2];
dot = (x[0]-corners[i][0])*face[iface][0] + dot = (x[0] - corners[i][0]) * face[iface][0] + (x[1] - corners[i][1]) * face[iface][1] +
(x[1]-corners[i][1])*face[iface][1] + (x[2] - corners[i][2]) * face[iface][2];
(x[2]-corners[i][2])*face[iface][2]; xproj[0] = x[0] - dot * face[iface][0];
xproj[0] = x[0] - dot*face[iface][0]; xproj[1] = x[1] - dot * face[iface][1];
xproj[1] = x[1] - dot*face[iface][1]; xproj[2] = x[2] - dot * face[iface][2];
xproj[2] = x[2] - dot*face[iface][2]; if (inside_tri(xproj, corners[i], corners[j], corners[k], face[iface])) {
if (inside_tri(xproj,corners[i],corners[j],corners[k],face[iface])) { distsq = closest(x, xproj, nearest, distsq);
distsq = closest(x,xproj,nearest,distsq); } else {
} point_on_line_segment(corners[i], corners[j], xproj, xline);
else { distsq = closest(x, xline, nearest, distsq);
point_on_line_segment(corners[i],corners[j],xproj,xline); point_on_line_segment(corners[j], corners[k], xproj, xline);
distsq = closest(x,xline,nearest,distsq); distsq = closest(x, xline, nearest, distsq);
point_on_line_segment(corners[j],corners[k],xproj,xline); point_on_line_segment(corners[i], corners[k], xproj, xline);
distsq = closest(x,xline,nearest,distsq); distsq = closest(x, xline, nearest, distsq);
point_on_line_segment(corners[i],corners[k],xproj,xline);
distsq = closest(x,xline,nearest,distsq);
} }
} }
@ -422,25 +470,24 @@ void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp)
if xproduct dot norm < 0.0 for any of 3 edges, then x is outside triangle if xproduct dot norm < 0.0 for any of 3 edges, then x is outside triangle
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int RegPrism::inside_tri(double *x, double *v1, double *v2, double *v3, int RegPrism::inside_tri(double *x, double *v1, double *v2, double *v3, double *norm)
double *norm)
{ {
double edge[3],pvec[3],xproduct[3]; double edge[3], pvec[3], xproduct[3];
MathExtra::sub3(v2,v1,edge); MathExtra::sub3(v2, v1, edge);
MathExtra::sub3(x,v1,pvec); MathExtra::sub3(x, v1, pvec);
MathExtra::cross3(edge,pvec,xproduct); MathExtra::cross3(edge, pvec, xproduct);
if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; if (MathExtra::dot3(xproduct, norm) < 0.0) return 0;
MathExtra::sub3(v3,v2,edge); MathExtra::sub3(v3, v2, edge);
MathExtra::sub3(x,v2,pvec); MathExtra::sub3(x, v2, pvec);
MathExtra::cross3(edge,pvec,xproduct); MathExtra::cross3(edge, pvec, xproduct);
if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; if (MathExtra::dot3(xproduct, norm) < 0.0) return 0;
MathExtra::sub3(v1,v3,edge); MathExtra::sub3(v1, v3, edge);
MathExtra::sub3(x,v3,pvec); MathExtra::sub3(x, v3, pvec);
MathExtra::cross3(edge,pvec,xproduct); MathExtra::cross3(edge, pvec, xproduct);
if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; if (MathExtra::dot3(xproduct, norm) < 0.0) return 0;
return 1; return 1;
} }
@ -452,7 +499,7 @@ double RegPrism::closest(double *x, double *near, double *nearest, double dsq)
double delx = x[0] - near[0]; double delx = x[0] - near[0];
double dely = x[1] - near[1]; double dely = x[1] - near[1];
double delz = x[2] - near[2]; double delz = x[2] - near[2];
double rsq = delx*delx + dely*dely + delz*delz; double rsq = delx * delx + dely * dely + delz * delz;
if (rsq >= dsq) return dsq; if (rsq >= dsq) return dsq;
nearest[0] = near[0]; nearest[0] = near[0];

119
src/region_sphere.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -23,52 +22,52 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{CONSTANT,VARIABLE}; enum { CONSTANT, VARIABLE };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr) Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr)
{ {
options(narg-6,&arg[6]); options(narg - 6, &arg[6]);
if (utils::strmatch(arg[2],"^v_")) { if (utils::strmatch(arg[2], "^v_")) {
xstr = utils::strdup(arg[2]+2); xstr = utils::strdup(arg[2] + 2);
xc = 0.0; xc = 0.0;
xstyle = VARIABLE; xstyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
xc = xscale*utils::numeric(FLERR,arg[2],false,lmp); xc = xscale * utils::numeric(FLERR, arg[2], false, lmp);
xstyle = CONSTANT; xstyle = CONSTANT;
} }
if (utils::strmatch(arg[3],"^v_")) { if (utils::strmatch(arg[3], "^v_")) {
ystr = utils::strdup(arg[3]+2); ystr = utils::strdup(arg[3] + 2);
yc = 0.0; yc = 0.0;
ystyle = VARIABLE; ystyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
yc = yscale*utils::numeric(FLERR,arg[3],false,lmp); yc = yscale * utils::numeric(FLERR, arg[3], false, lmp);
ystyle = CONSTANT; ystyle = CONSTANT;
} }
if (utils::strmatch(arg[4],"^v_")) { if (utils::strmatch(arg[4], "^v_")) {
zstr = utils::strdup(arg[4]+2); zstr = utils::strdup(arg[4] + 2);
zc = 0.0; zc = 0.0;
zstyle = VARIABLE; zstyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
zc = zscale*utils::numeric(FLERR,arg[4],false,lmp); zc = zscale * utils::numeric(FLERR, arg[4], false, lmp);
zstyle = CONSTANT; zstyle = CONSTANT;
} }
if (utils::strmatch(arg[5],"^v_")) { if (utils::strmatch(arg[5], "^v_")) {
rstr = utils::strdup(arg[5]+2); rstr = utils::strdup(arg[5] + 2);
radius = 0.0; radius = 0.0;
rstyle = VARIABLE; rstyle = VARIABLE;
varshape = 1; varshape = 1;
} else { } else {
radius = xscale*utils::numeric(FLERR,arg[5],false,lmp); radius = xscale * utils::numeric(FLERR, arg[5], false, lmp);
rstyle = CONSTANT; rstyle = CONSTANT;
} }
@ -79,7 +78,7 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
// error check // error check
if (radius < 0.0) error->all(FLERR,"Illegal region sphere radius: {}", radius); if (radius < 0.0) error->all(FLERR, "Illegal region sphere radius: {}", radius);
// extent of sphere // extent of sphere
// for variable radius, uses initial radius and origin for variable center // for variable radius, uses initial radius and origin for variable center
@ -92,7 +91,8 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
extent_yhi = yc + radius; extent_yhi = yc + radius;
extent_zlo = zc - radius; extent_zlo = zc - radius;
extent_zhi = zc + radius; extent_zhi = zc + radius;
} else bboxflag = 0; } else
bboxflag = 0;
cmax = 1; cmax = 1;
contact = new Contact[cmax]; contact = new Contact[cmax];
@ -103,11 +103,11 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
RegSphere::~RegSphere() RegSphere::~RegSphere()
{ {
delete [] xstr; delete[] xstr;
delete [] ystr; delete[] ystr;
delete [] zstr; delete[] zstr;
delete [] rstr; delete[] rstr;
delete [] contact; delete[] contact;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -128,7 +128,7 @@ int RegSphere::inside(double x, double y, double z)
double delx = x - xc; double delx = x - xc;
double dely = y - yc; double dely = y - yc;
double delz = z - zc; double delz = z - zc;
double r = sqrt(delx*delx + dely*dely + delz*delz); double r = sqrt(delx * delx + dely * dely + delz * delz);
if (r <= radius) return 1; if (r <= radius) return 1;
return 0; return 0;
@ -146,15 +146,15 @@ int RegSphere::surface_interior(double *x, double cutoff)
double delx = x[0] - xc; double delx = x[0] - xc;
double dely = x[1] - yc; double dely = x[1] - yc;
double delz = x[2] - zc; double delz = x[2] - zc;
double r = sqrt(delx*delx + dely*dely + delz*delz); double r = sqrt(delx * delx + dely * dely + delz * delz);
if (r > radius || r == 0.0) return 0; if (r > radius || r == 0.0) return 0;
double delta = radius - r; double delta = radius - r;
if (delta < cutoff) { if (delta < cutoff) {
contact[0].r = delta; contact[0].r = delta;
contact[0].delx = delx*(1.0-radius/r); contact[0].delx = delx * (1.0 - radius / r);
contact[0].dely = dely*(1.0-radius/r); contact[0].dely = dely * (1.0 - radius / r);
contact[0].delz = delz*(1.0-radius/r); contact[0].delz = delz * (1.0 - radius / r);
contact[0].radius = -radius; contact[0].radius = -radius;
contact[0].iwall = 0; contact[0].iwall = 0;
contact[0].varflag = 1; contact[0].varflag = 1;
@ -174,15 +174,15 @@ int RegSphere::surface_exterior(double *x, double cutoff)
double delx = x[0] - xc; double delx = x[0] - xc;
double dely = x[1] - yc; double dely = x[1] - yc;
double delz = x[2] - zc; double delz = x[2] - zc;
double r = sqrt(delx*delx + dely*dely + delz*delz); double r = sqrt(delx * delx + dely * dely + delz * delz);
if (r < radius) return 0; if (r < radius) return 0;
double delta = r - radius; double delta = r - radius;
if (delta < cutoff) { if (delta < cutoff) {
contact[0].r = delta; contact[0].r = delta;
contact[0].delx = delx*(1.0-radius/r); contact[0].delx = delx * (1.0 - radius / r);
contact[0].dely = dely*(1.0-radius/r); contact[0].dely = dely * (1.0 - radius / r);
contact[0].delz = delz*(1.0-radius/r); contact[0].delz = delz * (1.0 - radius / r);
contact[0].radius = radius; contact[0].radius = radius;
contact[0].iwall = 0; contact[0].iwall = 0;
contact[0].varflag = 1; contact[0].varflag = 1;
@ -197,19 +197,15 @@ int RegSphere::surface_exterior(double *x, double cutoff)
void RegSphere::shape_update() void RegSphere::shape_update()
{ {
if (xstyle == VARIABLE) if (xstyle == VARIABLE) xc = xscale * input->variable->compute_equal(xvar);
xc = xscale * input->variable->compute_equal(xvar);
if (ystyle == VARIABLE) if (ystyle == VARIABLE) yc = yscale * input->variable->compute_equal(yvar);
yc = yscale * input->variable->compute_equal(yvar);
if (zstyle == VARIABLE) if (zstyle == VARIABLE) zc = zscale * input->variable->compute_equal(zvar);
zc = zscale * input->variable->compute_equal(zvar);
if (rstyle == VARIABLE) { if (rstyle == VARIABLE) {
radius = xscale * input->variable->compute_equal(rvar); radius = xscale * input->variable->compute_equal(rvar);
if (radius < 0.0) if (radius < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
error->one(FLERR,"Variable evaluation in region gave bad value");
} }
} }
@ -221,34 +217,30 @@ void RegSphere::variable_check()
{ {
if (xstyle == VARIABLE) { if (xstyle == VARIABLE) {
xvar = input->variable->find(xstr); xvar = input->variable->find(xstr);
if (xvar < 0) if (xvar < 0) error->all(FLERR, "Variable name for region sphere does not exist");
error->all(FLERR,"Variable name for region sphere does not exist");
if (!input->variable->equalstyle(xvar)) if (!input->variable->equalstyle(xvar))
error->all(FLERR,"Variable for region sphere is invalid style"); error->all(FLERR, "Variable for region sphere is invalid style");
} }
if (ystyle == VARIABLE) { if (ystyle == VARIABLE) {
yvar = input->variable->find(ystr); yvar = input->variable->find(ystr);
if (yvar < 0) if (yvar < 0) error->all(FLERR, "Variable name for region sphere does not exist");
error->all(FLERR,"Variable name for region sphere does not exist");
if (!input->variable->equalstyle(yvar)) if (!input->variable->equalstyle(yvar))
error->all(FLERR,"Variable for region sphere is invalid style"); error->all(FLERR, "Variable for region sphere is invalid style");
} }
if (zstyle == VARIABLE) { if (zstyle == VARIABLE) {
zvar = input->variable->find(zstr); zvar = input->variable->find(zstr);
if (zvar < 0) if (zvar < 0) error->all(FLERR, "Variable name for region sphere does not exist");
error->all(FLERR,"Variable name for region sphere does not exist");
if (!input->variable->equalstyle(zvar)) if (!input->variable->equalstyle(zvar))
error->all(FLERR,"Variable for region sphere is invalid style"); error->all(FLERR, "Variable for region sphere is invalid style");
} }
if (rstyle == VARIABLE) { if (rstyle == VARIABLE) {
rvar = input->variable->find(rstr); rvar = input->variable->find(rstr);
if (rvar < 0) if (rvar < 0) error->all(FLERR, "Variable name for region sphere does not exist");
error->all(FLERR,"Variable name for region sphere does not exist");
if (!input->variable->equalstyle(rvar)) if (!input->variable->equalstyle(rvar))
error->all(FLERR,"Variable for region sphere is invalid style"); error->all(FLERR, "Variable for region sphere is invalid style");
} }
} }
@ -265,27 +257,26 @@ void RegSphere::set_velocity_shape()
xcenter[1] = yc; xcenter[1] = yc;
xcenter[2] = zc; xcenter[2] = zc;
forward_transform(xcenter[0], xcenter[1], xcenter[2]); forward_transform(xcenter[0], xcenter[1], xcenter[2]);
if (update->ntimestep > 0) rprev = prev[4]; if (update->ntimestep > 0)
else rprev = radius; rprev = prev[4];
else
rprev = radius;
prev[4] = radius; prev[4] = radius;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
add velocity due to shape change to wall velocity add velocity due to shape change to wall velocity
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void RegSphere::velocity_contact_shape(double *vwall, double *xc) void RegSphere::velocity_contact_shape(double *vwall, double *xc)
{ {
double delx, dely, delz; // Displacement of contact point in x,y,z double delx, dely, delz; // Displacement of contact point in x,y,z
delx = (xc[0] - xcenter[0])*(1 - rprev/radius); delx = (xc[0] - xcenter[0]) * (1 - rprev / radius);
dely = (xc[1] - xcenter[1])*(1 - rprev/radius); dely = (xc[1] - xcenter[1]) * (1 - rprev / radius);
delz = (xc[2] - xcenter[2])*(1 - rprev/radius); delz = (xc[2] - xcenter[2]) * (1 - rprev / radius);
vwall[0] += delx/update->dt; vwall[0] += delx / update->dt;
vwall[1] += dely/update->dt; vwall[1] += dely / update->dt;
vwall[2] += delz/update->dt; vwall[2] += delz / update->dt;
} }