Page MenuHomec4science
Files F91050847

compute_voronoi_atom.cpp
No OneTemporary
Actions

File Metadata

Created
Thu, Nov 7, 08:29

compute_voronoi_atom.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, 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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Daniel Schwen
------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "compute_voronoi_atom.h"
#include "atom.h"
#include "group.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include "comm.h"
#include "variable.h"
#include "input.h"
#include "force.h"
#include <vector>
using namespace LAMMPS_NS;
using namespace voro;
#define FACESDELTA 10000
/* ---------------------------------------------------------------------- */
ComputeVoronoi::ComputeVoronoi(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
int sgroup;
size_peratom_cols = 2;
peratom_flag = 1;
comm_forward = 1;
faces_flag = 0;
surface = VOROSURF_NONE;
maxedge = 0;
fthresh = ethresh = 0.0;
radstr = NULL;
onlyGroup = false;
occupation = false;
con_mono = NULL;
con_poly = NULL;
tags = NULL;
occvec = sendocc = lroot = lnext = NULL;
faces = NULL;
int iarg = 3;
while ( iarg<narg ) {
if (strcmp(arg[iarg], "occupation") == 0) {
occupation = true;
iarg++;
}
else if (strcmp(arg[iarg], "only_group") == 0) {
onlyGroup = true;
iarg++;
}
else if (strcmp(arg[iarg], "radius") == 0) {
if (iarg + 2 > narg || strstr(arg[iarg+1],"v_") != arg[iarg+1] )
error->all(FLERR,"Illegal compute voronoi/atom command");
int n = strlen(&arg[iarg+1][2]) + 1;
radstr = new char[n];
strcpy(radstr,&arg[iarg+1][2]);
iarg += 2;
}
else if (strcmp(arg[iarg], "surface") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
// group all is a special case where we just skip group testing
if(strcmp(arg[iarg+1], "all") == 0) {
surface = VOROSURF_ALL;
} else {
sgroup = group->find(arg[iarg+1]);
if (sgroup == -1) error->all(FLERR,"Could not find compute/voronoi surface group ID");
sgroupbit = group->bitmask[sgroup];
surface = VOROSURF_GROUP;
}
size_peratom_cols = 3;
iarg += 2;
} else if (strcmp(arg[iarg], "edge_histo") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
maxedge = force->inumeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg], "face_threshold") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
fthresh = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg], "edge_threshold") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
ethresh = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg], "neighbors") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
if (strcmp(arg[iarg+1],"yes") == 0) faces_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) faces_flag = 0;
else error->all(FLERR,"Illegal compute voronoi/atom command");
iarg += 2;
} else if (strcmp(arg[iarg], "peratom") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
if (strcmp(arg[iarg+1],"yes") == 0) peratom_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) peratom_flag = 0;
else error->all(FLERR,"Illegal compute voronoi/atom command");
iarg += 2;
}
else error->all(FLERR,"Illegal compute voronoi/atom command");
}
if (occupation && ( surface!=VOROSURF_NONE || maxedge>0 ) )
error->all(FLERR,"Illegal compute voronoi/atom command (occupation and (surface or edges))");
nmax = rmax = 0;
edge = rfield = sendvector = NULL;
voro = NULL;
if ( maxedge > 0 ) {
vector_flag = 1;
size_vector = maxedge+1;
memory->create(edge,maxedge+1,"voronoi/atom:edge");
memory->create(sendvector,maxedge+1,"voronoi/atom:sendvector");
vector = edge;
}
// store local face data: i, j, area
if (faces_flag) {
local_flag = 1;
size_local_cols = 3;
size_local_rows = 0;
nfacesmax = 0;
}
}
/* ---------------------------------------------------------------------- */
ComputeVoronoi::~ComputeVoronoi()
{
memory->destroy(edge);
memory->destroy(rfield);
memory->destroy(sendvector);
memory->destroy(voro);
delete[] radstr;
// voro++ container classes
delete con_mono;
delete con_poly;
// occupation analysis stuff
memory->destroy(lroot);
memory->destroy(lnext);
memory->destroy(occvec);
#ifdef NOTINPLACE
memory->destroy(sendocc);
#endif
memory->destroy(tags);
memory->destroy(faces);
}
/* ---------------------------------------------------------------------- */
void ComputeVoronoi::init()
{
}
/* ----------------------------------------------------------------------
gather compute vector data from other nodes
------------------------------------------------------------------------- */
void ComputeVoronoi::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow per atom array if necessary
int nlocal = atom->nlocal;
if (nlocal > nmax) {
memory->destroy(voro);
nmax = atom->nmax;
memory->create(voro,nmax,size_peratom_cols,"voronoi/atom:voro");
array_atom = voro;
}
// decide between occupation or per-frame tesselation modes
if (occupation) {
// build cells only once
int i, nall = nlocal + atom->nghost;
if (con_mono==NULL && con_poly==NULL) {
// generate the voronoi cell network for the initial structure
buildCells();
// save tags of atoms (i.e. of each voronoi cell)
memory->create(tags,nall,"voronoi/atom:tags");
for (i=0; i<nall; i++) tags[i] = atom->tag[i];
// linked list structure for cell occupation count on the atoms
oldnall= nall;
memory->create(lroot,nall,"voronoi/atom:lroot"); // point to first atom index in cell (or -1 for empty cell)
lnext = NULL;
lmax = 0;
// build the occupation buffer
oldnatoms = atom->natoms;
memory->create(occvec,oldnatoms,"voronoi/atom:occvec");
#ifdef NOTINPLACE
memory->create(sendocc,oldnatoms,"voronoi/atom:sendocc");
#endif
}
// get the occupation of each original voronoi cell
checkOccupation();
} else {
// build cells for each output
buildCells();
loopCells();
}
}
void ComputeVoronoi::buildCells()
{
int i;
const double e = 0.01;
int nlocal = atom->nlocal;
int dim = domain->dimension;
// in the onlyGroup mode we are not setting values for all atoms later in the voro loop
// initialize everything to zero here
if (onlyGroup) {
if (surface == VOROSURF_NONE)
for (i = 0; i < nlocal; i++) voro[i][0] = voro[i][1] = 0.0;
else
for (i = 0; i < nlocal; i++) voro[i][0] = voro[i][1] = voro[i][2] = 0.0;
}
double *sublo = domain->sublo, *sublo_lamda = domain->sublo_lamda, *boxlo = domain->boxlo;
double *subhi = domain->subhi, *subhi_lamda = domain->subhi_lamda;
double *cut = comm->cutghost;
double sublo_bound[3], subhi_bound[3], cut_bound[3];
double **x = atom->x;
// setup bounds for voro++ domain for orthogonal and triclinic simulation boxes
if( domain->triclinic ) {
// triclinic box: embed parallelepiped into orthogonal voro++ domain
// cutghost is in lamda coordinates for triclinic boxes, use subxx_lamda
double *h = domain->h;
for( i=0; i<3; ++i ) {
sublo_bound[i] = sublo[i]-cut[i]-e;
subhi_bound[i] = subhi[i]+cut[i]+e;
if (domain->periodicity[i]==0) {
sublo_bound[i] = MAX(sublo_bound[i],0.0);
subhi_bound[i] = MIN(subhi_bound[i],1.0);
}
}
if (dim == 2) {
sublo_bound[2] = 0.0;
subhi_bound[2] = 1.0;
}
sublo_bound[0] = h[0]*sublo_bound[0] + h[5]*sublo_bound[1] + h[4]*sublo_bound[2] + boxlo[0];
sublo_bound[1] = h[1]*sublo_bound[1] + h[3]*sublo_bound[2] + boxlo[1];
sublo_bound[2] = h[2]*sublo_bound[2] + boxlo[2];
subhi_bound[0] = h[0]*subhi_bound[0] + h[5]*subhi_bound[1] + h[4]*subhi_bound[2] + boxlo[0];
subhi_bound[1] = h[1]*subhi_bound[1] + h[3]*subhi_bound[2] + boxlo[1];
subhi_bound[2] = h[2]*subhi_bound[2] + boxlo[2];
} else {
// orthogonal box
for( i=0; i<3; ++i ) {
sublo_bound[i] = sublo[i]-cut[i]-e;
subhi_bound[i] = subhi[i]+cut[i]+e;
if (domain->periodicity[i]==0) {
sublo_bound[i] = MAX(sublo_bound[i],domain->boxlo[i]);
subhi_bound[i] = MIN(subhi_bound[i],domain->boxhi[i]);
}
}
if (dim == 2) {
sublo_bound[2] = sublo[2];
subhi_bound[2] = subhi[2];
}
}
// n = # of voro++ spatial hash cells (with approximately cubic cells)
int nall = nlocal + atom->nghost;
double n[3], V;
for( i=0; i<3; ++i ) n[i] = subhi_bound[i] - sublo_bound[i];
V = n[0]*n[1]*n[2];
for( i=0; i<3; ++i ) {
n[i] = round( n[i]*pow( double(nall)/(V*8.0), 0.333333 ) );
n[i] = n[i]==0 ? 1 : n[i];
}
// clear edge statistics
if ( maxedge > 0 )
for (i = 0; i <= maxedge; ++i) edge[i]=0;
// initialize voro++ container
// preallocates 8 atoms per cell
// voro++ allocates more memory if needed
int *mask = atom->mask;
if (radstr) {
// check and fetch atom style variable data
int radvar = input->variable->find(radstr);
if (radvar < 0)
error->all(FLERR,"Variable name for voronoi radius does not exist");
if (!input->variable->atomstyle(radvar))
error->all(FLERR,"Variable for voronoi radius is not atom style");
// prepare destination buffer for variable evaluation
if (nlocal > rmax) {
memory->destroy(rfield);
rmax = atom->nmax;
memory->create(rfield,rmax,"voronoi/atom:rfield");
}
// compute atom style radius variable
input->variable->compute_atom(radvar,0,rfield,1,0);
// communicate values to ghost atoms of neighboring nodes
comm->forward_comm_compute(this);
// polydisperse voro++ container
delete con_poly;
con_poly = new container_poly(sublo_bound[0],
subhi_bound[0],
sublo_bound[1],
subhi_bound[1],
sublo_bound[2],
subhi_bound[2],
int(n[0]),int(n[1]),int(n[2]),
false,false,false,8);
// pass coordinates for local and ghost atoms to voro++
for (i = 0; i < nall; i++) {
if( !onlyGroup || (mask[i] & groupbit) )
con_poly->put(i,x[i][0],x[i][1],x[i][2],rfield[i]);
}
} else {
// monodisperse voro++ container
delete con_mono;
con_mono = new container(sublo_bound[0],
subhi_bound[0],
sublo_bound[1],
subhi_bound[1],
sublo_bound[2],
subhi_bound[2],
int(n[0]),int(n[1]),int(n[2]),
false,false,false,8);
// pass coordinates for local and ghost atoms to voro++
for (i = 0; i < nall; i++)
if( !onlyGroup || (mask[i] & groupbit) )
con_mono->put(i,x[i][0],x[i][1],x[i][2]);
}
}
void ComputeVoronoi::checkOccupation()
{
// clear occupation vector
memset(occvec, 0, oldnatoms*sizeof(*occvec));
int i, j, k,
nlocal = atom->nlocal,
nall = atom->nghost + nlocal;
double rx, ry, rz,
**x = atom->x;
// prepare destination buffer for variable evaluation
if (nall > lmax) {
memory->destroy(lnext);
lmax = atom->nmax;
memory->create(lnext,lmax,"voronoi/atom:lnext");
}
// clear lroot
for (i=0; i<oldnall; ++i) lroot[i] = -1;
// clear lnext
for (i=0; i<nall; ++i) lnext[i] = -1;
// loop over all local atoms and find out in which of the local first frame voronoi cells the are in
// (need to loop over ghosts, too, to get correct occupation numbers for the second column)
for (i=0; i<nall; ++i) {
// again: find_voronoi_cell() should be in the common base class. Why it is not, I don't know. Ask the voro++ author.
if (( radstr && con_poly->find_voronoi_cell(x[i][0], x[i][1], x[i][2], rx, ry, rz, k)) ||
( !radstr && con_mono->find_voronoi_cell(x[i][0], x[i][1], x[i][2], rx, ry, rz, k) )) {
// increase occupation count of this particular cell
// only for local atoms, as we do an MPI reduce sum later
if (i<nlocal) occvec[tags[k]-1]++;
// add this atom to the linked list of cell j
if (lroot[k]<0)
lroot[k]=i;
else {
j = lroot[k];
while (lnext[j]>=0) j=lnext[j];
lnext[j] = i;
}
}
}
// MPI sum occupation
#ifdef NOTINPLACE
memcpy(sendocc, occvec, oldnatoms*sizeof(*occvec));
MPI_Allreduce(sendocc, occvec, oldnatoms, MPI_INT, MPI_SUM, world);
#else
MPI_Allreduce(MPI_IN_PLACE, occvec, oldnatoms, MPI_INT, MPI_SUM, world);
#endif
// determine the total number of atoms in this atom's currently occupied cell
int c;
for (i=0; i<oldnall; i++) { // loop over lroot (old voronoi cells)
// count
c = 0;
j = lroot[i];
while (j>=0) {
c++;
j = lnext[j];
}
// set
j = lroot[i];
while (j>=0) {
voro[j][1] = c;
j = lnext[j];
}
}
// cherry pick currently owned atoms
for (i=0; i<nlocal; i++) {
// set the new atom count in the atom's first frame voronoi cell
voro[i][0] = occvec[atom->tag[i]-1];
}
}
void ComputeVoronoi::loopCells()
{
// invoke voro++ and fetch results for owned atoms in group
voronoicell_neighbor c;
int i;
if (faces_flag) nfaces = 0;
if (radstr) {
c_loop_all cl(*con_poly);
if (cl.start()) do if (con_poly->compute_cell(c,cl)) {
i = cl.pid();
processCell(c,i);
} while (cl.inc());
} else {
c_loop_all cl(*con_mono);
if (cl.start()) do if (con_mono->compute_cell(c,cl)) {
i = cl.pid();
processCell(c,i);
} while (cl.inc());
}
if (faces_flag) size_local_rows = nfaces;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
void ComputeVoronoi::processCell(voronoicell_neighbor &c, int i)
{
int j,k, *mask = atom->mask;
std::vector<int> neigh, norder, vlist;
std::vector<double> narea, vcell;
bool have_narea = false;
// zero out surface area if surface computation was requested
if (surface != VOROSURF_NONE && !onlyGroup) voro[i][2] = 0.0;
if (i < atom->nlocal && (mask[i] & groupbit)) {
// cell volume
voro[i][0] = c.volume();
// number of cell faces
c.neighbors(neigh);
int neighs = neigh.size();
if (fthresh > 0) {
// count only faces above area threshold
c.face_areas(narea);
have_narea = true;
voro[i][1] = 0.0;
for (j=0; j<narea.size(); ++j)
if (narea[j] > fthresh) voro[i][1] += 1.0;
} else {
// unthresholded face count
voro[i][1] = neighs;
}
// cell surface area
if (surface == VOROSURF_ALL) {
voro[i][2] = c.surface_area();
} else if (surface == VOROSURF_GROUP) {
if (!have_narea) c.face_areas(narea);
voro[i][2] = 0.0;
// each entry in neigh should correspond to an entry in narea
if (neighs != narea.size())
error->one(FLERR,"Voro++ error: narea and neigh have a different size");
// loop over all faces (neighbors) and check if they are in the surface group
for (j=0; j<neighs; ++j)
if (neigh[j] >= 0 && mask[neigh[j]] & sgroupbit)
voro[i][2] += narea[j];
}
// histogram of number of face edges
if (maxedge>0) {
if (ethresh > 0) {
// count only edges above length threshold
c.vertices(vcell);
c.face_vertices(vlist); // for each face: vertex count followed list of vertex indices (n_1,v1_1,v2_1,v3_1,..,vn_1,n_2,v2_1,...)
double dx, dy, dz, r2, t2 = ethresh*ethresh;
for( j=0; j<vlist.size(); j+=vlist[j]+1 ) {
int a, b, nedge = 0;
// vlist[j] contains number of vertex indices for the current face
for( k=0; k<vlist[j]; ++k ) {
a = vlist[j+1+k]; // first vertex in edge
b = vlist[j+1+(k+1)%vlist[j]]; // second vertex in edge (possible wrap around to first vertex in list)
dx = vcell[a*3] - vcell[b*3];
dy = vcell[a*3+1] - vcell[b*3+1];
dz = vcell[a*3+2] - vcell[b*3+2];
r2 = dx*dx+dy*dy+dz*dz;
if (r2 > t2) nedge++;
}
// counted edges above threshold, now put into the correct bin
if (nedge>0) {
if (nedge<=maxedge)
edge[nedge-1]++;
else
edge[maxedge]++;
}
}
} else {
// unthresholded edge counts
c.face_orders(norder);
for (j=0; j<voro[i][1]; ++j)
if (norder[j]>0) {
if (norder[j]<=maxedge)
edge[norder[j]-1]++;
else
edge[maxedge]++;
}
}
}
// store info for local faces
if (faces_flag) {
if (nfaces+voro[i][1] > nfacesmax) {
while (nfacesmax < nfaces+voro[i][1]) nfacesmax += FACESDELTA;
memory->grow(faces,nfacesmax,size_local_cols,"compute/voronoi/atom:faces");
array_local = faces;
}
if (!have_narea) c.face_areas(narea);
if (neighs != narea.size())
error->one(FLERR,"Voro++ error: narea and neigh have a different size");
tagint itag, jtag;
tagint *tag = atom->tag;
itag = tag[i];
for (j=0; j<neighs; ++j)
if (narea[j] > fthresh) {
// external faces assigned the tag 0
int jj = neigh[j];
if (jj >= 0) jtag = tag[jj];
else jtag = 0;
faces[nfaces][0] = itag;
faces[nfaces][1] = jtag;
faces[nfaces][2] = narea[j];
nfaces++;
}
}
} else if (i < atom->nlocal) voro[i][0] = voro[i][1] = 0.0;
}
double ComputeVoronoi::memory_usage()
{
double bytes = size_peratom_cols * nmax * sizeof(double);
// estimate based on average coordination of 12
if (faces_flag) bytes += 12 * size_local_cols * nmax * sizeof(double);
return bytes;
}
void ComputeVoronoi::compute_vector()
{
invoked_vector = update->ntimestep;
if( invoked_peratom < invoked_vector ) compute_peratom();
for( int i=0; i<size_vector; ++i ) sendvector[i] = edge[i];
MPI_Allreduce(sendvector,edge,size_vector,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
void ComputeVoronoi::compute_local()
{
invoked_local = update->ntimestep;
if( invoked_peratom < invoked_local ) compute_peratom();
}
/* ---------------------------------------------------------------------- */
int ComputeVoronoi::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,m=0;
for (i = 0; i < n; ++i) buf[m++] = rfield[list[i]];
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeVoronoi::unpack_forward_comm(int n, int first, double *buf)
{
int i,last,m=0;
last = first + n;
for (i = first; i < last; ++i) rfield[i] = buf[m++];
}

Event Timeline

c4science · Help