diff --git a/src/vasp/automatedPreprocessing/automated_preprocessing.py b/src/vasp/automatedPreprocessing/automated_preprocessing.py index 85a39ff..7057081 100755 --- a/src/vasp/automatedPreprocessing/automated_preprocessing.py +++ b/src/vasp/automatedPreprocessing/automated_preprocessing.py @@ -24,7 +24,8 @@ from vampy.simulation.simulation_common import print_mesh_information from vasp.automatedPreprocessing.preprocessing_common import generate_mesh, distance_to_spheres_solid_thickness, \ - dist_sphere_spheres, convert_xml_mesh_to_hdf5, convert_vtu_mesh_to_xdmf, edge_length_evaluator + dist_sphere_spheres, convert_xml_mesh_to_hdf5, convert_vtu_mesh_to_xdmf, edge_length_evaluator, \ + check_flatten_boundary from vasp.simulations.simulation_common import load_mesh_and_data @@ -178,6 +179,7 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f # Get centerlines print("--- Get centerlines\n") inlet, outlets = get_centers_for_meshing(surface, has_multiple_inlets, str(base_path)) + num_inlets_outlets = len(inlet) // 3 + len(outlets) // 3 # Get point the furthest away from the inlet when only one boundary has_outlet = len(outlets) != 0 @@ -510,6 +512,10 @@ def try_generate_mesh(distance_to_sphere, number_of_sublayers_fluid, number_of_s print("--- Converting XML mesh to HDF5\n") convert_xml_mesh_to_hdf5(file_name_xml_mesh, scale_factor_h5) + # NOTE: stdev threshold is currently hardcoded, not sure if this should be a command line argument + print("--- Flattening the inlet/outlet if needed\n") + check_flatten_boundary(num_inlets_outlets, file_name_hdf5_mesh, threshold_stdev=0.001) + # Evaluate edge length for inspection print("--- Evaluating edge length\n") edge_length_evaluator(file_name_xml_mesh, file_name_edge_length_xdmf) diff --git a/src/vasp/automatedPreprocessing/check_flatten_in_out.py b/src/vasp/automatedPreprocessing/check_flatten_in_out.py deleted file mode 100755 index e5024ac..0000000 --- a/src/vasp/automatedPreprocessing/check_flatten_in_out.py +++ /dev/null @@ -1,144 +0,0 @@ -import argparse -import h5py -from pathlib import Path -from shutil import copyfile -import numpy as np -import os - - -def parse_arguments() -> argparse.Namespace: - """ - Parse command line arguments. - - Returns: - argparse.Namespace: Parsed command-line arguments. - """ - parser = argparse.ArgumentParser( - description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter - ) - parser.add_argument( - "--mesh-path", type=Path, default=None, help="Path to the mesh file" - ) - parser.add_argument( - "--num-inlets-outlets", - type=int, - help="Combined number of inlets and outlets (i.e, input 2 for 1 inlet and 1 outlet)", - ) - return parser.parse_args() - - -def check_flatten_boundary(num_inlets_outlets, mesh_path, threshold_stdev=0.001): - """ - Check whether inlets and outlets are flat, then flatten them if necessary - - Returns: - .h5 file of mesh with flattened outlets (flat_in_out_mesh_path) - """ - flat_in_out_mesh_path = Path(str(mesh_path).replace(".h5", "_flat_outlet.h5")) - copyfile(mesh_path, flat_in_out_mesh_path) - delete_mesh = True # For later, if outlets are already flat delete the mesh at - # flat_in_out_mesh_path - - vectorData = h5py.File(flat_in_out_mesh_path, "a") - facet_ids = np.array(vectorData["boundaries/values"]) - facet_topology = vectorData["boundaries/topology"] - - for inlet_id in range(2, 2 + num_inlets_outlets): - inlet_facet_ids = [i for i, x in enumerate(facet_ids) if x == inlet_id] - inlet_facet_topology = facet_topology[inlet_facet_ids, :] - inlet_nodes = np.unique(inlet_facet_topology.flatten()) - # pre-allocate arrays - inlet_facet_normals = np.zeros((len(inlet_facet_ids), 3)) - - # From: https://stackoverflow.com/questions/53698635/ - # how-to-define-a-plane-with-3-points-and-plot-it-in-3d - for idx, facet in enumerate(inlet_facet_topology): - p0 = vectorData["boundaries/coordinates"][facet[0]] - p1 = vectorData["boundaries/coordinates"][facet[1]] - p2 = vectorData["boundaries/coordinates"][facet[2]] - - x0, y0, z0 = p0 - x1, y1, z1 = p1 - x2, y2, z2 = p2 - - ux, uy, uz = [x1 - x0, y1 - y0, z1 - z0] # Vectors - vx, vy, vz = [x2 - x0, y2 - y0, z2 - z0] - - # cross product of vectors defines the plane normal - u_cross_v = [uy * vz - uz * vy, uz * vx - ux * vz, ux * vy - uy * vx] - normal = np.array(u_cross_v) - - # Facet unit normal vector (u_normal) - u_normal = normal / np.sqrt( - normal[0] ** 2 + normal[1] ** 2 + normal[2] ** 2 - ) - - # check if facet unit normal vector has opposite - # direction and reverse the vector if necessary - if idx == 0: - u_normal_baseline = u_normal - else: - angle = np.arccos( - np.clip(np.dot(u_normal_baseline, u_normal), -1.0, 1.0) - ) - if angle > np.pi / 2: - u_normal = -u_normal - - # record u_normal - inlet_facet_normals[idx, :] = u_normal - - # Average normal and d (we will assign this later to all facets) - normal_avg = np.mean(inlet_facet_normals, axis=0) - inlet_coords = np.array(vectorData["boundaries/coordinates"][inlet_nodes]) - point_avg = np.mean(inlet_coords, axis=0) - d_avg = -point_avg.dot(normal_avg) # plane coefficient - - # Standard deviation of components of normal vector - normal_stdev = np.std(inlet_facet_normals, axis=0) - if np.max(normal_stdev) > threshold_stdev: # if surface is not flat - print( - "Surface with ID {} is not flat: Standard deviation of facet unit\ -normals is {}, greater than threshold of {}".format( - inlet_id, np.max(normal_stdev), threshold_stdev - ) - ) - - # Move the inlet nodes into the average inlet plane (do same for outlets) - ArrayNames = [ - "boundaries/coordinates", - "mesh/coordinates", - "domains/coordinates", - ] - print("Moving nodes into a flat plane") - for ArrayName in ArrayNames: - vectorArray = vectorData[ArrayName] - for node_id in range(len(vectorArray)): - if node_id in inlet_nodes: - # from https://stackoverflow.com/questions/9605556/ - # how-to-project-a-point-onto-a-plane-in-3d (bobobobo) - node = vectorArray[node_id, :] - scalar_distance = node.dot(normal_avg) + d_avg - node_inplane = node - scalar_distance * normal_avg - vectorArray[node_id, :] = node_inplane - delete_mesh = False - - vectorData.close() - if delete_mesh is True: - print( - "Outlets and Inlets are all flat (Standard deviation of facet unit\ -normals is less than {})".format( - threshold_stdev - ) - ) - os.remove(flat_in_out_mesh_path) - - -def main(): - args = parse_arguments() - check_flatten_boundary( - args.num_inlets_outlets, args.mesh_path, threshold_stdev=0.001 - ) - - -if __name__ == "__main__": - main() diff --git a/src/vasp/automatedPreprocessing/preprocessing_common.py b/src/vasp/automatedPreprocessing/preprocessing_common.py index a003e15..1f6c432 100644 --- a/src/vasp/automatedPreprocessing/preprocessing_common.py +++ b/src/vasp/automatedPreprocessing/preprocessing_common.py @@ -4,6 +4,7 @@ from pathlib import Path from typing import Union +import h5py import numpy as np import meshio from dolfin import Mesh, MeshFunction, File, HDF5File, FunctionSpace, Function, XDMFFile, cells, Edge @@ -11,6 +12,7 @@ from morphman import vmtkscripts, write_polydata from vasp.automatedPreprocessing.vmtkmeshgeneratorfsi import vmtkMeshGeneratorFsi +from vasp.simulations.simulation_common import load_mesh_and_data from vtk import vtkPolyData @@ -326,3 +328,132 @@ def edge_length_evaluator(file_name_mesh: Union[str, Path], file_name_edge_lengt with XDMFFile(str(file_name_edge_length_xdmf)) as xdmf: xdmf.write_checkpoint(u, "edge_length", 0, append=False) + + +def check_flatten_boundary(num_inlets_outlets: int, mesh_path: Union[str, Path], + threshold_stdev: float = 0.001) -> None: + """ + Check whether inlets and outlets are flat, then flatten them if necessary + + Args: + num_inlets_outlets (int): Number of inlets and outlets in the mesh + mesh_path (Union[str, Path]): Path to the mesh file + threshold_stdev (float): Threshold for standard deviation of facet unit normals + + Returns: + None + """ + mesh_path = Path(mesh_path) + flat_in_out_mesh_path = Path(str(mesh_path).replace(".h5", "_flat_outlet.h5")) + # copy the mesh to a new file + flat_in_out_mesh_path.write_bytes(mesh_path.read_bytes()) + + vectorData = h5py.File(str(flat_in_out_mesh_path), "a") + facet_ids = np.array(vectorData["boundaries/values"]) + facet_topology = vectorData["boundaries/topology"] + + fix = False + for inlet_id in range(2, 2 + num_inlets_outlets): + inlet_facet_ids = [i for i, x in enumerate(facet_ids) if x == inlet_id] + inlet_facet_topology = facet_topology[inlet_facet_ids, :] + inlet_nodes = np.unique(inlet_facet_topology.flatten()) + # pre-allocate arrays + inlet_facet_normals = np.zeros((len(inlet_facet_ids), 3)) + + # From: https://stackoverflow.com/questions/53698635/ + # how-to-define-a-plane-with-3-points-and-plot-it-in-3d + for idx, facet in enumerate(inlet_facet_topology): + p0 = vectorData["boundaries/coordinates"][facet[0]] + p1 = vectorData["boundaries/coordinates"][facet[1]] + p2 = vectorData["boundaries/coordinates"][facet[2]] + + x0, y0, z0 = p0 + x1, y1, z1 = p1 + x2, y2, z2 = p2 + + ux, uy, uz = [x1 - x0, y1 - y0, z1 - z0] # Vectors + vx, vy, vz = [x2 - x0, y2 - y0, z2 - z0] + + # cross product of vectors defines the plane normal + u_cross_v = [uy * vz - uz * vy, uz * vx - ux * vz, ux * vy - uy * vx] + normal = np.array(u_cross_v) + + # Facet unit normal vector (u_normal) + u_normal = normal / np.sqrt( + normal[0] ** 2 + normal[1] ** 2 + normal[2] ** 2 + ) + + # check if facet unit normal vector has opposite + # direction and reverse the vector if necessary + if idx == 0: + u_normal_baseline = u_normal + else: + angle = np.arccos( + np.clip(np.dot(u_normal_baseline, u_normal), -1.0, 1.0) + ) + if angle > np.pi / 2: + u_normal = -u_normal + + # record u_normal + inlet_facet_normals[idx, :] = u_normal + + # Average normal and d (we will assign this later to all facets) + normal_avg = np.mean(inlet_facet_normals, axis=0) + inlet_coords = np.array(vectorData["boundaries/coordinates"][inlet_nodes]) + point_avg = np.mean(inlet_coords, axis=0) + d_avg = -point_avg.dot(normal_avg) # plane coefficient + + # Standard deviation of components of normal vector + normal_stdev = np.std(inlet_facet_normals, axis=0) + if np.max(normal_stdev) > threshold_stdev: # if surface is not flat + print( + "Surface with ID {} is not flat: Standard deviation of facet unit\ +normals is {}, greater than threshold of {}".format( + inlet_id, np.max(normal_stdev), threshold_stdev + ) + ) + + # Move the inlet nodes into the average inlet plane (do same for outlets) + ArrayNames = [ + "boundaries/coordinates", + "mesh/coordinates", + "domains/coordinates", + ] + print("Moving nodes into a flat plane") + for ArrayName in ArrayNames: + vectorArray = vectorData[ArrayName] + for node_id in range(len(vectorArray)): + if node_id in inlet_nodes: + # from https://stackoverflow.com/questions/9605556/ + # how-to-project-a-point-onto-a-plane-in-3d (bobobobo) + node = vectorArray[node_id, :] + scalar_distance = node.dot(normal_avg) + d_avg + node_inplane = node - scalar_distance * normal_avg + vectorArray[node_id, :] = node_inplane + fix = True + + if fix: + vectorData.close() + # Replace the original mesh file with the modified one + mesh_path.unlink() + mesh_path.write_bytes(flat_in_out_mesh_path.read_bytes()) + print("Changes made to the mesh file") + flat_in_out_mesh_path.unlink() + + # overwrite Paraview files for domains and boundaries + boundary_file = File(str(mesh_path.with_name(mesh_path.stem + "_boundaries.pvd"))) + domain_file = File(str(mesh_path.with_name(mesh_path.stem + "_domains.pvd"))) + + mesh, boundaries, domains = load_mesh_and_data(mesh_path) + boundary_file << boundaries + domain_file << domains + + else: + print( + "Surface with ID {} is flat: Standard deviation of facet unit\ +normals is {}, less than threshold of {}".format( + inlet_id, np.max(normal_stdev), threshold_stdev + )) + vectorData.close() + flat_in_out_mesh_path.unlink() + print("No changes made to the mesh file")