diff --git a/debris_flow.mp4 b/debris_flow.mp4 new file mode 100644 index 0000000..ed81c4c Binary files /dev/null and b/debris_flow.mp4 differ diff --git a/pydebflow.py b/pydebflow.py index 0333011..43b9ecb 100644 --- a/pydebflow.py +++ b/pydebflow.py @@ -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, @@ -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 ) @@ -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', diff --git a/run_simulation.py b/run_simulation.py index a29789b..edfebc7 100644 --- a/run_simulation.py +++ b/run_simulation.py @@ -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: """ @@ -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 """ @@ -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( @@ -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 @@ -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') @@ -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, @@ -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 ) diff --git a/src/core/terrain.py b/src/core/terrain.py index 9f6aa8b..4f29e64 100644 --- a/src/core/terrain.py +++ b/src/core/terrain.py @@ -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.""" diff --git a/src/gui/analysis_widgets.py b/src/gui/analysis_widgets.py new file mode 100644 index 0000000..3ab0e9e --- /dev/null +++ b/src/gui/analysis_widgets.py @@ -0,0 +1,477 @@ +""" +Analysis Widgets for PyDebFlow. + +Provides RAMMS-style post-simulation analysis tools: +- Cross-Section Profile: draw a line on terrain, see elevation + flow profile +- Hydrograph: click a point, see flow height / discharge over time +""" + +import numpy as np +from typing import Optional, List, Tuple + +from PyQt6.QtWidgets import ( + QWidget, QVBoxLayout, QHBoxLayout, QLabel, + QSlider, QPushButton, QSpinBox, QSplitter +) +from PyQt6.QtCore import Qt + +from matplotlib.backends.backend_qtagg import FigureCanvasQTAgg as FigureCanvas +from matplotlib.figure import Figure + + +class CrossSectionWidget(QWidget): + """ + RAMMS-style cross-section profile viewer. + + Click two points on the terrain map to define a transect line. + The lower plot shows terrain elevation + flow height along that line, + with a time slider to scrub through simulation timesteps. + """ + + def __init__(self, parent=None): + super().__init__(parent) + self.terrain = None + self.outputs = None # [(time, FlowState), ...] + self._pt1 = None + self._pt2 = None + self._current_frame = 0 + self._setup_ui() + + def _setup_ui(self): + layout = QVBoxLayout(self) + layout.setSpacing(4) + layout.setContentsMargins(4, 2, 4, 2) + + # Instructions + self.info_label = QLabel("Click two points on the terrain to define a cross-section line") + self.info_label.setStyleSheet("color: #aaa; font-size: 11px;") + layout.addWidget(self.info_label) + + # Splitter: top = terrain map, bottom = profile + splitter = QSplitter(Qt.Orientation.Vertical) + + # ── Top: Terrain Map ───────────────────────────── + self.map_figure = Figure(figsize=(5, 3), dpi=100, facecolor='#16213e') + self.map_canvas = FigureCanvas(self.map_figure) + self.map_ax = self.map_figure.add_subplot(111) + self.map_ax.set_facecolor('#1a1a2e') + self.map_ax.set_title("Click two points for cross-section", color='#e8e8e8', fontsize=9) + self._style_ax(self.map_ax) + self.map_figure.tight_layout() + self.map_canvas.mpl_connect('button_press_event', self._on_map_click) + splitter.addWidget(self.map_canvas) + + # ── Bottom: Profile Plot ───────────────────────── + self.profile_figure = Figure(figsize=(5, 2.5), dpi=100, facecolor='#16213e') + self.profile_canvas = FigureCanvas(self.profile_figure) + self.profile_ax = self.profile_figure.add_subplot(111) + self.profile_ax.set_facecolor('#1a1a2e') + self.profile_ax.set_title("Cross-Section Profile", color='#e8e8e8', fontsize=9) + self._style_ax(self.profile_ax) + self.profile_figure.tight_layout() + splitter.addWidget(self.profile_canvas) + + splitter.setStretchFactor(0, 3) + splitter.setStretchFactor(1, 2) + layout.addWidget(splitter, stretch=1) + + # ── Time Slider ────────────────────────────────── + slider_row = QHBoxLayout() + slider_row.addWidget(QLabel("Time:")) + self.time_slider = QSlider(Qt.Orientation.Horizontal) + self.time_slider.setMinimum(0) + self.time_slider.setMaximum(0) + self.time_slider.setValue(0) + self.time_slider.valueChanged.connect(self._on_slider_changed) + slider_row.addWidget(self.time_slider, stretch=1) + self.time_label = QLabel("t = 0.0 s") + self.time_label.setFixedWidth(100) + slider_row.addWidget(self.time_label) + + self.btn_reset = QPushButton("🔄 Reset Line") + self.btn_reset.setFixedWidth(90) + self.btn_reset.clicked.connect(self._reset_line) + slider_row.addWidget(self.btn_reset) + + layout.addLayout(slider_row) + + def set_data(self, terrain, outputs): + """Load terrain and simulation outputs.""" + self.terrain = terrain + self.outputs = outputs + self._pt1 = None + self._pt2 = None + self._current_frame = 0 + + if outputs: + self.time_slider.setMaximum(len(outputs) - 1) + self.time_slider.setValue(0) + + self._draw_map() + + # ── Drawing ────────────────────────────────────────── + + def _draw_map(self): + """Draw the terrain map with cross-section line.""" + self.map_ax.clear() + if self.terrain is None: + self.map_canvas.draw_idle() + return + + hillshade = self.terrain.get_hillshade() + self.map_ax.imshow(hillshade, cmap='gray', origin='upper', aspect='equal') + + # Show flow overlay for current frame + if self.outputs: + t, state = self.outputs[self._current_frame] + h = state.h_solid + state.h_fluid + masked = np.ma.masked_where(h < 0.01, h) + self.map_ax.imshow(masked, cmap='YlOrRd', alpha=0.5, origin='upper', aspect='equal') + + # Draw points and line + if self._pt1 is not None: + r, c = self._pt1 + self.map_ax.plot(c, r, 'o', color='#00ff88', markersize=8, markeredgecolor='white', markeredgewidth=1.5) + if self._pt2 is not None: + r, c = self._pt2 + self.map_ax.plot(c, r, 'o', color='#00ff88', markersize=8, markeredgecolor='white', markeredgewidth=1.5) + if self._pt1 is not None and self._pt2 is not None: + self.map_ax.plot( + [self._pt1[1], self._pt2[1]], [self._pt1[0], self._pt2[0]], + '--', color='#00ff88', linewidth=2 + ) + + self.map_ax.set_title("Terrain — click two points for cross-section", color='#e8e8e8', fontsize=9) + self._style_ax(self.map_ax) + self.map_figure.tight_layout() + self.map_canvas.draw_idle() + + def _draw_profile(self): + """Draw the cross-section profile.""" + self.profile_ax.clear() + + if self._pt1 is None or self._pt2 is None or self.terrain is None: + self.profile_ax.set_title("Select two points on the map above", color='#e8e8e8', fontsize=9) + self._style_ax(self.profile_ax) + self.profile_figure.tight_layout() + self.profile_canvas.draw_idle() + return + + # Sample along the line + n_samples = 200 + r1, c1 = self._pt1 + r2, c2 = self._pt2 + rows_line = np.linspace(r1, r2, n_samples) + cols_line = np.linspace(c1, c2, n_samples) + + # Distance along line in meters + dist = np.sqrt( + ((rows_line - r1) * self.terrain.cell_size) ** 2 + + ((cols_line - c1) * self.terrain.cell_size) ** 2 + ) + + # Sample terrain elevation + ri = np.clip(np.round(rows_line).astype(int), 0, self.terrain.rows - 1) + ci = np.clip(np.round(cols_line).astype(int), 0, self.terrain.cols - 1) + elev = self.terrain.elevation[ri, ci] + + # Plot terrain + self.profile_ax.fill_between(dist, elev.min() - 5, elev, color='#8B7355', alpha=0.4, label='_nolegend_') + self.profile_ax.plot(dist, elev, '-', color='#8B7355', linewidth=2, label='Terrain') + + # Plot flow at current timestep + if self.outputs: + t, state = self.outputs[self._current_frame] + h_total = (state.h_solid + state.h_fluid)[ri, ci] + flow_surface = elev + h_total + + mask = h_total > 0.01 + if mask.any(): + self.profile_ax.fill_between( + dist, elev, flow_surface, + where=mask, color='#ff4444', alpha=0.5, label='_nolegend_' + ) + self.profile_ax.plot(dist, flow_surface, '-', color='#ff4444', linewidth=1.5, label=f'Flow (t={t:.1f}s)') + + self.profile_ax.set_xlabel("Distance (m)", color='#aaa', fontsize=8) + self.profile_ax.set_ylabel("Elevation (m)", color='#aaa', fontsize=8) + self.profile_ax.set_title("Cross-Section Profile", color='#e8e8e8', fontsize=9) + self.profile_ax.legend(loc='upper right', fontsize=7, facecolor='#2a2a4a', edgecolor='#555', labelcolor='#ddd') + self._style_ax(self.profile_ax) + self.profile_figure.tight_layout() + self.profile_canvas.draw_idle() + + # ── Events ─────────────────────────────────────────── + + def _on_map_click(self, event): + if self.terrain is None or event.inaxes != self.map_ax: + return + col = int(round(event.xdata)) + row = int(round(event.ydata)) + row = max(0, min(row, self.terrain.rows - 1)) + col = max(0, min(col, self.terrain.cols - 1)) + + if self._pt1 is None: + self._pt1 = (row, col) + self.info_label.setText(f"Point A: ({row}, {col}) — click second point") + elif self._pt2 is None: + self._pt2 = (row, col) + length = np.sqrt( + ((self._pt2[0] - self._pt1[0]) * self.terrain.cell_size) ** 2 + + ((self._pt2[1] - self._pt1[1]) * self.terrain.cell_size) ** 2 + ) + self.info_label.setText(f"A({self._pt1[0]},{self._pt1[1]}) → B({self._pt2[0]},{self._pt2[1]}) | {length:.0f} m") + else: + # Reset and start new line + self._pt1 = (row, col) + self._pt2 = None + self.info_label.setText(f"Point A: ({row}, {col}) — click second point") + + self._draw_map() + self._draw_profile() + + def _on_slider_changed(self, value): + self._current_frame = value + if self.outputs and value < len(self.outputs): + t = self.outputs[value][0] + self.time_label.setText(f"t = {t:.1f} s") + self._draw_map() + self._draw_profile() + + def _reset_line(self): + self._pt1 = None + self._pt2 = None + self.info_label.setText("Click two points on the terrain to define a cross-section line") + self._draw_map() + self._draw_profile() + + def _style_ax(self, ax): + ax.tick_params(colors='#888', labelsize=7) + for spine in ax.spines.values(): + spine.set_color('#3d3d5c') + + +class HydrographWidget(QWidget): + """ + Hydrograph viewer — click a point on the terrain to see + flow height, velocity, and discharge over time at that location. + """ + + def __init__(self, parent=None): + super().__init__(parent) + self.terrain = None + self.outputs = None + self._monitor_point = None + self._setup_ui() + + def _setup_ui(self): + layout = QVBoxLayout(self) + layout.setSpacing(4) + layout.setContentsMargins(4, 2, 4, 2) + + # Info + self.info_label = QLabel("Click a point on the terrain to see hydrograph") + self.info_label.setStyleSheet("color: #aaa; font-size: 11px;") + layout.addWidget(self.info_label) + + # Manual entry row + entry_row = QHBoxLayout() + entry_row.addWidget(QLabel("Row:")) + self.row_spin = QSpinBox() + self.row_spin.setRange(0, 9999) + self.row_spin.setFixedWidth(70) + entry_row.addWidget(self.row_spin) + entry_row.addWidget(QLabel("Col:")) + self.col_spin = QSpinBox() + self.col_spin.setRange(0, 9999) + self.col_spin.setFixedWidth(70) + entry_row.addWidget(self.col_spin) + self.btn_set = QPushButton("📍 Set Point") + self.btn_set.setFixedWidth(90) + self.btn_set.clicked.connect(self._on_manual_set) + entry_row.addWidget(self.btn_set) + entry_row.addStretch() + layout.addLayout(entry_row) + + # Splitter: top = terrain, bottom = hydrograph plots + splitter = QSplitter(Qt.Orientation.Vertical) + + # ── Top: Terrain ───────────────────────────────── + self.map_figure = Figure(figsize=(5, 2.5), dpi=100, facecolor='#16213e') + self.map_canvas = FigureCanvas(self.map_figure) + self.map_ax = self.map_figure.add_subplot(111) + self.map_ax.set_facecolor('#1a1a2e') + self._style_ax(self.map_ax) + self.map_figure.tight_layout() + self.map_canvas.mpl_connect('button_press_event', self._on_map_click) + splitter.addWidget(self.map_canvas) + + # ── Bottom: Hydrograph Plots ───────────────────── + self.hydro_figure = Figure(figsize=(5, 3), dpi=100, facecolor='#16213e') + self.hydro_canvas = FigureCanvas(self.hydro_figure) + self.ax_height = self.hydro_figure.add_subplot(211) + self.ax_vel = self.hydro_figure.add_subplot(212) + for ax in [self.ax_height, self.ax_vel]: + ax.set_facecolor('#1a1a2e') + self._style_ax(ax) + self.hydro_figure.tight_layout() + splitter.addWidget(self.hydro_canvas) + + splitter.setStretchFactor(0, 2) + splitter.setStretchFactor(1, 3) + layout.addWidget(splitter, stretch=1) + + def set_data(self, terrain, outputs): + """Load terrain and simulation outputs.""" + self.terrain = terrain + self.outputs = outputs + self._monitor_point = None + + if terrain: + self.row_spin.setRange(0, terrain.rows - 1) + self.col_spin.setRange(0, terrain.cols - 1) + self.row_spin.setValue(terrain.rows // 2) + self.col_spin.setValue(terrain.cols // 2) + + self._draw_map() + self._draw_hydrograph() + + # ── Drawing ────────────────────────────────────────── + + def _draw_map(self): + self.map_ax.clear() + if self.terrain is None: + self.map_canvas.draw_idle() + return + + hillshade = self.terrain.get_hillshade() + self.map_ax.imshow(hillshade, cmap='gray', origin='upper', aspect='equal') + + # Show max flow extent + if self.outputs: + max_h = np.zeros((self.terrain.rows, self.terrain.cols)) + for _, state in self.outputs: + max_h = np.maximum(max_h, state.h_solid + state.h_fluid) + masked = np.ma.masked_where(max_h < 0.01, max_h) + self.map_ax.imshow(masked, cmap='YlOrRd', alpha=0.4, origin='upper', aspect='equal') + + # Draw monitor point + if self._monitor_point is not None: + r, c = self._monitor_point + self.map_ax.plot(c, r, 's', color='#00ccff', markersize=10, + markeredgecolor='white', markeredgewidth=2) + self.map_ax.annotate( + f'({r},{c})', (c, r), color='#00ccff', fontsize=8, + xytext=(8, -8), textcoords='offset points' + ) + + self.map_ax.set_title("Click to place monitor point", color='#e8e8e8', fontsize=9) + self._style_ax(self.map_ax) + self.map_figure.tight_layout() + self.map_canvas.draw_idle() + + def _draw_hydrograph(self): + self.ax_height.clear() + self.ax_vel.clear() + + if self._monitor_point is None or not self.outputs: + self.ax_height.set_title("Select a monitor point", color='#e8e8e8', fontsize=9) + self.ax_vel.set_title("", color='#e8e8e8', fontsize=9) + for ax in [self.ax_height, self.ax_vel]: + self._style_ax(ax) + self.hydro_figure.tight_layout() + self.hydro_canvas.draw_idle() + return + + r, c = self._monitor_point + times = [] + h_solid_series = [] + h_fluid_series = [] + h_total_series = [] + speed_s_series = [] + speed_f_series = [] + discharge_series = [] + + for t, state in self.outputs: + times.append(t) + hs = state.h_solid[r, c] + hf = state.h_fluid[r, c] + h_solid_series.append(hs) + h_fluid_series.append(hf) + h_total_series.append(hs + hf) + + vs = np.sqrt(state.u_solid[r, c]**2 + state.v_solid[r, c]**2) + vf = np.sqrt(state.u_fluid[r, c]**2 + state.v_fluid[r, c]**2) + speed_s_series.append(vs) + speed_f_series.append(vf) + + # Discharge = h * v * cell_width (per unit width) + v_avg = (hs * vs + hf * vf) / max(hs + hf, 1e-6) + discharge_series.append((hs + hf) * v_avg * self.terrain.cell_size) + + times = np.array(times) + + # ── Flow Height ────────────────────────────────── + self.ax_height.fill_between(times, 0, h_solid_series, color='#cc6633', alpha=0.6, label='Solid') + self.ax_height.fill_between(times, h_solid_series, h_total_series, color='#3399cc', alpha=0.6, label='Fluid') + self.ax_height.plot(times, h_total_series, '-', color='white', linewidth=1, label='Total') + self.ax_height.set_ylabel("Height (m)", color='#aaa', fontsize=8) + self.ax_height.set_title(f"Hydrograph at ({r}, {c})", color='#e8e8e8', fontsize=9) + self.ax_height.legend(loc='upper right', fontsize=6, facecolor='#2a2a4a', edgecolor='#555', labelcolor='#ddd') + self._style_ax(self.ax_height) + + # ── Velocity / Discharge ───────────────────────── + ax2 = self.ax_vel + ax2.plot(times, speed_s_series, '-', color='#cc6633', linewidth=1.2, label='v_solid') + ax2.plot(times, speed_f_series, '-', color='#3399cc', linewidth=1.2, label='v_fluid') + ax2.set_ylabel("Velocity (m/s)", color='#aaa', fontsize=8) + ax2.set_xlabel("Time (s)", color='#aaa', fontsize=8) + + # Discharge on secondary axis + ax2r = ax2.twinx() + ax2r.plot(times, discharge_series, '--', color='#ffcc00', linewidth=1, alpha=0.7, label='Q') + ax2r.set_ylabel("Discharge (m³/s)", color='#ffcc00', fontsize=8) + ax2r.tick_params(axis='y', colors='#ffcc00', labelsize=7) + ax2r.spines['right'].set_color('#ffcc00') + + # Combined legend + lines1, labels1 = ax2.get_legend_handles_labels() + lines2, labels2 = ax2r.get_legend_handles_labels() + ax2.legend(lines1 + lines2, labels1 + labels2, + loc='upper right', fontsize=6, facecolor='#2a2a4a', edgecolor='#555', labelcolor='#ddd') + self._style_ax(ax2) + + self.hydro_figure.tight_layout() + self.hydro_canvas.draw_idle() + + # ── Events ─────────────────────────────────────────── + + def _on_map_click(self, event): + if self.terrain is None or event.inaxes != self.map_ax: + return + col = int(round(event.xdata)) + row = int(round(event.ydata)) + row = max(0, min(row, self.terrain.rows - 1)) + col = max(0, min(col, self.terrain.cols - 1)) + + self._monitor_point = (row, col) + self.row_spin.setValue(row) + self.col_spin.setValue(col) + self.info_label.setText(f"Monitor point: ({row}, {col})") + self._draw_map() + self._draw_hydrograph() + + def _on_manual_set(self): + if self.terrain is None: + return + row = self.row_spin.value() + col = self.col_spin.value() + self._monitor_point = (row, col) + self.info_label.setText(f"Monitor point: ({row}, {col})") + self._draw_map() + self._draw_hydrograph() + + def _style_ax(self, ax): + ax.tick_params(colors='#888', labelsize=7) + for spine in ax.spines.values(): + spine.set_color('#3d3d5c') diff --git a/src/gui/main_window.py b/src/gui/main_window.py index ddd13fc..07dba0b 100644 --- a/src/gui/main_window.py +++ b/src/gui/main_window.py @@ -27,11 +27,12 @@ class SimulationWorker(QThread): finished = pyqtSignal(object) # outputs list error = pyqtSignal(str) - def __init__(self, terrain, params, t_end): + def __init__(self, terrain, params, t_end, release_zone=None): super().__init__() self.terrain = terrain self.params = params self.t_end = t_end + self.release_zone = release_zone self._is_cancelled = False def run(self): @@ -47,11 +48,15 @@ def run(self): # Initial state state = FlowState.zeros((self.terrain.rows, self.terrain.cols)) - release = self.terrain.create_release_zone( - self.terrain.rows // 5, - self.terrain.cols // 2, - 10, 5.0 - ) + if self.release_zone is not None: + release = self.release_zone + else: + # Fallback: auto-assign circular release zone + release = self.terrain.create_release_zone( + self.terrain.rows // 5, + self.terrain.cols // 2, + 10, 5.0 + ) state.h_solid = release * 0.7 state.h_fluid = release * 0.3 @@ -84,6 +89,7 @@ def __init__(self): self.terrain = None self.outputs = None self.worker = None + self.release_zone = None self._setup_ui() self._create_menu() @@ -270,6 +276,7 @@ def _create_viz_panel(self) -> QWidget:

Quick Start:

  1. Load DEM - Click "Load DEM..." or use "Synthetic" for testing
    2. +
  2. Mark Release Zone - Use the 🎯 Release Zone tab to click a point or draw a polygon on the terrain
  3. Configure - Adjust flow parameters or use a preset
  4. Simulate - Click "Run Simulation"
  5. Visualize - View 3D animation or export video
    6. @@ -283,6 +290,12 @@ def _create_viz_panel(self) -> QWidget: self.viz_tabs.addTab(info_widget, "ℹ️ Info") + # Release Zone tab + from src.gui.release_zone_widget import ReleaseZoneWidget + self.release_zone_widget = ReleaseZoneWidget() + self.release_zone_widget.release_zone_changed.connect(self._on_release_zone_changed) + self.viz_tabs.addTab(self.release_zone_widget, "🎯 Release Zone") + # Log tab self.log_text = QTextEdit() self.log_text.setReadOnly(True) @@ -298,6 +311,15 @@ def _create_viz_panel(self) -> QWidget: results_layout.addWidget(self.results_label) self.viz_tabs.addTab(results_widget, "📊 Results") + # Cross-Section Profile tab + from src.gui.analysis_widgets import CrossSectionWidget, HydrographWidget + self.cross_section_widget = CrossSectionWidget() + self.viz_tabs.addTab(self.cross_section_widget, "📏 Profile") + + # Hydrograph tab + self.hydrograph_widget = HydrographWidget() + self.viz_tabs.addTab(self.hydrograph_widget, "📈 Hydrograph") + layout.addWidget(self.viz_tabs) return panel @@ -575,9 +597,13 @@ def _load_dem(self): f"Elev: {self.terrain.elevation.min():.0f} - {self.terrain.elevation.max():.0f}m" ) - self.btn_run.setEnabled(True) + # Show the release zone widget with loaded terrain + self.release_zone_widget.set_terrain(self.terrain) + self.viz_tabs.setCurrentWidget(self.release_zone_widget) + self._log(f"✓ Loaded: {filepath}") - self.statusBar.showMessage(f"Loaded: {Path(filepath).name}") + self._log(" → Mark a release zone on the terrain to enable simulation") + self.statusBar.showMessage(f"Loaded: {Path(filepath).name} — Mark release zone") except Exception as e: QMessageBox.critical(self, "Error", f"Failed to load DEM:\n{e}") @@ -597,9 +623,13 @@ def _create_synthetic(self): f"Slope: 25°" ) - self.btn_run.setEnabled(True) + # Show the release zone widget with synthetic terrain + self.release_zone_widget.set_terrain(self.terrain) + self.viz_tabs.setCurrentWidget(self.release_zone_widget) + self._log("✓ Created synthetic terrain") - self.statusBar.showMessage("Synthetic terrain created") + self._log(" → Mark a release zone on the terrain to enable simulation") + self.statusBar.showMessage("Synthetic terrain created — Mark release zone") def _apply_preset(self, preset: str): """Apply parameter preset.""" @@ -618,12 +648,33 @@ def _apply_preset(self, preset: str): self.xi_spin.setValue(xi) self._log(f"Applied preset: {preset}") + def _on_release_zone_changed(self, release_zone): + """Handle release zone update from the interactive widget.""" + self.release_zone = release_zone + + if release_zone is not None and self.terrain is not None: + self.btn_run.setEnabled(True) + vol = float(release_zone.sum() * self.terrain.cell_size**2) + self._log(f"✓ Release zone set (volume ≈ {vol:.0f} m³)") + self.statusBar.showMessage("Release zone set — Ready to simulate") + else: + self.btn_run.setEnabled(False) + self.statusBar.showMessage("Mark a release zone to enable simulation") + def _run_simulation(self): """Run simulation.""" if self.terrain is None: QMessageBox.warning(self, "No Terrain", "Please load a DEM first.") return + if self.release_zone is None: + QMessageBox.warning( + self, "No Release Zone", + "Please mark a release zone on the terrain first.\n\n" + "Use the 🎯 Release Zone tab to click a point or draw a polygon." + ) + return + params = { 'solid_density': self.density_spin.value(), 'fluid_density': 1100.0, @@ -633,7 +684,8 @@ def _run_simulation(self): } self.worker = SimulationWorker( - self.terrain, params, self.time_spin.value() + self.terrain, params, self.time_spin.value(), + release_zone=self.release_zone ) self.worker.progress.connect(self._on_progress) self.worker.finished.connect(self._on_finished) @@ -681,6 +733,11 @@ def _on_finished(self, outputs): f"

    Max flow height: {max_h:.2f} m

    " f"

    Final volume: {(final.h_solid.sum() + final.h_fluid.sum()) * self.terrain.cell_size**2:.0f} m³

    " ) + + # Feed data to analysis widgets + self.cross_section_widget.set_data(self.terrain, outputs) + self.hydrograph_widget.set_data(self.terrain, outputs) + self.viz_tabs.setCurrentIndex(2) # Results tab def _on_error(self, error: str): @@ -691,7 +748,7 @@ def _on_error(self, error: str): QMessageBox.critical(self, "Simulation Error", error) def _show_3d_view(self): - """Show 3D visualization.""" + """Show 3D visualization with time slider.""" if self.outputs is None: QMessageBox.warning(self, "No Results", "Run a simulation first.") return @@ -707,12 +764,8 @@ def _show_3d_view(self): viewer.load_snapshots(snapshots, times) - # Show max height - max_h = np.zeros_like(self.terrain.elevation) - for snap in snapshots: - max_h = np.maximum(max_h, snap) - - viewer.show_static(max_h, "PyDebFlow - 3D View") + # Show interactive 3D view with time slider + viewer.show_with_slider("PyDebFlow - 3D View") except Exception as e: QMessageBox.critical(self, "Error", f"3D view failed:\n{e}") diff --git a/src/gui/release_zone_widget.py b/src/gui/release_zone_widget.py new file mode 100644 index 0000000..c17e326 --- /dev/null +++ b/src/gui/release_zone_widget.py @@ -0,0 +1,432 @@ +""" +Release Zone Interactive Widget for PyDebFlow. + +Provides an interactive matplotlib canvas embedded in PyQt6 for marking +release zones on loaded terrain. Supports: +- Point mode: Click to place a circular release zone +- Polygon mode: Click to add vertices, right-click to close +- Manual coordinate entry via spinboxes +""" + +import numpy as np +from typing import Optional, List, Tuple + +from PyQt6.QtWidgets import ( + QWidget, QVBoxLayout, QHBoxLayout, QGridLayout, + QPushButton, QLabel, QSpinBox, QDoubleSpinBox, + QGroupBox, QButtonGroup, QRadioButton, QMessageBox +) +from PyQt6.QtCore import pyqtSignal, Qt + +from matplotlib.backends.backend_qtagg import FigureCanvasQTAgg as FigureCanvas +from matplotlib.figure import Figure + + +class ReleaseZoneWidget(QWidget): + """ + Interactive widget for marking release zones on terrain. + + Embeds a matplotlib canvas showing the terrain hillshade/elevation + and allows users to mark point or polygon release zones interactively. + + Signals: + release_zone_changed(np.ndarray): Emitted when release zone is updated. + The array is the release height field (same shape as terrain). + """ + + release_zone_changed = pyqtSignal(object) # np.ndarray + + def __init__(self, parent=None): + super().__init__(parent) + + self.terrain = None + self.release_zone = None + + # Mode state + self._mode = 'point' # 'point' or 'polygon' + + # Point mode state + self._point_marker = None # (row, col) + + # Polygon mode state + self._polygon_vertices = [] # list of (row, col) + self._polygon_closed = False + + # Plot artists for overlay + self._overlay_img = None + self._marker_artists = [] + self._polygon_line = None + self._polygon_fill = None + + self._setup_ui() + self._connect_signals() + + def _setup_ui(self): + """Build the widget UI.""" + layout = QVBoxLayout(self) + layout.setSpacing(4) + layout.setContentsMargins(5, 2, 5, 2) + + # ── Top Row: Mode + Parameters (compact horizontal) ── + top_row = QHBoxLayout() + top_row.setSpacing(12) + + # Mode radios + self.btn_group = QButtonGroup(self) + + self.radio_point = QRadioButton("🎯 Point") + self.radio_point.setChecked(True) + self.radio_point.setToolTip("Click on terrain to place a circular release zone") + self.btn_group.addButton(self.radio_point) + top_row.addWidget(self.radio_point) + + self.radio_polygon = QRadioButton("📐 Polygon") + self.radio_polygon.setToolTip("Click to add vertices. Right-click to close polygon.") + self.btn_group.addButton(self.radio_polygon) + top_row.addWidget(self.radio_polygon) + + top_row.addSpacing(10) + + # Height + top_row.addWidget(QLabel("H:")) + self.height_spin = QDoubleSpinBox() + self.height_spin.setRange(0.5, 100.0) + self.height_spin.setValue(5.0) + self.height_spin.setSingleStep(0.5) + self.height_spin.setSuffix(" m") + self.height_spin.setFixedWidth(90) + top_row.addWidget(self.height_spin) + + # Radius + top_row.addWidget(QLabel("R:")) + self.radius_spin = QSpinBox() + self.radius_spin.setRange(1, 50) + self.radius_spin.setValue(10) + self.radius_spin.setToolTip("Radius for point mode") + self.radius_spin.setFixedWidth(60) + top_row.addWidget(self.radius_spin) + + top_row.addStretch() + layout.addLayout(top_row) + + # ── Second Row: Manual Coordinate Entry (compact) ──── + manual_row = QHBoxLayout() + manual_row.setSpacing(6) + + manual_row.addWidget(QLabel("Row:")) + self.row_spin = QSpinBox() + self.row_spin.setRange(0, 9999) + self.row_spin.setValue(0) + self.row_spin.setFixedWidth(70) + manual_row.addWidget(self.row_spin) + + manual_row.addWidget(QLabel("Col:")) + self.col_spin = QSpinBox() + self.col_spin.setRange(0, 9999) + self.col_spin.setValue(0) + self.col_spin.setFixedWidth(70) + manual_row.addWidget(self.col_spin) + + self.btn_add_point = QPushButton("➕ Add") + self.btn_add_point.setToolTip( + "Point mode: sets release center. Polygon mode: adds a vertex." + ) + self.btn_add_point.setFixedWidth(70) + manual_row.addWidget(self.btn_add_point) + + self.btn_close_polygon = QPushButton("✅ Close") + self.btn_close_polygon.setToolTip("Close polygon and compute release zone") + self.btn_close_polygon.setEnabled(False) + self.btn_close_polygon.setFixedWidth(70) + manual_row.addWidget(self.btn_close_polygon) + + manual_row.addStretch() + + self.status_label = QLabel("No terrain loaded") + self.status_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) + manual_row.addWidget(self.status_label) + + layout.addLayout(manual_row) + + # ── Matplotlib Canvas (takes all remaining space) ──── + self.figure = Figure(figsize=(6, 5), dpi=100, facecolor='#16213e') + self.canvas = FigureCanvas(self.figure) + self.canvas.setMinimumHeight(350) + self.ax = self.figure.add_subplot(111) + self.ax.set_facecolor('#1a1a2e') + self.ax.set_title("Load terrain to begin", color='#e8e8e8', fontsize=10) + self.ax.tick_params(colors='#888888', labelsize=8) + for spine in self.ax.spines.values(): + spine.set_color('#3d3d5c') + self.figure.tight_layout() + layout.addWidget(self.canvas, stretch=1) + + def _connect_signals(self): + """Wire up signals.""" + self.radio_point.toggled.connect(self._on_mode_changed) + self.radio_polygon.toggled.connect(self._on_mode_changed) + self.btn_add_point.clicked.connect(self._on_manual_add) + self.btn_close_polygon.clicked.connect(self._on_close_polygon) + self.canvas.mpl_connect('button_press_event', self._on_canvas_click) + + # ── Public API ────────────────────────────────────────── + + def set_terrain(self, terrain): + """ + Set the terrain to display. + + Args: + terrain: A Terrain instance with elevation, rows, cols, cell_size. + """ + self.terrain = terrain + self.release_zone = None + self._point_marker = None + self._polygon_vertices = [] + self._polygon_closed = False + + # Update spinbox ranges + self.row_spin.setRange(0, terrain.rows - 1) + self.col_spin.setRange(0, terrain.cols - 1) + self.row_spin.setValue(terrain.rows // 5) + self.col_spin.setValue(terrain.cols // 2) + + self._draw_terrain() + self.status_label.setText("Click on terrain to mark release zone") + + def get_release_zone(self) -> Optional[np.ndarray]: + """Return the current release zone array, or None if not set.""" + return self.release_zone + + def clear_zone(self): + """Clear the current release zone.""" + self.release_zone = None + self._point_marker = None + self._polygon_vertices = [] + self._polygon_closed = False + self.btn_close_polygon.setEnabled(False) + + if self.terrain is not None: + self._draw_terrain() + self.status_label.setText("Cleared — click to mark new zone") + + self.release_zone_changed.emit(None) + + # ── Internal: Drawing ─────────────────────────────────── + + def _draw_terrain(self): + """Redraw the terrain base image.""" + self.ax.clear() + + if self.terrain is None: + self.ax.set_title("No terrain loaded", color='#e8e8e8', fontsize=10) + self.canvas.draw_idle() + return + + # Show hillshade + hillshade = self.terrain.get_hillshade() + self.ax.imshow( + hillshade, cmap='gray', origin='upper', + aspect='equal', interpolation='bilinear' + ) + + # Overlay elevation contours + self.ax.contour( + self.terrain.elevation, levels=15, + colors='#00d9ff', linewidths=0.3, alpha=0.4 + ) + + self.ax.set_title("Terrain — Click to mark release zone", color='#e8e8e8', fontsize=10) + self.ax.set_xlabel("Column (j)", color='#888', fontsize=8) + self.ax.set_ylabel("Row (i)", color='#888', fontsize=8) + self.ax.tick_params(colors='#888888', labelsize=7) + for spine in self.ax.spines.values(): + spine.set_color('#3d3d5c') + + # Redraw overlays if any + self._draw_overlays() + + self.figure.tight_layout() + self.canvas.draw_idle() + + def _draw_overlays(self): + """Draw release zone overlays on top of terrain.""" + # Draw release zone heatmap + if self.release_zone is not None: + masked = np.ma.masked_where(self.release_zone < 0.01, self.release_zone) + self.ax.imshow( + masked, cmap='hot', alpha=0.55, + origin='upper', aspect='equal', interpolation='bilinear' + ) + + # Draw point marker + if self._point_marker is not None: + r, c = self._point_marker + radius = self.radius_spin.value() + circle = plt_Circle( + (c, r), radius, + fill=False, edgecolor='#00ff88', linewidth=2, linestyle='--' + ) + self.ax.add_patch(circle) + self.ax.plot(c, r, 'x', color='#00ff88', markersize=10, markeredgewidth=2) + + # Draw polygon vertices/edges + if self._polygon_vertices: + verts = self._polygon_vertices + rows = [v[0] for v in verts] + cols = [v[1] for v in verts] + + # Plot vertices + self.ax.plot(cols, rows, 'o', color='#ff6b6b', markersize=6, markeredgecolor='white', markeredgewidth=1) + + # Plot edges + if len(verts) > 1: + plot_cols = cols + ([cols[0]] if self._polygon_closed else []) + plot_rows = rows + ([rows[0]] if self._polygon_closed else []) + self.ax.plot(plot_cols, plot_rows, '-', color='#ff6b6b', linewidth=1.5) + + # Fill if closed + if self._polygon_closed and len(verts) >= 3: + from matplotlib.patches import Polygon as MplPolygon + xy = np.array([[c, r] for r, c in verts]) + poly = MplPolygon( + xy, closed=True, + facecolor='#ff6b6b', alpha=0.2, + edgecolor='#ff6b6b', linewidth=2 + ) + self.ax.add_patch(poly) + + # ── Internal: Event Handlers ──────────────────────────── + + def _on_mode_changed(self, checked): + """Handle mode radio button toggle.""" + if not checked: + return + + old_mode = self._mode + self._mode = 'point' if self.radio_point.isChecked() else 'polygon' + + if old_mode != self._mode: + # Clear current markers when switching modes + self._point_marker = None + self._polygon_vertices = [] + self._polygon_closed = False + self.btn_close_polygon.setEnabled(False) + + if self._mode == 'point': + self.radius_spin.setEnabled(True) + self.status_label.setText("Point mode — click to place release center") + else: + self.radius_spin.setEnabled(False) + self.status_label.setText("Polygon mode — click to add vertices, right-click to close") + + if self.terrain is not None: + self._draw_terrain() + + def _on_canvas_click(self, event): + """Handle mouse click on the matplotlib canvas.""" + if self.terrain is None: + return + if event.inaxes != self.ax: + return + + col = int(round(event.xdata)) + row = int(round(event.ydata)) + + # Clamp to terrain bounds + row = max(0, min(row, self.terrain.rows - 1)) + col = max(0, min(col, self.terrain.cols - 1)) + + if self._mode == 'point': + self._place_point(row, col) + elif self._mode == 'polygon': + if event.button == 3: # Right-click closes polygon + self._close_polygon() + else: + self._add_polygon_vertex(row, col) + + def _on_manual_add(self): + """Handle manual 'Add Point' button click.""" + if self.terrain is None: + QMessageBox.warning(self, "No Terrain", "Load a DEM first.") + return + + row = self.row_spin.value() + col = self.col_spin.value() + + if self._mode == 'point': + self._place_point(row, col) + else: + self._add_polygon_vertex(row, col) + + def _on_close_polygon(self): + """Handle 'Close Polygon' button click.""" + self._close_polygon() + + # ── Internal: Release Zone Logic ──────────────────────── + + def _place_point(self, row: int, col: int): + """Place a point release zone at (row, col).""" + self._point_marker = (row, col) + + height = self.height_spin.value() + radius = self.radius_spin.value() + + self.release_zone = self.terrain.create_release_zone( + center_i=row, center_j=col, + radius=radius, height=height + ) + + self.status_label.setText(f"Point at ({row}, {col}) — r={radius}, h={height:.1f}m") + self._draw_terrain() + self.release_zone_changed.emit(self.release_zone) + + def _add_polygon_vertex(self, row: int, col: int): + """Add a vertex to the polygon being drawn.""" + if self._polygon_closed: + # Start a new polygon + self._polygon_vertices = [] + self._polygon_closed = False + self.release_zone = None + + self._polygon_vertices.append((row, col)) + n = len(self._polygon_vertices) + + self.btn_close_polygon.setEnabled(n >= 3) + self.status_label.setText( + f"Polygon: {n} vertices — " + + ("right-click or Close to finish" if n >= 3 else f"need {3 - n} more") + ) + + self._draw_terrain() + + def _close_polygon(self): + """Close the current polygon and compute the release zone.""" + if len(self._polygon_vertices) < 3: + QMessageBox.warning( + self, "Not Enough Vertices", + "A polygon needs at least 3 vertices." + ) + return + + self._polygon_closed = True + height = self.height_spin.value() + + self.release_zone = self.terrain.create_polygon_release_zone( + vertices=self._polygon_vertices, + height=height, + smooth=True + ) + + n = len(self._polygon_vertices) + self.status_label.setText(f"Polygon ({n} vertices) — h={height:.1f}m ✓") + self.btn_close_polygon.setEnabled(False) + + self._draw_terrain() + self.release_zone_changed.emit(self.release_zone) + + +# Helper import for drawing circles (deferred to avoid import-time cost) +def plt_Circle(center, radius, **kwargs): + """Create a matplotlib Circle patch.""" + from matplotlib.patches import Circle + return Circle(center, radius, **kwargs) diff --git a/src/visualization/dem_viewer.py b/src/visualization/dem_viewer.py index 543a115..684b9a4 100644 --- a/src/visualization/dem_viewer.py +++ b/src/visualization/dem_viewer.py @@ -180,6 +180,177 @@ def show_static(self, flow_height: Optional[np.ndarray] = None, # Show plotter.show() + def show_with_slider(self, title: str = "PyDebFlow - 3D Flow Viewer") -> None: + """ + Show interactive 3D view with a time slider. + + Users can drag the slider to view the flow at any timestep. + The flow mesh updates in real-time as the slider moves. + + Keyboard controls: + Space: Play/Pause auto-playback + Left Arrow: Previous frame + Right Arrow: Next frame + + Requires snapshots to be loaded via load_snapshots() first. + """ + if not self.flow_snapshots: + raise ValueError("No snapshots loaded. Use load_snapshots() first.") + + pv = self._ensure_pyvista() + + n_frames = len(self.flow_snapshots) + max_height = max(s.max() for s in self.flow_snapshots) + + # Create plotter + plotter = pv.Plotter(title=title, window_size=(1400, 900)) + plotter.set_background('white', top='lightblue') + + # ── Terrain ───────────────────────────────────────── + terrain_mesh = self._create_terrain_mesh(pv) + hillshade = self._compute_hillshade() + terrain_mesh['hillshade'] = hillshade.flatten(order='F') + + plotter.add_mesh( + terrain_mesh, + scalars='hillshade', + cmap='gray', + show_scalar_bar=False, + smooth_shading=True, + name='terrain' + ) + + # ── Initial flow frame ────────────────────────────── + initial_flow = self.flow_snapshots[0] + flow_mesh = self._create_flow_mesh(initial_flow, pv) + + if flow_mesh is not None: + plotter.add_mesh( + flow_mesh, + scalars='flow_height', + cmap='YlOrRd', + clim=[0, max_height], + opacity=0.85, + smooth_shading=True, + scalar_bar_args={ + 'title': 'Flow Height (m)', + 'vertical': True, + 'position_x': 0.9, + 'position_y': 0.3, + 'width': 0.08, + 'height': 0.4, + }, + name='flow' + ) + + # ── Time text ─────────────────────────────────────── + plotter.add_text( + f"t = {self.snapshot_times[0]:.1f} s [Frame 1/{n_frames}]", + position='upper_left', + font_size=12, + color='black', + name='time_label' + ) + + # ── Camera ────────────────────────────────────────── + plotter.camera_position = 'iso' + plotter.camera.zoom(0.8) + plotter.show_axes() + + # ── State for callbacks ───────────────────────────── + viewer_self = self # Capture self for closures + state = {'frame': 0, 'playing': False} + + def update_to_frame(frame_idx): + """Update the 3D view to show a specific frame.""" + frame_idx = int(round(frame_idx)) + frame_idx = max(0, min(frame_idx, n_frames - 1)) + state['frame'] = frame_idx + + snapshot = viewer_self.flow_snapshots[frame_idx] + t = viewer_self.snapshot_times[frame_idx] + + # Remove old flow and add new one + new_flow = viewer_self._create_flow_mesh(snapshot, pv) + if new_flow is not None: + plotter.add_mesh( + new_flow, + scalars='flow_height', + cmap='YlOrRd', + clim=[0, max_height], + opacity=0.85, + smooth_shading=True, + name='flow' + ) + else: + # No flow in this frame — remove the actor + try: + plotter.remove_actor('flow') + except Exception: + pass + + # Update time text + plotter.add_text( + f"t = {t:.1f} s [Frame {frame_idx + 1}/{n_frames}]", + position='upper_left', + font_size=12, + color='black', + name='time_label' + ) + + # ── Slider widget ─────────────────────────────────── + plotter.add_slider_widget( + update_to_frame, + rng=[0, n_frames - 1], + value=0, + title="Time Step", + pointa=(0.15, 0.07), + pointb=(0.85, 0.07), + style='modern', + fmt="%.0f", + ) + + # ── Keyboard controls (optional, may not be available in older PyVista) ── + try: + def on_key_space(): + """Toggle auto-play.""" + state['playing'] = not state['playing'] + + def on_key_right(): + """Next frame.""" + new_idx = min(state['frame'] + 1, n_frames - 1) + update_to_frame(new_idx) + plotter.render() + + def on_key_left(): + """Previous frame.""" + new_idx = max(state['frame'] - 1, 0) + update_to_frame(new_idx) + plotter.render() + + plotter.add_key_event('space', on_key_space) + plotter.add_key_event('Right', on_key_right) + plotter.add_key_event('Left', on_key_left) + except AttributeError: + pass # Older PyVista — slider still works + + # ── Auto-play callback (optional, may not be available in older PyVista) ── + try: + def auto_play_tick(): + """Advance frame if playing.""" + if state['playing']: + new_idx = state['frame'] + 1 + if new_idx >= n_frames: + new_idx = 0 # Loop + update_to_frame(new_idx) + + plotter.add_callback(auto_play_tick, interval=300) + except AttributeError: + pass # Older PyVista — slider still works without auto-play + + # ── Show ──────────────────────────────────────────── + plotter.show() + def show_animation(self, title: str = "PyDebFlow - Debris Flow Animation") -> None: """ Show animated debris flow visualization. diff --git a/tests/test_api.py b/tests/test_api.py index e1c482a..e43199d 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -63,6 +63,34 @@ def test_create_release_zone(self): assert release.shape == (50, 40) assert release.max() <= 3.0 assert release.max() > 0.0 + + def test_create_polygon_release_zone(self): + """Test creating a polygon release zone.""" + from src import Terrain + + terrain = Terrain.create_synthetic_slope(rows=50, cols=40) + vertices = [(10, 10), (10, 30), (25, 30), (25, 10)] + release = terrain.create_polygon_release_zone(vertices, height=4.0) + + assert release.shape == (50, 40) + assert release.max() <= 4.0 + assert release.max() > 0.0 + # Center of polygon should have material + assert release[17, 20] > 0 + + def test_create_mask_release_zone(self): + """Test creating a mask-based release zone.""" + from src import Terrain + import numpy as np + + terrain = Terrain.create_synthetic_slope(rows=50, cols=40) + mask = np.zeros((50, 40), dtype=bool) + mask[5:15, 10:25] = True + release = terrain.create_mask_release_zone(mask, height=5.0) + + assert release.shape == (50, 40) + assert release.max() <= 5.0 + assert release.max() > 0.0 class TestFlowModelAPI: diff --git a/tests/test_release_zone.py b/tests/test_release_zone.py new file mode 100644 index 0000000..133c5c6 --- /dev/null +++ b/tests/test_release_zone.py @@ -0,0 +1,225 @@ +""" +Tests for release zone creation functionality. + +Tests the new polygon and mask-based release zone methods +added to the Terrain class. +""" + +import pytest +import numpy as np +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).parent.parent)) + +from src.core.terrain import Terrain + + +@pytest.fixture +def terrain(): + """Create a standard synthetic terrain for testing.""" + return Terrain.create_synthetic(rows=50, cols=40, cell_size=10.0, slope_angle=25.0) + + +class TestCircularReleaseZone: + """Test existing circular release zone (regression tests).""" + + def test_basic_release(self, terrain): + """Test basic circular release zone.""" + release = terrain.create_release_zone(10, 20, 5, 3.0) + assert release.shape == (50, 40) + assert release.max() <= 3.0 + assert release.max() > 0 + assert release[10, 20] == pytest.approx(3.0) # Center should be max + + def test_zero_outside_radius(self, terrain): + """Test that release is zero outside radius.""" + release = terrain.create_release_zone(25, 20, 5, 3.0) + assert release[0, 0] == 0.0 + assert release[49, 39] == 0.0 + + def test_edge_release(self, terrain): + """Test release zone at terrain edge doesn't crash.""" + release = terrain.create_release_zone(0, 0, 5, 3.0) + assert release.shape == (50, 40) + assert release.max() > 0 + + +class TestPolygonReleaseZone: + """Test polygon-based release zone creation.""" + + def test_triangle_release(self, terrain): + """Test triangular polygon release zone.""" + vertices = [(10, 15), (10, 25), (20, 20)] + release = terrain.create_polygon_release_zone(vertices, height=5.0) + + assert release.shape == (50, 40) + assert release.max() <= 5.0 + assert release.max() > 0 + # Outside polygon should be zero + assert release[0, 0] == 0.0 + assert release[49, 39] == 0.0 + + def test_rectangle_release(self, terrain): + """Test rectangular polygon release zone.""" + vertices = [(5, 10), (5, 30), (15, 30), (15, 10)] + release = terrain.create_polygon_release_zone(vertices, height=4.0) + + assert release.shape == (50, 40) + assert release.max() <= 4.0 + assert release.max() > 0 + # Center of rectangle should have material + assert release[10, 20] > 0 + + def test_polygon_uniform_height(self, terrain): + """Test polygon with uniform (non-smooth) height.""" + vertices = [(5, 10), (5, 30), (15, 30), (15, 10)] + release = terrain.create_polygon_release_zone(vertices, height=4.0, smooth=False) + + # All interior cells should have exactly height=4.0 + interior_values = release[release > 0] + assert np.allclose(interior_values, 4.0) + + def test_polygon_smooth_height(self, terrain): + """Test polygon with smooth (parabolic) height falloff.""" + vertices = [(5, 10), (5, 30), (15, 30), (15, 10)] + release = terrain.create_polygon_release_zone(vertices, height=4.0, smooth=True) + + # Should have some variation (not all the same) + interior_values = release[release > 0] + assert interior_values.std() > 0, "Smooth height should vary across polygon" + + def test_polygon_too_few_vertices(self, terrain): + """Test that < 3 vertices raises ValueError.""" + with pytest.raises(ValueError, match="at least 3"): + terrain.create_polygon_release_zone([(10, 10), (20, 20)], height=5.0) + + def test_large_polygon(self, terrain): + """Test polygon covering a large portion of terrain.""" + vertices = [(2, 2), (2, 37), (47, 37), (47, 2)] + release = terrain.create_polygon_release_zone(vertices, height=3.0) + + assert release.max() <= 3.0 + affected_cells = np.sum(release > 0) + total_cells = terrain.rows * terrain.cols + assert affected_cells > total_cells * 0.5 # Should cover > 50% + + def test_polygon_outside_bounds(self, terrain): + """Test polygon partially outside terrain bounds.""" + vertices = [(-5, 10), (-5, 30), (10, 30), (10, 10)] + release = terrain.create_polygon_release_zone(vertices, height=5.0) + + # Should still produce valid output, just clipped + assert release.shape == (50, 40) + assert not np.isnan(release).any() + + def test_pentagon_release(self, terrain): + """Test a 5-sided polygon.""" + vertices = [(10, 20), (15, 28), (25, 25), (25, 15), (15, 12)] + release = terrain.create_polygon_release_zone(vertices, height=6.0) + + assert release.shape == (50, 40) + assert release.max() <= 6.0 + assert release.max() > 0 + + +class TestMaskReleaseZone: + """Test mask-based release zone creation.""" + + def test_basic_mask(self, terrain): + """Test basic boolean mask release zone.""" + mask = np.zeros((terrain.rows, terrain.cols), dtype=bool) + mask[10:20, 15:25] = True + + release = terrain.create_mask_release_zone(mask, height=5.0) + + assert release.shape == (50, 40) + assert release.max() <= 5.0 + assert release.max() > 0 + assert release[0, 0] == 0.0 + + def test_mask_uniform(self, terrain): + """Test mask with uniform height (no smooth).""" + mask = np.zeros((terrain.rows, terrain.cols), dtype=bool) + mask[10:20, 15:25] = True + + release = terrain.create_mask_release_zone(mask, height=5.0, smooth=False) + + # All masked cells should be exactly 5.0 + assert np.allclose(release[mask], 5.0) + assert np.allclose(release[~mask], 0.0) + + def test_empty_mask(self, terrain): + """Test empty mask returns all zeros.""" + mask = np.zeros((terrain.rows, terrain.cols), dtype=bool) + release = terrain.create_mask_release_zone(mask, height=5.0) + + assert np.allclose(release, 0.0) + + def test_wrong_mask_shape(self, terrain): + """Test that wrong mask shape raises ValueError.""" + mask = np.zeros((10, 10), dtype=bool) + with pytest.raises(ValueError, match="doesn't match"): + terrain.create_mask_release_zone(mask, height=5.0) + + def test_single_cell_mask(self, terrain): + """Test single-cell mask.""" + mask = np.zeros((terrain.rows, terrain.cols), dtype=bool) + mask[25, 20] = True + + release = terrain.create_mask_release_zone(mask, height=3.0) + assert release[25, 20] > 0 + + def test_mask_smooth_has_max_at_centroid(self, terrain): + """Test that smooth mask has maximum near centroid.""" + mask = np.zeros((terrain.rows, terrain.cols), dtype=bool) + mask[10:20, 15:25] = True + + release = terrain.create_mask_release_zone(mask, height=5.0, smooth=True) + + # Centroid of 10:20, 15:25 is roughly (14.5, 19.5) + # Maximum should be near the center + max_pos = np.unravel_index(release.argmax(), release.shape) + assert 10 <= max_pos[0] <= 20 + assert 15 <= max_pos[1] <= 25 + + +class TestReleaseZoneIntegration: + """Integration tests verifying release zones work in simulations.""" + + def test_polygon_zone_in_simulation(self): + """Test that a polygon release zone can drive a simulation.""" + from src.core.flow_model import TwoPhaseFlowModel, FlowState, FlowParameters + from src.core.noc_tvd_solver import NOCTVDSolver, SolverConfig + + terrain = Terrain.create_synthetic(rows=30, cols=25, cell_size=10.0, slope_angle=25.0) + + # Create polygon release zone + vertices = [(3, 8), (3, 17), (8, 17), (8, 8)] + release = terrain.create_polygon_release_zone(vertices, height=3.0) + + # Setup simulation + params = FlowParameters(solid_density=2500.0, fluid_density=1100.0) + model = TwoPhaseFlowModel(params) + config = SolverConfig(cfl_number=0.4, max_timestep=0.5) + solver = NOCTVDSolver(terrain, model, config) + + state = FlowState.zeros((terrain.rows, terrain.cols)) + state.h_solid = release * 0.7 + state.h_fluid = release * 0.3 + + initial_volume = (state.h_solid.sum() + state.h_fluid.sum()) * terrain.cell_size**2 + + # Run short simulation + outputs = solver.run_simulation(state, t_end=2.0, output_interval=1.0) + + assert len(outputs) >= 2 + _, final = outputs[-1] + final_volume = (final.h_solid.sum() + final.h_fluid.sum()) * terrain.cell_size**2 + + # Volume should not increase significantly + assert final_volume <= initial_volume * 1.05 + + +if __name__ == '__main__': + pytest.main([__file__, '-v'])