DLR-AMR · jmark · Jul 17, 2023 · Jul 3, 2023 · Jul 3, 2023 · Jul 3, 2023
diff --git a/configure.ac b/configure.ac
@@ -82,6 +82,9 @@ AC_CONFIG_LINKS([test/testfiles/test_msh_file_vers4_ascii.msh:test/testfiles/tes
 AC_CONFIG_LINKS([test/testfiles/test_msh_file_vers4_bin.msh:test/testfiles/test_msh_file_vers4_bin.msh])
 AC_CONFIG_LINKS([test/testfiles/test_vtk_tri.vtu:test/testfiles/test_vtk_tri.vtu])
 AC_CONFIG_LINKS([test/testfiles/test_vtk_cube.vtp:test/testfiles/test_vtk_cube.vtp])
+AC_CONFIG_LINKS([test/testfiles/test_parallel_file.pvtu:test/testfiles/test_parallel_file.pvtu])
+AC_CONFIG_LINKS([test/testfiles/test_parallel_file_0.vtu:test/testfiles/test_parallel_file_0.vtu])
+AC_CONFIG_LINKS([test/testfiles/test_parallel_file_1.vtu:test/testfiles/test_parallel_file_1.vtu])
 AC_CONFIG_LINKS([example/IO/cmesh/gmsh/circlesquare_hybrid_hole.msh:example/IO/cmesh/gmsh/circlesquare_hybrid_hole.msh])
 AC_CONFIG_FILES([tutorials/features/t8_features_curved_meshes_generate_cmesh.geo:tutorials/features/t8_features_curved_meshes_generate_cmesh.geo])
 AC_CONFIG_FILES([tutorials/features/t8_features_curved_meshes_generate_cmesh_2d.geo:tutorials/features/t8_features_curved_meshes_generate_cmesh_2d.geo])

diff --git a/example/IO/cmesh/vtk/t8_cmesh_read_from_vtk.cxx b/example/IO/cmesh/vtk/t8_cmesh_read_from_vtk.cxx
@@ -162,11 +162,13 @@ main (int argc, char **argv)
                          "The prefix of the output-file.");
   sc_options_add_int (opt, 'c', "num_cell_values", &num_keys, 0,
                       "Number of values per cell stored in the vtk-file.");
-  sc_options_add_int (opt, 't', "type_of_file", &vtk_file_type, -1,
-                      " Set the type of the data in the file.\n"
-                      "0 for vtkUnstructuredGrid \n" "1 for vtkPolyData");
   sc_options_add_bool (opt, 'p', "partition", &partition, 0,
-                       "If set, partition the cmesh uniformly.\n");
+                       "If set, partition the cmesh uniformly.");
+  sc_options_add_int (opt, 't', "type_of_file", &vtk_file_type, -1,
+                      "Set the type of the data in the file.\n"
+                      "\t\t\t\t\t0 for vtkUnstructuredGrid \n"
+                      "\t\t\t\t\t1 for vtkPolyData\n"
+                      "\t\t\t\t\t2 for pvtu. Currently not working with -p");
   parsed =
     sc_options_parse (t8_get_package_id (), SC_LP_ERROR, opt, argc, argv);
 

diff --git a/src/Makefile.am b/src/Makefile.am
@@ -87,6 +87,7 @@ libt8_internal_headers = \
   src/t8_cmesh/t8_cmesh_offset.h \
   src/t8_vtk/t8_vtk_polydata.hxx \
   src/t8_vtk/t8_vtk_unstructured.hxx \
+  src/t8-vtk/t8_vtk_parallel.hxx \
   src/t8_forest/t8_forest_cxx.h  \
   src/t8_forest/t8_forest_ghost.h \
   src/t8_forest/t8_forest_balance.h src/t8_forest/t8_forest_types.h \
@@ -129,6 +130,7 @@ libt8_compiled_sources = \
   src/t8_cmesh/t8_cmesh_testcases.c \
   src/t8_vtk/t8_vtk_polydata.cxx \
   src/t8_vtk/t8_vtk_unstructured.cxx \
+  src/t8_vtk/t8_vtk_parallel.cxx \
   src/t8_vtk/t8_vtk_reader.cxx
 
 

diff --git a/src/t8_vtk/t8_vtk_parallel.cxx b/src/t8_vtk/t8_vtk_parallel.cxx
@@ -0,0 +1,104 @@
+/*
+This file is part of t8code.
+t8code is a C library to manage a collection (a forest) of multiple
+connected adaptive space-trees of general element classes in parallel.
+
+Copyright (C) 2023 the developers
+
+t8code is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+t8code is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with t8code; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+*/
+
+#include "t8_vtk_parallel.hxx"
+
+#if T8_WITH_VTK
+#include <vtkUnstructuredGrid.h>
+#include <vtkXMLPUnstructuredGridReader.h>
+#include <vtkAppendFilter.h>
+
+vtk_read_success_t
+t8_read_parallel (const char *filename, vtkSmartPointer < vtkDataSet > grid,
+                  sc_MPI_Comm comm)
+{
+  /* Check if we can open the parallel file */
+  FILE               *first_check;
+  first_check = fopen (filename, "r");
+  if (first_check == NULL) {
+    t8_errorf ("Can not find the file %s\n", filename);
+    return read_failure;
+  }
+  fclose (first_check);
+  /* Setup parallel reader. */
+  vtkSmartPointer < vtkXMLPUnstructuredGridReader > reader =
+    vtkSmartPointer < vtkXMLPUnstructuredGridReader >::New ();
+  if (!reader->CanReadFile (filename)) {
+    t8_errorf ("Unable to read file.\n");
+    return read_failure;
+  }
+
+  reader->SetFileName (filename);
+  reader->UpdateInformation ();
+
+  /* Get mpi size and rank */
+  const int           total_num_pieces = reader->GetNumberOfPieces ();
+  int                 mpiret;
+  int                 mpisize;
+  mpiret = sc_MPI_Comm_size (comm, &mpisize);
+  SC_CHECK_MPI (mpiret);
+
+  int                 mpirank;
+  mpiret = sc_MPI_Comm_rank (comm, &mpirank);
+  SC_CHECK_MPI (mpiret);
+
+  /* Setup number of pieces to read on this proc. */
+  int                 last_piece = -1;
+  int                 first_piece = 0;
+  if (mpisize >= total_num_pieces) {
+    /* The first n-procs read a piece each. */
+    first_piece = mpirank;
+    last_piece = first_piece + ((mpirank < mpisize) ? 1 : 0);
+  }
+  else {
+    first_piece = total_num_pieces / mpisize * mpirank;
+    last_piece = first_piece;
+    const int           prev_proc_first_piece =
+      total_num_pieces / mpisize * (mpirank - 1);
+    last_piece +=
+      first_piece == prev_proc_first_piece ? 0 : total_num_pieces / mpisize;
+    if (first_piece == total_num_pieces / mpisize * (mpisize - 1)) {
+      /* Read the last chunk of data */
+      last_piece = total_num_pieces - first_piece;
+    }
+  }
+
+  /* Read the pieces if there are any pieces to read on this proc. */
+  if (first_piece < last_piece) {
+    vtkNew < vtkAppendFilter > append;
+    for (int ipiece = first_piece; ipiece < last_piece; ipiece++) {
+      reader->UpdatePiece (ipiece, total_num_pieces, 0);
+      append->AddInputData (reader->GetOutput ());
+    }
+    /* Merge all read grids together */
+    append->Update ();
+    append->MergePointsOn ();
+    grid->ShallowCopy (append->GetOutput ());
+    t8_debugf ("[D] read %lli cells\n", grid->GetNumberOfCells ());
+  }
+  else {
+    t8_debugf ("[D] dont read any file on this proc\n");
+  }
+  return read_success;
+}
+
+#endif
diff --git a/src/t8_vtk/t8_vtk_parallel.hxx b/src/t8_vtk/t8_vtk_parallel.hxx
@@ -0,0 +1,47 @@
+/*
+This file is part of t8code.
+t8code is a C library to manage a collection (a forest) of multiple
+connected adaptive space-trees of general element classes in parallel.
+
+Copyright (C) 2023 the developers
+
+t8code is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+t8code is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with t8code; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+*/
+
+#ifndef T8_VTK_PARALLEL_HXX
+#define T8_VTK_PARALLEL_HXX
+
+#include <t8.h>
+#include "t8_vtk_types.h"
+
+#if T8_WITH_VTK
+#include <vtkDataSet.h>
+#include <vtkSmartPointer.h>
+
+/**
+ * Given a filename to a parallel vtk file (for example .pvtu) and its data files, 
+ * read a piece of the data files (like .vtu, .vtp, ...).
+ * 
+ * \param[in] filename  The name of a parallel vtk file (.pvtu for example)
+ * \param[in] grid On input a vtkSmartPointer, that will hold the grid described
+ *                 by the pieces read on this proc. 
+ * \returns        non-zero on success, zero if the reading failed. 
+ */
+vtk_read_success_t  t8_read_parallel (const char *filename,
+                                      vtkSmartPointer < vtkDataSet > grid,
+                                      sc_MPI_Comm comm);
+
+#endif /* T8_WITH_VTK */
+#endif /* T8_VTK_PARALLEL_HXX */
diff --git a/src/t8_vtk/t8_vtk_reader.cxx b/src/t8_vtk/t8_vtk_reader.cxx
@@ -3,7 +3,7 @@ This file is part of t8code.
 t8code is a C library to manage a collection (a forest) of multiple
 connected adaptive space-trees of general element classes in parallel.
 
-Copyright (C) 2015 the developers
+Copyright (C) 2023 the developers
 
 t8code is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -25,6 +25,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc.,
 #include "t8_vtk_reader.hxx"
 #include "t8_vtk_unstructured.hxx"
 #include "t8_vtk_polydata.hxx"
+#include "t8_vtk_parallel.hxx"
 #include "t8_vtk_types.h"
 #include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.h>
 
@@ -135,6 +136,16 @@ t8_file_to_vtkGrid (const char *filename,
     case VTK_POLYDATA_FILE:
       main_proc_read_successful = t8_read_poly (filename, vtkGrid);
       break;
+    case VTK_PARALLEL_FILE:
+      if (!partition) {
+        main_proc_read_successful = t8_read_unstructured (filename, vtkGrid);
+      }
+      else {
+        t8_errorf ("Distributed reading is not supported yet.\n");
+        vtkGrid = NULL;
+        break;
+      }
+      break;
     default:
       vtkGrid = NULL;
       t8_errorf ("Filetype not supported.\n");
@@ -402,7 +413,7 @@ t8_vtk_reader (const char *filename, const int partition,
   SC_CHECK_MPI (mpiret);
   mpiret = sc_MPI_Comm_rank (comm, &mpirank);
   SC_CHECK_MPI (mpiret);
-
+  t8_debugf ("[D] file_type: %i\n", (int) vtk_file_type);
   /* Ensure that the main-proc is a valid proc. */
   T8_ASSERT (0 <= main_proc && main_proc < mpisize);
   T8_ASSERT (filename != NULL);
@@ -416,6 +427,9 @@ t8_vtk_reader (const char *filename, const int partition,
   case VTK_POLYDATA_FILE:
     vtkGrid = vtkSmartPointer < vtkPolyData >::New ();
     break;
+  case VTK_PARALLEL_FILE:
+    vtkGrid = vtkSmartPointer < vtkUnstructuredGrid >::New ();
+    break;
   default:
     t8_errorf ("Filetype is not supported.\n");
     break;

diff --git a/src/t8_vtk/t8_vtk_reader.hxx b/src/t8_vtk/t8_vtk_reader.hxx
@@ -3,7 +3,7 @@ This file is part of t8code.
 t8code is a C library to manage a collection (a forest) of multiple
 connected adaptive space-trees of general element classes in parallel.
 
-Copyright (C) 2015 the developers
+Copyright (C) 2023 the developers
 
 t8code is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by

diff --git a/src/t8_vtk/t8_vtk_types.h b/src/t8_vtk/t8_vtk_types.h
@@ -60,9 +60,10 @@ const t8_eclass_t   t8_cmesh_vtk_type_to_t8_type[82] = {
  */
 typedef enum vtk_file_type
 {
-  VTK_FILE_ERROR = -1,          /*For Testing purpose. */
+  VTK_FILE_ERROR = -1,          /* For Testing purpose. */
   VTK_UNSTRUCTURED_FILE = 0,
   VTK_POLYDATA_FILE = 1,
+  VTK_PARALLEL_FILE = 2,        /* All parallel vtk files. */
   VTK_NUM_TYPES
 } vtk_file_type_t;
 

diff --git a/src/t8_vtk/t8_vtk_unstructured.cxx b/src/t8_vtk/t8_vtk_unstructured.cxx
@@ -26,6 +26,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc.,
 #include <vtkUnstructuredGrid.h>
 #include <vtkUnstructuredGridReader.h>
 #include <vtkXMLUnstructuredGridReader.h>
+#include <vtkXMLPUnstructuredGridReader.h>
 
 vtk_read_success_t
 t8_read_unstructured (const char *filename,
@@ -70,6 +71,18 @@ t8_read_unstructured (const char *filename,
     grid->ShallowCopy (vtkDataSet::SafeDownCast (reader->GetOutput ()));
     return read_success;
   }
+  else if (strcmp (extension, "pvtu") == 0) {
+    vtkSmartPointer < vtkXMLPUnstructuredGridReader > reader =
+      vtkSmartPointer < vtkXMLPUnstructuredGridReader >::New ();
+    if (!reader->CanReadFile (filename)) {
+      t8_errorf ("Unable to read file.\n");
+      return read_failure;
+    }
+    reader->SetFileName (filename);
+    reader->Update ();
+    grid->ShallowCopy (vtkDataSet::SafeDownCast (reader->GetOutput ()));
+    return read_success;
+  }
   else {
     /* Return failure if the reader is not used correctly */
     t8_global_errorf ("Please use .vtk or .vtu file\n");

diff --git a/test/t8_IO/t8_gtest_vtk_reader.cxx b/test/t8_IO/t8_gtest_vtk_reader.cxx
@@ -46,23 +46,28 @@ class vtk_reader : public testing::TestWithParam<std::tuple<vtk_file_type_t, int
       file = (int)file_type + 1;
       partition = std::get<1>(GetParam());
       main_proc = std::get<2>(GetParam());
+      if(file_type == VTK_PARALLEL_FILE && partition){
+        GTEST_SKIP();
+      }
     }
     int file;
     vtk_file_type_t file_type;
     int partition;
     int main_proc;
-    const char* failing_files[3] = {
+    const char* failing_files[4] = {
       "no_file",
       "non-existing-file.vtu",
-      "non-existing-file.vtp"
+      "non-existing-file.vtp",
+      "non-existing-file.pvtu"
     };
-    const char* test_files[3] = {
+    const char* test_files[4] = {
       "no_file",
       "test/testfiles/test_vtk_tri.vtu",
-      "test/testfiles/test_vtk_cube.vtp"
+      "test/testfiles/test_vtk_cube.vtp",
+      "test/testfiles/test_parallel_file.pvtu"
     };
-    const int num_points[3] = {0, 121, 24};
-    const int num_trees[3] = {0, 200, 12};
+    const int num_points[4] = {0, 121, 24, 6144};
+    const int num_trees[4] = {0, 200, 12, 1024};
 };
 /* *INDENT-ON* */
 
@@ -119,6 +124,8 @@ TEST_P (vtk_reader, vtk_to_pointSet)
 }
 
 /* *INDENT-OFF* */
+/* Currently does not work for parallel files. Replace with VTK_NUM_TYPES as soon
+ * as reading and constructing cmeshes from parallel files is enabled. */
 INSTANTIATE_TEST_SUITE_P (t8_gtest_vtk_reader, vtk_reader,testing::Combine (
                           testing::Range (VTK_FILE_ERROR, VTK_NUM_TYPES),
                           testing::Values (0, 1),

diff --git a/test/testfiles/test_parallel_file.pvtu b/test/testfiles/test_parallel_file.pvtu
@@ -0,0 +1,15 @@
+<VTKFile type="PUnstructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <PUnstructuredGrid GhostLevel="0">
+    <PCellData>
+      <PDataArray type="Int64" Name="treeid"/>
+      <PDataArray type="Int64" Name="mpirank"/>
+      <PDataArray type="Int64" Name="level"/>
+      <PDataArray type="Int64" Name="element_id"/>
+    </PCellData>
+    <PPoints>
+      <PDataArray type="Float32" Name="Points" NumberOfComponents="3"/>
+    </PPoints>
+    <Piece Source="test_parallel_file_0.vtu"/>
+    <Piece Source="test_parallel_file_1.vtu"/>
+  </PUnstructuredGrid>
+</VTKFile>
diff --git a/test/testfiles/test_parallel_file_0.vtu b/test/testfiles/test_parallel_file_0.vtu
diff --git a/test/testfiles/test_parallel_file_1.vtu b/test/testfiles/test_parallel_file_1.vtu