From a81c69a74b2ad41e82a41a860a0ce9a0aa36f51c Mon Sep 17 00:00:00 2001 From: aplowman Date: Thu, 28 Mar 2024 20:24:52 +0000 Subject: [PATCH] fix: use pyvista to resample VTU files onto a grid --- cipher_parse/cipher_output.py | 154 +++++----------------------------- 1 file changed, 20 insertions(+), 134 deletions(-) diff --git a/cipher_parse/cipher_output.py b/cipher_parse/cipher_output.py index 9e8a479..027b303 100644 --- a/cipher_parse/cipher_output.py +++ b/cipher_parse/cipher_output.py @@ -22,7 +22,6 @@ ) from cipher_parse.derived_outputs import num_voxels_per_phase -DEFAULT_PARAVIEW_EXE = "pvbatch" INC_DATA_NON_ARRAYS = ( "increment", "time", @@ -109,53 +108,6 @@ def parse_cipher_stdout(path_or_string, is_string=False): return out -def generate_VTI_files_from_VTU_files( - sampling_dimensions, - paraview_exe=DEFAULT_PARAVIEW_EXE, -): - """Generate a 'ParaView-python' script for generating VTI files from VTU files and - execute that script.""" - - script_name = "vtu2vti.py" - if len(sampling_dimensions) == 2: - sampling_dimensions += [1] - - with Path(script_name).open("wt") as fp: - fp.write( - dedent( - f""" - import os - - from paraview.simple import * - - vtu_files = [] - for root, dirs, files in os.walk("."): - for f in files: - if f.endswith(".vtu"): - vtu_files.append(f) - - for file_i_path in vtu_files: - file_i_base_name = file_i_path.split(".")[0] - vtu_data_i = XMLUnstructuredGridReader( - FileName=[os.getcwd() + os.path.sep + file_i_path] - ) - resampleToImage1 = ResampleToImage(Input=vtu_data_i) - resampleToImage1.SamplingDimensions = {sampling_dimensions!r} - SetActiveSource(resampleToImage1) - SaveData(file_i_base_name + ".vti", resampleToImage1) - """ - ) - ) - - proc = run(f"{paraview_exe} {script_name}", shell=True, stdout=PIPE, stderr=PIPE) - stdout = proc.stdout.decode() - stderr = proc.stderr.decode() - if stdout: - print(stdout) - if stderr: - print(stderr) - - class CIPHEROutput: """Class to hold output information from a CIPHER simulation.""" @@ -175,12 +127,6 @@ def __init__( cipher_input=None, ): default_options = { - "paraview_exe": DEFAULT_PARAVIEW_EXE, - "delete_VTIs": True, - "delete_VTUs": False, - "use_existing_VTIs": False, - "num_VTU_files": None, - "VTU_files_time_interval": None, "derive_outputs": None, "save_outputs": None, } @@ -201,14 +147,6 @@ def __init__( self._cipher_stdout = None self._geometries = None # assigned by set_geometries - if ( - options.get("VTU_files_time_interval") is not None - and options.get("num_VTU_files") is not None - ): - raise ValueError( - "Specify at most one of 'num_VTU_files' and 'VTU_files_time_interval'." - ) - for idx, i in enumerate(options["save_outputs"]): if i.get("number") is not None and i.get("time_interval") is not None: raise ValueError( @@ -282,12 +220,6 @@ def get_incremental_data(self): inp_dat = self.get_input_YAML_data() grid_size = inp_dat["grid_size"] - if not self.options["use_existing_VTIs"]: - generate_VTI_files_from_VTU_files( - sampling_dimensions=grid_size, - paraview_exe=self.options["paraview_exe"], - ) - outfile_base = inp_dat["solution_parameters"]["outfile"] output_lookup = { i: f"{outfile_base} output.{idx}" @@ -297,85 +229,46 @@ def get_incremental_data(self): list(self.directory.glob(f"{outfile_base}_*.vtu")), key=lambda x: int(re.search(r"\d+", x.name).group()), ) - vti_file_list = sorted( - list(self.directory.glob(f"{outfile_base}_*.vti")), - key=lambda x: int(re.search(r"\d+", x.name).group()), - ) - - # Move all VTU files to a sub-directory: - viz_dir = self.directory / "original_viz" - - vtu_orig_file_list = [] - if not viz_dir.is_dir(): - viz_dir.mkdir() - for viz_file_i in vtu_file_list: - dst_i = viz_dir.joinpath(viz_file_i.name).with_suffix( - ".viz" + viz_file_i.suffix - ) - shutil.move(viz_file_i, dst_i) - vtu_orig_file_list.append(dst_i) - else: - vtu_orig_file_list = sorted( - list(viz_dir.glob("*")), - key=lambda x: int(re.search(r"\d+", x.name).group()), - ) - - # Copy back to the root directory VTU files that we want to keep: - if self.options["num_VTU_files"]: - viz_files_keep_idx = get_subset_indices( - len(vti_file_list), - self.options["num_VTU_files"], - ) - elif self.options["VTU_files_time_interval"]: - viz_files_keep_idx = self._get_time_linear_subset_indices( - time_interval=self.options["VTU_files_time_interval"] - ) - else: - viz_files_keep_idx = [] - - for i in viz_files_keep_idx: - viz_file_i = vtu_orig_file_list[i] - dst_i = ( - self.directory.joinpath(viz_file_i.name) - .with_suffix("") - .with_suffix(".vtu") - ) - shutil.copy(viz_file_i, dst_i) - - if self.options["delete_VTUs"]: - print(f"Deleting original VTU files in directory: {viz_dir}") - shutil.rmtree(viz_dir) # get which files to include for each output/derived output outputs_keep_idx = {} for save_out_i in self.options["save_outputs"]: if "number" in save_out_i: - keep_idx = get_subset_indices(len(vti_file_list), save_out_i["number"]) + keep_idx = get_subset_indices(len(vtu_file_list), save_out_i["number"]) elif "time_interval" in save_out_i: keep_idx = self._get_time_linear_subset_indices( time_interval=save_out_i["time_interval"] ) else: - keep_idx = list(range(len(vti_file_list))) + keep_idx = list(range(len(vtu_file_list))) outputs_keep_idx[save_out_i["name"]] = keep_idx incremental_data = [] - for file_i_idx, file_i in enumerate(vti_file_list): + for file_i_idx, file_i in enumerate(vtu_file_list): + print(f"Reading VTU file {file_i.name}...", flush=True) mesh = pv.get_reader(file_i).read() - vtu_file_name = file_i.name.replace("vti", "vtu") + vtu_file_name = file_i.name + + img_data = pv.ImageData(dimensions=grid_size) + print( + f"Resampling VTU file {file_i.name} onto an image-data mesh...", + flush=True, + ) + img_mesh = img_data.sample(mesh) + inc_data_i = { - "increment": int(re.search(r"\d+", file_i.name).group()), + "increment": int(re.search(r"\d+", vtu_file_name).group()), "time": self.cipher_stdout["outputs"][vtu_file_name], - "dimensions": list(mesh.dimensions), - "spacing": list(mesh.spacing), - "number_VTI_cells": mesh.number_of_cells, - "number_VTI_points": mesh.number_of_points, + "dimensions": list(img_mesh.dimensions), + "spacing": list(img_mesh.spacing), + "number_VTI_cells": img_mesh.number_of_cells, + "number_VTI_points": img_mesh.number_of_points, } standard_outputs = {} for name in output_lookup: - arr_flat = mesh.get_array(output_lookup[name]) - arr = arr_flat.reshape(mesh.dimensions, order="F") + arr_flat = img_mesh.get_array(output_lookup[name]) + arr = arr_flat.reshape(img_mesh.dimensions, order="F") if name in STANDARD_OUTPUTS_TYPES: arr = arr.astype(STANDARD_OUTPUTS_TYPES[name]) standard_outputs[name] = np.array(arr) # convert from pyvista_ndarray @@ -401,13 +294,6 @@ def get_incremental_data(self): incremental_data.append(inc_data_i) - if self.options["delete_VTIs"] and not self.options["use_existing_VTIs"]: - for file_i in vti_file_list: - print(f"Deleting temporary VTI file: {file_i}") - os.remove(file_i) - - outputs_keep_idx["VTU_files"] = viz_files_keep_idx - return incremental_data, outputs_keep_idx @property