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

Merge pull request #3553 from evoyiatzis/patch-3

Browse Source 
Fixing bug #3545
This commit is contained in:
Axel Kohlmeyer
2022-12-21 12:05:44 -05:00
committed by GitHub
parent fa6251d83b 12a23b0ef1
commit d8620fc34c
1 changed files with 52 additions and 64 deletions
Download Patch File Download Diff File Expand all files Collapse all files

116
src/region_ellipsoid.cpp
View File

@ -197,39 +197,33 @@ int RegEllipsoid::surface_interior(double *x, double cutoff)
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 axes[3] = {a,b,c};
double coords[3] = {fabs(x[0] - xc), fabs(x[1] - yc), fabs(x[2] - zc)};
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;
}
int min, max;
if (axes[0] < axes[1]) {
min = 0;
max = 1;
} else {
min = 1;
max = 0;
}
if (axes[min] > axes[2]) {
min = 2;
}
if (axes[max] < axes[2]) {
max = 2;
}
int mid = 3 - min - max;
int sorting[3] = {min, mid, max};
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;
DistancePointEllipsoid(axes[sorting[2]], axes[sorting[1]], axes[sorting[0]], coords[sorting[2]], coords[sorting[1]],
coords[sorting[0]], x0[sorting[2]], x0[sorting[1]], x0[sorting[0]]);
contact[0].delx = x[0] - (copysign(x0[sorting[2]], x[0] - xc) + xc);
contact[0].dely = x[1] - (copysign(x0[sorting[1]], x[1] - yc) + yc);
contact[0].delz = x[2] - (copysign(x0[sorting[0]], x[2] - zc) + zc);
// contact[0].radius = -radius;
contact[0].iwall = 0;
contact[0].varflag = 1;
@ -254,12 +248,12 @@ int RegEllipsoid::surface_interior(double *x, double cutoff)
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].delx = x[0] - (copysign(x0, x[0] - xc) + xc);
contact[0].dely = x[1] - (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].delx = x[0] - (copysign(x1, x[0] - xc) + xc);
contact[0].dely = x[1] - (copysign(x0, x[1] - yc) + yc);
}
contact[0].delz = 0;
// contact[0].radius = -radius;
@ -298,39 +292,33 @@ int RegEllipsoid::surface_exterior(double *x, double cutoff)
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 axes[3] = {a,b,c};
double coords[3] = {fabs(x[0] - xc), fabs(x[1] - yc), fabs(x[2] - zc)};
if (axes[1] < axes[0]) {
sorting[0] = 1;
sorting[1] = 0;
axes[0] = b;
axes[1] = c;
int min, max;
if (axes[0] < axes[1]) {
min = 0;
max = 1;
} else {
min = 1;
max = 0;
}
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;
if (axes[min] > axes[2]) {
min = 2;
}
if (axes[max] < axes[2]) {
max = 2;
}
int mid = 3 - min - max;
int sorting[3] = {min, mid, max};
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;
DistancePointEllipsoid(axes[sorting[2]], axes[sorting[1]], axes[sorting[0]], coords[sorting[2]], coords[sorting[1]],
coords[sorting[0]], x0[sorting[2]], x0[sorting[1]], x0[sorting[0]]);
contact[0].delx = x[0] - (copysign(x0[sorting[2]], x[0] - xc) + xc);
contact[0].dely = x[1] - (copysign(x0[sorting[1]], x[1] - yc) + yc);
contact[0].delz = x[2] - (copysign(x0[sorting[0]], x[2] - zc) + zc);
// contact[0].radius = radius;
contact[0].iwall = 0;
contact[0].varflag = 1;
@ -355,12 +343,12 @@ int RegEllipsoid::surface_exterior(double *x, double cutoff)
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].delx = x[0] - (copysign(x0, x[0] - xc) + xc);
contact[0].dely = x[1] - (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].delx = x[0] - (copysign(x1, x[0] - xc) + xc);
contact[0].dely = x[1] - (copysign(x0, x[1] - yc) + yc);
}
contact[0].delz = 0;
// contact[0].radius = radius;