Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
/* ----------------------------------------------------------------------
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 "math.h"
#include "string.h"
#include "stdlib.h"
#include "compute_voronoi_atom.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include "comm.h"
#include <vector>
#include "voro++.hh"
using namespace LAMMPS_NS;
using namespace voro;
/* ---------------------------------------------------------------------- */
ComputeVoronoi::ComputeVoronoi(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute voronoi/atom command");
peratom_flag = 1;
size_peratom_cols = 2;
nmax = 0;
voro = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeVoronoi::~ComputeVoronoi()
{
memory->destroy(voro);
}
/* ---------------------------------------------------------------------- */
void ComputeVoronoi::init()
{
if (domain->dimension != 3)
error->all(FLERR,"Compute voronoi/atom not allowed for 2d systems");
if (domain->triclinic != 0)
error->all(FLERR,"Compute voronoi/atom not allowed for triclinic boxes");
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"voronoi/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute voronoi/atom command");
}
/* ----------------------------------------------------------------------
gather compute vector data from other nodes
------------------------------------------------------------------------- */
void ComputeVoronoi::compute_peratom()
{
int i;
invoked_peratom = update->ntimestep;
// grow per atom array if necessary
int local = atom->nlocal;
if (nlocal > nmax) {
memory->destroy(voro);
nmax = atom->nmax;
memory->create(voro,nmax,size_peratom_cols,"voronoi/atom:voro");
array_atom = voro;
}
// n = # of voro++ spatial hash cells
// TODO: make square
int nall = nlocal + atom->nghost;
int n = int(floor( pow( double(nall)/8.0, 0.333333 ) ));
n = (n==0) ? 1 : n;
// initialize voro++ container
// preallocates 8 atoms per cell
// voro++ allocates more memory if needed
double *sublo = domain->sublo;
double *subhi = domain->subhi;
double *cut = comm->cutghost;
container con(sublo[0]-cut[0],subhi[0]+cut[0],
sublo[1]-cut[1],subhi[1]+cut[1],
sublo[2]-cut[2],subhi[2]+cut[2],
n,n,n,false,false,false,8);
// pass coordinates for local and ghost atoms to voro++
double **x = atom->x;
for (i = 0; i < nall; i++)
con.put(i,x[i][0],x[i][1],x[i][2]);
// invoke voro++ and fetch results for owned atoms in group
std::vector<int> neigh;
voronoicell_neighbor c;
c_loop_all cl(con);
if (cl.start()) do if (con.compute_cell(c,cl)) {
i = cl.pid();
if (i < nlocal && (mask[i] & groupbit)) {
voro[i][0] = c.volume();
c.neighbors(neigh);
voro[i][1] = neigh.size();
} else if (i < nlocal) voro[i][0] = voro[i][1] = 0.0;
} while (cl.inc());
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeVoronoi::memory_usage()
{
double bytes = size_peratom_cols * nmax * sizeof(double);
return bytes;
}