Merge pull request #188 from keiyamamo/merge_flat_function

Merge flat function

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)

src/vasp/automatedPreprocessing/preprocessing_common.py

@@ -4,13 +4,15 @@
 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
 from vmtk import vmtkdistancetospheres, vmtkdijkstradistancetopoints
 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")