Chapter 9: Dynamic Response Workflows#

This chapter focuses on the parts of femlabpy that are easiest to misuse if you treat them as black boxes: modal output, frequency response, seismic loading, time-history arrays, and the plotting helpers that sit on top of them. The examples below follow the actual implementation in src/femlabpy/dynamics.py, src/femlabpy/modal.py, src/femlabpy/damping.py, src/femlabpy/plotting.py, and src/femlabpy/postprocess.py.

9.2 Frequency response workflow#

compute_frf() evaluates the complex dynamic stiffness

\[ \mathbf{Z}(\omega) = \mathbf{K} - \omega^2 \mathbf{M} + i \omega \mathbf{C} \]

at a set of frequencies and solves one linear system per frequency. The function returns two arrays:

  • freq_hz: the sampled frequency axis in Hz.

  • H: the complex transfer function between the selected input and output DOFs.

This is a scalar FRF, not a full matrix. The input and output DOF indices are global and 0-based.

plot_frf() then converts H into magnitude and phase plots. If mark_peaks=True, the helper uses scipy.signal.find_peaks to mark likely resonance peaks on the magnitude plot.

Example: one input, one output#

import numpy as np
import matplotlib.pyplot as plt

from femlabpy.dynamics import compute_frf, plot_frf

# Assume K, M, and C are already assembled and boundary conditions are applied.
freq_hz, H = compute_frf(
    M,
    C,
    K,
    input_dof=0,
    output_dof=0,
    freq_range=(0.1, 50.0),
    n_points=800,
)

fig = plot_frf(freq_hz, H, log_scale=True, mark_peaks=True)
plt.show()

If you only need the resonance locations, the sampled magnitude is usually enough. If you need a transfer function for later analysis, keep H and postprocess it yourself.

9.3 Seismic loading workflow#

seismic_load() is a load-function factory. It does not solve the dynamic problem by itself. Instead, it precomputes the effective base-excitation vector

\[ \mathbf{p}(t) = -\mathbf{M} \mathbf{r} a_g(t) \]

and returns a callable p(t) that interpolates the recorded ground acceleration at any requested time.

What matters in the code is that the function accepts either a dense or sparse mass matrix, multiplies it by the influence vector once, and then uses np.interp for the ground-motion history. That keeps the time integrator simple: the solver only needs to call p_func(t_next) at each step.

Example: ground-motion input#

import numpy as np
import matplotlib.pyplot as plt

from femlabpy.dynamics import seismic_load, solve_newmark, plot_time_history, plot_energy

dt_record = 0.01
n_points = 1000
t_record = np.arange(n_points) * dt_record
ag_g = 0.5 * np.sin(2.0 * np.pi * 2.5 * t_record) * np.exp(-0.2 * t_record)
ag = ag_g * 9.80665

# Reuse the same assembled model and boundary-condition table from the modal example.
C_damp = np.zeros_like(K)

# Influence vector: 1 for x-direction DOFs, 0 for all others.
influence = np.zeros(ndof, dtype=float)
influence[0::2] = 1.0

p_func = seismic_load(M, influence, ag, dt_record)

u0 = np.zeros((ndof, 1))
v0 = np.zeros((ndof, 1))

history = solve_newmark(
    M,
    C_damp,
    K,
    p_func,
    u0,
    v0,
    dt=dt_record,
    nsteps=n_points,
    C_bc=data["C"],
    dof=dof,
    compute_energy=True,
)

solve_newmark() returns a TimeHistory object. Its arrays are stored time-first:

  • history.t has shape (nsteps + 1,).

  • history.u, history.v, and history.a have shape (nsteps + 1, ndof).

  • Row 0 is the initial state.

  • Column j is the history of global DOF j.

That means a single DOF history is extracted as history.u[:, target_dof], not by reshaping the array.

Example: extracting one roof DOF#

roof_node = np.argmax(X[:, 1])
roof_dof_x = roof_node * dof
u_roof = history.u[:, roof_dof_x]

plt.figure(figsize=(10, 4))
plt.plot(history.t, u_roof * 1000.0, color="tab:blue", linewidth=1.5)
plt.xlabel("Time (s)")
plt.ylabel("Displacement (mm)")
plt.title("Roof displacement history")
plt.grid(True, alpha=0.3)
plt.show()

plot_time_history() is the library helper for the same task. It accepts one DOF index or a list of indices and plots displacement, velocity, or acceleration directly from the TimeHistory object.

9.4 Energy and postprocessing workflow#

When compute_energy=True, the solvers populate history.energy with a dictionary. In the current implementation the keys are:

  • kinetic

  • strain

  • total

plot_energy() expects exactly that dictionary. If the solver was not run with energy tracking enabled, the helper raises a ValueError instead of guessing.

Example: energy plot#

fig, ax = plt.subplots(figsize=(10, 4))
plot_energy(history, ax=ax)
plt.show()

The implementation is intentionally simple. It does not attempt to reconstruct external work or damping work histories. If you need those, you should compute them explicitly from your loading history and state vectors.

Example: using the built-in time-history plotter#

fig, ax = plt.subplots(figsize=(10, 4))
plot_time_history(history, [roof_dof_x], quantity="displacement", ax=ax)
plt.show()

This helper is useful when you want the standard library behavior and do not want to write array slicing and axis labeling by hand.