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

Merge pull request #2871 from evoyiatzis/master

Browse Source 
implementation of an "ellipsoidal" region option
This commit is contained in:
Axel Kohlmeyer
2022-04-29 19:46:15 -04:00
committed by GitHub
parent 44c4d0a4d6 a05bdfa76c
commit 87b0939fe7
3 changed files with 716 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

41
doc/src/region.rst
View File

@ -11,7 +11,7 @@ Syntax
region ID style args keyword arg ...
* ID = user-assigned name for the region
* style = *delete* or *block* or *cone* or *cylinder* or *plane* or *prism* or *sphere* or *union* or *intersect*
* style = *delete* or *block* or *cone* or *cylinder* or *ellipsoid* or *plane* or *prism* or *sphere* or *union* or *intersect*
.. parsed-literal::
@ -29,6 +29,10 @@ Syntax
radius = cylinder radius (distance units)
c1,c2, and radius can be a variable (see below)
lo,hi = bounds of cylinder in dim (distance units)
*ellipsoid* args = x y z a b c
x,y,z = center of ellipsoid (distance units)
a,b,c = half the length of the principal axes of the ellipsoid (distance units)
x,y,z,a,b and c can be a variable (see below)
*plane* args = px py pz nx ny nz
px,py,pz = point on the plane (distance units)
nx,ny,nz = direction normal to plane (distance units)
@ -60,7 +64,7 @@ Syntax
*lattice* = the geometry is defined in lattice units
*box* = the geometry is defined in simulation box units
*move* args = v_x v_y v_z
v_x,v_y,v_z = equal-style variables for x,y,z displacement of region over time
v_x,v_y,v_z = equal-style variables for x,y,z displacement of region over time (distance units)
*rotate* args = v_theta Px Py Pz Rx Ry Rz
v_theta = equal-style variable for rotaton of region over time (in radians)
Px,Py,Pz = origin for axis of rotation (distance units)
@ -158,6 +162,12 @@ Thus the third example above specifies a cylinder with its axis in the
y-direction located at x = 2.0 and z = 3.0, with a radius of 5.0, and
extending in the y-direction from -5.0 to the upper box boundary.
For style *ellipsoid*, an axis-aligned ellipsoid is defined. The
ellipsoid has its center at (x,y,z) and is defined by 3 axis-aligned
vectors given by A = (a,0,0); B = (0,b,0); C = (0,0,c). Note that
although the ellipsoid is specified as axis-aligned it can be rotated
via the optional *rotate* keyword.
For style *plane*, a plane is defined which contain the point
(px,py,pz) and has a normal vector (nx,ny,nz). The normal vector does
not have to be of unit length. The "inside" of the plane is the
@ -184,15 +194,21 @@ since if the maximum tilt factor is 5 (as in this example), then
configurations with tilt = ..., -15, -5, 5, 15, 25, ... are all
geometrically equivalent.
The *radius* value for style *sphere* and *cylinder* can be specified
as an equal-style :doc:`variable <variable>`. If the value 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 determine the radius of the region. For style *sphere*
also the x-, y-, and z- coordinate of the center of the sphere and for
style *cylinder* the two center positions c1 and c2 for the location of
the cylinder axes can be a variable with the same kind of effect and
requirements than for the radius.
For style *sphere*, a sphere is defined with its center at (x,y,z)
and with radius as its radius.
The *radius* value for styles *sphere* and *cylinder*, and the
parameters a,b,c for style *ellipsoid*, can each be specified as an
equal-style :doc:`variable <variable>`. Likewise, for style *sphere*
and *ellipsoid* the x-, y-, and z- coordinates of the center of the
sphere/ellipsoid can be specified as an equal-style variable. And for
style *cylinder* the two center positions c1 and c2 for the location
of the cylinder axes can be specified as a equal-style variable.
If the value 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 determine the radius of
the region.
Equal-style variables can specify formulas with various mathematical
functions, and include :doc:`thermo_style <thermo_style>` command
@ -250,6 +266,9 @@ define the lattice spacings which are used as follows:
to lo and hi. The spacings in the two radial dimensions are applied
to c1 and c2. The cylinder radius is scaled by the lattice
spacing in the dimension corresponding to c1.
* For style *ellipsoid*, the lattice spacing in dimensions x,y,z are
applied to the ellipsoid center x,y,z. The spacing in dimensions
x,y,z are applied to the ellipsoid radii a,b,c respectively.
* For style *plane*, the lattice spacing in dimension x is applied to
px and nx, similarly the spacings in dimensions y,z are applied to
py/ny and pz/nz.

625
src/region_ellipsoid.cpp Normal file
View File

@ -0,0 +1,625 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "region_ellipsoid.h"
#include "domain.h"
#include "error.h"
#include "input.h"
#include "update.h"
#include "variable.h"
#include <cmath>
#include <limits>
using namespace LAMMPS_NS;
enum { CONSTANT, VARIABLE };
/* ---------------------------------------------------------------------- */
RegEllipsoid::RegEllipsoid(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), astr(nullptr),
bstr(nullptr), cstr(nullptr)
{
options(narg - 8, &arg[8]);
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);
xstyle = CONSTANT;
}
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);
ystyle = CONSTANT;
}
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);
zstyle = CONSTANT;
}
if (utils::strmatch(arg[5], "^v_")) {
astr = utils::strdup(arg[5] + 2);
a = 0.0;
astyle = VARIABLE;
varshape = 1;
} else {
a = xscale * utils::numeric(FLERR, arg[5], false, lmp);
astyle = CONSTANT;
}
if (utils::strmatch(arg[6], "^v_")) {
bstr = utils::strdup(arg[6] + 2);
b = 0.0;
bstyle = VARIABLE;
varshape = 1;
} else {
b = yscale * utils::numeric(FLERR, arg[6], false, lmp);
bstyle = CONSTANT;
}
if (utils::strmatch(arg[7], "^v_")) {
cstr = utils::strdup(arg[7] + 2);
c = 0.0;
cstyle = VARIABLE;
varshape = 1;
} else {
c = zscale * utils::numeric(FLERR, arg[7], false, lmp);
cstyle = CONSTANT;
}
if (varshape) {
variable_check();
shape_update();
}
// error check
if (a < 0.0 || b < 0.0 || c < 0.0) error->all(FLERR, "Illegal region ellipsoid command");
// extent of ellipsoid
// for variable axes, uses initial axes and origin for variable center
if (interior) {
bboxflag = 1;
extent_xlo = xc - a;
extent_xhi = xc + a;
extent_ylo = yc - b;
extent_yhi = yc + b;
extent_zlo = zc - c;
extent_zhi = zc + c;
} else
bboxflag = 0;
cmax = 1;
contact = new Contact[cmax];
tmax = 1;
}
/* ---------------------------------------------------------------------- */
RegEllipsoid::~RegEllipsoid()
{
delete[] xstr;
delete[] ystr;
delete[] zstr;
delete[] astr;
delete[] bstr;
delete[] cstr;
delete[] contact;
}
/* ---------------------------------------------------------------------- */
void RegEllipsoid::init()
{
Region::init();
if (varshape) variable_check();
}
/* ----------------------------------------------------------------------
inside = 1 if x,y,z is inside or on surface
inside = 0 if x,y,z is outside and not on surface
------------------------------------------------------------------------- */
int RegEllipsoid::inside(double x, double y, double z)
{
if (domain->dimension == 3) {
double delx = b * c * (x - xc);
double dely = a * c * (y - yc);
double delz = a * b * (z - zc);
double r = delx * delx + dely * dely + delz * delz;
double rc = a * a * b * b * c * c;
if (r <= rc) return 1;
} else {
double delx = b * (x - xc);
double dely = a * (y - yc);
double r = delx * delx + dely * dely;
double rc = a * a * b * b;
if (r <= rc) return 1;
}
return 0;
}
/* ----------------------------------------------------------------------
one contact if 0 <= x < cutoff from inner surface of ellipsoid
no contact if outside (possible if called from union/intersect)
delxyz = vector from nearest point on ellipsoid to x
special case: no contact if x is at center of ellipsoid
------------------------------------------------------------------------- */
int RegEllipsoid::surface_interior(double *x, double cutoff)
{
if (domain->dimension == 3) {
double delx = b * c * (x[0] - xc);
double dely = a * c * (x[1] - yc);
double delz = a * b * (x[2] - zc);
double r = delx * delx + dely * dely + delz * delz;
double rc = a * a * b * b * c * c;
if (r > rc || r == 0.0) return 0;
double a_r = a - cutoff;
double b_r = b - cutoff;
double c_r = c - cutoff;
double delx_r = b_r * c_r * (x[0] - xc);
double dely_r = a_r * c_r * (x[1] - xc);
double delz_r = a_r * b_r * (x[2] - xc);
double r_r = delx_r * delx_r + dely_r * dely_r + delz_r * delz_r;
double rc_r = a_r * a_r * b_r * b_r * c_r * c_r;
if (r_r > rc_r) {
// sort the values
int sorting[3] = {0, 1, 2};
double axes[3] = {c, b, a};
double coords[3] = {fabs(x[2] - zc), fabs(x[1] - yc), fabs(x[0] - xc)};
if (axes[1] < axes[0]) {
sorting[0] = 1;
sorting[1] = 0;
axes[0] = b;
axes[1] = c;
}
if (axes[2] < axes[1]) {
int ti = sorting[2];
sorting[2] = sorting[1];
sorting[1] = ti;
double td = axes[2];
axes[2] = axes[1];
axes[1] = td;
if (axes[1] < axes[0]) ti = sorting[1];
sorting[1] = sorting[0];
sorting[0] = ti;
td = axes[1];
axes[1] = axes[0];
axes[0] = td;
}
double x0[3];
contact[0].r =
DistancePointEllipsoid(axes[2], axes[1], axes[0], coords[sorting[2]], coords[sorting[1]],
coords[sorting[0]], x0[2], x0[1], x0[0]);
contact[0].delx = copysign(x0[sorting[2]], x[0] - xc) + xc;
contact[0].dely = copysign(x0[sorting[1]], x[1] - yc) + yc;
contact[0].delz = copysign(x0[sorting[0]], x[2] - zc) + zc;
// contact[0].radius = -radius;
contact[0].iwall = 0;
contact[0].varflag = 1;
return 1;
}
return 0;
} else {
double delx = b * (x[0] - xc);
double dely = a * (x[1] - yc);
double r = delx * delx + dely * dely;
double rc = a * a * b * b;
if (r > rc || r == 0.0) return 0;
double a_r = a - cutoff;
double b_r = b - cutoff;
double delx_r = b_r * (x[0] - xc);
double dely_r = a_r * (x[1] - xc);
double r_r = delx_r * delx_r + dely_r * dely_r;
double rc_r = a_r * a_r * b_r * b_r;
if (r_r > rc_r) {
double x0, x1;
if (a >= b) {
contact[0].r = DistancePointEllipse(a, b, fabs(x[0] - xc), fabs(x[1] - yc), x0, x1);
contact[0].delx = copysign(x0, x[0] - xc) + xc;
contact[0].dely = copysign(x1, x[1] - yc) + yc;
} else {
contact[0].r = DistancePointEllipse(b, a, fabs(x[1] - yc), fabs(x[0] - xc), x0, x1);
contact[0].delx = copysign(x1, x[0] - xc) + xc;
contact[0].dely = copysign(x0, x[1] - yc) + yc;
}
contact[0].delz = 0;
// contact[0].radius = -radius;
contact[0].iwall = 0;
contact[0].varflag = 1;
return 1;
}
return 0;
}
}
/* ----------------------------------------------------------------------
one contact if 0 <= x < cutoff from outer surface of ellipsoid
no contact if inside (possible if called from union/intersect)
delxyz = vector from nearest point on ellipsoid to x
------------------------------------------------------------------------- */
int RegEllipsoid::surface_exterior(double *x, double cutoff)
{
if (domain->dimension == 3) {
double delx = b * c * (x[0] - xc);
double dely = a * c * (x[1] - yc);
double delz = a * b * (x[2] - zc);
double r = delx * delx + dely * dely + delz * delz;
double rc = a * a * b * b * c * c;
if (r < rc) return 0;
double a_r = a + cutoff;
double b_r = b + cutoff;
double c_r = c + cutoff;
double delx_r = b_r * c_r * (x[0] - xc);
double dely_r = a_r * c_r * (x[1] - xc);
double delz_r = a_r * b_r * (x[2] - xc);
double r_r = delx_r * delx_r + dely_r * dely_r + delz_r * delz_r;
double rc_r = a_r * a_r * b_r * b_r * c_r * c_r;
if (r_r < rc_r) {
// sort the values
int sorting[3] = {0, 1, 2};
double axes[3] = {c, b, a};
double coords[3] = {fabs(x[2] - zc), fabs(x[1] - yc), fabs(x[0] - xc)};
if (axes[1] < axes[0]) {
sorting[0] = 1;
sorting[1] = 0;
axes[0] = b;
axes[1] = c;
}
if (axes[2] < axes[1]) {
int ti = sorting[2];
sorting[2] = sorting[1];
sorting[1] = ti;
double td = axes[2];
axes[2] = axes[1];
axes[1] = td;
if (axes[1] < axes[0]) ti = sorting[1];
sorting[1] = sorting[0];
sorting[0] = ti;
td = axes[1];
axes[1] = axes[0];
axes[0] = td;
}
double x0[3];
contact[0].r =
DistancePointEllipsoid(axes[2], axes[1], axes[0], coords[sorting[2]], coords[sorting[1]],
coords[sorting[0]], x0[2], x0[1], x0[0]);
contact[0].delx = copysign(x0[sorting[2]], x[0] - xc) + xc;
contact[0].dely = copysign(x0[sorting[1]], x[1] - yc) + yc;
contact[0].delz = copysign(x0[sorting[0]], x[2] - zc) + zc;
// contact[0].radius = radius;
contact[0].iwall = 0;
contact[0].varflag = 1;
return 1;
}
return 0;
} else {
double delx = b * c * (x[0] - xc);
double dely = a * c * (x[1] - yc);
double r = delx * delx + dely * dely;
double rc = a * a * b * b;
if (r < rc) return 0;
double a_r = a + cutoff;
double b_r = b + cutoff;
double delx_r = b_r * (x[0] - xc);
double dely_r = a_r * (x[1] - xc);
double r_r = delx_r * delx_r + dely_r * dely_r;
double rc_r = a_r * a_r * b_r * b_r;
if (r_r < rc_r) {
double x0, x1;
if (a >= b) {
contact[0].r = DistancePointEllipse(a, b, fabs(x[0] - xc), fabs(x[1] - yc), x0, x1);
contact[0].delx = copysign(x0, x[0] - xc) + xc;
contact[0].dely = copysign(x1, x[1] - yc) + yc;
} else {
contact[0].r = DistancePointEllipse(b, a, fabs(x[1] - yc), fabs(x[0] - xc), x0, x1);
contact[0].delx = copysign(x1, x[0] - xc) + xc;
contact[0].dely = copysign(x0, x[1] - yc) + yc;
}
contact[0].delz = 0;
// contact[0].radius = radius;
contact[0].iwall = 0;
contact[0].varflag = 1;
return 1;
}
return 0;
}
}
/* ----------------------------------------------------------------------
change region shape via variable evaluation
------------------------------------------------------------------------- */
void RegEllipsoid::shape_update()
{
if (xstyle == VARIABLE) xc = xscale * input->variable->compute_equal(xvar);
if (ystyle == VARIABLE) yc = yscale * input->variable->compute_equal(yvar);
if (zstyle == VARIABLE) zc = zscale * input->variable->compute_equal(zvar);
if (astyle == VARIABLE) {
a = xscale * input->variable->compute_equal(avar);
if (a < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
}
if (bstyle == VARIABLE) {
b = yscale * input->variable->compute_equal(bvar);
if (b < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
}
if (cstyle == VARIABLE) {
c = zscale * input->variable->compute_equal(cvar);
if (c < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
}
}
/* ----------------------------------------------------------------------
error check on existence of variable
------------------------------------------------------------------------- */
void RegEllipsoid::variable_check()
{
if (xstyle == VARIABLE) {
xvar = input->variable->find(xstr);
if (xvar < 0) error->all(FLERR, "Variable name for region ellipsoid does not exist");
if (!input->variable->equalstyle(xvar))
error->all(FLERR, "Variable for region ellipsoid is invalid style");
}
if (ystyle == VARIABLE) {
yvar = input->variable->find(ystr);
if (yvar < 0) error->all(FLERR, "Variable name for region ellipsoid does not exist");
if (!input->variable->equalstyle(yvar))
error->all(FLERR, "Variable for region ellipsoid is invalid style");
}
if (zstyle == VARIABLE) {
zvar = input->variable->find(zstr);
if (zvar < 0) error->all(FLERR, "Variable name for region ellipsoid does not exist");
if (!input->variable->equalstyle(zvar))
error->all(FLERR, "Variable for region ellipsoid is invalid style");
}
if (astyle == VARIABLE) {
avar = input->variable->find(astr);
if (avar < 0) error->all(FLERR, "Variable name for region ellipsoid does not exist");
if (!input->variable->equalstyle(avar))
error->all(FLERR, "Variable for region ellipsoid is invalid style");
}
if (bstyle == VARIABLE) {
bvar = input->variable->find(bstr);
if (bvar < 0) error->all(FLERR, "Variable name for region ellipsoid does not exist");
if (!input->variable->equalstyle(bvar))
error->all(FLERR, "Variable for region ellipsoid is invalid style");
}
if (cstyle == VARIABLE) {
cvar = input->variable->find(cstr);
if (cvar < 0) error->all(FLERR, "Variable name for region ellipsoid does not exist");
if (!input->variable->equalstyle(cvar))
error->all(FLERR, "Variable for region ellipsoid is invalid style");
}
}
// ------------------------------------------------------------------
// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2021
// Based on https://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// Version: 4.0.2021.08.01
// ------------------------------------------------------------------
/* ----------------------------------------------------------------------
functions for the 2D case
------------------------------------------------------------------------- */
double RegEllipsoid::GetRoot2D(double r0, double z0, double z1, double g)
{
int maxIterations =
std::numeric_limits<double>::digits - std::numeric_limits<double>::min_exponent;
double n0 = r0 * z0;
double s0 = z1 - 1;
double s1 = (g < 0 ? 0 : sqrt(n0 * n0 + z1 * z1) - 1);
double s = 0;
for (int i = 0; i < maxIterations; ++i) {
s = (s0 + s1) / 2;
if (s == s0 || s == s1) { break; }
double ratio0 = n0 / (s + r0);
double ratio1 = z1 / (s + 1);
g = ratio0 * ratio0 + ratio1 * ratio1 - 1;
if (g > 0) {
s0 = s;
} else if (g < 0) {
s1 = s;
} else {
break;
}
}
return s;
}
double RegEllipsoid::DistancePointEllipse(double e0, double e1, double y0, double y1, double &x0,
double &x1)
{
double distance;
if (y1 > 0) {
if (y0 > 0) {
double z0 = y0 / e0;
double z1 = y1 / e1;
double g = z0 * z0 + z1 * z1 - 1;
if (g != 0) {
double r0 = (e0 * e0) / (e1 * e1);
double sbar = GetRoot2D(r0, z0, z1, g);
x0 = r0 * y0 / (sbar + r0);
x1 = y1 / (sbar + 1);
distance = sqrt((x0 - y0) * (x0 - y0) + (x1 - y1) * (x1 - y1));
} else {
x0 = y0;
x1 = y1;
distance = 0;
}
} else {
x0 = 0;
x1 = e1;
distance = fabs(y1 - e1);
}
} else {
double numer0 = e0 * y0;
double denom0 = e0 * e0 - e1 * e1;
if (numer0 < denom0) {
double xde0 = numer0 / denom0;
x0 = e0 * xde0;
x1 = e1 * sqrt(1 - xde0 * xde0);
distance = sqrt((x0 - y0) * (x0 - y0) + x1 * x1);
} else {
x0 = e0;
x1 = 0;
distance = fabs(y0 - e0);
}
}
return distance;
}
/* ----------------------------------------------------------------------
functions for the 3D case
------------------------------------------------------------------------- */
double RegEllipsoid::GetRoot3D(double r0, double r1, double z0, double z1, double z2, double g)
{
int maxIterations =
std::numeric_limits<double>::digits - std::numeric_limits<double>::min_exponent;
double n0 = r0 * z0;
double n1 = r1 * z1;
double s0 = z2 - 1;
double s1 = (g < 0 ? 0 : sqrt(n0 * n0 + n1 * n1 + z2 * z2) - 1);
double s = 0;
for (int i = 0; i < maxIterations; ++i) {
s = (s0 + s1) / 2;
if (s == s0 || s == s1) { break; }
double ratio0 = n0 / (s + r0);
double ratio1 = n1 / (s + r1);
double ratio2 = z2 / (s + 1);
g = ratio0 * ratio0 + ratio1 * ratio1 + ratio2 * ratio2 - 1;
if (g > 0) {
s0 = s;
} else if (g < 0) {
s1 = s;
} else {
break;
}
}
return s;
}
double RegEllipsoid::DistancePointEllipsoid(double e0, double e1, double e2, double y0, double y1,
double y2, double &x0, double &x1, double &x2)
{
double distance;
if (y2 > 0) {
if (y1 > 0) {
if (y0 > 0) {
double z0 = y0 / e0;
double z1 = y1 / e1;
double z2 = y2 / e2;
double g = z0 * z0 + z1 * z1 + z2 * z2 - 1;
if (g != 0) {
double r0 = e0 * e0 / (e2 * e2);
double r1 = e1 * e1 / (e2 * e2);
double sbar = GetRoot3D(r0, r1, z0, z1, z2, g);
x0 = r0 * y0 / (sbar + r0);
x1 = r1 * y1 / (sbar + r1);
x2 = y2 / (sbar + 1);
distance = sqrt((x0 - y0) * (x0 - y0) + (x1 - y1) * (x1 - y1) + (x2 - y2) * (x2 - y2));
} else {
x0 = y0;
x1 = y1;
x2 = y2;
distance = 0;
}
} else {
x0 = 0;
distance = DistancePointEllipse(e1, e2, y1, y2, x1, x2);
}
} else {
if (y0 > 0) {
x1 = 0;
distance = DistancePointEllipse(e0, e2, y0, y2, x0, x2);
} else {
x0 = 0;
x1 = 0;
x2 = e2;
distance = fabs(y2 - e2);
}
}
} else {
double denom0 = e0 * e0 - e2 * e2;
double denom1 = e1 * e1 - e2 * e2;
double numer0 = e0 * y0;
double numer1 = e1 * y1;
bool computed = false;
if (numer0 < denom0 && numer1 < denom1) {
double xde0 = numer0 / denom0;
double xde1 = numer1 / denom1;
double xde0sqr = xde0 * xde0;
double xde1sqr = xde1 * xde1;
double discr = 1 - xde0sqr - xde1sqr;
if (discr > 0) {
x0 = e0 * xde0;
x1 = e1 * xde1;
x2 = e2 * sqrt(discr);
distance = sqrt((x0 - y0) * (x0 - y0) + (x1 - y1) * (x1 - y1) + x2 * x2);
computed = true;
}
}
if (!computed) {
x2 = 0;
distance = DistancePointEllipse(e0, e1, y0, y1, x0, x1);
}
}
return distance;
}

61
src/region_ellipsoid.h Normal file
View File

@ -0,0 +1,61 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef REGION_CLASS
// clang-format off
RegionStyle(ellipsoid,RegEllipsoid);
// clang-format on
#else
#ifndef LMP_REGION_ELLIPSOID_H
#define LMP_REGION_ELLIPSOID_H
#include "region.h"
namespace LAMMPS_NS {
class RegEllipsoid : public Region {
public:
RegEllipsoid(class LAMMPS *, int, char **);
~RegEllipsoid() override;
void init() override;
int inside(double, double, double) override;
int surface_interior(double *, double) override;
int surface_exterior(double *, double) override;
void shape_update() override;
private:
double xc, yc, zc;
double a, b, c;
int xstyle, xvar;
int ystyle, yvar;
int zstyle, zvar;
int astyle, avar;
int bstyle, bvar;
int cstyle, cvar;
char *xstr, *ystr, *zstr;
char *astr, *bstr, *cstr;
void variable_check();
double GetRoot2D(double r0, double z0, double z1, double g);
double GetRoot3D(double r0, double r1, double z0, double z1, double z2, double g);
double DistancePointEllipse(double e0, double e1, double y0, double y1, double &x0, double &x1);
double DistancePointEllipsoid(double e0, double e1, double e2, double y0, double y1, double y2,
double &x0, double &x1, double &x2);
};
} // namespace LAMMPS_NS
#endif
#endif