1+ #include < limits>
2+ #include < cmath>
3+ #include < algorithm>
4+ #include " bin_manager.h"
5+
6+
7+
8+ void BinManager::init_bins (
9+ double sr,
10+ const std::vector<NeighborAtom>& inside_atoms,
11+ const std::vector<NeighborAtom>& ghost_atoms
12+ )
13+ {
14+ sradius = sr;
15+ if (inside_atoms.empty () && ghost_atoms.empty ())
16+ {
17+ x_min=y_min=z_min=0 ;
18+ x_max=y_max=z_max=0 ;
19+ nbinx=nbiny=nbinz=1 ;
20+ bins.clear ();
21+ bins.resize (1 );
22+ return ;
23+ }
24+
25+ x_min = y_min = z_min = std::numeric_limits<double >::max ();
26+
27+ x_max = y_max = z_max = std::numeric_limits<double >::lowest ();
28+
29+
30+ auto update_bounds = [&](const std::vector<NeighborAtom>& atoms)
31+ {
32+ for (const auto & atom : atoms)
33+ {
34+ x_min = std::min (x_min, atom.position_x );
35+ x_max = std::max (x_max, atom.position_x );
36+
37+ y_min = std::min (y_min, atom.position_y );
38+ y_max = std::max (y_max, atom.position_y );
39+
40+ z_min = std::min (z_min, atom.position_z );
41+ z_max = std::max (z_max, atom.position_z );
42+ }
43+ };
44+
45+ update_bounds (inside_atoms);
46+ update_bounds (ghost_atoms);
47+
48+ bin_sizex = bin_sizey = bin_sizez = sradius;
49+
50+ nbinx = std::ceil ((x_max - x_min) / bin_sizex);
51+ nbiny = std::ceil ((y_max - y_min) / bin_sizey);
52+ nbinz = std::ceil ((z_max - z_min) / bin_sizez);
53+
54+ nbinx = std::max (1 , nbinx);
55+ nbiny = std::max (1 , nbiny);
56+ nbinz = std::max (1 , nbinz);
57+
58+ int nbins = nbinx * nbiny * nbinz;
59+
60+ bins.clear ();
61+
62+ bins.resize (nbins);
63+
64+ for (int ix = 0 ; ix < nbinx; ++ix)
65+ {
66+ for (int iy = 0 ; iy < nbiny; ++iy)
67+ {
68+ for (int iz = 0 ; iz < nbinz; ++iz)
69+ {
70+ int idx = ix * nbiny * nbinz + iy * nbinz + iz;
71+
72+ bins[idx].id_x = ix;
73+ bins[idx].id_y = iy;
74+ bins[idx].id_z = iz;
75+
76+ bins[idx].atoms .clear ();
77+ }
78+ }
79+ }
80+ }
81+
82+ void BinManager::do_binning (
83+ const std::vector<NeighborAtom>& inside_atoms,
84+ const std::vector<NeighborAtom>& ghost_atoms
85+ )
86+ {
87+ auto bin_atom = [&](const NeighborAtom& atom)
88+ {
89+ int ix = std::min (
90+ std::max (int ((atom.position_x - x_min) / bin_sizex), 0 ),
91+ nbinx - 1
92+ );
93+
94+ int iy = std::min (
95+ std::max (int ((atom.position_y - y_min) / bin_sizey), 0 ),
96+ nbiny - 1
97+ );
98+
99+ int iz = std::min (
100+ std::max (int ((atom.position_z - z_min) / bin_sizez), 0 ),
101+ nbinz - 1
102+ );
103+
104+ int idx = ix * nbiny * nbinz + iy * nbinz + iz;
105+
106+ bins[idx].atoms .push_back (atom);
107+ };
108+
109+ for (const auto & atom : inside_atoms) bin_atom (atom);
110+
111+ for (const auto & atom : ghost_atoms) bin_atom (atom);
112+ }
113+
114+ void BinManager::build_atom_neighbors (
115+ NeighborList& neighbor_list,
116+ std::vector<NeighborAtom>& atoms
117+ )
118+ {
119+ assert (atoms.size () == neighbor_list.numneigh .size ());
120+
121+ double sradius2 = sradius * sradius;
122+
123+ neighbor_list.reset ();
124+
125+ for (int i = 0 ; i < atoms.size (); i++)
126+ {
127+ std::vector<int > neigh_tmp;
128+
129+ int ix = std::min (
130+ std::max (int ((atoms[i].position_x - x_min) / bin_sizex), 0 ),
131+ nbinx - 1
132+ );
133+
134+ int iy = std::min (
135+ std::max (int ((atoms[i].position_y - y_min) / bin_sizey), 0 ),
136+ nbiny - 1
137+ );
138+
139+ int iz = std::min (
140+ std::max (int ((atoms[i].position_z - z_min) / bin_sizez), 0 ),
141+ nbinz - 1
142+ );
143+
144+ for (int dx = -1 ; dx <= 1 ; dx++)
145+ {
146+ for (int dy = -1 ; dy <= 1 ; dy++)
147+ {
148+ for (int dz = -1 ; dz <= 1 ; dz++)
149+ {
150+ int jx = ix + dx;
151+ int jy = iy + dy;
152+ int jz = iz + dz;
153+
154+ if (jx < 0 || jx >= nbinx ||
155+ jy < 0 || jy >= nbiny ||
156+ jz < 0 || jz >= nbinz)
157+ continue ;
158+
159+ int nidx = jx * nbiny * nbinz + jy * nbinz + jz;
160+
161+ for (const NeighborAtom& natom : bins[nidx].atoms )
162+ {
163+ double dx = atoms[i].position_x - natom.position_x ;
164+ double dy = atoms[i].position_y - natom.position_y ;
165+ double dz = atoms[i].position_z - natom.position_z ;
166+
167+ double dist2 = dx * dx + dy * dy + dz * dz;
168+
169+ if (dist2 <= sradius2 && dist2 != 0 )
170+ {
171+ neigh_tmp.push_back (natom.atom_id );
172+ }
173+ }
174+ }
175+ }
176+ }
177+ int n = neigh_tmp.size ();
178+
179+ // std::cout<<n<<std::endl;
180+
181+ int * ptr = neighbor_list.allocator .allocate (n);
182+
183+ for (int k = 0 ; k < n; k++)
184+ {
185+ assert (ptr != nullptr );
186+ ptr[k] = neigh_tmp[k];
187+ }
188+
189+ neighbor_list.firstneigh [i] = ptr;
190+ neighbor_list.numneigh [i] = n;
191+ }
192+ }
193+
194+ void BinManager::clear ()
195+ {
196+ for (auto & bin : bins)
197+ {
198+ bin.atoms .clear ();
199+ }
200+
201+ bins.clear ();
202+ }
0 commit comments