8.2  Systematic Truss Analysis

Element Stiffness Matrix in Local Coordinates

Now, consider some arbitrary rod in a two dimensional setting.

Each node can be translated in the x- and y-direction and each node can carry loads in x- and y-directions.

To formulate the element equations systematically, we introduce the concept of the element stiffness matrix. In local coordinates, each element has four degrees of freedom: two displacement components at each end node. However, for a rod that carries only axial loads, the transverse components (perpendicular to the rod axis) do not participate in the force-displacement relationship.

Each rod has a local coordinate system \(x'\) and \(y'\) and local node numbers 1 and 2 . The local coordinate system is rotated compared to the global coordinate system, and the angle is denoted by \(\theta\). In the local coordinates we have

\[ \begin{cases}k\left(u_1'-u_3'\right) & =f_1' \\ 0 & =f_2^{\prime} \\ k\left(u_3^{\prime}-u_1^{\prime}\right) & =f_3^{\prime} \\ 0 & =f_4^{\prime}\end{cases} \]

where indices 1 and 2 refer to the first node (axial and transverse directions), and indices 3 and 4 refer to the second node. The above relation can be written on matrix form:

\[ k\left[\begin{array}{cccc} 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{array}\right]\left[\begin{array}{c} u_1^{\prime} \\ u_2^{\prime} \\ u_3^{\prime} \\ u_4^{\prime} \end{array}\right]=\left[\begin{array}{c} f_1^{\prime} \\ f_2^{\prime} \\ f_3^{\prime} \\ f_4^{\prime} \end{array}\right] \]

or

\[ \boxed{\mathbf{K}'_e \mathbf{u}' = \mathbf{f}'} \tag{8.2.1}\]

where the local element stiffness matrix is:

\[ \mathbf{K}'_e = k \begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \]

Note that this matrix is singular (has zero determinant) because two degrees of freedom are unconstrained. This is expected for a rod element that has no stiffness in the transverse direction.

The Transformation Matrix

The displacement components in global coordinates \((u, v)\) relate to local coordinates \((u', v')\) through the rotation transformation:

\[ \begin{bmatrix} u' \\ v' \end{bmatrix} = \begin{bmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix} \]

This rotation matrix transforms displacement vectors from global to local coordinates. For the complete element with two nodes (four degrees of freedom), we construct the transformation matrix \(\mathbf{L}\) by applying the rotation to each node:

\[ \boxed{\mathbf{L} = \begin{bmatrix} \cos\theta & \sin\theta & 0 & 0 \\ -\sin\theta & \cos\theta & 0 & 0 \\ 0 & 0 & \cos\theta & \sin\theta \\ 0 & 0 & -\sin\theta & \cos\theta \end{bmatrix}} \]

An important property of this matrix is that it is orthogonal, meaning:

\[ \mathbf{L}^\mathsf{T} \mathbf{L} = \mathbf{L} \mathbf{L}^\mathsf{T} = \mathbf{I} \quad \text{and} \quad \mathbf{L}^{-1} = \mathbf{L}^\mathsf{T} \]

This orthogonality has profound physical significance. The strain energy stored in a deformed element is given by \(U = \frac{1}{2}\mathbf{u}^\mathsf{T}\mathbf{K}\mathbf{u}\). Since a rotation is merely a change in perspective (how we describe the same physical configuration), the energy must remain unchanged regardless of which coordinate system we use. Mathematically, the orthogonality of \(\mathbf{L}\) ensures that the quadratic form \(\mathbf{u}^\mathsf{T}\mathbf{K}\mathbf{u}\) is preserved under the transformation, meaning the calculated strain energy is identical whether computed in local or global coordinates. This consistency is essential: if the energy changed with coordinate system choice, our analysis would yield different results depending on an arbitrary modeling decision rather than the underlying physics.

Using the transformation matrix we can transform local quantities to global

\[ \mathbf{u}_e^{\prime}=\mathbf{L}^T \mathbf{u}_e^{\prime} \Leftrightarrow \mathbf{u}_e^{\prime}=\mathbf{L} \mathbf{u}_e \tag{8.2.2}\]

and

\[ \mathbf{f}_e=\mathbf{L}^T \mathbf{f}_e^{\prime} \Leftrightarrow \mathbf{f}_e^{\prime}=\mathbf{L} \mathbf{f}_e \tag{8.2.3}\]

Inserting (8.2.2) and (8.2.3) into the global equilibrium equation (8.2.1) gives

\[ \mathbf{L}^T \mathbf{K}_e^{\prime} \mathbf{L} \mathbf{u}_e=\mathbf{f}_e \]

or

\[ \mathbf{S}_e \mathbf{u}_e=\mathbf{f}_e \]

where \(\mathbf{S}_e\) is known as the element stiffness matrix in global coordinates.

\[ \boxed{\mathbf{S}_e = \mathbf{L}^\mathsf{T} \mathbf{K}'_e \mathbf{L}} \]

Element Deformation

The deformation of an element is the change in its length, which corresponds to the relative displacement of its end nodes in the axial direction. From the figure we can see that the deformation of the rod is

\[ \delta=\Delta u=\left(u_j^x \cos \theta+u_j^y \sin \theta\right)-\left(u_i^x \cos \theta+u_i^y \sin \theta\right) \]

This can be written om matrix form as

\[ \delta=\underbrace{\left[\begin{array}{llll} -c & -s & c & s \end{array}\right]}_{\mathbf{T}}\left[\begin{array}{c} u_i^x \\ u_i^y \\ u_j^x \\ u_j^y \end{array}\right]=\mathbf{T} \cdot \mathbf{u}_e \]

where \(\mathbf T\) is known as the transformation vector that projects the nodal displacements onto the rod axis direction. We thus have

\[ \boxed{\delta = \mathbf{T} \cdot \mathbf{u}_e} \]

Member Force

The axial force in the member follows directly from the deformation:

\[ \boxed{N = k \delta = \frac{EA}{L} \mathbf{T} \cdot \mathbf{u}_e} \]

A positive value indicates tension, while a negative value indicates compression.

Strain and Stress

The engineering strain is the deformation divided by the original length:

\[ \boxed{\varepsilon = \frac{\delta}{L}} \]

Note that the length of the element \(L\) can be computed from the nodal coordinates as

\[ L = |\mathbf r| = | \mathbf p_j - \mathbf p_i | \]

The axial stress follows from Hooke’s law or directly from the force:

\[ \boxed{\sigma = \frac{N}{A} = E\varepsilon} \]

These relationships allow us to assess whether the structure remains within acceptable stress limits and to understand the load paths through the truss.

Element load vector (or reaction force vector)

\[ \mathbf{f}^e=\mathbf{T}^T N=k \mathbf{T}_{4 \times 1}^T \mathbf{T}_{1 \times 4} \mathbf{u}_{4 \times 1}^e=k\left[\begin{array}{c} -c \\ -s \\ c \\ s \end{array}\right]\left[\begin{array}{llll} -c & -s & c & s \end{array}\right]\left[\begin{array}{c} u_i^x \\ u_i^y \\ u_j^x \\ u_j^y \end{array}\right] \]

or

\[ \left[\begin{array}{l} f_i^x \\ f_i^y \\ f_j^x \\ f_j^y \end{array}\right]=\underbrace{\frac{E A}{L}\left[\begin{array}{cccc} c^2 & c s & -c^2 & -c s \\ c s & s^2 & -c s & -s^2 \\ -c^2 & -c s & c^2 & c s \\ c s & -s^2 & c s & s^2 \end{array}\right]}_{\mathbf{S}_e}\left[\begin{array}{c} u_i^x \\ u_i^y \\ u_j^x \\ u_j^y \end{array}\right] \]

Symbolic Derivation of Element Stiffness

Let us first verify the element stiffness matrix transformation symbolically. This provides confidence in our implementation and reveals the mathematical structure.

import sympy as sp

# Define symbolic variables
theta, k = sp.symbols("theta k", real=True, positive=True)
c, s = sp.cos(theta), sp.sin(theta)

# Local element stiffness matrix (template)
K_prime = k * sp.Matrix([[1, 0, -1, 0], [0, 0, 0, 0], [-1, 0, 1, 0], [0, 0, 0, 0]])

# Transformation matrix (L^T maps global to local)
L_T = sp.Matrix([[c, -s, 0, 0], [s, c, 0, 0], [0, 0, c, -s], [0, 0, s, c]])

Let’s verify the orthogonality property of the transformation matrix:

# Verify L^T * L = I
identity_check = sp.trigsimp(L_T * L_T.T)
identity_check

As expected, we obtain the identity matrix, confirming orthogonality. Now we compute the global element stiffness:

# Global element stiffness: S_e = L^T * K'_e * L
S_e = sp.trigsimp(L_T.T * K_prime * L_T)

Verification at Special Angles

To build confidence in our derivation, let’s evaluate the global stiffness matrix at specific angles where we can verify the result by inspection.

# At theta = 0 (horizontal rod)
S_e_0 = S_e.subs(theta, 0)
S_e_0

For \(\theta = 0\), the rod is horizontal, so the axial direction aligns with the global \(x\)-axis. The stiffness matrix shows coupling only in the \(x\)-direction (indices 0 and 2), which is exactly what we expect.

# At theta = 90 degrees (vertical rod)
S_e_90 = S_e.subs(theta, sp.pi / 2)
S_e_90

For \(\theta = 90°\), the rod is vertical, and the stiffness couples only the \(y\)-direction displacements (indices 1 and 3). This confirms our transformation is working correctly.

Step by step implementation

We begin with some preliminaries and initialize the global system matrix and vectors using the number of degrees of freedom.

nele = elements.shape[0]  # number of elements
nnod = nodes.shape[0]  # number of nodes
ndofs = nnod * 2  # 2 DOF per node (u_x, u_y)

S = np.zeros((ndofs, ndofs))  # Global stiffness matrix
f = np.zeros(ndofs)  # Global force vector
u = np.zeros(ndofs)  # Global displacement vector
Sprime = np.array([[1, 0, -1, 0], [0, 0, 0, 0], [-1, 0, 1, 0], [0, 0, 0, 0]])
Ltot = 0  # Total length of truss
E = 210000 * np.ones(nele)  # Young's modulus in MPa
A = 4 * 6 * np.ones(nele)  # Cross-sectional area in mm^2

print(f"Number of elements: {nele}")
print(f"Number of nodes: {nnod}")

Number of elements: 5
Number of nodes: 4

Element1

We begin by computing which global degrees of freedom correspond to our first element.

iel = 0  # Element index (0-based)
n1, n2 = elements[iel]  # Nodes for this element

The node numbers correspond to degrees of freedom, for our problem we have two, i.e., \(x-\) and \(y-\) translation. The data structure of dofs is the same as for the displacements and loads above.

# DOF indices for the first element
dofs = np.array([2 * n1, 2 * n1 + 1, 2 * n2, 2 * n2 + 1])

Computing the element stiffness matrix

p1 = nodes[n1]  # node coordinates
p2 = nodes[n2]  # node coordinates

r = p2 - p1  # element vector

L = np.linalg.norm(r)  # element length

er = r / L  # unit vector along element

c, s = er[0], er[1]  # direction cosines

R = np.array([[c, -s], [s, c]])

LT = np.block([[R, np.zeros((2, 2))], [np.zeros((2, 2)), R]])

Ai = A[iel]
Ei = E[iel]
Ki = Ei * Ai / L
print(f"Element {iel} stiffness: Ki = {Ki} MPa")
print(f"Element {iel} area: Ai = {Ai} mm^2")
print(f"Element {iel} Young's modulus: Ei = {Ei} MPa")
print(f"Element {iel} length: L = {L} mm")

Element 0 stiffness: Ki = 10080.0 MPa
Element 0 area: Ai = 24.0 mm^2
Element 0 Young's modulus: Ei = 210000.0 MPa
Element 0 length: L = 500.0 mm

Se = Ki * LT * Sprime * LT.T

Assembly process

The element stiffness matrix \(\mathbf S_e\) needs to be added to the global stiffness matrix \(\mathbf S\) at the correct locations corresponding to the global degrees of freedom of the element. This process is known as assembly. We append (add) the stiffness contributions using the following code:

np.ix_(dofs, dofs) creates a tuple of index arrays that can be used to index into a 2D array. This is particularly useful for selecting submatrices from a larger matrix based on the rows dofs and columns dofs. We take the rows and columns of the global stiffness matrix S that correspond to the degrees of freedom of the current element and add the local element stiffness matrix Se to this submatrix. This effectively assembles the contributions of the current element into the global stiffness matrix.

S[np.ix_(dofs, dofs)] += Se

This is the simplest form of assembly in code. It is easy to read and understand. It is however not very efficient for very large problems with hundreds of thousands of elements. So there exists other methods of dealing with this. But that is a topic for a course in high performance computing.

The rest of the elements are dealt with in an automated manner by executing the above lines of code for each element number, iel. This is done using the for loop:

nele = elements.shape[0]  # number of elements
nnod = nodes.shape[0]  # number of nodes
dofs = nnod * 2  # 2 DOF per node (u_x, u_y)

S = np.zeros((dofs, dofs))  # Global stiffness matrix
f = np.zeros(dofs)  # Global force vector
u = np.zeros(dofs)  # Global displacement vector
Sprime = np.array([[1, 0, -1, 0], [0, 0, 0, 0], [-1, 0, 1, 0], [0, 0, 0, 0]])
Ltot = 0  # Total length of truss
E = 210000 * np.ones(nele)  # Young's modulus in MPa
A = 4 * 6 * np.ones(nele)  # Cross-sectional area in mm^2


for iel, (n1, n2) in enumerate(elements):
    r = nodes[n2] - nodes[n1]
    L = np.linalg.norm(r)
    Ltot += L
    c, s = r / L
    R = np.array([[c, -s], [s, c]])
    LT = np.block([[R, np.zeros((2, 2))], [np.zeros((2, 2)), R]])
    Ai = A[iel]
    Ei = E[iel]
    Ki = Ei * Ai / L
    Se = Ki * LT @ Sprime @ LT.T
    dofs = np.array([2 * n1, 2 * n1 + 1, 2 * n2, 2 * n2 + 1])
    # print(Se)
    S[np.ix_(dofs, dofs)] += Se

This code block is the cornerstone of the systematic truss analysis. It automates the assembly of the global stiffness matrix by iterating over each element, computing its local stiffness matrix, transforming it to global coordinates, and adding it to the appropriate locations in the global matrix. This loop encapsulates the essence of the direct stiffness method, allowing us to handle complex truss structures with ease.

Applying Boundary Conditions and Solving the System

The resulting system \(\mathbf {S u = f}\), however, is singular, meaning it cannot be solved directly. In mechanical terms, the structure contains rigid body rotations or translations that are not constrained by supports. To resolve this, we apply boundary conditions by modifying the global stiffness matrix and load vector to account for fixed supports. This typically involves removing rows and columns corresponding to constrained degrees of freedom.

Note that special treatment of the right-hand side is required in cases where the displacements are non-zero. See [TODO] for details.

We assign the known boundary conditions and loads. Take a look at the figure below for reference.

presc = np.array([0, 1, 3])  # Prescribed DOFs (0-based)
u[presc] = 0.0  # Prescribed displacements in mm

# Loads
f[2 * 4 - 1] = -10000.0  # Load in N at node 4 in negative y-direction

free = np.setdiff1d(np.arange(nnod * 2), presc)

S[np.ix_(free, free)]


u_free = np.linalg.solve(S[np.ix_(free, free)], f[free])
u[free] = u_free

We can separate the x- and y-components of the displacements into a two column matrix for easier handling in the visualization step below.

U = u.reshape((nnod, 2))

Visualization of Results

Visualization is essential for understanding structural behavior. We plot the deformed shape (with amplification for visibility) and color-code the members according to their stress levels.

UR = np.sum(U**2, axis=1) ** 0.5

Code

plt.style.use("default")
fig, ax = plt.subplots(figsize=(10, 8))
ax.set_aspect("equal")

# Amplification factor for visualization
scale = 80

# Plot undeformed structure
for n1, n2 in elements:
    p1 = nodes[n1]
    p2 = nodes[n2]
    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], "k--", linewidth=1, alpha=0.5)

# Plot deformed structure
for n1, n2 in elements:
    p1 = nodes[n1] + scale * U[n1]
    p2 = nodes[n2] + scale * U[n2]
    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], "k-", linewidth=2)

Computing Member Quantities

Once the displacements have been computed we can compute the normal forces, resultant forces, deformations, strains and stresses.

React = np.zeros(nnod * 2)  # Reaction forces in the nodes, two DOF per node
delta = np.zeros(nele)  # Element deformations
N = delta.copy()  # Normal forces in the elements
epsilon = delta.copy()  # Strains in the elements
sigma = delta.copy()  # Stresses in the elements

for iel, (n1, n2) in enumerate(elements):
    r = nodes[n2] - nodes[n1]
    L = np.linalg.norm(r)
    c, s = r / L
    T = np.array([-c, -s, c, s])

    Ai = A[iel]
    Ei = E[iel]
    Ki = Ei * Ai / L

    dofs = np.array([2 * n1, 2 * n1 + 1, 2 * n2, 2 * n2 + 1])

    ue = u[dofs]
    delta_i = T @ ue
    N_i = Ki * delta_i
    fe = N_i * T
    epsilon_i = delta_i / L
    sigma_i = N_i / Ai

    React[dofs] += fe
    N[iel] = N_i
    delta[iel] = delta_i
    epsilon[iel] = epsilon_i
    sigma[iel] = sigma_i


# Split reactions into x and y components for clarity
Rx = React[0::2]
Ry = React[1::2]
R = np.array([Rx, Ry]).T

Visualize element results on deformed shape

# Plot the member stresses on the deformed shape
# | code-fold: true
# | label: fig-deformed-shape-stress
plt.style.use("default")
fig, ax = plt.subplots(figsize=(10, 8))
ax.set_aspect("equal")

# Amplification factor for visualization
scale = 80

# Plot undeformed structure
for n1, n2 in elements:
    p1 = nodes[n1]
    p2 = nodes[n2]
    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], "k--", linewidth=1, alpha=0.5)

# Plot deformed structure with color-coded stress
deformed_nodes = nodes + scale * u.reshape(-1, 2)
element_field = sigma  # in MPa
field = np.array(element_field)
norm = plt.Normalize(vmin=field.min(), vmax=field.max())
cmap = plt.cm.jet

for elem_idx, (n1, n2) in enumerate(elements):
    p1 = deformed_nodes[n1]
    p2 = deformed_nodes[n2]
    color = cmap(norm(element_field[elem_idx]))
    ax.plot(
        [p1[0], p2[0]], [p1[1], p2[1]], color=color, linewidth=3, solid_capstyle="round"
    )

# Plot nodes
# ax.scatter(nodes[:, 0], nodes[:, 1], c='blue', s=100, zorder=5, label='Undeformed')
ax.scatter(deformed_nodes[:, 0], deformed_nodes[:, 1], c="k", s=100, zorder=5)

# Add node labels
# for i in range(nnod):
#    ax.annotate(f'{i+1}', nodes[i] + np.array([-20, 10]), fontsize=12)

# Add element labels at midpoints of deformed elements
for i, (n1, n2) in enumerate(elements):
    mid_point = (deformed_nodes[n1] + deformed_nodes[n2]) / 2
    ax.text(
        mid_point[0],
        mid_point[1],
        str(f"$\\sigma={element_field[i]:.2f}$ MPa"),
        backgroundcolor="w",
        horizontalalignment="center",
        verticalalignment="center",
        fontsize=10,
        fontweight="normal",
    )

# Add colorbar
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, shrink=0.4)
cbar.set_label("Stress [MPa]", fontsize=12)

ax.set_xlabel("x [mm]", fontsize=12)
ax.set_ylabel("y [mm]", fontsize=12)
ax.set_title(f"Deformed Shape (Scale Factor: {scale}×)", fontsize=14)
# ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
plt.tight_layout()
plt.show()

# Plot the member force on the deformed shape
# | code-fold: true
# | label: fig-deformed-shape-force
plt.style.use("default")
fig, ax = plt.subplots(figsize=(10, 8))
ax.set_aspect("equal")

# Amplification factor for visualization
scale = 80

# Plot undeformed structure
for n1, n2 in elements:
    p1 = nodes[n1]
    p2 = nodes[n2]
    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], "k--", linewidth=1, alpha=0.5)

# Plot deformed structure with color-coded stress
deformed_nodes = nodes + scale * u.reshape(-1, 2)
element_field = N
field = np.array(element_field)
norm = plt.Normalize(vmin=field.min(), vmax=field.max())
cmap = plt.cm.jet

for elem_idx, (n1, n2) in enumerate(elements):
    p1 = deformed_nodes[n1]
    p2 = deformed_nodes[n2]
    color = cmap(norm(element_field[elem_idx]))
    ax.plot(
        [p1[0], p2[0]], [p1[1], p2[1]], color=color, linewidth=3, solid_capstyle="round"
    )

# Plot nodes
# ax.scatter(nodes[:, 0], nodes[:, 1], c='blue', s=100, zorder=5, label='Undeformed')
ax.scatter(deformed_nodes[:, 0], deformed_nodes[:, 1], c="k", s=100, zorder=5)

# Add node labels
# for i in range(nnod):
#    ax.annotate(f'{i+1}', nodes[i] + np.array([-20, 10]), fontsize=12)

# Add element labels at midpoints of deformed elements
for i, (n1, n2) in enumerate(elements):
    mid_point = (deformed_nodes[n1] + deformed_nodes[n2]) / 2
    ax.text(
        mid_point[0],
        mid_point[1],
        str(f"$N={element_field[i]:.0f}$ N"),
        backgroundcolor="w",
        horizontalalignment="center",
        verticalalignment="center",
        fontsize=10,
        fontweight="normal",
    )

# Add colorbar
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, shrink=0.4)
cbar.set_label("Element force [N]", fontsize=12)

ax.set_xlabel("x [mm]", fontsize=12)
ax.set_ylabel("y [mm]", fontsize=12)
ax.set_title(f"Deformed Shape (Scale Factor: {scale}×)", fontsize=14)
# ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
plt.tight_layout()
plt.show()

# Plot the member deformation on the deformed shape
# | code-fold: true
# | label: fig-deformed-shape-deformation
plt.style.use("default")
fig, ax = plt.subplots(figsize=(10, 8))
ax.set_aspect("equal")

# Amplification factor for visualization
scale = 80

# Plot undeformed structure
for n1, n2 in elements:
    p1 = nodes[n1]
    p2 = nodes[n2]
    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], "k--", linewidth=1, alpha=0.5)

# Plot deformed structure with color-coded stress
deformed_nodes = nodes + scale * u.reshape(-1, 2)
element_field = delta  # in mm
field = np.array(element_field)
norm = plt.Normalize(vmin=field.min(), vmax=field.max())
cmap = plt.cm.jet

for elem_idx, (n1, n2) in enumerate(elements):
    p1 = deformed_nodes[n1]
    p2 = deformed_nodes[n2]
    color = cmap(norm(element_field[elem_idx]))
    ax.plot(
        [p1[0], p2[0]], [p1[1], p2[1]], color=color, linewidth=3, solid_capstyle="round"
    )

# Plot nodes
# ax.scatter(nodes[:, 0], nodes[:, 1], c='blue', s=100, zorder=5, label='Undeformed')
ax.scatter(deformed_nodes[:, 0], deformed_nodes[:, 1], c="k", s=100, zorder=5)

# Add node labels
# for i in range(nnod):
#    ax.annotate(f'{i+1}', nodes[i] + np.array([-20, 10]), fontsize=12)

# Add element labels at midpoints of deformed elements
for i, (n1, n2) in enumerate(elements):
    mid_point = (deformed_nodes[n1] + deformed_nodes[n2]) / 2
    ax.text(
        mid_point[0],
        mid_point[1],
        str(f"$\\delta={element_field[i]:.3f}$ mm"),
        backgroundcolor="w",
        horizontalalignment="center",
        verticalalignment="center",
        fontsize=10,
        fontweight="normal",
    )

# Add colorbar
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, shrink=0.4)
cbar.set_label("Element deformation [mm]", fontsize=12)

ax.set_xlabel("x [mm]", fontsize=12)
ax.set_ylabel("y [mm]", fontsize=12)
ax.set_title(f"Deformed Shape (Scale Factor: {scale}×)", fontsize=14)
# ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
plt.tight_layout()
plt.show()

# Plot the member strain on the deformed shape
# | code-fold: true
# | label: fig-deformed-shape-strain
plt.style.use("default")
fig, ax = plt.subplots(figsize=(10, 8))
ax.set_aspect("equal")

# Amplification factor for visualization
scale = 80

# Plot undeformed structure
for n1, n2 in elements:
    p1 = nodes[n1]
    p2 = nodes[n2]
    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], "k--", linewidth=1, alpha=0.5)

# Plot deformed structure with color-coded stress
deformed_nodes = nodes + scale * u.reshape(-1, 2)
element_field = epsilon  # in mm
field = np.array(element_field)
norm = plt.Normalize(vmin=field.min(), vmax=field.max())
cmap = plt.cm.jet

for elem_idx, (n1, n2) in enumerate(elements):
    p1 = deformed_nodes[n1]
    p2 = deformed_nodes[n2]
    color = cmap(norm(element_field[elem_idx]))
    ax.plot(
        [p1[0], p2[0]], [p1[1], p2[1]], color=color, linewidth=3, solid_capstyle="round"
    )

# Plot nodes
# ax.scatter(nodes[:, 0], nodes[:, 1], c='blue', s=100, zorder=5, label='Undeformed')
ax.scatter(deformed_nodes[:, 0], deformed_nodes[:, 1], c="k", s=100, zorder=5)

# Add node labels
# for i in range(nnod):
#    ax.annotate(f'{i+1}', nodes[i] + np.array([-20, 10]), fontsize=12)

# Add element labels at midpoints of deformed elements
for i, (n1, n2) in enumerate(elements):
    mid_point = (deformed_nodes[n1] + deformed_nodes[n2]) / 2
    ax.text(
        mid_point[0],
        mid_point[1],
        str(f"$\\varepsilon={element_field[i]*100:.3f}$ %"),
        backgroundcolor="w",
        horizontalalignment="center",
        verticalalignment="center",
        fontsize=10,
        fontweight="normal",
    )

# Add colorbar
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, shrink=0.4)
cbar.set_label("Element strain [%]", fontsize=12)

ax.set_xlabel("x [mm]", fontsize=12)
ax.set_ylabel("y [mm]", fontsize=12)
ax.set_title(f"Deformed Shape (Scale Factor: {scale}×)", fontsize=14)
# ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
plt.tight_layout()
plt.show()