tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

enable and apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-04-27 18:35:44 -04:00
parent 74afc04580
commit cd3ce8459b
1 changed files with 211 additions and 237 deletions
Download Patch File Download Diff File Expand all files Collapse all files

448
src/region_ellipsoid.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -21,75 +20,77 @@
#include "variable.h"
#include <cmath>
#include <limits>
using namespace LAMMPS_NS;
enum{CONSTANT,VARIABLE};
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)
Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), astr(nullptr),
bstr(nullptr), cstr(nullptr)
{
options(narg-8,&arg[8]);
options(narg - 8, &arg[8]);
if (utils::strmatch(arg[2],"^v_")) {
xstr = utils::strdup(arg[2]+2);
if (utils::strmatch(arg[2], "^v_")) {
xstr = utils::strdup(arg[2] + 2);
xc = 0.0;
xstyle = VARIABLE;
varshape = 1;
} else {
xc = xscale*utils::numeric(FLERR,arg[2],false,lmp);
xc = xscale * utils::numeric(FLERR, arg[2], false, lmp);
xstyle = CONSTANT;
}
if (utils::strmatch(arg[3],"^v_")) {
ystr = utils::strdup(arg[3]+2);
if (utils::strmatch(arg[3], "^v_")) {
ystr = utils::strdup(arg[3] + 2);
yc = 0.0;
ystyle = VARIABLE;
varshape = 1;
} else {
yc = yscale*utils::numeric(FLERR,arg[3],false,lmp);
yc = yscale * utils::numeric(FLERR, arg[3], false, lmp);
ystyle = CONSTANT;
}
if (utils::strmatch(arg[4],"^v_")) {
zstr = utils::strdup(arg[4]+2);
if (utils::strmatch(arg[4], "^v_")) {
zstr = utils::strdup(arg[4] + 2);
zc = 0.0;
zstyle = VARIABLE;
varshape = 1;
} else {
zc = zscale*utils::numeric(FLERR,arg[4],false,lmp);
zc = zscale * utils::numeric(FLERR, arg[4], false, lmp);
zstyle = CONSTANT;
}
if (utils::strmatch(arg[5],"^v_")) {
astr = utils::strdup(arg[5]+2);
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);
a = xscale * utils::numeric(FLERR, arg[5], false, lmp);
astyle = CONSTANT;
}
if (utils::strmatch(arg[6],"^v_")) {
bstr = utils::strdup(arg[6]+2);
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);
b = yscale * utils::numeric(FLERR, arg[6], false, lmp);
bstyle = CONSTANT;
}
if (utils::strmatch(arg[7],"^v_")) {
cstr = utils::strdup(arg[7]+2);
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);
c = zscale * utils::numeric(FLERR, arg[7], false, lmp);
cstyle = CONSTANT;
}
@ -100,7 +101,7 @@ RegEllipsoid::RegEllipsoid(LAMMPS *lmp, int narg, char **arg) :
// error check
if (a < 0.0 || b < 0.0 || c < 0.0) error->all(FLERR,"Illegal region ellipsoid command");
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
@ -113,7 +114,8 @@ RegEllipsoid::RegEllipsoid(LAMMPS *lmp, int narg, char **arg) :
extent_yhi = yc + b;
extent_zlo = zc - c;
extent_zhi = zc + c;
} else bboxflag = 0;
} else
bboxflag = 0;
cmax = 1;
contact = new Contact[cmax];
@ -124,13 +126,13 @@ RegEllipsoid::RegEllipsoid(LAMMPS *lmp, int narg, char **arg) :
RegEllipsoid::~RegEllipsoid()
{
delete [] xstr;
delete [] ystr;
delete [] zstr;
delete [] astr;
delete [] bstr;
delete [] cstr;
delete [] contact;
delete[] xstr;
delete[] ystr;
delete[] zstr;
delete[] astr;
delete[] bstr;
delete[] cstr;
delete[] contact;
}
/* ---------------------------------------------------------------------- */
@ -149,17 +151,17 @@ void RegEllipsoid::init()
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;
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;
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;
@ -176,91 +178,90 @@ 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;
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;
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)};
double coords[3] = {fabs(x[2] - zc), fabs(x[1] - yc), fabs(x[0] - xc)};
if (axes[1] < axes[0])
{
if (axes[1] < axes[0]) {
sorting[0] = 1;
sorting[1] = 0;
axes[0] = b;
axes[1] = c;
}
if (axes[2] < axes[1])
{
if (axes[2] < axes[1]) {
int ti = sorting[2];
sorting[2] = sorting[1];
sorting[1]= ti;
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;
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].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;
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;
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;
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].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].radius = -radius;
contact[0].iwall = 0;
contact[0].varflag = 1;
return 1;
@ -278,91 +279,90 @@ int RegEllipsoid::surface_interior(double *x, double cutoff)
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;
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;
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)};
double coords[3] = {fabs(x[2] - zc), fabs(x[1] - yc), fabs(x[0] - xc)};
if (axes[1] < axes[0])
{
if (axes[1] < axes[0]) {
sorting[0] = 1;
sorting[1] = 0;
axes[0] = b;
axes[1] = c;
}
if (axes[2] < axes[1])
{
if (axes[2] < axes[1]) {
int ti = sorting[2];
sorting[2] = sorting[1];
sorting[1]= ti;
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;
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].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;
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;
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;
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].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].radius = radius;
contact[0].iwall = 0;
contact[0].varflag = 1;
return 1;
@ -377,31 +377,21 @@ int RegEllipsoid::surface_exterior(double *x, double cutoff)
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 (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 (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 (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");
if (c < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value");
}
}
@ -413,50 +403,44 @@ 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 (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");
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 (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");
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 (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");
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 (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");
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 (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");
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 (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");
error->all(FLERR, "Variable for region ellipsoid is invalid style");
}
}
@ -476,18 +460,18 @@ void RegEllipsoid::variable_check()
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;
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 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;}
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;
g = ratio0 * ratio0 + ratio1 * ratio1 - 1;
if (g > 0) {
s0 = s;
} else if (g < 0) {
@ -499,23 +483,21 @@ double RegEllipsoid::GetRoot2D(double r0, double z0, double z1, double g)
return s;
}
double RegEllipsoid::DistancePointEllipse(double e0, double e1, double y0, double y1, double& x0, double& x1)
double RegEllipsoid::DistancePointEllipse(double e0, double e1, double y0, double y1, double &x0,
double &x1)
{
double distance;
if (y1 > 0)
{
if (y0 > 0)
{
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 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);
x0 = r0 * y0 / (sbar + r0);
x1 = y1 / (sbar + 1);
distance = sqrt((x0-y0)*(x0-y0) + (x1-y1)*(x1-y1));
distance = sqrt((x0 - y0) * (x0 - y0) + (x1 - y1) * (x1 - y1));
} else {
x0 = y0;
x1 = y1;
@ -527,18 +509,17 @@ double RegEllipsoid::DistancePointEllipse(double e0, double e1, double y0, doubl
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);
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);
distance = fabs(y0 - e0);
}
}
return distance;
@ -550,20 +531,20 @@ double RegEllipsoid::DistancePointEllipse(double e0, double e1, double y0, doubl
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;
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 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;
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) {
@ -575,28 +556,25 @@ double RegEllipsoid::GetRoot3D(double r0, double r1, double z0, double z1, doubl
return s;
}
double RegEllipsoid::DistancePointEllipsoid(double e0, double e1, double e2, double y0, double y1, double y2, double& x0, double& x1, double& x2)
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);
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));
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;
@ -608,8 +586,7 @@ double RegEllipsoid::DistancePointEllipsoid(double e0, double e1, double e2, dou
distance = DistancePointEllipse(e1, e2, y1, y2, x1, x2);
}
} else {
if (y0 > 0)
{
if (y0 > 0) {
x1 = 0;
distance = DistancePointEllipse(e0, e2, y0, y2, x0, x2);
} else {
@ -620,29 +597,26 @@ double RegEllipsoid::DistancePointEllipsoid(double e0, double e1, double e2, dou
}
}
} else {
double denom0 = e0*e0 - e2*e2;
double denom1 = e1*e1 - e2*e2;
double numer0 = e0*y0;
double numer1 = e1*y1;
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;
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);
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)
{
if (!computed) {
x2 = 0;
distance = DistancePointEllipse(e0, e1, y0, y1, x0, x1);
}