Plotting and Postprocessing#

This chapter covers stress recovery, reaction extraction, and visualization in femlabpy.

1. Stress Recovery Mathematics#

1.1 Gauss point to node extrapolation#

For isoparametric elements, stresses are most accurate at Gauss points. To obtain nodal values, we extrapolate using the inverse of the shape function matrix evaluated at Gauss points.

For a Q4 element with \(2 \times 2\) Gauss integration, the Gauss point coordinates are \((\pm 1/\sqrt{3}, \pm 1/\sqrt{3})\). The extrapolation uses the shape functions evaluated at corner nodes (\(\xi, \eta = \pm 1\)):

\[ \boldsymbol{\sigma}^{nodal} = \mathbf{E} \boldsymbol{\sigma}^{GP} \]

where \(\mathbf{E}\) is the \(4 \times 4\) extrapolation matrix:

\[ E_{ij} = N_i(\xi_j^{GP}, \eta_j^{GP}) \]

evaluated with corner shape functions at Gauss point locations.

1.2 Nodal averaging#

When multiple elements share a node, the extrapolated values are averaged:

\[ \sigma_i^{avg} = \frac{\sum_{e \in \text{elements at } i} \sigma_i^{(e)}}{n_{\text{elements at } i}} \]

This smoothing reduces discontinuities between elements.

1.3 Stress components#

The Voigt notation stress vector contains:

2D (plane stress/strain): $\( \boldsymbol{\sigma} = [\sigma_{xx}, \sigma_{yy}, \sigma_{xy}]^T \)$

3D: $\( \boldsymbol{\sigma} = [\sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \sigma_{xy}, \sigma_{yz}, \sigma_{zx}]^T \)$

2. Reaction Forces#

2.1 Theory#

At constrained DOFs, the reaction force equals the internal force:

\[ R_i = F_i^{int} = \sum_e \int_{\Omega^e} \mathbf{B}^T \boldsymbol{\sigma} \, dV \bigg|_{\text{DOF } i} \]

For equilibrium: \(\mathbf{R} + \mathbf{F}^{ext} = \mathbf{F}^{int}\)

At free DOFs: \(R_i = 0\) (by definition)

At constrained DOFs: \(R_i = F_i^{int} - F_i^{ext}\)

2.2 Implementation#

The reaction() function extracts reactions from the assembled internal force vector:

R = fp.reaction(q, C, dof)

Returns [node, local_dof, reaction] for each constrained DOF.

3. Mesh and Field Plotting#

The plotting helpers in plotting.py work directly on the topology and coordinate arrays used elsewhere in the package.

plotelem#

plotelem(T, X, ...) draws the undeformed mesh. It reads each element row, subtracts one from the one-based node ids, and plots the resulting polygon or polyline.

The helper also supports optional node and element labels. That is useful when you are checking topology ordering or debugging a bad mesh import.

plotforces#

plotforces(T, X, P, ...) draws nodal loads as arrows. The element table is not used directly; it is accepted only for a consistent plotting signature.

The arrow size is scaled from the mesh span so a small model and a large model remain readable with the same function call.

plotbc#

plotbc(T, X, C, ...) visualizes prescribed displacements. Zero-valued constraints are shown as markers, while nonzero prescribed values are drawn as arrows.

That distinction matches the actual solver behavior:

  • zero values are fixed degrees of freedom;

  • nonzero values are enforced displacements.

plotu#

plotu(T, X, u, dof=None, component=0, ...) draws a nodal contour field. It accepts either:

  • one scalar value per node,

  • a nodal array with several components per node, or

  • a flattened global vector with length nn * dof.

When several components are available, component=0 plots the nodal magnitude and component=1, 2, ... selects a specific component. Flat Gmsh coordinate arrays of the form (x, y, z) are reduced automatically to a 2D plotting view when all out-of-plane coordinates are constant.

The contour still works the same way internally: each element is colored by the mean of its nodal samples, using PolyCollection in 2D and Poly3DCollection in 3D.

Element Contours#

plotq4#

plotq4(T, X, S, scomp, ...) reconstructs a nodal contour from Q4 Gauss-point results. The code expects the Q4 storage layout used by the Q4 recovery functions: each component is stored at four Gauss points per element.

The implementation does three things:

  1. it extracts the requested component with one-based scomp;

  2. it maps the four Gauss-point values to corner values with a small extrapolation matrix;

  3. it averages contributions from adjacent elements at each node.

The final contour is displayed with tripcolor(..., shading="gouraud") so the result looks smooth across the mesh.

plott3#

plott3(T, X, S, scomp, ...) follows the same idea for T3 output. The recovery is simpler because the T3 storage is already element-centered in the form used by the rest of the codebase.

Reactions#

The reaction() function in postprocess.py reads the internal force vector and extracts the support reactions associated with constrained degrees of freedom.

It does not recompute the finite element solution. It assumes you already have the assembled internal force vector q and the constraint table C.

Inputs and outputs#

  • q is the global internal force vector.

  • C is the constraint table.

  • dof tells the helper how to flatten node and local DOF ids.

  • comp optionally filters the returned rows by a single constrained component.

The output is either:

  • [node, local_dof, reaction], or

  • [constraint_row, reaction] when comp is supplied.

Why this helper is small#

Support reactions are a recovery task, not a solve task. The helper exists so the solve path can stay focused on matrix assembly while the postprocess layer handles the bookkeeping required to report reactions cleanly.

Typical Workflow#

import numpy as np
import femlabpy as fp

K, p, q = fp.init(nn, dof)
K = fp.kq4e(K, T, X, G)
p = fp.setload(p, P)
K, p, _ = fp.setbc(K, p, C, dof)
u = np.linalg.solve(K, p)
q, S, E = fp.qq4e(q, T, X, G, u)
R = fp.reaction(q, C, dof)

Then use the plotting helpers to inspect the result:

  • plotelem for geometry,

  • plotforces for loads,

  • plotbc for constraints,

  • plotu for scalar nodal fields or displacement-component contours,

  • plotq4 or plott3 for stress or strain contours.