OctreeLogger.h

Go to the documentation of this file.

 
#pragma once
 
#include <array>
#include <cstring>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
 
#include "autopas/containers/octree/OctreeNodeInterface.h"
#include "autopas/containers/octree/OctreeNodeWrapper.h"
 
namespace autopas {
template <typename Particle_T>
class OctreeLogger {
 public:
  explicit OctreeLogger() = default;
 
  ~OctreeLogger() = default;
 
  void logTree(OctreeNodeWrapper<Particle_T> *wrapper) { logTree(wrapper->getRaw()); }
 
  void logTree(OctreeNodeInterface<Particle_T> *root) {
    // Load the leaf boxes
    using Position = std::array<double, 3>;
    using Box = std::pair<Position, Position>;
    std::vector<Box> boxes;
    root->appendAllLeafBoxes(boxes);
    auto boxCount = boxes.size();
    auto pointCount = 8 * boxCount;
 
    // Open the VTK file
    char filename[256] = {0};
    snprintf(filename, sizeof(filename), "octree_%d.vtk", iteration++);
    std::ofstream vtkFile;
    vtkFile.open(filename);
 
    if (not vtkFile.is_open()) {
      throw std::runtime_error("OctreeLogger::logTree(): Failed to open file \"" + std::string(filename) + "\".");
    }
 
    // Write the header
    vtkFile << "# vtk DataFile Version 2.0\n"
            << "Octree boxes\n"
            << "ASCII\n"
 
            << "DATASET UNSTRUCTURED_GRID\n"
            << "\n";
 
    // Write points
    vtkFile << "POINTS " << pointCount << " float\n";  // Points header
    for (Box box : boxes) {
      Position min = box.first;
      Position max = box.second;
 
      auto [minX, minY, minZ] = min;
      auto [maxX, maxY, maxZ] = max;
 
      // Write the points in the order of the VTK_HEXAHEDRON. Each point is
      // written on its own line.
 
      vtkFile << minX << " " << minY << " " << minZ << "\n";  // 0 ---
      vtkFile << maxX << " " << minY << " " << minZ << "\n";  // 1 +--
      vtkFile << maxX << " " << maxY << " " << minZ << "\n";  // 2 ++-
      vtkFile << minX << " " << maxY << " " << minZ << "\n";  // 3 -+-
      vtkFile << minX << " " << minY << " " << maxZ << "\n";  // 4 --+
      vtkFile << maxX << " " << minY << " " << maxZ << "\n";  // 5 +-+
      vtkFile << maxX << " " << maxY << " " << maxZ << "\n";  // 6 +++
      vtkFile << minX << " " << maxY << " " << maxZ << "\n";  // 7 -++
    }
    vtkFile << "\n";
 
    // Write cells
    auto cellListSize = pointCount + boxCount;
    vtkFile << "CELLS " << boxCount << " " << cellListSize << "\n";
    for (auto boxIndex = 0; boxIndex < boxCount; ++boxIndex) {
      vtkFile << "8 ";  // Output # of elements in the following line.
      for (auto pointIndex = 0; pointIndex < 8; ++pointIndex) {
        // Generate an index that points to the corresponding point in the points list.
        auto offset = 8 * boxIndex + pointIndex;
        vtkFile << offset;
        if (pointIndex < 7) {
          vtkFile << " ";
        }
      }
      vtkFile << "\n";
    }
    vtkFile << "\n";
 
    // Write cell types
    vtkFile << "CELL_TYPES " << boxCount << "\n";
    for (int i = 0; i < boxes.size(); ++i) {
      vtkFile << "12\n";  // Write VTK_HEXAHEDRON type for each cell
    }
 
    // Cleanup
    vtkFile.close();
 
    // TODO(johannes): Enclose with macro
    // #ifdef AUTOPAS_LOG_OCTREE
    // #endif
  }
 
  static FILE *leavesToJSON(FILE *out, std::vector<OctreeLeafNode<Particle_T> *> &leaves) {
    if (out) {
      fprintf(out, "[\n");
      for (int leafIndex = 0; leafIndex < leaves.size(); ++leafIndex) {
        auto leaf = leaves[leafIndex];
        fprintf(out, "{\"minmax\": ");
        outBoxCoordinatesJSON(out, leaf);
 
        // Print face neighbors
        fprintf(out, ", \"fn\": [");
        bool first = true;
        for (auto face : octree::Tables::faces) {
          auto neighbor = leaf->GTEQ_FACE_NEIGHBOR(face);
          if (neighbor) {
            if (!first) {
              fprintf(out, ", ");
            }
            first = false;
            outBoxCoordinatesJSON(out, neighbor);
          }
        }
 
        // Print face neighbor leaves
        fprintf(out, "], \"fnl\": [");
        first = true;
        for (auto face : octree::Tables::faces) {
          auto neighbor = leaf->GTEQ_FACE_NEIGHBOR(face);
          if (neighbor) {
            auto neighborLeaves = neighbor->getNeighborLeaves(face);
            for (auto neighborLeaf : neighborLeaves) {
              if (!first) {
                fprintf(out, ", ");
              }
              first = false;
              outBoxCoordinatesJSON(out, neighborLeaf);
            }
          }
        }
 
        // Print edge neighbors
        fprintf(out, "], \"en\": [");
        first = true;
        for (auto edge : octree::Tables::edges) {
          auto neighbor = leaf->GTEQ_EDGE_NEIGHBOR(edge);
          if (neighbor) {
            if (!first) {
              fprintf(out, ", ");
            }
            first = false;
            outBoxCoordinatesJSON(out, neighbor);
          }
        }
 
        // Print edge neighbor leaves
        fprintf(out, "], \"enl\": [");
        first = true;
        for (auto edge : octree::Tables::edges) {
          auto neighbor = leaf->GTEQ_EDGE_NEIGHBOR(edge);
          if (neighbor) {
            auto neighborLeaves = neighbor->getNeighborLeaves(edge);
            for (auto neighborLeaf : neighborLeaves) {
              if (!first) {
                fprintf(out, ", ");
              }
              first = false;
              outBoxCoordinatesJSON(out, neighborLeaf);
            }
          }
        }
 
        // Print vertex neighbors
        fprintf(out, "], \"vn\": [");
        first = true;
        for (auto vertex : octree::Tables::vertices) {
          auto neighbor = leaf->GTEQ_VERTEX_NEIGHBOR(vertex);
          if (neighbor) {
            if (!first) {
              fprintf(out, ", ");
            }
            first = false;
            outBoxCoordinatesJSON(out, neighbor);
          }
        }
 
        // Print vertex neighbor leaves
        fprintf(out, "], \"vnl\": [");
        first = true;
        for (auto vertex : octree::Tables::vertices) {
          auto neighbor = leaf->GTEQ_VERTEX_NEIGHBOR(vertex);
          if (neighbor) {
            auto neighborLeaves = neighbor->getNeighborLeaves(vertex);
 
            for (auto neighborLeaf : neighborLeaves) {
              if (!first) {
                fprintf(out, ", ");
              }
              first = false;
              outBoxCoordinatesJSON(out, neighborLeaf);
            }
          }
        }
 
        fprintf(out, "]}");
        if (leafIndex < (leaves.size() - 1)) {
          fprintf(out, ",");
        }
        fprintf(out, "\n");
      }
      fprintf(out, "]\n");
    } else {
      fprintf(stderr, "ERROR: Dump file is nullptr.\n");
    }
 
    return out;
  }
 
  static FILE *particlesToJSON(FILE *out, OctreeNodeInterface<Particle_T> *root) {
    if (out) {
      // Get all leaves
      std::vector<OctreeLeafNode<Particle_T> *> leaves;
      root->appendAllLeaves(leaves);
 
      fprintf(out, "[");
      for (int leafIndex = 0; leafIndex < leaves.size(); ++leafIndex) {
        OctreeLeafNode<Particle_T> *leaf = leaves[leafIndex];
 
        auto n = leaf->size();
        for (int particleIndex = 0; particleIndex < n; ++particleIndex) {
          Particle_T &particle = leaf->at(particleIndex);
          auto p = particle.getR();
          fputc('[', out);
          out3(out, p[0], p[1], p[2]);
          fputc(']', out);
          if ((particleIndex < (n - 1)) || (leafIndex < (leaves.size() - 1))) {
            fprintf(out, ",");
          }
        }
      }
      fprintf(out, "]");
    } else {
      fprintf(stderr, "ERROR: Dump file is nullptr.\n");
    }
 
    return out;
  }
 
  static void octreeToJSON(OctreeNodeInterface<Particle_T> *owned, OctreeNodeInterface<Particle_T> *halo,
                           std::vector<OctreeLeafNode<Particle_T> *> &ownedLeaves,
                           std::vector<OctreeLeafNode<Particle_T> *> &haloLeaves) {
#if AUTOPAS_LOG_OCTREE
    // Log all owned leaves for this octree
    fclose(OctreeLogger::leavesToJSON(fopen("owned.json", "w"), ownedLeaves));
    // Log all halo leaves for this octree
    fclose(OctreeLogger::leavesToJSON(fopen("halo.json", "w"), haloLeaves));
    FILE *particles = fopen("particles.json", "w");
    fprintf(particles, "{\"owned\": ");
    OctreeLogger::particlesToJSON(particles, owned);
    fprintf(particles, ", \"halo\": \n");
    OctreeLogger::particlesToJSON(particles, halo);
    fprintf(particles, "}");
    fclose(particles);
#endif
  }
 
 private:
  static void out3(FILE *out, double x, double y, double z) { fprintf(out, "%.3f,%.3f,%.3f", x, y, z); }
 
  static void outLocationArrayJSON(FILE *out, const std::array<double, 3> &min, const std::array<double, 3> &max) {
    fputc('[', out);
    out3(out, min[0], min[1], min[2]);
    fputc(',', out);
    out3(out, max[0], max[1], max[2]);
    fputc(']', out);
  }
 
  static void outBoxCoordinatesJSON(FILE *out, OctreeNodeInterface<Particle_T> *node) {
    const auto min = node->getBoxMin();
    const auto max = node->getBoxMax();
    outLocationArrayJSON(out, min, max);
  }
 
  int unsigned iteration = 0;
};
}  // namespace autopas