Advanced Results Processing#

This notebook demonstrates advanced techniques for extracting and processing ospgrillage results using xarray. It builds on the Super-T Tutorial.

Setup#

We first recreate the Super-T bridge model from the tutorial. See that notebook for detailed explanations.

[1]:

import ospgrillage as og
import numpy as np

# Units
kilo, milli, N, m = 1e3, 1e-3, 1, 1
mm = milli * m
m2, m3, m4 = m**2, m**3, m**4
kN = kilo * N
MPa = N / mm**2
GPa = kilo * MPa

[2]:

concrete = og.create_material(material="concrete", code="AS5100-2017", grade="65MPa")

longitudinal_section = og.create_section(A=1.025*m2, J=0.1878*m3, Iz=0.3694*m4, Iy=0.3634*m4, Az=0.4979*m2, Ay=0.309*m2)
edge_longitudinal_section = og.create_section(A=0.934*m2, J=0.1857*m3, Iz=0.3478*m4, Iy=0.213602*m4, Az=0.444795*m2, Ay=0.258704*m2)
transverse_section = og.create_section(A=0.504*m2, J=5.22303e-3*m3, Iy=0.32928*m4, Iz=1.3608e-3*m4, Ay=0.42*m2, Az=0.42*m2, unit_width=True)
end_transverse_section = og.create_section(A=0.504/2*m2, J=2.5e-3*m3, Iy=2.73e-2*m4, Iz=6.8e-4*m4, Ay=0.21*m2, Az=0.21*m2, unit_width=True)

longitudinal_beam = og.create_member(section=longitudinal_section, material=concrete)
edge_longitudinal_beam = og.create_member(section=edge_longitudinal_section, material=concrete)
transverse_slab = og.create_member(section=transverse_section, material=concrete)
end_transverse_slab = og.create_member(section=end_transverse_section, material=concrete)

[3]:

L = 33.5 * m
w = 11.565 * m

model = og.create_grillage(
    bridge_name="Super-T 33_5m", long_dim=L, width=w, skew=0,
    num_long_grid=7, num_trans_grid=11, edge_beam_dist=1.05 * m,
)
model.set_member(longitudinal_beam, member="interior_main_beam")
model.set_member(longitudinal_beam, member="exterior_main_beam_1")
model.set_member(longitudinal_beam, member="exterior_main_beam_2")
model.set_member(edge_longitudinal_beam, member="edge_beam")
model.set_member(transverse_slab, member="transverse_slab")
model.set_member(end_transverse_slab, member="start_edge")
model.set_member(end_transverse_slab, member="end_edge")
model.create_osp_model(pyfile=False)

[4]:

# Dead loads
beam_mag = 22.4 * kN / m
DL_self_weight = og.create_load_case(name="Super T self weight")
for z_pos in model.Mesh_obj.noz[1:-1]:
    p1 = og.create_load_vertex(x=0, z=z_pos, p=beam_mag)
    p2 = og.create_load_vertex(x=L, z=z_pos, p=beam_mag)
    DL_self_weight.add_load(og.create_load(loadtype="line", point1=p1, point2=p2, name="SW"))
model.add_load_case(DL_self_weight)

# Overlay slab
overlay = og.create_load(loadtype="patch", name="overlay",
    point1=og.create_load_vertex(x=0, z=0, p=4.32*kN/m2),
    point2=og.create_load_vertex(x=L, z=0, p=4.32*kN/m2),
    point3=og.create_load_vertex(x=L, z=w, p=4.32*kN/m2),
    point4=og.create_load_vertex(x=0, z=w, p=4.32*kN/m2))
DL_overlay = og.create_load_case(name="Overlay self weight")
DL_overlay.add_load(overlay)
model.add_load_case(DL_overlay)

# SIDL
asphalt = og.create_load(loadtype="patch", name="asphalt",
    point1=og.create_load_vertex(x=0, z=0, p=1.43*kN/m2),
    point2=og.create_load_vertex(x=L, z=0, p=1.43*kN/m2),
    point3=og.create_load_vertex(x=L, z=w, p=1.43*kN/m2),
    point4=og.create_load_vertex(x=0, z=w, p=1.43*kN/m2))
left_barrier = og.create_load(loadtype="line", name="left barrier",
    point1=og.create_load_vertex(x=0, z=0, p=6.54*kN/m),
    point2=og.create_load_vertex(x=L, z=0, p=6.54*kN/m))
right_barrier = og.create_load(loadtype="line", name="right barrier",
    point1=og.create_load_vertex(x=0, z=w, p=6.54*kN/m),
    point2=og.create_load_vertex(x=L, z=w, p=6.54*kN/m))
SIDL = og.create_load_case(name="SIDL")
SIDL.add_load(asphalt)
SIDL.add_load(left_barrier)
SIDL.add_load(right_barrier)
model.add_load_case(SIDL)

# M1600 traffic
gap = 6.25 * m
z_coord = [w/2, w/2 - 3.2*m, w/2 + 3.2*m]
alf = [1.0, 0.8, 0.4]
dla = 1.3
udl_width = 3.2 * m
udl_mag = 6*kN/m / udl_width

m1600_lc_list = [og.create_load_case(name=f"M1600 L{i+1}") for i in range(3)]
M1600_moving_loads = []
for i in range(3):
    vehicle = og.create_load_model(model_type="M1600", gap=gap).create()
    vehicle.set_global_coord(og.Point(0, 0, z_coord[i]))
    m1600_lc_list[i].add_load(vehicle, load_factor=alf[i])
    M1600_moving_loads.append(vehicle)
    udl = og.create_load(loadtype="patch", name="UDL",
        point1=og.create_load_vertex(x=-L, z=z_coord[i]-udl_width/2, p=udl_mag),
        point2=og.create_load_vertex(x=2*L, z=z_coord[i]-udl_width/2, p=udl_mag),
        point3=og.create_load_vertex(x=2*L, z=z_coord[i]+udl_width/2, p=udl_mag),
        point4=og.create_load_vertex(x=-L, z=z_coord[i]+udl_width/2, p=udl_mag))
    m1600_lc_list[i].add_load(udl, load_factor=alf[i])
    model.add_load_case(m1600_lc_list[i], load_factor=dla)

# Moving loads
path = og.create_moving_path(
    start_point=og.create_point(x=-25, y=0, z=0),
    end_point=og.Point(L, 0, 0))
for i in range(3):
    ml = og.create_moving_load(name=f"Moving M1600 L{i+1}")
    ml.set_path(path)
    ml.add_load(M1600_moving_loads[i])
    model.add_load_case(ml)

[5]:

model.analyze()

Understanding the Results Dataset#

The get_results() method returns an xarray Dataset containing multiple DataArray objects. Each DataArray is indexed along named dimensions.

[6]:

# Retrieve static load case results only (fast).
# Moving load results have hundreds of increments — we fetch those
# separately in the Envelopes section below.
static_load_cases = ["Super T self weight", "Overlay self weight", "SIDL",
                     "M1600 L1", "M1600 L2", "M1600 L3"]
results = model.get_results(load_case=static_load_cases)
results

[6]:

<xarray.Dataset> Size: 546kB
Dimensions:        (Loadcase: 6, Node: 77, Component: 30, Element: 136, Nodes: 2)
Coordinates:
  * Loadcase       (Loadcase) <U53 1kB 'Super T self weight' ... 'M1600 L3'
  * Node           (Node) int64 616B 1 2 3 4 5 6 7 8 ... 70 71 72 73 74 75 76 77
  * Component      (Component) <U9 1kB 'Mx_i' 'Mx_j' 'My_i' ... 'x' 'y' 'z'
  * Element        (Element) int64 1kB 1 2 3 4 5 6 7 ... 131 132 133 134 135 136
  * Nodes          (Nodes) <U1 8B 'i' 'j'
Data variables:
    displacements  (Loadcase, Node, Component) object 111kB nan nan ... 0.0
    velocity       (Loadcase, Node, Component) object 111kB nan nan ... nan nan
    acceleration   (Loadcase, Node, Component) object 111kB nan nan ... nan nan
    forces         (Loadcase, Element, Component) object 196kB -49532.6262150...
    ele_nodes      (Loadcase, Element, Nodes) object 13kB 1 2 2 3 ... 13 77 14

The Dataset contains several DataArrays, each indexed along named dimensions:

  • ``displacements`` — (Loadcase, Node, Component) with translations (x, y, z), rotations (theta_x, …), velocities (dx, dy, dz), and accelerations (ddx, …)

  • ``forces`` — (Loadcase, Element, Component) with i-end and j-end forces (Vx_i, Vy_i, …, Mz_j)

  • ``ele_nodes`` — maps each Element tag to its i and j node tags

Above we retrieved only the six static load cases. Moving loads generate hundreds of incremental positions — we’ll fetch those later with load_case="Moving M1600 L2" when we need them for envelopes.

[7]:

# Summarise the static load cases in the results
loadcases = list(results.coords["Loadcase"].values)
print(f"Static load cases ({len(loadcases)}):")
for lc in loadcases:
    print(f"  - {lc}")

# Show the available component names (cast to plain str for tidy output)
disp_comps = [str(c) for c in results.displacements.coords["Component"].values]
force_comps = [str(c) for c in results.forces.coords["Component"].values]

# Displacement: translation + rotation + velocity + acceleration
translations = [c for c in disp_comps if c in ("x", "y", "z")]
rotations    = [c for c in disp_comps if c.startswith("theta")]
print(f"\nDisplacement components ({len(disp_comps)} total):")
print(f"  Translations: {translations}")
print(f"  Rotations:    {rotations}")

# Forces: i-end and j-end
force_i = [c for c in force_comps if c.endswith("_i")]
force_j = [c for c in force_comps if c.endswith("_j")]
print(f"\nForce components ({len(force_comps)} total):")
print(f"  i-end: {force_i}")
print(f"  j-end: {force_j}")

Static load cases (6):
  - Super T self weight
  - Overlay self weight
  - SIDL
  - M1600 L1
  - M1600 L2
  - M1600 L3

Displacement components (30 total):
  Translations: ['x', 'y', 'z']
  Rotations:    ['theta_x', 'theta_y', 'theta_z']

Force components (30 total):
  i-end: ['Mx_i', 'My_i', 'Mz_i', 'Vx_i', 'Vy_i', 'Vz_i']
  j-end: ['Mx_j', 'My_j', 'Mz_j', 'Vx_j', 'Vy_j', 'Vz_j']

[8]:

# Access a specific DataArray
displacements = results.displacements
forces = results.forces
print(f"Displacements shape: {displacements.shape}")
print(f"Forces shape: {forces.shape}")

Displacements shape: (6, 77, 30)
Forces shape: (6, 136, 30)

Selecting Results with .sel()#

xarray’s .sel() method is the primary way to extract specific results. You can select along any dimension.

[9]:

# Select vertical displacement for a specific load case
dy_sw = results.displacements.sel(Loadcase="Super T self weight", Component="y")
print(f"Self-weight deflections (first 5 nodes): {dy_sw.values[:5]}")

Self-weight deflections (first 5 nodes): [0.00035442798327501944 0.0 0.0 0.0 0.0]

[10]:

# Select bending moment for specific elements
member_name = "exterior_main_beam_1"
ext_beam_elements = model.get_element(member=member_name, options="elements")
ext_beam_nodes = model.get_element(member=member_name, options="nodes")

mz = results.forces.sel(Loadcase="M1600 L1", Component="Mz_i", Element=ext_beam_elements)
print(f"M1600 L1 bending moments on Beam 1 (kNm):\n{mz.values / 1000}")

M1600 L1 bending moments on Beam 1 (kNm):
[2.5831676077568555 314.3865886359853 671.1786019300126 1009.9015807302987
 1213.742805394697 1247.0461404027237 1117.6406533931802 852.967891701035
 511.0095849386099 203.6161309343516]

[11]:

# Select multiple load cases at once
multi_lc = results.forces.sel(
    Loadcase=["M1600 L1", "M1600 L2", "M1600 L3"],
    Component="Mz_i",
    Element=ext_beam_elements,
)
print(f"Shape: {multi_lc.shape}  (3 load cases x {len(ext_beam_elements)} elements)")

Shape: (3, 10)  (3 load cases x 10 elements)

Moving Load Results#

Moving load analyses generate one load case per incremental position. Use the load_case keyword to retrieve a specific moving load — this avoids fetching the entire (large) dataset.

[12]:

# Moving load results for Lane 2
moving_results = model.get_results(load_case="Moving M1600 L2")
print(f"Moving M1600 L2 has {len(moving_results.coords['Loadcase'])} incremental positions")

# Peek at first and last increment names
lc_names = list(moving_results.coords["Loadcase"].values)
print(f"  First: {lc_names[0]}")
print(f"  Last:  {lc_names[-1]}")

Moving M1600 L2 has 50 incremental positions
  First: Moving M1600 L2 at global position [-25.00,0.00,0.00]
  Last:  Moving M1600 L2 at global position [33.50,0.00,0.00]

Load Combinations#

Pass a dictionary of {load_case_name: factor} to compute factored combinations on the fly.

[13]:

combination = {
    "Super T self weight": 1.0,
    "M1600 L1": 1.0,
    "M1600 L2": 1.0,
    "M1600 L3": 1.0,
    "SIDL": 1.3,
    "Overlay self weight": 1.0,
}
comb_results = model.get_results(combinations=combination)
comb_results

[13]:

<xarray.Dataset> Size: 93kB
Dimensions:        (Node: 77, Component: 30, Element: 136, Nodes: 2)
Coordinates:
  * Node           (Node) int64 616B 1 2 3 4 5 6 7 8 ... 70 71 72 73 74 75 76 77
  * Component      (Component) <U9 1kB 'Mx_i' 'Mx_j' 'My_i' ... 'x' 'y' 'z'
  * Element        (Element) int64 1kB 1 2 3 4 5 6 7 ... 131 132 133 134 135 136
  * Nodes          (Nodes) <U1 8B 'i' 'j'
Data variables:
    displacements  (Node, Component) object 18kB nan nan ... 0.0
    velocity       (Node, Component) object 18kB nan nan nan nan ... nan nan nan
    acceleration   (Node, Component) object 18kB nan nan nan nan ... nan nan nan
    forces         (Element, Component) object 33kB -190594.99977261 ... nan
    ele_nodes      (Element, Nodes) object 2kB 6.3 12.6 12.6 ... 81.9 485.1 88.2
[14]:

# Extract maximum combined bending moment
comb_mz = comb_results.forces.sel(Component="Mz_i", Element=ext_beam_elements)
print(f"Max combined BM on Beam 1 = {max(comb_mz.values / 1000):.2f} kN m")

Max combined BM on Beam 1 = 8056.09 kN m

Envelopes#

The create_envelope() function finds the maximum (or minimum) response across all load cases for each element. This is essential for design — it tells you the worst-case response regardless of which load case caused it.

[15]:

# Envelope of static load cases
envelope = og.create_envelope(ds=results, load_effect="Mz_i", array="forces")
bending_env = envelope.get()
bending_env

[15]:

<xarray.DataArray 'forces' (Element: 136, Component: 30)> Size: 33kB
array([[-5766.9712876072845, -1301.6871548769068, -0.0, ..., nan, nan,
        nan],
       [19379.122833009842, 20228.82109868856, -0.0, ..., nan, nan, nan],
       [25026.028246444723, 28974.427908193567, -0.0, ..., nan, nan, nan],
       ...,
       [49399.213424556474, 40113.33714174717, 0.0, ..., nan, nan, nan],
       [57392.499131045755, -11412.184752349387, 0.0, ..., nan, nan, nan],
       [49532.54937515611, -5732.40598089605, 0.0, ..., nan, nan, nan]],
      shape=(136, 30), dtype=object)
Coordinates:
  * Element    (Element) int64 1kB 1 2 3 4 5 6 7 ... 130 131 132 133 134 135 136
  * Component  (Component) <U9 1kB 'Mx_i' 'Mx_j' 'My_i' 'My_j' ... 'x' 'y' 'z'
[16]:

# Envelope of moving load positions
moving_envelope = og.create_envelope(ds=moving_results, load_effect="Mz_i", array="forces")
moving_env = moving_envelope.get()
moving_env.sel(Element=ext_beam_elements, Component="Mz_i")

[16]:

<xarray.DataArray 'forces' (Element: 10)> Size: 80B
array([2.773708512602406e-11, 1110897.064136581, 1320696.3477219746,
       1562899.1039277394, 1754332.398995135, 1793432.9329442543,
       1776146.0418655646, 1616824.7536388119, 1386801.6559747239,
       1192240.6062729026], dtype=object)
Coordinates:
  * Element    (Element) int64 80B 20 33 46 59 72 85 98 111 124 131
    Component  <U9 36B 'Mz_i'

Query Mode — Which Load Case Governs?#

Set query_mode=True to find out which load case produced the maximum at each element. The result is a DataArray whose values are load case name strings — one per element — telling you which position governs.

[17]:

query_envelope = og.create_envelope(
    ds=moving_results, load_effect="Mz_i", array="forces", query_mode=True
)
query_env = query_envelope.get()

# Which moving load position gives max BM on each element of Beam 1?
governing = query_env.sel(Element=ext_beam_elements, Component="Mz_i")
print("Governing load case per element:")
for elem, lc_name in zip(governing.Element.values, governing.values):
    print(f"  Element {elem}: {lc_name}")

Governing load case per element:
  Element 20: Moving M1600 L2 at global position [33.50,0.00,0.00]
  Element 33: Moving M1600 L2 at global position [2.46,0.00,0.00]
  Element 46: Moving M1600 L2 at global position [4.85,0.00,0.00]
  Element 59: Moving M1600 L2 at global position [3.65,0.00,0.00]
  Element 72: Moving M1600 L2 at global position [6.04,0.00,0.00]
  Element 85: Moving M1600 L2 at global position [4.85,0.00,0.00]
  Element 98: Moving M1600 L2 at global position [7.23,0.00,0.00]
  Element 111: Moving M1600 L2 at global position [8.43,0.00,0.00]
  Element 124: Moving M1600 L2 at global position [9.62,0.00,0.00]
  Element 131: Moving M1600 L2 at global position [10.82,0.00,0.00]

Custom Post-Processing#

Extract numpy arrays for custom calculations.

[18]:

# Support nodes sit on the bridge edges (edge_node_recorder maps node tag -> edge group)
support_nodes = list(model.Mesh_obj.edge_node_recorder.keys())
print(f"Support node tags ({len(support_nodes)}): {support_nodes}")

# Vertical displacement at supports under self-weight (should be near zero)
dy_supports = results.displacements.sel(
    Loadcase="Super T self weight", Component="y", Node=support_nodes
)
print(f"Support deflections: {dy_supports.values}")

Support node tags (14): [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
Support deflections: [0.00035442798327501944 0.0 0.0 0.0 0.0 0.0 0.00035441352875530157
 0.00035442798327510477 0.0 0.0 0.0 0.0 0.0 0.00035441352875520356]

[19]:

# Compare maximum bending across all beams
for member in ["exterior_main_beam_1", "interior_main_beam", "exterior_main_beam_2"]:
    elements = model.get_element(member=member, options="elements")
    mz = results.forces.sel(Loadcase="M1600 L1", Component="Mz_i", Element=elements)
    print(f"{member}: max Mz = {max(mz.values) / 1000:.1f} kN m")

exterior_main_beam_1: max Mz = 1247.0 kN m
interior_main_beam: max Mz = 1676.3 kN m
exterior_main_beam_2: max Mz = 1247.0 kN m

Filtering Member Groups in Plots#

The convenience plotting functions (plot_bmd, plot_sfd, plot_def) accept a member parameter that controls which member groups are plotted. You can pass:

  • A string — a single member name (e.g. "interior_main_beam")

  • An og.Members bitflag — one or more groups combined with |

  • None (default) — all members for plotly, longitudinal only for matplotlib

The og.Members enum provides individual flags for each member type and pre-defined composites for common selections:

[20]:

# Available flags
print("Individual flags:")
for flag in og.Members:
    print(f"  og.Members.{flag.name}")

print("\nComposites:")
print(f"  og.Members.LONGITUDINAL = {og.Members.LONGITUDINAL}")
print(f"  og.Members.TRANSVERSE   = {og.Members.TRANSVERSE}")
print(f"  og.Members.ALL          = {og.Members.ALL}")

Individual flags:
  og.Members.EDGE_BEAM
  og.Members.EXTERIOR_MAIN_BEAM_1
  og.Members.INTERIOR_MAIN_BEAM
  og.Members.EXTERIOR_MAIN_BEAM_2
  og.Members.START_EDGE
  og.Members.END_EDGE
  og.Members.TRANSVERSE_SLAB

Composites:
  og.Members.LONGITUDINAL = Members.LONGITUDINAL
  og.Members.TRANSVERSE   = Members.TRANSVERSE
  og.Members.ALL          = Members.ALL

[21]:

# Plot BMD for longitudinal members only
og.plot_bmd(model, results, members=og.Members.LONGITUDINAL,
            loadcase="M1600 L1", backend="plotly", figsize=(14, 8));

[22]:

# Combine individual flags to plot just edge beam and interior main beam
og.plot_sfd(model, results,
            members=og.Members.EDGE_BEAM | og.Members.INTERIOR_MAIN_BEAM,
            loadcase="M1600 L1", backend="plotly", figsize=(14, 8));

[23]:

# Plot deflection for transverse members only
og.plot_def(model, results,
            loadcase="M1600 L1", backend="plotly", figsize=(14, 8));

Saving and Loading Results#

The save_filename keyword writes the Dataset to a NetCDF file — useful for archiving results or sharing them without re-running the analysis. The saved file is self-contained: it includes node coordinates and member-element connectivity, so you can plot from it without the original model object.

[ ]:

# Save results to NetCDF
results.to_netcdf("super_t_results.nc")
print("Saved to super_t_results.nc")

# Reload from disk — no model object needed
import xarray as xr
reloaded = xr.open_dataset("super_t_results.nc")
print(f"Reloaded dataset: {len(reloaded.data_vars)} variables, "
      f"{len(reloaded.coords['Loadcase'])} load cases")

# Create a lightweight proxy for plotting
proxy = og.model_proxy_from_results(reloaded)

# Verify round-trip fidelity
np.testing.assert_array_equal(
    results.forces.values, reloaded.forces.values
)
print("Round-trip verified: saved and reloaded forces match exactly.")

[ ]:

# Plot directly from the reloaded file — no original model needed
og.plot_bmd(proxy, reloaded, members=og.Members.LONGITUDINAL,
            loadcase="M1600 L1", backend="plotly", figsize=(14, 8));

[ ]:

# Clean up temporary file
import os
os.remove("super_t_results.nc")
print("Cleaned up temporary file.")

Shell-Beam Results and plot_srf#

The examples above use a beam-only model. For shell_beam models the slab is modelled with shell elements and the results Dataset contains additional DataArrays:

  • ``forces_shell`` — (Loadcase, Element, Component) with shell element nodal forces at 4 corner nodes (Vx_i…Vz_l, Mx_i…Mz_l) plus nodal displacements.

  • ``stresses_shell`` — (Loadcase, Element, Stress) with 8 stress resultants at 4 Gauss points per element (32 values): N11, N22, N12, M11, M22, M12, Q13, Q23.

plot_srf() renders contour plots over the shell mesh for shell forces (Mx, My, Mz, Vx, Vy, Vz), nodal displacements (Dx, Dy, Dz), and stress resultants (N11–Q23).

[ ]:

# Build a small shell-beam model
concrete_shell = og.create_material(
    material="concrete", code="AS5100-2017", grade="50MPa", rho=2400
)
shell_sec = og.create_section(h=0.2)
slab_shell = og.create_member(section=shell_sec, material=concrete_shell)

# Reuse the beam section from earlier
shell_long_sec = og.create_section(
    A=0.034, J=2.08e-3, Iz=6.77e-3, Iy=2.04e-3, Az=6.10e-3, Ay=3.99e-3,
)
shell_long_beam = og.create_member(section=shell_long_sec, material=concrete)

Ls, ws = 10.0 * m, 7.0 * m
shell_model = og.create_grillage(
    bridge_name="Shell Beam Demo", long_dim=Ls, width=ws, skew=0,
    num_long_grid=5, num_trans_grid=11, edge_beam_dist=1.0 * m,
    mesh_type="Orth", model_type="shell_beam",
    max_mesh_size_z=1.0, max_mesh_size_x=1.0,
    offset_beam_y_dist=0.499, beam_width=0.89,
)
shell_model.set_member(shell_long_beam, member="interior_main_beam")
shell_model.set_member(shell_long_beam, member="exterior_main_beam_1")
shell_model.set_member(shell_long_beam, member="exterior_main_beam_2")
shell_model.set_shell_members(slab_shell)
shell_model.create_osp_model(pyfile=False)

# Apply a patch load
patch = og.create_load(
    loadtype="patch", name="UDL",
    point1=og.create_load_vertex(x=0, z=0, p=10 * kN / m**2),
    point2=og.create_load_vertex(x=Ls, z=0, p=10 * kN / m**2),
    point3=og.create_load_vertex(x=Ls, z=ws, p=10 * kN / m**2),
    point4=og.create_load_vertex(x=0, z=ws, p=10 * kN / m**2),
)
udl_lc = og.create_load_case(name="UDL")
udl_lc.add_load(patch)
shell_model.add_load_case(udl_lc)

shell_model.analyze()
shell_results = shell_model.get_results()
shell_results

[ ]:

# Inspect the shell-specific DataArrays
print("forces_shell:")
print(f"  shape:      {shell_results.forces_shell.shape}")
print(f"  dimensions: {shell_results.forces_shell.dims}")
print(f"  components: {list(shell_results.forces_shell.coords['Component'].values)}")

print("\nstresses_shell:")
print(f"  shape:      {shell_results.stresses_shell.shape}")
print(f"  dimensions: {shell_results.stresses_shell.dims}")

# The Stress coordinate contains 32 labels: 8 resultants x 4 Gauss points
stress_labels = list(shell_results.stresses_shell.coords["Stress"].values)
print(f"  stress labels ({len(stress_labels)}): {stress_labels[:8]} ...")
print(f"  (8 resultants x 4 Gauss points = {len(stress_labels)} values per element)")

Shell contour plots#

plot_srf supports three families of component:

Family

Components

Source DataArray

Shell forces

Vx, Vy, Vz, Mx, My, Mz

forces_shell

Displacements

Dx, Dy, Dz

displacements

Stress resultants

N11, N22, N12, M11, M22, M12, Q13, Q23

stresses_shell

 
[ ]:

# Shell bending moment contour (Mx)
og.plot_srf(shell_results, component="Mx", backend="plotly");

[ ]:

# Vertical displacement contour (Dy) — use a sequential colorscale
og.plot_srf(shell_results, component="Dy", colorscale="Viridis", backend="plotly");

[ ]:

# Stress resultant contour (M11 — plate bending about local 2-axis)
og.plot_srf(shell_results, component="M11", backend="plotly");

Composing shell contours with beam diagrams#

Pass the Plotly figure returned by plot_srf as ax= to plot_bmd, plot_sfd, plot_tmd, or plot_def to overlay beam diagrams on the shell contour.

[ ]:

# Shell Mx contour with BMD overlay
fig = og.plot_srf(shell_results, component="Mx", backend="plotly", show=False)
og.plot_bmd(shell_model, shell_results, ax=fig, backend="plotly");

Summary#

Key techniques covered:

Technique

Method

Get static results

model.get_results(load_case=[...])

Get moving results

model.get_results(load_case="Moving M1600 L2")

Select specific results

results.forces.sel(Loadcase=..., Component=..., Element=...)

Load combinations

model.get_results(combinations={...})

Envelopes

og.create_envelope(ds=..., load_effect=..., array=...)

Query governing case

og.create_envelope(..., query_mode=True)

Extract numpy arrays

.values on any DataArray

Filter member groups

og.plot_bmd(..., members=og.Members.LONGITUDINAL)

Save results to file

model.get_results(save_filename="out.nc")

Reload from file

xr.open_dataset("out.nc")

Plot from saved file

proxy = og.model_proxy_from_results(ds) then og.plot_bmd(proxy, ds)

Shell contour plot

og.plot_srf(results, component="Mx")

Stress resultant contour

og.plot_srf(results, component="M11")

Compose contour + beam

fig = og.plot_srf(..., show=False) then og.plot_bmd(..., ax=fig)