jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14480 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2016-01-23 00:49:16 +00:00
parent 5d99bf664e
commit 3621171480
4 changed files with 132 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
src/atom_vec_line.cpp
View File

@ -21,10 +21,12 @@
#include "modify.h" #include "modify.h"
#include "force.h" #include "force.h"
#include "fix.h" #include "fix.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define EPSILON 0.001 #define EPSILON 0.001
@ -37,7 +39,7 @@ AtomVecLine::AtomVecLine(LAMMPS *lmp) : AtomVec(lmp)
comm_x_only = comm_f_only = 0; comm_x_only = comm_f_only = 0;
size_forward = 4; size_forward = 4;
size_reverse = 6; size_reverse = 6;
size_border = 10; size_border = 12;
size_velocity = 6; size_velocity = 6;
size_data_atom = 8; size_data_atom = 8;
size_data_vel = 7; size_data_vel = 7;
@ -46,7 +48,8 @@ AtomVecLine::AtomVecLine(LAMMPS *lmp) : AtomVec(lmp)
atom->line_flag = 1; atom->line_flag = 1;
atom->molecule_flag = atom->rmass_flag = 1; atom->molecule_flag = atom->rmass_flag = 1;
atom->omega_flag = atom->torque_flag = 1; atom->radius_flag = atom->omega_flag = atom->torque_flag = 1;
atom->sphere_flag = 1;
nlocal_bonus = nghost_bonus = nmax_bonus = 0; nlocal_bonus = nghost_bonus = nmax_bonus = 0;
bonus = NULL; bonus = NULL;
@ -93,6 +96,7 @@ void AtomVecLine::grow(int n)
molecule = memory->grow(atom->molecule,nmax,"atom:molecule"); molecule = memory->grow(atom->molecule,nmax,"atom:molecule");
rmass = memory->grow(atom->rmass,nmax,"atom:rmass"); rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
radius = memory->grow(atom->radius,nmax,"atom:radius");
omega = memory->grow(atom->omega,nmax,3,"atom:omega"); omega = memory->grow(atom->omega,nmax,3,"atom:omega");
torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque"); torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque");
line = memory->grow(atom->line,nmax,"atom:line"); line = memory->grow(atom->line,nmax,"atom:line");
@ -112,7 +116,7 @@ void AtomVecLine::grow_reset()
mask = atom->mask; image = atom->image; mask = atom->mask; image = atom->image;
x = atom->x; v = atom->v; f = atom->f; x = atom->x; v = atom->v; f = atom->f;
molecule = atom->molecule; rmass = atom->rmass; molecule = atom->molecule; rmass = atom->rmass;
omega = atom->omega; torque = atom->torque; radius = atom->radius; omega = atom->omega; torque = atom->torque;
line = atom->line; line = atom->line;
} }
@ -149,6 +153,7 @@ void AtomVecLine::copy(int i, int j, int delflag)
molecule[j] = molecule[i]; molecule[j] = molecule[i];
rmass[j] = rmass[i]; rmass[j] = rmass[i];
radius[j] = radius[i];
omega[j][0] = omega[i][0]; omega[j][0] = omega[i][0];
omega[j][1] = omega[i][1]; omega[j][1] = omega[i][1];
omega[j][2] = omega[i][2]; omega[j][2] = omega[i][2];
@ -212,6 +217,12 @@ void AtomVecLine::set_length(int i, double value)
nlocal_bonus--; nlocal_bonus--;
line[i] = -1; line[i] = -1;
} else bonus[line[i]].length = value; } else bonus[line[i]].length = value;
// also set radius = half of length
// unless value = 0.0, then set diameter = 1.0
radius[i] = 0.5 * value;
if (value == 0.0) radius[i] = 0.5;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -482,6 +493,8 @@ int AtomVecLine::pack_border(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (line[j] < 0) buf[m++] = ubuf(0).d; if (line[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -508,6 +521,8 @@ int AtomVecLine::pack_border(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (line[j] < 0) buf[m++] = ubuf(0).d; if (line[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -543,6 +558,8 @@ int AtomVecLine::pack_border_vel(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (line[j] < 0) buf[m++] = ubuf(0).d; if (line[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -576,6 +593,8 @@ int AtomVecLine::pack_border_vel(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (line[j] < 0) buf[m++] = ubuf(0).d; if (line[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -602,6 +621,8 @@ int AtomVecLine::pack_border_vel(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (line[j] < 0) buf[m++] = ubuf(0).d; if (line[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -641,6 +662,8 @@ int AtomVecLine::pack_border_hybrid(int n, int *list, double *buf)
for (i = 0; i < n; i++) { for (i = 0; i < n; i++) {
j = list[i]; j = list[i];
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (line[j] < 0) buf[m++] = ubuf(0).d; if (line[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -668,6 +691,8 @@ void AtomVecLine::unpack_border(int n, int first, double *buf)
type[i] = (int) ubuf(buf[m++]).i; type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i; mask[i] = (int) ubuf(buf[m++]).i;
molecule[i] = (tagint) ubuf(buf[m++]).i; molecule[i] = (tagint) ubuf(buf[m++]).i;
radius[i] = buf[m++];
rmass[i] = buf[m++];
line[i] = (int) ubuf(buf[m++]).i; line[i] = (int) ubuf(buf[m++]).i;
if (line[i] == 0) line[i] = -1; if (line[i] == 0) line[i] = -1;
else { else {
@ -704,6 +729,8 @@ void AtomVecLine::unpack_border_vel(int n, int first, double *buf)
type[i] = (int) ubuf(buf[m++]).i; type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i; mask[i] = (int) ubuf(buf[m++]).i;
molecule[i] = (tagint) ubuf(buf[m++]).i; molecule[i] = (tagint) ubuf(buf[m++]).i;
radius[i] = buf[m++];
rmass[i] = buf[m++];
line[i] = (int) ubuf(buf[m++]).i; line[i] = (int) ubuf(buf[m++]).i;
if (line[i] == 0) line[i] = -1; if (line[i] == 0) line[i] = -1;
else { else {
@ -739,6 +766,8 @@ int AtomVecLine::unpack_border_hybrid(int n, int first, double *buf)
last = first + n; last = first + n;
for (i = first; i < last; i++) { for (i = first; i < last; i++) {
molecule[i] = (tagint) ubuf(buf[m++]).i; molecule[i] = (tagint) ubuf(buf[m++]).i;
radius[i] = buf[m++];
rmass[i] = buf[m++];
line[i] = (int) ubuf(buf[m++]).i; line[i] = (int) ubuf(buf[m++]).i;
if (line[i] == 0) line[i] = -1; if (line[i] == 0) line[i] = -1;
else { else {
@ -775,6 +804,7 @@ int AtomVecLine::pack_exchange(int i, double *buf)
buf[m++] = ubuf(molecule[i]).d; buf[m++] = ubuf(molecule[i]).d;
buf[m++] = rmass[i]; buf[m++] = rmass[i];
buf[m++] = radius[i];
buf[m++] = omega[i][0]; buf[m++] = omega[i][0];
buf[m++] = omega[i][1]; buf[m++] = omega[i][1];
buf[m++] = omega[i][2]; buf[m++] = omega[i][2];
@ -816,6 +846,7 @@ int AtomVecLine::unpack_exchange(double *buf)
molecule[nlocal] = (tagint) ubuf(buf[m++]).i; molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
rmass[nlocal] = buf[m++]; rmass[nlocal] = buf[m++];
radius[nlocal] = buf[m++];
omega[nlocal][0] = buf[m++]; omega[nlocal][0] = buf[m++];
omega[nlocal][1] = buf[m++]; omega[nlocal][1] = buf[m++];
omega[nlocal][2] = buf[m++]; omega[nlocal][2] = buf[m++];
@ -851,8 +882,8 @@ int AtomVecLine::size_restart()
int n = 0; int n = 0;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (line[i] >= 0) n += 19; if (line[i] >= 0) n += 20;
else n += 17; else n += 18;
if (atom->nextra_restart) if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++) for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
@ -884,6 +915,7 @@ int AtomVecLine::pack_restart(int i, double *buf)
buf[m++] = ubuf(molecule[i]).d; buf[m++] = ubuf(molecule[i]).d;
buf[m++] = rmass[i]; buf[m++] = rmass[i];
buf[m++] = radius[i];
buf[m++] = omega[i][0]; buf[m++] = omega[i][0];
buf[m++] = omega[i][1]; buf[m++] = omega[i][1];
buf[m++] = omega[i][2]; buf[m++] = omega[i][2];
@ -931,6 +963,7 @@ int AtomVecLine::unpack_restart(double *buf)
molecule[nlocal] = (tagint) ubuf(buf[m++]).i; molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
rmass[nlocal] = buf[m++]; rmass[nlocal] = buf[m++];
radius[nlocal] = buf[m++];
omega[nlocal][0] = buf[m++]; omega[nlocal][0] = buf[m++];
omega[nlocal][1] = buf[m++]; omega[nlocal][1] = buf[m++];
omega[nlocal][2] = buf[m++]; omega[nlocal][2] = buf[m++];
@ -978,7 +1011,8 @@ void AtomVecLine::create_atom(int itype, double *coord)
v[nlocal][2] = 0.0; v[nlocal][2] = 0.0;
molecule[nlocal] = 0; molecule[nlocal] = 0;
rmass[nlocal] = 1.0; radius[nlocal] = 0.5;
rmass[nlocal] = 4.0*MY_PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal];
omega[nlocal][0] = 0.0; omega[nlocal][0] = 0.0;
omega[nlocal][1] = 0.0; omega[nlocal][1] = 0.0;
omega[nlocal][2] = 0.0; omega[nlocal][2] = 0.0;
@ -1012,6 +1046,12 @@ void AtomVecLine::data_atom(double *coord, imageint imagetmp, char **values)
if (rmass[nlocal] <= 0.0) if (rmass[nlocal] <= 0.0)
error->one(FLERR,"Invalid density in Atoms section of data file"); error->one(FLERR,"Invalid density in Atoms section of data file");
if (line[nlocal] < 0) {
radius[nlocal] = 0.5;
rmass[nlocal] *= 4.0*MY_PI/3.0 *
radius[nlocal]*radius[nlocal]*radius[nlocal];
} else radius[nlocal] = 0.0;
x[nlocal][0] = coord[0]; x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1]; x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2]; x[nlocal][2] = coord[2];
@ -1047,6 +1087,12 @@ int AtomVecLine::data_atom_hybrid(int nlocal, char **values)
if (rmass[nlocal] <= 0.0) if (rmass[nlocal] <= 0.0)
error->one(FLERR,"Invalid density in Atoms section of data file"); error->one(FLERR,"Invalid density in Atoms section of data file");
if (line[nlocal] < 0) {
radius[nlocal] = 0.5;
rmass[nlocal] *= 4.0*MY_PI/3.0 *
radius[nlocal]*radius[nlocal]*radius[nlocal];
} else radius[nlocal] = 0.0;
return 3; return 3;
} }
@ -1084,9 +1130,10 @@ void AtomVecLine::data_atom_bonus(int m, char **values)
x[m][0] = xc; x[m][0] = xc;
x[m][1] = yc; x[m][1] = yc;
// reset line mass // reset line radius and mass
// previously stored density in rmass // rmass currently holds density
radius[m] = 0.5 * length;
rmass[m] *= length; rmass[m] *= length;
bonus[nlocal_bonus].ilocal = m; bonus[nlocal_bonus].ilocal = m;
@ -1132,7 +1179,8 @@ void AtomVecLine::pack_data(double **buf)
buf[i][2] = ubuf(type[i]).d; buf[i][2] = ubuf(type[i]).d;
if (line[i] < 0) buf[i][3] = ubuf(0).d; if (line[i] < 0) buf[i][3] = ubuf(0).d;
else buf[i][3] = ubuf(1).d; else buf[i][3] = ubuf(1).d;
if (line[i] < 0) buf[i][4] = rmass[i]; if (line[i] < 0)
buf[i][4] = rmass[i] / (4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i]);
else buf[i][4] = rmass[i]/bonus[line[i]].length; else buf[i][4] = rmass[i]/bonus[line[i]].length;
buf[i][5] = x[i][0]; buf[i][5] = x[i][0];
buf[i][6] = x[i][1]; buf[i][6] = x[i][1];
@ -1152,7 +1200,8 @@ int AtomVecLine::pack_data_hybrid(int i, double *buf)
buf[0] = ubuf(molecule[i]).d; buf[0] = ubuf(molecule[i]).d;
if (line[i] < 0) buf[1] = ubuf(0).d; if (line[i] < 0) buf[1] = ubuf(0).d;
else buf[1] = ubuf(1).d; else buf[1] = ubuf(1).d;
if (line[i] < 0) buf[2] = rmass[i]; if (line[i] < 0)
buf[2] = rmass[i] / (4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i]);
else buf[2] = rmass[i]/bonus[line[i]].length; else buf[2] = rmass[i]/bonus[line[i]].length;
return 3; return 3;
} }
@ -1255,6 +1304,7 @@ bigint AtomVecLine::memory_usage()
if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax); if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax); if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
if (atom->memcheck("radius")) bytes += memory->usage(radius,nmax);
if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3); if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3);
if (atom->memcheck("torque")) if (atom->memcheck("torque"))
bytes += memory->usage(torque,nmax*comm->nthreads,3); bytes += memory->usage(torque,nmax*comm->nthreads,3);

2
src/atom_vec_line.h
View File

@ -89,7 +89,7 @@ class AtomVecLine : public AtomVec {
imageint *image; imageint *image;
double **x,**v,**f; double **x,**v,**f;
tagint *molecule; tagint *molecule;
double *rmass; double *rmass,*radius;
double **omega,**torque; double **omega,**torque;
int *line; int *line;

80
src/atom_vec_tri.cpp
View File

@ -22,10 +22,12 @@
#include "modify.h" #include "modify.h"
#include "force.h" #include "force.h"
#include "fix.h" #include "fix.h"
#include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst;
#define EPSILON 0.001 #define EPSILON 0.001
@ -38,7 +40,7 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
comm_x_only = comm_f_only = 0; comm_x_only = comm_f_only = 0;
size_forward = 7; size_forward = 7;
size_reverse = 6; size_reverse = 6;
size_border = 24; size_border = 26;
size_velocity = 6; size_velocity = 6;
size_data_atom = 8; size_data_atom = 8;
size_data_vel = 7; size_data_vel = 7;
@ -47,10 +49,13 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
atom->tri_flag = 1; atom->tri_flag = 1;
atom->molecule_flag = atom->rmass_flag = 1; atom->molecule_flag = atom->rmass_flag = 1;
atom->angmom_flag = atom->torque_flag = 1; atom->radius_flag = atom->angmom_flag = atom->torque_flag = 1;
nlocal_bonus = nghost_bonus = nmax_bonus = 0; nlocal_bonus = nghost_bonus = nmax_bonus = 0;
bonus = NULL; bonus = NULL;
if (domain->dimension != 3)
error->all(FLERR,"Atom_style tri can only be used in 3d simulations");
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -94,6 +99,7 @@ void AtomVecTri::grow(int n)
molecule = memory->grow(atom->molecule,nmax,"atom:molecule"); molecule = memory->grow(atom->molecule,nmax,"atom:molecule");
rmass = memory->grow(atom->rmass,nmax,"atom:rmass"); rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
radius = memory->grow(atom->radius,nmax,"atom:radius");
angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom"); angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom");
torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque"); torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque");
tri = memory->grow(atom->tri,nmax,"atom:tri"); tri = memory->grow(atom->tri,nmax,"atom:tri");
@ -113,7 +119,7 @@ void AtomVecTri::grow_reset()
mask = atom->mask; image = atom->image; mask = atom->mask; image = atom->image;
x = atom->x; v = atom->v; f = atom->f; x = atom->x; v = atom->v; f = atom->f;
molecule = atom->molecule; rmass = atom->rmass; molecule = atom->molecule; rmass = atom->rmass;
angmom = atom->angmom; torque = atom->torque; radius = atom->radius; angmom = atom->angmom; torque = atom->torque;
tri = atom->tri; tri = atom->tri;
} }
@ -151,6 +157,7 @@ void AtomVecTri::copy(int i, int j, int delflag)
molecule[j] = molecule[i]; molecule[j] = molecule[i];
rmass[j] = rmass[i]; rmass[j] = rmass[i];
radius[j] = radius[i];
angmom[j][0] = angmom[i][0]; angmom[j][0] = angmom[i][0];
angmom[j][1] = angmom[i][1]; angmom[j][1] = angmom[i][1];
angmom[j][2] = angmom[i][2]; angmom[j][2] = angmom[i][2];
@ -202,6 +209,9 @@ void AtomVecTri::clear_bonus()
void AtomVecTri::set_equilateral(int i, double size) void AtomVecTri::set_equilateral(int i, double size)
{ {
// also set radius = distance from center to corner-pt = len(c1)
// unless size = 0.0, then set diameter = 1.0
if (tri[i] < 0) { if (tri[i] < 0) {
if (size == 0.0) return; if (size == 0.0) return;
if (nlocal_bonus == nmax_bonus) grow_bonus(); if (nlocal_bonus == nmax_bonus) grow_bonus();
@ -226,9 +236,11 @@ void AtomVecTri::set_equilateral(int i, double size)
inertia[0] = sqrt(3.0)/96.0 * size*size*size*size; inertia[0] = sqrt(3.0)/96.0 * size*size*size*size;
inertia[1] = sqrt(3.0)/96.0 * size*size*size*size; inertia[1] = sqrt(3.0)/96.0 * size*size*size*size;
inertia[2] = sqrt(3.0)/48.0 * size*size*size*size; inertia[2] = sqrt(3.0)/48.0 * size*size*size*size;
radius[i] = MathExtra::len3(c1);
bonus[nlocal_bonus].ilocal = i; bonus[nlocal_bonus].ilocal = i;
tri[i] = nlocal_bonus++; tri[i] = nlocal_bonus++;
} else if (size == 0.0) { } else if (size == 0.0) {
radius[i] = 0.5;
copy_bonus(nlocal_bonus-1,tri[i]); copy_bonus(nlocal_bonus-1,tri[i]);
nlocal_bonus--; nlocal_bonus--;
tri[i] = -1; tri[i] = -1;
@ -249,6 +261,7 @@ void AtomVecTri::set_equilateral(int i, double size)
inertia[0] = sqrt(3.0)/96.0 * size*size*size*size; inertia[0] = sqrt(3.0)/96.0 * size*size*size*size;
inertia[1] = sqrt(3.0)/96.0 * size*size*size*size; inertia[1] = sqrt(3.0)/96.0 * size*size*size*size;
inertia[2] = sqrt(3.0)/48.0 * size*size*size*size; inertia[2] = sqrt(3.0)/48.0 * size*size*size*size;
radius[i] = MathExtra::len3(c1);
} }
} }
@ -581,6 +594,8 @@ int AtomVecTri::pack_border(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (tri[j] < 0) buf[m++] = ubuf(0).d; if (tri[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -626,6 +641,8 @@ int AtomVecTri::pack_border(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (tri[j] < 0) buf[m++] = ubuf(0).d; if (tri[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -681,6 +698,8 @@ int AtomVecTri::pack_border_vel(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (tri[j] < 0) buf[m++] = ubuf(0).d; if (tri[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -733,6 +752,8 @@ int AtomVecTri::pack_border_vel(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (tri[j] < 0) buf[m++] = ubuf(0).d; if (tri[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -778,6 +799,8 @@ int AtomVecTri::pack_border_vel(int n, int *list, double *buf,
buf[m++] = ubuf(type[j]).d; buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d; buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (tri[j] < 0) buf[m++] = ubuf(0).d; if (tri[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -837,6 +860,8 @@ int AtomVecTri::pack_border_hybrid(int n, int *list, double *buf)
for (i = 0; i < n; i++) { for (i = 0; i < n; i++) {
j = list[i]; j = list[i];
buf[m++] = ubuf(molecule[j]).d; buf[m++] = ubuf(molecule[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (tri[j] < 0) buf[m++] = ubuf(0).d; if (tri[j] < 0) buf[m++] = ubuf(0).d;
else { else {
buf[m++] = ubuf(1).d; buf[m++] = ubuf(1).d;
@ -884,6 +909,8 @@ void AtomVecTri::unpack_border(int n, int first, double *buf)
type[i] = (int) ubuf(buf[m++]).i; type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i; mask[i] = (int) ubuf(buf[m++]).i;
molecule[i] = (tagint) ubuf(buf[m++]).i; molecule[i] = (tagint) ubuf(buf[m++]).i;
radius[i] = buf[m++];
rmass[i] = buf[m++];
tri[i] = (int) ubuf(buf[m++]).i; tri[i] = (int) ubuf(buf[m++]).i;
if (tri[i] == 0) tri[i] = -1; if (tri[i] == 0) tri[i] = -1;
else { else {
@ -940,6 +967,8 @@ void AtomVecTri::unpack_border_vel(int n, int first, double *buf)
type[i] = (int) ubuf(buf[m++]).i; type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i; mask[i] = (int) ubuf(buf[m++]).i;
molecule[i] = (tagint) ubuf(buf[m++]).i; molecule[i] = (tagint) ubuf(buf[m++]).i;
radius[i] = buf[m++];
rmass[i] = buf[m++];
tri[i] = (int) ubuf(buf[m++]).i; tri[i] = (int) ubuf(buf[m++]).i;
if (tri[i] == 0) tri[i] = -1; if (tri[i] == 0) tri[i] = -1;
else { else {
@ -995,6 +1024,8 @@ int AtomVecTri::unpack_border_hybrid(int n, int first, double *buf)
last = first + n; last = first + n;
for (i = first; i < last; i++) { for (i = first; i < last; i++) {
molecule[i] = (tagint) ubuf(buf[m++]).i; molecule[i] = (tagint) ubuf(buf[m++]).i;
radius[i] = buf[m++];
rmass[i] = buf[m++];
tri[i] = (int) ubuf(buf[m++]).i; tri[i] = (int) ubuf(buf[m++]).i;
if (tri[i] == 0) tri[i] = -1; if (tri[i] == 0) tri[i] = -1;
else { else {
@ -1050,6 +1081,7 @@ int AtomVecTri::pack_exchange(int i, double *buf)
buf[m++] = ubuf(molecule[i]).d; buf[m++] = ubuf(molecule[i]).d;
buf[m++] = rmass[i]; buf[m++] = rmass[i];
buf[m++] = radius[i];
buf[m++] = angmom[i][0]; buf[m++] = angmom[i][0];
buf[m++] = angmom[i][1]; buf[m++] = angmom[i][1];
buf[m++] = angmom[i][2]; buf[m++] = angmom[i][2];
@ -1110,6 +1142,7 @@ int AtomVecTri::unpack_exchange(double *buf)
molecule[nlocal] = (tagint) ubuf(buf[m++]).i; molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
rmass[nlocal] = buf[m++]; rmass[nlocal] = buf[m++];
radius[nlocal] = buf[m++];
angmom[nlocal][0] = buf[m++]; angmom[nlocal][0] = buf[m++];
angmom[nlocal][1] = buf[m++]; angmom[nlocal][1] = buf[m++];
angmom[nlocal][2] = buf[m++]; angmom[nlocal][2] = buf[m++];
@ -1164,8 +1197,8 @@ int AtomVecTri::size_restart()
int n = 0; int n = 0;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (tri[i] >= 0) n += 33; if (tri[i] >= 0) n += 34;
else n += 17; else n += 18;
if (atom->nextra_restart) if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++) for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
@ -1197,6 +1230,7 @@ int AtomVecTri::pack_restart(int i, double *buf)
buf[m++] = ubuf(molecule[i]).d; buf[m++] = ubuf(molecule[i]).d;
buf[m++] = rmass[i]; buf[m++] = rmass[i];
buf[m++] = radius[i];
buf[m++] = angmom[i][0]; buf[m++] = angmom[i][0];
buf[m++] = angmom[i][1]; buf[m++] = angmom[i][1];
buf[m++] = angmom[i][2]; buf[m++] = angmom[i][2];
@ -1263,6 +1297,7 @@ int AtomVecTri::unpack_restart(double *buf)
molecule[nlocal] = (tagint) ubuf(buf[m++]).i; molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
rmass[nlocal] = buf[m++]; rmass[nlocal] = buf[m++];
radius[nlocal] = buf[m++];
angmom[nlocal][0] = buf[m++]; angmom[nlocal][0] = buf[m++];
angmom[nlocal][1] = buf[m++]; angmom[nlocal][1] = buf[m++];
angmom[nlocal][2] = buf[m++]; angmom[nlocal][2] = buf[m++];
@ -1329,7 +1364,8 @@ void AtomVecTri::create_atom(int itype, double *coord)
v[nlocal][2] = 0.0; v[nlocal][2] = 0.0;
molecule[nlocal] = 0; molecule[nlocal] = 0;
rmass[nlocal] = 1.0; radius[nlocal] = 0.5;
rmass[nlocal] = 4.0*MY_PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal];
angmom[nlocal][0] = 0.0; angmom[nlocal][0] = 0.0;
angmom[nlocal][1] = 0.0; angmom[nlocal][1] = 0.0;
angmom[nlocal][2] = 0.0; angmom[nlocal][2] = 0.0;
@ -1363,6 +1399,12 @@ void AtomVecTri::data_atom(double *coord, imageint imagetmp, char **values)
if (rmass[nlocal] <= 0.0) if (rmass[nlocal] <= 0.0)
error->one(FLERR,"Invalid density in Atoms section of data file"); error->one(FLERR,"Invalid density in Atoms section of data file");
if (tri[nlocal] < 0) {
radius[nlocal] = 0.5;
rmass[nlocal] *= 4.0*MY_PI/3.0 *
radius[nlocal]*radius[nlocal]*radius[nlocal];
} else radius[nlocal] = 0.0;
x[nlocal][0] = coord[0]; x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1]; x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2]; x[nlocal][2] = coord[2];
@ -1398,6 +1440,12 @@ int AtomVecTri::data_atom_hybrid(int nlocal, char **values)
if (rmass[nlocal] <= 0.0) if (rmass[nlocal] <= 0.0)
error->one(FLERR,"Invalid density in Atoms section of data file"); error->one(FLERR,"Invalid density in Atoms section of data file");
if (tri[nlocal] < 0) {
radius[nlocal] = 0.5;
rmass[nlocal] *= 4.0*MY_PI/3.0 *
radius[nlocal]*radius[nlocal]*radius[nlocal];
} else radius[nlocal] = 0.0;
return 3; return 3;
} }
@ -1457,10 +1505,19 @@ void AtomVecTri::data_atom_bonus(int m, char **values)
x[m][1] = centroid[1]; x[m][1] = centroid[1];
x[m][2] = centroid[2]; x[m][2] = centroid[2];
// reset tri mass // reset tri radius and mass
// previously stored density in rmass // rmass currently holds density
// tri area = 0.5 len(U x V), where U,V are edge vectors from one vertex // tri area = 0.5 len(U x V), where U,V are edge vectors from one vertex
double c4[3];
MathExtra::sub3(c1,centroid,c4);
radius[m] = MathExtra::lensq3(c4);
MathExtra::sub3(c2,centroid,c4);
radius[m] = MAX(radius[m],MathExtra::lensq3(c4));
MathExtra::sub3(c3,centroid,c4);
radius[m] = MAX(radius[m],MathExtra::lensq3(c4));
radius[m] = sqrt(radius[m]);
double norm[3]; double norm[3];
MathExtra::cross3(c2mc1,c3mc1,norm); MathExtra::cross3(c2mc1,c3mc1,norm);
double area = 0.5 * MathExtra::len3(norm); double area = 0.5 * MathExtra::len3(norm);
@ -1567,7 +1624,8 @@ void AtomVecTri::pack_data(double **buf)
buf[i][2] = ubuf(type[i]).d; buf[i][2] = ubuf(type[i]).d;
if (tri[i] < 0) buf[i][3] = ubuf(0).d; if (tri[i] < 0) buf[i][3] = ubuf(0).d;
else buf[i][3] = ubuf(1).d; else buf[i][3] = ubuf(1).d;
if (tri[i] < 0) buf[i][4] = rmass[i]; if (tri[i] < 0)
buf[i][4] = rmass[i] / (4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i]);
else { else {
MathExtra::sub3(bonus[tri[i]].c2,bonus[tri[i]].c1,c2mc1); MathExtra::sub3(bonus[tri[i]].c2,bonus[tri[i]].c1,c2mc1);
MathExtra::sub3(bonus[tri[i]].c3,bonus[tri[i]].c1,c3mc1); MathExtra::sub3(bonus[tri[i]].c3,bonus[tri[i]].c1,c3mc1);
@ -1593,7 +1651,8 @@ int AtomVecTri::pack_data_hybrid(int i, double *buf)
buf[0] = ubuf(molecule[i]).d; buf[0] = ubuf(molecule[i]).d;
if (tri[i] < 0) buf[1] = ubuf(0).d; if (tri[i] < 0) buf[1] = ubuf(0).d;
else buf[1] = ubuf(1).d; else buf[1] = ubuf(1).d;
if (tri[i] < 0) buf[2] = rmass[i]; if (tri[i] < 0)
buf[2] = rmass[i] / (4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i]);
else { else {
double c2mc1[3],c3mc1[3],norm[3]; double c2mc1[3],c3mc1[3],norm[3];
MathExtra::sub3(bonus[tri[i]].c2,bonus[tri[i]].c1,c2mc1); MathExtra::sub3(bonus[tri[i]].c2,bonus[tri[i]].c1,c2mc1);
@ -1703,6 +1762,7 @@ bigint AtomVecTri::memory_usage()
if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax); if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax); if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
if (atom->memcheck("radius")) bytes += memory->usage(radius,nmax);
if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3); if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3);
if (atom->memcheck("torque")) bytes += if (atom->memcheck("torque")) bytes +=
memory->usage(torque,nmax*comm->nthreads,3); memory->usage(torque,nmax*comm->nthreads,3);

2
src/atom_vec_tri.h
View File

@ -91,7 +91,7 @@ class AtomVecTri : public AtomVec {
imageint *image; imageint *image;
double **x,**v,**f; double **x,**v,**f;
tagint *molecule; tagint *molecule;
double *rmass; double *rmass,*radius;
double **angmom,**torque; double **angmom,**torque;
int *tri; int *tri;