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

new fix wall/gran/region command, REBO bug fix, new example log files

Browse Source
This commit is contained in:
Steve Plimpton
2016-10-06 15:47:41 -06:00
parent c213457550
commit b3d2fb91bb
63 changed files with 6354 additions and 490 deletions
Download Patch File Download Diff File Expand all files Collapse all files

207
src/region_cone.cpp
View File

@ -34,6 +34,11 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) :
{
options(narg-9,&arg[9]);
// check open face settings
if (openflag && (open_faces[3] || open_faces[4] || open_faces[5]))
error->all(FLERR,"Invalid region cone open setting");
if (strcmp(arg[2],"x") && strcmp(arg[2],"y") && strcmp(arg[2],"z"))
error->all(FLERR,"Illegal region cylinder command");
axis = arg[2][0];
@ -118,7 +123,6 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) :
if (interior) {
bboxflag = 1;
if (axis == 'x') {
extent_xlo = lo;
extent_xhi = hi;
@ -145,10 +149,13 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) :
}
} else bboxflag = 0;
// particle could be contact cone surface and 2 ends
// particle could be close to cone surface and 2 ends
// particle can only touch surface and 1 end
cmax = 3;
contact = new Contact[cmax];
if (interior) tmax = 2;
else tmax = 1;
}
/* ---------------------------------------------------------------------- */
@ -167,7 +174,7 @@ int RegCone::inside(double x, double y, double z)
{
double del1,del2,dist;
double currentradius;
int inside = 0;
int inside;
if (axis == 'x') {
del1 = y - c1;
@ -226,7 +233,7 @@ int RegCone::surface_interior(double *x, double cutoff)
// surflo = pt on outer circle of bottom end plane, same dir as x vs axis
// surfhi = pt on outer circle of top end plane, same dir as x vs axis
if (r > 0.0) {
if (r > 0.0 && !open_faces[2]) {
surflo[0] = lo;
surflo[1] = c1 + del1*radiuslo/r;
surflo[2] = c2 + del2*radiuslo/r;
@ -243,22 +250,29 @@ int RegCone::surface_interior(double *x, double cutoff)
contact[n].delx = delx;
contact[n].dely = dely;
contact[n].delz = delz;
contact[n].radius = -2.0*(radiuslo + (xs[0]-lo)*
(radiushi-radiuslo)/(hi-lo));
contact[n].iwall = 2;
n++;
}
}
delta = x[0] - lo;
if (delta < cutoff) {
if (delta < cutoff && !open_faces[0]) {
contact[n].r = delta;
contact[n].delx = delta;
contact[n].dely = contact[n].delz = 0.0;
contact[n].radius = 0;
contact[n].iwall = 0;
n++;
}
delta = hi - x[0];
if (delta < cutoff) {
if (delta < cutoff && !open_faces[1]) {
contact[n].r = delta;
contact[n].delx = -delta;
contact[n].dely = contact[n].delz = 0.0;
contact[n].radius = 0;
contact[n].iwall = 1;
n++;
}
@ -266,7 +280,8 @@ int RegCone::surface_interior(double *x, double cutoff)
del1 = x[0] - c1;
del2 = x[2] - c2;
r = sqrt(del1*del1 + del2*del2);
currentradius = radiuslo + (x[1]-lo)*(radiushi-radiuslo)/(hi-lo);
currentradius = radiuslo + (x[1]-lo)*
(radiushi-radiuslo)/(hi-lo);
// y is exterior to cone
@ -276,7 +291,7 @@ int RegCone::surface_interior(double *x, double cutoff)
// surflo = pt on outer circle of bottom end plane, same dir as y vs axis
// surfhi = pt on outer circle of top end plane, same dir as y vs axis
if (r > 0.0) {
if (r > 0.0 && !open_faces[2]) {
surflo[0] = c1 + del1*radiuslo/r;
surflo[1] = lo;
surflo[2] = c2 + del2*radiuslo/r;
@ -293,22 +308,29 @@ int RegCone::surface_interior(double *x, double cutoff)
contact[n].delx = delx;
contact[n].dely = dely;
contact[n].delz = delz;
contact[n].iwall = 2;
contact[n].radius = -2.0*(radiuslo + (xs[1]-lo)*
(radiushi-radiuslo)/(hi-lo));
n++;
}
}
delta = x[1] - lo;
if (delta < cutoff) {
if (delta < cutoff && !open_faces[0]) {
contact[n].r = delta;
contact[n].delz = delta;
contact[n].delx = contact[n].dely = 0.0;
contact[n].iwall = 0;
contact[n].radius = 0;
n++;
}
delta = hi - x[1];
if (delta < cutoff) {
if (delta < cutoff && !open_faces[1]) {
contact[n].r = delta;
contact[n].delz = -delta;
contact[n].delx = contact[n].dely = 0.0;
contact[n].iwall = 1;
contact[n].radius = 0;
n++;
}
@ -326,7 +348,7 @@ int RegCone::surface_interior(double *x, double cutoff)
// surflo = pt on outer circle of bottom end plane, same dir as z vs axis
// surfhi = pt on outer circle of top end plane, same dir as z vs axis
if (r > 0.0) {
if (r > 0.0 && !open_faces[2]) {
surflo[0] = c1 + del1*radiuslo/r;
surflo[1] = c2 + del2*radiuslo/r;
surflo[2] = lo;
@ -343,22 +365,29 @@ int RegCone::surface_interior(double *x, double cutoff)
contact[n].delx = delx;
contact[n].dely = dely;
contact[n].delz = delz;
contact[n].iwall = 2;
contact[n].radius = -2.0*(radiuslo + (xs[2]-lo)*
(radiushi-radiuslo)/(hi-lo));
n++;
}
}
delta = x[2] - lo;
if (delta < cutoff) {
if (delta < cutoff && !open_faces[0]) {
contact[n].r = delta;
contact[n].delz = delta;
contact[n].delx = contact[n].dely = 0.0;
contact[n].iwall = 0;
contact[n].radius = 0;
n++;
}
delta = hi - x[2];
if (delta < cutoff) {
if (delta < cutoff && !open_faces[1]) {
contact[n].r = delta;
contact[n].delz = -delta;
contact[n].delx = contact[n].dely = 0.0;
contact[n].iwall = 1;
contact[n].radius = 0;
n++;
}
}
@ -374,7 +403,7 @@ int RegCone::surface_interior(double *x, double cutoff)
int RegCone::surface_exterior(double *x, double cutoff)
{
double del1,del2,r,currentradius,distsq;
double del1,del2,r,currentradius,distsq,distsqprev,crad;
double corner1[3],corner2[3],corner3[3],corner4[3],xp[3],nearest[3];
if (axis == 'x') {
@ -383,11 +412,15 @@ int RegCone::surface_exterior(double *x, double cutoff)
r = sqrt(del1*del1 + del2*del2);
currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo);
// radius of curvature, only used for granular walls
double crad = 0.0;
// x is far enough from cone that there is no contact
// x is interior to cone
if (r >= maxradius+cutoff ||
x[0] <= lo-cutoff || x[0] >= hi+cutoff) return 0;
if (r >= maxradius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff)
return 0;
if (r < currentradius && x[0] > lo && x[0] < hi) return 0;
// x is exterior to cone or on its surface
@ -411,14 +444,31 @@ int RegCone::surface_exterior(double *x, double cutoff)
corner4[2] = c2;
distsq = BIG;
point_on_line_segment(corner1,corner2,x,xp);
distsq = closest(x,xp,nearest,distsq);
point_on_line_segment(corner1,corner3,x,xp);
distsq = closest(x,xp,nearest,distsq);
point_on_line_segment(corner2,corner4,x,xp);
distsq = closest(x,xp,nearest,distsq);
if (!open_faces[2]) {
point_on_line_segment(corner1,corner2,x,xp);
distsq = closest(x,xp,nearest,distsq);
crad = -2.0*(radiuslo + (nearest[0]-lo)*(radiushi-radiuslo)/(hi-lo));
}
if (!open_faces[0]) {
point_on_line_segment(corner1,corner3,x,xp);
distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq);
if (distsq < distsqprev) crad = 0.0;
}
if (!open_faces[1]) {
point_on_line_segment(corner2,corner4,x,xp);
distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq);
if (distsq < distsqprev) crad = 0.0;
}
if (distsq == BIG) return 0;
add_contact(0,x,nearest[0],nearest[1],nearest[2]);
contact[0].radius = crad;
contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1;
return 0;
@ -456,14 +506,30 @@ int RegCone::surface_exterior(double *x, double cutoff)
corner4[2] = c2;
distsq = BIG;
point_on_line_segment(corner1,corner2,x,xp);
distsq = closest(x,xp,nearest,distsq);
point_on_line_segment(corner1,corner3,x,xp);
distsq = closest(x,xp,nearest,distsq);
point_on_line_segment(corner2,corner4,x,xp);
distsq = closest(x,xp,nearest,distsq);
if (!open_faces[2]){
point_on_line_segment(corner1,corner2,x,xp);
distsq = closest(x,xp,nearest,distsq);
crad = -2.0*(radiuslo + (nearest[1]-lo)*(radiushi-radiuslo)/(hi-lo));
}
if (!open_faces[0]) {
point_on_line_segment(corner1,corner3,x,xp);
distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq);
if (distsq < distsqprev) crad = 0;
}
if (!open_faces[1]) {
point_on_line_segment(corner2,corner4,x,xp);
distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq);
if (distsq < distsqprev) crad = 0;
}
add_contact(0,x,nearest[0],nearest[1],nearest[2]);
contact[0].radius = crad;
contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1;
return 0;
@ -476,8 +542,8 @@ int RegCone::surface_exterior(double *x, double cutoff)
// z is far enough from cone that there is no contact
// z is interior to cone
if (r >= maxradius+cutoff ||
x[2] <= lo-cutoff || x[2] >= hi+cutoff) return 0;
if (r >= maxradius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff)
return 0;
if (r < currentradius && x[2] > lo && x[2] < hi) return 0;
// z is exterior to cone or on its surface
@ -501,52 +567,35 @@ int RegCone::surface_exterior(double *x, double cutoff)
corner4[2] = hi;
distsq = BIG;
point_on_line_segment(corner1,corner2,x,xp);
distsq = closest(x,xp,nearest,distsq);
point_on_line_segment(corner1,corner3,x,xp);
distsq = closest(x,xp,nearest,distsq);
point_on_line_segment(corner2,corner4,x,xp);
distsq = closest(x,xp,nearest,distsq);
if (!open_faces[2]){
point_on_line_segment(corner1,corner2,x,xp);
distsq = closest(x,xp,nearest,distsq);
crad = -2.0*(radiuslo + (nearest[2]-lo)*(radiushi-radiuslo)/(hi-lo));
}
if (!open_faces[0]) {
point_on_line_segment(corner1,corner3,x,xp);
distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq);
if (distsq < distsqprev) crad = 0;
}
if (!open_faces[1]) {
point_on_line_segment(corner2,corner4,x,xp);
distsqprev = distsq;
distsq = closest(x,xp,nearest,distsq);
if (distsq < distsqprev) crad = 0;
}
add_contact(0,x,nearest[0],nearest[1],nearest[2]);
contact[0].radius = crad;
contact[0].iwall = 0;
if (contact[0].r < cutoff) return 1;
return 0;
}
}
/* ----------------------------------------------------------------------
find nearest point to C on line segment A,B and return it as D
project (C-A) onto (B-A)
t = length of that projection, normalized by length of (B-A)
t <= 0, C is closest to A
t >= 1, C is closest to B
else closest point is between A and B
------------------------------------------------------------------------- */
void RegCone::point_on_line_segment(double *a, double *b,
double *c, double *d)
{
double ba[3],ca[3];
subtract(a,b,ba);
subtract(a,c,ca);
double t = dotproduct(ca,ba) / dotproduct(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];
}
}
/* ---------------------------------------------------------------------- */
double RegCone::closest(double *x, double *near, double *nearest, double dsq)
@ -562,23 +611,3 @@ double RegCone::closest(double *x, double *near, double *nearest, double dsq)
nearest[2] = near[2];
return rsq;
}
/* ----------------------------------------------------------------------
v3 = v2 - v1
------------------------------------------------------------------------- */
void RegCone::subtract(double *v1, double *v2, double *v3)
{
v3[0] = v2[0] - v1[0];
v3[1] = v2[1] - v1[1];
v3[2] = v2[2] - v1[2];
}
/* ----------------------------------------------------------------------
return dotproduct = v1 dot v2
------------------------------------------------------------------------- */
double RegCone::dotproduct(double *v1, double *v2)
{
return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
}