From 24b34062e46dde9a887518a81bfa6bac0d655c0e Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Mon, 10 Feb 2025 14:50:00 +0100 Subject: [PATCH 1/6] improve refined and smooth mesh script --- .../refine_and_smooth_mesh.py | 66 ++++++++++++++++++- 1 file changed, 65 insertions(+), 1 deletion(-) diff --git a/creating_geometric_models/refine_and_smooth_mesh.py b/creating_geometric_models/refine_and_smooth_mesh.py index 635ad31..e88d5d8 100755 --- a/creating_geometric_models/refine_and_smooth_mesh.py +++ b/creating_geometric_models/refine_and_smooth_mesh.py @@ -3,6 +3,7 @@ import argparse import os import trimesh +import numpy as np parser = argparse.ArgumentParser(description="refine and smooth") parser.add_argument("input_file", help="surface in ts file format") @@ -22,19 +23,82 @@ required=True, help="number of and refine/smoothing steps", ) + +parser.add_argument( + "--remesh_first", + nargs=1, + metavar=("mesh_size"), + type=float, + required=True, + help="remesh first the surface with pygalmesh", +) + +parser.add_argument( + "--fix_boundary", + dest="fix_boundary", + action="store_true", + help="fix boundary edge when interpolating", +) + + args = parser.parse_args() fid = open(args.input_file) myFace = Face.from_ts(fid) -a = trimesh.Trimesh(vertices=myFace.vertex, faces=myFace.connect) +if args.remesh_first: + """ + # pygalmesh installed with: + sudo apt install libcgal-dev libeigen3-dev + pip install pygalmesh + """ + import pygalmesh + import meshio + + mesh = meshio.Mesh(points=myFace.vertex, cells=[("triangle", myFace.connect)]) + mesh.write("trash_me.vtk") + mesh_size = args.remesh_first[0] + print(f"remeshing with pygalmesh aiming for mesh size {mesh_size}") + mesh = pygalmesh.remesh_surface( + "trash_me.vtk", + max_edge_size_at_feature_edges=mesh_size, + min_facet_angle=25, + max_radius_surface_delaunay_ball=mesh_size, + max_facet_distance=mesh_size, + verbose=False, + ) + os.remove("trash_me.vtk") + print("done remeshing") + a = trimesh.Trimesh(vertices=mesh.points, faces=mesh.cells[0].data) +else: + a = trimesh.Trimesh(vertices=myFace.vertex, faces=myFace.connect) + for i in range(args.N[0]): a = a.subdivide() for i in range(args.P[0]): a = a.subdivide() + + if args.fix_boundary: + # Identify boundary edges using grouping + unique_edges = a.edges[ + trimesh.grouping.group_rows(a.edges_sorted, require_count=1) + ] + # Extract unique vertices on boundary + boundary_vertices = np.unique(unique_edges) + a_before_smoothing = a.copy() + a = trimesh.smoothing.filter_taubin(a) + if args.fix_boundary: + # Replace only the internal vertices + a.vertices[np.isin(np.arange(len(a.vertices)), boundary_vertices)] = ( + a_before_smoothing.vertices[ + np.isin(np.arange(len(a.vertices)), boundary_vertices) + ] + ) + + myFace = Face(a.vertices, a.faces) basename, ext = os.path.splitext(args.input_file) myFace.write(f"{basename}_refined_smooth_{args.N[0]}{ext}") From c6a90339027a94698273c39bda46586e65c5addf Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Mon, 10 Feb 2025 14:57:09 +0100 Subject: [PATCH 2/6] remove required --- creating_geometric_models/refine_and_smooth_mesh.py | 1 - 1 file changed, 1 deletion(-) diff --git a/creating_geometric_models/refine_and_smooth_mesh.py b/creating_geometric_models/refine_and_smooth_mesh.py index e88d5d8..a935e4b 100755 --- a/creating_geometric_models/refine_and_smooth_mesh.py +++ b/creating_geometric_models/refine_and_smooth_mesh.py @@ -29,7 +29,6 @@ nargs=1, metavar=("mesh_size"), type=float, - required=True, help="remesh first the surface with pygalmesh", ) From d5f11e4b86bfec5675e957fda991343614d373e2 Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Wed, 28 Jan 2026 10:17:33 +0100 Subject: [PATCH 3/6] small improvements to Grid --- creating_geometric_models/Face.py | 75 ++++++++++++++----- creating_geometric_models/Grid.py | 25 +++++-- .../create_surface_from_rectilinear_grid.py | 20 ++--- 3 files changed, 87 insertions(+), 33 deletions(-) diff --git a/creating_geometric_models/Face.py b/creating_geometric_models/Face.py index d6cf998..44e5ec3 100644 --- a/creating_geometric_models/Face.py +++ b/creating_geometric_models/Face.py @@ -1,15 +1,24 @@ from collections import defaultdict import numpy as np +from tqdm import tqdm, trange +import pyvista as pv class Face: def __init__(self, vertex, connect): - self.connect = connect self.vertex = vertex self.local_vid_lookup = {} self.ntriangles = self.connect.shape[0] + @classmethod + def from_pyvista(self, fname): + mesh = pv.read(fname) + mesh = mesh.clean(tolerance=0.1) + faces = mesh.faces.reshape((-1, 4))[:, 1:] + myFace = Face(mesh.points, faces) + return myFace + @classmethod def from_ts(self, fid): surface_name = "undefined" @@ -72,7 +81,9 @@ def proj(self, sProj): transformer = Transformer.from_crs("epsg:4326", sProj, always_xy=True) print("projecting the node coordinates") - self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(self.vertex[:, 0], self.vertex[:, 1]) + self.vertex[:, 0], self.vertex[:, 1] = transformer.transform( + self.vertex[:, 0], self.vertex[:, 1] + ) print(self.vertex) print("done projecting") @@ -82,7 +93,9 @@ def convert_projected_to_latlon(self, sProj): transformer = Transformer.from_crs(sProj, "epsg:4326", always_xy=True) print("convert the node coordinates to lat/lon") - self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(self.vertex[:, 0], self.vertex[:, 1]) + self.vertex[:, 0], self.vertex[:, 1] = transformer.transform( + self.vertex[:, 0], self.vertex[:, 1] + ) print(self.vertex) print("done converting") @@ -109,11 +122,15 @@ def enforce_min_depth(self, depth_min): mesh = trimesh.Trimesh(self.vertex, self.connect) assert depth_min > 0 # list vertex of the face boundary - unique_edges = mesh.edges[trimesh.grouping.group_rows(mesh.edges_sorted, require_count=1)] - if len(unique_edges)==0: - raise ValueError("the surface mesh has no boundary edges. " + unique_edges = mesh.edges[ + trimesh.grouping.group_rows(mesh.edges_sorted, require_count=1) + ] + if len(unique_edges) == 0: + raise ValueError( + "the surface mesh has no boundary edges. " "This means you extracted an incorrect set of surfaces " - "(the surface mesh of a closed region).") + "(the surface mesh of a closed region)." + ) boundary_vid = np.array(list(set(list(unique_edges.flatten())))) # mask shallow water region shallow = np.where(mesh.vertices[:, 2] > -depth_min) @@ -142,17 +159,38 @@ def __writeTs(self, fname, write_full_vertex_array=True, append=False): mode = "a" if append else "w" with open(fname, mode) as fout: - fout.write("GOCAD TSURF 1\nHEADER {\nname:%s\nborder: true\nmesh: false\n*border*bstone: true\n}\nTFACE\n" % (fname)) - for ivx in vertex_id_2_write: - fout.write("VRTX %s %s %s %s\n" % (ivx, self.vertex[ivx - 1, 0], self.vertex[ivx - 1, 1], self.vertex[ivx - 1, 2])) - - for i in range(self.ntriangles): - fout.write("TRGL %d %d %d\n" % (self.connect[i, 0] + 1, self.connect[i, 1] + 1, self.connect[i, 2] + 1)) + fout.write( + "GOCAD TSURF 1\nHEADER {\nname:%s\nborder: true\nmesh: false\n*border*bstone: true\n}\nTFACE\n" + % (fname) + ) + for ivx in tqdm(vertex_id_2_write): + fout.write( + "VRTX %s %s %s %s\n" + % ( + ivx, + self.vertex[ivx - 1, 0], + self.vertex[ivx - 1, 1], + self.vertex[ivx - 1, 2], + ) + ) + + for i in trange(self.ntriangles): + fout.write( + "TRGL %d %d %d\n" + % ( + self.connect[i, 0] + 1, + self.connect[i, 1] + 1, + self.connect[i, 2] + 1, + ) + ) fout.write("END\n") def __computeNormal(self): "compute efficiently the normals" - normal = np.cross(self.vertex[self.connect[:, 1], :] - self.vertex[self.connect[:, 0], :], self.vertex[self.connect[:, 2], :] - self.vertex[self.connect[:, 0], :]).astype(float) + normal = np.cross( + self.vertex[self.connect[:, 1], :] - self.vertex[self.connect[:, 0], :], + self.vertex[self.connect[:, 2], :] - self.vertex[self.connect[:, 0], :], + ).astype(float) norm = np.linalg.norm(normal, axis=1) self.normal = np.divide(normal, norm[:, None]) @@ -161,11 +199,14 @@ def __writeStl(self, fname, append=False): mode = "a" if append else "w" with open(fname, mode) as fout: fout.write(f"solid {fname}\n") - for k in range(self.ntriangles): + for k in trange(self.ntriangles): fout.write("facet normal %e %e %e\n" % tuple(self.normal[k, :])) fout.write("outer loop\n") for i in range(0, 3): - fout.write("vertex %.10e %.10e %.10e\n" % tuple(self.vertex[self.connect[k, i], :])) + fout.write( + "vertex %.10e %.10e %.10e\n" + % tuple(self.vertex[self.connect[k, i], :]) + ) fout.write("endloop\n") fout.write("endfacet\n") fout.write(f"endsolid {fname}\n") @@ -177,7 +218,7 @@ def __writebStl(self, fname): fout.seek(80) fout.write(struct.pack(" Date: Wed, 25 Feb 2026 15:16:42 +0100 Subject: [PATCH 4/6] add rotate option --- creating_geometric_models/generate_box.py | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/creating_geometric_models/generate_box.py b/creating_geometric_models/generate_box.py index f01f7c1..fb1bc60 100755 --- a/creating_geometric_models/generate_box.py +++ b/creating_geometric_models/generate_box.py @@ -46,6 +46,13 @@ default=[10e3], ) +parser.add_argument( + "--rotate", + metavar=("angle_in_degree"), + help="rotate box around it center and Uz", + type=float, +) + parser.add_argument( "--shrink", nargs=1, @@ -107,12 +114,26 @@ def read_lat_lon_range_from_netcdf(fname): xwidth = xmax - xmin ywidth = ymax - ymin height = args.zdim[1] - zmin + box_origin = [ xmin + 0.5 * (1 - args.shrink[0]) * xwidth, ymin + 0.5 * (1 - args.shrink[0]) * ywidth, args.zdim[0], ] box_dimension = [args.shrink[0] * xwidth, args.shrink[0] * ywidth, height] - geom.add_box(box_origin, box_dimension, args.meshSize[0]) + box = geom.add_box(box_origin, box_dimension, args.meshSize[0]) + + if args.rotate: + center = [ + box_origin[0] + 0.5 * box_dimension[0], + box_origin[1] + 0.5 * box_dimension[1], + box_origin[2] + 0.5 * box_dimension[2], + ] + + angle_rad = np.radians(args.rotate) + axis_vector = [0.0, 0.0, 1.0] + + geom.rotate(box, center, angle_rad, axis_vector) + mesh = geom.generate_mesh(dim=2) mesh.write(args.output_file) From e6d86614acafbe38f58d430696a7098486f8bc8a Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Thu, 26 Feb 2026 09:22:05 +0100 Subject: [PATCH 5/6] allow writing to vtk --- creating_geometric_models/Face.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/creating_geometric_models/Face.py b/creating_geometric_models/Face.py index 44e5ec3..625993c 100644 --- a/creating_geometric_models/Face.py +++ b/creating_geometric_models/Face.py @@ -224,6 +224,15 @@ def __writebStl(self, fname): fout.write(struct.pack("<3f", *self.vertex[self.connect[k, i], :])) fout.write(struct.pack(" Date: Thu, 26 Feb 2026 09:26:13 +0100 Subject: [PATCH 6/6] add inputCRS option --- creating_geometric_models/projTs.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/creating_geometric_models/projTs.py b/creating_geometric_models/projTs.py index 678ceed..40007c9 100755 --- a/creating_geometric_models/projTs.py +++ b/creating_geometric_models/projTs.py @@ -3,14 +3,17 @@ parser = argparse.ArgumentParser(description="project ts file") parser.add_argument("ts_file", help="ts filename") -parser.add_argument("--proj", nargs=1, metavar=("projname"), help="transform vertex array to projected system.\ - projname: name of the (projected) Coordinate Reference System (CRS) (e.g. EPSG:32646 for UTM46N)") +parser.add_argument("--inputCRS", metavar=("projname"), help="input CRS", default = "epsg:4326") +parser.add_argument("--proj", metavar=("projname"), help="transform vertex array to projected system.\ + projname: name of the (projected) Coordinate Reference System (CRS) (e.g. EPSG:32646 for UTM46N)", required=True) + + args = parser.parse_args() # set projection from pyproj import Transformer -transformer = Transformer.from_crs("epsg:4326", args.proj[0], always_xy=True) +transformer = Transformer.from_crs(args.inputCRS, args.proj, always_xy=True) # read Ts file fid = open(args.ts_file)