Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added debris_flow.mp4
Binary file not shown.
13 changes: 13 additions & 0 deletions pydebflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,15 @@ def cmd_simulate(args):
visualize=not args.no_viz
)
elif args.dem:
# Parse polygon vertices if provided
release_vertices = None
if hasattr(args, 'release_polygon') and args.release_polygon:
coords = [int(x.strip()) for x in args.release_polygon.split(',')]
if len(coords) % 2 != 0:
print("Error: --release-polygon must have even number of values")
sys.exit(1)
release_vertices = [(coords[i], coords[i+1]) for i in range(0, len(coords), 2)]

run_dem_simulation(
dem_file=args.dem,
output_dir=args.output,
Expand All @@ -45,6 +54,7 @@ def cmd_simulate(args):
release_j=args.release_col,
release_radius=args.release_radius,
release_height=args.release_height,
release_vertices=release_vertices,
animate_3d=args.animate and not args.no_viz,
export_video=args.video
)
Expand Down Expand Up @@ -162,6 +172,9 @@ def main():
help='Release zone radius in cells (default: 10)')
sim_release.add_argument('--release-height', type=float, default=5.0,
help='Release zone height in meters (default: 5.0)')
sim_release.add_argument('--release-polygon', type=str, metavar='VERTICES',
help='Polygon release zone as comma-separated row,col pairs '
'(e.g. "10,20,10,40,30,40,30,20")')

sim_viz = sim_parser.add_argument_group('Visualization')
sim_viz.add_argument('--animate', '-a', action='store_true',
Expand Down
47 changes: 38 additions & 9 deletions run_simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -272,6 +272,7 @@ def run_dem_simulation(dem_file: str,
release_j: int = None,
release_radius: int = 10,
release_height: float = 5.0,
release_vertices: list = None,
animate_3d: bool = True,
export_video: bool = False) -> None:
"""
Expand All @@ -284,6 +285,7 @@ def run_dem_simulation(dem_file: str,
release_i, release_j: Release zone center (auto if None)
release_radius: Release zone radius in cells
release_height: Initial release height (m)
release_vertices: List of (row, col) tuples for polygon release zone
animate_3d: Show 3D animation
export_video: Export animation to MP4
"""
Expand All @@ -298,12 +300,6 @@ def run_dem_simulation(dem_file: str,
print(f" Cell size: {terrain.cell_size} m")
print(f" Elevation range: {terrain.elevation.min():.1f} to {terrain.elevation.max():.1f} m")

# Auto-detect release zone if not specified
if release_i is None:
release_i = terrain.rows // 5 # Upper portion
if release_j is None:
release_j = terrain.cols // 2 # Center

# Setup model
print("\n[2/6] Configuring flow model...")
params = FlowParameters(
Expand All @@ -319,9 +315,26 @@ def run_dem_simulation(dem_file: str,
solver = NOCTVDSolver(terrain, model, solver_config)

# Initialize release
print(f"\n[3/6] Creating release zone at ({release_i}, {release_j})...")
print(f"\n[3/6] Creating release zone...")
state = FlowState.zeros((terrain.rows, terrain.cols))
release = terrain.create_release_zone(release_i, release_j, release_radius, release_height)

if release_vertices is not None and len(release_vertices) >= 3:
# Polygon release zone
release = terrain.create_polygon_release_zone(
vertices=release_vertices,
height=release_height,
smooth=True
)
print(f" Polygon release with {len(release_vertices)} vertices")
else:
# Circular release zone (original behavior)
if release_i is None:
release_i = terrain.rows // 5
if release_j is None:
release_j = terrain.cols // 2
release = terrain.create_release_zone(release_i, release_j, release_radius, release_height)
print(f" Circular release at ({release_i}, {release_j}), radius={release_radius}")

state.h_solid = release * 0.7
state.h_fluid = release * 0.3

Expand Down Expand Up @@ -457,7 +470,10 @@ def main():
release_group.add_argument('--release-radius', type=int, default=10,
help='Release zone radius in cells (default: 10)')
release_group.add_argument('--release-height', type=float, default=5.0,
help='Release zone height in meters (default: 5.0)')
help='Release zone height in meters (default: 5.0)')
release_group.add_argument('--release-polygon', type=str, metavar='VERTICES',
help='Polygon release zone as comma-separated row,col pairs '
'(e.g. "10,20,10,40,30,40,30,20")')

# Visualization
viz_group = parser.add_argument_group('Visualization')
Expand Down Expand Up @@ -485,6 +501,18 @@ def main():

# DEM simulation
if args.dem_file:
# Parse polygon vertices if provided
release_vertices = None
if args.release_polygon:
coords = [int(x.strip()) for x in args.release_polygon.split(',')]
if len(coords) % 2 != 0:
print("Error: --release-polygon must have even number of values (row,col pairs)")
sys.exit(1)
release_vertices = [(coords[i], coords[i+1]) for i in range(0, len(coords), 2)]
if len(release_vertices) < 3:
print("Error: --release-polygon requires at least 3 vertex pairs")
sys.exit(1)

run_dem_simulation(
dem_file=args.dem_file,
output_dir=args.output_dir,
Expand All @@ -493,6 +521,7 @@ def main():
release_j=args.release_col,
release_radius=args.release_radius,
release_height=args.release_height,
release_vertices=release_vertices,
animate_3d=args.animate_3d and not args.no_viz,
export_video=args.export_video
)
Expand Down
77 changes: 77 additions & 0 deletions src/core/terrain.py
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,83 @@ def create_release_zone(self, center_i: int, center_j: int,

return release

def create_polygon_release_zone(self, vertices: list, height: float,
smooth: bool = True) -> np.ndarray:
"""
Create a release zone from a polygon defined by vertices.

Args:
vertices: List of (row, col) tuples defining the polygon boundary.
At least 3 vertices required.
height: Maximum release height in meters.
smooth: If True, apply parabolic falloff from centroid.
If False, uniform height within polygon.

Returns:
Release height array with shape (rows, cols).
"""
if len(vertices) < 3:
raise ValueError("Polygon requires at least 3 vertices")

from matplotlib.path import Path as MplPath

# Create grid of all cell coordinates
rows_idx, cols_idx = np.mgrid[0:self.rows, 0:self.cols]
points = np.column_stack([rows_idx.ravel(), cols_idx.ravel()])

# Create polygon path and test containment
poly_path = MplPath(vertices)
mask = poly_path.contains_points(points).reshape(self.rows, self.cols)

return self.create_mask_release_zone(mask, height, smooth=smooth)

def create_mask_release_zone(self, mask: np.ndarray, height: float,
smooth: bool = True) -> np.ndarray:
"""
Create a release zone from a boolean mask.

Args:
mask: Boolean array of shape (rows, cols). True = release zone.
height: Maximum release height in meters.
smooth: If True, apply parabolic falloff from centroid.
If False, uniform height within mask.

Returns:
Release height array with shape (rows, cols).
"""
if mask.shape != (self.rows, self.cols):
raise ValueError(
f"Mask shape {mask.shape} doesn't match terrain ({self.rows}, {self.cols})"
)

release = np.zeros_like(self.elevation)

if not mask.any():
return release

if smooth:
# Compute centroid of masked region
masked_coords = np.argwhere(mask)
centroid_i = masked_coords[:, 0].mean()
centroid_j = masked_coords[:, 1].mean()

# Max distance from centroid to any masked cell
dists = np.sqrt(
(masked_coords[:, 0] - centroid_i)**2 +
(masked_coords[:, 1] - centroid_j)**2
)
max_dist = dists.max() if dists.max() > 0 else 1.0

# Apply parabolic falloff
for idx in range(len(masked_coords)):
i, j = masked_coords[idx]
d = dists[idx]
release[i, j] = height * (1 - (d / max_dist)**2)
else:
release[mask] = height

return release

@classmethod
def from_geotiff(cls, filepath: str) -> 'Terrain':
"""Load terrain from GeoTIFF file."""
Expand Down
Loading
Loading