Added script splitvtufiles.py into the postprocess_scripts directory

postprocess_scripts/README.md

@@ -67,3 +67,16 @@ This script calculates the following quantities for each of the time states:
 Results for all time states are outputted in the file **domain_stats_vs_t.txt**.
 
 For each row, the first column corresponds to the frame number, the second to the time,  the third column to the number of domains, the fourth column to the average size, and the fifth column to the standard deviation of the size.
+
+#### splitvtufiles.py
+
+This script splits the output .vtu files into several files, allowing for parallel visualization in VisIt or ParaView. To run it simply 
+1. Open a terminal
+
+2. Go to the directory where the series of output files from a single simulation are located (these files must be generated using the single file per process option, see the **Output** section from [this page](https://prisms-center.github.io/phaseField/doxygen_files/input_file.html))
+
+3. Type:
+
+  \$ python splitvtufiles.py <no. files>
+
+  where <no. files> is the number of files into which you want to split each solution\*.vtu file.

postprocess_scripts/splitvtufiles.py

@@ -0,0 +1,77 @@
+import pyvista as pv
+import numpy as np
+import os
+import glob
+import sys
+
+def split_mesh(input_file, output_prefix, number_of_splits):
+    # Read the .vtu file
+    mesh = pv.read(input_file)
+
+    # Number of cells in each split
+    cells_per_split = mesh.n_cells // number_of_splits
+
+    split_filenames = []
+
+    for i in range(number_of_splits):
+        start = i * cells_per_split
+        end = start + cells_per_split
+        if i == number_of_splits - 1:
+            end = mesh.n_cells
+
+        # Extract a subset of cells
+        submesh = mesh.extract_cells(np.arange(start, end))
+
+        # Write the submesh to a new .vtu file
+        split_filename = f"{output_prefix}.{i}.vtu"
+        submesh.save(split_filename)
+        split_filenames.append(split_filename)
+
+    # Create the .pvtu file
+    create_pvtu_file(output_prefix, split_filenames, mesh)
+
+def create_pvtu_file(base_name, split_filenames, reference_mesh):
+    with open(f"{base_name}.pvtu", 'w') as file:
+        file.write('<?xml version="1.0"?>\n')
+        file.write('<VTKFile type="PUnstructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32">\n')
+        file.write('  <PUnstructuredGrid GhostLevel="0">\n')
+
+        # Write point data arrays
+        for array in reference_mesh.point_data.keys():
+            num_components = reference_mesh.point_data[array].shape[1] if len(reference_mesh.point_data[array].shape) > 1 else 1
+            data_type = pv.get_array_association(reference_mesh, array)
+            file.write(f'    <PDataArray type="{data_type}" Name="{array}" NumberOfComponents="{num_components}"/>\n')
+
+        # Write piece references
+        for filename in split_filenames:
+            file.write(f'    <Piece Source="{filename}"/>\n')
+
+        file.write('  </PUnstructuredGrid>\n')
+        file.write('</VTKFile>\n')
+
+def process_vtu_files(no_parts):
+    # Construct the search pattern for .vtu files
+    search_pattern = os.path.join('solution-*.vtu')
+
+    # Find all .vtu files in the specified directory
+    vtu_files = glob.glob(search_pattern)
+
+    # Process each .vtu file sequentially
+    for file in vtu_files:
+        print(f"Splitting {file}")
+        # Add your code here to process each file
+        # For example, reading the file, performing some analysis, etc.
+        filenm_noext=file[:-4]
+        split_mesh(file, filenm_noext, int(no_parts))
+
+# Main
+# Call the function
+try:
+    process_vtu_files(sys.argv[1])
+except Exception as e:
+    print(f"An error occurred: {e}")
+    print("Error: A single integer argument must be provided specifying the number files in which the system will be split")
+
+
+
+