From 24b34062e46dde9a887518a81bfa6bac0d655c0e Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Mon, 10 Feb 2025 14:50:00 +0100 Subject: [PATCH 1/9] 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/9] 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/9] 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/9] 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/9] 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/9] 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) From 43391d388ab1163a1162b0f0a62767f237a25f0d Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Mon, 23 Mar 2026 10:58:05 +0100 Subject: [PATCH 7/9] fix imprt form py_vista in case not a PolyData --- creating_geometric_models/Face.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/creating_geometric_models/Face.py b/creating_geometric_models/Face.py index 625993c..dfa84b3 100644 --- a/creating_geometric_models/Face.py +++ b/creating_geometric_models/Face.py @@ -15,8 +15,11 @@ def __init__(self, vertex, connect): def from_pyvista(self, fname): mesh = pv.read(fname) mesh = mesh.clean(tolerance=0.1) + if mesh.faces is None: + mesh = mesh.extract_surface() faces = mesh.faces.reshape((-1, 4))[:, 1:] myFace = Face(mesh.points, faces) + return myFace @classmethod @@ -232,7 +235,6 @@ def __writeVtk(self, fname): mesh = pv.PolyData(self.vertex, faces) mesh.save(fname) - def write(self, fname, write_full_vertex_array=True, append=False): import os From c7e755bdd68ce12d257be52d46e2f8cf3aea862e Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Mon, 23 Mar 2026 11:20:04 +0100 Subject: [PATCH 8/9] refactoring --- creating_geometric_models/Face.py | 20 ++++++++++++++++--- .../refine_and_smooth_mesh.py | 14 ++++--------- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/creating_geometric_models/Face.py b/creating_geometric_models/Face.py index dfa84b3..3fd815a 100644 --- a/creating_geometric_models/Face.py +++ b/creating_geometric_models/Face.py @@ -1,7 +1,8 @@ -from collections import defaultdict +import os + import numpy as np -from tqdm import tqdm, trange import pyvista as pv +from tqdm import tqdm, trange class Face: @@ -11,6 +12,20 @@ def __init__(self, vertex, connect): self.local_vid_lookup = {} self.ntriangles = self.connect.shape[0] + @classmethod + def from_file(cls, fname): + ext = os.path.splitext(fname)[1].lower() + + if ext in [".ts"]: + with open(fname, "r") as fid: + return cls.from_ts(fid) + + elif ext in [".vtk", ".vtu", ".xdmf", ".ply", ".stl", ".obj"]: + return cls.from_pyvista(fname) + + else: + raise ValueError(f"Unsupported file format: {ext}") + @classmethod def from_pyvista(self, fname): mesh = pv.read(fname) @@ -28,7 +43,6 @@ def from_ts(self, fid): xyzl = [] vid = {} trl = [] - prev_vert = -1 ivertex = 0 line = fid.readline() if not line: diff --git a/creating_geometric_models/refine_and_smooth_mesh.py b/creating_geometric_models/refine_and_smooth_mesh.py index a935e4b..b338e91 100755 --- a/creating_geometric_models/refine_and_smooth_mesh.py +++ b/creating_geometric_models/refine_and_smooth_mesh.py @@ -9,16 +9,12 @@ parser.add_argument("input_file", help="surface in ts file format") parser.add_argument( "--N", - nargs=1, - metavar=("N"), type=int, required=True, help="number of initial refine steps", ) parser.add_argument( "--P", - nargs=1, - metavar=("P"), type=int, required=True, help="number of and refine/smoothing steps", @@ -42,8 +38,7 @@ args = parser.parse_args() -fid = open(args.input_file) -myFace = Face.from_ts(fid) +myFace = Face.from_file(args.input_file) if args.remesh_first: """ @@ -72,10 +67,9 @@ else: a = trimesh.Trimesh(vertices=myFace.vertex, faces=myFace.connect) - -for i in range(args.N[0]): +for i in range(args.N): a = a.subdivide() -for i in range(args.P[0]): +for i in range(args.P): a = a.subdivide() if args.fix_boundary: @@ -100,4 +94,4 @@ myFace = Face(a.vertices, a.faces) basename, ext = os.path.splitext(args.input_file) -myFace.write(f"{basename}_refined_smooth_{args.N[0]}{ext}") +myFace.write(f"{basename}_refined_smooth_{args.N}.ts") From 148cafb64f3ca563882420863c6de6de591a7bee Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Mon, 23 Mar 2026 11:29:12 +0100 Subject: [PATCH 9/9] make name clearer --- .../create_surface_from_rectilinear_grid.py | 4 ++-- creating_geometric_models/{projTs.py => project_surface.py} | 0 2 files changed, 2 insertions(+), 2 deletions(-) rename creating_geometric_models/{projTs.py => project_surface.py} (100%) diff --git a/creating_geometric_models/create_surface_from_rectilinear_grid.py b/creating_geometric_models/create_surface_from_rectilinear_grid.py index 916b731..76da909 100755 --- a/creating_geometric_models/create_surface_from_rectilinear_grid.py +++ b/creating_geometric_models/create_surface_from_rectilinear_grid.py @@ -9,7 +9,7 @@ parser = argparse.ArgumentParser(description="create surface from a (possibly sparse, e.g. Slab2.0 dataset) structured dataset (e.g. netcdf)") parser.add_argument("input_file", help="netcdf file") parser.add_argument("output_file", help="output file (ext in stl, bstl, ts)") -parser.add_argument("--subsample" type=int, metavar=("onesample_every"), default=[1], help="use only one value every onesample_every in both direction") +parser.add_argument("--downsample", type=int, metavar=("onesample_every"), default=1, help="use only one value every onesample_every in both direction") parser.add_argument("--objectname", metavar=("objectname"), help="name of the surface in gocad") parser.add_argument("--hole", nargs=4, metavar=(("x0"), ("x1"), ("y0"), ("y1")), help="isolate a hole in surface defined by x0<=x<=x1 and y0<=y<=y1 (stl and ts output only)", type=float) parser.add_argument("--crop", nargs=4, metavar=(("x0"), ("x1"), ("y0"), ("y1")), help="select only surfaces in x0<=x<=x1 and y0<=y<=y1", type=float) @@ -31,7 +31,7 @@ else: args.objectname = args.objectname -structured_grid = Grid(args.input_file, args.subsample) +structured_grid = Grid(args.input_file, args.downsample) structured_grid.crop(args.crop) if args.change_zero_elevation: diff --git a/creating_geometric_models/projTs.py b/creating_geometric_models/project_surface.py similarity index 100% rename from creating_geometric_models/projTs.py rename to creating_geometric_models/project_surface.py