180 std::vector<T> local_info(info_length);
186 for (
auto i = 0; i < 2; ++i)
188 local_info.at(d) =
vertices.at(i)[d];
192 MPI_Allgather(local_info.data(), info_length, MPIDataTypeT,
193 generator_points.data(), info_length, MPIDataTypeT,
198 voro::container con_old(
199 this->
sysSize.at(0), this->sysSize.at(1), this->sysSize.at(2),
200 this->sysSize.at(3), this->sysSize.at(4), this->sysSize.at(5),
202 false,
false,
false, nDomains + 10);
205 for (
auto d = 0; d < nDomains; ++d) {
206 con_old.put(d, generator_points.at(info_length * d),
207 generator_points.at(info_length * d + 1),
208 generator_points.at(info_length * d + 2));
211#ifdef ALL_DEBUG_ENABLED
213 std::ostringstream ss_local_gp_pov;
214 ss_local_gp_pov <<
"voronoi/generator_points_" << std::setw(7)
215 << std::setfill(
'0') << loadbalancing_step <<
".pov";
216 std::ostringstream ss_local_gp_gnu;
217 ss_local_gp_gnu <<
"voronoi/generator_points_" << std::setw(7)
218 << std::setfill(
'0') << loadbalancing_step <<
".gnu";
219 std::ostringstream ss_local_vc_pov;
220 ss_local_vc_pov <<
"voronoi/voronoi_cells_" << std::setw(7)
221 << std::setfill(
'0') << loadbalancing_step <<
".pov";
222 std::ostringstream ss_local_vc_gnu;
223 ss_local_vc_gnu <<
"voronoi/voronoi_cells_" << std::setw(7)
224 << std::setfill(
'0') << loadbalancing_step <<
".gnu";
226 loadbalancing_step++;
229 con_old.draw_particles(ss_local_gp_gnu.str().c_str());
230 con_old.draw_cells_gnuplot(ss_local_vc_gnu.str().c_str());
231 con_old.draw_particles_pov(ss_local_gp_pov.str().c_str());
232 con_old.draw_cells_pov(ss_local_vc_pov.str().c_str());
237 std::vector<double> cell_vertices;
238 voro::voronoicell_neighbor c;
241 voro::c_loop_all cl_old(con_old);
248 std::map<int, double> neig_area;
250 for (
auto i = 0; i < 1; ++i) {
252 std::vector<std::vector<int>> next_neighbors(
253 std::max((
int)neighbors.size(), 1));
254 std::vector<std::vector<double>> neighbors_area(
255 std::max((
int)neighbors.size(), 1));
257 if (cl_old.start()) {
260 cl_old.pos(pid, x, y, z, r);
263 con_old.compute_cell(c, cl_old);
264 c.neighbors(next_neighbors.at(idx));
265 c.face_areas(neighbors_area.at(idx));
266 for (
int j = 0; j < (int)next_neighbors.at(idx).size(); ++j)
267 neig_area.insert(std::pair<int, double>(
268 next_neighbors.at(idx).at(j), neighbors_area.at(idx).at(j)));
270 if (idx == (
int)neighbors.size())
274 if (std::count(neighbors.begin(), neighbors.end(), pid) == 1) {
275 con_old.compute_cell(c, cl_old);
276 c.neighbors(next_neighbors.at(idx));
278 if (idx == (
int)neighbors.size())
282 }
while (cl_old.inc());
284 for (
auto it = next_neighbors.begin(); it != next_neighbors.end(); ++it) {
285 neighbors.insert(neighbors.begin(), it->begin(), it->end());
288 std::sort(neighbors.begin(), neighbors.end());
289 auto uniq_it = std::unique(neighbors.begin(), neighbors.end());
290 neighbors.resize(std::distance(neighbors.begin(), uniq_it));
292 std::remove(neighbors.begin(), neighbors.end(), this->localRank),
294 neighbors.erase(std::remove_if(neighbors.begin(), neighbors.end(),
295 [](
int x) { return x < 0; }),
299 std::vector<double> volumes(neighbors.size());
300 std::vector<double> surfaces(neighbors.size());
302 for (
int i = 0; i < (int)neighbors.size(); ++i) {
303 auto it = neig_area.find(neighbors.at(i));
304 surfaces.at(i) = it->second;
310 if (cl_old.start()) {
312 cl_old.pos(pid, x, y, z, r);
313 auto it = std::find(neighbors.begin(), neighbors.end(), pid);
315 if (it != neighbors.end()) {
316 int idx = std::distance(neighbors.begin(), it);
317 con_old.compute_cell(c, cl_old);
318 volumes.at(idx) = c.volume();
321 con_old.compute_cell(c, cl_old);
322 local_volume = c.volume();
324 }
while (cl_old.inc());
326#ifdef ALL_DEBUG_ENABLED
329 std::cout <<
"ALL::Voronoi_LB<T,W>::balance begin checking..." << std::endl;
336 T work_load = (T)local_info.at(info_length - 1);
337 for (
auto it = neighbors.begin(); it != neighbors.end(); ++it) {
339 work_load += generator_points.at(info_length * neighbor + info_length - 1);
342#ifdef ALL_DEBUG_ENABLED
345 std::cout <<
"ALL::Voronoi_LB<T,W>::balance find neighbors..." << std::endl;
350 for (
auto it = neighbors.begin(); it != neighbors.end(); ++it) {
354 for (
int d = 0; d <
dimension && correct; ++d) {
355 diff.at(d) = generator_points.at(info_length * neighbor +
dimension + d) -
359 0.5 * (this->sysSize[2 * d + 1] - this->sysSize[2 * d])) {
362 }
else if (diff.at(d) <
363 -0.5 * (this->sysSize[2 * d + 1] - this->sysSize[2 * d])) {
368 max_diff = std::min(max_diff,
369 sqrt(diff.at(0) * diff.at(0) + diff.at(1) * diff.at(1) +
370 diff.at(2) * diff.at(2)));
373 T work_diff = (T)0.0;
375 (generator_points.at(info_length * neighbor + info_length - 1) -
376 local_info.at(info_length - 1));
378 (generator_points.at(info_length * neighbor + info_length - 1) +
379 local_info.at(info_length - 1));
381 if (work_density < 1e-6)
384 work_diff /= work_density;
386 shift.at(d) += 0.5 * work_diff * diff.at(d) / this->
gamma;
391#ifdef ALL_DEBUG_ENABLED
394 std::cout <<
"ALL::Voronoi_LB<T,W>::balance after neighbors found..."
398 norm = sqrt(shift.at(0) * shift.at(0) + shift.at(1) * shift.at(1) +
399 shift.at(2) * shift.at(2));
403 if (norm > 0.45 * max_diff)
404 scale = 0.45 * max_diff / norm;
406#ifdef ALL_DEBUG_ENABLED
409 std::cout <<
"ALL::Voronoi_LB<T,W>::balance find new neighbors..." << std::endl;
414 local_info.at(d) += scale * shift.at(d);
417 local_info.at(d) = (local_info.at(d) < this->
sysSize.at(2 * d))
418 ? this->
sysSize.at(2 * d) + 1.0
419 : ((local_info.at(d) >= this->
sysSize.at(2 * d + 1))
420 ? this->sysSize.at(2 * d + 1) - 1.0
424 MPI_Allgather(local_info.data(), info_length, MPIDataTypeT,
425 generator_points.data(), info_length, MPIDataTypeT,
428#ifdef ALL_DEBUG_ENABLED
431 std::cout <<
"ALL::Voronoi_LB<T,W>::balance create new containers..."
437 voro::container con_new(
438 this->
sysSize.at(0), this->sysSize.at(1), this->sysSize.at(2),
439 this->sysSize.at(3), this->sysSize.at(4), this->sysSize.at(5),
441 false,
false,
false, nDomains + 10);
444 for (
int d = 0; d < nDomains; ++d) {
445 con_new.put(d, generator_points.at(info_length * d),
446 generator_points.at(info_length * d + 1),
447 generator_points.at(info_length * d + 2));
450#ifdef ALL_DEBUG_ENABLED
451 std::ostringstream ss_local_gp_pov2;
452 ss_local_gp_pov2 <<
"voronoi/generator_points_" << std::setw(7)
453 << std::setfill(
'0') << loadbalancing_step <<
".pov";
454 std::ostringstream ss_local_gp_gnu2;
455 ss_local_gp_gnu2 <<
"voronoi/generator_points_" << std::setw(7)
456 << std::setfill(
'0') << loadbalancing_step <<
".gnu";
457 std::ostringstream ss_local_vc_pov2;
458 ss_local_vc_pov2 <<
"voronoi/voronoi_cells_" << std::setw(7)
459 << std::setfill(
'0') << loadbalancing_step <<
".pov";
460 std::ostringstream ss_local_vc_gnu2;
461 ss_local_vc_gnu2 <<
"voronoi/voronoi_cells_" << std::setw(7)
462 << std::setfill(
'0') << loadbalancing_step <<
".gnu";
464 loadbalancing_step++;
467 con_new.draw_particles(ss_local_gp_gnu2.str().c_str());
468 con_new.draw_cells_gnuplot(ss_local_vc_gnu2.str().c_str());
469 con_new.draw_particles_pov(ss_local_gp_pov2.str().c_str());
470 con_new.draw_cells_pov(ss_local_vc_pov2.str().c_str());
474#ifdef ALL_DEBUG_ENABLED
477 std::cout <<
"ALL::Voronoi_LB<T,W>::balance storing shifted vertices..."
486 voro::c_loop_all cl_new(con_new);
491#ifdef ALL_DEBUG_ENABLED
494 std::cout <<
"ALL::Voronoi_LB<T,W>::balance searching neighboring vertices..."
500 std::vector<std::vector<int>> next_neighbors(
501 std::max((
int)neighbors.size(), 1));
503 if (cl_new.start()) {
506 cl_new.pos(pid, x, y, z, r);
509 con_new.compute_cell(c, cl_new);
510 c.neighbors(next_neighbors.at(idx));
512 if (idx == (
int)neighbors.size())
516 if (std::count(neighbors.begin(), neighbors.end(), pid) == 1) {
517 con_new.compute_cell(c, cl_new);
518 c.neighbors(next_neighbors.at(idx));
520 if (idx == (
int)neighbors.size())
524 }
while (cl_new.inc());
526 for (
auto it = next_neighbors.begin(); it != next_neighbors.end(); ++it) {
527 neighbors.insert(neighbors.begin(), it->begin(), it->end());
530 std::sort(neighbors.begin(), neighbors.end());
531 auto uniq_it = std::unique(neighbors.begin(), neighbors.end());
532 neighbors.resize(std::distance(neighbors.begin(), uniq_it));
534 std::remove(neighbors.begin(), neighbors.end(), this->localRank),
536 neighbors.erase(std::remove_if(neighbors.begin(), neighbors.end(),
537 [](
int x) { return x < 0; }),
542 nNeighbors[0] = neighbors.size();
549 for (
auto n : neighbors) {