diff --git a/src/region_block.cpp b/src/region_block.cpp index c20d3abb17..8b50521d7c 100644 --- a/src/region_block.cpp +++ b/src/region_block.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -24,100 +23,125 @@ 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) : - 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; - 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[2],"INF") == 0) xlo = -BIG; - else if (domain->triclinic == 0) xlo = domain->boxlo[0]; - else xlo = domain->boxlo_bound[0]; - } else if (utils::strmatch(arg[2],"^v_")) { - xlostr = utils::strdup(arg[2]+2); - xlo = 0.0; - xlostyle = VARIABLE; - varshape = 1; - } else xlo = xscale*utils::numeric(FLERR,arg[2],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[2], "INF") == 0) + xlo = -BIG; + else if (domain->triclinic == 0) + xlo = domain->boxlo[0]; + else + xlo = domain->boxlo_bound[0]; + } else if (utils::strmatch(arg[2], "^v_")) { + xlostr = utils::strdup(arg[2] + 2); + xlo = 0.0; + xlostyle = VARIABLE; + varshape = 1; + } else + xlo = xscale * utils::numeric(FLERR, arg[2], false, lmp); 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[3],"INF") == 0) xhi = BIG; - else if (domain->triclinic == 0) xhi = domain->boxhi[0]; - else xhi = domain->boxhi_bound[0]; - } else if (utils::strmatch(arg[3],"^v_")) { - xhistr = utils::strdup(arg[3]+2); - xhi = 0.0; - xhistyle = VARIABLE; - varshape = 1; - } else xhi = xscale*utils::numeric(FLERR,arg[3],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[3], "INF") == 0) + xhi = BIG; + else if (domain->triclinic == 0) + xhi = domain->boxhi[0]; + else + xhi = domain->boxhi_bound[0]; + } else if (utils::strmatch(arg[3], "^v_")) { + xhistr = utils::strdup(arg[3] + 2); + xhi = 0.0; + xhistyle = VARIABLE; + varshape = 1; + } else + xhi = xscale * utils::numeric(FLERR, arg[3], false, lmp); 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[4],"INF") == 0) ylo = -BIG; - else if (domain->triclinic == 0) ylo = domain->boxlo[1]; - else ylo = domain->boxlo_bound[1]; - } else if (utils::strmatch(arg[4],"^v_")) { - ylostr = utils::strdup(arg[4]+2); - ylo = 0.0; - ylostyle = VARIABLE; - varshape = 1; - } else ylo = yscale*utils::numeric(FLERR,arg[4],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[4], "INF") == 0) + ylo = -BIG; + else if (domain->triclinic == 0) + ylo = domain->boxlo[1]; + else + ylo = domain->boxlo_bound[1]; + } else if (utils::strmatch(arg[4], "^v_")) { + ylostr = utils::strdup(arg[4] + 2); + ylo = 0.0; + ylostyle = VARIABLE; + varshape = 1; + } else + ylo = yscale * utils::numeric(FLERR, arg[4], false, lmp); 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[5],"INF") == 0) yhi = BIG; - else if (domain->triclinic == 0) yhi = domain->boxhi[1]; - else yhi = domain->boxhi_bound[1]; - } else if (utils::strmatch(arg[5],"^v_")) { - yhistr = utils::strdup(arg[5]+2); - yhi = 0.0; - yhistyle = VARIABLE; - varshape = 1; - } else yhi = yscale*utils::numeric(FLERR,arg[5],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[5], "INF") == 0) + yhi = BIG; + else if (domain->triclinic == 0) + yhi = domain->boxhi[1]; + else + yhi = domain->boxhi_bound[1]; + } else if (utils::strmatch(arg[5], "^v_")) { + yhistr = utils::strdup(arg[5] + 2); + yhi = 0.0; + yhistyle = VARIABLE; + varshape = 1; + } else + yhi = yscale * utils::numeric(FLERR, arg[5], false, lmp); 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[6],"INF") == 0) zlo = -BIG; - else if (domain->triclinic == 0) zlo = domain->boxlo[2]; - else zlo = domain->boxlo_bound[2]; - } else if (utils::strmatch(arg[6],"^v_")) { - zlostr = utils::strdup(arg[6]+2); - zlo = 0.0; - zlostyle = VARIABLE; - varshape = 1; - } else zlo = zscale*utils::numeric(FLERR,arg[6],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[6], "INF") == 0) + zlo = -BIG; + else if (domain->triclinic == 0) + zlo = domain->boxlo[2]; + else + zlo = domain->boxlo_bound[2]; + } else if (utils::strmatch(arg[6], "^v_")) { + zlostr = utils::strdup(arg[6] + 2); + zlo = 0.0; + zlostyle = VARIABLE; + varshape = 1; + } else + zlo = zscale * utils::numeric(FLERR, arg[6], false, lmp); 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[7],"INF") == 0) zhi = BIG; - else if (domain->triclinic == 0) zhi = domain->boxhi[2]; - else zhi = domain->boxhi_bound[2]; - } else if (utils::strmatch(arg[7],"^v_")) { - zhistr = utils::strdup(arg[7]+2); - zhi = 0.0; - zhistyle = VARIABLE; - varshape = 1; - } else zhi = zscale*utils::numeric(FLERR,arg[7],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[7], "INF") == 0) + zhi = BIG; + else if (domain->triclinic == 0) + zhi = domain->boxhi[2]; + else + zhi = domain->boxhi_bound[2]; + } else if (utils::strmatch(arg[7], "^v_")) { + zhistr = utils::strdup(arg[7] + 2); + zhi = 0.0; + zhistyle = VARIABLE; + varshape = 1; + } else + zhi = zscale * utils::numeric(FLERR, arg[7], false, lmp); if (varshape) { variable_check(); @@ -126,9 +150,9 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : // error check - 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 (zlo > zhi) error->all(FLERR,"Illegal region block zlo: {} >= zhi: {}", zlo, zhi); + 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 (zlo > zhi) error->all(FLERR, "Illegal region block zlo: {} >= zhi: {}", zlo, zhi); // extent of block @@ -140,15 +164,18 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : extent_yhi = yhi; extent_zlo = zlo; extent_zhi = zhi; - } else bboxflag = 0; + } else + bboxflag = 0; // particle could be close to all 6 planes // particle can only touch 3 planes cmax = 6; contact = new Contact[cmax]; - if (interior) tmax = 3; - else tmax = 1; + if (interior) + tmax = 3; + else + tmax = 1; // open face data structs @@ -235,13 +262,13 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : RegBlock::~RegBlock() { if (copymode) return; - delete [] xlostr; - delete [] xhistr; - delete [] ylostr; - delete [] yhistr; - delete [] zlostr; - delete [] zhistr; - delete [] contact; + delete[] xlostr; + delete[] xhistr; + delete[] ylostr; + delete[] yhistr; + delete[] zlostr; + delete[] zhistr; + delete[] contact; } /* ---------------------------------------------------------------------- */ @@ -259,8 +286,7 @@ void RegBlock::init() int RegBlock::inside(double x, double y, double z) { - if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) - return 1; + if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) return 1; return 0; } @@ -277,8 +303,7 @@ int RegBlock::surface_interior(double *x, double cutoff) // x is exterior to block - if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || - x[2] < zlo || x[2] > zhi) return 0; + if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || x[2] < zlo || x[2] > zhi) return 0; // 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) { - double xp,yp,zp; - double xc,yc,zc,dist,mindist; + double xp, yp, zp; + double xc, yc, zc, dist, mindist; // x is far enough from block that there is no contact // x is interior to block - if (x[0] <= xlo-cutoff || x[0] >= xhi+cutoff || - x[1] <= ylo-cutoff || x[1] >= yhi+cutoff || - x[2] <= zlo-cutoff || x[2] >= zhi+cutoff) return 0; - if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && - x[2] > zlo && x[2] < zhi) return 0; + if (x[0] <= xlo - cutoff || x[0] >= xhi + cutoff || x[1] <= ylo - cutoff || + x[1] >= yhi + cutoff || x[2] <= zlo - cutoff || x[2] >= zhi + cutoff) + return 0; + if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && x[2] > zlo && x[2] < zhi) return 0; // x is exterior to block or on its surface // 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 if (!openflag) { - if (x[0] < xlo) xp = xlo; - else if (x[0] > xhi) xp = xhi; - else xp = x[0]; - if (x[1] < ylo) yp = ylo; - 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]; + if (x[0] < xlo) + xp = xlo; + else if (x[0] > xhi) + xp = xhi; + else + xp = x[0]; + if (x[1] < ylo) + yp = ylo; + 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 { mindist = BIG; for (int i = 0; i < 6; i++) { 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) { xp = xc; 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; if (contact[0].r < cutoff) return 1; return 0; @@ -403,28 +436,17 @@ int RegBlock::surface_exterior(double *x, double cutoff) change region shape via variable evaluation ------------------------------------------------------------------------- */ -void RegBlock::shape_update() // addition +void RegBlock::shape_update() // addition { - if (xlostyle == VARIABLE) - 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 (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 (xlostyle == VARIABLE) 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 (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) - error->one(FLERR,"Variable evaluation in region gave bad value"); + error->one(FLERR, "Variable evaluation in region gave bad value"); // face[0] @@ -489,54 +511,48 @@ void RegBlock::shape_update() // addition error check on existence of variable ------------------------------------------------------------------------- */ -void RegBlock::variable_check() // addition +void RegBlock::variable_check() // addition { if (xlostyle == VARIABLE) { xlovar = input->variable->find(xlostr); - if (xlovar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (xlovar < 0) error->all(FLERR, "Variable name for region block does not exist"); 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) { xhivar = input->variable->find(xhistr); - if (xhivar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (xhivar < 0) error->all(FLERR, "Variable name for region block does not exist"); 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) { ylovar = input->variable->find(ylostr); - if (ylovar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (ylovar < 0) error->all(FLERR, "Variable name for region block does not exist"); 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) { yhivar = input->variable->find(yhistr); - if (yhivar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (yhivar < 0) error->all(FLERR, "Variable name for region block does not exist"); 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) { zlovar = input->variable->find(zlostr); - if (zlovar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (zlovar < 0) error->all(FLERR, "Variable name for region block does not exist"); 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) { zhivar = input->variable->find(zhistr); - if (zhivar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (zhivar < 0) error->all(FLERR, "Variable name for region block does not exist"); 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 --------------------------------------------------------------------------*/ -double RegBlock::find_closest_point(int i, double *x, - double &xc, double &yc, double &zc) +double RegBlock::find_closest_point(int i, double *x, double &xc, double &yc, double &zc) { - double dot,d2,d2min; - double xr[3],xproj[3],p[3]; + double dot, d2, d2min; + double xr[3], xproj[3], p[3]; xr[0] = x[0] - corners[i][0][0]; xr[1] = x[1] - corners[i][0][1]; xr[2] = x[2] - corners[i][0][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[1] = xr[1] - dot*face[i][1]; - xproj[2] = xr[2] - dot*face[i][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[1] = xr[1] - dot * face[i][1]; + xproj[2] = xr[2] - dot * face[i][2]; d2min = BIG; // check if point projects inside of face if (inside_face(xproj, i)) { - d2 = d2min = dot*dot; + d2 = d2min = dot * dot; xc = xproj[0] + corners[i][0][0]; yc = xproj[1] + corners[i][0][1]; zc = xproj[2] + corners[i][0][2]; - // check each edge + // check each edge } else { - 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]) + - (p[2]-x[2])*(p[2]-x[2]); + 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]) + + (p[2] - x[2]) * (p[2] - x[2]); if (d2 < d2min) { d2min = d2; xc = p[0]; @@ -582,9 +597,9 @@ double RegBlock::find_closest_point(int i, double *x, zc = p[2]; } - 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]) + - (p[2]-x[2])*(p[2]-x[2]); + 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]) + + (p[2] - x[2]) * (p[2] - x[2]); if (d2 < d2min) { d2min = d2; xc = p[0]; @@ -592,9 +607,9 @@ double RegBlock::find_closest_point(int i, double *x, zc = p[2]; } - 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]) + - (p[2]-x[2])*(p[2]-x[2]); + 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]) + + (p[2] - x[2]) * (p[2] - x[2]); if (d2 < d2min) { d2min = d2; xc = p[0]; @@ -602,9 +617,9 @@ double RegBlock::find_closest_point(int i, double *x, zc = p[2]; } - 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]) + - (p[2]-x[2])*(p[2]-x[2]); + 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]) + + (p[2] - x[2]) * (p[2] - x[2]); if (d2 < d2min) { d2min = d2; xc = p[0]; @@ -623,14 +638,12 @@ double RegBlock::find_closest_point(int i, double *x, int RegBlock::inside_face(double *xproj, int iface) { if (iface < 2) { - if (xproj[1] > 0 && (xproj[1] < yhi-ylo) && - xproj[2] > 0 && (xproj[2] < zhi-zlo)) return 1; + if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1; } else if (iface < 4) { - if (xproj[0] > 0 && (xproj[0] < (xhi-xlo)) && - xproj[2] > 0 && (xproj[2] < (zhi-zlo))) return 1; + if (xproj[0] > 0 && (xproj[0] < (xhi - xlo)) && xproj[2] > 0 && (xproj[2] < (zhi - zlo))) + return 1; } else { - if (xproj[0] > 0 && xproj[0] < (xhi-xlo) && - xproj[1] > 0 && xproj[1] < (yhi-ylo)) return 1; + if (xproj[0] > 0 && xproj[0] < (xhi - xlo) && xproj[1] > 0 && xproj[1] < (yhi - ylo)) return 1; } return 0; diff --git a/src/region_cone.cpp b/src/region_cone.cpp index e71a6a72a3..76bedf4c92 100644 --- a/src/region_cone.cpp +++ b/src/region_cone.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -26,102 +25,120 @@ using namespace LAMMPS_NS; -#define BIG 1.0e20 +static constexpr double BIG = 1.0e20; /* ---------------------------------------------------------------------- */ -RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : - Region(lmp, narg, arg), lo(0.0), hi(0.0) +RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo(0.0), hi(0.0) { - options(narg-9,&arg[9]); + options(narg - 9, &arg[9]); // check open face settings if (openflag) - for (int i=3; i<6; i++) - if (open_faces[i]) error->all(FLERR,"Illegal region cone open face: {}", i+1); + for (int i = 3; i < 6; i++) + 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) - error->all(FLERR,"Illegal region cone axis: {}", arg[2]); + 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]); axis = arg[2][0]; if (axis == 'x') { - c1 = yscale*utils::numeric(FLERR,arg[3],false,lmp); - c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); - radiuslo = yscale*utils::numeric(FLERR,arg[5],false,lmp); - radiushi = yscale*utils::numeric(FLERR,arg[6],false,lmp); + c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp); + c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp); + radiuslo = yscale * utils::numeric(FLERR, arg[5], false, lmp); + radiushi = yscale * utils::numeric(FLERR, arg[6], false, lmp); } else if (axis == 'y') { - c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); - c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); - radiuslo = xscale*utils::numeric(FLERR,arg[5],false,lmp); - radiushi = xscale*utils::numeric(FLERR,arg[6],false,lmp); + c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp); + c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp); + radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp); + radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp); } else if (axis == 'z') { - c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); - c2 = yscale*utils::numeric(FLERR,arg[4],false,lmp); - radiuslo = xscale*utils::numeric(FLERR,arg[5],false,lmp); - radiushi = xscale*utils::numeric(FLERR,arg[6],false,lmp); + c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp); + c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp); + radiuslo = xscale * utils::numeric(FLERR, arg[5], 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) - 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 (strcmp(arg[7],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[0]; - else lo = domain->boxlo_bound[0]; + if (strcmp(arg[7], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[0]; + else + lo = domain->boxlo_bound[0]; } if (axis == 'y') { - if (strcmp(arg[7],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[1]; - else lo = domain->boxlo_bound[1]; + if (strcmp(arg[7], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[1]; + else + lo = domain->boxlo_bound[1]; } if (axis == 'z') { - if (strcmp(arg[7],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[2]; - else lo = domain->boxlo_bound[2]; + if (strcmp(arg[7], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[2]; + else + lo = domain->boxlo_bound[2]; } } else { - 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 == 'z') lo = zscale*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 == '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) - 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 (strcmp(arg[8],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[0]; - else hi = domain->boxhi_bound[0]; + if (strcmp(arg[8], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[0]; + else + hi = domain->boxhi_bound[0]; } if (axis == 'y') { - if (strcmp(arg[8],"INF") == 0) hi = BIG; - if (domain->triclinic == 0) hi = domain->boxhi[1]; - else hi = domain->boxhi_bound[1]; + if (strcmp(arg[8], "INF") == 0) hi = BIG; + if (domain->triclinic == 0) + hi = domain->boxhi[1]; + else + hi = domain->boxhi_bound[1]; } if (axis == 'z') { - if (strcmp(arg[8],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[2]; - else hi = domain->boxhi_bound[2]; + if (strcmp(arg[8], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[2]; + else + hi = domain->boxhi_bound[2]; } } else { - 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 == 'z') hi = zscale*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 == 'z') hi = zscale * utils::numeric(FLERR, arg[8], false, lmp); } // error check - 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 (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 (radiuslo == 0.0 && radiushi == 0.0) - error->all(FLERR,"Illegal radius in region cone command"); - if (hi == lo) error->all(FLERR,"Illegal cone length 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"); // extent of cone - if (radiuslo > radiushi) maxradius = radiuslo; - else maxradius = radiushi; + if (radiuslo > radiushi) + maxradius = radiuslo; + else + maxradius = radiushi; if (interior) { bboxflag = 1; @@ -149,22 +166,25 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : extent_zlo = lo; extent_zhi = hi; } - } else bboxflag = 0; + } else + bboxflag = 0; // particle could be close to cone surface and 2 ends // particle can only touch surface and 1 end cmax = 3; contact = new Contact[cmax]; - if (interior) tmax = 2; - else tmax = 1; + if (interior) + tmax = 2; + else + tmax = 1; } /* ---------------------------------------------------------------------- */ RegCone::~RegCone() { - delete [] contact; + delete[] contact; } /* ---------------------------------------------------------------------- @@ -174,30 +194,36 @@ RegCone::~RegCone() int RegCone::inside(double x, double y, double z) { - double del1,del2,dist; + double del1, del2, dist; double currentradius; if (axis == 'x') { del1 = y - c1; del2 = z - c2; - dist = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x-lo)*(radiushi-radiuslo)/(hi-lo); - if (dist <= currentradius && x >= lo && x <= hi) return 1; - else return 0; + dist = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x - lo) * (radiushi - radiuslo) / (hi - lo); + if (dist <= currentradius && x >= lo && x <= hi) + return 1; + else + return 0; } else if (axis == 'y') { del1 = x - c1; del2 = z - c2; - dist = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (y-lo)*(radiushi-radiuslo)/(hi-lo); - if (dist <= currentradius && y >= lo && y <= hi) return 1; - else return 0; + dist = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (y - lo) * (radiushi - radiuslo) / (hi - lo); + if (dist <= currentradius && y >= lo && y <= hi) + return 1; + else + return 0; } else if (axis == 'z') { del1 = x - c1; del2 = y - c2; - dist = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (z-lo)*(radiushi-radiuslo)/(hi-lo); - if (dist <= currentradius && z >= lo && z <= hi) return 1; - else return 0; + dist = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (z - lo) * (radiushi - radiuslo) / (hi - lo); + if (dist <= currentradius && z >= lo && z <= hi) + return 1; + else + 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) { - double del1,del2,r,currentradius,delx,dely,delz,dist,delta; - double surflo[3],surfhi[3],xs[3]; + double del1, del2, r, currentradius, delx, dely, delz, dist, delta; + double surflo[3], surfhi[3], xs[3]; int n = 0; if (axis == 'x') { del1 = x[1] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[0] - lo) * (radiushi - radiuslo) / (hi - lo); // x is exterior to cone @@ -234,23 +260,22 @@ int RegCone::surface_interior(double *x, double cutoff) if (r > 0.0 && !open_faces[2]) { surflo[0] = lo; - surflo[1] = c1 + del1*radiuslo/r; - surflo[2] = c2 + del2*radiuslo/r; + surflo[1] = c1 + del1 * radiuslo / r; + surflo[2] = c2 + del2 * radiuslo / r; surfhi[0] = hi; - surfhi[1] = c1 + del1*radiushi/r; - surfhi[2] = c2 + del2*radiushi/r; - point_on_line_segment(surflo,surfhi,x,xs); + surfhi[1] = c1 + del1 * radiushi / r; + surfhi[2] = c2 + del2 * radiushi / r; + point_on_line_segment(surflo, surfhi, x, xs); delx = x[0] - xs[0]; dely = x[1] - xs[1]; delz = x[2] - xs[2]; - dist = sqrt(delx*delx + dely*dely + delz*delz); + dist = sqrt(delx * delx + dely * dely + delz * delz); if (dist < cutoff) { contact[n].r = dist; contact[n].delx = delx; contact[n].dely = dely; contact[n].delz = delz; - contact[n].radius = -2.0*(radiuslo + (xs[0]-lo)* - (radiushi-radiuslo)/(hi-lo)); + contact[n].radius = -2.0 * (radiuslo + (xs[0] - lo) * (radiushi - radiuslo) / (hi - lo)); contact[n].iwall = 2; n++; } @@ -278,9 +303,8 @@ int RegCone::surface_interior(double *x, double cutoff) } else if (axis == 'y') { del1 = x[0] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[1]-lo)* - (radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[1] - lo) * (radiushi - radiuslo) / (hi - lo); // 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 if (r > 0.0 && !open_faces[2]) { - surflo[0] = c1 + del1*radiuslo/r; + surflo[0] = c1 + del1 * radiuslo / r; surflo[1] = lo; - surflo[2] = c2 + del2*radiuslo/r; - surfhi[0] = c1 + del1*radiushi/r; + surflo[2] = c2 + del2 * radiuslo / r; + surfhi[0] = c1 + del1 * radiushi / r; surfhi[1] = hi; - surfhi[2] = c2 + del2*radiushi/r; - point_on_line_segment(surflo,surfhi,x,xs); + surfhi[2] = c2 + del2 * radiushi / r; + point_on_line_segment(surflo, surfhi, x, xs); delx = x[0] - xs[0]; dely = x[1] - xs[1]; delz = x[2] - xs[2]; - dist = sqrt(delx*delx + dely*dely + delz*delz); + dist = sqrt(delx * delx + dely * dely + delz * delz); if (dist < cutoff) { contact[n].r = dist; contact[n].delx = delx; contact[n].dely = dely; contact[n].delz = delz; contact[n].iwall = 2; - contact[n].radius = -2.0*(radiuslo + (xs[1]-lo)* - (radiushi-radiuslo)/(hi-lo)); + contact[n].radius = -2.0 * (radiuslo + (xs[1] - lo) * (radiushi - radiuslo) / (hi - lo)); n++; } } @@ -336,8 +359,8 @@ int RegCone::surface_interior(double *x, double cutoff) } else { del1 = x[0] - c1; del2 = x[1] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[2] - lo) * (radiushi - radiuslo) / (hi - lo); // 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 if (r > 0.0 && !open_faces[2]) { - surflo[0] = c1 + del1*radiuslo/r; - surflo[1] = c2 + del2*radiuslo/r; + surflo[0] = c1 + del1 * radiuslo / r; + surflo[1] = c2 + del2 * radiuslo / r; surflo[2] = lo; - surfhi[0] = c1 + del1*radiushi/r; - surfhi[1] = c2 + del2*radiushi/r; + surfhi[0] = c1 + del1 * radiushi / r; + surfhi[1] = c2 + del2 * radiushi / r; surfhi[2] = hi; - point_on_line_segment(surflo,surfhi,x,xs); + point_on_line_segment(surflo, surfhi, x, xs); delx = x[0] - xs[0]; dely = x[1] - xs[1]; delz = x[2] - xs[2]; - dist = sqrt(delx*delx + dely*dely + delz*delz); + dist = sqrt(delx * delx + dely * dely + delz * delz); if (dist < cutoff) { contact[n].r = dist; contact[n].delx = delx; contact[n].dely = dely; contact[n].delz = delz; contact[n].iwall = 2; - contact[n].radius = -2.0*(radiuslo + (xs[2]-lo)* - (radiushi-radiuslo)/(hi-lo)); + contact[n].radius = -2.0 * (radiuslo + (xs[2] - lo) * (radiushi - radiuslo) / (hi - lo)); n++; } } delta = x[2] - lo; - if (delta < cutoff && !open_faces[0]) { + if (delta < cutoff && !open_faces[0]) { contact[n].r = delta; contact[n].delz = delta; contact[n].delx = contact[n].dely = 0.0; @@ -381,7 +403,7 @@ int RegCone::surface_interior(double *x, double cutoff) n++; } delta = hi - x[2]; - if (delta < cutoff && !open_faces[1]) { + if (delta < cutoff && !open_faces[1]) { contact[n].r = delta; contact[n].delz = -delta; 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) { - double del1,del2,r,currentradius,distsq,distsqprev,crad; - double corner1[3],corner2[3],corner3[3],corner4[3],xp[3],nearest[3]; + double del1, del2, r, currentradius, distsq, distsqprev, crad; + double corner1[3], corner2[3], corner3[3], corner4[3], xp[3], nearest[3]; if (axis == 'x') { del1 = x[1] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[0] - lo) * (radiushi - radiuslo) / (hi - lo); // 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 interior to cone - if (r >= maxradius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff) - return 0; + if (r >= maxradius + cutoff || x[0] <= lo - cutoff || x[0] >= hi + cutoff) return 0; if (r < currentradius && x[0] > lo && x[0] < hi) return 0; // 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 corner1[0] = lo; - corner1[1] = c1 + del1*radiuslo/r; - corner1[2] = c2 + del2*radiuslo/r; + corner1[1] = c1 + del1 * radiuslo / r; + corner1[2] = c2 + del2 * radiuslo / r; corner2[0] = hi; - corner2[1] = c1 + del1*radiushi/r; - corner2[2] = c2 + del2*radiushi/r; + corner2[1] = c1 + del1 * radiushi / r; + corner2[2] = c2 + del2 * radiushi / r; corner3[0] = lo; corner3[1] = c1; corner3[2] = c2; @@ -445,27 +466,27 @@ int RegCone::surface_exterior(double *x, double cutoff) distsq = BIG; if (!open_faces[2]) { - point_on_line_segment(corner1,corner2,x,xp); - distsq = closest(x,xp,nearest,distsq); - crad = -2.0*(radiuslo + (nearest[0]-lo)*(radiushi-radiuslo)/(hi-lo)); + point_on_line_segment(corner1, corner2, x, xp); + distsq = closest(x, xp, nearest, distsq); + crad = -2.0 * (radiuslo + (nearest[0] - lo) * (radiushi - radiuslo) / (hi - lo)); } if (!open_faces[0]) { - point_on_line_segment(corner1,corner3,x,xp); + point_on_line_segment(corner1, corner3, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0.0; } if (!open_faces[1]) { - point_on_line_segment(corner2,corner4,x,xp); + point_on_line_segment(corner2, corner4, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0.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].iwall = 0; if (contact[0].r < cutoff) return 1; @@ -474,8 +495,8 @@ int RegCone::surface_exterior(double *x, double cutoff) } else if (axis == 'y') { del1 = x[0] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[1]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[1] - lo) * (radiushi - radiuslo) / (hi - lo); // 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 interior to cone - if (r >= maxradius+cutoff || - x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0; + if (r >= maxradius + cutoff || x[1] <= lo - cutoff || x[1] >= hi + cutoff) return 0; if (r < currentradius && x[1] > lo && x[1] < hi) return 0; // 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 // do not add contact point if r >= cutoff - corner1[0] = c1 + del1*radiuslo/r; + corner1[0] = c1 + del1 * radiuslo / r; corner1[1] = lo; - corner1[2] = c2 + del2*radiuslo/r; - corner2[0] = c1 + del1*radiushi/r; + corner1[2] = c2 + del2 * radiuslo / r; + corner2[0] = c1 + del1 * radiushi / r; corner2[1] = hi; - corner2[2] = c2 + del2*radiushi/r; + corner2[2] = c2 + del2 * radiushi / r; corner3[0] = c1; corner3[1] = lo; corner3[2] = c2; @@ -511,26 +531,26 @@ int RegCone::surface_exterior(double *x, double cutoff) distsq = BIG; if (!open_faces[2]) { - point_on_line_segment(corner1,corner2,x,xp); - distsq = closest(x,xp,nearest,distsq); - crad = -2.0*(radiuslo + (nearest[1]-lo)*(radiushi-radiuslo)/(hi-lo)); + point_on_line_segment(corner1, corner2, x, xp); + distsq = closest(x, xp, nearest, distsq); + crad = -2.0 * (radiuslo + (nearest[1] - lo) * (radiushi - radiuslo) / (hi - lo)); } if (!open_faces[0]) { - point_on_line_segment(corner1,corner3,x,xp); + point_on_line_segment(corner1, corner3, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0; } if (!open_faces[1]) { - point_on_line_segment(corner2,corner4,x,xp); + point_on_line_segment(corner2, corner4, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); 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].iwall = 0; if (contact[0].r < cutoff) return 1; @@ -539,8 +559,8 @@ int RegCone::surface_exterior(double *x, double cutoff) } else { del1 = x[0] - c1; del2 = x[1] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[2] - lo) * (radiushi - radiuslo) / (hi - lo); // 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 interior to cone - if (r >= maxradius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff) - return 0; + if (r >= maxradius + cutoff || x[2] <= lo - cutoff || x[2] >= hi + cutoff) return 0; if (r < currentradius && x[2] > lo && x[2] < hi) return 0; // 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 // do not add contact point if r >= cutoff - corner1[0] = c1 + del1*radiuslo/r; - corner1[1] = c2 + del2*radiuslo/r; + corner1[0] = c1 + del1 * radiuslo / r; + corner1[1] = c2 + del2 * radiuslo / r; corner1[2] = lo; - corner2[0] = c1 + del1*radiushi/r; - corner2[1] = c2 + del2*radiushi/r; + corner2[0] = c1 + del1 * radiushi / r; + corner2[1] = c2 + del2 * radiushi / r; corner2[2] = hi; corner3[0] = c1; corner3[1] = c2; @@ -576,26 +595,26 @@ int RegCone::surface_exterior(double *x, double cutoff) distsq = BIG; if (!open_faces[2]) { - point_on_line_segment(corner1,corner2,x,xp); - distsq = closest(x,xp,nearest,distsq); - crad = -2.0*(radiuslo + (nearest[2]-lo)*(radiushi-radiuslo)/(hi-lo)); + point_on_line_segment(corner1, corner2, x, xp); + distsq = closest(x, xp, nearest, distsq); + crad = -2.0 * (radiuslo + (nearest[2] - lo) * (radiushi - radiuslo) / (hi - lo)); } if (!open_faces[0]) { - point_on_line_segment(corner1,corner3,x,xp); + point_on_line_segment(corner1, corner3, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0; } if (!open_faces[1]) { - point_on_line_segment(corner2,corner4,x,xp); + point_on_line_segment(corner2, corner4, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); 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].iwall = 0; 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 dely = x[1] - near[1]; 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; nearest[0] = near[0]; diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp index 259b855a66..c6119c0c3b 100644 --- a/src/region_cylinder.cpp +++ b/src/region_cylinder.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -25,94 +24,97 @@ using namespace LAMMPS_NS; -#define BIG 1.0e20 -enum{CONSTANT,VARIABLE}; +static constexpr double BIG = 1.0e20; + +enum { CONSTANT, VARIABLE }; /* ---------------------------------------------------------------------- */ 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 if (openflag) - for (int i=3; i<6; i++) - if (open_faces[i]) error->all(FLERR,"Illegal region cylinder open face: {}", i+1); + for (int i = 3; i < 6; i++) + 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) - error->all(FLERR,"Illegal region cylinder axis: {}", arg[2]); + 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]); axis = arg[2][0]; if (axis == 'x') { - if (utils::strmatch(arg[3],"^v_")) { - c1str = utils::strdup(arg[3]+2); + if (utils::strmatch(arg[3], "^v_")) { + c1str = utils::strdup(arg[3] + 2); c1 = 0.0; c1style = VARIABLE; varshape = 1; } else { - c1 = yscale*utils::numeric(FLERR,arg[3],false,lmp); + c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp); c1style = CONSTANT; } - if (utils::strmatch(arg[4],"^v_")) { - c2str = utils::strdup(arg[4]+2); + if (utils::strmatch(arg[4], "^v_")) { + c2str = utils::strdup(arg[4] + 2); c2 = 0.0; c2style = VARIABLE; varshape = 1; } else { - c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); + c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp); c2style = CONSTANT; } } else if (axis == 'y') { - if (utils::strmatch(arg[3],"^v_")) { - c1str = utils::strdup(arg[3]+2); + if (utils::strmatch(arg[3], "^v_")) { + c1str = utils::strdup(arg[3] + 2); c1 = 0.0; c1style = VARIABLE; varshape = 1; } else { - c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); + c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp); c1style = CONSTANT; } - if (utils::strmatch(arg[4],"^v_")) { - c2str = utils::strdup(arg[4]+2); + if (utils::strmatch(arg[4], "^v_")) { + c2str = utils::strdup(arg[4] + 2); c2 = 0.0; c2style = VARIABLE; varshape = 1; } else { - c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); + c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp); c2style = CONSTANT; } } else if (axis == 'z') { - if (utils::strmatch(arg[3],"^v_")) { - c1str = utils::strdup(arg[3]+2); + if (utils::strmatch(arg[3], "^v_")) { + c1str = utils::strdup(arg[3] + 2); c1 = 0.0; c1style = VARIABLE; varshape = 1; } else { - c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); + c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp); c1style = CONSTANT; } - if (utils::strmatch(arg[4],"^v_")) { - c2str = utils::strdup(arg[4]+2); + if (utils::strmatch(arg[4], "^v_")) { + c2str = utils::strdup(arg[4] + 2); c2 = 0.0; c2style = VARIABLE; varshape = 1; } else { - c2 = yscale*utils::numeric(FLERR,arg[4],false,lmp); + c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp); c2style = CONSTANT; } } - if (utils::strmatch(arg[5],"^v_")) { - rstr = utils::strdup(arg[5]+2); + if (utils::strmatch(arg[5], "^v_")) { + rstr = utils::strdup(arg[5] + 2); radius = 0.0; rstyle = VARIABLE; varshape = 1; } else { - radius = utils::numeric(FLERR,arg[5],false,lmp); - if (axis == 'x') radius *= yscale; - else radius *= xscale; + radius = utils::numeric(FLERR, arg[5], false, lmp); + if (axis == 'x') + radius *= yscale; + else + radius *= xscale; rstyle = CONSTANT; } @@ -121,57 +123,75 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) : 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) - 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 (strcmp(arg[6],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[0]; - else lo = domain->boxlo_bound[0]; + if (strcmp(arg[6], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[0]; + else + lo = domain->boxlo_bound[0]; } if (axis == 'y') { - if (strcmp(arg[6],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[1]; - else lo = domain->boxlo_bound[1]; + if (strcmp(arg[6], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[1]; + else + lo = domain->boxlo_bound[1]; } if (axis == 'z') { - if (strcmp(arg[6],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[2]; - else lo = domain->boxlo_bound[2]; + if (strcmp(arg[6], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[2]; + else + lo = domain->boxlo_bound[2]; } } else { - 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 == 'z') lo = zscale*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 == '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) - 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 (strcmp(arg[7],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[0]; - else hi = domain->boxhi_bound[0]; + if (strcmp(arg[7], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[0]; + else + hi = domain->boxhi_bound[0]; } if (axis == 'y') { - if (strcmp(arg[7],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[1]; - else hi = domain->boxhi_bound[1]; + if (strcmp(arg[7], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[1]; + else + hi = domain->boxhi_bound[1]; } if (axis == 'z') { - if (strcmp(arg[7],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[2]; - else hi = domain->boxhi_bound[2]; + if (strcmp(arg[7], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[2]; + else + hi = domain->boxhi_bound[2]; } } else { - 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 == 'z') hi = zscale*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 == 'z') hi = zscale * utils::numeric(FLERR, arg[7], false, lmp); } // 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 // for variable radius, uses initial radius @@ -202,25 +222,28 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) : extent_zlo = lo; extent_zhi = hi; } - } else bboxflag = 0; + } else + bboxflag = 0; // particle could be close to cylinder surface and 2 ends // particle can only touch surface and 1 end cmax = 3; contact = new Contact[cmax]; - if (interior) tmax = 2; - else tmax = 1; + if (interior) + tmax = 2; + else + tmax = 1; } /* ---------------------------------------------------------------------- */ RegCylinder::~RegCylinder() { - delete [] c1str; - delete [] c2str; - delete [] rstr; - delete [] contact; + delete[] c1str; + delete[] c2str; + delete[] rstr; + delete[] contact; } /* ---------------------------------------------------------------------- */ @@ -238,27 +261,33 @@ void RegCylinder::init() int RegCylinder::inside(double x, double y, double z) { - double del1,del2,dist; + double del1, del2, dist; int inside; if (axis == 'x') { del1 = y - c1; del2 = z - c2; - dist = sqrt(del1*del1 + del2*del2); - if (dist <= radius && x >= lo && x <= hi) inside = 1; - else inside = 0; + dist = sqrt(del1 * del1 + del2 * del2); + if (dist <= radius && x >= lo && x <= hi) + inside = 1; + else + inside = 0; } else if (axis == 'y') { del1 = x - c1; del2 = z - c2; - dist = sqrt(del1*del1 + del2*del2); - if (dist <= radius && y >= lo && y <= hi) inside = 1; - else inside = 0; + dist = sqrt(del1 * del1 + del2 * del2); + if (dist <= radius && y >= lo && y <= hi) + inside = 1; + else + inside = 0; } else { del1 = x - c1; del2 = y - c2; - dist = sqrt(del1*del1 + del2*del2); - if (dist <= radius && z >= lo && z <= hi) inside = 1; - else inside = 0; + dist = sqrt(del1 * del1 + del2 * del2); + if (dist <= radius && z >= lo && z <= hi) + inside = 1; + else + inside = 0; } return inside; @@ -274,14 +303,14 @@ int RegCylinder::inside(double x, double y, double z) int RegCylinder::surface_interior(double *x, double cutoff) { - double del1,del2,r,delta; + double del1, del2, r, delta; int n = 0; if (axis == 'x') { del1 = x[1] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); + r = sqrt(del1 * del1 + del2 * del2); // 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]) { contact[n].r = delta; contact[n].delx = 0.0; - contact[n].dely = del1*(1.0-radius/r); - contact[n].delz = del2*(1.0-radius/r); - contact[n].radius = -2.0*radius; + contact[n].dely = del1 * (1.0 - radius / r); + contact[n].delz = del2 * (1.0 - radius / r); + contact[n].radius = -2.0 * radius; contact[n].iwall = 2; contact[n].varflag = 1; n++; @@ -324,7 +353,7 @@ int RegCylinder::surface_interior(double *x, double cutoff) } else if (axis == 'y') { del1 = x[0] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); + r = sqrt(del1 * del1 + del2 * del2); // y is exterior to cylinder @@ -335,10 +364,10 @@ int RegCylinder::surface_interior(double *x, double cutoff) delta = radius - r; if (delta < cutoff && r > 0.0 && !open_faces[2]) { 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].delz = del2*(1.0-radius/r); - contact[n].radius = -2.0*radius; + contact[n].delz = del2 * (1.0 - radius / r); + contact[n].radius = -2.0 * radius; contact[n].iwall = 2; contact[n].varflag = 1; n++; @@ -367,7 +396,7 @@ int RegCylinder::surface_interior(double *x, double cutoff) } else { del1 = x[0] - c1; del2 = x[1] - c2; - r = sqrt(del1*del1 + del2*del2); + r = sqrt(del1 * del1 + del2 * del2); // z is exterior to cylinder @@ -378,10 +407,10 @@ int RegCylinder::surface_interior(double *x, double cutoff) delta = radius - r; if (delta < cutoff && r > 0.0 && !open_faces[2]) { contact[n].r = delta; - contact[n].delx = del1*(1.0-radius/r); - contact[n].dely = del2*(1.0-radius/r); + contact[n].delx = del1 * (1.0 - radius / r); + contact[n].dely = del2 * (1.0 - radius / r); contact[n].delz = 0.0; - contact[n].radius = -2.0*radius; + contact[n].radius = -2.0 * radius; contact[n].iwall = 2; contact[n].varflag = 1; n++; @@ -419,12 +448,12 @@ int RegCylinder::surface_interior(double *x, double cutoff) int RegCylinder::surface_exterior(double *x, double cutoff) { - double del1,del2,r; - double xp,yp,zp; + double del1, del2, r; + double xp, yp, zp; double dx, dr, dr2, d2, d2prev; - // radius of curvature for granular - // 0 for flat surfaces (infinite case), 2*radius for curved portion + // radius of curvature for granular + // 0 for flat surfaces (infinite case), 2*radius for curved portion double crad = 0.0; int varflag = 0; @@ -432,12 +461,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (axis == 'x') { del1 = x[1] - c1; 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 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; // x is exterior to cylinder or on its surface @@ -448,48 +477,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff) d2prev = BIG; if (!openflag) { if (r > radius) { - yp = c1 + del1*radius/r; - zp = c2 + del2*radius/r; - crad = 2.0*radius; + yp = c1 + del1 * radius / r; + zp = c2 + del2 * radius / r; + crad = 2.0 * radius; varflag = 1; } else { yp = x[1]; zp = x[2]; } - if (x[0] < lo) xp = lo; - else if (x[0] > hi) xp = hi; - else xp = x[0]; + if (x[0] < lo) + xp = lo; + else if (x[0] > hi) + xp = hi; + else + xp = x[0]; } else { - // closest point on curved surface + // closest point on curved surface dr = r - radius; - dr2 = dr*dr; + dr2 = dr * dr; if (!open_faces[2]) { - yp = c1 + del1*radius/r; - zp = c2 + del2*radius/r; + yp = c1 + del1 * radius / r; + zp = c2 + del2 * radius / r; if (x[0] < lo) { - dx = lo-x[0]; + dx = lo - x[0]; xp = lo; - } - else if (x[0] > hi) { - dx = x[0]-hi; + } else if (x[0] > hi) { + dx = x[0] - hi; xp = hi; - } - else { + } else { dx = 0; xp = x[0]; } - d2 = d2prev = dr2 + dx*dx; + d2 = d2prev = dr2 + dx * dx; } // closest point on bottom cap if (!open_faces[0]) { dx = lo - x[0]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { xp = lo; if (r < radius) { @@ -504,8 +536,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (!open_faces[1]) { dx = hi - x[0]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { xp = hi; 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].varflag = varflag; contact[0].iwall = 0; @@ -526,12 +560,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff) } else if (axis == 'y') { del1 = x[0] - c1; 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 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; // y is exterior to cylinder or on its surface @@ -542,48 +576,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff) d2prev = BIG; if (!openflag) { if (r > radius) { - xp = c1 + del1*radius/r; - zp = c2 + del2*radius/r; - crad = 2.0*radius; + xp = c1 + del1 * radius / r; + zp = c2 + del2 * radius / r; + crad = 2.0 * radius; varflag = 1; } else { xp = x[0]; zp = x[2]; } - if (x[1] < lo) yp = lo; - else if (x[1] > hi) yp = hi; - else yp = x[1]; + if (x[1] < lo) + yp = lo; + else if (x[1] > hi) + yp = hi; + else + yp = x[1]; } else { - // closest point on curved surface + // closest point on curved surface dr = r - radius; - dr2 = dr*dr; + dr2 = dr * dr; if (!open_faces[2]) { - xp = c1 + del1*radius/r; - zp = c2 + del2*radius/r; + xp = c1 + del1 * radius / r; + zp = c2 + del2 * radius / r; if (x[1] < lo) { - dx = lo-x[1]; + dx = lo - x[1]; yp = lo; - } - else if (x[1] > hi) { - dx = x[1]-hi; + } else if (x[1] > hi) { + dx = x[1] - hi; yp = hi; - } - else { + } else { dx = 0; yp = x[1]; } - d2 = d2prev = dr2 + dx*dx; + d2 = d2prev = dr2 + dx * dx; } // closest point on bottom cap if (!open_faces[0]) { dx = lo - x[1]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { yp = lo; if (r < radius) { @@ -598,8 +635,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (!open_faces[1]) { dx = hi - x[1]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { yp = hi; 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].varflag = varflag; contact[0].iwall = 0; @@ -620,12 +659,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff) } else { del1 = x[0] - c1; 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 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; // z is exterior to cylinder or on its surface @@ -636,46 +675,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff) d2prev = BIG; if (!openflag) { if (r > radius) { - xp = c1 + del1*radius/r; - yp = c2 + del2*radius/r; - crad = 2.0*radius; + xp = c1 + del1 * radius / r; + yp = c2 + del2 * radius / r; + crad = 2.0 * radius; varflag = 1; } else { xp = x[0]; yp = x[1]; } - if (x[2] < lo) zp = lo; - else if (x[2] > hi) zp = hi; - else zp = x[2]; + if (x[2] < lo) + zp = lo; + else if (x[2] > hi) + zp = hi; + else + zp = x[2]; } else { - // closest point on curved surface + // closest point on curved surface dr = r - radius; - dr2 = dr*dr; + dr2 = dr * dr; if (!open_faces[2]) { - xp = c1 + del1*radius/r; - yp = c2 + del2*radius/r; + xp = c1 + del1 * radius / r; + yp = c2 + del2 * radius / r; if (x[2] < lo) { - dx = lo-x[2]; + dx = lo - x[2]; zp = lo; } else if (x[2] > hi) { - dx = x[2]-hi; + dx = x[2] - hi; zp = hi; } else { dx = 0; zp = x[2]; } - d2prev = dr2 + dx*dx; + d2prev = dr2 + dx * dx; } // closest point on bottom cap if (!open_faces[0]) { dx = lo - x[2]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { zp = lo; if (r < radius) { @@ -690,8 +734,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (!open_faces[1]) { dx = hi - x[2]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { zp = hi; 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].varflag = varflag; contact[0].iwall = 0; @@ -721,22 +767,21 @@ void RegCylinder::shape_update() if (c2style == VARIABLE) c2 = input->variable->compute_equal(c2var); if (rstyle == VARIABLE) { radius = input->variable->compute_equal(rvar); - if (radius < 0.0) - error->one(FLERR,"Variable evaluation in region gave bad value"); + if (radius < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value"); } if (axis == 'x') { if (c1style == VARIABLE) c1 *= yscale; if (c2style == VARIABLE) c2 *= zscale; - if (rstyle == VARIABLE) radius *= yscale; + if (rstyle == VARIABLE) radius *= yscale; } else if (axis == 'y') { if (c1style == VARIABLE) c1 *= xscale; if (c2style == VARIABLE) c2 *= zscale; - if (rstyle == VARIABLE) radius *= xscale; - } else { // axis == 'z' + if (rstyle == VARIABLE) radius *= xscale; + } else { // axis == 'z' if (c1style == VARIABLE) c1 *= xscale; 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) { c1var = input->variable->find(c1str); - if (c1var < 0) - error->all(FLERR,"Variable name for region cylinder does not exist"); + if (c1var < 0) error->all(FLERR, "Variable name for region cylinder does not exist"); 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) { c2var = input->variable->find(c2str); - if (c2var < 0) - error->all(FLERR,"Variable name for region cylinder does not exist"); + if (c2var < 0) error->all(FLERR, "Variable name for region cylinder does not exist"); 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) { rvar = input->variable->find(rstr); - if (rvar < 0) - error->all(FLERR,"Variable name for region cylinder does not exist"); + if (rvar < 0) error->all(FLERR, "Variable name for region cylinder does not exist"); 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. These values do not depend on the contact, so this function is @@ -795,35 +836,34 @@ void RegCylinder::set_velocity_shape() xcenter[2] = 0; } forward_transform(xcenter[0], xcenter[1], xcenter[2]); - if (update->ntimestep > 0) rprev = prev[4]; - else rprev = radius; + if (update->ntimestep > 0) + rprev = prev[4]; + else + rprev = radius; prev[4] = radius; } - - /* ---------------------------------------------------------------------- add velocity due to shape change to wall velocity ------------------------------------------------------------------------- */ 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') { delx = 0; - dely = (xc[1] - xcenter[1])*(1 - rprev/radius); - delz = (xc[2] - xcenter[2])*(1 - rprev/radius); + dely = (xc[1] - xcenter[1]) * (1 - rprev / radius); + delz = (xc[2] - xcenter[2]) * (1 - rprev / radius); } else if (axis == 'y') { - delx = (xc[0] - xcenter[0])*(1 - rprev/radius); + delx = (xc[0] - xcenter[0]) * (1 - rprev / radius); dely = 0; - delz = (xc[2] - xcenter[2])*(1 - rprev/radius); + delz = (xc[2] - xcenter[2]) * (1 - rprev / radius); } else { - delx = (xc[0] - xcenter[0])*(1 - rprev/radius); - dely = (xc[1] - xcenter[1])*(1 - rprev/radius); + delx = (xc[0] - xcenter[0]) * (1 - rprev / radius); + dely = (xc[1] - xcenter[1]) * (1 - rprev / radius); delz = 0; } - vwall[0] += delx/update->dt; - vwall[1] += dely/update->dt; - vwall[2] += delz/update->dt; + vwall[0] += delx / update->dt; + vwall[1] += dely / update->dt; + vwall[2] += delz / update->dt; } - diff --git a/src/region_prism.cpp b/src/region_prism.cpp index 7c06c0aa7d..da2048cce0 100644 --- a/src/region_prism.cpp +++ b/src/region_prism.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -26,105 +25,126 @@ 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) { - 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[2],"INF") == 0) xlo = -BIG; - else xlo = domain->boxlo[0]; - } else xlo = xscale*utils::numeric(FLERR,arg[2],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[2], "INF") == 0) + xlo = -BIG; + 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[3],"INF") == 0) xhi = BIG; - else xhi = domain->boxhi[0]; - } else xhi = xscale*utils::numeric(FLERR,arg[3],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[3], "INF") == 0) + xhi = BIG; + 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[4],"INF") == 0) ylo = -BIG; - else ylo = domain->boxlo[1]; - } else ylo = yscale*utils::numeric(FLERR,arg[4],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[4], "INF") == 0) + ylo = -BIG; + 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[5],"INF") == 0) yhi = BIG; - else yhi = domain->boxhi[1]; - } else yhi = yscale*utils::numeric(FLERR,arg[5],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[5], "INF") == 0) + yhi = BIG; + 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[6],"INF") == 0) zlo = -BIG; - else zlo = domain->boxlo[2]; - } else zlo = zscale*utils::numeric(FLERR,arg[6],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[6], "INF") == 0) + zlo = -BIG; + 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) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[7],"INF") == 0) zhi = BIG; - else zhi = domain->boxhi[2]; - } else zhi = zscale*utils::numeric(FLERR,arg[7],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[7], "INF") == 0) + zhi = BIG; + else + zhi = domain->boxhi[2]; + } else + zhi = zscale * utils::numeric(FLERR, arg[7], false, lmp); - xy = xscale*utils::numeric(FLERR,arg[8],false,lmp); - xz = xscale*utils::numeric(FLERR,arg[9],false,lmp); - yz = yscale*utils::numeric(FLERR,arg[10],false,lmp); + xy = xscale * utils::numeric(FLERR, arg[8], false, lmp); + xz = xscale * utils::numeric(FLERR, arg[9], false, lmp); + yz = yscale * utils::numeric(FLERR, arg[10], false, lmp); // error check // 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 - 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 (zlo >= zhi) error->all(FLERR,"Illegal region prism zlo: {} >= zhi: {}", zlo ,zhi); + 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 (zlo >= zhi) error->all(FLERR, "Illegal region prism zlo: {} >= zhi: {}", zlo, zhi); 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) - 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) - 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) - 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) - 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) - 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 if (interior) { bboxflag = 1; - extent_xlo = MIN(xlo,xlo+xy); - extent_xlo = MIN(extent_xlo,extent_xlo+xz); - extent_ylo = MIN(ylo,ylo+yz); + extent_xlo = MIN(xlo, xlo + xy); + extent_xlo = MIN(extent_xlo, extent_xlo + xz); + extent_ylo = MIN(ylo, ylo + yz); extent_zlo = zlo; - extent_xhi = MAX(xhi,xhi+xy); - extent_xhi = MAX(extent_xhi,extent_xhi+xz); - extent_yhi = MAX(yhi,yhi+yz); + extent_xhi = MAX(xhi, xhi + xy); + extent_xhi = MAX(extent_xhi, extent_xhi + xz); + extent_yhi = MAX(yhi, yhi + yz); extent_zhi = zhi; - } else bboxflag = 0; + } else + bboxflag = 0; // particle could be close to all 6 planes // particle can only touch 3 planes cmax = 6; contact = new Contact[cmax]; - if (interior) tmax = 3; - else tmax = 1; + if (interior) + tmax = 3; + else + tmax = 1; // h = transformation matrix from tilt coords (0-1) to box coords (xyz) // 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[2][2] = zhi - zlo; - hinv[0][0] = 1.0/h[0][0]; - 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[1][1] = 1.0/h[1][1]; - hinv[1][2] = -h[1][2] / (h[1][1]*h[2][2]); - hinv[2][2] = 1.0/h[2][2]; + hinv[0][0] = 1.0 / h[0][0]; + 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[1][1] = 1.0 / h[1][1]; + hinv[1][2] = -h[1][2] / (h[1][1] * h[2][2]); + hinv[2][2] = 1.0 / h[2][2]; // corners = 8 corner points of prism // order = x varies fastest, then y, finally z // clo/chi = lo and hi corner pts of prism - a[0] = xhi-xlo; + a[0] = xhi - xlo; a[1] = 0.0; a[2] = 0.0; b[0] = xy; - b[1] = yhi-ylo; + b[1] = yhi - ylo; b[2] = 0.0; c[0] = xz; c[1] = yz; - c[2] = zhi-zlo; + c[2] = zhi - zlo; clo[0] = corners[0][0] = xlo; 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 // order = xy plane, xz plane, yz plane - MathExtra::cross3(a,b,face[0]); - MathExtra::cross3(b,a,face[1]); - MathExtra::cross3(c,a,face[2]); - MathExtra::cross3(a,c,face[3]); - MathExtra::cross3(b,c,face[4]); - MathExtra::cross3(c,b,face[5]); + MathExtra::cross3(a, b, face[0]); + MathExtra::cross3(b, a, face[1]); + MathExtra::cross3(c, a, face[2]); + MathExtra::cross3(a, c, face[3]); + MathExtra::cross3(b, c, face[4]); + MathExtra::cross3(c, b, face[5]); // remap open face indices to be consistent if (openflag) { int temp[6]; - for (int i = 0; i < 6; i++) - temp[i] = open_faces[i]; + for (int i = 0; i < 6; i++) temp[i] = open_faces[i]; open_faces[0] = temp[4]; open_faces[1] = temp[5]; 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 // order = xy plane, xz plane, yz plane - tri[0][0] = 0; tri[0][1] = 1; tri[0][2] = 3; - 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[0][0] = 0; + tri[0][1] = 1; + tri[0][2] = 3; + 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[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[4][0] = 0; + tri[4][1] = 4; + tri[4][2] = 5; + 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[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; + tri[8][0] = 2; + tri[8][1] = 6; + tri[8][2] = 4; + 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() { - delete [] contact; + delete[] contact; } /* ---------------------------------------------------------------------- @@ -258,12 +301,11 @@ RegPrism::~RegPrism() 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 b = hinv[1][1]*(y-ylo) + hinv[1][2]*(z-zlo); - double c = hinv[2][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 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) - return 1; + if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0) return 1; return 0; } @@ -283,10 +325,12 @@ int RegPrism::surface_interior(double *x, double cutoff) // x is exterior to prism for (i = 0; i < 6; i++) { - if (i % 2) corner = chi; - else 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 (i % 2) + corner = chi; + else + 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; } @@ -296,15 +340,17 @@ int RegPrism::surface_interior(double *x, double cutoff) for (i = 0; i < 6; i++) { if (open_faces[i]) continue; - if (i % 2) corner = chi; - else 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 (i % 2) + corner = chi; + else + 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) { contact[n].r = dot; - contact[n].delx = dot*face[i][0]; - contact[n].dely = dot*face[i][1]; - contact[n].delz = dot*face[i][2]; + contact[n].delx = dot * face[i][0]; + contact[n].dely = dot * face[i][1]; + contact[n].delz = dot * face[i][2]; contact[n].radius = 0; contact[n].iwall = i; n++; @@ -325,25 +371,29 @@ int RegPrism::surface_exterior(double *x, double cutoff) int i; double dot; double *corner; - double xp,yp,zp; + double xp, yp, zp; // x is far enough from prism that there is no contact for (i = 0; i < 6; i++) { - if (i % 2) corner = chi; - else 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 (i % 2) + corner = chi; + else + 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; } // x is interior to prism for (i = 0; i < 6; i++) { - if (i % 2) corner = chi; - else 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 (i % 2) + corner = chi; + else + 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; } @@ -354,8 +404,8 @@ int RegPrism::surface_exterior(double *x, double cutoff) // could be edge or corner pt of prism // do not add contact point if r >= cutoff - find_nearest(x,xp,yp,zp); - add_contact(0,x,xp,yp,zp); + find_nearest(x, xp, yp, zp); + add_contact(0, x, xp, yp, zp); contact[0].radius = 0; contact[0].iwall = 0; 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) { - int i,j,k,iface; - double xproj[3],xline[3],nearest[3]; + int i, j, k, iface; + double xproj[3], xline[3], nearest[3]; double dot; // 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; for (int itri = 0; itri < 12; itri++) { - iface = itri/2; + iface = itri / 2; if (open_faces[iface]) continue; i = tri[itri][0]; j = tri[itri][1]; k = tri[itri][2]; - dot = (x[0]-corners[i][0])*face[iface][0] + - (x[1]-corners[i][1])*face[iface][1] + - (x[2]-corners[i][2])*face[iface][2]; - xproj[0] = x[0] - dot*face[iface][0]; - xproj[1] = x[1] - dot*face[iface][1]; - xproj[2] = x[2] - dot*face[iface][2]; - if (inside_tri(xproj,corners[i],corners[j],corners[k],face[iface])) { - distsq = closest(x,xproj,nearest,distsq); - } - else { - point_on_line_segment(corners[i],corners[j],xproj,xline); - distsq = closest(x,xline,nearest,distsq); - point_on_line_segment(corners[j],corners[k],xproj,xline); - distsq = closest(x,xline,nearest,distsq); - point_on_line_segment(corners[i],corners[k],xproj,xline); - distsq = closest(x,xline,nearest,distsq); + dot = (x[0] - corners[i][0]) * face[iface][0] + (x[1] - corners[i][1]) * face[iface][1] + + (x[2] - corners[i][2]) * face[iface][2]; + xproj[0] = x[0] - dot * face[iface][0]; + xproj[1] = x[1] - dot * face[iface][1]; + xproj[2] = x[2] - dot * face[iface][2]; + if (inside_tri(xproj, corners[i], corners[j], corners[k], face[iface])) { + distsq = closest(x, xproj, nearest, distsq); + } else { + point_on_line_segment(corners[i], corners[j], xproj, xline); + distsq = closest(x, xline, nearest, distsq); + point_on_line_segment(corners[j], corners[k], xproj, xline); + 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 ------------------------------------------------------------------------- */ -int RegPrism::inside_tri(double *x, double *v1, double *v2, double *v3, - double *norm) +int RegPrism::inside_tri(double *x, double *v1, double *v2, double *v3, double *norm) { - double edge[3],pvec[3],xproduct[3]; + double edge[3], pvec[3], xproduct[3]; - MathExtra::sub3(v2,v1,edge); - MathExtra::sub3(x,v1,pvec); - MathExtra::cross3(edge,pvec,xproduct); - if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; + MathExtra::sub3(v2, v1, edge); + MathExtra::sub3(x, v1, pvec); + MathExtra::cross3(edge, pvec, xproduct); + if (MathExtra::dot3(xproduct, norm) < 0.0) return 0; - MathExtra::sub3(v3,v2,edge); - MathExtra::sub3(x,v2,pvec); - MathExtra::cross3(edge,pvec,xproduct); - if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; + MathExtra::sub3(v3, v2, edge); + MathExtra::sub3(x, v2, pvec); + MathExtra::cross3(edge, pvec, xproduct); + if (MathExtra::dot3(xproduct, norm) < 0.0) return 0; - MathExtra::sub3(v1,v3,edge); - MathExtra::sub3(x,v3,pvec); - MathExtra::cross3(edge,pvec,xproduct); - if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; + MathExtra::sub3(v1, v3, edge); + MathExtra::sub3(x, v3, pvec); + MathExtra::cross3(edge, pvec, xproduct); + if (MathExtra::dot3(xproduct, norm) < 0.0) return 0; return 1; } @@ -452,7 +499,7 @@ double RegPrism::closest(double *x, double *near, double *nearest, double dsq) double delx = x[0] - near[0]; double dely = x[1] - near[1]; 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; nearest[0] = near[0]; diff --git a/src/region_sphere.cpp b/src/region_sphere.cpp index 03eb762266..b67e805ac4 100644 --- a/src/region_sphere.cpp +++ b/src/region_sphere.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -23,52 +22,52 @@ using namespace LAMMPS_NS; -enum{CONSTANT,VARIABLE}; +enum { CONSTANT, VARIABLE }; /* ---------------------------------------------------------------------- */ 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_")) { - xstr = utils::strdup(arg[2]+2); + if (utils::strmatch(arg[2], "^v_")) { + xstr = utils::strdup(arg[2] + 2); xc = 0.0; xstyle = VARIABLE; varshape = 1; } else { - xc = xscale*utils::numeric(FLERR,arg[2],false,lmp); + xc = xscale * utils::numeric(FLERR, arg[2], false, lmp); xstyle = CONSTANT; } - if (utils::strmatch(arg[3],"^v_")) { - ystr = utils::strdup(arg[3]+2); + if (utils::strmatch(arg[3], "^v_")) { + ystr = utils::strdup(arg[3] + 2); yc = 0.0; ystyle = VARIABLE; varshape = 1; } else { - yc = yscale*utils::numeric(FLERR,arg[3],false,lmp); + yc = yscale * utils::numeric(FLERR, arg[3], false, lmp); ystyle = CONSTANT; } - if (utils::strmatch(arg[4],"^v_")) { - zstr = utils::strdup(arg[4]+2); + if (utils::strmatch(arg[4], "^v_")) { + zstr = utils::strdup(arg[4] + 2); zc = 0.0; zstyle = VARIABLE; varshape = 1; } else { - zc = zscale*utils::numeric(FLERR,arg[4],false,lmp); + zc = zscale * utils::numeric(FLERR, arg[4], false, lmp); zstyle = CONSTANT; } - if (utils::strmatch(arg[5],"^v_")) { - rstr = utils::strdup(arg[5]+2); + if (utils::strmatch(arg[5], "^v_")) { + rstr = utils::strdup(arg[5] + 2); radius = 0.0; rstyle = VARIABLE; varshape = 1; } else { - radius = xscale*utils::numeric(FLERR,arg[5],false,lmp); + radius = xscale * utils::numeric(FLERR, arg[5], false, lmp); rstyle = CONSTANT; } @@ -79,7 +78,7 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : // 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 // 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_zlo = zc - radius; extent_zhi = zc + radius; - } else bboxflag = 0; + } else + bboxflag = 0; cmax = 1; contact = new Contact[cmax]; @@ -103,11 +103,11 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : RegSphere::~RegSphere() { - delete [] xstr; - delete [] ystr; - delete [] zstr; - delete [] rstr; - delete [] contact; + delete[] xstr; + delete[] ystr; + delete[] zstr; + delete[] rstr; + delete[] contact; } /* ---------------------------------------------------------------------- */ @@ -128,7 +128,7 @@ int RegSphere::inside(double x, double y, double z) double delx = x - xc; double dely = y - yc; 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; return 0; @@ -146,15 +146,15 @@ int RegSphere::surface_interior(double *x, double cutoff) double delx = x[0] - xc; double dely = x[1] - yc; 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; double delta = radius - r; if (delta < cutoff) { contact[0].r = delta; - contact[0].delx = delx*(1.0-radius/r); - contact[0].dely = dely*(1.0-radius/r); - contact[0].delz = delz*(1.0-radius/r); + contact[0].delx = delx * (1.0 - radius / r); + contact[0].dely = dely * (1.0 - radius / r); + contact[0].delz = delz * (1.0 - radius / r); contact[0].radius = -radius; contact[0].iwall = 0; contact[0].varflag = 1; @@ -174,15 +174,15 @@ int RegSphere::surface_exterior(double *x, double cutoff) double delx = x[0] - xc; double dely = x[1] - yc; 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; double delta = r - radius; if (delta < cutoff) { contact[0].r = delta; - contact[0].delx = delx*(1.0-radius/r); - contact[0].dely = dely*(1.0-radius/r); - contact[0].delz = delz*(1.0-radius/r); + contact[0].delx = delx * (1.0 - radius / r); + contact[0].dely = dely * (1.0 - radius / r); + contact[0].delz = delz * (1.0 - radius / r); contact[0].radius = radius; contact[0].iwall = 0; contact[0].varflag = 1; @@ -197,19 +197,15 @@ int RegSphere::surface_exterior(double *x, double cutoff) void RegSphere::shape_update() { - if (xstyle == VARIABLE) - xc = xscale * input->variable->compute_equal(xvar); + if (xstyle == VARIABLE) xc = xscale * input->variable->compute_equal(xvar); - if (ystyle == VARIABLE) - yc = yscale * input->variable->compute_equal(yvar); + if (ystyle == VARIABLE) yc = yscale * input->variable->compute_equal(yvar); - if (zstyle == VARIABLE) - zc = zscale * input->variable->compute_equal(zvar); + if (zstyle == VARIABLE) zc = zscale * input->variable->compute_equal(zvar); if (rstyle == VARIABLE) { radius = xscale * input->variable->compute_equal(rvar); - if (radius < 0.0) - error->one(FLERR,"Variable evaluation in region gave bad value"); + if (radius < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value"); } } @@ -221,34 +217,30 @@ void RegSphere::variable_check() { if (xstyle == VARIABLE) { xvar = input->variable->find(xstr); - if (xvar < 0) - error->all(FLERR,"Variable name for region sphere does not exist"); + if (xvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); 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) { yvar = input->variable->find(ystr); - if (yvar < 0) - error->all(FLERR,"Variable name for region sphere does not exist"); + if (yvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); 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) { zvar = input->variable->find(zstr); - if (zvar < 0) - error->all(FLERR,"Variable name for region sphere does not exist"); + if (zvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); 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) { rvar = input->variable->find(rstr); - if (rvar < 0) - error->all(FLERR,"Variable name for region sphere does not exist"); + if (rvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); 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[2] = zc; forward_transform(xcenter[0], xcenter[1], xcenter[2]); - if (update->ntimestep > 0) rprev = prev[4]; - else rprev = radius; + if (update->ntimestep > 0) + rprev = prev[4]; + else + rprev = radius; prev[4] = radius; } - - /* ---------------------------------------------------------------------- add velocity due to shape change to wall velocity ------------------------------------------------------------------------- */ 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); - dely = (xc[1] - xcenter[1])*(1 - rprev/radius); - delz = (xc[2] - xcenter[2])*(1 - rprev/radius); + delx = (xc[0] - xcenter[0]) * (1 - rprev / radius); + dely = (xc[1] - xcenter[1]) * (1 - rprev / radius); + delz = (xc[2] - xcenter[2]) * (1 - rprev / radius); - vwall[0] += delx/update->dt; - vwall[1] += dely/update->dt; - vwall[2] += delz/update->dt; + vwall[0] += delx / update->dt; + vwall[1] += dely / update->dt; + vwall[2] += delz / update->dt; } -