diff --git a/doc/src/fix_indent.rst b/doc/src/fix_indent.rst index 15790e15d02..e041f9f29b1 100644 --- a/doc/src/fix_indent.rst +++ b/doc/src/fix_indent.rst @@ -8,33 +8,44 @@ 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 *cylinder* or *plane* or *side* or *units* +* gstyle = *sphere* or *cylinder* or *cone* or *plane* .. parsed-literal:: *sphere* args = x y z R - x,y,z = position of center of indenter (distance units) + 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) + 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) + 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) *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 + *in* = the indenter acts on particles inside the sphere or cylinder or cone + *out* = the indenter acts on particles outside the sphere or cylinder or cone *units* value = *lattice* or *box* lattice = the geometry is defined in lattice units box = the geometry is defined in simulation box units @@ -53,12 +64,12 @@ Description Insert an indenter within a simulation box. The indenter repels all atoms in the group that touch it, so it can be used to push into a -material or as an obstacle in a flow. Or it can be used as a +material or as an obstacle in a flow. Alternatively, 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 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 @@ -75,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 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 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 @@ -97,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. @@ -110,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 @@ -141,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 @@ -166,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 9adb337dd66..87ed5091bd9 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -23,6 +23,7 @@ #include "error.h" #include "input.h" #include "lattice.h" +#include "math_extra.h" #include "modify.h" #include "respa.h" #include "update.h" @@ -34,14 +35,15 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{NONE,SPHERE,CYLINDER,PLANE}; -enum{INSIDE,OUTSIDE}; +enum{NONE, SPHERE, CYLINDER, PLANE, CONE}; +enum{INSIDE, OUTSIDE}; /* ---------------------------------------------------------------------- */ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr), pstr(nullptr) + xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr), pstr(nullptr), + rlostr(nullptr), rhistr(nullptr), lostr(nullptr), histr(nullptr) { if (narg < 4) utils::missing_cmd_args(FLERR, "fix indent", error); @@ -56,21 +58,19 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : ilevel_respa = 0; k = utils::numeric(FLERR,arg[3],false,lmp); + if (k < 0.0) error->all(FLERR, "Illegal fix indent force constant: {}", k); k3 = k/3.0; - // read options from end of input line + // read geometry of indenter and optional args - options(narg-4,&arg[4]); + int iarg = geometry(narg-4,&arg[4]) + 4; + options(narg-iarg,&arg[iarg]); // setup scaling - double xscale,yscale,zscale; - if (scaleflag) { - xscale = domain->lattice->xlattice; - yscale = domain->lattice->ylattice; - zscale = domain->lattice->zlattice; - } - else xscale = yscale = zscale = 1.0; + const double xscale { scaleflag ? domain->lattice->xlattice : 1.0}; + const double yscale { scaleflag ? domain->lattice->ylattice : 1.0}; + const double zscale { scaleflag ? domain->lattice->zlattice : 1.0}; // apply scaling factors to geometry @@ -79,14 +79,39 @@ 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; + + double scaling_factor = 1.0; + switch (cdim) { + case 0: + scaling_factor = xscale; + break; + case 1: + scaling_factor = yscale; + break; + case 2: + scaling_factor = zscale; + break; + } + + if (!rlostr) rlovalue *= scaling_factor; + if (!rhistr) rhivalue *= scaling_factor; + if (!lostr) lovalue *= scaling_factor; + if (!histr) hivalue *= scaling_factor; + } else if (istyle == PLANE) { 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; - if (xstr || ystr || zstr || rstr || pstr) varflag = 1; + if (xstr || ystr || zstr || rstr || pstr || rlostr || rhistr || lostr || histr) varflag = 1; indenter_flag = 0; indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0; @@ -101,6 +126,10 @@ FixIndent::~FixIndent() delete [] zstr; delete [] rstr; delete [] pstr; + delete [] rlostr; + delete [] rhistr; + delete [] lostr; + delete [] histr; } /* ---------------------------------------------------------------------- */ @@ -153,6 +182,34 @@ 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) + error->all(FLERR,"Variable {} for fix indent does not exist", rlostr); + 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) + error->all(FLERR,"Variable {} for fix indent does not exist", rhistr); + 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) + error->all(FLERR,"Variable {} for fix indent does not exist", lostr); + 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) + error->all(FLERR,"Variable {} for fix indent does not exist", histr); + if (!input->variable->equalstyle(hivar)) + error->all(FLERR,"Variable {} for fix indent is invalid style", histr); + } if (utils::strmatch(update->integrate_style,"^respa")) { ilevel_respa = (dynamic_cast(update->integrate))->nlevels-1; @@ -192,32 +249,30 @@ void FixIndent::post_force(int /*vflag*/) indenter_flag = 0; 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); + if (zstr) ctr[2] = input->variable->compute_equal(zvar); + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double delx, dely, delz, r, dr, fmag, fx, fy, fz; + // spherical indenter if (istyle == SPHERE) { - // ctr = current indenter center - // remap into periodic box + // remap indenter center into periodic box - double ctr[3]; - if (xstr) ctr[0] = input->variable->compute_equal(xvar); - else ctr[0] = xvalue; - if (ystr) ctr[1] = input->variable->compute_equal(yvar); - else ctr[1] = yvalue; - if (zstr) ctr[2] = input->variable->compute_equal(zvar); - else ctr[2] = zvalue; domain->remap(ctr); - double radius; - if (rstr) radius = input->variable->compute_equal(rvar); - else radius = rvalue; - - double **x = atom->x; - double **f = atom->f; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - double delx,dely,delz,r,dr,fmag,fx,fy,fz; + double radius { rstr ? input->variable->compute_equal(rvar) : rvalue}; + if (radius < 0.0) error->all(FLERR, "Illegal fix indent sphere radius: {}", radius); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -254,56 +309,18 @@ void FixIndent::post_force(int /*vflag*/) // remap into periodic box // 3rd coord is just near box for remap(), since isn't used - double ctr[3]; - if (cdim == 0) { - ctr[0] = domain->boxlo[0]; - if (ystr) ctr[1] = input->variable->compute_equal(yvar); - else ctr[1] = yvalue; - if (zstr) ctr[2] = input->variable->compute_equal(zvar); - else ctr[2] = zvalue; - } else if (cdim == 1) { - if (xstr) ctr[0] = input->variable->compute_equal(xvar); - else ctr[0] = xvalue; - ctr[1] = domain->boxlo[1]; - if (zstr) ctr[2] = input->variable->compute_equal(zvar); - else ctr[2] = zvalue; - } else { - if (xstr) ctr[0] = input->variable->compute_equal(xvar); - else ctr[0] = xvalue; - if (ystr) ctr[1] = input->variable->compute_equal(yvar); - else ctr[1] = yvalue; - ctr[2] = domain->boxlo[2]; - } + ctr[cdim] = domain->boxlo[cdim]; domain->remap(ctr); - double radius; - if (rstr) radius = input->variable->compute_equal(rvar); - else radius = rvalue; - - double **x = atom->x; - double **f = atom->f; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - double delx,dely,delz,r,dr,fmag,fx,fy,fz; + double radius { rstr ? input->variable->compute_equal(rvar) : rvalue}; + if (radius < 0.0) error->all(FLERR, "Illegal fix indent cylinder radius: {}", radius); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (cdim == 0) { - delx = 0; - dely = x[i][1] - ctr[1]; - delz = x[i][2] - ctr[2]; - } else if (cdim == 1) { - delx = x[i][0] - ctr[0]; - dely = 0; - delz = x[i][2] - ctr[2]; - } else { - delx = x[i][0] - ctr[0]; - dely = x[i][1] - ctr[1]; - delz = 0; - } - domain->minimum_image(delx,dely,delz); - r = sqrt(delx*delx + dely*dely + delz*delz); + double del[3] {x[i][0] - ctr[0], x[i][1] - ctr[1], x[i][2] - ctr[2]}; + del[cdim] = 0; + domain->minimum_image(del[0], del[1], del[2]); + r = sqrt(del[0]*del[0] + del[1]*del[1] + del[2]*del[2]); if (side == OUTSIDE) { dr = r - radius; fmag = k*dr*dr; @@ -312,17 +329,83 @@ void FixIndent::post_force(int /*vflag*/) fmag = -k*dr*dr; } if (dr >= 0.0) continue; + fx = del[0]*fmag/r; + fy = del[1]*fmag/r; + fz = del[2]*fmag/r; + f[i][0] += fx; + f[i][1] += fy; + f[i][2] += fz; + indenter[0] -= k3 * dr*dr*dr; + indenter[1] -= fx; + indenter[2] -= fy; + indenter[3] -= fz; + } + + // conical indenter + + } else if (istyle == CONE) { + + double radiuslo { rlostr ? input->variable->compute_equal(rlovar) : rlovalue }; + if (radiuslo < 0.0) error->all(FLERR, "Illegal fix indent cone lower radius: {}", radiuslo); + double radiushi { rhistr ? input->variable->compute_equal(rhivar) : rhivalue }; + if (radiushi < 0.0) error->all(FLERR, "Illegal fix indent cone high radius: {}", radiushi); + + double initial_lo { lostr ? input->variable->compute_equal(lovar) : lovalue }; + double initial_hi { histr ? input->variable->compute_equal(hivar) : hivalue }; + + ctr[cdim] = 0.5 * (initial_hi + initial_lo); + + domain->remap(ctr); + + double hi = ctr[cdim] + 0.5 * (initial_hi - initial_lo); + double lo = ctr[cdim] - 0.5 * (initial_hi - initial_lo); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + + delx = x[i][0] - ctr[0]; + dely = x[i][1] - ctr[1]; + delz = x[i][2] - ctr[2]; + domain->minimum_image(delx, dely, delz); + + double x0[3] {delx + ctr[0], dely + ctr[1], delz + ctr[2]}; + r = sqrt(delx * delx + dely * dely + delz * delz); + + // 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 + // 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 }; + fmag = force_sign * k * dr * dr; + fx = delx*fmag/r; fy = dely*fmag/r; fz = delz*fmag/r; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; - indenter[0] -= k3 * dr*dr*dr; + indenter[0] -= k3 * dr * dr * dr; indenter[1] -= fx; indenter[2] -= fy; indenter[3] -= fz; } + } // planar indenter @@ -330,25 +413,16 @@ void FixIndent::post_force(int /*vflag*/) // plane = current plane position - double plane; - if (pstr) plane = input->variable->compute_equal(pvar); - else plane = pvalue; - - double **x = atom->x; - double **f = atom->f; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - double dr,fatom; + double plane { pstr ? input->variable->compute_equal(pvar) : pvalue}; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { dr = planeside * (plane - x[i][cdim]); if (dr >= 0.0) continue; - fatom = -planeside * k*dr*dr; - f[i][cdim] += fatom; - indenter[0] -= k3 * dr*dr*dr; - indenter[cdim+1] -= fatom; + fmag = -planeside * k * dr * dr; + f[i][cdim] += fmag; + indenter[0] -= k3 * dr * dr * dr; + indenter[cdim+1] -= fmag; } } @@ -400,94 +474,173 @@ 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); istyle = NONE; xstr = ystr = zstr = rstr = pstr = nullptr; xvalue = yvalue = zvalue = rvalue = pvalue = 0.0; + + // 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) +{ scaleflag = 1; side = OUTSIDE; 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],"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; @@ -500,6 +653,186 @@ 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 +------------------------------------------------------------------------- */ + +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; + + double del[3] {x[0] - center[0], x[1] - center[1], x[2] - center[2]}; + del[dir] = 0.0; + + double dist = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); + double currentradius = rlo + (x[dir] - lo) * (rhi - rlo) / (hi - lo); + + if (dist > currentradius) return false; + + return true; +} + +/* ---------------------------------------------------------------------- + distance between an exterior point and a cone +------------------------------------------------------------------------- */ + +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]; + + double point[3] {x, y, z}; + + double del[3] {x - center[0], y - center[1], z - center[2]}; + del[dir] = 0.0; + + double r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); + + corner1[0] = center[0] + del[0] * rlo / r; + corner1[1] = center[1] + del[1] * rlo / r; + corner1[2] = center[2] + del[2] * rlo / r; + corner1[dir] = lo; + + corner2[0] = center[0] + del[0] * rhi / r; + corner2[1] = center[1] + del[1] * rhi / r; + corner2[2] = center[2] + del[2] * rhi / r; + corner2[dir] = hi; + + double corner3[3] {center[0], center[1], center[2]}; + corner3[dir] = lo; + + double corner4[3] {center[0], center[1], center[2]}; + 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); + + x -= nearest[0]; + y -= nearest[1]; + z -= nearest[2]; + + return; +} + +/* ---------------------------------------------------------------------- + distance between an interior point and a cone +------------------------------------------------------------------------- */ + +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]; + double initial_point[3] {x, y, z}; + 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; + } else { + dist_disk = (hi - initial_point[dir]) * (hi - initial_point[dir]); + point[dir] = initial_point[dir] - hi; + } + + // 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]); + + surflo[0] = center[0] + del[0] * rlo / r; + surflo[1] = center[1] + del[1] * rlo / r; + surflo[2] = center[2] + del[2] * rlo / r; + surflo[dir] = lo; + + surfhi[0] = center[0] + del[0] * rhi / r; + surfhi[1] = center[1] + del[1] * rhi / r; + surfhi[2] = center[2] + del[2] * rhi / r; + surfhi[dir] = hi; + + point_on_line_segment(surflo, surfhi, initial_point, xs); + + double dx[3] {initial_point[0] - xs[0], initial_point[1] - xs[1], initial_point[2] - xs[2]}; + dist_surf = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; + if (dist_surf < dist_disk) { + x = dx[0]; + y = dx[1]; + z = dx[2]; + } else { + x = point[0]; + y = point[1]; + z = point[2]; + } + + return; +} + +/* ---------------------------------------------------------------------- + helper function extracted from region.cpp +------------------------------------------------------------------------- */ + +void FixIndent::point_on_line_segment(double *a, double *b, double *c, double *d) +{ + double ba[3], ca[3]; + + MathExtra::sub3(b, a, ba); + MathExtra::sub3(c, a, ca); + double t = MathExtra::dot3(ca, ba) / MathExtra::dot3(ba, ba); + if (t <= 0.0) { + d[0] = a[0]; + d[1] = a[1]; + d[2] = a[2]; + } else if (t >= 1.0) { + d[0] = b[0]; + d[1] = b[1]; + d[2] = b[2]; + } else { + d[0] = a[0] + t * ba[0]; + d[1] = a[1] + t * ba[1]; + d[2] = a[2] + t * ba[2]; + } +} + +/* ---------------------------------------------------------------------- + helper function extracted from region_cone.cpp +------------------------------------------------------------------------- */ + +double FixIndent::closest(double *x, double *near, double *nearest, double dsq) +{ + double dx = x[0] - near[0]; + double dy = x[1] - near[1]; + double dz = x[2] - near[2]; + double rsq = dx * dx + dy * dy + dz * dz; + if (rsq >= dsq) return dsq; + + nearest[0] = near[0]; + nearest[1] = near[1]; + nearest[2] = near[2]; + return rsq; +} diff --git a/src/fix_indent.h b/src/fix_indent.h index 527e9ec2772..37e1623df9d 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -49,7 +49,24 @@ class FixIndent : public Fix { int cdim, varflag; int ilevel_respa; + char *rlostr, *rhistr, *lostr, *histr; + int rlovar, rhivar, lovar, hivar; + double rlovalue, rhivalue, lovalue, hivalue; + + // methods for argument parsing + + 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