From 8deee3d954fc2514f28f1c4de7861dbf5a799bf8 Mon Sep 17 00:00:00 2001
From: Rafael Ravedutti <rafaelravedutti@gmail.com>
Date: Tue, 8 Feb 2022 00:11:10 +0100
Subject: [PATCH] Add cluster edges in VTK output

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
---
 gromacs/includes/vtk.h |  1 +
 gromacs/main.c         |  2 ++
 gromacs/vtk.c          | 41 +++++++++++++++++++++++++++++++++++++++++
 3 files changed, 44 insertions(+)

diff --git a/gromacs/includes/vtk.h b/gromacs/includes/vtk.h
index 66b3b32..f839835 100644
--- a/gromacs/includes/vtk.h
+++ b/gromacs/includes/vtk.h
@@ -26,4 +26,5 @@
 #define __VTK_H_
 extern int write_local_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep);
 extern int write_ghost_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep);
+extern int write_cluster_edges_to_vtk_file(const char* filename, Atom* atom, int timestep);
 #endif
diff --git a/gromacs/main.c b/gromacs/main.c
index 4694741..6b2c5e1 100644
--- a/gromacs/main.c
+++ b/gromacs/main.c
@@ -273,6 +273,7 @@ int main(int argc, char** argv) {
     if(param.vtk_file != NULL) {
         write_local_atoms_to_vtk_file(param.vtk_file, &atom, 0);
         write_ghost_atoms_to_vtk_file(param.vtk_file, &atom, 0);
+        write_cluster_edges_to_vtk_file(param.vtk_file, &atom, 0);
     }
 
     for(int n = 0; n < param.ntimes; n++) {
@@ -303,6 +304,7 @@ int main(int argc, char** argv) {
         if(param.vtk_file != NULL) {
             write_local_atoms_to_vtk_file(param.vtk_file, &atom, n + 1);
             write_ghost_atoms_to_vtk_file(param.vtk_file, &atom, n + 1);
+            write_cluster_edges_to_vtk_file(param.vtk_file, &atom, n + 1);
         }
     }
 
diff --git a/gromacs/vtk.c b/gromacs/vtk.c
index 4b788ea..a9b0b34 100644
--- a/gromacs/vtk.c
+++ b/gromacs/vtk.c
@@ -88,3 +88,44 @@ int write_ghost_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep
     fclose(fp);
     return 0;
 }
+
+int write_cluster_edges_to_vtk_file(const char* filename, Atom* atom, int timestep) {
+    char timestep_filename[128];
+    snprintf(timestep_filename, sizeof timestep_filename, "%s_edges_%d.vtk", filename, timestep);
+    FILE* fp = fopen(timestep_filename, "wb");
+    int Nclusters_all = atom->Nclusters_local + atom->Nclusters_ghost;
+    int tot_lines = 0;
+    int i = 0;
+
+    if(fp == NULL) {
+        fprintf(stderr, "Could not open VTK file for writing!\n");
+        return -1;
+    }
+
+    fprintf(fp, "# vtk DataFile Version 2.0\n");
+    fprintf(fp, "Particle data\n");
+    fprintf(fp, "ASCII\n");
+    fprintf(fp, "DATASET POLYDATA\n");
+    fprintf(fp, "POINTS %d double\n", atom->Nlocal + atom->Nghost);
+    for(int ci = 0; ci < Nclusters_all; ++ci) {
+        MD_FLOAT *cptr = cluster_pos_ptr(ci);
+        for(int cii = 0; cii < atom->clusters[ci].natoms; ++cii) {
+            fprintf(fp, "%.4f %.4f %.4f\n", cluster_x(cptr, cii), cluster_y(cptr, cii), cluster_z(cptr, cii));
+        }
+
+        tot_lines += atom->clusters[ci].natoms;
+    }
+    fprintf(fp, "\n\n");
+    fprintf(fp, "LINES %d %d\n", Nclusters_all, Nclusters_all + tot_lines);
+    for(int ci = 0; ci < Nclusters_all; ++ci) {
+        fprintf(fp, "%d ", atom->clusters[ci].natoms);
+        for(int cii = 0; cii < atom->clusters[ci].natoms; ++cii) {
+            fprintf(fp, "%d ", i++);
+        }
+
+        fprintf(fp, "\n");
+    }
+    fprintf(fp, "\n\n");
+    fclose(fp);
+    return 0;
+}