From 2d13cdee386ad5bbaceead4c2c774e2072a5d934 Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Fri, 2 Feb 2024 21:55:07 +0100
Subject: [PATCH] First working verison

---
 .../NonUniformGridCPU/NonUniformGridCPU.cpp   |   4 +-
 .../simulation_setup/benchmark_configs.py     |   9 +-
 .../NonUniformGridGPU/NonUniformGridGPU.cpp   |   2 +-
 .../UniformGridCPU/UniformGridCPU.cpp         |   2 +-
 .../UniformGridGPU/UniformGridGPU.cpp         |   2 +-
 src/vtk/VTKOutput.cpp                         | 336 ++++++++++--------
 src/vtk/VTKOutput.h                           |   9 +-
 7 files changed, 202 insertions(+), 162 deletions(-)

diff --git a/apps/benchmarks/NonUniformGridCPU/NonUniformGridCPU.cpp b/apps/benchmarks/NonUniformGridCPU/NonUniformGridCPU.cpp
index 6476d82ec..ad4450721 100644
--- a/apps/benchmarks/NonUniformGridCPU/NonUniformGridCPU.cpp
+++ b/apps/benchmarks/NonUniformGridCPU/NonUniformGridCPU.cpp
@@ -174,8 +174,8 @@ int main(int argc, char** argv)
       if (vtkWriteFrequency > 0)
       {
          auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, "vtk_out",
-                                                         "simulation_step", false, true, true, false, 0);
-         auto velWriter = make_shared< field::VTKWriter< VelocityField_T > >(velFieldID, "vel");
+                                                         "simulation_step", false, false, true, false, 0);
+         auto velWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(velFieldID, "vel");
          vtkOutput->addCellDataWriter(velWriter);
 
          vtkOutput->addBeforeFunction([&]() {
diff --git a/apps/benchmarks/NonUniformGridCPU/simulation_setup/benchmark_configs.py b/apps/benchmarks/NonUniformGridCPU/simulation_setup/benchmark_configs.py
index 8c99c62c2..e4a391491 100644
--- a/apps/benchmarks/NonUniformGridCPU/simulation_setup/benchmark_configs.py
+++ b/apps/benchmarks/NonUniformGridCPU/simulation_setup/benchmark_configs.py
@@ -6,6 +6,7 @@ import sys
 
 DB_FILE = os.environ.get('DB_FILE', "cpu_benchmark.sqlite3")
 
+
 class Scenario:
     def __init__(self,
                  domain_size=(64, 64, 64),
@@ -63,7 +64,7 @@ class Scenario:
             }
         }
 
-        if(print_dict):
+        if (print_dict):
             wlb.log_info_on_root("Scenario:\n" + pformat(config_dict))
 
         return config_dict
@@ -100,8 +101,8 @@ def validation_run():
     """Run with full periodic shear flow or boundary scenario (ldc) to check if the code works"""
     wlb.log_info_on_root("Validation run")
 
-    domain_size = (96, 96, 96)
-    cells_per_block = (32, 32, 32)
+    domain_size = (2, 4, 2)
+    cells_per_block = (2, 2, 2)
 
     root_blocks = tuple([d // c for d, c in zip(domain_size, cells_per_block)])
 
@@ -117,6 +118,7 @@ def validation_run():
                         write_setup_vtk=True)
     scenarios.add(scenario)
 
+
 def scaling():
     wlb.log_info_on_root("Running scaling benchmark...")
 
@@ -134,5 +136,6 @@ def scaling():
                         timesteps=10)
     scenarios.add(scenario)
 
+
 validation_run()
 # scaling()
diff --git a/apps/benchmarks/NonUniformGridGPU/NonUniformGridGPU.cpp b/apps/benchmarks/NonUniformGridGPU/NonUniformGridGPU.cpp
index 919755d6d..e76c0b511 100644
--- a/apps/benchmarks/NonUniformGridGPU/NonUniformGridGPU.cpp
+++ b/apps/benchmarks/NonUniformGridGPU/NonUniformGridGPU.cpp
@@ -243,7 +243,7 @@ int main(int argc, char** argv)
       {
          auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, "vtk_out",
                                                          "simulation_step", false, true, true, false, 0);
-         auto velWriter = make_shared< field::VTKWriter< VelocityField_T > >(velFieldCpuID, "vel");
+         auto velWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(velFieldCpuID, "vel");
          vtkOutput->addCellDataWriter(velWriter);
 
          vtkOutput->addBeforeFunction([&]() {
diff --git a/apps/benchmarks/UniformGridCPU/UniformGridCPU.cpp b/apps/benchmarks/UniformGridCPU/UniformGridCPU.cpp
index dfcd22a87..4674cfae9 100644
--- a/apps/benchmarks/UniformGridCPU/UniformGridCPU.cpp
+++ b/apps/benchmarks/UniformGridCPU/UniformGridCPU.cpp
@@ -174,7 +174,7 @@ int main(int argc, char** argv)
       {
          auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, "vtk_out",
                                                          "simulation_step", false, true, true, false, 0);
-         auto velWriter = make_shared< field::VTKWriter< VelocityField_T > >(velFieldId, "vel");
+         auto velWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(velFieldId, "vel");
          vtkOutput->addCellDataWriter(velWriter);
 
          vtkOutput->addBeforeFunction([&]() {
diff --git a/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp b/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp
index fdc8969d6..91b7a0210 100644
--- a/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp
+++ b/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp
@@ -205,7 +205,7 @@ int main(int argc, char** argv)
       {
          auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, "vtk_out",
                                                          "simulation_step", false, true, true, false, 0);
-         auto velWriter = make_shared< field::VTKWriter< VelocityField_T > >(velFieldCpuID, "vel");
+         auto velWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(velFieldCpuID, "vel");
          vtkOutput->addCellDataWriter(velWriter);
 
          vtkOutput->addBeforeFunction([&]() {
diff --git a/src/vtk/VTKOutput.cpp b/src/vtk/VTKOutput.cpp
index 345d4fda3..c7b48af63 100644
--- a/src/vtk/VTKOutput.cpp
+++ b/src/vtk/VTKOutput.cpp
@@ -983,14 +983,6 @@ void VTKOutput::computeVTUCells( const IBlock& block, CellVector & cellsOut ) co
 void VTKOutput::writeBlocks( const std::string& path, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates )
 {
    WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
-
-   std::vector< const IBlock* > blocks;
-   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block )
-   {
-      if( selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
-         blocks.push_back( block.get() );
-   }
-
    if( !configured_ ) {
       if( !forcePVTU_ && cellInclusionFunctions_.empty() && cellExclusionFunctions_.empty() &&
           blockStorage_->getNumberOfLevels() == 1 && ghostLayers_ == 0 ) // uniform data -> vti
@@ -1000,41 +992,42 @@ void VTKOutput::writeBlocks( const std::string& path, const Set<SUID>& requiredS
       configured_ = true;
    }
 
-   for( auto it = blocks.begin(); it != blocks.end(); ++it )
+   if(uniformGrid_)
    {
-      WALBERLA_ASSERT_NOT_NULLPTR( *it );
-      const IBlock& block = **it;
-
-      std::ostringstream file;
-      file << path << "/block [" << block.getId() << "].";
-
-      if( uniformGrid_ ) // uniform data -> vti
+      for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block )
       {
-         file << "vti";
-         std::ofstream ofs( file.str().c_str()  );
-         if( samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0) )
-            writeVTI( ofs, block );
-         else
-            writeVTI_sampling( ofs, block );
-         ofs.close();
-      }
-      else // unstructured data -> vtu
-      {
-         CellVector cells; // cells to be written to file
-         computeVTUCells( block, cells );
+         if( !selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
+            continue;
+
+         WALBERLA_ASSERT_NOT_NULLPTR(*block);
+         std::ostringstream file;
+         file << path << "/block [" << block->getId() << "].";
 
-         if( !cells.empty() )
+         if (uniformGrid_) // uniform data -> vti
          {
-            file << "vtu";
-            std::ofstream ofs( file.str().c_str()  );
-            if( samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0) )
-               writeVTU( ofs, block, cells );
+            file << "vti";
+            std::ofstream ofs(file.str().c_str());
+            if (samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0))
+               writeVTI(ofs, *block);
             else
-               writeVTU_sampling( ofs, block, cells );
+               writeVTI_sampling(ofs, *block);
             ofs.close();
          }
       }
    }
+   else
+   {
+      const int process = MPIManager::instance()->rank();
+      std::ostringstream file;
+      file << path << "/dataRank[" << process << "].";
+      file << "vtu";
+      std::ofstream ofs( file.str().c_str()  );
+      if( samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0) )
+         writeVTU( ofs, requiredStates, incompatibleStates );
+      else
+         writeVTU_sampling( ofs, requiredStates, incompatibleStates );
+      ofs.close();
+   }
 }
 
 
@@ -1042,49 +1035,36 @@ void VTKOutput::writeBlocks( const std::string& path, const Set<SUID>& requiredS
 void VTKOutput::writeBlockPieces( std::ostream & oss, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates )
 {
    WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
-
-   std::vector< const IBlock* > blocks;
-   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block )
-   {
-      if( selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
-         blocks.push_back( block.get() );
-   }
-
    if( !configured_ ) {
       if( !forcePVTU_ && cellInclusionFunctions_.empty() && cellExclusionFunctions_.empty() &&
-         blockStorage_->getNumberOfLevels() == 1 && ghostLayers_ == 0 ) // uniform data -> vti
+          blockStorage_->getNumberOfLevels() == 1 && ghostLayers_ == 0 ) // uniform data -> vti
       {
          uniformGrid_ = true;
       }
       configured_ = true;
    }
 
-   for( auto it = blocks.begin(); it != blocks.end(); ++it )
+   if(uniformGrid_)
    {
-      WALBERLA_ASSERT_NOT_NULLPTR( *it );
-      const IBlock& block = **it;
-
-      if( uniformGrid_ ) // uniform data -> vti
+      for (auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block)
       {
-         if( samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0) )
-            writeVTIPiece( oss, block );
-         else
-            writeVTIPiece_sampling( oss, block );
-      }
-      else // unstructured data -> vtu
-      {
-         CellVector cells; // cells to be written to file
-         computeVTUCells( block, cells );
+         if (!selectable::isSetSelected(uid::globalState() + block->getState(), requiredStates, incompatibleStates))
+            continue;
+         WALBERLA_ASSERT_NOT_NULLPTR(*it);
 
-         if( !cells.empty() )
-         {
-            if( samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0) )
-               writeVTUPiece( oss, block, cells );
-            else
-               writeVTUPiece_sampling( oss, block, cells );
-         }
+         if (samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0))
+            writeVTIPiece(oss, *block);
+         else
+            writeVTIPiece_sampling(oss, *block);
       }
    }
+   else // unstructured data -> vtu
+   {
+      if( samplingDx_ <= real_c(0) || samplingDy_ <= real_c(0) || samplingDz_ <= real_c(0) )
+         writeVTUPiece( oss, requiredStates, incompatibleStates );
+      else
+         writeVTUPiece_sampling( oss, requiredStates, incompatibleStates );
+   }
 }
 
 
@@ -1196,13 +1176,13 @@ void VTKOutput::writeVTIPiece_sampling( std::ostream& ofs, const IBlock& block )
 
 
 
-void VTKOutput::writeVTU( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const
+void VTKOutput::writeVTU( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const
 {
    ofs << "<?xml version=\"1.0\"?>\n"
        << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
        << " <UnstructuredGrid>\n";
 
-   writeVTUPiece( ofs, block, cells );
+   writeVTUPiece( ofs, requiredStates, incompatibleStates );
 
    ofs << " </UnstructuredGrid>\n"
        << "</VTKFile>" << std::endl;
@@ -1210,77 +1190,87 @@ void VTKOutput::writeVTU( std::ostream& ofs, const IBlock& block, const CellVect
 
 
 
-void VTKOutput::writeVTUPiece( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const
+void VTKOutput::writeVTUPiece( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const
 {
    WALBERLA_ASSERT_NOT_NULLPTR(blockStorage_);
-
-   // setting up vertex-index mapping --->
+   const uint_t finestLevel = blockStorage_->getNumberOfLevels() - 1;
 
    std::map< Vertex, Index, VertexCompare > vimap; // vertex<->index map
    std::vector< VertexCoord >               vc;    // vertex coordinates
    std::vector< Index >                     ci;    // ci[0] to ci[7]: indices for cell number one, ci[8] to ci[15]: ...
+   uint_t numberOfCells = 0;
 
-   for (auto cell = cells.begin(); cell != cells.end(); ++cell)
+   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block )
    {
-      AABB aabb;
-      blockStorage_->getBlockLocalCellAABB(block, *cell, aabb);
+      if( !selectable::isSetSelected( uid::globalState() + block->getState(), requiredStates, incompatibleStates ) )
+         continue;
 
-      for (cell_idx_t z = 0; z != 2; ++z) {
-         for (cell_idx_t y = 0; y != 2; ++y) {
-            for (cell_idx_t x = 0; x != 2; ++x)
-            {
-               Vertex v(cell->x() + x, cell->y() + y, cell->z() + z);
+      // TODO why ????
+      // CellVector cells; // cells to be written to file
+      // computeVTUCells( *block, cells );
 
-               auto mapping = vimap.find(v);
-               if (mapping != vimap.end()) // vertex already exists
-               {
-                  ci.push_back(mapping->second);
-               }
-               else // new vertex
+      const uint_t level = blockStorage_->getLevel(*block);
+      const cell_idx_t factorToFinest = 1 << (finestLevel - level);
+      const CellInterval cells = blockStorage_->getBlockCellBB(*block); //  These are global cells
+
+      for (auto cell = cells.begin(); cell != cells.end(); ++cell)
+      {
+         numberOfCells++;
+         const AABB aabb = blockStorage_->getCellAABB(*cell, level);
+         for (cell_idx_t z = 0; z != 2; ++z) {
+            for (cell_idx_t y = 0; y != 2; ++y) {
+               for (cell_idx_t x = 0; x != 2; ++x)
                {
-                  vimap[v] = numeric_cast< Index >(vc.size());
-                  ci.push_back(numeric_cast< Index >(vc.size()));
-                  vc.emplace_back((x == 0) ? aabb.xMin() : aabb.xMax(),
-                     (y == 0) ? aabb.yMin() : aabb.yMax(),
-                     (z == 0) ? aabb.zMin() : aabb.zMax());
+                  const Vertex v((cell->x() + x) * factorToFinest, (cell->y() + y) * factorToFinest, (cell->z() + z) * factorToFinest);
+                  auto mapping = vimap.find(v);
+                  if (mapping != vimap.end()) // vertex already exists
+                  {
+                     ci.push_back(mapping->second);
+                  }
+                  else // new vertex
+                  {
+                     vimap[v] = numeric_cast< Index >(vc.size());
+                     ci.push_back(numeric_cast< Index >(vc.size()));
+                     vc.emplace_back((x == 0) ? aabb.xMin() : aabb.xMax(),
+                                     (y == 0) ? aabb.yMin() : aabb.yMax(),
+                                     (z == 0) ? aabb.zMin() : aabb.zMax());
+                  }
                }
             }
          }
       }
    }
-
-   WALBERLA_ASSERT_EQUAL(ci.size(), 8 * cells.size());
-
    // <--- setting up vertex-index mapping
 
-   writeVTUHeaderPiece(ofs, uint_c(cells.size()), vc, ci);
+   writeVTUHeaderPiece(ofs, numberOfCells, vc, ci);
 
    ofs << "   <CellData>\n";
 
-   if (ghostLayers_ > 0)
-   {
-      ofs << "    <DataArray type=\"" << vtk::typeToString< uint8_t >()
-          << "\" Name=\"vtkGhostLevels\" NumberOfComponents=\"1\" format=\"" << format_ << "\">\n"
-          << "     ";
-
-      if (binary_)
-      {
-         Base64Writer base64;
-         for (auto cell = cells.begin(); cell != cells.end(); ++cell)
-            base64 << ghostLayerNr(block, cell->x(), cell->y(), cell->z());
-         base64.toStream(ofs);
-      }
-      else
-      {
-         for (auto cell = cells.begin(); cell != cells.end(); ++cell)
-            ofs << uint_c(ghostLayerNr(block, cell->x(), cell->y(), cell->z())) << " ";
-         ofs << "\n";
-      }
-
-      ofs << "    </DataArray>\n";
-   }
+   // TODO think about ghostlayer
+//   if (ghostLayers_ > 0)
+//   {
+//      ofs << "    <DataArray type=\"" << vtk::typeToString< uint8_t >()
+//          << "\" Name=\"vtkGhostLevels\" NumberOfComponents=\"1\" format=\"" << format_ << "\">\n"
+//          << "     ";
+//
+//      if (binary_)
+//      {
+//         Base64Writer base64;
+//         for (auto cell = cells.begin(); cell != cells.end(); ++cell)
+//            base64 << ghostLayerNr(block, cell->x(), cell->y(), cell->z());
+//         base64.toStream(ofs);
+//      }
+//      else
+//      {
+//         for (auto cell = cells.begin(); cell != cells.end(); ++cell)
+//            ofs << uint_c(ghostLayerNr(block, cell->x(), cell->y(), cell->z())) << " ";
+//         ofs << "\n";
+//      }
+//
+//      ofs << "    </DataArray>\n";
+//   }
 
-   writeCellData(ofs, block, cells);
+   writeCellData(ofs, requiredStates, incompatibleStates);
 
    ofs << "   </CellData>\n"
        << "  </Piece>\n";
@@ -1288,13 +1278,13 @@ void VTKOutput::writeVTUPiece( std::ostream& ofs, const IBlock& block, const Cel
 
 
 
-void VTKOutput::writeVTU_sampling( std::ostream& ofs, const IBlock& block, const CellVector& localCells ) const
+void VTKOutput::writeVTU_sampling( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const
 {
    ofs << "<?xml version=\"1.0\"?>\n"
        << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"" << endianness_ << "\">\n"
        << " <UnstructuredGrid>\n";
 
-   writeVTUPiece_sampling( ofs, block, localCells );
+   writeVTUPiece_sampling( ofs, requiredStates, incompatibleStates );
 
    ofs << " </UnstructuredGrid>\n"
        << "</VTKFile>" << std::endl;
@@ -1303,51 +1293,56 @@ void VTKOutput::writeVTU_sampling( std::ostream& ofs, const IBlock& block, const
 
 
 
-void VTKOutput::writeVTUPiece_sampling(std::ostream& ofs, const IBlock& block, const CellVector& localCells) const
+void VTKOutput::writeVTUPiece_sampling(std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const
 {
-   std::vector< SamplingCell > cells = getSamplingCells(block, localCells);
-
-   // setting up vertex-index mapping --->
 
    std::map< Vertex, Index, VertexCompare > vimap; // vertex<->index map
    std::vector< VertexCoord >               vc;    // vertex coordinates
    std::vector< Index >                     ci;    // ci[0] to ci[7]: indices for cell number one, ci[8] to ci[15]: ...
+   uint_t numberOfCells = 0;
 
-   for (auto cell = cells.begin(); cell != cells.end(); ++cell) {
-      for (cell_idx_t z = 0; z != 2; ++z) {
-         for (cell_idx_t y = 0; y != 2; ++y) {
-            for (cell_idx_t x = 0; x != 2; ++x)
-            {
-               Vertex v(cell->coordinates_.x() + x, cell->coordinates_.y() + y, cell->coordinates_.z() + z);
+   for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block )
+   {
+      if (!selectable::isSetSelected(uid::globalState() + block->getState(), requiredStates, incompatibleStates))
+         continue;
 
-               auto mapping = vimap.find(v);
-               if (mapping != vimap.end()) // vertex already exists
-               {
-                  ci.push_back(mapping->second);
-               }
-               else // new vertex
+      CellVector localCells; // cells to be written to file
+      computeVTUCells( *block, localCells );
+      // setting up vertex-index mapping --->
+      std::vector< SamplingCell > cells = getSamplingCells(*block, localCells);
+      for (auto cell = cells.begin(); cell != cells.end(); ++cell)
+      {
+         numberOfCells++;
+         for (cell_idx_t z = 0; z != 2; ++z)
+         {
+            for (cell_idx_t y = 0; y != 2; ++y)
+            {
+               for (cell_idx_t x = 0; x != 2; ++x)
                {
-                  vimap[v] = numeric_cast< Index >(vc.size());
-                  ci.push_back(numeric_cast< Index >(vc.size()));
-                  vc.emplace_back((x == 0) ? cell->aabb_.xMin() : cell->aabb_.xMax(),
-                     (y == 0) ? cell->aabb_.yMin() : cell->aabb_.yMax(),
-                     (z == 0) ? cell->aabb_.zMin() : cell->aabb_.zMax());
+                  Vertex v(cell->coordinates_.x() + x, cell->coordinates_.y() + y, cell->coordinates_.z() + z);
+
+                  auto mapping = vimap.find(v);
+                  if (mapping != vimap.end()) // vertex already exists
+                  {
+                     ci.push_back(mapping->second);
+                  }
+                  else // new vertex
+                  {
+                     vimap[v] = numeric_cast< Index >(vc.size());
+                     ci.push_back(numeric_cast< Index >(vc.size()));
+                     vc.emplace_back((x == 0) ? cell->aabb_.xMin() : cell->aabb_.xMax(),
+                                     (y == 0) ? cell->aabb_.yMin() : cell->aabb_.yMax(),
+                                     (z == 0) ? cell->aabb_.zMin() : cell->aabb_.zMax());
+                  }
                }
             }
          }
       }
    }
 
-   WALBERLA_ASSERT_EQUAL(ci.size(), 8 * cells.size());
-
-   // <--- setting up vertex-index mapping
-
-   writeVTUHeaderPiece(ofs, uint_c(cells.size()), vc, ci);
-
+   writeVTUHeaderPiece(ofs, numberOfCells, vc, ci);
    ofs << "   <CellData>\n";
-
-   writeCellData(ofs, block, cells);
-
+   writeCellData(ofs, requiredStates, incompatibleStates);
    ofs << "   </CellData>\n"
       << "  </Piece>\n";
 }
@@ -1523,8 +1518,6 @@ std::vector< VTKOutput::SamplingCell > VTKOutput::getSamplingCells( const IBlock
    return samplingCells;
 }
 
-
-
 void VTKOutput::writeCellData( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const
 {
    WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
@@ -1534,7 +1527,7 @@ void VTKOutput::writeCellData( std::ostream& ofs, const IBlock& block, const Cel
       (*writer)->configure( block, *blockStorage_ );
 
       ofs << "    <DataArray type=\"" << (*writer)->typeString() << "\" Name=\"" << (*writer)->identifier()
-                                      << "\" NumberOfComponents=\"" << (*writer)->fSize() << "\" format=\"" << format_ << "\">\n";
+          << "\" NumberOfComponents=\"" << (*writer)->fSize() << "\" format=\"" << format_ << "\">\n";
 
       if( binary_ )
       {
@@ -1561,6 +1554,49 @@ void VTKOutput::writeCellData( std::ostream& ofs, const IBlock& block, const Cel
 }
 
 
+void VTKOutput::writeCellData( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const
+{
+   WALBERLA_ASSERT_NOT_NULLPTR( blockStorage_ );
+
+   for( auto writer = cellDataWriter_.begin(); writer != cellDataWriter_.end(); ++writer )
+   {
+      ofs << "    <DataArray type=\"" << (*writer)->typeString() << "\" Name=\"" << (*writer)->identifier()
+          << "\" NumberOfComponents=\"" << (*writer)->fSize() << "\" format=\"" << format_ << "\">\n";
+
+      for( auto block = blockStorage_->begin(); block != blockStorage_->end(); ++block )
+      {
+         if (!selectable::isSetSelected(uid::globalState() + block->getState(), requiredStates, incompatibleStates))
+            continue;
+
+         CellVector cells; // cells to be written to file
+         computeVTUCells(*block, cells);
+         (*writer)->configure( *block, *blockStorage_ );
+
+         if( binary_ )
+         {
+            Base64Writer base64;
+            for( auto cell = cells.begin(); cell != cells.end(); ++cell )
+               for( uint_t f = 0; f != (*writer)->fSize(); ++f )
+                  (*writer)->push( base64, cell->x(), cell->y(), cell->z(), cell_idx_c(f) );
+            ofs << "     "; base64.toStream( ofs );
+         }
+         else
+         {
+            for( auto cell = cells.begin(); cell != cells.end(); ++cell ) {
+               ofs << "     ";
+               for( uint_t f = 0; f != (*writer)->fSize(); ++f )
+               {
+                  (*writer)->push( ofs, cell->x(), cell->y(), cell->z(), cell_idx_c(f) );
+                  ofs << ( ( f == (*writer)->fSize() - 1 ) ? "\n" : " " );
+               }
+            }
+         }
+      }
+      ofs << "    </DataArray>\n";
+   }
+}
+
+
 
 void VTKOutput::writeCellData( std::ostream& ofs, const IBlock& block, const std::vector< SamplingCell >& cells ) const
 {
diff --git a/src/vtk/VTKOutput.h b/src/vtk/VTKOutput.h
index d90586b5c..03dea5375 100644
--- a/src/vtk/VTKOutput.h
+++ b/src/vtk/VTKOutput.h
@@ -241,11 +241,11 @@ private:
    void writeVTIPiece( std::ostream& ofs, const IBlock& block ) const;
    void writeVTIPiece_sampling( std::ostream& ofs, const IBlock& block ) const;
 
-   void writeVTU( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const;
-   void writeVTU_sampling( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const;
+   void writeVTU( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const;
+   void writeVTU_sampling( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const;
 
-   void writeVTUPiece(std::ostream& ofs, const IBlock& block, const CellVector& cells) const;
-   void writeVTUPiece_sampling(std::ostream& ofs, const IBlock& block, const CellVector& cells) const;
+   void writeVTUPiece( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const;
+   void writeVTUPiece_sampling(std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates) const;
 
    void writeVTUHeader( std::ofstream& ofs, const uint_t numberOfCells, const std::vector< VertexCoord > & vc, const std::vector< Index > & ci ) const;
    void writeVTUHeaderPiece (std::ostream& ofs, const uint_t numberOfCells, const std::vector< VertexCoord > & vc, const std::vector< Index > & ci) const;
@@ -255,6 +255,7 @@ private:
    std::vector< SamplingCell > getSamplingCells( const IBlock& block, const CellVector& cells ) const;
 
    void writeCellData( std::ostream& ofs, const IBlock& block, const CellVector& cells ) const;
+   void writeCellData( std::ostream& ofs, const Set<SUID>& requiredStates, const Set<SUID>& incompatibleStates ) const;
    void writeCellData( std::ostream& ofs, const IBlock& block, const std::vector< SamplingCell >& cells ) const;
 
    void writePVD();
-- 
GitLab