diff --git a/docs/source/notebooks/influence_lines_and_surfaces.ipynb b/docs/source/notebooks/influence_lines_and_surfaces.ipynb new file mode 100644 index 00000000..2524366d --- /dev/null +++ b/docs/source/notebooks/influence_lines_and_surfaces.ipynb @@ -0,0 +1,360 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Influence Lines and Surfaces\n", + "\n", + "This notebook demonstrates the recommended influence-analysis workflow in `ospgrillage`:\n", + "\n", + "- multi-lane influence lines in one result object/file,\n", + "- 2D and path-overlay influence-line plotting,\n", + "- station-based influence-surface extraction,\n", + "- physical-coordinate influence-surface plotting for skewed/curved geometry,\n", + "- semantic NetCDF output naming (`.il.nc` / `.is.nc`) and CSV export." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "from pathlib import Path\n", + "\n", + "cwd = Path.cwd().resolve()\n", + "for candidate in [cwd, *cwd.parents]:\n", + " src_dir = candidate / \"src\"\n", + " if (src_dir / \"ospgrillage\").exists():\n", + " sys.path.insert(0, str(src_dir))\n", + " break\n", + "\n", + "import ospgrillage as og\n", + "\n", + "try:\n", + " import plotly.io as pio\n", + " pio.renderers.default = \"notebook_connected\"\n", + "except Exception:\n", + " pass" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Technical context\n", + "\n", + "- `shape_function=\"hermite\"` uses Hermite interpolation on quadrilateral\n", + " regions and a DKT-style condensed triangular distributor on 3-node skew\n", + " regions.\n", + "- Influence lines often use `load_coord=\"station\"` because station is a\n", + " stable abscissa along a lane path even when the bridge is skewed/curved.\n", + "- Influence surfaces are built from admissible mesh station points and retain\n", + " mapped physical (`x`,`z`) coordinates for plotting/export." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Build a compact bridge model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "concrete = og.create_material(material=\"concrete\", code=\"AS5100-2017\", grade=\"50MPa\")\n", + "main_section = og.create_section(A=0.896, J=0.133, Iy=0.213, Iz=0.259, Ay=0.233, Az=0.58)\n", + "slab_section = og.create_section(A=0.04428, J=2.6e-4, Iy=1.1e-4, Iz=2.42e-4, Ay=3.69e-1, Az=3.69e-1, unit_width=True)\n", + "edge_section = og.create_section(A=0.044625, J=2.28e-3, Iy=2.23e-1, Iz=1.2e-3, Ay=3.72e-2, Az=3.72e-2)\n", + "\n", + "main_beam = og.create_member(member_name=\"Main Beam\", section=main_section, material=concrete)\n", + "slab = og.create_member(member_name=\"Slab\", section=slab_section, material=concrete)\n", + "edge_beam = og.create_member(member_name=\"Edge Beam\", section=edge_section, material=concrete)\n", + "\n", + "bridge = og.create_grillage(\n", + " bridge_name=\"Influence Demo\",\n", + " long_dim=10,\n", + " width=7,\n", + " skew=0,\n", + " num_long_grid=7,\n", + " num_trans_grid=5,\n", + " edge_beam_dist=1,\n", + " mesh_type=\"Ortho\",\n", + ")\n", + "bridge.set_member(main_beam, member=\"interior_main_beam\")\n", + "bridge.set_member(edge_beam, member=\"exterior_main_beam_1\")\n", + "bridge.set_member(edge_beam, member=\"exterior_main_beam_2\")\n", + "bridge.set_member(edge_beam, member=\"edge_beam\")\n", + "bridge.set_member(slab, member=\"transverse_slab\")\n", + "bridge.set_member(edge_beam, member=\"start_edge\")\n", + "bridge.set_member(edge_beam, member=\"end_edge\")\n", + "bridge.create_osp_model(pyfile=False)\n", + "\n", + "target_element = bridge.get_element(member=\"interior_main_beam\", options=\"elements\")[0];" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Influence lines from multiple lane paths" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "paths = {\n", + " \"Lane 1\": og.Path(start_point=og.Point(0, 0, 1.5), end_point=og.Point(10, 0, 1.5), increments=11),\n", + " \"Lane 2\": og.Path(start_point=og.Point(0, 0, 3.5), end_point=og.Point(10, 0, 3.5), increments=11),\n", + " \"Lane 3\": og.Path(start_point=og.Point(0, 0, 5.5), end_point=og.Point(10, 0, 5.5), increments=11),\n", + "}\n", + "\n", + "ils = bridge.analyze_influence_lines(\n", + " paths=paths,\n", + " axle_load=1.0,\n", + " shape_function=\"hermite\",\n", + ");" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ils.plot(\n", + " array=\"forces\",\n", + " component=\"Mz_j\",\n", + " element=target_element,\n", + " load_coord=\"station\",\n", + " title=\"Influence lines by lane (station abscissa)\",\n", + " ylabel=\"Mz_j ordinate\",\n", + ");" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig_il_path = ils.plot(\n", + " array=\"forces\",\n", + " component=\"Mz_j\",\n", + " element=target_element,\n", + " backend=\"plotly\",\n", + " view=\"path\",\n", + " load_coord=\"station\",\n", + " title=\"Influence lines overlaid on axle paths\",\n", + " show=False,\n", + ");\n", + "# Optional in an interactive notebook frontend:\n", + "# fig_il_path.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Influence surface from admissible station grid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iss = bridge.analyze_influence_surfaces(\n", + " name=\"Deck IS\",\n", + " point_load=1.0,\n", + " shape_function=\"hermite\",\n", + ");" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iss.plot(\n", + " array=\"forces\",\n", + " component=\"Mz_j\",\n", + " element=target_element,\n", + " x_coord=\"longitudinal_station\",\n", + " y_coord=\"transverse_station\",\n", + " coordinate_space=\"physical\",\n", + " title=\"Influence surface contour (physical x-z)\",\n", + ");" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig_is_3d = iss.plot(\n", + " array=\"forces\",\n", + " component=\"Mz_j\",\n", + " element=target_element,\n", + " x_coord=\"longitudinal_station\",\n", + " y_coord=\"transverse_station\",\n", + " coordinate_space=\"physical\",\n", + " backend=\"plotly\",\n", + " view=\"surface3d\",\n", + " title=\"Influence surface mapped to physical deck coordinates\",\n", + " show=False,\n", + ");\n", + "# Optional in an interactive notebook frontend:\n", + "# fig_is_3d.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Curved Deck Example\n", + "\n", + "Curved meshes follow the same API as the dedicated `curve_mesh`\n", + "notebook. The compact check below uses a moderate highway-style\n", + "curvature (`long_dim=30`, `mesh_radius=200`) and plots the influence\n", + "surface in physical `x-z` space." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "curved_bridge = og.create_grillage(\n", + " bridge_name=\"Curved Influence Demo\",\n", + " long_dim=30,\n", + " width=7,\n", + " skew=0,\n", + " num_long_grid=7,\n", + " num_trans_grid=15,\n", + " edge_beam_dist=1,\n", + " mesh_type=\"Ortho\",\n", + " mesh_radius=200,\n", + ")\n", + "curved_bridge.set_member(main_beam, member=\"interior_main_beam\")\n", + "curved_bridge.set_member(edge_beam, member=\"exterior_main_beam_1\")\n", + "curved_bridge.set_member(edge_beam, member=\"exterior_main_beam_2\")\n", + "curved_bridge.set_member(edge_beam, member=\"edge_beam\")\n", + "curved_bridge.set_member(slab, member=\"transverse_slab\")\n", + "curved_bridge.set_member(edge_beam, member=\"start_edge\")\n", + "curved_bridge.set_member(edge_beam, member=\"end_edge\")\n", + "curved_bridge.create_osp_model(pyfile=False)\n", + "\n", + "curved_nodes = curved_bridge.get_nodes()\n", + "curved_target_node = min(\n", + " curved_nodes,\n", + " key=lambda tag: (\n", + " (curved_nodes[tag][\"coordinate\"][0] - 15.0) ** 2\n", + " + (curved_nodes[tag][\"coordinate\"][2] - 3.5) ** 2\n", + " ),\n", + ")\n", + "curved_iss = curved_bridge.analyze_influence_surfaces(\n", + " name=\"Curved Deck IS\",\n", + " point_load=1.0,\n", + " shape_function=\"hermite\",\n", + ");\n", + "ax_curved = curved_iss.plot(\n", + " array=\"displacements\",\n", + " component=\"y\",\n", + " node=curved_target_node,\n", + " x_coord=\"longitudinal_station\",\n", + " y_coord=\"transverse_station\",\n", + " coordinate_space=\"physical\",\n", + " title=\"Curved-bridge IS: deck y displacement near midspan (L = 30 m, R = 200 m)\",\n", + ");\n", + "ax_curved.set_aspect(\"equal\", adjustable=\"box\");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Save combined influence outputs\n", + "\n", + "Both result objects retain save metadata and can be written directly to\n", + "single NetCDF files. If a stem is given, semantic suffixes are applied:\n", + "`*.il.nc` for influence lines and `*.is.nc` for influence surfaces." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "il_nc = ils.save(\"lane_ils\")\n", + "is_nc = iss.save(\"deck_is\")\n", + "print(f\"Saved {il_nc}\")\n", + "print(f\"Saved {is_nc}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Export CSV\n", + "\n", + "- IL CSV is long-form with one row per path station.\n", + "- IS CSV can be exported as station-grid (`layout=\"grid\"`) and optional\n", + " companion point map with physical coordinates." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "il_csv = ils.to_csv(\n", + " \"lane_ils.csv\",\n", + " array=\"forces\",\n", + " component=\"Mz_j\",\n", + " element=target_element,\n", + " load_coord=\"station\",\n", + ")\n", + "\n", + "is_csv = iss.to_csv(\n", + " \"deck_is_grid.csv\",\n", + " array=\"forces\",\n", + " component=\"Mz_j\",\n", + " element=target_element,\n", + " include_physical_coords=True,\n", + ")\n", + "\n", + "print(f\"Saved {il_csv}\")\n", + "print(f\"Saved {is_csv}\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/source/pages/api_analysis.rst b/docs/source/pages/api_analysis.rst index b28eb3d1..8cc0a7e2 100644 --- a/docs/source/pages/api_analysis.rst +++ b/docs/source/pages/api_analysis.rst @@ -3,6 +3,16 @@ Analysis and results Running the analysis and extracting results. +For influence analyses, `shape_function="hermite"` denotes the higher-order +load-distribution path used by *ospgrillage*: the existing Hermite +quadrilateral distributor is used in four-node regions, while three-node skew +regions use the DKT-style condensed triangular distributor. + +For skewed/curved bridges, influence-line station abscissa +(``load_coord="station"``) and influence-surface station axes +(``longitudinal_station``, ``transverse_station``) provide geometry-robust +reduction coordinates independent of global Cartesian layout. + OspGrillage methods ------------------- @@ -10,7 +20,20 @@ OspGrillage methods :toctree: generated/ ~ospgrillage.osp_grillage.OspGrillage.analyze + ~ospgrillage.osp_grillage.OspGrillage.analyze_influence_lines + ~ospgrillage.osp_grillage.OspGrillage.analyze_influence_surfaces + ~ospgrillage.osp_grillage.OspGrillage.analyze_il + ~ospgrillage.osp_grillage.OspGrillage.analyze_is + ~ospgrillage.osp_grillage.OspGrillage.analyze_influence_line + ~ospgrillage.osp_grillage.OspGrillage.analyze_influence_surface + ~ospgrillage.osp_grillage.OspGrillage.get_il + ~ospgrillage.osp_grillage.OspGrillage.get_is ~ospgrillage.osp_grillage.OspGrillage.get_results + ~ospgrillage.osp_grillage.OspGrillage.get_influence_results + ~ospgrillage.osp_grillage.OspGrillage.plot_il + ~ospgrillage.osp_grillage.OspGrillage.plot_is + ~ospgrillage.osp_grillage.OspGrillage.export_il + ~ospgrillage.osp_grillage.OspGrillage.export_is Class reference --------------- @@ -26,3 +49,25 @@ Results .. autoclass:: ospgrillage.osp_grillage.Results :show-inheritance: + +InfluenceResultSet +~~~~~~~~~~~~~~~~~~ + +.. autoclass:: ospgrillage.osp_grillage.InfluenceResultSet + :show-inheritance: + +InfluenceLineResults +~~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: ospgrillage.osp_grillage.InfluenceLineResults + :show-inheritance: + + Includes ``plot()``, ``save()/to_netcdf()``, and ``to_csv()`` helpers. + +InfluenceSurfaceResults +~~~~~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: ospgrillage.osp_grillage.InfluenceSurfaceResults + :show-inheritance: + + Includes ``plot()``, ``save()/to_netcdf()``, and ``to_csv()`` helpers. diff --git a/docs/source/pages/api_postprocessing.rst b/docs/source/pages/api_postprocessing.rst index b0ee2bd2..9b2ed252 100644 --- a/docs/source/pages/api_postprocessing.rst +++ b/docs/source/pages/api_postprocessing.rst @@ -1,16 +1,35 @@ -Post-processing -=============== +Post-processing +=============== + +Influence-line and influence-surface post-processing preserves the +`shape_function` metadata from the stored influence study. When a study was run +with `shape_function="hermite"`, that corresponds to the higher-order load +distribution path: Hermite on quadrilateral regions and the DKT-style condensed +distributor on three-node triangular regions. + +For influence surfaces, :func:`~ospgrillage.postprocessing.plot_is` supports +``coordinate_space="station"`` and ``coordinate_space="physical"``. Physical +space uses mapped deck coordinates and triangulated rendering so curved/skewed +decks are plotted as contiguous surfaces (Matplotlib and Plotly backends). + +For influence lines, :func:`~ospgrillage.postprocessing.plot_il` supports both +``view="ordinate"`` and ``view="path"``. Path view overlays IL ordinates in +3D along the load trajectory using model geometry embedded in the results file. + +Factory functions +----------------- -Factory functions ------------------ - -.. autosummary:: - :toctree: generated/ - - ~ospgrillage.postprocessing.create_envelope - ~ospgrillage.postprocessing.plot_force - ~ospgrillage.postprocessing.plot_bmd - ~ospgrillage.postprocessing.plot_sfd +.. autosummary:: + :toctree: generated/ + + ~ospgrillage.postprocessing.create_envelope + ~ospgrillage.postprocessing.create_influence_line + ~ospgrillage.postprocessing.create_influence_surface + ~ospgrillage.postprocessing.plot_il + ~ospgrillage.postprocessing.plot_is + ~ospgrillage.postprocessing.plot_force + ~ospgrillage.postprocessing.plot_bmd + ~ospgrillage.postprocessing.plot_sfd ~ospgrillage.postprocessing.plot_tmd ~ospgrillage.postprocessing.plot_def ~ospgrillage.postprocessing.plot_model @@ -20,14 +39,26 @@ Factory functions Class reference --------------- -Envelope -~~~~~~~~ - -.. autoclass:: ospgrillage.postprocessing.Envelope - :show-inheritance: - -PostProcessor -~~~~~~~~~~~~~ +Envelope +~~~~~~~~ + +.. autoclass:: ospgrillage.postprocessing.Envelope + :show-inheritance: + +InfluenceLine +~~~~~~~~~~~~~ + +.. autoclass:: ospgrillage.postprocessing.InfluenceLine + :show-inheritance: + +InfluenceSurface +~~~~~~~~~~~~~~~~ + +.. autoclass:: ospgrillage.postprocessing.InfluenceSurface + :show-inheritance: + +PostProcessor +~~~~~~~~~~~~~ .. autoclass:: ospgrillage.postprocessing.PostProcessor :show-inheritance: diff --git a/docs/source/pages/ospgui.md b/docs/source/pages/ospgui.md index b76fd5a8..03ed1264 100644 --- a/docs/source/pages/ospgui.md +++ b/docs/source/pages/ospgui.md @@ -51,13 +51,21 @@ post-processing — no Python scripting required. Use **File > Open Results (.nc)** (or Ctrl+O) to load a self-contained NetCDF file produced by -{meth}`~ospgrillage.osp_grillage.OspGrillage.get_results`: +{meth}`~ospgrillage.osp_grillage.OspGrillage.get_results` or +{meth}`~ospgrillage.osp_grillage.OspGrillage.get_influence_results`: ```python results = bridge.get_results(save_filename="bridge_results.nc") ``` -The interface switches to a results view with five interactive tabs: +Recommended naming is semantic NetCDF: + +- ordinary analysis results: `bridge.res.nc` +- influence-line results: `bridge.il.nc` +- influence-surface results: `bridge.is.nc` + +For ordinary analysis results, the interface switches to a results view with +five interactive tabs: | Tab | Contents | |---|---| @@ -68,7 +76,43 @@ The interface switches to a results view with five interactive tabs: | **Shell Contour** | Shell element contour plot (shell_beam models only) | A left-hand panel provides **loadcase** and **member filter** controls -that apply to all tabs. +that apply to all ordinary-result tabs. + +Influence result files open in a dedicated query mode instead: + +| Dataset type | GUI tab | Key controls | +|---|---|---| +| **Influence line** | **Influence Line** | response array, component, node/element target, load-path axis | +| **Influence surface** | **Influence Surface** | response array, component, node/element target, surface axes, contour/3-D surface view | + +Influence plots are rendered with the Plotly backend in the GUI, while the +same datasets remain available through Python with +{meth}`~ospgrillage.osp_grillage.OspGrillage.get_il`, +{meth}`~ospgrillage.osp_grillage.OspGrillage.get_is`, +{func}`~ospgrillage.postprocessing.plot_il`, and +{func}`~ospgrillage.postprocessing.plot_is`. + +The GUI classifies each loaded dataset and enables only relevant tabs: + +- ordinary results (`*.res.nc`) enable Deflection/BMD/SFD/TMD (and Shell Contour for shell models) +- influence-line results (`*.il.nc`) enable only **Influence Line** +- influence-surface results (`*.is.nc`) enable only **Influence Surface** + +Influence-surface plots in the GUI are always rendered in physical deck +coordinates (`x`,`z`) using triangulated contiguous surfaces, which is +more reliable for skewed and curved geometry. + +Practical interpretation of controls: + +- **IL Axis = `station`** uses cumulative distance along each path (best default + for skewed/curved lane paths). +- **IS X/Y Axis** selects the reduction coordinates (for example station axes); + rendering still appears in mapped physical `x,z` coordinates. + +For influence-line overlays, save a combined influence-line Dataset with an +`InfluenceLine` dimension and open that single file in the GUI. The +`IL Study` selector can then show one stored lane/path study or `All` on the +same plot. ```{image} ../images/gui_results_deflection.png :width: 100% diff --git a/docs/source/pages/tutorials.md b/docs/source/pages/tutorials.md index 15a6fe5c..ccff55e6 100644 --- a/docs/source/pages/tutorials.md +++ b/docs/source/pages/tutorials.md @@ -9,6 +9,7 @@ or viewed as rendered HTML. ../notebooks/super_t_tutorial ../notebooks/advanced_results +../notebooks/influence_lines_and_surfaces ../notebooks/mesh_showcase ../notebooks/multi_span ../notebooks/curve_mesh diff --git a/docs/source/pages/working_with_results.md b/docs/source/pages/working_with_results.md index 1679e9c1..3b3fd5ec 100644 --- a/docs/source/pages/working_with_results.md +++ b/docs/source/pages/working_with_results.md @@ -117,6 +117,97 @@ by_name = force_array.sel(Loadcase="patch load case at global position [0,0,0]" by_index = force_array.isel(Loadcase=0) ``` +## Influence lines and surfaces + +Influence studies are created through +{meth}`~ospgrillage.osp_grillage.OspGrillage.analyze_influence_lines` and +{meth}`~ospgrillage.osp_grillage.OspGrillage.analyze_influence_surfaces`, +which return durable result objects. Those objects carry the compiled +Dataset, can be plotted directly, and can be saved later without recomputing: + +```python +ils = bridge.analyze_influence_lines( + paths={"Lane 1": lane_1_path, "Lane 2": lane_2_path} +) +ils.plot(array="forces", component="Mz_i", element=42) +ils.plot(array="forces", component="Mz_i", element=42, backend="plotly", view="path") +ils.save("lane_ils.nc") + +iss = bridge.analyze_influence_surfaces(x=[2, 4], z=[2, 3], name="Deck IS") +iss.plot(array="forces", component="Mz_i", element=42, view="surface3d") + +# Station-space reduction (recommended for skewed/curved decks) +iss.plot( + array="forces", + component="Mz_i", + element=42, + x_coord="longitudinal_station", + y_coord="transverse_station", +) +``` + +CSV export is available directly from the result objects: + +```python +# IL: one long-form CSV containing all lanes/paths +ils.to_csv( + "lane_ils.csv", + array="forces", + component="Mz_i", + element=42, + load_coord="x", +) + +# IS: rectangular station-grid CSV of ordinates +# (+ optional companion point map with physical x/z) +iss.to_csv( + "deck_is_grid.csv", + array="forces", + component="Mz_i", + element=42, + include_physical_coords=True, +) +``` + +By default, influence surfaces are generated from the mesh station grid +(admissible deck points). For curved or skewed bridges, reduce in station +space with `x_coord="longitudinal_station"` and +`y_coord="transverse_station"` to avoid rectangular-grid assumptions. + +When plotting a reduced surface, `plot_is(..., coordinate_space="physical")` +uses mapped deck coordinates (`x`, `z`) and triangulates valid points so the +surface remains contiguous on curved geometry in both Matplotlib and Plotly. +Use `coordinate_space="station"` to keep station axes. + +The lower-level helpers +{func}`~ospgrillage.postprocessing.create_influence_line`, +{func}`~ospgrillage.postprocessing.create_influence_surface`, +{func}`~ospgrillage.postprocessing.plot_il`, and +{func}`~ospgrillage.postprocessing.plot_is` remain available when direct +`xarray` manipulation is preferred. + +Influence lines support both the traditional `view="ordinate"` plot and a +`view="path"` mode that draws the ordinates along the load path on top of the +bridge model. + +### Technical notes (DKT and station mapping) + +- `shape_function="hermite"` uses the higher-order consistent-load path: + Hermite interpolation on quadrilateral regions, and a DKT-style condensed + triangular distributor where skew geometry creates 3-node regions. +- Influence-line `load_coord="station"` is cumulative distance along each axle + path; this is typically preferred for skewed/curved paths because global `x` + is not a uniform path-abscissa there. +- Influence-surface default generation uses admissible mesh station points + (not a global rectangular deck bounding box), then stores both station + coordinates and mapped physical (`x`,`z`) coordinates. +- Physical plotting (`coordinate_space="physical"`) triangulates valid mapped + points so curved/skewed surfaces remain contiguous without assuming a + rectangular Cartesian deck. + +For a complete workflow, see the +{doc}`../notebooks/influence_lines_and_surfaces` tutorial notebook. + ```{note} For information on the full range of indexing and selection operations available on DataArrays, see the @@ -200,6 +291,11 @@ comb = example_bridge.get_results( ) ``` +If you omit the NetCDF extension and pass a stem (for example +`save_filename="my_results"`), *ospgrillage* writes semantic names: +`*.res.nc` for ordinary results, `*.il.nc` for influence-line exports, and +`*.is.nc` for influence-surface exports. + This writes the full xarray Dataset — including node coordinates and member-element connectivity — to a ``.nc`` file in the current working directory. diff --git a/src/ospgrillage/__init__.py b/src/ospgrillage/__init__.py index cada9832..c01a9959 100644 --- a/src/ospgrillage/__init__.py +++ b/src/ospgrillage/__init__.py @@ -22,6 +22,9 @@ "OspGrillageBeam", "OspGrillageShell", "create_grillage", + "InfluenceResultSet", + "InfluenceLineResults", + "InfluenceSurfaceResults", # Members & sections "GrillageMember", "Section", @@ -59,10 +62,16 @@ "create_point", # Post-processing & plotting "Envelope", + "InfluenceLine", + "InfluenceSurface", "Members", "PostProcessor", "create_envelope", "model_proxy_from_results", + "create_influence_line", + "create_influence_surface", + "plot_il", + "plot_is", "plot_force", "plot_bmd", "plot_sfd", diff --git a/src/ospgrillage/load.py b/src/ospgrillage/load.py index 2a2d7081..4755228f 100644 --- a/src/ospgrillage/load.py +++ b/src/ospgrillage/load.py @@ -1923,5 +1923,4 @@ def _edge_data(xi, zi, xj, zj): ) rhs = np.array([1.0, x - sum(Nmz), z - sum(Nmx)], dtype=float) Nv = np.linalg.solve(coeff_matrix, rhs).tolist() - return Nv, Nmx, Nmz diff --git a/src/ospgrillage/osp_grillage.py b/src/ospgrillage/osp_grillage.py index d401da31..43fbd178 100644 --- a/src/ospgrillage/osp_grillage.py +++ b/src/ospgrillage/osp_grillage.py @@ -5,6 +5,7 @@ This module also handles all load case assignment, analysis, and results by wrapping `OpenSeesPy` command for analysis """ import ast +import csv import dataclasses from dataclasses import dataclass from copy import deepcopy @@ -13,6 +14,7 @@ import json import logging import math +from pathlib import Path as FilePath from typing import List, Tuple, Union, TYPE_CHECKING import numpy as np @@ -25,6 +27,8 @@ LoadVertex, MovingLoad, NodalLoad, + # ospgrillage moving-load path definition, not pathlib.Path + Path, PatchLoading, PointLoad, ShapeFunction, @@ -43,7 +47,11 @@ Envelope, PostProcessor, create_envelope, + create_influence_line, + create_influence_surface, plot_force, + plot_il as plot_influence_line, + plot_is as plot_influence_surface, ) from ospgrillage.utils import ( calculate_area_given_vertices, @@ -73,10 +81,33 @@ "create_grillage", "Analysis", "GrillageElement", + "InfluenceLineResults", + "InfluenceResultSet", + "InfluenceSurfaceResults", "Results", ] +def _normalise_netcdf_filename(filename: str, file_tag: str = None) -> str: + """Return a semantic NetCDF filename. + + Rules: + - If ``filename`` already ends with ``.nc`` (case-insensitive), keep it. + - If it ends with ``.res/.il/.is``, append ``.nc``. + - Otherwise append ``.{file_tag}.nc`` when ``file_tag`` is provided, + else append ``.nc``. + """ + filename = str(filename) + lower = filename.lower() + if lower.endswith(".nc"): + return filename + if lower.endswith((".res", ".il", ".is")): + return f"{filename}.nc" + if file_tag: + return f"{filename}.{file_tag}.nc" + return f"{filename}.nc" + + def _format_ops_cmd(name: str, args: tuple, kwargs: dict) -> str: """Serialise an ops.name(*args, **kwargs) call to Python source code.""" parts = [repr(a) for a in args] @@ -226,6 +257,325 @@ class GrillageElement: a: list +@dataclass +class InfluenceResultSet: + """Separate result container for influence-line or influence-surface analyses.""" + + name: str + kind: str + results: "Results" + load_positions: list + station_positions: list = None + shape_function: str = "linear" + + def compile_data_array(self, local_force_option=True, main_ele_tags=None): + ds = self.results.compile_data_array( + local_force_option=local_force_option, + main_ele_tags=main_ele_tags, + ) + if self.load_positions: + xs = [float(pos[0]) for pos in self.load_positions] + ys = [float(pos[1]) for pos in self.load_positions] + zs = [float(pos[2]) for pos in self.load_positions] + ds = ds.assign_coords( + load_position_x=("Loadcase", xs), + load_position_y=("Loadcase", ys), + load_position_z=("Loadcase", zs), + ) + if self.station_positions: + longitudinal = [float(pos[0]) for pos in self.station_positions] + transverse = [float(pos[1]) for pos in self.station_positions] + ds = ds.assign_coords( + load_position_longitudinal_station=("Loadcase", longitudinal), + load_position_transverse_station=("Loadcase", transverse), + ) + ds.attrs["influence_name"] = self.name + ds.attrs["influence_type"] = self.kind + ds.attrs["shape_function"] = self.shape_function + return ds + + def compile_named_data_array(self, local_force_option=True, main_ele_tags=None): + """Return the compiled Dataset with a leading InfluenceLine/Surface dimension.""" + ds = self.compile_data_array( + local_force_option=local_force_option, + main_ele_tags=main_ele_tags, + ) + index_name = "InfluenceLine" if self.kind == "line" else "InfluenceSurface" + ds = ds.assign_coords( + loadcase_position_index=("Loadcase", np.arange(ds.sizes["Loadcase"], dtype=int)) + ) + ds = ds.assign_coords( + loadcase_label=("Loadcase", [str(label) for label in ds.coords["Loadcase"].values]) + ) + ds = ds.assign_coords(Loadcase=ds.coords["loadcase_position_index"].values) + ds = ds.expand_dims({index_name: [self.name]}) + return ds + + +class _BaseInfluenceResults: + """Durable user-facing container for compiled influence results.""" + + def __init__(self, dataset, *, save_filename=None, file_tag="nc"): + self.dataset = dataset + self.data = dataset + self._file_tag = file_tag + self._save_filename = ( + _normalise_netcdf_filename(save_filename, file_tag=file_tag) + if save_filename + else None + ) + + def save(self, filename: str = None) -> str: + """Save the compiled influence Dataset to NetCDF.""" + target = filename or self._save_filename + if not target: + raise ValueError("filename= is required because no save target was stored") + target = _normalise_netcdf_filename(target, file_tag=self._file_tag) + self.dataset.to_netcdf(target) + self._save_filename = target + return target + + def to_netcdf(self, filename: str = None) -> str: + """Alias for :meth:`save`.""" + return self.save(filename) + + +class InfluenceLineResults(_BaseInfluenceResults): + """User-facing influence-line result container.""" + + def __init__(self, dataset, *, save_filename=None): + super().__init__( + dataset, + save_filename=save_filename, + file_tag="il", + ) + + def get_line(self, **kwargs): + """Return one reduced influence line or a dict of overlaid lines.""" + component = kwargs.get("component", None) + array = kwargs.get("array", None) + if component is None or array is None: + raise ValueError("array= and component= are required to extract an influence line") + + selector_kwargs = dict( + ds=self.dataset, + array=array, + component=component, + node=kwargs.get("node", None), + element=kwargs.get("element", None), + load_coord=kwargs.get("load_coord", "x"), + ) + influence_line = kwargs.get("influence_line", None) + if "InfluenceLine" in self.dataset.dims: + if influence_line is None or influence_line == "All": + return { + str(name): create_influence_line( + influence_line=str(name), + **selector_kwargs, + ).get() + for name in self.dataset.coords["InfluenceLine"].values + } + selector_kwargs["influence_line"] = influence_line + return create_influence_line(**selector_kwargs).get() + + def plot(self, **kwargs): + """Reduce and plot one or more influence lines.""" + plot_kwargs = dict(kwargs) + plot_kwargs.setdefault("dataset", self.dataset) + return plot_influence_line(self.get_line(**kwargs), **plot_kwargs) + + def to_csv(self, filename: str, **kwargs) -> str: + """Export one or more reduced influence lines to a single CSV file. + + The CSV contains long-form rows with columns: + ``influence_line, abscissa, x, y, z, station, ordinate``. + """ + line_data = self.get_line(**kwargs) + if isinstance(line_data, dict): + labelled_lines = list(line_data.items()) + else: + label = kwargs.get("influence_line", None) + if label in (None, "All"): + label = self.dataset.attrs.get("influence_name", "Influence Line") + labelled_lines = [(str(label), line_data)] + + with open(filename, "w", newline="") as csv_file: + writer = csv.writer(csv_file) + writer.writerow( + [ + "influence_line", + "abscissa", + "x", + "y", + "z", + "station", + "ordinate", + ] + ) + for line_name, line in labelled_lines: + axis_name = line.dims[0] + axis_values = np.asarray(line.coords[axis_name].values, dtype=float) + x_values = np.asarray(line.coords["x"].values, dtype=float) + y_values = np.asarray(line.coords["y"].values, dtype=float) + z_values = np.asarray(line.coords["z"].values, dtype=float) + station_values = np.asarray(line.coords["station"].values, dtype=float) + ordinate_values = np.asarray(line.values, dtype=float) + for idx in range(len(ordinate_values)): + writer.writerow( + [ + str(line_name), + float(axis_values[idx]), + float(x_values[idx]), + float(y_values[idx]), + float(z_values[idx]), + float(station_values[idx]), + float(ordinate_values[idx]), + ] + ) + return filename + + +class InfluenceSurfaceResults(_BaseInfluenceResults): + """User-facing influence-surface result container.""" + + def __init__(self, dataset, *, save_filename=None): + super().__init__( + dataset, + save_filename=save_filename, + file_tag="is", + ) + + def get_surface(self, **kwargs): + """Return one reduced influence surface.""" + component = kwargs.get("component", None) + array = kwargs.get("array", None) + if component is None or array is None: + raise ValueError( + "array= and component= are required to extract an influence surface" + ) + return create_influence_surface( + ds=self.dataset, + array=array, + component=component, + node=kwargs.get("node", None), + element=kwargs.get("element", None), + x_coord=kwargs.get("x_coord", "x"), + y_coord=kwargs.get("y_coord", "z"), + influence_surface=kwargs.get("influence_surface", None), + ).get() + + def plot(self, **kwargs): + """Reduce and plot one influence surface.""" + return plot_influence_surface(self.get_surface(**kwargs), **kwargs) + + def _default_surface_axes(self, kwargs): + """Return selector kwargs with station axes preferred when available.""" + selector_kwargs = dict(kwargs) + has_station_coords = ( + "load_position_longitudinal_station" in self.dataset.coords + and "load_position_transverse_station" in self.dataset.coords + ) + if has_station_coords: + selector_kwargs.setdefault("x_coord", "longitudinal_station") + selector_kwargs.setdefault("y_coord", "transverse_station") + return selector_kwargs + + def to_csv(self, filename: str, **kwargs): + """Export a reduced influence surface to CSV. + + Parameters + ---------- + filename : str + Output CSV path. + layout : {"grid", "long"}, default "grid" + ``"grid"`` writes a rectangular station/axis matrix of ordinates. + ``"long"`` writes one row per grid point. + include_physical_coords : bool, default False + If ``True`` and mapped physical coordinates are available, include + physical ``x,z`` coordinates in long-form output, or write a + companion ``*_points.csv`` file for grid output. + """ + layout = kwargs.pop("layout", "grid") + include_physical_coords = kwargs.pop("include_physical_coords", False) + selector_kwargs = self._default_surface_axes(kwargs) + isurface = self.get_surface(**selector_kwargs) + x_name, y_name = isurface.dims + x_values = np.asarray(isurface.coords[x_name].values, dtype=float) + y_values = np.asarray(isurface.coords[y_name].values, dtype=float) + ordinate = np.asarray(isurface.values, dtype=float) + + has_physical = ( + "x" in isurface.coords + and "z" in isurface.coords + and isurface.coords["x"].dims == isurface.dims + and isurface.coords["z"].dims == isurface.dims + ) + + if layout == "long": + with open(filename, "w", newline="") as csv_file: + writer = csv.writer(csv_file) + header = [x_name, y_name] + if include_physical_coords and has_physical: + header += ["x", "z"] + header += ["ordinate"] + writer.writerow(header) + x_phys = ( + np.asarray(isurface.coords["x"].values, dtype=float) + if has_physical + else None + ) + z_phys = ( + np.asarray(isurface.coords["z"].values, dtype=float) + if has_physical + else None + ) + for ix, x_coord in enumerate(x_values): + for iy, y_coord in enumerate(y_values): + row = [float(x_coord), float(y_coord)] + if include_physical_coords and has_physical: + row += [float(x_phys[ix, iy]), float(z_phys[ix, iy])] + row += [float(ordinate[ix, iy])] + writer.writerow(row) + return filename + + if layout != "grid": + raise ValueError("layout must be either 'grid' or 'long'") + + with open(filename, "w", newline="") as csv_file: + writer = csv.writer(csv_file) + writer.writerow([x_name] + [float(val) for val in y_values]) + for ix, x_coord in enumerate(x_values): + row = [float(x_coord)] + [float(val) for val in ordinate[ix, :]] + writer.writerow(row) + + if include_physical_coords and has_physical: + base = FilePath(filename) + suffix = base.suffix or ".csv" + points_file = base.with_name(f"{base.stem}_points{suffix}") + x_phys = np.asarray(isurface.coords["x"].values, dtype=float) + z_phys = np.asarray(isurface.coords["z"].values, dtype=float) + with open(points_file, "w", newline="") as csv_file: + writer = csv.writer(csv_file) + writer.writerow([x_name, y_name, "x", "z", "ordinate"]) + for ix, x_coord in enumerate(x_values): + for iy, y_coord in enumerate(y_values): + writer.writerow( + [ + float(x_coord), + float(y_coord), + float(x_phys[ix, iy]), + float(z_phys[ix, iy]), + float(ordinate[ix, iy]), + ] + ) + return { + "grid": filename, + "points": points_file.as_posix(), + } + + return filename + + class OspGrillage: """ Base class representing an OpenSees grillage structural model. @@ -420,6 +770,7 @@ def __init__( # Mesh objects, pyfile flag, and verbose flag self.pyfile = None self.results = None + self.influence_result_set = dict() self.diagnostics = kwargs.get( "diagnostics", False ) # flag for diagnostics printed to terminal @@ -2341,6 +2692,761 @@ def _embed_geometry(self, ds): return ds + def _run_analysis_case_dicts( + self, + load_case_dict_list: list, + results_obj: "Results", + analysis_type: str = "Static", + step: int = 1, + time_increment: float = 0.01, + ) -> None: + """Evaluate a list of prepared load-case dictionaries into a target Results object.""" + for load_case_dict in load_case_dict_list: + load_case_obj = load_case_dict["loadcase"] + load_command = load_case_dict["load_command"] + load_factor = load_case_dict["load_factor"] + load_case_analysis = Analysis( + analysis_name=load_case_obj.name, + analysis_type=analysis_type, + step=step, + ops_grillage_name=self.model_name, + pyfile=self.pyfile, + time_series_counter=self.global_time_series_counter, + pattern_counter=self.global_pattern_counter, + node_counter=self.Mesh_obj.node_counter, + ele_counter=self.Mesh_obj.element_counter, + constraint_type=self.constraint_type, + load_case=load_case_obj, + time_increment=time_increment, + ) + load_case_analysis.add_load_command(load_command, load_factor=load_factor) + ( + self.global_time_series_counter, + self.global_pattern_counter, + node_disp, + ele_force, + self.analysis_command, + ) = load_case_analysis.evaluate_analysis() + results_obj.extract_analysis(analysis_obj=load_case_analysis) + + def analyze_influence_line(self, **kwargs) -> InfluenceResultSet: + """ + Create and run a dedicated 100 kN axle influence-line analysis. + + This analysis is stored separately from ordinary load-case results and + is retrieved via :func:`get_influence_results`. + + :param name: Influence analysis name. + :type name: str + :param path: Optional explicit Path object describing the axle centre-line movement. + :type path: Path, optional + :param start_point: Start point of the axle path. + :type start_point: Point + :param end_point: End point of the axle path. + :type end_point: Point + :param increments: Number of axle positions along the path. Defaults to ``50``. + :type increments: int, optional + :param step: Target spacing between axle positions along the path. + If provided, overrides ``increments``. + :type step: float, optional + :param axle_load: Axle load magnitude. Defaults to ``100e3``. + :type axle_load: float, optional + """ + name = kwargs.get("name", "Influence Line") + path = kwargs.get("path", None) + start_point = kwargs.get("start_point", None) + end_point = kwargs.get("end_point", None) + increments = kwargs.get("increments", 50) + step_size = kwargs.get("step", None) + axle_load = kwargs.get("axle_load", 100e3) + shape_function = kwargs.get("shape_function", "linear") + analysis_type = kwargs.get("analysis_type", "Static") + analysis_step = kwargs.get("analysis_step", 1) + time_increment = kwargs.get("time_increment", 0.01) + store_results = kwargs.get("store_results", True) + + if path is None and (start_point is None or end_point is None): + raise ValueError("Either path= or both start_point= and end_point= are required for influence-line analysis") + + if path is None: + if step_size is not None: + total_length = get_distance(start_point, end_point) + increments = max(int(np.floor(total_length / step_size)) + 1, 2) + path = Path(start_point=start_point, end_point=end_point, increments=increments) + axle = PointLoad( + name=f"{name} axle", + point1=LoadVertex(0, 0, 0, axle_load), + shape_function=shape_function, + ) + moving_load = MovingLoad(name=name) + moving_load.set_path(path) + moving_load.add_load(axle) + moving_load.parse_moving_load_cases() + + load_case_dict_list = [] + for moving_load_case_list in moving_load.moving_load_case: + for increment_load_case in moving_load_case_list: + load_str = self._distribute_load_types_to_model( + load_case_obj=increment_load_case + ) + load_case_dict_list.append( + { + "name": increment_load_case.name, + "loadcase": increment_load_case, + "load_command": load_str, + "load_factor": 1, + } + ) + + influence_results = Results(self.Mesh_obj) + self._run_analysis_case_dicts( + load_case_dict_list=load_case_dict_list, + results_obj=influence_results, + analysis_type=analysis_type, + step=analysis_step, + time_increment=time_increment, + ) + result_set = InfluenceResultSet( + name=name, + kind="line", + results=influence_results, + load_positions=path.get_path_points(), + shape_function=shape_function, + ) + if store_results: + self.influence_result_set[name] = result_set + return result_set + + return self.get_il( + result_set=result_set, + array=kwargs.get("array", None), + component=kwargs.get("component", None), + node=kwargs.get("node", None), + element=kwargs.get("element", None), + load_coord=kwargs.get("load_coord", "x"), + ) + + def analyze_influence_surface(self, **kwargs) -> InfluenceResultSet: + """ + Create and run a dedicated 100 kN roving-point influence-surface analysis. + + This analysis is stored separately from ordinary load-case results and + is retrieved via :func:`get_influence_results`. + + :param name: Influence analysis name. + :type name: str + :param x: Iterable of global x coordinates for the roving point + (legacy rectangular-grid mode). + :type x: iterable + :param z: Iterable of global z coordinates for the roving point + (legacy rectangular-grid mode). + :type z: iterable + :param longitudinal_stations: Optional iterable of longitudinal mesh + station values (from ``Mesh_obj.nox``) for admissible deck sampling. + :type longitudinal_stations: iterable, optional + :param transverse_stations: Optional iterable of transverse mesh + station values (from ``Mesh_obj.noz``) for admissible deck sampling. + :type transverse_stations: iterable, optional + :param x_step: Default spacing used to generate x coordinates from the deck extent. + :type x_step: float, optional + :param z_step: Default spacing used to generate z coordinates from the deck extent. + :type z_step: float, optional + :param y: Model-plane y coordinate. Defaults to ``0``. + :type y: float, optional + :param point_load: Point-load magnitude. Defaults to ``100e3``. + :type point_load: float, optional + + Notes + ----- + When neither ``x``/``z`` nor station lists are provided, the default + surface is generated from the mesh station grid (admissible deck + points), not from a global rectangular bounding box. + """ + name = kwargs.get("name", "Influence Surface") + x_points = kwargs.get("x", None) + z_points = kwargs.get("z", None) + longitudinal_stations = kwargs.get("longitudinal_stations", None) + transverse_stations = kwargs.get("transverse_stations", None) + x_step = kwargs.get("x_step", 1.0) + z_step = kwargs.get("z_step", 1.0) + y = kwargs.get("y", 0) + point_load = kwargs.get("point_load", 100e3) + shape_function = kwargs.get("shape_function", "linear") + analysis_type = kwargs.get("analysis_type", "Static") + analysis_step = kwargs.get("analysis_step", 1) + time_increment = kwargs.get("time_increment", 0.01) + store_results = kwargs.get("store_results", True) + + station_positions = None + if longitudinal_stations is not None or transverse_stations is not None: + load_positions, station_positions = self._get_surface_station_positions( + longitudinal_stations=longitudinal_stations, + transverse_stations=transverse_stations, + ) + elif x_points is None and z_points is None: + load_positions, station_positions = self._get_surface_station_positions() + else: + if x_points is None or z_points is None: + all_nodes = self.get_nodes() + x_coords = [node_data["coordinate"][0] for node_data in all_nodes.values()] + z_coords = [node_data["coordinate"][2] for node_data in all_nodes.values()] + if x_points is None: + x_points = np.arange(min(x_coords), max(x_coords) + x_step * 0.5, x_step) + if z_points is None: + z_points = np.arange(min(z_coords), max(z_coords) + z_step * 0.5, z_step) + load_positions = [(float(x), float(y), float(z)) for x in x_points for z in z_points] + + load_case_dict_list = [] + for x, y_val, z in load_positions: + point_load_case = LoadCase( + name=f"{name} at global position [{x:.2f},{y_val:.2f},{z:.2f}]" + ) + point = PointLoad( + name=f"{name} point", + point1=LoadVertex(x, y_val, z, point_load), + shape_function=shape_function, + ) + point_load_case.add_load(point) + load_str = self._distribute_load_types_to_model(load_case_obj=point_load_case) + load_case_dict_list.append( + { + "name": point_load_case.name, + "loadcase": point_load_case, + "load_command": load_str, + "load_factor": 1, + } + ) + + influence_results = Results(self.Mesh_obj) + self._run_analysis_case_dicts( + load_case_dict_list=load_case_dict_list, + results_obj=influence_results, + analysis_type=analysis_type, + step=analysis_step, + time_increment=time_increment, + ) + result_set = InfluenceResultSet( + name=name, + kind="surface", + results=influence_results, + load_positions=load_positions, + station_positions=station_positions, + shape_function=shape_function, + ) + if store_results: + self.influence_result_set[name] = result_set + return result_set + + return self.get_is( + result_set=result_set, + array=kwargs.get("array", None), + component=kwargs.get("component", None), + node=kwargs.get("node", None), + element=kwargs.get("element", None), + x_coord=kwargs.get("x_coord", "x"), + y_coord=kwargs.get("y_coord", "z"), + ) + + def analyze_influence_lines(self, **kwargs) -> InfluenceLineResults: + """Run one or more influence-line studies and return a result object.""" + save_filename = kwargs.get("save_filename", None) + local_force_flag = kwargs.get("local_forces", False) + + paths = kwargs.get("paths", None) + if paths is None: + single_name = kwargs.get("name", "Influence Line") + result_set = self.analyze_influence_line(**kwargs) + ds = self._compile_influence_dataset( + result_set, + local_force_flag=local_force_flag, + ) + else: + if not isinstance(paths, dict) or not paths: + raise ValueError("paths= must be a non-empty dict of {name: Path}") + common_kwargs = dict(kwargs) + common_kwargs.pop("paths", None) + common_kwargs.pop("path", None) + common_kwargs.pop("name", None) + common_kwargs.pop("save_filename", None) + common_kwargs.pop("local_forces", None) + line_names = [] + for line_name, path in paths.items(): + self.analyze_influence_line( + name=str(line_name), + path=path, + **common_kwargs, + ) + line_names.append(str(line_name)) + ds = self._compile_combined_influence_line_dataset( + line_names, + local_force_flag=local_force_flag, + ) + + result = InfluenceLineResults(ds, save_filename=save_filename) + if save_filename: + result.save() + return result + + def analyze_influence_surfaces(self, **kwargs) -> InfluenceSurfaceResults: + """Run an influence-surface study and return a result object.""" + save_filename = kwargs.get("save_filename", None) + local_force_flag = kwargs.get("local_forces", False) + result_set = self.analyze_influence_surface(**kwargs) + ds = self._compile_influence_dataset( + result_set, + local_force_flag=local_force_flag, + ) + result = InfluenceSurfaceResults(ds, save_filename=save_filename) + if save_filename: + result.save() + return result + + def _compile_influence_dataset(self, result_set: InfluenceResultSet, local_force_flag=False): + ds = result_set.compile_data_array( + local_force_option=local_force_flag, + main_ele_tags=self.Mesh_obj.element_counter, + ) + return self._embed_geometry(ds) + + def _get_surface_station_positions( + self, + longitudinal_stations=None, + transverse_stations=None, + ): + """Map logical deck stations to admissible physical load positions.""" + x_station_values = [float(val) for val in np.asarray(self.Mesh_obj.nox).tolist()] + z_station_values = [float(val) for val in np.asarray(self.Mesh_obj.noz).tolist()] + + def _select_indices(requested, available): + if requested is None: + return list(range(len(available))) + indices = [] + for value in requested: + matches = [ + idx for idx, candidate in enumerate(available) + if np.isclose(candidate, float(value)) + ] + if not matches: + raise ValueError( + f"Requested station {value!r} is not available in the mesh station grid" + ) + indices.append(matches[0]) + return indices + + x_indices = _select_indices(longitudinal_stations, x_station_values) + z_indices = _select_indices(transverse_stations, z_station_values) + + y_ref = getattr(self.Mesh_obj, "y_elevation", 0.0) + station_to_coord = {} + for node_data in self.Mesh_obj.node_spec.values(): + x_group = node_data.get("x_group") + z_group = node_data.get("z_group") + if not isinstance(x_group, (int, np.integer)) or not isinstance(z_group, (int, np.integer)): + continue + coord = node_data["coordinate"] + score = abs(coord[1] - y_ref) + key = (int(x_group), int(z_group)) + current = station_to_coord.get(key) + if current is None or score < current[1]: + station_to_coord[key] = (coord, score) + + if not station_to_coord: + raise ValueError("No mesh station groups were found for influence-surface sampling") + + # Build longitudinal adjacency by row (constant z-group) from beam connectivity. + row_adjacency = {} + long_like_elements = list(getattr(self.Mesh_obj, "long_ele", [])) + list( + getattr(self.Mesh_obj, "connect_ele", []) + ) + for element in long_like_elements: + if len(element) < 3: + continue + node_i = self.Mesh_obj.node_spec.get(element[1], {}) + node_j = self.Mesh_obj.node_spec.get(element[2], {}) + x_group_i = node_i.get("x_group") + x_group_j = node_j.get("x_group") + z_group_i = node_i.get("z_group") + z_group_j = node_j.get("z_group") + if ( + not isinstance(x_group_i, (int, np.integer)) + or not isinstance(x_group_j, (int, np.integer)) + or not isinstance(z_group_i, (int, np.integer)) + or not isinstance(z_group_j, (int, np.integer)) + ): + continue + x_group_i = int(x_group_i) + x_group_j = int(x_group_j) + z_group_i = int(z_group_i) + z_group_j = int(z_group_j) + if z_group_i != z_group_j or x_group_i == x_group_j: + continue + row_adj = row_adjacency.setdefault(z_group_i, {}) + row_adj.setdefault(x_group_i, set()).add(x_group_j) + row_adj.setdefault(x_group_j, set()).add(x_group_i) + + edge_node_recorder = getattr(self.Mesh_obj, "edge_node_recorder", None) + first_edge = ( + min(edge_node_recorder.values()) + if isinstance(edge_node_recorder, dict) and edge_node_recorder + else None + ) + + def _row_start_group(z_group, row_groups, adjacency, x_coords): + if first_edge is not None: + counts = {} + for node_tag, edge_id in edge_node_recorder.items(): + if edge_id != first_edge: + continue + node_data = self.Mesh_obj.node_spec.get(node_tag, {}) + node_z_group = node_data.get("z_group") + node_x_group = node_data.get("x_group") + if ( + not isinstance(node_z_group, (int, np.integer)) + or not isinstance(node_x_group, (int, np.integer)) + ): + continue + if int(node_z_group) != z_group: + continue + node_x_group = int(node_x_group) + if node_x_group not in row_groups: + continue + counts[node_x_group] = counts.get(node_x_group, 0) + 1 + if counts: + return max(counts, key=counts.get) + + end_groups = [group for group in row_groups if len(adjacency.get(group, set())) <= 1] + if end_groups: + return min(end_groups, key=lambda group: x_coords[group]) + return min(row_groups, key=lambda group: x_coords[group]) + + def _order_row_groups(z_group): + row_groups = sorted( + { + group_x + for (group_x, group_z) in station_to_coord + if group_z == z_group + } + ) + if not row_groups: + return [] + x_coords = { + group_x: float(station_to_coord[(group_x, z_group)][0][0]) + for group_x in row_groups + } + adjacency = {group: set() for group in row_groups} + for group, neighbours in row_adjacency.get(z_group, {}).items(): + if group not in adjacency: + continue + adjacency[group].update([n for n in neighbours if n in adjacency]) + + start_group = _row_start_group(z_group, row_groups, adjacency, x_coords) + ordered = [] + visited = set() + previous = None + current = start_group + while current is not None and current not in visited: + ordered.append(current) + visited.add(current) + neighbours = [ + group + for group in adjacency.get(current, set()) + if group != previous and group not in visited + ] + if not neighbours: + break + if len(neighbours) == 1: + next_group = neighbours[0] + else: + next_group = min(neighbours, key=lambda group: x_coords[group]) + previous = current + current = next_group + + if len(ordered) < len(row_groups): + remaining = [group for group in row_groups if group not in visited] + remaining.sort(key=lambda group: x_coords[group]) + ordered.extend(remaining) + return ordered + + row_data = {} + for z_group in sorted({group_z for _, group_z in station_to_coord}): + ordered_groups = _order_row_groups(z_group) + if not ordered_groups: + continue + points = np.asarray( + [station_to_coord[(group_x, z_group)][0] for group_x in ordered_groups], + dtype=float, + ) + if len(points) == 1: + cumulative = np.asarray([0.0], dtype=float) + else: + segment_lengths = np.linalg.norm( + np.diff(points[:, [0, 2]], axis=0), + axis=1, + ) + cumulative = np.concatenate([[0.0], np.cumsum(segment_lengths)]) + row_data[z_group] = { + "points": points, + "cumulative": cumulative, + "z_ref": float(np.mean(points[:, 2])), + } + + requested_z_indices = sorted(set(z_indices)) + if len(row_data) < len(requested_z_indices): + raise ValueError("No admissible surface load positions were found for the requested station grid") + + if all(z_idx in row_data for z_idx in requested_z_indices): + # In standard meshes the node z_group is the transverse station index. + # Keep that direct mapping to avoid cyclic row reorder on curved decks. + z_index_to_group = {z_idx: z_idx for z_idx in requested_z_indices} + else: + z_index_to_group = {} + remaining_groups = set(row_data.keys()) + for z_idx in requested_z_indices: + target_z = z_station_values[z_idx] + chosen_group = min( + remaining_groups, + key=lambda group: abs(row_data[group]["z_ref"] - target_z), + ) + z_index_to_group[z_idx] = chosen_group + remaining_groups.remove(chosen_group) + + x_station_start = x_station_values[0] + x_station_end = x_station_values[-1] + x_station_span = x_station_end - x_station_start + + def _station_ratio(station): + if np.isclose(x_station_span, 0.0): + return 0.0 + return (station - x_station_start) / x_station_span + + def _interpolate_row(row_points, cumulative, station): + ratio = float(np.clip(_station_ratio(station), 0.0, 1.0)) + if len(row_points) == 1: + coord = row_points[0] + return float(coord[0]), float(coord[1]), float(coord[2]) + total = cumulative[-1] + if np.isclose(total, 0.0): + coord = row_points[0] + return float(coord[0]), float(coord[1]), float(coord[2]) + target = ratio * total + segment_index = int(np.searchsorted(cumulative, target, side="right") - 1) + segment_index = max(0, min(segment_index, len(cumulative) - 2)) + segment_start = cumulative[segment_index] + segment_end = cumulative[segment_index + 1] + denom = segment_end - segment_start + local_ratio = 0.0 if np.isclose(denom, 0.0) else (target - segment_start) / denom + p0 = row_points[segment_index] + p1 = row_points[segment_index + 1] + interp = p0 + local_ratio * (p1 - p0) + return float(interp[0]), float(interp[1]), float(interp[2]) + + load_positions = [] + station_positions = [] + for x_idx in x_indices: + for z_idx in z_indices: + z_group = z_index_to_group[z_idx] + row_points = row_data[z_group]["points"] + cumulative = row_data[z_group]["cumulative"] + if len(row_points) == len(x_station_values): + # Use exact mesh row points when row station count matches. + point = row_points[x_idx] + coord = (float(point[0]), float(point[1]), float(point[2])) + else: + coord = _interpolate_row(row_points, cumulative, x_station_values[x_idx]) + load_positions.append(coord) + station_positions.append((x_station_values[x_idx], z_station_values[z_idx])) + + if not load_positions: + raise ValueError("No admissible surface load positions were found for the requested station grid") + + return load_positions, station_positions + + def _compile_combined_influence_line_dataset(self, names, local_force_flag=False): + """Combine multiple stored influence-line result sets into one Dataset.""" + if not names: + raise ValueError("At least one influence-line name is required") + + datasets = [] + expected_kind = "line" + for line_name in names: + if line_name not in self.influence_result_set: + raise ValueError(f"No influence analysis named {line_name!r} has been run") + result_set = self.influence_result_set[line_name] + if result_set.kind != expected_kind: + raise ValueError("Combined influence-line export only supports line studies") + ds = result_set.compile_named_data_array( + local_force_option=local_force_flag, + main_ele_tags=self.Mesh_obj.element_counter, + ) + ds = self._embed_geometry(ds) + datasets.append(ds) + + base = datasets[0] + for ds in datasets[1:]: + if ds.attrs.get("model_type") != base.attrs.get("model_type"): + raise ValueError("All combined influence-line studies must come from the same model type") + if not ds["node_coordinates"].equals(base["node_coordinates"]): + raise ValueError("All combined influence-line studies must share the same node layout") + + combined = xr.concat( + datasets, + dim="InfluenceLine", + join="outer", + compat="override", + coords="all", + data_vars="all", + ) + combined.attrs["influence_name"] = ", ".join(names) + combined.attrs["influence_type"] = "line" + combined.attrs["influence_overlay"] = "multi" + combined.attrs["shape_function"] = ",".join( + [self.influence_result_set[line_name].shape_function for line_name in names] + ) + return combined + + def get_influence_results(self, name: str = None, **kwargs): + """ + Return a stored influence-line or influence-surface result Dataset. + + :param name: Influence analysis name. + :type name: str + :param names: Optional list of influence-line analysis names to combine + into one Dataset for overlay plotting. + :type names: list[str], optional + :param save_filename: Optional NetCDF output path. + :type save_filename: str, optional + :param local_forces: If ``True``, use local element forces. Defaults to ``False``. + :type local_forces: bool, optional + """ + local_force_flag = kwargs.get("local_forces", False) + save_filename = kwargs.get("save_filename", None) + + names = kwargs.get("names", None) + if names is not None: + ds = self._compile_combined_influence_line_dataset( + list(names), + local_force_flag=local_force_flag, + ) + else: + if name is None: + raise ValueError("name= or names= is required") + if name not in self.influence_result_set: + raise ValueError(f"No influence analysis named {name!r} has been run") + ds = self._compile_influence_dataset( + self.influence_result_set[name], + local_force_flag=local_force_flag, + ) + if save_filename: + save_target = _normalise_netcdf_filename( + save_filename, + file_tag="il" if ds.attrs.get("influence_type") == "line" else "is", + ) + ds.to_netcdf(save_target) + return ds + + def clear_influence_results(self, name: str = None) -> None: + """Clear one named influence analysis or all stored influence analyses.""" + if name is None: + self.influence_result_set = dict() + else: + self.influence_result_set.pop(name, None) + + def analyze_il(self, **kwargs): + """Short alias for :func:`analyze_influence_line`.""" + return self.analyze_influence_line(**kwargs) + + def analyze_is(self, **kwargs): + """Short alias for :func:`analyze_influence_surface`.""" + return self.analyze_influence_surface(**kwargs) + + def get_il(self, name: str = None, result_set: InfluenceResultSet = None, **kwargs): + """Extract a reduced influence line using xarray-style ``array`` and ``component`` selectors.""" + if result_set is None: + if name is None: + raise ValueError("name= or result_set= is required") + ds = self.get_influence_results(name, local_forces=kwargs.get("local_forces", False)) + else: + ds = result_set.compile_data_array( + local_force_option=kwargs.get("local_forces", False), + main_ele_tags=self.Mesh_obj.element_counter, + ) + + component = kwargs.get("component", None) + array = kwargs.get("array", None) + if component is None or array is None: + raise ValueError("array= and component= are required to extract an influence line") + return create_influence_line( + ds=ds, + array=array, + component=component, + node=kwargs.get("node", None), + element=kwargs.get("element", None), + load_coord=kwargs.get("load_coord", "x"), + ).get() + + def get_is(self, name: str = None, result_set: InfluenceResultSet = None, **kwargs): + """Extract a reduced influence surface using xarray-style ``array`` and ``component`` selectors.""" + if result_set is None: + if name is None: + raise ValueError("name= or result_set= is required") + ds = self.get_influence_results(name, local_forces=kwargs.get("local_forces", False)) + else: + ds = result_set.compile_data_array( + local_force_option=kwargs.get("local_forces", False), + main_ele_tags=self.Mesh_obj.element_counter, + ) + + component = kwargs.get("component", None) + array = kwargs.get("array", None) + if component is None or array is None: + raise ValueError("array= and component= are required to extract an influence surface") + return create_influence_surface( + ds=ds, + array=array, + component=component, + node=kwargs.get("node", None), + element=kwargs.get("element", None), + x_coord=kwargs.get("x_coord", "x"), + y_coord=kwargs.get("y_coord", "z"), + ).get() + + def plot_il(self, il, **kwargs): + """Plot a reduced influence line.""" + return plot_influence_line(il, **kwargs) + + def plot_is(self, isurface, **kwargs): + """Plot a reduced influence surface as a filled contour.""" + return plot_influence_surface(isurface, **kwargs) + + def export_il(self, il, filename: str) -> None: + """Export a reduced influence line to CSV.""" + axis_name = il.dims[0] + data = np.column_stack( + [ + il.coords["x"].values, + il.coords["y"].values, + il.coords["z"].values, + il.values, + ] + ) + header = "x,y,z,ordinate" + np.savetxt(filename, data, delimiter=",", header=header, comments="") + + def export_is(self, isurface, filename: str) -> None: + """Export a reduced influence surface to CSV in long-form.""" + x_name, y_name = isurface.dims[0], isurface.dims[1] + rows = [] + for x in isurface.coords[x_name].values: + for y in isurface.coords[y_name].values: + rows.append([x, y, isurface.sel({x_name: x, y_name: y}).item()]) + np.savetxt( + filename, + np.asarray(rows, dtype=float), + delimiter=",", + header=f"{x_name},{y_name},ordinate", + comments="", + ) + def get_results(self, **kwargs): """ Return analysis results as an xarray ``Dataset``. @@ -2381,6 +3487,11 @@ def get_results(self, **kwargs): # get kwargs comb = kwargs.get("combinations", False) # if Boolean true save_filename = kwargs.get("save_filename", None) # str of file name + save_target = ( + _normalise_netcdf_filename(save_filename, file_tag="res") + if save_filename + else None + ) specific_load_case = kwargs.get("load_case", None) # str of fil local_force_flag = kwargs.get("local_forces", False) basic_da = self.results.compile_data_array( @@ -2507,14 +3618,14 @@ def get_results(self, **kwargs): combination_array = factored_array.assign_coords( Loadcase=coordinate_name_list ) - if save_filename: - combination_array.to_netcdf(save_filename) + if save_target: + combination_array.to_netcdf(save_target) return combination_array else: # return raw data array for manual post processing - if save_filename: - basic_da.to_netcdf(save_filename) + if save_target: + basic_da.to_netcdf(save_target) return basic_da def get_element(self, **kwargs) -> Union[List[float]]: diff --git a/src/ospgrillage/ospgui.py b/src/ospgrillage/ospgui.py index ba2bec10..eb89a945 100644 --- a/src/ospgrillage/ospgui.py +++ b/src/ospgrillage/ospgui.py @@ -1,6 +1,7 @@ import sys import os import logging +import numpy as np logger = logging.getLogger(__name__) @@ -59,6 +60,88 @@ class QApplication: pass # type: ignore[assignment] +def _classify_results_kind(ds) -> str: + """Classify a loaded results Dataset for GUI tab enablement. + + Returns one of ``"ordinary"``, ``"influence_line"``, or + ``"influence_surface"``. + """ + influence_type = str(ds.attrs.get("influence_type", "")).strip().lower() + if influence_type == "line": + return "influence_line" + if influence_type == "surface": + return "influence_surface" + + # Combined influence exports encode study dimensions explicitly. + if "InfluenceLine" in ds.dims or "InfluenceLine" in ds.coords: + return "influence_line" + if "InfluenceSurface" in ds.dims or "InfluenceSurface" in ds.coords: + return "influence_surface" + + coord_names = set(ds.coords) + if { + "load_position_longitudinal_station", + "load_position_transverse_station", + }.issubset(coord_names): + return "influence_surface" + + # Legacy influence exports may only carry influence_name + load positions. + if "influence_name" in ds.attrs and {"load_position_x", "load_position_z"}.issubset(coord_names): + try: + x_vals = [float(v) for v in ds.coords["load_position_x"].values.tolist()] + z_vals = [float(v) for v in ds.coords["load_position_z"].values.tolist()] + ux = len(set(x_vals)) + uz = len(set(z_vals)) + if ux > 1 and uz > 1 and ux * uz <= len(x_vals): + return "influence_surface" + except Exception: + pass + return "influence_line" + + return "ordinary" + + +def _resolve_influence_surface_axes(ds, mode: str = "physical"): + """Resolve influence-surface reduction axes and display coordinate-space. + + Reduction is performed on station coordinates when available because that + produces a contiguous 2D grid for skewed/curved decks. Display space can + still be physical x-z or station, depending on user mode. + """ + mode_value = str(mode or "physical").strip().lower() + has_station_coords = { + "load_position_longitudinal_station", + "load_position_transverse_station", + }.issubset(set(ds.coords)) + + if has_station_coords: + if mode_value in {"station", "stations"}: + return "longitudinal_station", "transverse_station", "station" + return "longitudinal_station", "transverse_station", "physical" + + # Fallback for legacy/non-station datasets. + return "x", "z", "physical" + + +def _nonempty_components(da): + """Return Component labels that contain at least one finite/non-null value.""" + if "Component" not in da.dims: + return [] + mask = da.notnull() + reduce_dims = [dim for dim in da.dims if dim != "Component"] + if reduce_dims: + mask = mask.any(dim=reduce_dims) + valid = [] + for comp in da.coords["Component"].values: + try: + has_data = bool(mask.sel(Component=comp).item()) + except Exception: + has_data = False + if has_data: + valid.append(str(comp)) + return valid + + class BridgeInputWidget(QWidget): """Tabbed input form for bridge geometry, materials, sections, and members. @@ -769,6 +852,73 @@ def __init__(self): self.contour_group.setVisible(False) # shown only for shell_beam layout.addWidget(self.contour_group) + # --- Influence controls --- + self.influence_group = QGroupBox("Influence Query") + influence_layout = QVBoxLayout() + + self.influence_common_group = QGroupBox("Target Response") + influence_common_layout = QFormLayout() + self.influence_array_combo = QComboBox() + self.influence_component_combo = QComboBox() + self._influence_target_type = None + self.influence_target_id_combo = QComboBox() + influence_common_layout.addRow("Response:", self.influence_array_combo) + influence_common_layout.addRow("Component:", self.influence_component_combo) + self.influence_target_id_combo.setToolTip( + "Identifier of the selected response location." + ) + self.influence_target_id_label = QLabel("Location ID") + influence_common_layout.addRow(self.influence_target_id_label, self.influence_target_id_combo) + self.influence_common_group.setLayout(influence_common_layout) + influence_layout.addWidget(self.influence_common_group) + + self.influence_line_group = QGroupBox("Influence Line") + influence_line_layout = QFormLayout() + self.influence_line_set_combo = QComboBox() + self.influence_line_view_combo = QComboBox() + for mode in ("ordinate", "path"): + self.influence_line_view_combo.addItem(mode) + self.influence_line_view_combo.setCurrentText("path") + self.influence_line_axis_combo = QComboBox() + for axis in ("station", "x", "y", "z"): + self.influence_line_axis_combo.addItem(axis) + self.influence_line_axis_combo.setToolTip( + "Influence-line abscissa. Use 'station' for cumulative path distance." + ) + influence_line_layout.addRow("IL Study:", self.influence_line_set_combo) + influence_line_layout.addRow("IL View:", self.influence_line_view_combo) + influence_line_layout.addRow("IL Abscissa:", self.influence_line_axis_combo) + self.influence_line_group.setLayout(influence_line_layout) + influence_layout.addWidget(self.influence_line_group) + + self.influence_surface_group = QGroupBox("Influence Surface") + influence_surface_layout = QFormLayout() + self.influence_surface_mode_combo = QComboBox() + self.influence_surface_mode_combo.addItem("Physical (x-z)", "physical") + self.influence_surface_mode_combo.addItem( + "Stations (longitudinal-transverse)", "stations" + ) + self.influence_surface_mode_combo.setToolTip( + "Single IS coordinate mode. Physical uses x-z; stations uses longitudinal-transverse station grid." + ) + self.influence_surface_view_combo = QComboBox() + for mode in ("contour", "surface3d"): + self.influence_surface_view_combo.addItem(mode) + self.influence_surface_show_layer_check = QCheckBox("Show IS Layer") + self.influence_surface_show_layer_check.setChecked(True) + self.influence_surface_hover_check = QCheckBox("Enable IS Hover") + self.influence_surface_hover_check.setChecked(False) + + influence_surface_layout.addRow("IS Coordinates:", self.influence_surface_mode_combo) + influence_surface_layout.addRow("IS View:", self.influence_surface_view_combo) + influence_surface_layout.addRow("", self.influence_surface_show_layer_check) + influence_surface_layout.addRow("", self.influence_surface_hover_check) + self.influence_surface_group.setLayout(influence_surface_layout) + influence_layout.addWidget(self.influence_surface_group) + self.influence_group.setLayout(influence_layout) + self.influence_group.setVisible(False) + layout.addWidget(self.influence_group) + # --- Back button --- self.btn_back = QPushButton("Back to Wizard") layout.addWidget(self.btn_back) @@ -786,16 +936,53 @@ def populate_loadcases(self, loadcase_names): def update_available_members(self, proxy): """Enable/disable checkboxes based on which members have elements.""" + def _norm_member_key(name): + text = str(name or "") + return "".join(ch for ch in text.lower() if ch.isalnum()) + + def _has_elements(info): + if isinstance(info, dict): + values = info.get("elements", None) + if values is None: + values = info.get("element", None) + else: + values = info + if values is None: + return False + if isinstance(values, np.ndarray): + return values.size > 0 + if isinstance(values, (list, tuple, set)): + if len(values) == 0: + return False + for item in values: + if isinstance(item, np.ndarray): + if item.size > 0: + return True + elif isinstance(item, (list, tuple, set)): + if len(item) > 0: + return True + elif item is not None: + return True + return False + return bool(values) + + raw_members = getattr(proxy, "_members", {}) or {} + members_by_norm = {_norm_member_key(name): info for name, info in raw_members.items()} + has_member_catalog = bool(members_by_norm) + if not has_member_catalog: + for cb in self.member_checkboxes.values(): + cb.setEnabled(True) + cb.setChecked(True) + return + for flag_name, cb in self.member_checkboxes.items(): - member_name = flag_name.lower() - has_elements = False - info = proxy._members.get(member_name, {}) - for group in info.get("elements", []): - if group: - has_elements = True - break - cb.setEnabled(has_elements) - cb.setChecked(has_elements) + info = members_by_norm.get(_norm_member_key(flag_name), None) + has_elements = _has_elements(info) + # Keep filters interactive when the file includes member metadata + # but uses naming that does not exactly match our enum labels. + unknown_mapping = has_member_catalog and info is None + cb.setEnabled(has_elements or unknown_mapping) + cb.setChecked(has_elements or unknown_mapping) def set_file_info(self, filename, summary): """Update the file info label.""" @@ -835,6 +1022,214 @@ def set_shell_contour_enabled(self, enabled): "" if enabled else "QGroupBox { color: #999; }" ) + def set_influence_mode(self, enabled, result_kind=None): + """Show/hide influence-query controls and relevant IL/IS subsets.""" + self.influence_group.setVisible(enabled) + self.loadcase_combo.setEnabled(not enabled) + show_line = enabled and result_kind == "influence_line" + show_surface = enabled and result_kind == "influence_surface" + if enabled and result_kind not in {"influence_line", "influence_surface"}: + show_line = True + show_surface = True + self.influence_line_group.setVisible(show_line) + self.influence_surface_group.setVisible(show_surface) + + def populate_influence_controls(self, ds): + """Populate array/component/target selectors for an influence Dataset.""" + arrays = [] + excluded_arrays = {"velocity", "velocities", "acceleration", "accelerations"} + for name, da in ds.data_vars.items(): + if str(name).lower() in excluded_arrays: + continue + if "Loadcase" in da.dims and "Component" in da.dims: + if "Node" in da.dims or "Element" in da.dims: + if _nonempty_components(da): + arrays.append(name) + preferred_order = {"displacements": 0, "forces": 1} + arrays.sort(key=lambda array_name: (preferred_order.get(array_name, 99), str(array_name))) + + self.influence_array_combo.blockSignals(True) + self.influence_array_combo.clear() + for name in arrays: + self.influence_array_combo.addItem(name) + self.influence_array_combo.blockSignals(False) + + self.influence_line_set_combo.blockSignals(True) + self.influence_line_set_combo.clear() + if "InfluenceLine" in ds.coords: + self.influence_line_set_combo.addItem("All") + for name in ds.coords["InfluenceLine"].values: + self.influence_line_set_combo.addItem(str(name)) + self.influence_line_set_combo.setEnabled(True) + else: + self.influence_line_set_combo.addItem("Current") + self.influence_line_set_combo.setEnabled(False) + self.influence_line_set_combo.blockSignals(False) + + has_station_coords = { + "load_position_longitudinal_station", + "load_position_transverse_station", + }.issubset(set(ds.coords)) + self.influence_surface_mode_combo.blockSignals(True) + self.influence_surface_mode_combo.clear() + self.influence_surface_mode_combo.addItem("Physical (x-z)", "physical") + if has_station_coords: + self.influence_surface_mode_combo.addItem( + "Stations (longitudinal-transverse)", "stations" + ) + self.influence_surface_mode_combo.blockSignals(False) + + if arrays: + self._update_influence_component_and_target(ds, arrays[0]) + else: + self.influence_component_combo.blockSignals(True) + self.influence_component_combo.clear() + self.influence_component_combo.blockSignals(False) + + def _update_influence_component_and_target(self, ds, array_name): + """Refresh component and target-id selectors from the selected array.""" + da = ds[array_name] + + target_types = [] + if "Node" in da.dims: + target_types.append("Node") + if "Element" in da.dims: + target_types.append("Element") + + if not target_types: + self._influence_target_type = None + self.influence_target_id_label.setText("Location ID") + self.influence_target_id_combo.blockSignals(True) + self.influence_target_id_combo.clear() + self.influence_target_id_combo.blockSignals(False) + self.influence_target_id_combo.setEnabled(False) + self.influence_target_id_combo.setToolTip( + f"No valid response location IDs are available for response '{array_name}'." + ) + self.influence_component_combo.blockSignals(True) + self.influence_component_combo.clear() + self.influence_component_combo.blockSignals(False) + return + + preferred = self._influence_target_type + if preferred not in target_types: + preferred = "Node" if "Node" in target_types else target_types[0] + self._influence_target_type = preferred + self.influence_target_id_label.setText(f"{preferred} ID") + self.influence_target_id_combo.setEnabled(True) + self.influence_target_id_combo.setToolTip( + f"Identifier of the selected {preferred.lower()} response location." + ) + + self._update_influence_target_ids(ds, array_name) + self._update_influence_components(ds, array_name) + + def _update_influence_target_ids(self, ds, array_name): + """Refresh target IDs when the target type changes.""" + da = ds[array_name] + target_type = self._influence_target_type + if target_type is None: + self.influence_target_id_combo.blockSignals(True) + self.influence_target_id_combo.clear() + self.influence_target_id_combo.blockSignals(False) + return + coord_name = "Node" if target_type == "Node" else "Element" + target_ids = da.coords[coord_name].values if coord_name in da.coords else [] + + self.influence_target_id_combo.blockSignals(True) + self.influence_target_id_combo.clear() + for target_id in target_ids: + self.influence_target_id_combo.addItem(str(int(target_id))) + self.influence_target_id_combo.blockSignals(False) + + def _update_influence_components(self, ds, array_name): + """Refresh component options for selected response array/location.""" + da = ds[array_name] + target_type = self._influence_target_type + target_id = self.influence_target_id_combo.currentText() + + da_for_components = da + if target_type == "Node" and "Node" in da.dims and target_id: + try: + da_for_components = da.sel(Node=int(target_id)) + except Exception: + da_for_components = da + elif target_type == "Element" and "Element" in da.dims and target_id: + try: + da_for_components = da.sel(Element=int(target_id)) + except Exception: + da_for_components = da + + components = _nonempty_components(da_for_components) + if not components: + components = _nonempty_components(da) + + self.influence_component_combo.blockSignals(True) + self.influence_component_combo.clear() + for comp in components: + self.influence_component_combo.addItem(str(comp)) + self.influence_component_combo.blockSignals(False) + + def selected_influence_array(self): + return self.influence_array_combo.currentText() or None + + def selected_influence_component(self): + component = self.influence_component_combo.currentText() + if not component: + return None + return component + + def selected_influence_target(self): + target_type = self._influence_target_type + target_id = self.influence_target_id_combo.currentText() + if not target_type or not target_id: + return None, None + return target_type, int(target_id) + + def selected_influence_line_axis(self): + return self.influence_line_axis_combo.currentText() or "x" + + def selected_influence_line_study(self): + text = self.influence_line_set_combo.currentText() + if not text or text == "Current": + return None + if text == "All": + return "All" + return text + + def selected_influence_line_view(self): + return self.influence_line_view_combo.currentText() or "path" + + def set_influence_line_axis_enabled(self, enabled: bool): + """Enable IL abscissa selection only when it affects the plot.""" + self.influence_line_axis_combo.setEnabled(enabled) + if enabled: + self.influence_line_axis_combo.setToolTip( + "Influence-line abscissa. Use 'station' for cumulative path distance." + ) + else: + self.influence_line_axis_combo.setToolTip( + "IL view='path' always uses station order along the load path." + ) + + def selected_influence_surface_mode(self): + mode = self.influence_surface_mode_combo.currentData() + if mode: + return str(mode) + text = self.influence_surface_mode_combo.currentText().lower() + if "station" in text: + return "stations" + return "physical" + + def selected_influence_surface_view(self): + return self.influence_surface_view_combo.currentText() or "contour" + + def is_influence_surface_layer_visible(self): + return self.influence_surface_show_layer_check.isChecked() + + def is_influence_surface_hover_enabled(self): + return self.influence_surface_hover_check.isChecked() + class BridgeAnalysisGUI(QMainWindow): """Main window for the *ospgui* bridge geometry generator. @@ -976,6 +1371,7 @@ def __init__(self): self._mode = "wizard" # "wizard" or "results" self._model_proxy = None # _ModelProxy from loaded results self._results = None # xarray Dataset + self._results_kind = "ordinary" # ordinary, influence_line, influence_surface self._stale_tabs = set() # result tab names needing re-render # Create UI components @@ -1008,7 +1404,15 @@ def show_about(self): "About ospgui", "ospgui — GUI for ospgrillage\n\n" "Wizard mode: define bridge grillage geometry\n" - "Results mode: view BMD/SFD/TMD/Deflection from .nc files\n\n" + "Results mode: view BMD/SFD/TMD/Deflection or Influence plots\n" + "from NetCDF result files.\n\n" + "Recommended file naming:\n" + "- ordinary: *.res.nc\n" + "- influence lines: *.il.nc\n" + "- influence surfaces: *.is.nc\n\n" + "Tabs are enabled automatically by dataset type.\n" + "Influence-surface plots default to physical x-z space " + "(or station grid if selected).\n\n" "ospgrillage v" + getattr(__import__("ospgrillage"), "__version__", "unknown"), ) @@ -1125,6 +1529,36 @@ def create_main_content(self): self.results_panel.contour_overlay_combo.currentIndexChanged.connect( self._on_results_control_changed ) + self.results_panel.influence_array_combo.currentIndexChanged.connect( + self._on_influence_array_changed + ) + self.results_panel.influence_component_combo.currentIndexChanged.connect( + self._on_results_control_changed + ) + self.results_panel.influence_target_id_combo.currentIndexChanged.connect( + self._on_influence_target_id_changed + ) + self.results_panel.influence_line_set_combo.currentIndexChanged.connect( + self._on_results_control_changed + ) + self.results_panel.influence_line_view_combo.currentIndexChanged.connect( + self._on_influence_line_view_changed + ) + self.results_panel.influence_line_axis_combo.currentIndexChanged.connect( + self._on_results_control_changed + ) + self.results_panel.influence_surface_mode_combo.currentIndexChanged.connect( + self._on_results_control_changed + ) + self.results_panel.influence_surface_view_combo.currentIndexChanged.connect( + self._on_results_control_changed + ) + self.results_panel.influence_surface_show_layer_check.stateChanged.connect( + self._on_results_control_changed + ) + self.results_panel.influence_surface_hover_check.stateChanged.connect( + self._on_results_control_changed + ) self.results_panel.btn_back.clicked.connect(self._switch_to_wizard) self.left_stack.addWidget(self.results_panel) # index 1 @@ -1167,7 +1601,15 @@ def create_main_content(self): "diagrams here.

" ) self._result_tab_widgets = {} - for label in ("Deflection", "BMD", "SFD", "TMD", "Shell Contour"): + for label in ( + "Deflection", + "BMD", + "SFD", + "TMD", + "Shell Contour", + "Influence Line", + "Influence Surface", + ): if _WEBENGINE_AVAILABLE: tab = QWebEngineView() tab.setHtml(_placeholder) @@ -1648,35 +2090,76 @@ def _open_results_file(self): self._model_proxy = proxy self._results = ds + self._results_kind = _classify_results_kind(ds) + summary_name = os.path.basename(file_name) + loadcase_count = len(ds.coords["Loadcase"].values) if "Loadcase" in ds.coords else 0 + kind_label = { + "ordinary": "Ordinary results", + "influence_line": "Influence line results", + "influence_surface": "Influence surface results", + }.get(self._results_kind, "Results") + summary_text = ( + f"{kind_label} | {len(ds.data_vars)} variables, {loadcase_count} load cases" + ) # Populate controls loadcase_names = [] - if "Loadcase" in ds.coords: - loadcase_names = [str(lc) for lc in ds.coords["Loadcase"].values] + if "Loadcase" in self._results.coords: + loadcase_names = [str(lc) for lc in self._results.coords["Loadcase"].values] self.results_panel.populate_loadcases(loadcase_names) - self.results_panel.update_available_members(proxy) - - n_vars = len(ds.data_vars) - n_lc = len(loadcase_names) - self.results_panel.set_file_info( - os.path.basename(file_name), - f"{n_vars} variables, {n_lc} load cases", + self.results_panel.update_available_members(self._model_proxy) + self.results_panel.set_influence_mode( + self._results_kind != "ordinary", + self._results_kind, ) + if self._results_kind != "ordinary": + self.results_panel.populate_influence_controls(self._results) + self.results_panel.set_influence_line_axis_enabled( + self.results_panel.selected_influence_line_view() != "path" + ) + + self.results_panel.set_file_info(summary_name, summary_text) # Shell contour visibility - is_shell = proxy.model_type == "shell_beam" + is_shell = self._model_proxy.model_type == "shell_beam" self.results_panel.set_shell_contour_visible(is_shell) - # Enable/disable the Shell Contour tab + # Enable/disable tabs based on dataset kind for i in range(self.results_tabs.count()): - if self.results_tabs.tabText(i) == "Shell Contour": - self.results_tabs.setTabEnabled(i, is_shell) - break - - # Mark all result tabs stale and switch mode - self._stale_tabs = {"BMD", "SFD", "TMD", "Deflection"} - if is_shell: - self._stale_tabs.add("Shell Contour") + label = self.results_tabs.tabText(i) + if self._results_kind == "ordinary": + enabled = label in {"Deflection", "BMD", "SFD", "TMD"} + if label == "Shell Contour": + enabled = is_shell + elif self._results_kind == "influence_line": + enabled = label == "Influence Line" + else: + enabled = label == "Influence Surface" + self.results_tabs.setTabEnabled(i, enabled) + + if self._results_kind == "ordinary": + for i in range(self.results_tabs.count()): + if self.results_tabs.tabText(i) == "Deflection": + self.results_tabs.setCurrentIndex(i) + break + + # Mark all relevant tabs stale and switch mode + if self._results_kind == "ordinary": + self._stale_tabs = {"BMD", "SFD", "TMD", "Deflection"} + if is_shell: + self._stale_tabs.add("Shell Contour") + elif self._results_kind == "influence_line": + self._stale_tabs = {"Influence Line"} + for i in range(self.results_tabs.count()): + if self.results_tabs.tabText(i) == "Influence Line": + self.results_tabs.setCurrentIndex(i) + break + else: + self._stale_tabs = {"Influence Surface"} + for i in range(self.results_tabs.count()): + if self.results_tabs.tabText(i) == "Influence Surface": + self.results_tabs.setCurrentIndex(i) + break self._switch_to_results() # Grey out contour controls unless Shell Contour tab is active @@ -1690,7 +2173,12 @@ def _open_results_file(self): # ------------------------------------------------------------------ def _on_results_control_changed(self, _=None): """Slot: loadcase or member filter changed — debounced refresh.""" - self._stale_tabs = {"BMD", "SFD", "TMD", "Deflection", "Shell Contour"} + if self._results_kind == "ordinary": + self._stale_tabs = {"BMD", "SFD", "TMD", "Deflection", "Shell Contour"} + elif self._results_kind == "influence_line": + self._stale_tabs = {"Influence Line"} + else: + self._stale_tabs = {"Influence Surface"} # Debounce: multiple checkbox changes in quick succession are # collapsed into a single render via a short single-shot timer. if not hasattr(self, "_debounce_timer"): @@ -1708,6 +2196,32 @@ def _on_result_tab_changed(self, index): self.results_panel.set_shell_contour_enabled(label == "Shell Contour") self._refresh_current_result_tab() + def _on_influence_array_changed(self, _=None): + """Refresh target selectors when the influence response array changes.""" + if self._results is None: + return + array_name = self.results_panel.selected_influence_array() + if array_name: + self.results_panel._update_influence_component_and_target( + self._results, array_name + ) + self._on_results_control_changed() + + def _on_influence_target_id_changed(self, _=None): + """Refresh influence components when target ID changes.""" + if self._results is None: + return + array_name = self.results_panel.selected_influence_array() + if array_name: + self.results_panel._update_influence_components(self._results, array_name) + self._on_results_control_changed() + + def _on_influence_line_view_changed(self, _=None): + """Keep IL abscissa control aligned with selected IL view mode.""" + view_mode = self.results_panel.selected_influence_line_view() + self.results_panel.set_influence_line_axis_enabled(view_mode != "path") + self._on_results_control_changed() + def _refresh_current_result_tab(self): """Render only the currently visible result tab if it is stale.""" if self._model_proxy is None or self._results is None: @@ -1730,7 +2244,131 @@ def _refresh_current_result_tab(self): return try: - if label == "Shell Contour": + if label == "Influence Line": + array_name = self.results_panel.selected_influence_array() + component = self.results_panel.selected_influence_component() + target_type, target_id = self.results_panel.selected_influence_target() + if not array_name or not component or target_id is None: + raise ValueError( + "Influence line selection is incomplete. " + "Choose response array/component and a valid location." + ) + line_view = self.results_panel.selected_influence_line_view() + il_kwargs = dict( + array=array_name, + component=component, + load_coord=( + "station" + if line_view == "path" + else self.results_panel.selected_influence_line_axis() + ), + ) + if target_type == "Node": + il_kwargs["node"] = target_id + elif target_type == "Element": + il_kwargs["element"] = target_id + + selected_study = self.results_panel.selected_influence_line_study() + if "InfluenceLine" in self._results.dims and selected_study == "All": + overlay_ils = {} + for influence_line in self._results.coords["InfluenceLine"].values: + overlay_ils[str(influence_line)] = og.create_influence_line( + ds=self._results, + influence_line=str(influence_line), + **il_kwargs, + ).get() + plot_data = overlay_ils + title = f"{component} influence lines" + elif "InfluenceLine" in self._results.dims: + plot_data = og.create_influence_line( + ds=self._results, + influence_line=selected_study, + **il_kwargs, + ).get() + title = f"{component} influence line" + else: + plot_data = og.create_influence_line( + ds=self._results, + **il_kwargs, + ).get() + title = f"{component} influence line" + + fig = og.plot_il( + plot_data, + backend="plotly", + show=False, + title=title, + view=line_view, + dataset=self._results, + members=members, + show_nodes=True, + show_node_labels=False, + ordinate_aspect=0.9, + ) + elif label == "Influence Surface": + array_name = self.results_panel.selected_influence_array() + component = self.results_panel.selected_influence_component() + target_type, target_id = self.results_panel.selected_influence_target() + if not array_name or not component or target_id is None: + raise ValueError( + "Influence surface selection is incomplete. " + "Choose response array/component and a valid location." + ) + show_surface_layer = self.results_panel.is_influence_surface_layer_visible() + enable_surface_hover = self.results_panel.is_influence_surface_hover_enabled() + selected_mode = self.results_panel.selected_influence_surface_mode() + x_coord, y_coord, coordinate_space = _resolve_influence_surface_axes( + self._results, selected_mode + ) + if selected_mode in {"station", "stations"} and coordinate_space != "station": + self.statusbar.showMessage( + "Station coordinates unavailable in this file; using physical x-z.", + 4000, + ) + is_kwargs = { + "ds": self._results, + "array": array_name, + "component": component, + "x_coord": x_coord, + "y_coord": y_coord, + } + if target_type == "Node": + is_kwargs["node"] = target_id + elif target_type == "Element": + is_kwargs["element"] = target_id + isurface = og.create_influence_surface(**is_kwargs).get() + finite_count = int(np.isfinite(np.asarray(isurface.values, dtype=float)).sum()) + if finite_count < 3: + raise ValueError( + "Selected component/location has too few valid points " + "to build an influence surface. Choose another component." + ) + model_fig = og.plot_model( + self._model_proxy, + backend="plotly", + show=False, + members=members, + show_nodes=True, + show_supports=True, + show_rigid_links=True, + title=None, + ) + if show_surface_layer: + fig = og.plot_is( + isurface, + backend="plotly", + show=False, + ax=model_fig, + title=f"{component} influence surface", + view=self.results_panel.selected_influence_surface_view(), + coordinate_space=coordinate_space, + opacity=0.78, + surface_hover=enable_surface_hover, + ) + else: + fig = model_fig + fig.update_layout(title=f"{component} influence surface (model only)") + elif label == "Shell Contour": # Shell contour — reads directly from Dataset comp = self.results_panel.contour_component_combo.currentText() cs = self.results_panel.contour_colorscale_combo.currentText() diff --git a/src/ospgrillage/postprocessing.py b/src/ospgrillage/postprocessing.py index 9cb6a58b..00608d4a 100644 --- a/src/ospgrillage/postprocessing.py +++ b/src/ospgrillage/postprocessing.py @@ -12,7 +12,9 @@ import matplotlib.pyplot as plt import opsvis as opsv import numpy as np +import re from typing import TYPE_CHECKING, Union +import xarray as xr # if TYPE_CHECKING: from ospgrillage.load import ShapeFunction @@ -20,10 +22,16 @@ __all__ = [ "Envelope", + "InfluenceLine", + "InfluenceSurface", "Members", "PostProcessor", "create_envelope", "model_proxy_from_results", + "create_influence_line", + "create_influence_surface", + "plot_il", + "plot_is", "plot_force", "plot_bmd", "plot_sfd", @@ -54,6 +62,11 @@ def __init__(self, ds): # Build node_spec dict: {tag: {"coordinate": [x, y, z]}} coords_da = ds["node_coordinates"] + # Combined influence datasets may carry an extra leading study dimension. + # Geometry is identical across studies, so take the first slice. + for dim in list(coords_da.dims): + if dim not in {"Node", "Axis"}: + coords_da = coords_da.isel({dim: 0}) self._node_spec = {} for tag in coords_da.coords["Node"].values: self._node_spec[int(tag)] = { @@ -63,6 +76,56 @@ def __init__(self, ds): # Build member-element mapping from JSON attr self._members = json.loads(ds.attrs.get("member_elements", "{}")) + # Optional beam element connectivity (element -> [node_i, node_j]) for + # robust model rendering from saved files even when member node lists are + # not ordered along the geometric line. + self._element_nodes = {} + ele_nodes = ds.get("ele_nodes", None) + if ele_nodes is not None: + ele_da = ele_nodes + for dim in list(ele_da.dims): + if dim not in {"Element", "Nodes"}: + ele_da = ele_da.isel({dim: 0}) + if {"Element", "Nodes"}.issubset(set(ele_da.dims)): + for ele_tag in ele_da.coords["Element"].values: + node_vals = np.asarray(ele_da.sel(Element=ele_tag).values).tolist() + node_tags = [] + for node_val in node_vals: + try: + node_f = float(node_val) + except (TypeError, ValueError): + continue + if np.isnan(node_f): + continue + node_tags.append(int(node_f)) + if len(node_tags) >= 2: + self._element_nodes[int(ele_tag)] = tuple(node_tags[:2]) + + # Optional shell connectivity (for shell-beam model visualisation when + # plotting from saved results files without a live OspGrillage object). + self._shell_element_nodes = [] + ele_nodes_shell = ds.get("ele_nodes_shell", None) + if ele_nodes_shell is not None: + shell_da = ele_nodes_shell + for dim in list(shell_da.dims): + if dim not in {"Element", "Nodes"}: + shell_da = shell_da.isel({dim: 0}) + if {"Element", "Nodes"}.issubset(set(shell_da.dims)): + for node_row in np.asarray(shell_da.values): + node_tags = [] + for node_tag in np.asarray(node_row).tolist(): + if node_tag is None: + continue + try: + node_val = float(node_tag) + except (TypeError, ValueError): + continue + if np.isnan(node_val): + continue + node_tags.append(int(node_val)) + if len(node_tags) >= 3: + self._shell_element_nodes.append(tuple(node_tags)) + # Reconstruct common_grillage_element_z_group (member → list of # z-group indices, e.g. {"interior_main_beam": [0, 1, 2]}). self.common_grillage_element_z_group = {} @@ -130,6 +193,9 @@ def model_proxy_from_results(ds): # default"; ``title=None`` means "no title"; ``title="..."`` is a custom # override. _AUTO = object() +_LOADCASE_POSITION_RE = re.compile( + r"^(?P.+) at global position \[(?P[-+\d.eE]+),(?P[-+\d.eE]+),(?P[-+\d.eE]+)\]$" +) # --------------------------------------------------------------------------- @@ -203,6 +269,157 @@ def create_envelope(**kwargs): return Envelope(**kwargs) +def create_influence_line(**kwargs): + """ + Create an influence line object from stored moving-load or grid-load results. + + The helper extracts one response quantity from an ``xarray.Dataset`` and + reindexes it against one load-position coordinate, typically ``x`` for a + driving-lane influence line. + + :param ds: Result DataSet from :func:`~ospgrillage.osp_grillage.OspGrillage.get_results`. + :type ds: xarray.Dataset + :param component: Specific response component to extract. + :type component: str + :param array: Data array to query, e.g. ``"displacements"``, ``"forces"``, + ``"forces_beam"``, or ``"forces_shell"``. + :type array: str, optional + :param load_coord: Load-position coordinate to use as the influence-line axis. + One of ``"x"``, ``"y"``, ``"z"``, or cumulative path ``"station"``. + Defaults to ``"x"``. + :type load_coord: str, optional + :param loadcase: Optional load case name or list of names to include. + :type loadcase: str or list[str], optional + :param values: Optional explicit ``(x, y, z)`` load positions corresponding to + the selected load cases. If omitted, positions are parsed from the + ``Loadcase`` coordinate strings. + :type values: list[tuple[float, float, float]], optional + :param node: Optional node tag or list of node tags to select. + :type node: int or list[int], optional + :param element: Optional element tag or list of element tags to select. + :type element: int or list[int], optional + :returns: :class:`InfluenceLine` object. + """ + return InfluenceLine(**kwargs) + + +def create_influence_surface(**kwargs): + """ + Create an influence surface object from stored results. + + The helper extracts one response quantity from an ``xarray.Dataset`` and + reshapes it onto a 2D load-position grid, typically ``x`` by ``z`` for + bridge deck influence surfaces. + + :param ds: Result DataSet from :func:`~ospgrillage.osp_grillage.OspGrillage.get_results`. + :type ds: xarray.Dataset + :param component: Specific response component to extract. + :type component: str + :param array: Data array to query, e.g. ``"displacements"``, ``"forces"``, + ``"forces_beam"``, or ``"forces_shell"``. + :type array: str, optional + :param x_coord: Load-position coordinate to use on the first surface axis. + Defaults to ``"x"``. Station-based datasets may also use + ``"longitudinal_station"`` or ``"transverse_station"``. + :type x_coord: str, optional + :param y_coord: Load-position coordinate to use on the second surface axis. + Defaults to ``"z"``. Station-based datasets may also use + ``"longitudinal_station"`` or ``"transverse_station"``. + :type y_coord: str, optional + :param loadcase: Optional load case name or list of names to include. + :type loadcase: str or list[str], optional + :param values: Optional explicit ``(x, y, z)`` load positions corresponding to + the selected load cases. If omitted, positions are parsed from the + ``Loadcase`` coordinate strings. + :type values: list[tuple[float, float, float]], optional + :param node: Optional node tag or list of node tags to select. + :type node: int or list[int], optional + :param element: Optional element tag or list of element tags to select. + :type element: int or list[int], optional + :returns: :class:`InfluenceSurface` object. + """ + return InfluenceSurface(**kwargs) + + +def _parse_loadcase_position(loadcase_name): + match = _LOADCASE_POSITION_RE.match(str(loadcase_name)) + if not match: + raise ValueError( + "Unable to determine load position from Loadcase={!r}. " + "Expected names like ' at global position [x,y,z]' " + "or pass values=[(x, y, z), ...].".format(loadcase_name) + ) + return ( + float(match.group("x")), + float(match.group("y")), + float(match.group("z")), + ) + + +def _normalise_loadcase_selection(loadcase): + if loadcase is None: + return None + if isinstance(loadcase, str): + return [loadcase] + return list(loadcase) + + +def _select_response_data(ds, array, component, node=None, element=None, loadcase=None): + da = getattr(ds, array) + sel_kwargs = {"Component": component} + if "Node" in da.dims and node is not None: + sel_kwargs["Node"] = node + if "Element" in da.dims and element is not None: + sel_kwargs["Element"] = element + if loadcase is not None: + sel_kwargs["Loadcase"] = _normalise_loadcase_selection(loadcase) + return da.sel(**sel_kwargs) + + +def _get_position_index(da, values=None): + loadcases = da.coords["Loadcase"].values.tolist() + if values is None: + if all( + coord in da.coords + for coord in ("load_position_x", "load_position_y", "load_position_z") + ): + positions = list( + zip( + da.coords["load_position_x"].values.tolist(), + da.coords["load_position_y"].values.tolist(), + da.coords["load_position_z"].values.tolist(), + ) + ) + else: + positions = [_parse_loadcase_position(name) for name in loadcases] + else: + positions = [tuple(map(float, position)) for position in values] + if len(positions) != len(loadcases): + raise ValueError( + "values= must have the same length as the selected Loadcase coordinate" + ) + return xr.DataArray( + np.asarray(positions, dtype=float), + dims=("Loadcase", "position_component"), + coords={ + "Loadcase": da.coords["Loadcase"].values, + "position_component": ["x", "y", "z"], + }, + name="load_position", + ) + + +def _cumulative_path_station(position_index): + """Return cumulative path distance for each load position.""" + coords = np.asarray(position_index.values, dtype=float) + if len(coords) == 0: + return np.asarray([], dtype=float) + if len(coords) == 1: + return np.asarray([0.0], dtype=float) + deltas = np.diff(coords, axis=0) + return np.concatenate([[0.0], np.cumsum(np.linalg.norm(deltas, axis=1))]) + + class Envelope: """ Class for defining envelope of :class:`~ospgrillage.osp_grillage.OspGrillage`'s @@ -304,6 +521,814 @@ def get(self): return getattr(da, self.selected_xarray_command)(dim="Loadcase") +class _BaseInfluence: + """Shared implementation for influence-line and influence-surface helpers.""" + + def __init__(self, ds, component: str = None, **kwargs): + self.ds = ds + self.component = component if component is not None else kwargs.get("load_effect", None) + self.array = kwargs.get("array", "displacements") + self.node = kwargs.get("node", None) + self.element = kwargs.get("element", None) + self.loadcase = kwargs.get("loadcase", None) + self.values = kwargs.get("values", None) + self.influence_line = kwargs.get("influence_line", None) + self.influence_surface = kwargs.get("influence_surface", None) + + if ds is None: + raise ValueError("Missing ds=: an xarray Dataset is required") + if self.component is None: + raise ValueError("Missing component=: specify a Component label to extract") + + def _get_base_response(self): + da = _select_response_data( + ds=self.ds, + array=self.array, + component=self.component, + node=self.node, + element=self.element, + loadcase=self.loadcase, + ) + if "InfluenceLine" in da.dims: + if self.influence_line is None: + raise ValueError( + "Dataset contains multiple InfluenceLine studies; specify " + "influence_line= or extract each line explicitly." + ) + da = da.sel(InfluenceLine=self.influence_line) + if "InfluenceSurface" in da.dims: + if self.influence_surface is None: + raise ValueError( + "Dataset contains multiple InfluenceSurface studies; specify " + "influence_surface= explicitly." + ) + da = da.sel(InfluenceSurface=self.influence_surface) + if "Loadcase" not in da.dims: + raise ValueError( + "Selected response does not have a Loadcase dimension. Influence queries " + "require results for multiple loading positions." + ) + return da, _get_position_index(da, values=self.values) + + +class InfluenceLine(_BaseInfluence): + """ + Class for extracting a response history against one load-position coordinate. + + Call :func:`InfluenceLine.get` to obtain an ``xarray.DataArray`` indexed by + the selected load-position axis. + """ + + def __init__(self, ds, component: str = None, **kwargs): + super().__init__(ds=ds, component=component, **kwargs) + self.load_coord = kwargs.get("load_coord", "x") + if self.load_coord not in {"x", "y", "z", "station"}: + raise ValueError("load_coord must be one of 'x', 'y', 'z', or 'station'") + + def get(self): + da, position_index = self._get_base_response() + station = _cumulative_path_station(position_index) + da = da.assign_coords( + x=("Loadcase", position_index.sel(position_component="x").values), + y=("Loadcase", position_index.sel(position_component="y").values), + z=("Loadcase", position_index.sel(position_component="z").values), + station=("Loadcase", station), + ) + da = da.swap_dims({"Loadcase": self.load_coord}).sortby(self.load_coord) + return da.drop_vars("Loadcase", errors="ignore") + + +class InfluenceSurface(_BaseInfluence): + """ + Class for extracting a response field across a 2D load-position grid. + + Call :func:`InfluenceSurface.get` to obtain an ``xarray.DataArray`` indexed by + two load-position coordinates. + """ + + def __init__(self, ds, component: str = None, **kwargs): + super().__init__(ds=ds, component=component, **kwargs) + self.x_coord = kwargs.get("x_coord", "x") + self.y_coord = kwargs.get("y_coord", "z") + valid_coords = {"x", "y", "z", "longitudinal_station", "transverse_station"} + if self.x_coord not in valid_coords or self.y_coord not in valid_coords: + raise ValueError( + "x_coord and y_coord must each be one of 'x', 'y', 'z', " + "'longitudinal_station', or 'transverse_station'" + ) + if self.x_coord == self.y_coord: + raise ValueError("x_coord and y_coord must be different coordinates") + + def get(self): + da, position_index = self._get_base_response() + da = da.assign_coords( + x=("Loadcase", position_index.sel(position_component="x").values), + y=("Loadcase", position_index.sel(position_component="y").values), + z=("Loadcase", position_index.sel(position_component="z").values), + ) + if "load_position_longitudinal_station" in da.coords: + da = da.assign_coords( + longitudinal_station=( + "Loadcase", + da.coords["load_position_longitudinal_station"].values, + ) + ) + if "load_position_transverse_station" in da.coords: + da = da.assign_coords( + transverse_station=( + "Loadcase", + da.coords["load_position_transverse_station"].values, + ) + ) + if da.indexes.get("Loadcase", None) is not None and not da.indexes["Loadcase"].is_unique: + raise ValueError("Loadcase coordinates must be unique before building an influence surface") + da = da.set_index(Loadcase=[self.x_coord, self.y_coord]).unstack("Loadcase") + return da.sortby(self.x_coord).sortby(self.y_coord) + + +def _normalise_influence_line_input(il): + """Return ``[(label, data_array), ...]`` for one or more influence lines.""" + if isinstance(il, xr.DataArray): + label = il.name or "Influence Line" + return [(label, il)] + if isinstance(il, dict): + return [(str(label), da) for label, da in il.items()] + if isinstance(il, (list, tuple)): + lines = [] + for idx, da in enumerate(il, start=1): + label = getattr(da, "name", None) or f"Influence Line {idx}" + lines.append((label, da)) + return lines + raise TypeError("plot_il expects a DataArray, list/tuple of DataArrays, or dict of labelled DataArrays") + + +def _plot_il_path_plotly(lines, **kwargs): + """Plot one or more influence lines along their load paths on the bridge model.""" + go = _import_plotly() + dataset = kwargs.get("dataset", None) + if dataset is None: + raise ValueError("plot_il(view='path') requires dataset= with model geometry metadata") + + proxy = model_proxy_from_results(dataset) + fig = kwargs.get("ax", None) + if fig is None: + fig = plot_model( + proxy, + backend="plotly", + show=False, + title=kwargs.get("title", None), + members=kwargs.get("members", None), + show_nodes=kwargs.get("show_nodes", True), + show_node_labels=kwargs.get("show_node_labels", False), + show_supports=kwargs.get("show_supports", True), + show_rigid_links=kwargs.get("show_rigid_links", True), + ) + + scale = kwargs.get("scale", 1.0) + fill = kwargs.get("fill", True) + fill_alpha = kwargs.get("fill_alpha", 0.35) + marker = kwargs.get("marker", "circle") + show_top_markers = kwargs.get("show_top_markers", False) + show_path_markers = kwargs.get("show_path_markers", True) + color = kwargs.get("color", None) + ordinate_aspect = float(kwargs.get("ordinate_aspect", 0.8)) + default_colours = ["black", "blue", "red", "green", "orange", "purple", "brown", "grey"] + + def _build_segmented_ribbon_mesh(x_vals, z_vals, base_vals, top_vals, eps=1e-12): + """Build a non-self-intersecting ribbon mesh between top and baseline. + + The mesh is built per segment and explicitly split at ordinate + sign-changes (top crossing baseline) to avoid folded quads. + """ + vertices = [] + i_idx, j_idx, k_idx = [], [], [] + + def add_triangle(p0, p1, p2): + start = len(vertices) + vertices.extend([p0, p1, p2]) + i_idx.append(start) + j_idx.append(start + 1) + k_idx.append(start + 2) + + for seg_idx in range(len(x_vals) - 1): + x0, x1 = float(x_vals[seg_idx]), float(x_vals[seg_idx + 1]) + y0, y1 = float(z_vals[seg_idx]), float(z_vals[seg_idx + 1]) + b0, b1 = float(base_vals[seg_idx]), float(base_vals[seg_idx + 1]) + t0, t1 = float(top_vals[seg_idx]), float(top_vals[seg_idx + 1]) + o0 = t0 - b0 + o1 = t1 - b1 + + if not np.isfinite([x0, x1, y0, y1, b0, b1, t0, t1]).all(): + continue + + if abs(o0) <= eps: + o0 = 0.0 + t0 = b0 + if abs(o1) <= eps: + o1 = 0.0 + t1 = b1 + if o0 == 0.0 and o1 == 0.0: + continue + + p0_top = (x0, y0, t0) + p1_top = (x1, y1, t1) + p0_base = (x0, y0, b0) + p1_base = (x1, y1, b1) + + if o0 == 0.0: + # Segment starts on baseline: one wedge triangle. + add_triangle(p0_base, p1_top, p1_base) + continue + if o1 == 0.0: + # Segment ends on baseline: one wedge triangle. + add_triangle(p0_top, p1_base, p0_base) + continue + + if o0 * o1 > 0.0: + # Same sign: one quad = two triangles. + add_triangle(p0_top, p1_top, p0_base) + add_triangle(p1_top, p1_base, p0_base) + continue + + # Ordinates change sign -> split at zero crossing. + alpha = abs(o0) / (abs(o0) + abs(o1)) + xc = x0 + alpha * (x1 - x0) + yc = y0 + alpha * (y1 - y0) + bc = b0 + alpha * (b1 - b0) + p_cross = (xc, yc, bc) # top == base at crossing + + add_triangle(p0_top, p_cross, p0_base) + add_triangle(p1_top, p1_base, p_cross) + + if not vertices: + return None + vertices = np.asarray(vertices, dtype=float) + return { + "x": vertices[:, 0], + "y": vertices[:, 1], + "z": vertices[:, 2], + "i": np.asarray(i_idx, dtype=int), + "j": np.asarray(j_idx, dtype=int), + "k": np.asarray(k_idx, dtype=int), + } + + def _iter_valid_runs(mask): + start = None + for idx, is_valid in enumerate(mask): + if is_valid and start is None: + start = idx + elif not is_valid and start is not None: + if idx - start >= 2: + yield start, idx + start = None + if start is not None and len(mask) - start >= 2: + yield start, len(mask) + + for idx, (label, da) in enumerate(lines): + x_vals = np.asarray(da.coords["x"].values, dtype=float) + y_vals = ( + np.asarray(da.coords["y"].values, dtype=float) + if "y" in da.coords + else np.zeros_like(x_vals, dtype=float) + ) + z_vals = np.asarray(da.coords["z"].values, dtype=float) + ord_vals = np.asarray(da.values, dtype=float) * scale + base_vals = -y_vals + top_vals = base_vals + ord_vals + valid = ( + np.isfinite(x_vals) + & np.isfinite(z_vals) + & np.isfinite(base_vals) + & np.isfinite(top_vals) + ) + if np.count_nonzero(valid) < 2: + continue + + x_plot = x_vals.copy() + z_plot = z_vals.copy() + base_plot = base_vals.copy() + top_plot = top_vals.copy() + x_plot[~valid] = np.nan + z_plot[~valid] = np.nan + base_plot[~valid] = np.nan + top_plot[~valid] = np.nan + + if color is not None and not isinstance(color, (list, tuple)): + line_colour = color + elif isinstance(color, (list, tuple)): + line_colour = color[idx] + else: + line_colour = default_colours[idx % len(default_colours)] + line_dict = None + if color is not None and not isinstance(color, (list, tuple)): + line_dict = dict(color=line_colour, width=6) + elif isinstance(color, (list, tuple)): + line_dict = dict(color=line_colour, width=6) + else: + line_dict = dict(color=line_colour, width=6) + top_mode = "lines+markers" if show_top_markers else "lines" + top_marker_dict = dict(symbol=marker, size=4) if show_top_markers and marker else None + + fig.add_trace( + go.Scatter3d( + x=x_plot, + y=z_plot, + z=top_plot, + mode=top_mode, + line=line_dict, + marker=top_marker_dict, + name=label, + ) + ) + base_mode = "lines+markers" if show_path_markers else "lines" + base_marker = dict(symbol=marker, size=4) if show_path_markers and marker else None + fig.add_trace( + go.Scatter3d( + x=x_plot, + y=z_plot, + z=base_plot, + mode=base_mode, + line=dict( + color=line_colour, + width=3, + ), + marker=base_marker, + name=f"{label} baseline", + showlegend=False, + ) + ) + if fill: + finite_ord = np.abs(top_vals[valid] - base_vals[valid]) + ord_span = float(np.nanmax(finite_ord)) if finite_ord.size else 0.0 + zero_tol = float(kwargs.get("zero_tol", max(1e-10, ord_span * 1e-12))) + for run_start, run_end in _iter_valid_runs(valid): + ribbon = _build_segmented_ribbon_mesh( + x_vals=x_vals[run_start:run_end], + z_vals=z_vals[run_start:run_end], + base_vals=base_vals[run_start:run_end], + top_vals=top_vals[run_start:run_end], + eps=zero_tol, + ) + if ribbon is None: + continue + fig.add_trace( + go.Mesh3d( + x=ribbon["x"], + y=ribbon["y"], + z=ribbon["z"], + i=ribbon["i"], + j=ribbon["j"], + k=ribbon["k"], + color=line_colour, + opacity=fill_alpha, + flatshading=True, + showlegend=False, + hoverinfo="skip", + ) + ) + + aspect = _spatial_aspect_ratio(fig) + aspect["z"] = ordinate_aspect + fig.update_layout( + scene=dict( + xaxis_title="x (m)", + yaxis_title="z (m)", + zaxis_title=kwargs.get("ylabel", "ordinate"), + aspectmode="manual", + aspectratio=aspect, + ) + ) + return fig + + +def plot_il(il, **kwargs): + """Plot a 1D influence line DataArray. + + Parameters + ---------- + il : xarray.DataArray + Influence line with exactly one dimension. + backend : {"matplotlib", "plotly"}, default "matplotlib" + Rendering backend. + show : bool, default True + Whether to display the figure immediately. + ax : matplotlib.axes.Axes or plotly.graph_objects.Figure, optional + Existing target to draw into. + title, xlabel, ylabel, color, marker : optional + Standard plotting customisation keywords. + """ + backend = kwargs.get("backend", "matplotlib") + show = kwargs.get("show", True) + title = kwargs.get("title", None) + xlabel = kwargs.get("xlabel", None) + ylabel = kwargs.get("ylabel", "ordinate") + color = kwargs.get("color", None) + marker = kwargs.get("marker", None) + view = kwargs.get("view", "ordinate") + lines = _normalise_influence_line_input(il) + for _, da in lines: + if len(da.dims) != 1: + raise ValueError("plot_il requires 1D DataArray inputs") + + x_name = lines[0][1].dims[0] + if xlabel is None: + xlabel = x_name + + if view == "path": + if backend != "plotly": + raise ValueError("plot_il(view='path') currently requires backend='plotly'") + fig = _plot_il_path_plotly(lines, **kwargs) + if title is not None: + fig.update_layout(title=title) + if show: + fig.show() + return fig + + if backend == "matplotlib": + ax = kwargs.get("ax", None) + if ax is None: + figsize = kwargs.get("figsize", None) + _, ax = plt.subplots(figsize=figsize) + for idx, (label, da) in enumerate(lines): + line_kwargs = {} + if color is not None and not isinstance(color, (list, tuple)): + line_kwargs["color"] = color + elif isinstance(color, (list, tuple)): + line_kwargs["color"] = color[idx] + if marker is not None: + line_kwargs["marker"] = marker + ax.plot( + da.coords[da.dims[0]].values, + da.values, + label=label, + **line_kwargs, + ) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + if title is not None: + ax.set_title(title) + if len(lines) > 1: + ax.legend() + if show: + plt.show() + return ax + + if backend == "plotly": + import plotly.graph_objects as go + + fig = kwargs.get("ax", None) + if fig is None: + fig = go.Figure() + for idx, (label, da) in enumerate(lines): + line_dict = None + if color is not None and not isinstance(color, (list, tuple)): + line_dict = dict(color=color) + elif isinstance(color, (list, tuple)): + line_dict = dict(color=color[idx]) + marker_dict = dict(symbol=marker) if marker else None + fig.add_trace( + go.Scatter( + x=da.coords[da.dims[0]].values, + y=da.values, + mode="lines+markers" if marker else "lines", + line=line_dict, + marker=marker_dict, + name=label, + ) + ) + fig.update_layout( + title=title, + xaxis_title=xlabel, + yaxis_title=ylabel, + ) + if show: + fig.show() + return fig + + raise ValueError("Unknown backend: choose 'matplotlib' or 'plotly'") + + +def _grid_triangles_from_mask(valid_mask): + """Build triangle connectivity from a structured grid validity mask.""" + if valid_mask.ndim != 2: + raise ValueError("valid_mask must be a 2D array") + + rows, cols = valid_mask.shape + compact_index = -np.ones_like(valid_mask, dtype=int) + compact_index[valid_mask] = np.arange(np.count_nonzero(valid_mask)) + + triangles = [] + for row in range(rows - 1): + for col in range(cols - 1): + v00 = (row, col) + v01 = (row, col + 1) + v11 = (row + 1, col + 1) + v10 = (row + 1, col) + + if valid_mask[v00] and valid_mask[v01] and valid_mask[v11]: + triangles.append( + [ + compact_index[v00], + compact_index[v01], + compact_index[v11], + ] + ) + if valid_mask[v00] and valid_mask[v11] and valid_mask[v10]: + triangles.append( + [ + compact_index[v00], + compact_index[v11], + compact_index[v10], + ] + ) + + if triangles: + return np.asarray(triangles, dtype=int) + return np.empty((0, 3), dtype=int) + + +def _prepare_influence_surface_plot_data(isurface, coordinate_space): + """Prepare coordinate arrays and optional triangulation metadata for ``plot_is``.""" + x_name, y_name = isurface.dims + ordinate = np.asarray(isurface.values.T, dtype=float) + + use_physical_coords = ( + coordinate_space in {"auto", "physical"} + and "x" in isurface.coords + and "z" in isurface.coords + and isurface.coords["x"].dims == isurface.dims + and isurface.coords["z"].dims == isurface.dims + ) + if coordinate_space == "station": + use_physical_coords = False + + if use_physical_coords: + x_grid = np.asarray(isurface.coords["x"].values.T, dtype=float) + y_grid = np.asarray(isurface.coords["z"].values.T, dtype=float) + xlabel = "x" + ylabel = "z" + else: + x_axis = np.asarray(isurface.coords[x_name].values, dtype=float) + y_axis = np.asarray(isurface.coords[y_name].values, dtype=float) + x_grid, y_grid = np.meshgrid(x_axis, y_axis) + xlabel = x_name + ylabel = y_name + + valid_mask = ( + np.isfinite(x_grid) + & np.isfinite(y_grid) + & np.isfinite(ordinate) + ) + use_triangulation = bool(use_physical_coords or np.any(~valid_mask)) + + tri_data = None + if use_triangulation: + triangles = _grid_triangles_from_mask(valid_mask) + x_points = x_grid[valid_mask] + y_points = y_grid[valid_mask] + ordinates = ordinate[valid_mask] + + if len(x_points) < 3 or len(triangles) == 0: + raise ValueError( + "Influence surface plotting requires at least three connected " + "valid points when using physical coordinates" + ) + tri_data = { + "x": x_points, + "y": y_points, + "z": ordinates, + "triangles": triangles, + } + + return { + "x_name": x_name, + "y_name": y_name, + "ordinate": ordinate, + "x_grid": x_grid, + "y_grid": y_grid, + "xlabel": xlabel, + "ylabel": ylabel, + "use_triangulation": use_triangulation, + "tri_data": tri_data, + } + + +def plot_is(isurface, **kwargs): + """Plot a 2D influence surface DataArray. + + Parameters + ---------- + isurface : xarray.DataArray + Influence surface with exactly two dimensions. + backend : {"matplotlib", "plotly"}, default "matplotlib" + Rendering backend. + show : bool, default True + Whether to display the figure immediately. + ax : matplotlib.axes.Axes or plotly.graph_objects.Figure, optional + Existing target to draw into. + title, xlabel, ylabel, colorscale : optional + Standard plotting customisation keywords. + """ + if len(isurface.dims) != 2: + raise ValueError("plot_is requires a 2D DataArray") + + backend = kwargs.get("backend", "matplotlib") + show = kwargs.get("show", True) + title = kwargs.get("title", None) + colorscale = kwargs.get("colorscale", "RdBu_r") + view = kwargs.get("view", "contour") + coordinate_space = kwargs.get("coordinate_space", "auto") + contour_levels = kwargs.get("levels", 20) + opacity = kwargs.get("opacity", None) + ordinate_aspect = float(kwargs.get("ordinate_aspect", 0.8)) + + plot_data = _prepare_influence_surface_plot_data(isurface, coordinate_space) + z = plot_data["ordinate"] + x_grid = plot_data["x_grid"] + y_grid = plot_data["y_grid"] + use_triangulation = plot_data["use_triangulation"] + tri_data = plot_data["tri_data"] + xlabel = kwargs.get("xlabel", plot_data["xlabel"]) + ylabel = kwargs.get("ylabel", plot_data["ylabel"]) + + if backend == "matplotlib": + ax = kwargs.get("ax", None) + if ax is None: + figsize = kwargs.get("figsize", None) + if view == "surface3d": + fig = plt.figure(figsize=figsize) + ax = fig.add_subplot(111, projection="3d") + else: + _, ax = plt.subplots(figsize=figsize) + if view == "surface3d": + if use_triangulation: + surface = ax.plot_trisurf( + tri_data["x"], + tri_data["y"], + tri_data["z"], + triangles=tri_data["triangles"], + cmap=colorscale, + ) + else: + surface = ax.plot_surface(x_grid, y_grid, z, cmap=colorscale) + ax.set_zlabel("ordinate") + plt.colorbar(surface, ax=ax, label="ordinate") + else: + if use_triangulation: + import matplotlib.tri as mtri + + tri = mtri.Triangulation( + tri_data["x"], + tri_data["y"], + tri_data["triangles"], + ) + contour = ax.tricontourf( + tri, + tri_data["z"], + levels=contour_levels, + cmap=colorscale, + ) + else: + contour = ax.contourf(x_grid, y_grid, z, levels=contour_levels, cmap=colorscale) + plt.colorbar(contour, ax=ax, label="ordinate") + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + if title is not None: + ax.set_title(title) + if show: + plt.show() + return ax + + if backend == "plotly": + import plotly.graph_objects as go + + fig = kwargs.get("ax", None) + has_existing_traces = fig is not None and len(getattr(fig, "data", ())) > 0 + surface_hover = bool(kwargs.get("surface_hover", not has_existing_traces)) + if fig is None: + fig = go.Figure() + colorbar_x = kwargs.get("colorbar_x", -0.08 if has_existing_traces else 1.02) + colorbar_kw = dict(title="ordinate", x=colorbar_x, xanchor="right") + if view == "surface3d": + if use_triangulation: + fig.add_trace( + go.Mesh3d( + x=tri_data["x"], + y=tri_data["y"], + z=tri_data["z"], + i=tri_data["triangles"][:, 0], + j=tri_data["triangles"][:, 1], + k=tri_data["triangles"][:, 2], + intensity=tri_data["z"], + colorscale=colorscale, + colorbar=colorbar_kw, + showscale=True, + name="Influence Surface", + hovertemplate=( + f"{xlabel}: %{{x:.3f}}
" + f"{ylabel}: %{{y:.3f}}
" + "ordinate: %{intensity:.3g}" + ) if surface_hover else None, + hoverinfo="skip" if not surface_hover else None, + ) + ) + else: + fig.add_trace( + go.Surface( + x=x_grid, + y=y_grid, + z=z, + colorscale=colorscale, + colorbar=colorbar_kw, + hoverinfo="skip" if not surface_hover else None, + ) + ) + fig.update_layout( + title=title, + scene=dict( + xaxis_title=xlabel, + yaxis_title=ylabel, + zaxis_title="ordinate", + aspectmode="manual", + aspectratio=dict(x=1, y=1, z=ordinate_aspect), + ), + ) + else: + if use_triangulation: + mesh_kwargs = {} + if opacity is not None: + mesh_kwargs["opacity"] = opacity + fig.add_trace( + go.Mesh3d( + x=tri_data["x"], + y=tri_data["y"], + z=np.zeros_like(tri_data["z"]), + i=tri_data["triangles"][:, 0], + j=tri_data["triangles"][:, 1], + k=tri_data["triangles"][:, 2], + intensity=tri_data["z"], + colorscale=colorscale, + colorbar=colorbar_kw, + showscale=True, + name="Influence Surface", + hovertemplate=( + f"{xlabel}: %{{x:.3f}}
" + f"{ylabel}: %{{y:.3f}}
" + "ordinate: %{intensity:.3g}" + ) if surface_hover else None, + hoverinfo="skip" if not surface_hover else None, + **mesh_kwargs, + ) + ) + fig.update_layout( + title=title, + scene=dict( + xaxis_title=xlabel, + yaxis_title=ylabel, + zaxis=dict( + title="", + showticklabels=False, + visible=False, + ), + aspectmode="data", + camera=dict( + projection=dict(type="orthographic"), + eye=dict(x=0.0, y=0.0, z=2.2), + up=dict(x=0, y=1, z=0), + ), + ), + ) + aspect = _spatial_aspect_ratio(fig) + aspect["z"] = ordinate_aspect + fig.update_layout( + scene=dict( + aspectmode="manual", + aspectratio=aspect, + ) + ) + else: + fig.add_trace( + go.Heatmap( + x=x_grid[0, :], + y=y_grid[:, 0], + z=z, + colorscale=colorscale, + colorbar=colorbar_kw, + hoverinfo="skip" if not surface_hover else None, + ) + ) + fig.update_layout( + title=title, + xaxis_title=xlabel, + yaxis_title=ylabel, + ) + if show: + fig.show() + return fig + + raise ValueError("Unknown backend: choose 'matplotlib' or 'plotly'") + + # --------------------------------------------------------------------------- # Force component lookup tables # --------------------------------------------------------------------------- @@ -2186,7 +3211,7 @@ def plot_srf( show : bool Display the plot immediately. Default ``True`` for plotly, ``False`` for matplotlib. - **kwargs + kwargs : dict, optional Forwarded to the backend renderer. Accepted keys include *figsize*, *title*, *opacity*, *show_colorbar*, *averaging*, and *ax* / *fig* for composing onto an existing figure. @@ -2319,9 +3344,11 @@ def _extract_support_data(grillage_obj): ``fixity_type`` (``"pin"`` / ``"roller"`` / ``"fixed"`` / ``"other"``). """ + if not hasattr(grillage_obj, "Mesh_obj"): + return [] mesh = grillage_obj.Mesh_obj node_spec = mesh.node_spec - support_type_dict = grillage_obj.edge_support_type_dict + support_type_dict = getattr(grillage_obj, "edge_support_type_dict", {}) # Determine which attribute stores the node → edge-group mapping. edge_recorder = getattr(mesh, "edge_node_recorder", {}) @@ -2388,79 +3415,189 @@ def _extract_mesh_data(grillage_obj): Returns ------- data : dict[str, list[tuple]] - ``{member_name: [(node_i_xyz, node_j_xyz, ele_tag), ...]}`` + ``{member_name: [(node_i_xyz, node_j_xyz, ele_tag, node_i_tag, node_j_tag), ...]}`` all_nodes : dict[int, list] ``{node_tag: [x, y, z]}`` quads : list[tuple] ``[(c1, c2, c3, c4), ...]`` where each ``ci`` is ``[x, y, z]``. Non-empty only for shell-type meshes. """ - mesh = grillage_obj.Mesh_obj - node_spec = mesh.node_spec - z_group_map = grillage_obj.common_grillage_element_z_group - - data = {} # member_name -> list of (coord_i, coord_j, ele_tag) + data = {} # member_name -> list of (coord_i, coord_j, ele_tag, node_i_tag, node_j_tag) all_nodes = {} # node_tag -> [x, y, z] + if hasattr(grillage_obj, "Mesh_obj"): + mesh = grillage_obj.Mesh_obj + node_spec = mesh.node_spec + z_group_map = grillage_obj.common_grillage_element_z_group + + def _add_elements(member_name, ele_list): + entries = data.setdefault(member_name, []) + for ele in ele_list: + tag, ni, nj = ele[0], ele[1], ele[2] + ci = node_spec[ni]["coordinate"] + cj = node_spec[nj]["coordinate"] + entries.append((ci, cj, tag, int(ni), int(nj))) + all_nodes[ni] = ci + all_nodes[nj] = cj + + for member_name in [ + "edge_beam", + "exterior_main_beam_1", + "interior_main_beam", + "exterior_main_beam_2", + ]: + if member_name not in z_group_map: + continue + for z_idx in z_group_map[member_name]: + if z_idx in mesh.z_group_to_ele: + _add_elements(member_name, mesh.z_group_to_ele[z_idx]) + + for member_name in ["start_edge", "end_edge"]: + if member_name not in z_group_map: + continue + for edge_idx in z_group_map[member_name]: + if edge_idx in mesh.edge_group_to_ele: + _add_elements(member_name, mesh.edge_group_to_ele[edge_idx]) + + _add_elements("transverse_slab", mesh.trans_ele) + + quads = [] + links = [] + if hasattr(grillage_obj, "shell_element_command_list"): + grid = getattr(mesh, "grid_number_dict", {}) + for node_list in grid.values(): + valid = [n for n in node_list if isinstance(n, (int, float))] + if len(valid) >= 3: + coords = [node_spec[n]["coordinate"] for n in valid] + quads.append(tuple(coords)) + for n in valid: + all_nodes[n] = node_spec[n]["coordinate"] + + link_dict = getattr(mesh, "link_dict", {}) + for beam_node, slab_nodes in link_dict.items(): + bc = node_spec[beam_node]["coordinate"] + all_nodes[beam_node] = bc + for sn in slab_nodes: + sc = node_spec[sn]["coordinate"] + links.append((sc, bc)) + all_nodes[sn] = sc + + return data, all_nodes, quads, links + + if isinstance(grillage_obj, _ModelProxy): + node_spec = grillage_obj._node_spec + element_node_map = getattr(grillage_obj, "_element_nodes", {}) + + def _flatten_numeric(values): + flat = [] + stack = list(values) if isinstance(values, list) else [values] + while stack: + item = stack.pop(0) + if isinstance(item, list): + stack = list(item) + stack + continue + try: + flat.append(int(item)) + except (TypeError, ValueError): + continue + return flat + + for member_name, info in grillage_obj._members.items(): + entries = data.setdefault(member_name, []) + element_groups = info.get("elements", []) + node_groups = info.get("nodes", []) + for group_idx, nodes in enumerate(node_groups): + elements = ( + element_groups[group_idx] if group_idx < len(element_groups) else [] + ) + element_ids = _flatten_numeric(elements) + + added_from_elements = False + if element_ids and element_node_map: + for etag in element_ids: + node_pair = element_node_map.get(int(etag)) + if not node_pair or len(node_pair) < 2: + continue + ni, nj = int(node_pair[0]), int(node_pair[1]) + if ni not in node_spec or nj not in node_spec: + continue + ci = node_spec[ni]["coordinate"] + cj = node_spec[nj]["coordinate"] + entries.append((ci, cj, int(etag), ni, nj)) + all_nodes[ni] = ci + all_nodes[nj] = cj + added_from_elements = True + if added_from_elements: + continue + + if nodes and isinstance(nodes, list) and len(nodes) == 1 and isinstance(nodes[0], list): + nodes = nodes[0] + if not nodes or len(nodes) < 2: + continue + + if ( + isinstance(nodes[0], list) + and all(isinstance(seg, list) and len(seg) >= 2 for seg in nodes) + ): + for seg_idx, seg in enumerate(nodes): + ni = int(seg[0]) + nj = int(seg[1]) + if ni not in node_spec or nj not in node_spec: + continue + ci = node_spec[ni]["coordinate"] + cj = node_spec[nj]["coordinate"] + tag = int(elements[seg_idx]) if seg_idx < len(elements) else seg_idx + 1 + entries.append((ci, cj, tag, ni, nj)) + all_nodes[ni] = ci + all_nodes[nj] = cj + continue + + for seg_idx, (ni, nj) in enumerate(zip(nodes[:-1], nodes[1:])): + ni = int(ni) + nj = int(nj) + if ni not in node_spec or nj not in node_spec: + continue + ci = node_spec[ni]["coordinate"] + cj = node_spec[nj]["coordinate"] + tag = int(elements[seg_idx]) if seg_idx < len(elements) else seg_idx + 1 + entries.append((ci, cj, tag, ni, nj)) + all_nodes[ni] = ci + all_nodes[nj] = cj + quads = [] + for shell_nodes in getattr(grillage_obj, "_shell_element_nodes", []): + coords = [] + valid = True + for node_tag in shell_nodes: + node_info = node_spec.get(int(node_tag)) + if node_info is None: + valid = False + break + coord = node_info["coordinate"] + coords.append(coord) + all_nodes[int(node_tag)] = coord + if valid and len(coords) >= 3: + quads.append(tuple(coords)) + return data, all_nodes, quads, [] - def _add_elements(member_name, ele_list): - entries = data.setdefault(member_name, []) - for ele in ele_list: - tag, ni, nj = ele[0], ele[1], ele[2] - ci = node_spec[ni]["coordinate"] - cj = node_spec[nj]["coordinate"] - entries.append((ci, cj, tag)) - all_nodes[ni] = ci - all_nodes[nj] = cj - - # Longitudinal members (edge_beam, exterior_main_beam_1, interior, exterior_2) - for member_name in [ - "edge_beam", - "exterior_main_beam_1", - "interior_main_beam", - "exterior_main_beam_2", - ]: - if member_name not in z_group_map: - continue - for z_idx in z_group_map[member_name]: - if z_idx in mesh.z_group_to_ele: - _add_elements(member_name, mesh.z_group_to_ele[z_idx]) + raise TypeError("Unsupported grillage object for mesh extraction") - # Start / end edges - for member_name in ["start_edge", "end_edge"]: - if member_name not in z_group_map: - continue - for edge_idx in z_group_map[member_name]: - if edge_idx in mesh.edge_group_to_ele: - _add_elements(member_name, mesh.edge_group_to_ele[edge_idx]) - - # Transverse slab - _add_elements("transverse_slab", mesh.trans_ele) - - # Shell quad elements and rigid links (only present for ShellLinkMesh) - quads = [] - links = [] # list of (slab_coord, beam_coord) pairs - if hasattr(grillage_obj, "shell_element_command_list"): - grid = getattr(mesh, "grid_number_dict", {}) - for node_list in grid.values(): - # Each entry is [n1, n2, n3, n4] (may contain [] for degenerate grids) - valid = [n for n in node_list if isinstance(n, (int, float))] - if len(valid) >= 3: - coords = [node_spec[n]["coordinate"] for n in valid] - quads.append(tuple(coords)) - for n in valid: - all_nodes[n] = node_spec[n]["coordinate"] - # Rigid links: link_dict maps offset beam node → [slab_node1, slab_node2] - link_dict = getattr(mesh, "link_dict", {}) - for beam_node, slab_nodes in link_dict.items(): - bc = node_spec[beam_node]["coordinate"] - all_nodes[beam_node] = bc - for sn in slab_nodes: - sc = node_spec[sn]["coordinate"] - links.append((sc, bc)) - all_nodes[sn] = sc +def _filter_member_geometry(data, members): + """Filter extracted mesh geometry to selected member groups.""" + if members is None: + return data + selected_names = set(_resolve_members(members, backend="plotly")) + return {name: elements for name, elements in data.items() if name in selected_names} + - return data, all_nodes, quads, links +def _visible_node_ids_from_member_data(data): + """Return node tags referenced by currently visible member line segments.""" + node_ids = set() + for elements in data.values(): + for element_entry in elements: + if len(element_entry) >= 5: + node_ids.add(int(element_entry[3])) + node_ids.add(int(element_entry[4])) + return node_ids def _plot_model_matplotlib( @@ -2468,6 +3605,7 @@ def _plot_model_matplotlib( *, figsize=None, ax=None, + members=None, title=_AUTO, show_nodes=True, show_node_labels=False, @@ -2479,6 +3617,14 @@ def _plot_model_matplotlib( ): """Render a 2-D plan view (x vs z) of the grillage mesh.""" data, all_nodes, quads, links = _extract_mesh_data(grillage_obj) + data = _filter_member_geometry(data, members) + if members is not None: + visible_nodes = _visible_node_ids_from_member_data(data) + all_nodes = { + node_tag: coord + for node_tag, coord in all_nodes.items() + if node_tag in visible_nodes + } if ax is None: fig, ax = plt.subplots(figsize=figsize) @@ -2511,7 +3657,8 @@ def _plot_model_matplotlib( if color_by_member else "k" ) - for ci, cj, etag in elements: + for element_entry in elements: + ci, cj, etag = element_entry[:3] ax.plot([ci[0], cj[0]], [ci[2], cj[2]], "-", color=colour, linewidth=1) if show_element_labels: mx = 0.5 * (ci[0] + cj[0]) @@ -2639,6 +3786,7 @@ def _plot_model_plotly( *, fig=None, figsize=None, + members=None, title=_AUTO, show_nodes=True, show_node_labels=False, @@ -2654,6 +3802,14 @@ def _plot_model_plotly( fig = go.Figure() data, all_nodes, quads, links = _extract_mesh_data(grillage_obj) + data = _filter_member_geometry(data, members) + if members is not None: + visible_nodes = _visible_node_ids_from_member_data(data) + all_nodes = { + node_tag: coord + for node_tag, coord in all_nodes.items() + if node_tag in visible_nodes + } colours = ["grey", "blue", "green", "blue", "red", "red", "orange"] colour_map = dict(zip(_MEMBER_COLOURS.keys(), colours)) @@ -2661,11 +3817,19 @@ def _plot_model_plotly( # Draw elements grouped by member for member_name, elements in data.items(): colour = colour_map.get(member_name, "black") if color_by_member else "black" - xs, ys, zs = [], [], [] - for ci, cj, etag in elements: + xs, ys, zs, hover_text = [], [], [], [] + for element_entry in elements: + ci, cj, etag = element_entry[:3] + ni = int(element_entry[3]) if len(element_entry) > 3 else None + nj = int(element_entry[4]) if len(element_entry) > 4 else None xs.extend([ci[0], cj[0], None]) ys.extend([ci[2], cj[2], None]) zs.extend([-ci[1], -cj[1], None]) + if ni is not None and nj is not None: + text = f"member: {member_name}
element: {etag}
nodes: {ni}-{nj}" + else: + text = f"member: {member_name}
element: {etag}" + hover_text.extend([text, text, None]) fig.add_trace( go.Scatter3d( x=xs, @@ -2674,6 +3838,13 @@ def _plot_model_plotly( mode="lines", line=dict(color=colour, width=3), name=member_name, + text=hover_text, + hovertemplate=( + "%{text}
" + "x: %{x:.3f}
" + "z: %{y:.3f}
" + "plot y: %{z:.3f}" + ), ) ) @@ -2737,7 +3908,8 @@ def _plot_model_plotly( nx = [all_nodes[n][0] for n in ntags] ny = [all_nodes[n][2] for n in ntags] nz = [-all_nodes[n][1] for n in ntags] - text = [str(n) for n in ntags] if show_node_labels else None + text = [str(n) for n in ntags] + node_y = [all_nodes[n][1] for n in ntags] mode = "markers+text" if show_node_labels else "markers" fig.add_trace( go.Scatter3d( @@ -2747,10 +3919,17 @@ def _plot_model_plotly( mode=mode, marker=dict(size=2, color="black"), text=text, + customdata=np.asarray(node_y, dtype=float), textposition="top center", textfont=dict(size=8), name="nodes", showlegend=False, + hovertemplate=( + "node: %{text}
" + "x: %{x:.3f}
" + "z: %{y:.3f}
" + "y: %{customdata:.3f}" + ), ) ) @@ -2758,7 +3937,8 @@ def _plot_model_plotly( if show_element_labels: ex, ey, ez, etexts = [], [], [], [] for member_name, elements in data.items(): - for ci, cj, etag in elements: + for element_entry in elements: + ci, cj, etag = element_entry[:3] ex.append(0.5 * (ci[0] + cj[0])) ey.append(0.5 * (ci[2] + cj[2])) ez.append(-0.5 * (ci[1] + cj[1])) @@ -2882,6 +4062,7 @@ def plot_model(grillage_obj, *, backend="matplotlib", **kwargs): in ( "fig", "figsize", + "members", "title", "show_nodes", "show_node_labels", @@ -2906,6 +4087,7 @@ def plot_model(grillage_obj, *, backend="matplotlib", **kwargs): in ( "figsize", "ax", + "members", "title", "show_nodes", "show_node_labels", diff --git a/tests/test_load.py b/tests/test_load.py index 45c0a872..ad32d4d4 100644 --- a/tests/test_load.py +++ b/tests/test_load.py @@ -1411,8 +1411,6 @@ def test_dkt_triangle_shape_function_centroid_partition_and_equilibrium(): assert np.isclose(sum(Nv), 1.0) assert np.isclose(sum(Nmx), 0.0) assert np.isclose(sum(Nmz), 0.0) - - def test_dkt_triangle_shape_function_matches_1d_hermite_on_edge(): """On a triangle edge, the DKT distributor should reduce to the 1D Hermite beam form.""" Nv, Nmx, Nmz = og.ShapeFunction.dkt_triangle_shape_function( diff --git a/tests/test_ospgui.py b/tests/test_ospgui.py index d2efe51b..32073f69 100644 --- a/tests/test_ospgui.py +++ b/tests/test_ospgui.py @@ -9,14 +9,15 @@ 2. main() reports the missing dependency and exits with code 1. """ -import importlib -import sys -import os - -import pytest - -sys.path.insert(0, os.path.abspath(".")) -sys.path.insert(0, os.path.abspath("../")) +import importlib +import sys +import os + +import pytest +import xarray as xr + +sys.path.insert(0, os.path.abspath(".")) +sys.path.insert(0, os.path.abspath("../")) def test_ospgui_importable_without_pyqt6(): @@ -29,7 +30,7 @@ def test_ospgui_importable_without_pyqt6(): pytest.skip("PyQt6 is installed; cannot verify the absent-flag path") -def test_ospgui_main_exits_when_pyqt6_absent(capsys): +def test_ospgui_main_exits_when_pyqt6_absent(capsys): """main() must call sys.exit(1) and write a helpful message to stderr.""" from ospgrillage.ospgui import main, _PYQT6_AVAILABLE @@ -40,5 +41,102 @@ def test_ospgui_main_exits_when_pyqt6_absent(capsys): main() assert exc_info.value.code == 1 - captured = capsys.readouterr() - assert "pip install" in captured.err + captured = capsys.readouterr() + assert "pip install" in captured.err + + +def test_ospgui_results_kind_classifier(): + """Classification should identify ordinary, IL, and IS datasets robustly.""" + mod = importlib.import_module("ospgrillage.ospgui") + + ordinary = xr.Dataset( + coords={"Loadcase": ["LC1"], "Node": [1], "Component": ["y"]}, + data_vars={"displacements": (["Loadcase", "Node", "Component"], [[[0.0]]])}, + ) + assert mod._classify_results_kind(ordinary) == "ordinary" + + il_attr = xr.Dataset(coords={"Loadcase": ["LC1"]}, attrs={"influence_type": "line"}) + assert mod._classify_results_kind(il_attr) == "influence_line" + + is_attr = xr.Dataset(coords={"Loadcase": ["LC1"]}, attrs={"influence_type": "surface"}) + assert mod._classify_results_kind(is_attr) == "influence_surface" + + il_dim = xr.Dataset(coords={"InfluenceLine": ["Lane 1"], "Loadcase": [0, 1]}) + assert mod._classify_results_kind(il_dim) == "influence_line" + + is_station = xr.Dataset( + coords={ + "Loadcase": [0, 1, 2, 3], + "load_position_longitudinal_station": ("Loadcase", [0, 0, 1, 1]), + "load_position_transverse_station": ("Loadcase", [0, 1, 0, 1]), + } + ) + assert mod._classify_results_kind(is_station) == "influence_surface" + + legacy_line = xr.Dataset( + coords={ + "Loadcase": [0, 1, 2], + "load_position_x": ("Loadcase", [0.0, 1.0, 2.0]), + "load_position_z": ("Loadcase", [0.0, 0.5, 1.0]), + }, + attrs={"influence_name": "Lane IL"}, + ) + assert mod._classify_results_kind(legacy_line) == "influence_line" + + +def test_ospgui_resolve_influence_surface_axes(): + mod = importlib.import_module("ospgrillage.ospgui") + + ds_with_station = xr.Dataset( + coords={ + "Loadcase": [0, 1, 2, 3], + "load_position_longitudinal_station": ("Loadcase", [0, 0, 1, 1]), + "load_position_transverse_station": ("Loadcase", [0, 1, 0, 1]), + } + ) + assert mod._resolve_influence_surface_axes(ds_with_station, "stations") == ( + "longitudinal_station", + "transverse_station", + "station", + ) + assert mod._resolve_influence_surface_axes(ds_with_station, "physical") == ( + "longitudinal_station", + "transverse_station", + "physical", + ) + + ds_without_station = xr.Dataset( + coords={ + "Loadcase": [0, 1], + "load_position_x": ("Loadcase", [0.0, 1.0]), + "load_position_z": ("Loadcase", [2.0, 2.0]), + } + ) + assert mod._resolve_influence_surface_axes(ds_without_station, "stations") == ( + "x", + "z", + "physical", + ) + assert mod._resolve_influence_surface_axes(ds_without_station, "physical") == ( + "x", + "z", + "physical", + ) + + +def test_ospgui_nonempty_components_filters_nan_components(): + mod = importlib.import_module("ospgrillage.ospgui") + + da = xr.DataArray( + data=[ + [[1.0, float("nan"), 2.0], [3.0, float("nan"), 4.0]], + [[5.0, float("nan"), 6.0], [7.0, float("nan"), 8.0]], + ], + dims=("Loadcase", "Node", "Component"), + coords={ + "Loadcase": [0, 1], + "Node": [10, 11], + "Component": ["x", "theta", "y"], + }, + ) + assert mod._nonempty_components(da) == ["x", "y"] diff --git a/tests/test_postprocessing.py b/tests/test_postprocessing.py index 56e58eea..5216d609 100644 --- a/tests/test_postprocessing.py +++ b/tests/test_postprocessing.py @@ -1,9 +1,1169 @@ +import csv +import json + from ospgrillage import __version__ as version +import ospgrillage.osp_grillage as og_model from fixtures import * +import xarray as xr sys.path.insert(0, os.path.abspath("../")) +def test_influence_line_from_loadcase_positions(): + ds = xr.Dataset( + data_vars={ + "displacements": ( + ["Loadcase", "Node", "Component"], + np.array( + [ + [[1.0, 10.0], [2.0, 20.0]], + [[3.0, 30.0], [4.0, 40.0]], + [[5.0, 50.0], [6.0, 60.0]], + ] + ), + ) + }, + coords={ + "Loadcase": [ + "100 kN axle at global position [0.00,0.00,2.00]", + "100 kN axle at global position [1.00,0.00,2.00]", + "100 kN axle at global position [2.00,0.00,2.00]", + ], + "Node": [1, 2], + "Component": ["x", "y"], + }, + ) + + influence_line = og.create_influence_line( + ds=ds, + array="displacements", + load_effect="y", + node=2, + ).get() + + assert influence_line.dims == ("x",) + assert np.allclose(influence_line.coords["x"].values, [0.0, 1.0, 2.0]) + assert np.allclose(influence_line.coords["z"].values, [2.0, 2.0, 2.0]) + assert np.allclose(influence_line.values, [20.0, 40.0, 60.0]) + + +def test_influence_surface_from_explicit_positions(): + ds = xr.Dataset( + data_vars={ + "forces": ( + ["Loadcase", "Element", "Component"], + np.array( + [ + [[10.0]], + [[11.0]], + [[12.0]], + [[13.0]], + ] + ), + ) + }, + coords={ + "Loadcase": ["case_a", "case_b", "case_c", "case_d"], + "Element": [42], + "Component": ["Mz_i"], + }, + ) + + influence_surface = og.create_influence_surface( + ds=ds, + array="forces", + load_effect="Mz_i", + element=42, + values=[ + (0.0, 0.0, 2.0), + (1.0, 0.0, 2.0), + (0.0, 0.0, 3.0), + (1.0, 0.0, 3.0), + ], + ).get() + + assert influence_surface.dims == ("x", "z") + assert np.allclose(influence_surface.coords["x"].values, [0.0, 1.0]) + assert np.allclose(influence_surface.coords["z"].values, [2.0, 3.0]) + assert influence_surface.sel(x=0.0, z=2.0).item() == 10.0 + assert influence_surface.sel(x=1.0, z=2.0).item() == 11.0 + assert influence_surface.sel(x=0.0, z=3.0).item() == 12.0 + assert influence_surface.sel(x=1.0, z=3.0).item() == 13.0 + + +def test_influence_surface_sorts_unsorted_positions(): + ds = xr.Dataset( + data_vars={ + "forces": ( + ["Loadcase", "Element", "Component"], + np.array( + [ + [[13.0]], + [[10.0]], + [[11.0]], + [[12.0]], + ] + ), + ) + }, + coords={ + "Loadcase": ["case_d", "case_a", "case_b", "case_c"], + "Element": [42], + "Component": ["Mz_i"], + }, + ) + + influence_surface = og.create_influence_surface( + ds=ds, + array="forces", + load_effect="Mz_i", + element=42, + values=[ + (1.0, 0.0, 3.0), + (0.0, 0.0, 2.0), + (1.0, 0.0, 2.0), + (0.0, 0.0, 3.0), + ], + ).get() + + assert np.allclose(influence_surface.coords["x"].values, [0.0, 1.0]) + assert np.allclose(influence_surface.coords["z"].values, [2.0, 3.0]) + assert influence_surface.sel(x=0.0, z=2.0).item() == 10.0 + assert influence_surface.sel(x=1.0, z=2.0).item() == 11.0 + assert influence_surface.sel(x=0.0, z=3.0).item() == 12.0 + assert influence_surface.sel(x=1.0, z=3.0).item() == 13.0 + + +def test_influence_line_station_uses_cumulative_path_distance(): + ds = xr.Dataset( + data_vars={ + "displacements": ( + ["Loadcase", "Node", "Component"], + np.array([[[1.0]], [[2.0]], [[3.0]]]), + ) + }, + coords={ + "Loadcase": ["a", "b", "c"], + "Node": [25], + "Component": ["y"], + }, + ) + + il = og.create_influence_line( + ds=ds, + array="displacements", + component="y", + node=25, + load_coord="station", + values=[ + (0.0, 0.0, 0.0), + (3.0, 0.0, 4.0), + (6.0, 0.0, 8.0), + ], + ).get() + + assert il.dims == ("station",) + assert np.allclose(il.coords["station"].values, [0.0, 5.0, 10.0]) + assert np.allclose(il.coords["x"].values, [0.0, 3.0, 6.0]) + assert np.allclose(il.coords["z"].values, [0.0, 4.0, 8.0]) + + +def test_plot_influence_line_matplotlib(): + il = xr.DataArray( + data=np.array([1.0, 2.0, 3.0]), + dims=("x",), + coords={"x": [0.0, 1.0, 2.0], "z": ("x", [2.0, 2.0, 2.0])}, + ) + + ax = og.plot_il(il, title="IL", show=False) + + assert ax.get_title() == "IL" + assert len(ax.lines) == 1 + + +def test_plot_influence_line_overlay_matplotlib(): + il_a = xr.DataArray(data=np.array([1.0, 2.0, 3.0]), dims=("x",), coords={"x": [0.0, 1.0, 2.0]}) + il_b = xr.DataArray(data=np.array([3.0, 2.0, 1.0]), dims=("x",), coords={"x": [0.0, 1.0, 2.0]}) + + ax = og.plot_il({"Lane 1": il_a, "Lane 2": il_b}, title="Overlay", show=False) + + assert ax.get_title() == "Overlay" + assert len(ax.lines) == 2 + assert ax.get_legend() is not None + + +def test_plot_influence_surface_matplotlib(): + isurface = xr.DataArray( + data=np.array([[10.0, 11.0], [12.0, 13.0]]), + dims=("x", "z"), + coords={"x": [0.0, 1.0], "z": [2.0, 3.0]}, + ) + + ax = og.plot_is(isurface, title="IS", show=False) + + assert ax.get_title() == "IS" + assert len(ax.collections) > 0 + + +def test_plot_influence_surface_matplotlib_surface3d(): + isurface = xr.DataArray( + data=np.array([[10.0, 11.0], [12.0, 13.0]]), + dims=("x", "z"), + coords={"x": [0.0, 1.0], "z": [2.0, 3.0]}, + ) + + ax = og.plot_is(isurface, title="IS 3D", view="surface3d", show=False) + + assert ax.get_title() == "IS 3D" + assert hasattr(ax, "plot_surface") + + +def _build_curved_influence_surface(): + data = np.array( + [ + [10.0, 11.0, 12.0], + [13.0, 14.0, 15.0], + [16.0, 17.0, np.nan], + ] + ) + x_phys = np.array( + [ + [0.0, 0.0, 0.0], + [1.0, 1.1, 1.2], + [2.0, 2.2, np.nan], + ] + ) + z_phys = np.array( + [ + [0.0, 1.0, 2.0], + [0.2, 1.2, 2.2], + [0.5, 1.5, np.nan], + ] + ) + return xr.DataArray( + data=data, + dims=("longitudinal_station", "transverse_station"), + coords={ + "longitudinal_station": [0.0, 1.0, 2.0], + "transverse_station": [0.0, 1.0, 2.0], + "x": (("longitudinal_station", "transverse_station"), x_phys), + "z": (("longitudinal_station", "transverse_station"), z_phys), + }, + ) + + +def test_plot_influence_surface_curved_physical_matplotlib(): + isurface = _build_curved_influence_surface() + + ax = og.plot_is( + isurface, + coordinate_space="physical", + title="Curved IS", + show=False, + ) + + assert ax.get_title() == "Curved IS" + assert len(ax.collections) > 0 + + +def test_plot_influence_line_plotly(): + go = pytest.importorskip("plotly.graph_objects") + il = xr.DataArray( + data=np.array([1.0, 2.0, 3.0]), + dims=("x",), + coords={"x": [0.0, 1.0, 2.0]}, + ) + + fig = og.plot_il(il, backend="plotly", title="IL", show=False) + + assert isinstance(fig, go.Figure) + assert fig.layout.title.text == "IL" + + +def test_plot_influence_line_overlay_plotly(): + go = pytest.importorskip("plotly.graph_objects") + il_a = xr.DataArray(data=np.array([1.0, 2.0, 3.0]), dims=("x",), coords={"x": [0.0, 1.0, 2.0]}) + il_b = xr.DataArray(data=np.array([3.0, 2.0, 1.0]), dims=("x",), coords={"x": [0.0, 1.0, 2.0]}) + + fig = og.plot_il({"Lane 1": il_a, "Lane 2": il_b}, backend="plotly", title="Overlay", show=False) + + assert isinstance(fig, go.Figure) + assert len(fig.data) == 2 + + +def test_plot_model_proxy_includes_shell_mesh_when_available(): + go = pytest.importorskip("plotly.graph_objects") + ds = xr.Dataset( + data_vars={ + "node_coordinates": ( + ("Node", "Axis"), + np.array( + [ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 0.0, 1.0], + [0.0, 0.0, 1.0], + ] + ), + ), + "ele_nodes_shell": ( + ("Element", "Nodes"), + np.array([[1.0, 2.0, 3.0, 4.0]]), + ), + }, + coords={ + "Node": [1, 2, 3, 4], + "Axis": ["x", "y", "z"], + "Element": [1], + "Nodes": [0, 1, 2, 3], + }, + attrs={ + "member_elements": "{}", + "model_type": "shell_beam", + }, + ) + proxy = og.model_proxy_from_results(ds) + fig = og.plot_model(proxy, backend="plotly", show=False) + assert isinstance(fig, go.Figure) + assert any(trace.type == "mesh3d" and trace.name == "shell_slab" for trace in fig.data) + + +def test_plot_model_proxy_prefers_ele_nodes_connectivity_over_node_order(): + go = pytest.importorskip("plotly.graph_objects") + member_elements = { + "interior_main_beam": { + "elements": [[1, 2]], + "nodes": [[[1, 3, 2]]], # intentionally misordered + } + } + ds = xr.Dataset( + data_vars={ + "node_coordinates": ( + ("Node", "Axis"), + np.array( + [ + [0.0, 0.0, 0.0], # node 1 + [1.0, 0.0, 0.0], # node 2 + [2.0, 0.0, 0.0], # node 3 + ] + ), + ), + "ele_nodes": ( + ("Element", "Nodes"), + np.array( + [ + [1, 2], # element 1 + [2, 3], # element 2 + ] + ), + ), + }, + coords={ + "Node": [1, 2, 3], + "Axis": ["x", "y", "z"], + "Element": [1, 2], + "Nodes": [0, 1], + }, + attrs={ + "member_elements": json.dumps(member_elements), + "model_type": "beam_only", + }, + ) + proxy = og.model_proxy_from_results(ds) + fig = og.plot_model(proxy, backend="plotly", show=False, show_nodes=True) + assert isinstance(fig, go.Figure) + trace = next(t for t in fig.data if getattr(t, "name", "") == "interior_main_beam") + assert list(trace.x[:6]) == [0.0, 1.0, None, 1.0, 2.0, None] + non_null_text = [txt for txt in trace.text if txt is not None] + assert any("element:" in txt for txt in non_null_text) + assert any("nodes:" in txt for txt in non_null_text) + node_trace = next(t for t in fig.data if getattr(t, "name", "") == "nodes") + assert "node:" in str(node_trace.hovertemplate) + assert set(node_trace.text) == {"1", "2", "3"} + + +def test_plot_model_plotly_member_filter_limits_visible_members(): + go = pytest.importorskip("plotly.graph_objects") + member_elements = { + "interior_main_beam": { + "elements": [[1]], + "nodes": [[[1, 2]]], + }, + "transverse_slab": { + "elements": [[2]], + "nodes": [[[2, 3]]], + }, + } + ds = xr.Dataset( + data_vars={ + "node_coordinates": ( + ("Node", "Axis"), + np.array( + [ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [2.0, 0.0, 0.0], + ] + ), + ), + "ele_nodes": ( + ("Element", "Nodes"), + np.array( + [ + [1, 2], + [2, 3], + ] + ), + ), + }, + coords={ + "Node": [1, 2, 3], + "Axis": ["x", "y", "z"], + "Element": [1, 2], + "Nodes": [0, 1], + }, + attrs={ + "member_elements": json.dumps(member_elements), + "model_type": "beam_only", + }, + ) + proxy = og.model_proxy_from_results(ds) + fig = og.plot_model( + proxy, + backend="plotly", + show=False, + show_nodes=True, + members=og.Members.INTERIOR_MAIN_BEAM, + ) + assert isinstance(fig, go.Figure) + names = {getattr(t, "name", "") for t in fig.data} + assert "interior_main_beam" in names + assert "transverse_slab" not in names + node_trace = next(t for t in fig.data if getattr(t, "name", "") == "nodes") + assert set(node_trace.text) == {"1", "2"} + + +def test_plot_influence_surface_plotly(): + go = pytest.importorskip("plotly.graph_objects") + isurface = xr.DataArray( + data=np.array([[10.0, 11.0], [12.0, 13.0]]), + dims=("x", "z"), + coords={"x": [0.0, 1.0], "z": [2.0, 3.0]}, + ) + + fig = og.plot_is(isurface, backend="plotly", title="IS", show=False) + + assert isinstance(fig, go.Figure) + assert fig.layout.title.text == "IS" + + +def test_plot_influence_surface_plotly_surface3d(): + go = pytest.importorskip("plotly.graph_objects") + isurface = xr.DataArray( + data=np.array([[10.0, 11.0], [12.0, 13.0]]), + dims=("x", "z"), + coords={"x": [0.0, 1.0], "z": [2.0, 3.0]}, + ) + + fig = og.plot_is(isurface, backend="plotly", title="IS 3D", view="surface3d", show=False) + + assert isinstance(fig, go.Figure) + assert len(fig.data) == 1 + assert fig.data[0].type == "surface" + + +def test_normalise_netcdf_filename_semantic_suffixes(): + assert og_model._normalise_netcdf_filename("results", file_tag="res") == "results.res.nc" + assert og_model._normalise_netcdf_filename("lane", file_tag="il") == "lane.il.nc" + assert og_model._normalise_netcdf_filename("surface", file_tag="is") == "surface.is.nc" + assert og_model._normalise_netcdf_filename("already.nc", file_tag="res") == "already.nc" + assert og_model._normalise_netcdf_filename("bridge.res", file_tag="res") == "bridge.res.nc" + + +def test_plot_influence_surface_curved_physical_plotly_contour(): + go = pytest.importorskip("plotly.graph_objects") + isurface = _build_curved_influence_surface() + + fig = og.plot_is( + isurface, + backend="plotly", + coordinate_space="physical", + title="Curved IS", + show=False, + ) + + assert isinstance(fig, go.Figure) + assert len(fig.data) == 1 + assert fig.data[0].type == "mesh3d" + + +def test_plot_influence_surface_curved_physical_plotly_surface3d(): + go = pytest.importorskip("plotly.graph_objects") + isurface = _build_curved_influence_surface() + + fig = og.plot_is( + isurface, + backend="plotly", + coordinate_space="physical", + title="Curved IS 3D", + view="surface3d", + show=False, + ) + + assert isinstance(fig, go.Figure) + assert len(fig.data) == 1 + assert fig.data[0].type == "mesh3d" + + +def test_plot_influence_surface_plotly_overlay_on_model(bridge_model_42_negative): + go = pytest.importorskip("plotly.graph_objects") + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + iss = example_bridge.analyze_influence_surfaces(name="Deck IS Default") + isurface = iss.get_surface( + array="forces", + component="Mz_i", + element=1, + x_coord="longitudinal_station", + y_coord="transverse_station", + ) + proxy = og.model_proxy_from_results(iss.dataset) + model_fig = og.plot_model(proxy, backend="plotly", show=False, show_nodes=False) + fig = og.plot_is( + isurface, + backend="plotly", + show=False, + ax=model_fig, + coordinate_space="physical", + opacity=0.6, + ) + assert isinstance(fig, go.Figure) + assert any(trace.type == "scatter3d" for trace in fig.data) + assert any(trace.type == "mesh3d" for trace in fig.data) + assert fig.layout.scene.aspectmode == "manual" + is_trace = next( + trace + for trace in fig.data + if trace.type == "mesh3d" and getattr(trace, "name", "") == "Influence Surface" + ) + assert float(is_trace.colorbar.x) < 0.0 + assert getattr(is_trace, "hoverinfo", None) == "skip" + + +def test_analyze_influence_line_separate_results(bridge_model_42_negative): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + result_set = example_bridge.analyze_influence_line( + name="Lane IL", + start_point=og.Point(2, 0, 2), + end_point=og.Point(4, 0, 2), + increments=3, + ) + influence_results = example_bridge.get_influence_results("Lane IL") + + assert result_set.kind == "line" + assert influence_results.attrs["influence_type"] == "line" + assert influence_results.attrs["shape_function"] == "linear" + assert "Lane IL" in example_bridge.influence_result_set + assert influence_results.sizes["Loadcase"] == 3 + assert np.allclose(influence_results.coords["load_position_x"].values, [2.0, 3.0, 4.0]) + assert "node_coordinates" in influence_results + assert influence_results.attrs["model_type"] == example_bridge.model_type + + +def test_analyze_influence_lines_returns_result_object(bridge_model_42_negative, tmp_path): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + ils = example_bridge.analyze_influence_lines( + name="Lane IL", + start_point=og.Point(2, 0, 2), + end_point=og.Point(4, 0, 2), + increments=3, + ) + + assert isinstance(ils, og.InfluenceLineResults) + assert ils.dataset.attrs["influence_type"] == "line" + assert ils.dataset.attrs["influence_name"] == "Lane IL" + assert np.allclose(ils.dataset.coords["load_position_x"].values, [2.0, 3.0, 4.0]) + + save_path = tmp_path / "lane_il.nc" + assert ils.save(save_path.as_posix()) == save_path.as_posix() + assert save_path.exists() + + stem_path = tmp_path / "lane_il_export" + saved = ils.to_netcdf(stem_path.as_posix()) + assert saved.endswith(".il.nc") + assert os.path.exists(saved) + + +def test_analyze_influence_lines_combines_named_paths(bridge_model_42_negative): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + path_1 = og.Path( + start_point=og.Point(2, 0, 2), + end_point=og.Point(4, 0, 2), + increments=3, + ) + path_2 = og.Path( + start_point=og.Point(2, 0, 4), + end_point=og.Point(4, 0, 4), + increments=3, + ) + + ils = example_bridge.analyze_influence_lines( + paths={"Lane 1": path_1, "Lane 2": path_2} + ) + + assert isinstance(ils, og.InfluenceLineResults) + assert "InfluenceLine" in ils.dataset.dims + assert list(ils.dataset.coords["InfluenceLine"].values) == ["Lane 1", "Lane 2"] + + overlay = ils.get_line(array="displacements", component="y", node=25) + assert sorted(overlay) == ["Lane 1", "Lane 2"] + ax = ils.plot(array="displacements", component="y", node=25, show=False) + assert len(ax.lines) == 2 + go = pytest.importorskip("plotly.graph_objects") + fig = ils.plot( + array="displacements", + component="y", + node=25, + backend="plotly", + view="path", + show=False, + ) + assert isinstance(fig, go.Figure) + assert any(trace.type == "scatter3d" for trace in fig.data) + assert any(trace.type == "mesh3d" for trace in fig.data) + lane_2_trace = next(trace for trace in fig.data if getattr(trace, "name", "") == "Lane 2") + lane_2_base = next(trace for trace in fig.data if getattr(trace, "name", "") == "Lane 2 baseline") + assert lane_2_trace.mode == "lines" + assert lane_2_base.mode == "lines+markers" + assert np.allclose(np.asarray(lane_2_trace.y, dtype=float), [4.0, 4.0, 4.0]) + assert np.allclose(np.asarray(lane_2_base.z, dtype=float), [0.0, 0.0, 0.0]) + + +def test_plot_il_path_fill_ignores_noise_and_gaps(bridge_model_42_negative): + go = pytest.importorskip("plotly.graph_objects") + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + ils = example_bridge.analyze_influence_lines( + name="Lane IL", + start_point=og.Point(2, 0, 2), + end_point=og.Point(8, 0, 2), + increments=7, + ) + il = ils.get_line(array="displacements", component="y", node=25) + il_with_noise = il.copy( + data=np.asarray([1e-13, -1e-13, 5e-14, np.nan, np.nan, 0.1, -0.2], dtype=float) + ) + x_vals = np.asarray(il_with_noise.coords["x"].values, dtype=float) + + fig = og.plot_il( + il_with_noise, + backend="plotly", + view="path", + dataset=ils.dataset, + show=False, + ) + + mesh_traces = [trace for trace in fig.data if isinstance(trace, go.Mesh3d)] + assert len(mesh_traces) == 1 + x_mesh = np.asarray(mesh_traces[0].x, dtype=float) + expected_min = min(x_vals[5], x_vals[6]) - 1e-9 + expected_max = max(x_vals[5], x_vals[6]) + 1e-9 + assert np.nanmin(x_mesh) >= expected_min + assert np.nanmax(x_mesh) <= expected_max + + +def test_get_combined_influence_line_results(bridge_model_42_negative): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + example_bridge.analyze_influence_line( + name="Lane 1", + start_point=og.Point(2, 0, 2), + end_point=og.Point(4, 0, 2), + increments=3, + ) + example_bridge.analyze_influence_line( + name="Lane 2", + start_point=og.Point(2, 0, 4), + end_point=og.Point(4, 0, 4), + increments=3, + ) + + combined = example_bridge.get_influence_results(names=["Lane 1", "Lane 2"]) + + assert combined.attrs["influence_type"] == "line" + assert combined.attrs["influence_overlay"] == "multi" + assert "InfluenceLine" in combined.dims + assert list(combined.coords["InfluenceLine"].values) == ["Lane 1", "Lane 2"] + assert "loadcase_label" in combined.coords + assert np.allclose( + combined.sel(InfluenceLine="Lane 1").coords["load_position_z"].values[:3], + [2.0, 2.0, 2.0], + equal_nan=False, + ) + + +def test_create_influence_line_from_combined_results(bridge_model_42_negative): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + example_bridge.analyze_influence_line( + name="Lane 1", + start_point=og.Point(2, 0, 2), + end_point=og.Point(4, 0, 2), + increments=3, + ) + example_bridge.analyze_influence_line( + name="Lane 2", + start_point=og.Point(2, 0, 4), + end_point=og.Point(4, 0, 4), + increments=3, + ) + + combined = example_bridge.get_influence_results(names=["Lane 1", "Lane 2"]) + + with pytest.raises(ValueError, match="multiple InfluenceLine studies"): + og.create_influence_line( + ds=combined, + array="displacements", + component="y", + node=25, + ).get() + + il_lane_1 = og.create_influence_line( + ds=combined, + array="displacements", + component="y", + node=25, + influence_line="Lane 1", + ).get() + il_lane_2 = og.create_influence_line( + ds=combined, + array="displacements", + component="y", + node=25, + influence_line="Lane 2", + ).get() + + assert list(il_lane_1.coords["x"].values) == [2.0, 3.0, 4.0] + assert list(il_lane_2.coords["z"].values) == [4.0, 4.0, 4.0] + + +def test_analyze_influence_line_preserves_hermite_shape_function(bridge_model_42_negative): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + example_bridge.analyze_influence_line( + name="Lane IL Hermite", + start_point=og.Point(2, 0, 2), + end_point=og.Point(4, 0, 2), + increments=3, + shape_function="hermite", + ) + influence_results = example_bridge.get_influence_results("Lane IL Hermite") + + assert influence_results.attrs["shape_function"] == "hermite" + + +def test_analyze_influence_surface_separate_results(bridge_model_42_negative): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + result_set = example_bridge.analyze_influence_surface( + name="Deck IS", + x=[2, 4], + z=[2, 3], + ) + influence_results = example_bridge.get_influence_results("Deck IS") + + assert result_set.kind == "surface" + assert influence_results.attrs["influence_type"] == "surface" + assert influence_results.attrs["shape_function"] == "linear" + assert "Deck IS" in example_bridge.influence_result_set + assert influence_results.sizes["Loadcase"] == 4 + assert np.allclose(influence_results.coords["load_position_x"].values, [2.0, 2.0, 4.0, 4.0]) + assert np.allclose(influence_results.coords["load_position_z"].values, [2.0, 3.0, 2.0, 3.0]) + + +def test_analyze_influence_surfaces_returns_result_object(bridge_model_42_negative, tmp_path): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + iss = example_bridge.analyze_influence_surfaces( + name="Deck IS", + x=[2, 4], + z=[2, 3], + ) + + assert isinstance(iss, og.InfluenceSurfaceResults) + assert iss.dataset.attrs["influence_type"] == "surface" + assert iss.dataset.attrs["influence_name"] == "Deck IS" + + isurface = iss.get_surface(array="forces", component="Mz_i", element=1) + assert isurface.dims == ("x", "z") + ax = iss.plot(array="forces", component="Mz_i", element=1, show=False) + assert len(ax.collections) > 0 + + save_path = tmp_path / "deck_is.nc" + assert iss.to_netcdf(save_path.as_posix()) == save_path.as_posix() + assert save_path.exists() + + stem_path = tmp_path / "deck_is_export" + saved = iss.to_netcdf(stem_path.as_posix()) + assert saved.endswith(".is.nc") + assert os.path.exists(saved) + + +def test_analyze_influence_surfaces_defaults_to_mesh_station_grid(bridge_model_42_negative): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + iss = example_bridge.analyze_influence_surfaces(name="Deck IS Default") + + longitudinal = iss.dataset.coords["load_position_longitudinal_station"].values + transverse = iss.dataset.coords["load_position_transverse_station"].values + + assert np.all(np.isin(np.unique(longitudinal), np.asarray(example_bridge.Mesh_obj.nox))) + assert np.all(np.isin(np.unique(transverse), np.asarray(example_bridge.Mesh_obj.noz))) + assert len(longitudinal) == iss.dataset.sizes["Loadcase"] + + station_surface = iss.get_surface( + array="forces", + component="Mz_i", + element=1, + x_coord="longitudinal_station", + y_coord="transverse_station", + ) + assert station_surface.dims == ("longitudinal_station", "transverse_station") + + +@pytest.mark.parametrize( + "bridge_fixture", + ["bridge_model_42_negative", "bridge_42_0_angle_mesh"], +) +def test_influence_surface_station_mapping_preserves_longitudinal_order(request, bridge_fixture): + og.ops.wipeAnalysis() + example_bridge = request.getfixturevalue(bridge_fixture) + + iss = example_bridge.analyze_influence_surfaces(name="Deck IS Station Ordering") + station_surface = iss.get_surface( + array="forces", + component="Mz_i", + element=1, + x_coord="longitudinal_station", + y_coord="transverse_station", + ) + + transverse_stations = np.asarray( + station_surface.coords["transverse_station"].values, + dtype=float, + ) + longitudinal_stations = np.asarray( + station_surface.coords["longitudinal_station"].values, + dtype=float, + ) + assert transverse_stations.size > 0 + assert longitudinal_stations.size > 1 + + first_transverse = float(transverse_stations[0]) + x_line = np.asarray( + station_surface.sel(transverse_station=first_transverse).coords["x"].values, + dtype=float, + ) + finite = np.isfinite(x_line) + assert np.count_nonzero(finite) == longitudinal_stations.size + assert np.all(np.diff(x_line[finite]) >= -1e-9) + + +@pytest.mark.parametrize("mesh_radius", [20.0, 1000.0]) +def test_analyze_influence_surfaces_curved_mesh_station_mapping(ref_bridge_properties, mesh_radius): + og.ops.wipeAnalysis() + I_beam, slab, exterior_I_beam, _ = ref_bridge_properties + + example_bridge = og.create_grillage( + bridge_name=f"Curved IS R={mesh_radius:g}", + long_dim=10, + width=7, + skew=0, + num_long_grid=7, + num_trans_grid=7, + mesh_type="Ortho", + mesh_radius=mesh_radius, + ) + example_bridge.set_member(I_beam, member="interior_main_beam") + example_bridge.set_member(exterior_I_beam, member="exterior_main_beam_1") + example_bridge.set_member(exterior_I_beam, member="exterior_main_beam_2") + example_bridge.set_member(exterior_I_beam, member="edge_beam") + example_bridge.set_member(slab, member="transverse_slab") + example_bridge.set_member(exterior_I_beam, member="start_edge") + example_bridge.set_member(exterior_I_beam, member="end_edge") + example_bridge.create_osp_model(pyfile=False) + + iss = example_bridge.analyze_influence_surfaces(name="Curved Deck IS") + station_surface = iss.get_surface( + array="forces", + component="Mz_i", + element=1, + x_coord="longitudinal_station", + y_coord="transverse_station", + ) + + longitudinal_stations = np.asarray( + station_surface.coords["longitudinal_station"].values, + dtype=float, + ) + transverse_stations = np.asarray( + station_surface.coords["transverse_station"].values, + dtype=float, + ) + assert station_surface.shape == (len(longitudinal_stations), len(transverse_stations)) + assert len(longitudinal_stations) == len(np.asarray(example_bridge.Mesh_obj.nox, dtype=float)) + assert len(transverse_stations) == len(np.asarray(example_bridge.Mesh_obj.noz, dtype=float)) + + for z_station in transverse_stations: + x_line = np.asarray( + station_surface.sel(transverse_station=float(z_station)).coords["x"].values, + dtype=float, + ) + finite = np.isfinite(x_line) + assert np.count_nonzero(finite) == len(longitudinal_stations) + assert np.all(np.diff(x_line[finite]) >= -1e-9) + + for longitudinal_station in longitudinal_stations: + z_line = np.asarray( + station_surface.sel(longitudinal_station=float(longitudinal_station)).coords["z"].values, + dtype=float, + ) + finite = np.isfinite(z_line) + assert np.count_nonzero(finite) == len(transverse_stations) + assert np.all(np.diff(z_line[finite]) >= -1e-9) + + +def test_influence_line_results_to_csv_combined_paths(bridge_model_42_negative, tmp_path): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + + path_1 = og.Path( + start_point=og.Point(2, 0, 2), + end_point=og.Point(4, 0, 2), + increments=3, + ) + path_2 = og.Path( + start_point=og.Point(2, 0, 4), + end_point=og.Point(4, 0, 4), + increments=3, + ) + ils = example_bridge.analyze_influence_lines(paths={"Lane 1": path_1, "Lane 2": path_2}) + + csv_file = tmp_path / "lane_ils.csv" + assert ( + ils.to_csv( + csv_file.as_posix(), + array="displacements", + component="y", + node=25, + load_coord="x", + ) + == csv_file.as_posix() + ) + assert csv_file.exists() + + with open(csv_file, newline="") as fh: + rows = list(csv.DictReader(fh)) + + assert len(rows) == 6 + assert {row["influence_line"] for row in rows} == {"Lane 1", "Lane 2"} + assert all("x" in row and "ordinate" in row for row in rows) + + lane_1_x = sorted( + [float(row["x"]) for row in rows if row["influence_line"] == "Lane 1"] + ) + assert lane_1_x == [2.0, 3.0, 4.0] + + +def test_influence_surface_results_to_csv_grid_and_points(bridge_model_42_negative, tmp_path): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + iss = example_bridge.analyze_influence_surfaces(name="Deck IS Default") + + grid_file = tmp_path / "deck_is_grid.csv" + exported = iss.to_csv( + grid_file.as_posix(), + array="forces", + component="Mz_i", + element=1, + include_physical_coords=True, + ) + + assert isinstance(exported, dict) + assert exported["grid"] == grid_file.as_posix() + points_file = tmp_path / "deck_is_grid_points.csv" + assert exported["points"] == points_file.as_posix() + assert grid_file.exists() + assert points_file.exists() + + with open(grid_file, newline="") as fh: + grid_rows = list(csv.reader(fh)) + assert grid_rows[0][0] == "longitudinal_station" + assert len(grid_rows) > 1 + + with open(points_file, newline="") as fh: + point_rows = list(csv.reader(fh)) + assert point_rows[0] == [ + "longitudinal_station", + "transverse_station", + "x", + "z", + "ordinate", + ] + assert len(point_rows) > 1 + + long_file = tmp_path / "deck_is_long.csv" + assert ( + iss.to_csv( + long_file.as_posix(), + array="forces", + component="Mz_i", + element=1, + layout="long", + include_physical_coords=True, + ) + == long_file.as_posix() + ) + with open(long_file, newline="") as fh: + long_rows = list(csv.reader(fh)) + assert long_rows[0] == [ + "longitudinal_station", + "transverse_station", + "x", + "z", + "ordinate", + ] + assert len(long_rows) > 1 + + +def test_influence_surface_results_plot_defaults_auto_space( + bridge_model_42_negative, + monkeypatch, +): + og.ops.wipeAnalysis() + example_bridge = bridge_model_42_negative + iss = example_bridge.analyze_influence_surfaces(name="Deck IS Default") + + captured = {} + + def _fake_plot(isurface, **kwargs): + captured["coordinate_space"] = kwargs.get("coordinate_space", None) + return kwargs.get("coordinate_space", None) + + monkeypatch.setattr("ospgrillage.osp_grillage.plot_influence_surface", _fake_plot) + + iss.plot( + array="forces", + component="Mz_i", + element=1, + x_coord="longitudinal_station", + y_coord="transverse_station", + show=False, + ) + assert captured["coordinate_space"] is None + + iss.plot( + array="forces", + component="Mz_i", + element=1, + x_coord="longitudinal_station", + y_coord="transverse_station", + coordinate_space="physical", + show=False, + ) + assert captured["coordinate_space"] == "physical" + + +def test_influence_line_midspan_moment_matches_l_over_4(ref_bridge_properties): + og.ops.wipeAnalysis() + I_beam, slab, exterior_I_beam, concrete = ref_bridge_properties + + L = 10.0 + bridge = og.create_grillage( + bridge_name="IL_simple_beam", + long_dim=L, + width=2.0, + skew=0, + num_long_grid=3, + num_trans_grid=3, + edge_beam_dist=1.0, + mesh_type="Ortho", + ) + bridge.set_member(I_beam, member="interior_main_beam") + bridge.set_member(exterior_I_beam, member="exterior_main_beam_1") + bridge.set_member(exterior_I_beam, member="exterior_main_beam_2") + bridge.set_member(exterior_I_beam, member="edge_beam") + bridge.set_member(slab, member="transverse_slab") + bridge.set_member(exterior_I_beam, member="start_edge") + bridge.set_member(exterior_I_beam, member="end_edge") + bridge.create_osp_model(pyfile=False) + + interior_elements = bridge.get_element(member="interior_main_beam", options="elements") + assert len(interior_elements) == 2 + + bridge.analyze_il( + name="unit_axle_il", + start_point=og.Point(0, 0, 1.0), + end_point=og.Point(L, 0, 1.0), + step=0.5, + axle_load=1.0, + ) + il = bridge.get_il( + name="unit_axle_il", + array="forces", + component="Mz_j", + element=interior_elements[0], + load_coord="x", + ) + + midspan_ordinate = float(il.sel(x=L / 2, method="nearest")) + assert np.isclose(midspan_ordinate, -L / 4, atol=0.15) + + +def test_hermite_triangle_region_uses_dkt_style_distribution(bridge_model_42_negative): + og.ops.wipeAnalysis() + bridge = bridge_model_42_negative + + tri_grid_nodes = None + for nodes in bridge.Mesh_obj.grid_number_dict.values(): + if len(nodes) == 3: + tri_grid_nodes = nodes + break + + assert tri_grid_nodes is not None + coords = [bridge.Mesh_obj.node_spec[node]["coordinate"] for node in tri_grid_nodes] + point = [ + sum(coord[0] for coord in coords) / 3, + 0, + sum(coord[2] for coord in coords) / 3, + ] + + load_cmd = bridge._assign_load_to_four_node(point=point, mag=1.0, shape_func="hermite") + + assert len(load_cmd) == 3 + vertical_sum = sum(command[1][2] for command in load_cmd) + mx_sum = sum(command[1][4] for command in load_cmd) + mz_sum = sum(command[1][6] for command in load_cmd) + assert any(abs(command[1][4]) > 0 for command in load_cmd) + assert any(abs(command[1][6]) > 0 for command in load_cmd) + assert np.isclose(vertical_sum, 1.0) + assert np.isclose(mx_sum, 0.0) + assert np.isclose(mz_sum, 0.0) + + def test_envelope(bridge_model_42_negative): # test functionality of Envelope class and output og.ops.wipeAnalysis()