openfoam-hydraulic-validation/generate_spillway_stl.py at main · sgoudarzi/openfoam-hydraulic-validation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#!/usr/bin/env python3
"""
generate_spillway_stl.py
Generates a watertight STL of a WES-standard ogee spillway,
2D-extruded for OpenFOAM single-cell-spanwise simulation.

Geometry (after USBR Design of Small Dams, Chapter 9):
  - Vertical upstream face from base (z=0) to crest (z=P)
  - WES downstream profile: z = P - 0.5*Hd*(x/Hd)**1.85, x from 0 to 2.5*Hd
  - Tangent extension from end of WES to base
  - Closed at base by horizontal segment back to upstream face

Important: the STL is extruded in y BEYOND the OpenFOAM domain
(y_stl: -0.05 to 0.15, domain y: 0 to 0.1).  This ensures the closed
body's y-normal faces do not coincide with the bounding-box
'frontAndBack' empty patches, so snappyHexMesh handles cell removal
cleanly without patch-type conflicts.

Usage:  pip install numpy-stl numpy
        python3 generate_spillway_stl.py
        mv spillway.stl constant/triSurface/
"""
import numpy as np
from stl import mesh

# --- USER PARAMETERS ----------------------------------------------------
Hd       = 1.0      # design head [m]
P        = 4.0      # spillway height [m]
n_crest  = 80       # discretization of WES downstream profile
y_min    = 0.0    # extrude beyond domain (domain y: 0 to 0.1)
y_max    = 0.1
out_file = 'spillway.stl'

# --- BUILD 2D POLYGON IN xz-PLANE (ordered CCW viewed from +y) ---------
poly = []
poly.append([0.0, 0.0])              # base of upstream face
poly.append([0.0, P])                # crest apex
for i in range(1, n_crest + 1):      # WES profile
    x = 2.5 * Hd * i / n_crest
    z = P - 0.5 * Hd * (x / Hd) ** 1.85
    poly.append([x, z])
# Tangent line from end of WES to base
x_end = 2.5 * Hd
z_end = P - 0.5 * Hd * (x_end / Hd) ** 1.85
slope = -0.5 * 1.85 * (x_end / Hd) ** 0.85 / Hd     # dz/dx at x_end
x_base_ds = x_end - z_end / slope                    # where tangent hits z=0
poly.append([x_base_ds, 0.0])

poly = np.array(poly)
N = len(poly)
print(f'Polygon: {N} vertices, x in [{poly[:,0].min():.3f}, {poly[:,0].max():.3f}], '
      f'z in [{poly[:,1].min():.3f}, {poly[:,1].max():.3f}]')

# --- TRIANGULATE FRONT (y_min) AND BACK (y_max) FACES ------------------
front, back = [], []
for i in range(1, N - 1):
    front.append([
        [poly[0,   0], y_min, poly[0,   1]],
        [poly[i+1, 0], y_min, poly[i+1, 1]],
        [poly[i,   0], y_min, poly[i,   1]],
    ])
    back.append([
        [poly[0,   0], y_max, poly[0,   1]],
        [poly[i,   0], y_max, poly[i,   1]],
        [poly[i+1, 0], y_max, poly[i+1, 1]],
    ])

# --- TRIANGULATE SIDE STRIPS (one quad per polygon edge) ---------------
sides = []
for i in range(N):
    j = (i + 1) % N
    p0 = [poly[i, 0], y_min, poly[i, 1]]
    p1 = [poly[j, 0], y_min, poly[j, 1]]
    p2 = [poly[j, 0], y_max, poly[j, 1]]
    p3 = [poly[i, 0], y_max, poly[i, 1]]
    sides.append([p0, p1, p2])
    sides.append([p0, p2, p3])

# --- SAVE STL -----------------------------------------------------------
all_tris = np.array(front + back + sides, dtype=np.float64)
sm = mesh.Mesh(np.zeros(len(all_tris), dtype=mesh.Mesh.dtype))
for k, tri in enumerate(all_tris):
    sm.vectors[k] = tri
sm.save(out_file)
print(f'Wrote {len(all_tris)} triangles to {out_file}')
print(f'STL extent: y in [{y_min}, {y_max}] (domain y: 0 to 0.1)')