ALL.hpp

Go to the documentation of this file.

/*
Copyright 2018-2026 Rene Halver, Forschungszentrum Juelich GmbH, Germany
Copyright 2018-2026 Godehard Sutmann, Forschungszentrum Juelich GmbH, Germany
Copyright 2026-2026 Jonas Kroschewski, Forschungszentrum Juelich GmbH, Germany
 
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
 
1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
 
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation and/or
   other materials provided with the distribution.
 
3. Neither the name of the copyright holder nor the names of its contributors
may be used to endorse or promote products derived from this software without
   specific prior written permission.
 
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
 
#ifndef ALL_MAIN_HEADER_INC
#define ALL_MAIN_HEADER_INC
 
#include "ALL_CustomExceptions.hpp"
#include "ALL_LB.hpp"
#include "ALL_Point.hpp"
#include "ALL_Tensor.hpp"
#ifdef ALL_VORONOI_ACTIVE
#include "ALL_Voronoi.hpp"
#endif
#include "ALL_ForceBased.hpp"
#include "ALL_Histogram.hpp"
#include "ALL_Staggered.hpp"
#include <algorithm>
#include <iomanip>
#include <memory>
#include <numeric>
#include <sstream>
#include <string>
#include <vector>
 
#include <cerrno>
#include <sys/stat.h>
 
#ifdef ALL_VTK_OUTPUT
#include <limits>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkFloatArray.h>
#include <vtkInformation.h>
#include <vtkIntArray.h>
#include <vtkMPIController.h>
#include <vtkPolyhedron.h>
#include <vtkProgrammableFilter.h>
#include <vtkSmartPointer.h>
#include <vtkUnstructuredGrid.h>
#include <vtkVoxel.h>
#include <vtkXMLPUnstructuredGridWriter.h>
#include <vtkXMLUnstructuredGridWriter.h>
#include <vtkUnstructuredGridWriter.h> //For vtkXMLUnstructuredGridWriter debug
 
#ifdef VTK_CELL_ARRAY_V2
#include <vtkNew.h>
#endif
#endif
 
namespace ALL {
 
#define ALL_ESTIMATION_LIMIT_DEFAULT 0

enum LB_t : int {
  STAGGERED = 0,
  TENSOR = 1,
  FORCEBASED = 2,
#ifdef ALL_VORONOI_ACTIVE
  VORONOI = 3,
#endif
  HISTOGRAM = 4,
  TENSOR_MAX = 5,
  UNIMPLEMENTED = 128
};

 
template <class T, class W> class ALL {
public:
  ALL()
      : loadbalancing_step(0)
#ifdef ALL_VTK_OUTPUT
        ,
        vtk_init(false)
#endif
  {
    // determine correct MPI data type for template T
    if (std::is_same<T, double>::value)
      MPIDataTypeT = MPI_DOUBLE;
    else if (std::is_same<T, float>::value)
      MPIDataTypeT = MPI_FLOAT;
    else if (std::is_same<T, int>::value)
      MPIDataTypeT = MPI_INT;
    else if (std::is_same<T, long>::value)
      MPIDataTypeT = MPI_LONG;
 
    // determine correct MPI data type for template W
    if (std::is_same<W, double>::value)
      MPIDataTypeW = MPI_DOUBLE;
    else if (std::is_same<W, float>::value)
      MPIDataTypeW = MPI_FLOAT;
    else if (std::is_same<W, int>::value)
      MPIDataTypeW = MPI_INT;
    else if (std::is_same<W, long>::value)
      MPIDataTypeW = MPI_LONG;
  }

  ALL(const LB_t m, const int d, const T g) : ALL() {
    method = m;
    switch (method) {
    case LB_t::TENSOR:
      balancer.reset(new Tensor_LB<T, W>(d, (W)0, g));
      break;
    case LB_t::STAGGERED:
      balancer.reset(new Staggered_LB<T, W>(d, (W)0, g));
      break;
    case LB_t::FORCEBASED:
      balancer.reset(new ForceBased_LB<T, W>(d, (W)0, g));
      break;
#ifdef ALL_VORONOI_ACTIVE
    case LB_t::VORONOI:
      balancer.reset(new Voronoi_LB<T, W>(d, (W)0, g));
      break;
#endif
    case LB_t::HISTOGRAM:
      balancer.reset(new Histogram_LB<T, W>(d, std::vector<W>(10), g));
      break;
    case LB_t::TENSOR_MAX:
      balancer.reset(new TensorMax_LB<T, W>(d, (W)0, g));
      break;
    default:
      throw InvalidArgumentException(__FILE__, __func__, __LINE__,
                                     "Unknown type of loadbalancing passed.");
    }
    balancer->setDimension(d);
    balancer->setGamma(g);
    estimation_limit = ALL_ESTIMATION_LIMIT_DEFAULT;
  }

  ALL(const LB_t m, const int d, const std::vector<Point<T>> &inp, const T g)
      : ALL(m, d, g) {
    balancer->setVertices(inp);
    calculate_outline();
  }

  ~ALL() = default;

  void setVertices(const std::vector<Point<T>> &inp) {
    int dim = inp.at(0).getDimension();
    if (dim != balancer->getDimension())
      throw PointDimensionMissmatchException(
          __FILE__, __func__, __LINE__,
          "Dimension of ALL::Points in input vector do not match dimension of "
          "ALL "
          "object.");
    balancer->setVertices(inp);
    calculate_outline();
  }

  void setCommunicator(const MPI_Comm comm) {
    int comm_type;
    MPI_Topo_test(comm, &comm_type);
    if (comm_type == MPI_CART) {
      communicator = comm;
      balancer->setCommunicator(communicator);
      // set procGridLoc and procGridDim
      const int dimension = balancer->getDimension();
      int location[dimension];
      int size[dimension];
      int periodicity[dimension];
      MPI_Cart_get(communicator, dimension, size, periodicity, location);
      procGridLoc.assign(location, location + dimension);
      procGridDim.assign(size, size + dimension);
      procGridSet = true;
    } else {
      int rank;
      MPI_Comm_rank(comm, &rank);
      if (procGridSet) {
        // compute 1D index from the location in the process grid (using z-y-x
        // ordering, as MPI_Dims_create)
        int idx;
        switch (balancer->getDimension()) {
        case 3:
          idx = procGridLoc.at(2) + procGridLoc.at(1) * procGridDim.at(2) +
                procGridLoc.at(0) * procGridDim.at(2) * procGridDim.at(1);
          break;
        case 2:
          idx = procGridLoc.at(1) + procGridLoc.at(0) * procGridDim.at(1);
          break;
        case 1:
          idx = procGridLoc.at(0);
          break;
        default:
          throw InternalErrorException(
              __FILE__, __func__, __LINE__,
              "ALL cannot deal with more then three dimensions (or less then "
              "one).");
          break;
        }
 
        // create new temporary communicator with correct ranks
        MPI_Comm temp_comm;
        MPI_Comm_split(comm, 0, idx, &temp_comm);
 
        std::vector<int> periods(3, 1);
 
        // transform temporary communicator to cartesian communicator
        MPI_Cart_create(temp_comm, balancer->getDimension(), procGridDim.data(),
                        periods.data(), 0, &communicator);
        balancer->setCommunicator(communicator);
      } else {
        throw InternalErrorException(
            __FILE__, __func__, __LINE__,
            "When using a non-cartesian communicator process grid parameters "
            "required (not set).");
      }
    }
  }

  T getGamma() { return balancer->getGamma(); };

  void setGamma(const T g) { balancer->setGamma(g); };
 
  // method to set the automatic calculation of the correction value
  // @param[in] b the boolean to activate/deactivate automatic calculation
  void setGammaAutoCalc(const bool b) { balancer->setGammaAutoCalc(b); }

  bool getGammaAutoCalc() { return balancer->getGammaAutoCalc(); }

  void setWork(const W work) {
    // set new value for work (single value for whole domain)
    balancer->setWork(work);
  }
 
  // method to set a multi-dimensional work for the local domain
  void setWork(const std::vector<W> &work) { balancer->setWork(work); }

  void setup() { balancer->setup(); }

  std::vector<Point<T>> &balance() {
    nVertices = balancer->getNVertices();
    switch (method) {
    case LB_t::TENSOR:
      balancer->balance(loadbalancing_step);
      calculate_outline();
      break;
    case LB_t::STAGGERED:
 
      /* ToDo: Check if repetition is useful at all and
       *       change it to be included for all methods,
       *       not only STAGGERED */
      /*
      // estimation to improve speed of convergence
      // changing vertices during estimation
      std::vector<Point<T>> tmp_vertices = balancer->getVertices();
      // saved original vertices
      std::vector<Point<T>> old_vertices = balancer->getVertices();
      // old temporary vertices
      std::vector<Point<T>> old_tmp_vertices = balancer->getVertices();
      // result temporary vertices
      std::vector<Point<T>> result_vertices = balancer->getVertices();
      // temporary work
      W tmp_work = balancer->getWork().at(0);
      // old temporary work
      W old_tmp_work = balancer->getWork().at(0);
 
      W max_work;
      W sum_work;
      double d_max;
      int n_ranks;
 
      // compute work density on local domain
      double V = 1.0;
      for (int i = 0; i < balancer->getDimension(); ++i)
          V *= ( outline.at(1)[i] - outline.at(0)[i] );
      double rho = workArray.at(0) / V;
 
      // collect maximum work in system
      MPI_Allreduce(workArray.data(), &max_work, 1, MPIDataTypeW, MPI_MAX,
      communicator); MPI_Allreduce(workArray.data(), &sum_work, 1,
      MPIDataTypeW, MPI_SUM, communicator);
      MPI_Comm_size(communicator,&n_ranks);
      d_max = (double)(max_work) * (double)n_ranks / (double)sum_work - 1.0;
 
 
      // internal needs to be readded to argument list
      const bool internal = false
      for (int i = 0; i < estimation_limit && d_max > 0.1 && !internal; ++i)
      {
          balancer->setWork(tmp_work);
          balancer->setup();
          balancer->balance(true);
          bool sane = true;
          // check if the estimated boundaries are not too deformed
          for (int j = 0; j < balancer->getDimension(); ++j)
              sane = sane && (result_vertices.at(0)[j] <=
      old_vertices.at(1)[j]);
          MPI_Allreduce(MPI_IN_PLACE,&sane,1,MPI_CXX_BOOL,MPI_LAND,communicator);
          if (sane)
          {
              old_tmp_vertices = tmp_vertices;
              tmp_vertices = result_vertices;
              V = 1.0;
              for (int i = 0; i < balancer->getDimension(); ++i)
                  V *= ( tmp_vertices.at(1)[i] - tmp_vertices.at(0)[i] );
              old_tmp_work = tmp_work;
              tmp_work = rho * V;
          }
          else if (!sane || i == estimation_limit - 1)
          {
              balancer->getVertices() = old_tmp_vertices;
              workArray->at(0) = old_tmp_work;
              calculate_outline();
              i = estimation_limit;
          }
      }
 
      ((Staggered_LB<T,W>*)balancer.get())->setVertices(outline->data());
      */
      balancer->balance(loadbalancing_step);
      calculate_outline();
      break;
    case LB_t::FORCEBASED:
      balancer->balance(loadbalancing_step);
      calculate_outline();
      break;
#ifdef ALL_VORONOI_ACTIVE
    case LB_t::VORONOI:
      balancer->balance(loadbalancing_step);
      break;
#endif
    case LB_t::HISTOGRAM:
      balancer->balance(loadbalancing_step);
      calculate_outline();
      break;
    case LB_t::TENSOR_MAX:
      balancer->balance(loadbalancing_step);
      calculate_outline();
      break;
 
    default:
      throw InvalidArgumentException(__FILE__, __func__, __LINE__,
                                     "Unknown type of loadbalancing passed.");
    }
    loadbalancing_step++;
    return balancer->getVertices();
  }
 
  /**********************/
  /*  getter functions  */
  /**********************/

  std::vector<Point<T>> &getPrevVertices() {
    return balancer->getPrevVertices();
  };

  std::vector<Point<T>> &getVertices() { return balancer->getVertices(); };

  int getDimension() { return balancer->getDimension(); }

  void getWork(std::vector<W> &result) { result = balancer->getWork(); }

  W getWork() { return balancer->getWork().at(0); }

  std::vector<int> &getNeighbors() { return balancer->getNeighbors(); }

  std::vector<T> &getNeighborVertices() {
    return balancer->getNeighborVertices();
  }

  W getEfficiency(){
    return balancer->getEfficiency();
  }

  W getEstimatedEfficiency() {
    // currently only implemented for HISTOGRAM
    return balancer->getEstimatedEfficiency();
  }
 
  /**********************/
  /*  setter functions  */
  /**********************/

  void setSysSize(const std::vector<T> &sysSize) {
    balancer.get()->setSysSize(sysSize);
  }

  void setMethodData(const void *data) {
    balancer.get()->setAdditionalData(data);
  }

  void setProcGridParams(const std::vector<int> &loc,
                         const std::vector<int> &size) {
    // todo(s.schulz): We should probably throw if the dimension does not match
    procGridLoc.insert(procGridLoc.begin(), loc.begin(), loc.end());
    procGridDim.insert(procGridDim.begin(), size.begin(), size.end());
    procGridSet = true;
  }

  void setProcTag(int tag) { procTag = tag; }

  void setMinDomainSize(const std::vector<T> &minSize) {
    (balancer.get())->setMinDomainSize(minSize);
  }
 
#ifdef ALL_VTK_OUTPUT
  void printVTKoutlines(const int step) {
    // define grid points, i.e. vertices of local domain
    auto points = vtkSmartPointer<vtkPoints>::New();
    // allocate space for eight points (describing the cuboid around the domain)
    points->Allocate(8 * balancer->getDimension());
 
    int n_ranks;
    int local_rank;
 
    if (!vtk_init) {
      controller = vtkMPIController::New();
      controller->Initialize();
      controller->SetGlobalController(controller);
      vtk_init = true;
    }
 
    MPI_Comm_rank(communicator, &local_rank);
    MPI_Comm_size(communicator, &n_ranks);
 
    std::vector<Point<W>> tmp_outline(2);
    std::vector<W> tmp_0(
        {outline.at(0)[0], outline.at(0)[1], outline.at(0)[2]});
    std::vector<W> tmp_1(
        {outline.at(1)[0], outline.at(1)[1], outline.at(1)[2]});
    tmp_outline.at(0) = Point<W>(tmp_0);
    tmp_outline.at(1) = Point<W>(tmp_1);
 
    for (auto z = 0; z <= 1; ++z)
      for (auto y = 0; y <= 1; ++y)
        for (auto x = 0; x <= 1; ++x) {
          points->InsertPoint(x + 2 * y + 4 * z, tmp_outline.at(x)[0],
                              tmp_outline.at(y)[1], tmp_outline.at(z)[2]);
        }
 
    auto hexa = vtkSmartPointer<vtkVoxel>::New();
    for (int i = 0; i < 8; ++i)
      hexa->GetPointIds()->SetId(i, i);
 
    auto cellArray = vtkSmartPointer<vtkCellArray>::New();
    cellArray->InsertNextCell(hexa);
 
    // define work array (length = 1)
    auto work = vtkSmartPointer<vtkFloatArray>::New();
    work->SetNumberOfComponents(1);
    work->SetNumberOfTuples(1);
    work->SetName("work");
    W total_work = std::accumulate(balancer->getWork().begin(),
                                   balancer->getWork().end(), (W)0);
    work->SetValue(0, total_work);
 
    // define rank array (length = 4, x,y,z, rank)
    int rank = 0;
    MPI_Comm_rank(communicator, &rank);
    auto coords = vtkSmartPointer<vtkFloatArray>::New();
    coords->SetNumberOfComponents(3);
    coords->SetNumberOfTuples(1);
    coords->SetName("coords");
    coords->SetValue(0, procGridLoc.at(0));
    coords->SetValue(1, procGridLoc.at(1));
    coords->SetValue(2, procGridLoc.at(2));
 
    auto rnk = vtkSmartPointer<vtkFloatArray>::New();
    rnk->SetNumberOfComponents(1);
    rnk->SetNumberOfTuples(1);
    rnk->SetName("MPI rank");
    rnk->SetValue(0, rank);
 
    // define tag array (length = 1)
    auto tag = vtkSmartPointer<vtkIntArray>::New();
    tag->SetNumberOfComponents(1);
    tag->SetNumberOfTuples(1);
    tag->SetName("tag");
    tag->SetValue(0, procTag);
 
    // determine extent of system
    W known_unused global_extent[6];
    W local_min[3];
    W local_max[3];
    W global_min[3];
    W global_max[3];
    W width[3];
    W volume = (W)1.0;
 
    for (int i = 0; i < 3; ++i) {
      local_min[i] = outline.at(0)[i];
      local_max[i] = outline.at(1)[i];
      width[i] = local_max[i] - local_min[i];
      volume *= width[i];
    }
 
    W surface = (W)2.0 * width[0] * width[1] + (W)2.0 * width[1] * width[2] +
                (W)2.0 * width[0] * width[2];
 
    auto expanse = vtkSmartPointer<vtkFloatArray>::New();
    expanse->SetNumberOfComponents(6);
    expanse->SetNumberOfTuples(1);
    expanse->SetName("expanse");
    expanse->SetValue(0, width[0]);
    expanse->SetValue(1, width[0]);
    expanse->SetValue(2, width[0]);
    expanse->SetValue(3, volume);
    expanse->SetValue(4, surface);
    expanse->SetValue(5, surface / volume);
 
    MPI_Allreduce(&local_min, &global_min, 3, MPIDataTypeW, MPI_MIN,
                  communicator);
    MPI_Allreduce(&local_max, &global_max, 3, MPIDataTypeW, MPI_MAX,
                  communicator);
 
    for (int i = 0; i < 3; ++i) {
      global_extent[2 * i] = global_min[i];
      global_extent[2 * i + 1] = global_max[i];
    }
 
    // create a structured grid and assign data to it
    auto unstructuredGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
    unstructuredGrid->SetPoints(points);
    unstructuredGrid->SetCells(VTK_VOXEL, cellArray);
    unstructuredGrid->GetCellData()->AddArray(work);
    unstructuredGrid->GetCellData()->AddArray(expanse);
    unstructuredGrid->GetCellData()->AddArray(rnk);
    unstructuredGrid->GetCellData()->AddArray(coords);
    unstructuredGrid->GetCellData()->AddArray(tag);
    
    //DEBUG 
    /*
    createDirectory("vtk_outline_debug");
    std::ostringstream ss_local_debug;
    ss_local_debug << "vtk_outline_debug/ALL_vtk_outline_" << std::setw(7)
             << std::setfill('0') << step << "_" << local_rank << ".vtk";
 
    auto debug_writer = vtkSmartPointer<vtkUnstructuredGridWriter>::New();
    debug_writer->SetFileName(ss_local_debug.str().c_str());
    debug_writer->SetInputData(unstructuredGrid);
    debug_writer->Write();    
    */
    //ENDDEBUG
 
    createDirectory("vtk_outline");
    std::ostringstream ss_local;
    ss_local << "vtk_outline/ALL_vtk_outline_" << std::setw(7)
             << std::setfill('0') << step << "_" << local_rank << ".vtu";
 
    //DEBUG 
    /*
    createDirectory("vtk_outline_debug");
    std::ostringstream ss_local_debug;
    ss_local_debug << "vtk_outline_debug/ALL_vtk_outline_" << std::setw(7)
             << std::setfill('0') << step << "_" << local_rank << ".vtk";
 
    auto debug_writer = vtkSmartPointer<vtkUnstructuredGridWriter>::New();
    debug_writer->SetFileName(ss_local_debug.str().c_str());
    debug_writer->SetInputData(unstructuredGrid);
    debug_writer->Write();    
    */
    //ENDDEBUG
 
    auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
    writer->SetInputData(unstructuredGrid);
    writer->SetFileName(ss_local.str().c_str());
    //writer->SetDataModeToAscii(); //Causes Segmentation fault. Was used for debugging and is replaced by vtkUnstructuredGridWriter
    writer->SetDataModeToBinary();
    writer->Write();
 
    // if (local_rank == 0)
    //{
    std::ostringstream ss_para;
    ss_para << "vtk_outline/ALL_vtk_outline_" << std::setw(7)
            << std::setfill('0') << step << ".pvtu";
 
    // create the parallel writer
    auto parallel_writer =
        vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New();
    parallel_writer->SetFileName(ss_para.str().c_str());
    parallel_writer->SetNumberOfPieces(n_ranks);
    parallel_writer->SetStartPiece(local_rank);
    parallel_writer->SetEndPiece(local_rank);
    parallel_writer->SetInputData(unstructuredGrid);
    //parallel_writer->SetDataModeToAscii(); //Causes Segmentation fault. Was used for debugging and is replaced by vtkUnstructuredGridWriter
    parallel_writer->SetDataModeToBinary();
    parallel_writer->Write();
    //}
  }

  void printVTKvertices(const int step) {
    int n_ranks;
    int local_rank;
 
    MPI_Comm_rank(communicator, &local_rank);
    MPI_Comm_size(communicator, &n_ranks);
 
    // local vertices
    // (vertices + work)
    T local_vertices[nVertices * balancer->getDimension() + 1];
 
    for (int v = 0; v < nVertices; ++v) {
      for (int d = 0; d < balancer->getDimension(); ++d) {
        local_vertices[v * balancer->getDimension() + d] =
            balancer->getVertices().at(v)[d];
      }
    }
    local_vertices[nVertices * balancer->getDimension()] =
        (T)balancer->getWork().at(0);
 
    /*
    T *global_vertices;
    if (local_rank == 0) {
      global_vertices =
          new T[n_ranks * (nVertices * balancer->getDimension() + 1)];
    }
    */
    T global_vertices[n_ranks * (nVertices * balancer->getDimension() + 1)];
 
    // collect all works and vertices on a single process
    MPI_Gather(local_vertices, nVertices * balancer->getDimension() + 1,
               MPIDataTypeT, global_vertices,
               nVertices * balancer->getDimension() + 1, MPIDataTypeT, 0,
               communicator);
 
    if (local_rank == 0) {
      auto vtkpoints = vtkSmartPointer<vtkPoints>::New();
#ifdef VTK_CELL_ARRAY_V2
      vtkNew<vtkUnstructuredGrid> unstructuredGrid;
      unstructuredGrid->Allocate(n_ranks + 10);
#else
      auto unstructuredGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
#endif
 
      // enter vertices into unstructured grid
      for (int i = 0; i < n_ranks; ++i) {
        for (int v = 0; v < nVertices; ++v) {
          vtkpoints->InsertNextPoint(
              global_vertices[i * (nVertices * balancer->getDimension() + 1) +
                              v * balancer->getDimension() + 0],
              global_vertices[i * (nVertices * balancer->getDimension() + 1) +
                              v * balancer->getDimension() + 1],
              global_vertices[i * (nVertices * balancer->getDimension() + 1) +
                              v * balancer->getDimension() + 2]);
        }
      }
      unstructuredGrid->SetPoints(vtkpoints);
 
      // data arrays for work and cell id
      auto work = vtkSmartPointer<vtkFloatArray>::New();
      auto cell = vtkSmartPointer<vtkFloatArray>::New();
      work->SetNumberOfComponents(1);
      work->SetNumberOfTuples(n_ranks);
      work->SetName("work");
      cell->SetNumberOfComponents(1);
      cell->SetNumberOfTuples(n_ranks);
      cell->SetName("cell id");
 
      for (int n = 0; n < n_ranks; ++n) {
 
#ifdef VTK_CELL_ARRAY_V2
        // define grid points, i.e. vertices of local domain
        vtkIdType pointIds[8] = {8 * n + 0, 8 * n + 1, 8 * n + 2, 8 * n + 3,
                                 8 * n + 4, 8 * n + 5, 8 * n + 6, 8 * n + 7};
 
        vtkIdType faces[48] = {
            3,         8 * n + 0, 8 * n + 2, 8 * n + 1, 
            3,         8 * n + 1, 8 * n + 2, 8 * n + 3, 
            3,         8 * n + 0, 8 * n + 4, 8 * n + 2, 
            3,         8 * n + 2, 8 * n + 4, 8 * n + 6, 
            3,         8 * n + 2, 8 * n + 6, 8 * n + 3, 
            3,         8 * n + 3, 8 * n + 6, 8 * n + 7, 
            3,         8 * n + 1, 8 * n + 5, 8 * n + 3, 
            3,         8 * n + 3, 8 * n + 5, 8 * n + 7, 
            3,         8 * n + 0, 8 * n + 4, 8 * n + 1, 
            3,         8 * n + 1, 8 * n + 4, 8 * n + 5, 
            3,         8 * n + 4, 8 * n + 6, 8 * n + 5, 
            3,         8 * n + 5, 8 * n + 6, 8 * n + 7};
 
        unstructuredGrid->InsertNextCell(VTK_POLYHEDRON, 8, pointIds, 12,
                                         faces);
#else
        // define grid points, i.e. vertices of local domain
        vtkIdType pointIds[8] = {8 * n + 0, 8 * n + 1, 8 * n + 2, 8 * n + 3,
                                 8 * n + 4, 8 * n + 5, 8 * n + 6, 8 * n + 7};
 
        auto faces = vtkSmartPointer<vtkCellArray>::New();
        // setup faces of polyhedron
        vtkIdType f0[3] = {8 * n + 0, 8 * n + 2, 8 * n + 1};
        vtkIdType f1[3] = {8 * n + 1, 8 * n + 2, 8 * n + 3};
 
        vtkIdType f2[3] = {8 * n + 0, 8 * n + 4, 8 * n + 2};
        vtkIdType f3[3] = {8 * n + 2, 8 * n + 4, 8 * n + 6};
 
        vtkIdType f4[3] = {8 * n + 2, 8 * n + 6, 8 * n + 3};
        vtkIdType f5[3] = {8 * n + 3, 8 * n + 6, 8 * n + 7};
 
        vtkIdType f6[3] = {8 * n + 1, 8 * n + 5, 8 * n + 3};
        vtkIdType f7[3] = {8 * n + 3, 8 * n + 5, 8 * n + 7};
 
        vtkIdType f8[3] = {8 * n + 0, 8 * n + 4, 8 * n + 1};
        vtkIdType f9[3] = {8 * n + 1, 8 * n + 4, 8 * n + 5};
 
        vtkIdType fa[3] = {8 * n + 4, 8 * n + 6, 8 * n + 5};
        vtkIdType fb[3] = {8 * n + 5, 8 * n + 6, 8 * n + 7};
 
        faces->InsertNextCell(3, f0);
        faces->InsertNextCell(3, f1);
        faces->InsertNextCell(3, f2);
        faces->InsertNextCell(3, f3);
        faces->InsertNextCell(3, f4);
        faces->InsertNextCell(3, f5);
        faces->InsertNextCell(3, f6);
        faces->InsertNextCell(3, f7);
        faces->InsertNextCell(3, f8);
        faces->InsertNextCell(3, f9);
        faces->InsertNextCell(3, fa);
        faces->InsertNextCell(3, fb);
 
        unstructuredGrid->InsertNextCell(VTK_POLYHEDRON, 8, pointIds, 12,
                                         faces->GetPointer());
#endif
        work->SetValue(
            n, global_vertices[n * (nVertices * balancer->getDimension() + 1) +
                               8 * balancer->getDimension()]);
        cell->SetValue(n, (T)n);
      }
      unstructuredGrid->GetCellData()->AddArray(work);
      unstructuredGrid->GetCellData()->AddArray(cell);
 
      createDirectory("vtk_vertices");
      std::ostringstream filename;
      filename << "vtk_vertices/ALL_vtk_vertices_" << std::setw(7)
               << std::setfill('0') << step << ".vtu";
      auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
      writer->SetInputData(unstructuredGrid);
      writer->SetFileName(filename.str().c_str());
      writer->SetDataModeToAscii();
      // writer->SetDataModeToBinary();
      writer->Write();
 
      // delete[] global_vertices;
    }
  }
#endif
 
private:
  LB_t method;

  std::vector<Point<T>> outline;

  std::unique_ptr<LB<T, W>> balancer;

  std::vector<W> workArray;

  MPI_Comm communicator;

  int nVertices;

  int nNeighbors;

  void calculate_outline() {
    // calculation only possible if there are vertices defining the domain
    if (balancer->getVertices().size() > 0) {
      outline.resize(2);
      // setup the outline with the first point
      for (int i = 0; i < balancer->getDimension(); ++i) {
        outline.at(0) = balancer->getVertices().at(0);
        outline.at(1) = balancer->getVertices().at(0);
      }
      // compare value of each outline point with all vertices to find the
      // maximum extend of the domain in each dimension
      for (int i = 1; i < (int)balancer->getVertices().size(); ++i) {
        Point<T> p = balancer->getVertices().at(i);
        for (int j = 0; j < balancer->getDimension(); ++j) {
          outline.at(0)[j] = std::min(outline.at(0)[j], p[j]);
          outline.at(1)[j] = std::max(outline.at(1)[j], p[j]);
        }
      }
    }
  }

  void createDirectory(const char *filename) {
    errno = 0;
    mode_t mode = S_IRWXU | S_IRWXG; // rwx for user and group
    if (mkdir(filename, mode)) {
      if (errno != EEXIST) {
        throw FilesystemErrorException(__FILE__, __func__, __LINE__,
                                       "Could not create output directory.");
      }
    }
    // check permission in case directory already existed, but had wrong
    // permissions
    struct stat attr;
    stat(filename, &attr);
    if ((attr.st_mode & S_IRWXU) != S_IRWXU) {
      throw FilesystemErrorException(
          __FILE__, __func__, __LINE__,
          "Possibly already existing directory does not have correct "
          "permissions (rwx) for user");
    }
  }

  int estimation_limit;

  std::vector<int> procGridLoc;
  std::vector<int> procGridDim;
  bool procGridSet;
  int procTag;

  MPI_Datatype MPIDataTypeT;

  MPI_Datatype MPIDataTypeW;

  int loadbalancing_step;
 
#ifdef ALL_VTK_OUTPUT
  bool vtk_init;

  vtkMPIController *controller;

  vtkMPIController *controller_vertices;
#endif
};
 
} // namespace ALL
 
#endif
 
// VIM modeline for indentation
// vim: sw=2 et