diff --git a/doc/src/fix_indent.rst b/doc/src/fix_indent.rst index d478d7dc50..5658c06373 100644 --- a/doc/src/fix_indent.rst +++ b/doc/src/fix_indent.rst @@ -8,13 +8,12 @@ Syntax .. code-block:: LAMMPS - fix ID group-ID indent K keyword values ... + fix ID group-ID indent K gstyle args keyword value ... * ID, group-ID are documented in :doc:`fix ` command * indent = style name of this fix command * K = force constant for indenter surface (force/distance\^2 units) -* one or more keyword/value pairs may be appended -* keyword = *sphere* or *cone* or *cylinder* or *plane* or *side* or *units* +* gstyle = *sphere* or *cylinder* or *cone* or *plane* .. parsed-literal:: @@ -22,22 +21,28 @@ Syntax x, y, z = position of center of indenter (distance units) R = sphere radius of indenter (distance units) any of x, y, z, R can be a variable (see below) + *cylinder* args = dim c1 c2 R + dim = *x* or *y* or *z* = axis of cylinder + c1, c2 = coords of cylinder axis in other 2 dimensions (distance units) + R = cylinder radius of indenter (distance units) + any of c1,c2,R can be a variable (see below) *cone* args = dim c1 c2 radlo radhi lo hi dim = *x* or *y* or *z* = axis of cone c1, c2 = coords of cone axis in other 2 dimensions (distance units) radlo,radhi = cone radii at lo and hi end (distance units) lo,hi = bounds of cone in dim (distance units) any of c1, c2, radlo, radhi, lo, hi can be a variable (see below) - *cylinder* args = dim c1 c2 R - dim = *x* or *y* or *z* = axis of cylinder - c1, c2 = coords of cylinder axis in other 2 dimensions (distance units) - R = cylinder radius of indenter (distance units) - any of c1,c2,R can be a variable (see below) *plane* args = dim pos side dim = *x* or *y* or *z* = plane perpendicular to this dimension pos = position of plane in dimension x, y, or z (distance units) pos can be a variable (see below) side = *lo* or *hi* + +* zero or more keyword/value pairs may be appended +* keyword = *side* or *units* + + .. parsed-literal:: + *side* value = *in* or *out* *in* = the indenter acts on particles inside the sphere or cylinder *out* = the indenter acts on particles outside the sphere or cylinder @@ -63,8 +68,8 @@ material or as an obstacle in a flow. Or it can be used as a constraining wall around a simulation; see the discussion of the *side* keyword below. -The indenter can either be spherical or conical or cylindrical or planar. You -must set one of those 3 keywords. +The *gstyle* geometry of the indenter can either be a sphere, a +cylinder, a cone, or a plane. A spherical indenter exerts a force of magnitude @@ -81,15 +86,20 @@ A cylindrical indenter exerts the same force, except that *r* is the distance from the atom to the center axis of the cylinder. The cylinder extends infinitely along its axis. -Spherical, conical and cylindrical indenters account for periodic boundaries in -two ways. First, the center point of a spherical indenter (x,y,z) or -axis of a conical/cylindrical indenter (c1,c2) is remapped back into the -simulation box, if the box is periodic in a particular dimension. -This occurs every timestep if the indenter geometry is specified with -a variable (see below), e.g. it is moving over time. Second, the -calculation of distance to the indenter center or axis accounts for -periodic boundaries. Both of these mean that an indenter can -effectively move through and straddle one or more periodic boundaries. +A conical indenter is similar to a cylindrical indenter except that it +has a finite length (between *lo* and *hi*), and that two different +radii (one at each end, *radlo* and *radhi*) can be defined. + +Spherical, cylindrical, and conical indenters account for periodic +boundaries in two ways. First, the center point of a spherical +indenter (x,y,z) or axis of a cylindrical/conical indenter (c1,c2) is +remapped back into the simulation box, if the box is periodic in a +particular dimension. This occurs every timestep if the indenter +geometry is specified with a variable (see below), e.g. it is moving +over time. Second, the calculation of distance to the indenter center +or axis accounts for periodic boundaries. Both of these mean that an +indenter can effectively move through and straddle one or more +periodic boundaries. A planar indenter is really an axis-aligned infinite-extent wall exerting the same force on atoms in the system, where *R* is the @@ -103,9 +113,13 @@ is specified as *hi*\ . Any of the 4 quantities defining a spherical indenter's geometry can be specified as an equal-style :doc:`variable `, namely *x*, -*y*, *z*, or *R*\ . Similarly, for a cylindrical indenter, any of *c1*, -*c2*, or *R*, can be a variable. For a planar indenter, *pos* can be -a variable. If the value is a variable, it should be specified as +*y*, *z*, or *R*\ . For a cylindrical indenter, any of the 3 +quantities *c1*, *c2*, or *R*, can be a variable. For a conical +indenter, any of the 6 quantities *c1*, *c2*, *radlo*, *radhi*, *lo*, +or *hi* can be a variable. For a planar indenter, the single value +*pos* can be a variable. + +If any of these values is a variable, it should be specified as v_name, where name is the variable name. In this case, the variable will be evaluated each timestep, and its value used to define the indenter geometry. @@ -116,7 +130,8 @@ command keywords for the simulation box parameters and timestep and elapsed time. Thus it is easy to specify indenter properties that change as a function of time or span consecutive runs in a continuous fashion. For the latter, see the *start* and *stop* keywords of the -:doc:`run ` command and the *elaplong* keyword of :doc:`thermo_style custom ` for details. +:doc:`run ` command and the *elaplong* keyword of +:doc:`thermo_style custom ` for details. For example, if a spherical indenter's x-position is specified as v_x, then this variable definition will keep it's center at a relative @@ -147,12 +162,13 @@ rate. If the *side* keyword is specified as *out*, which is the default, then particles outside the indenter are pushed away from its outer -surface, as described above. This only applies to spherical or -cylindrical indenters. If the *side* keyword is specified as *in*, -the action of the indenter is reversed. Particles inside the indenter -are pushed away from its inner surface. In other words, the indenter -is now a containing wall that traps the particles inside it. If the -radius shrinks over time, it will squeeze the particles. +surface, as described above. This only applies to spherical, +cylindrical, and conical indenters. If the *side* keyword is +specified as *in*, the action of the indenter is reversed. Particles +inside the indenter are pushed away from its inner surface. In other +words, the indenter is now a containing wall that traps the particles +inside it. If the radius shrinks over time, it will squeeze the +particles. The *units* keyword determines the meaning of the distance units used to define the indenter geometry. A *box* value selects standard @@ -172,10 +188,10 @@ lattice spacings in a variable formula. The force constant *K* is not affected by the *units* keyword. It is always in force/distance\^2 units where force and distance are defined -by the :doc:`units ` command. If you wish K to be scaled by the -lattice spacing, you can define K with a variable whose formula -contains *xlat*, *ylat*, *zlat* keywords of the -:doc:`thermo_style ` command, e.g. +by the :doc:`units ` command. If you wish K to be scaled by +the lattice spacing, you can define K with a variable whose formula +contains *xlat*, *ylat*, *zlat* keywords of the :doc:`thermo_style +` command, e.g. .. code-block:: LAMMPS diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 65cf91c2d2..8d450bee75 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -38,21 +38,6 @@ using namespace FixConst; enum{NONE, SPHERE, CYLINDER, PLANE, CONE}; enum{INSIDE, OUTSIDE}; -static bool PointInsideCone(int dir, double *center, double lo, - double hi, double rlo, double rhi, double *point); - -static void DistanceExteriorPoint(int dir, double *center, - double lo, double hi, double rlo, double rhi, double &x, - double &y, double &z); - -static void DistanceInteriorPoint(int dir, double *center, - double lo, double hi, double rlo, double rhi, double &x, - double &y, double &z); - -static void point_on_line_segment(double *a, double *b, double *c, double *d); - -static double closest(double *x, double *near, double *nearest, double dsq); - /* ---------------------------------------------------------------------- */ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : @@ -76,10 +61,11 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : if (k < 0.0) error->all(FLERR, "Illegal fix indent force constant: {}", k); k3 = k/3.0; - // read options from end of input line - - options(narg-4,&arg[4]); + // read geometry of indenter and optional args + int iarg = geometry(narg-4,&arg[4]) + 4; + options(narg-iarg,&arg[iarg]); + // setup scaling const double xscale { scaleflag ? domain->lattice->xlattice : 1.0}; @@ -93,8 +79,8 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : if (!ystr) yvalue *= yscale; if (!zstr) zvalue *= zscale; if (!rstr) rvalue *= xscale; + } else if (istyle == CONE) { - if (!xstr) xvalue *= xscale; if (!ystr) yvalue *= yscale; if (!zstr) zvalue *= zscale; @@ -121,6 +107,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : if (cdim == 0 && !pstr) pvalue *= xscale; else if (cdim == 1 && !pstr) pvalue *= yscale; else if (cdim == 2 && !pstr) pvalue *= zscale; + } else error->all(FLERR,"Unknown fix indent keyword: {}", istyle); varflag = 0; @@ -195,7 +182,6 @@ void FixIndent::init() if (!input->variable->equalstyle(pvar)) error->all(FLERR,"Variable {} for fix indent is invalid style", pstr); } - if (rlostr) { rlovar = input->variable->find(rlostr); if (rlovar < 0) @@ -203,7 +189,6 @@ void FixIndent::init() if (!input->variable->equalstyle(rlovar)) error->all(FLERR,"Variable {} for fix indent is invalid style", rlostr); } - if (rhistr) { rhivar = input->variable->find(rhistr); if (rhivar < 0) @@ -211,7 +196,6 @@ void FixIndent::init() if (!input->variable->equalstyle(rhivar)) error->all(FLERR,"Variable {} for fix indent is invalid style", rhistr); } - if (lostr) { lovar = input->variable->find(lostr); if (lovar < 0) @@ -219,7 +203,6 @@ void FixIndent::init() if (!input->variable->equalstyle(lovar)) error->all(FLERR,"Variable {} for fix indent is invalid style", lostr); } - if (histr) { hivar = input->variable->find(histr); if (hivar < 0) @@ -267,6 +250,7 @@ void FixIndent::post_force(int /*vflag*/) indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0; // ctr = current indenter centerz + double ctr[3] {xvalue, yvalue, zvalue}; if (xstr) ctr[0] = input->variable->compute_equal(xvar); if (ystr) ctr[1] = input->variable->compute_equal(yvar); @@ -284,6 +268,7 @@ void FixIndent::post_force(int /*vflag*/) if (istyle == SPHERE) { // remap indenter center into periodic box + domain->remap(ctr); double radius { rstr ? input->variable->compute_equal(rvar) : rvalue}; @@ -387,20 +372,24 @@ void FixIndent::post_force(int /*vflag*/) double x0[3] {delx + ctr[0], dely + ctr[1], delz + ctr[2]}; r = sqrt(delx * delx + dely * dely + delz * delz); - // find if the particle is inside or outside the cone + // check if particle is inside or outside the cone + bool point_inside_cone = PointInsideCone(cdim, ctr, lo, hi, radiuslo, radiushi, x0); if (side == INSIDE && point_inside_cone) continue; if (side == OUTSIDE && !point_inside_cone) continue; // find the distance between the point and the cone + if (point_inside_cone) { DistanceInteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]); } else { DistanceExteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]); } - // compute the force from the center of the cone - it is different from the approach of fix wall/region + // compute the force from the center of the cone + // this is different from how it is done in fix wall/region + dr = sqrt(x0[0] * x0[0] + x0[1] * x0[1] + x0[2] * x0[2]); int force_sign = { point_inside_cone ? 1 : -1 }; @@ -486,10 +475,10 @@ double FixIndent::compute_vector(int n) } /* ---------------------------------------------------------------------- - parse optional parameters at end of input line + parse input args for geometry of indenter ------------------------------------------------------------------------- */ -void FixIndent::options(int narg, char **arg) +int FixIndent::geometry(int narg, char **arg) { if (narg < 0) utils::missing_cmd_args(FLERR, "fix indent", error); @@ -499,139 +488,168 @@ void FixIndent::options(int narg, char **arg) scaleflag = 1; side = OUTSIDE; + // sphere + + if (strcmp(arg[0],"sphere") == 0) { + if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); + if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent sphere", error); + + if (utils::strmatch(arg[1],"^v_")) { + xstr = utils::strdup(arg[1]+2); + } else xvalue = utils::numeric(FLERR,arg[1],false,lmp); + if (utils::strmatch(arg[2],"^v_")) { + ystr = utils::strdup(arg[2]+2); + } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + if (utils::strmatch(arg[4],"^v_")) { + rstr = utils::strdup(arg[4]+2); + } else rvalue = utils::numeric(FLERR,arg[4],false,lmp); + + istyle = SPHERE; + return 5; + } + + // cylinder + + if (strcmp(arg[0],"cylinder") == 0) { + if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); + if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error); + + if (strcmp(arg[1],"x") == 0) { + cdim = 0; + if (utils::strmatch(arg[2],"^v_")) { + ystr = utils::strdup(arg[2]+2); + } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + } else if (strcmp(arg[1],"y") == 0) { + cdim = 1; + if (utils::strmatch(arg[2],"^v_")) { + xstr = utils::strdup(arg[2]+2); + } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + } else if (strcmp(arg[1],"z") == 0) { + cdim = 2; + if (utils::strmatch(arg[2],"^v_")) { + xstr = utils::strdup(arg[2]+2); + } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); + if (utils::strmatch(arg[3],"^v_")) { + ystr = utils::strdup(arg[3]+2); + } else yvalue = utils::numeric(FLERR,arg[3],false,lmp); + } else error->all(FLERR,"Unknown fix indent cylinder argument: {}", arg[1]); + + if (utils::strmatch(arg[4],"^v_")) { + rstr = utils::strdup(arg[4]+2); + } else rvalue = utils::numeric(FLERR,arg[4],false,lmp); + + istyle = CYLINDER; + return 5; + } + + // cone + + if (strcmp(arg[0],"cone") == 0) { + if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); + if (8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error); + + if (strcmp(arg[1],"x") == 0) { + cdim = 0; + + if (utils::strmatch(arg[2],"^v_")) { + ystr = utils::strdup(arg[2]+2); + } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); + + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + + } else if (strcmp(arg[1],"y") == 0) { + cdim = 1; + + if (utils::strmatch(arg[2],"^v_")) { + xstr = utils::strdup(arg[2]+2); + } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); + + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + + } else if (strcmp(arg[1],"z") == 0) { + cdim = 2; + + if (utils::strmatch(arg[2],"^v_")) { + xstr = utils::strdup(arg[2]+2); + } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); + + if (utils::strmatch(arg[3],"^v_")) { + ystr = utils::strdup(arg[3]+2); + } else yvalue = utils::numeric(FLERR,arg[3],false,lmp); + + } else error->all(FLERR,"Unknown fix indent cone argument: {}", arg[1]); + + if (utils::strmatch(arg[4],"^v_")) { + rlostr = utils::strdup(arg[4]+2); + } else rlovalue = utils::numeric(FLERR,arg[4],false,lmp); + + if (utils::strmatch(arg[5],"^v_")) { + rhistr = utils::strdup(arg[5]+2); + } else rhivalue = utils::numeric(FLERR,arg[5],false,lmp); + + if (utils::strmatch(arg[6],"^v_")) { + lostr = utils::strdup(arg[6]+2); + } else lovalue = utils::numeric(FLERR,arg[6],false,lmp); + + if (utils::strmatch(arg[7],"^v_")) { + histr = utils::strdup(arg[7]+2); + } else hivalue = utils::numeric(FLERR,arg[7],false,lmp); + + istyle = CONE; + return 8; + } + + // plane + + if (strcmp(arg[0],"plane") == 0) { + if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); + if (4 > narg) utils::missing_cmd_args(FLERR, "fix indent plane", error); + if (strcmp(arg[1],"x") == 0) cdim = 0; + else if (strcmp(arg[1],"y") == 0) cdim = 1; + else if (strcmp(arg[1],"z") == 0) cdim = 2; + else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[1]); + + if (utils::strmatch(arg[2],"^v_")) { + pstr = utils::strdup(arg[2]+2); + } else pvalue = utils::numeric(FLERR,arg[2],false,lmp); + + if (strcmp(arg[3],"lo") == 0) planeside = -1; + else if (strcmp(arg[3],"hi") == 0) planeside = 1; + else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[3]); + istyle = PLANE; + return 4; + } + + // invalid istyle arg + + error->all(FLERR,"Unknown fix indent argument: {}", arg[0]); + + return 0; +} + +/* ---------------------------------------------------------------------- + parse optional input args +------------------------------------------------------------------------- */ + +void FixIndent::options(int narg, char **arg) +{ int iarg = 0; + while (iarg < narg) { - if (strcmp(arg[iarg],"sphere") == 0) { - if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "fix indent sphere", error); - - if (utils::strmatch(arg[iarg+1],"^v_")) { - xstr = utils::strdup(arg[iarg+1]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (utils::strmatch(arg[iarg+2],"^v_")) { - ystr = utils::strdup(arg[iarg+2]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (utils::strmatch(arg[iarg+4],"^v_")) { - rstr = utils::strdup(arg[iarg+4]+2); - } else rvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); - - istyle = SPHERE; - iarg += 5; - - } else if (strcmp(arg[iarg],"cylinder") == 0) { - if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error); - - if (strcmp(arg[iarg+1],"x") == 0) { - cdim = 0; - if (utils::strmatch(arg[iarg+2],"^v_")) { - ystr = utils::strdup(arg[iarg+2]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - } else if (strcmp(arg[iarg+1],"y") == 0) { - cdim = 1; - if (utils::strmatch(arg[iarg+2],"^v_")) { - xstr = utils::strdup(arg[iarg+2]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - } else if (strcmp(arg[iarg+1],"z") == 0) { - cdim = 2; - if (utils::strmatch(arg[iarg+2],"^v_")) { - xstr = utils::strdup(arg[iarg+2]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) { - ystr = utils::strdup(arg[iarg+3]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - } else error->all(FLERR,"Unknown fix indent cylinder argument: {}", arg[iarg+1]); - - if (utils::strmatch(arg[iarg+4],"^v_")) { - rstr = utils::strdup(arg[iarg+4]+2); - } else rvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); - - istyle = CYLINDER; - iarg += 5; - - } else if (strcmp(arg[iarg],"plane") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "fix indent plane", error); - if (strcmp(arg[iarg+1],"x") == 0) cdim = 0; - else if (strcmp(arg[iarg+1],"y") == 0) cdim = 1; - else if (strcmp(arg[iarg+1],"z") == 0) cdim = 2; - else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[iarg+1]); - - if (utils::strmatch(arg[iarg+2],"^v_")) { - pstr = utils::strdup(arg[iarg+2]+2); - } else pvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - - if (strcmp(arg[iarg+3],"lo") == 0) planeside = -1; - else if (strcmp(arg[iarg+3],"hi") == 0) planeside = 1; - else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[iarg+3]); - istyle = PLANE; - iarg += 4; - - } else if (strcmp(arg[iarg],"cone") == 0) { - - if (iarg+8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error); - - if (strcmp(arg[iarg+1],"x") == 0) { - cdim = 0; - - if (utils::strmatch(arg[iarg+2],"^v_")) { - ystr = utils::strdup(arg[iarg+2]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - - } else if (strcmp(arg[iarg+1],"y") == 0) { - cdim = 1; - - if (utils::strmatch(arg[iarg+2],"^v_")) { - xstr = utils::strdup(arg[iarg+2]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - - } else if (strcmp(arg[iarg+1],"z") == 0) { - cdim = 2; - - if (utils::strmatch(arg[iarg+2],"^v_")) { - xstr = utils::strdup(arg[iarg+2]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - - if (utils::strmatch(arg[iarg+3],"^v_")) { - ystr = utils::strdup(arg[iarg+3]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - - } else error->all(FLERR,"Unknown fix indent cone argument: {}", arg[iarg+1]); - - if (utils::strmatch(arg[iarg+4],"^v_")) { - rlostr = utils::strdup(arg[iarg+4]+2); - } else rlovalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); - - if (utils::strmatch(arg[iarg+5],"^v_")) { - rhistr = utils::strdup(arg[iarg+5]+2); - } else rhivalue = utils::numeric(FLERR,arg[iarg+5],false,lmp); - - if (utils::strmatch(arg[iarg+6],"^v_")) { - lostr = utils::strdup(arg[iarg+6]+2); - } else lovalue = utils::numeric(FLERR,arg[iarg+6],false,lmp); - - if (utils::strmatch(arg[iarg+7],"^v_")) { - histr = utils::strdup(arg[iarg+7]+2); - } else hivalue = utils::numeric(FLERR,arg[iarg+7],false,lmp); - - istyle = CONE; - iarg += 8; - - } else if (strcmp(arg[iarg],"units") == 0) { + if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix indent units", error); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; @@ -644,16 +662,17 @@ void FixIndent::options(int narg, char **arg) else if (strcmp(arg[iarg+1],"out") == 0) side = OUTSIDE; else error->all(FLERR,"Unknown fix indent side argument: {}", arg[iarg+1]); iarg += 2; + } else error->all(FLERR,"Unknown fix indent argument: {}", arg[iarg]); } } /* ---------------------------------------------------------------------- - determines if a point is inside (true) or outside (false) of a cone + determines if a point is inside (true) or outside (false) of a cone ------------------------------------------------------------------------- */ -bool PointInsideCone(int dir, double *center, double lo, - double hi, double rlo, double rhi, double *x) +bool FixIndent::PointInsideCone(int dir, double *center, double lo, + double hi, double rlo, double rhi, double *x) { if ((x[dir] > hi) || (x[dir] < lo)) return false; @@ -669,10 +688,12 @@ bool PointInsideCone(int dir, double *center, double lo, } /* ---------------------------------------------------------------------- - distance between an exterior point and a cone + distance between an exterior point and a cone ------------------------------------------------------------------------- */ -void DistanceExteriorPoint(int dir, double *center, double lo, double hi, - double rlo, double rhi, double &x, double &y, double &z) + +void FixIndent::DistanceExteriorPoint(int dir, double *center, double lo, double hi, + double rlo, double rhi, + double &x, double &y, double &z) { double xp[3], nearest[3], corner1[3], corner2[3]; @@ -700,17 +721,21 @@ void DistanceExteriorPoint(int dir, double *center, double lo, double hi, corner4[dir] = hi; // initialize distance to a big number + double distsq = 1.0e20; // check the first triangle + point_on_line_segment(corner1, corner2, point, xp); distsq = closest(point, xp, nearest, distsq); // check the second triangle + point_on_line_segment(corner1, corner3, point, xp); distsq = closest(point, xp, nearest, distsq); // check the third triangle + point_on_line_segment(corner2, corner4, point, xp); distsq = closest(point, xp, nearest, distsq); @@ -722,12 +747,12 @@ void DistanceExteriorPoint(int dir, double *center, double lo, double hi, } /* ---------------------------------------------------------------------- - distance between an interior point and a cone + distance between an interior point and a cone ------------------------------------------------------------------------- */ -void DistanceInteriorPoint(int dir, double *center, - double lo, double hi, double rlo, double rhi, double &x, - double &y, double &z) +void FixIndent::DistanceInteriorPoint(int dir, double *center, + double lo, double hi, double rlo, double rhi, double &x, + double &y, double &z) { double r, dist_disk, dist_surf; double surflo[3], surfhi[3], xs[3]; @@ -735,6 +760,7 @@ void DistanceInteriorPoint(int dir, double *center, double point[3] {0.0, 0.0, 0.0}; // initial check with the two disks + if ( (initial_point[dir] - lo) < (hi - initial_point[dir]) ) { dist_disk = (initial_point[dir] - lo) * (initial_point[dir] - lo); point[dir] = initial_point[dir] - lo; @@ -744,6 +770,7 @@ void DistanceInteriorPoint(int dir, double *center, } // check with the points in the conical surface + double del[3] {x - center[0], y - center[1], z - center[2]}; del[dir] = 0.0; r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); @@ -776,10 +803,10 @@ void DistanceInteriorPoint(int dir, double *center, } /* ---------------------------------------------------------------------- - helper function extracted from region.cpp + helper function extracted from region.cpp ------------------------------------------------------------------------- */ -void point_on_line_segment(double *a, double *b, double *c, double *d) +void FixIndent::point_on_line_segment(double *a, double *b, double *c, double *d) { double ba[3], ca[3]; @@ -802,10 +829,10 @@ void point_on_line_segment(double *a, double *b, double *c, double *d) } /* ---------------------------------------------------------------------- - helper function extracted from region_cone.cpp + helper function extracted from region_cone.cpp ------------------------------------------------------------------------- */ -double closest(double *x, double *near, double *nearest, double dsq) +double FixIndent::closest(double *x, double *near, double *nearest, double dsq) { double dx = x[0] - near[0]; double dy = x[1] - near[1]; diff --git a/src/fix_indent.h b/src/fix_indent.h index 6f33f6fbb1..202a138729 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -53,7 +53,20 @@ class FixIndent : public Fix { int rlovar, rhivar, lovar, hivar; double rlovalue, rhivalue, lovalue, hivalue; + // methods for argument + + int geometry(int, char **); void options(int, char **); + + // methods for conical indenter + + bool PointInsideCone(int, double *, double, double, double, double, double *); + void DistanceExteriorPoint(int, double *, double, double, double, double, + double &, double &, double &); + void DistanceInteriorPoint(int, double *, double, double, double, double, + double &, double &, double &); + void point_on_line_segment(double *, double *, double *, double *); + double closest(double *, double *, double *, double); }; } // namespace LAMMPS_NS