diff --git a/doc/src/create_atoms.rst b/doc/src/create_atoms.rst index fd241ecf73..7b8280b07b 100644 --- a/doc/src/create_atoms.rst +++ b/doc/src/create_atoms.rst @@ -53,6 +53,8 @@ Syntax factor = scale factor for setting atom radius *radthresh* value = radius (distance units) radius = threshold value for *mesh* to determine when to split triangles + *quasirandom* value = number density (inverse distance squared units) + number density = minimum number density for particles placed on *mesh* triangles *rotate* values = theta Rx Ry Rz theta = rotation angle for single molecule (degrees) Rx,Ry,Rz = rotation vector for single molecule @@ -148,6 +150,12 @@ one sphere per triangle. If the atom style in use allows to set a per-atom radius this radius is set to the average distance of the triangle vertices from its center times the value of the *radscale* keyword (default: 1.0). +The *mesh* style also has a *quasirandom* option which uses a quasirandom +sequence to distribute particles on mesh triangles using an approach +by :ref:`(Roberts) `. This option takes a minimal number +density as an argument. Particles are added to the triangle until this +number density is met or exceeded such that every triangle will have +at least one particle. .. note:: @@ -509,3 +517,10 @@ assigned the argument *type* as their atom type (when single atoms are being created). The other defaults are *remap* = no, *rotate* = random, *radscale* = 1.0, *radthresh* = x-lattice spacing, *overlap* not checked, *maxtry* = 10, and *units* = lattice. + +---------- + +.. _Roberts2019: + +**(Roberts)** R. Roberts (2019) "Evenly Distributing Points in a Triangle." Extreme Learning. +``_ diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index ad75284498..b84e9c6710 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -48,12 +48,15 @@ using MathConst::THIRD; static constexpr double BIG = 1.0e30; static constexpr double EPSILON = 1.0e-6; static constexpr double LB_FACTOR = 1.1; +static constexpr double INV_P_CONST = 0.7548777; +static constexpr double INV_SQ_P_CONST = 0.5698403; static constexpr int DEFAULT_MAXTRY = 1000; enum { BOX, REGION, SINGLE, RANDOM, MESH }; enum { ATOM, MOLECULE }; enum { COUNT, INSERT, INSERT_SELECTED }; enum { NONE, RATIO, SUBSET }; +enum { TRICENTER, QUASIRANDOM }; /* ---------------------------------------------------------------------- */ @@ -135,6 +138,7 @@ void CreateAtoms::command(int narg, char **arg) int scaleflag = 1; remapflag = 0; + mesh_style = TRICENTER; mode = ATOM; int molseed; ranmol = nullptr; @@ -261,6 +265,8 @@ void CreateAtoms::command(int narg, char **arg) if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "create_atoms radscale", error); if (style != MESH) error->all(FLERR, "Create_atoms radscale can only be used with mesh style"); + if (mesh_style != TRICENTER) + error->all(FLERR, "Create_atoms radscale can only be used with tricenter mesh style"); if (!atom->radius_flag) error->all(FLERR, "Must have atom attribute radius to set radscale factor"); radscale = utils::numeric(FLERR, arg[iarg + 1], false, lmp); @@ -269,8 +275,17 @@ void CreateAtoms::command(int narg, char **arg) if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "create_atoms radthresh", error); if (style != MESH) error->all(FLERR, "Create_atoms radthresh can only be used with mesh style"); + if (mesh_style != TRICENTER) + error->all(FLERR, "Create_atoms radthresh can only be used with tricenter mesh style"); radthresh = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; + } else if (strcmp(arg[iarg], "quasirandom") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "create_atoms quasirandom", error); + if (style != MESH) + error->all(FLERR, "Create_atoms quasirandom can only be used with mesh style"); + mesh_style = QUASIRANDOM; + mesh_density = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; } else error->all(FLERR, "Illegal create_atoms command option {}", arg[iarg]); } @@ -925,6 +940,86 @@ int CreateAtoms::add_tricenter(const double vert[3][3], tagint molid) return ilocal; } +/* ---------------------------------------------------------------------- + add monodisperse atoms by mapping quasirandom sequence onto a triangle + method by Martin Roberts 2019 http://extremelearning.com.au/evenly-distributing-points-in-a-triangle/ +------------------------------------------------------------------------- */ + +int CreateAtoms::add_quasirandom(const double vert[3][3], tagint molid) +{ + double ab[3], ac[3], bc[3], temp[3], point[3], ref[3]; + double lab, lac, lbc, area, xi, yi; + double seed = 0.5; + + // Order vertices such that a has the largest angle + MathExtra::sub3(vert[1], vert[0], ab); + MathExtra::sub3(vert[2], vert[0], ac); + MathExtra::sub3(vert[2], vert[1], bc); + + lab = MathExtra::len3(ab); + lac = MathExtra::len3(ac); + lbc = MathExtra::len3(bc); + + if (lac > lab && lac > lbc) { + // B has the largest angle, relabel as A + MathExtra::scale3(-1.0, ab); + MathExtra::copy3(bc, ac); + MathExtra::copy3(vert[1], ref); + } else if (lab > lac && lab > lbc) { + // C has the largest angle, relabel as A + MathExtra::scale3(-1.0, ac); + MathExtra::copy3(bc, ab); + MathExtra::scale3(-1.0, ab); + MathExtra::copy3(vert[2], ref); + } else { + MathExtra::copy3(vert[0], ref); + } + + // Estimate number of particles from area, could scale by RCP + MathExtra::cross3(ab, ac, temp); + area = 0.5 * MathExtra::len3(temp); + int nparticles = ceil(mesh_density * area); + + for (int i = 0; i < nparticles; i++) { + // Define point in unit square + xi = (i + 1) * INV_P_CONST; + yi = (i + 1) * INV_SQ_P_CONST; // Add seed function of vertices + + xi += seed; + yi += seed; + + xi = std::fmod(xi, 1.0); + yi = std::fmod(yi, 1.0); + + // Map to tirangle using parallelogram method + if ((xi + yi) < 1) { + MathExtra::scale3(xi, ac, point); + MathExtra::scale3(yi, ab, temp); + MathExtra::add3(point, temp, point); + } else { + xi = 1.0 - xi; + yi = 1.0 - yi; + MathExtra::scale3(xi, ac, point); + MathExtra::scale3(yi, ab, temp); + MathExtra::add3(point, temp, point); + } + + MathExtra::add3(point, ref, point); + + // if atom/molecule is in my subbox, create it + if ((point[0] >= sublo[0]) && (point[0] < subhi[0]) && (point[1] >= sublo[1]) && + (point[1] < subhi[1]) && (point[2] >= sublo[2]) && (point[2] < subhi[2])) { + + atom->avec->create_atom(ntype, point); + int idx = atom->nlocal - 1; + if (atom->molecule_flag) atom->molecule[idx] = molid; + } + + } + + return nparticles; +} + /* ---------------------------------------------------------------------- add atoms at center of triangulated mesh ------------------------------------------------------------------------- */ @@ -990,10 +1085,14 @@ void CreateAtoms::add_mesh(const char *filename) if (!line || !utils::strmatch(line, "^ *endfacet")) throw TokenizerException("Error reading endfacet", ""); - // now we have three vertices ... proceed with adding atom in center of triangle - // or splitting recursively into halves as needed to get the desired density ++ntriangle; - atomlocal += add_tricenter(vert, molid); + if (mesh_style == TRICENTER) { + // now we have three vertices ... proceed with adding atom in center of triangle + // or splitting recursively into halves as needed to get the desired density + atomlocal += add_tricenter(vert, molid); + } else if (mesh_style == QUASIRANDOM) { + atomlocal += add_quasirandom(vert, molid); + } } } catch (std::exception &e) { error->all(FLERR, "Error reading triangles from file {}: {}", filename, e.what()); diff --git a/src/create_atoms.h b/src/create_atoms.h index f9d96f24a3..dcd0a87d94 100644 --- a/src/create_atoms.h +++ b/src/create_atoms.h @@ -41,7 +41,7 @@ class CreateAtoms : public Command { double subsetfrac; int *basistype; double xone[3], quatone[4]; - double radthresh, radscale; + double radthresh, radscale, mesh_density; int varflag, vvar, xvar, yvar, zvar; char *vstr, *xstr, *ystr, *zstr; @@ -54,6 +54,7 @@ class CreateAtoms : public Command { int *flag; // flag subset of particles to insert on lattice int *next; + int mesh_style; class Region *region; class Molecule *onemol; @@ -67,6 +68,7 @@ class CreateAtoms : public Command { void add_random(); void add_mesh(const char *); int add_tricenter(const double [3][3], tagint); + int add_quasirandom(const double [3][3], tagint); void add_lattice(); void loop_lattice(int); void add_molecule(double *);