Integration methods: polygon vs fiber¶
GenSec evaluates two very different families of integrals over the cross-section \(\Omega\). The geometric quantities (area, static moments, second-moments, elastic and plastic moduli) involve integrands that are polynomial in the coordinates and independent of the deformation state, and therefore admit closed-form evaluation via Green’s theorem on a polygon. The mechanical response quantities (axial force, bending moments, moment–curvature and interaction surfaces) involve integrands of the form \(\sigma(\varepsilon(x, y))\) that are non-linear in the coordinates, typically only \(C^0\) at the branch points of the constitutive law, and dependent on the current deformation state; they admit no closed form and must be evaluated by numerical quadrature, which in this context is universally done by fiber summation.
GenSec uses each method exclusively where it is best suited. This page documents the two methods, their accuracy, their cost, and the reasons they are not mixed.
The polygonal method (Green’s theorem)¶
Formulation¶
For any polynomial integrand \(P(x, y)\) and any simply connected region \(\Omega\) bounded by a closed curve \(\partial\Omega\), Green’s theorem converts the 2-D area integral into a 1-D line integral:
with \(Q\) any anti-derivative of \(P\) with respect to one of the coordinates. When \(\partial\Omega\) is a polygon with vertices \((x_i, y_i),\; i = 1, \dots, n\) (indexed modulo \(n\)), the line integral reduces to a sum of elementary integrals over straight segments, each of which is itself a closed-form polynomial expression.
For the six moment integrals needed to compute \(A, S_x, S_y, I_{xx,O}, I_{yy,O}, I_{xy,O}\) of the polygon about the origin, one obtains
with \(c_i = x_i y_{i+1} - x_{i+1} y_i\). Rings with holes are handled by summing the integrals of every ring, using counter-clockwise orientation for the exterior and clockwise for every hole, so that the signed contribution of a hole subtracts automatically.
Accuracy¶
The formulas above are exact: for any polygon with vertices given exactly, the result is the true value of the integral, up to floating-point round-off of order \(10^{-15}\) relative for standard double-precision arithmetic.
When the section boundary is a smooth curve that has been approximated by a polygon of \(N\) vertices — for instance a circle of radius \(R\) discretised as a regular \(N\)-gon — the only error is the boundary approximation error. For a regular \(N\)-gon inscribed in a circle,
so the relative error of the area decays as \(\mathcal{O}(N^{-2})\). The same rate governs all second-moment integrals. Numerically (see the table below), 64 vertices give an error of \(\approx 1.6\times 10^{-3}\), 256 vertices give \(\approx 10^{-4}\).
Cost¶
Each moment integral is a single pass over the
\(n\)-vertex boundary, i.e. \(\mathcal{O}(n)\) in time
and space. A typical RC cross-section has
\(n \lesssim 100\) even for complex boundaries, so all
geometric properties are computed in milliseconds. Full
computation of compute_section_properties() (area,
moments, principal axes, \(W\), kern) on a 100-vertex
polygon is \(\lesssim 3\,\mathrm{ms}\) on a 2024 laptop
CPU.
The fiber method¶
Formulation¶
When the integrand depends on the deformation state through a constitutive law \(\sigma = \hat\sigma(\varepsilon)\), and the strain is a linear function of the fiber coordinates,
the resultant axial force and bending moments are integrals of the form
with \(\mathbf{e} = (\varepsilon_0, \chi_y, \chi_z)\). The integrand is not polynomial in \(y, z\); it has piecewise-smooth branches separated by the pivot points of the constitutive law (e.g. \(\varepsilon_{cr}\) for concrete cracking, \(\varepsilon_{c2}, \varepsilon_{cu2}\) for EC2 parabola–rectangle, \(\varepsilon_{sy}, \varepsilon_{su}\) for elastoplastic steel). Green’s theorem cannot be applied.
The fiber method partitions \(\Omega\) into a finite set of cells \(\{\Omega_k\}\) with centroids \((y_k, z_k)\) and areas \(A_k\), and approximates each integral by a midpoint sum:
Point fibers (rebars, tendons, FRP strips) are added analytically: they are not quadrature points but exact point contributions \(A_s\,\hat\sigma_s(\varepsilon(y_s, z_s))\). Embedded rebars additionally subtract the bulk contribution at the same location to avoid double-counting the area they physically displace, so the total contribution is \(A_s\,[\hat\sigma_s - \hat\sigma_c(\varepsilon_s)]\).
Boundary fitting — what GenSec actually does¶
Naïve fiber meshes — a rectangular \(N \times N\) grid with every cell contributing either “fully in” or “fully out” based on whether the cell centre lies inside the polygon — suffer from a boundary-fitting pathology: the sign of the local quadrature error depends on how each boundary cell happens to straddle \(\partial\Omega\), and the pattern changes erratically with \(N\). Convergence drops to \(\mathcal{O}(h)\) or worse, and can be non-monotone (see the disc table below).
GenSec avoids this completely. The grid mesher in
gensec.geometry.geometry intersects every cell with the
polygon using Shapely, and each fiber inherits the centroid
and area of the intersection:
cell = box(x0, y0, x1, y1)
clipped = polygon.intersection(cell)
if not clipped.is_empty:
x_fib, y_fib = clipped.centroid.x, clipped.centroid.y
A_fib = clipped.area
This is exact up to the polygon’s own discretisation of the boundary. Every fiber sits on a polygonal sub-domain, and the midpoint rule recovers its theoretical \(\mathcal{O}(h^2)\) rate — monotonic, without the fluctuations seen in the naïve version. The clipping is a one-time cost at section assembly: the subsequent Newton iterations in the solver operate on the pre-computed fiber arrays \((x_k, y_k, A_k)\) and never invoke Shapely.
For sections where a structured grid is inadequate (complex
boundaries, locally refined regions), GenSec also supports a
constrained Delaunay mesh via the triangle backend; the
accuracy discussion below applies identically.
Accuracy — branch points of the constitutive law¶
Even with perfect boundary fitting, the constitutive law introduces a second source of error. Where \(\hat\sigma\) is only \(C^0\) — the parabola/plateau transition at \(\varepsilon = \varepsilon_{c2}\), the yield point at \(\varepsilon = \varepsilon_{sy}\), the tension cut-off at \(\varepsilon = 0\) for concrete — the midpoint rule has a local error that degrades to \(\mathcal{O}(h)\) inside a strip of width \(\mathcal{O}(h)\) around each iso-strain line \(\varepsilon(y, z) = \varepsilon_b\). Integrating the strip over its \(\mathcal{O}(h)\) width gives an \(\mathcal{O}(h^2)\) contribution to the total error, so the global rate is preserved; the constant is larger than for a smooth integrand, but the asymptotic behaviour is still quadratic.
The RC benchmark in RC benchmark — fiber method vs exact PR integration confirms this: convergence is \(\mathcal{O}(h^2)\) with a constant that would be no worse if the integrand were polynomial.
Accuracy — coupling with the non-linear solver¶
The fiber sum is evaluated many times inside a Newton or bisection loop on \(\mathbf{e}\). Residual errors in the quadrature propagate into the equilibrium solution, so over-coarse meshes manifest as spurious oscillation or failure of the non-linear solver, not just a slightly wrong answer. In practice, a mesh size of a few millimetres on an RC column gives \(10^{-5}\)–\(10^{-4}\) relative error on \((N, M_y, M_z)\), which is more than enough: the engineering uncertainty on material constants is several orders of magnitude larger.
Why not use higher-order quadrature (Gauss, Simpson)?¶
For a smooth integrand on a rectangular domain, quadrature rules of order \(p\) achieve error \(\mathcal{O}(h^{p})\) per cell with a fixed number of evaluation points per cell, vastly outperforming midpoint (\(p = 2\)). Gauss–Legendre with \(k\) points per axis, for instance, integrates polynomials of degree \(2k-1\) exactly. One might ask why GenSec does not use such a rule.
The answer is in the constitutive law.
Higher-order rules gain nothing across branch points. Their exactness proofs assume the integrand is a polynomial — or at least \(C^p\)-smooth — over the cell. Across an iso-strain line where \(\hat\sigma\) is only \(C^0\), every quadrature rule has local error \(\mathcal{O}(h^2)\) regardless of its nominal order. The total error then integrates as \(\mathcal{O}(h^2)\) globally, the same as midpoint. Investing in Gauss points buys nothing on the strips that dominate the error, and pays extra evaluations of \(\hat\sigma\) everywhere else.
Higher-order rules can be strictly worse. When \(\hat\sigma\) is non-smooth and non-monotonic — for instance the descending branch of a concrete law past \(\varepsilon_{cu2}\), or a damage model with load reversal — Gauss rules can produce spurious cancellations between points on opposite sides of the discontinuity. Midpoint, with a single sample per cell, cannot exhibit this pathology; it systematically over- or under-estimates each cell consistently with the iso-strain structure.
Mesh refinement near branch points is the right remedy. If the user truly needs higher accuracy near a transition line, the right tool is local mesh refinement (smaller cells in the strip where \(\varepsilon \approx \varepsilon_b\)), not a higher-order rule on the coarse mesh. GenSec’s triangular mesher supports non-uniform sizing for exactly this purpose; cells can be refined to align with the expected branch-point bands.
The short answer, then, is that midpoint on a conforming mesh is the optimal rule for this class of integrands: it is robust against non-smoothness, achieves the best asymptotic rate any quadrature can achieve here, and offers the most predictable behaviour inside the non-linear solver.
Cost¶
Each evaluation of \((N, M_y, M_z)\) is \(\mathcal{O}(n_{\mathrm{fib}})\) where \(n_{\mathrm{fib}}\) is the total fiber count. A full verification — build-up of a 3-D interaction surface with \(n_{\mathrm{pts}}\) sampled points, each requiring tens of Newton iterations — therefore costs \(\mathcal{O}(n_{\mathrm{pts}} \cdot n_{\mathrm{iter}} \cdot n_{\mathrm{fib}})\) constitutive-law evaluations. See RC benchmark — fiber method vs exact PR integration for per-call timings on a representative RC section, and Engineering cost budget for the full-surface cost budget.
Why the two methods are not mixed¶
Two consistency arguments and one efficiency argument.
Consistency. All the geometric quantities \(A, S, I, W, Z\) derive from polynomial integrands on the same domain \(\Omega\). If the area were computed by fiber summation and the moments of inertia by Green’s theorem, the relation \(\rho = \sqrt{I/A}\) would mix two incommensurable approximations; the reported radius of gyration would inherit two different error contributions with different rates.
Similarly, the plastic neutral axis in \(Z\) is found by splitting \(A_{\mathrm{id}}\) in half. If \(A\) came from a fiber sum with \(10^{-3}\) relative error, the PNA position would have matching error, and the resulting \(Z\) would lose the exactness property that Green’s theorem naturally provides on a polygon.
Regularity. The fiber method is designed to tolerate constitutive-law branch points. Applying it to a polynomial integrand is not wrong, just wasteful: the branch-point machinery is carrying weight it does not need. Conversely, Green’s theorem is exact for polynomial integrands and strictly inapplicable to non-polynomial ones — there is no way to bend the formulation to accept \(\sigma(\varepsilon(y, z))\).
Efficiency. A fiber mesh of \(10^{4}\) cells inside a rectangular polygon costs \(10^{4}\) floating-point operations per moment integral; Green’s theorem on the same polygon costs 4. For quantities that never depend on the deformation state, there is no argument for the more expensive option.
GenSec’s allocation of methods¶
Quantity |
Method |
Reason |
|---|---|---|
\(A, S_x, S_y, x_G, y_G\) |
Polygonal |
Polynomial integrand; exact on a polygon. |
\(I_x, I_y, I_{xy}, I_\xi, I_\eta, \alpha\) |
Polygonal |
Polynomial integrand; exact on a polygon. |
\(c_y^{\pm}, c_x^{\pm}, c_\xi^{\pm}, c_\eta^{\pm}\) |
Vertex enumeration |
Deterministic max over the vertex + rebar set. |
\(W_x^{\pm}, W_y^{\pm}, W_\xi^{\pm}, W_\eta^{\pm}\) |
Polygonal |
Derived from \(I / c\); inherits the exactness of the polygonal moments. |
\(Z_x, Z_y, Z_\xi, Z_\eta\) |
Polygonal + bisection |
Polygon split by a half-plane (exact); PNA located by bisection on a monotone objective; first-moment integrals computed in closed form on each half. |
\((N, M_y, M_z)\) for a given \(\mathbf{e}\) |
Fiber (midpoint, Shapely-clipped mesh) |
Non-polynomial integrand with constitutive-law branch points. |
Moment–curvature \(M(\chi)\), failure surface, utilisation ratio |
Fiber + non-linear solver |
Built on top of the fiber response. |
Torsional constant \(I_t\) |
(Future) St.-Venant FEM |
Requires solving \(\nabla^2 \psi = -2\) with \(\psi = 0\) on \(\partial\Omega\); neither polynomial nor expressible through fiber sums. |
Numerical illustration¶
Disc — polygonal vs fiber (with and without clipping)¶
Disc of radius \(R = 500\,\mathrm{mm}\) (exact \(A = \pi R^2\) and \(I = \pi R^4 / 4\)). For the polygonal method, \(N\) is the number of boundary vertices of the regular \(N\)-gon approximating the disc. For the fiber methods, \(N\) is the number of cells per side of the bounding square grid. The naïve column uses a cell-centre mask (in/out); the clipped column uses Shapely intersection to assign each fiber the area and centroid of its cell ∩ disc.
\(N\) |
Polygonal \(N\)-gon |
Fiber, naïve mask |
Fiber, Shapely clipping |
|---|---|---|---|
8 |
\(1.88\!\times\!10^{-1}\) |
\(5.94\!\times\!10^{-2}\) |
\(1.96\!\times\!10^{-2}\) |
16 |
\(5.02\!\times\!10^{-2}\) |
\(6.81\!\times\!10^{-2}\) |
\(5.06\!\times\!10^{-3}\) |
32 |
\(1.28\!\times\!10^{-2}\) |
\(1.90\!\times\!10^{-2}\) |
\(1.29\!\times\!10^{-3}\) |
64 |
\(3.21\!\times\!10^{-3}\) |
\(6.84\!\times\!10^{-3}\) |
\(3.26\!\times\!10^{-4}\) |
128 |
\(8.03\!\times\!10^{-4}\) |
\(3.73\!\times\!10^{-3}\) |
\(8.42\!\times\!10^{-5}\) |
256 |
\(2.01\!\times\!10^{-4}\) |
\(1.50\!\times\!10^{-4}\) |
\(2.34\!\times\!10^{-5}\) |
512 |
\(5.02\!\times\!10^{-5}\) |
\(4.45\!\times\!10^{-5}\) |
|
1024 |
\(1.26\!\times\!10^{-5}\) |
\(1.03\!\times\!10^{-4}\) |
Three observations.
The polygonal column shows a clean, monotone \(\mathcal{O}(N^{-2})\) decay: every doubling of \(N\) reduces the error by a factor close to 4 throughout. The only limit is boundary approximation; it is predictable and stable.
The naïve mask column shows non-monotone convergence with fluctuations caused by the boundary-fitting effect: the sign and magnitude of the local error depend on how each cell straddles the disc boundary, and that pattern changes erratically with \(N\). At \(N = 1024\) the error is actually larger than at \(N = 512\). This is the pathology GenSec avoids.
The Shapely-clipped column recovers the clean \(\mathcal{O}(h^2)\) rate of the clipped-mesh midpoint rule. The ratio between successive rows is again close to 4. For a given \(N\), the clipped fiber is about four times more accurate than the polygonal \(N\)-gon; this is because the \(N \times N\) grid produces \(\mathcal{O}(N^2)\) cells instead of \(\mathcal{O}(N)\) boundary vertices.
RC benchmark — fiber method vs exact PR integration¶
Rectangular section \(B \times H = 300 \times 500\, \mathrm{mm}\), centred on the origin, cover 40 mm. Two rebars φ12 at \(y = +210\,\mathrm{mm}\) (compression side) and three rebars φ20 at \(y = -210\,\mathrm{mm}\) (tension side), all embedded. Concrete EC2 parabola–rectangle with \(f_{cd} = 17\,\mathrm{MPa}\), \(\varepsilon_{c2} = 0.002\), \(\varepsilon_{cu2} = 0.0035\), no tension capacity. Steel B450C (\(f_{yd} = 391.3\,\mathrm{MPa}\), \(E_s = 200\,000\,\mathrm{MPa}\)), symmetric elastoplastic law. Prescribed ultimate strain profile: \(\varepsilon_{\mathrm{top}} = \varepsilon_{cu2} = 0.0035\), neutral axis at \(y = 50\,\mathrm{mm}\) (NA depth from compressed edge \(x = 200\,\mathrm{mm}\)), giving \(\varepsilon_0 = -8.75\!\times\!10^{-4}\), \(\chi = 1.75\!\times\!10^{-5}\,\mathrm{mm}^{-1}\).
The concrete integrand contains all three regions in the compressed zone — tension cut-off from \(y = -250\) to \(y = 50\), parabola from \(y = 50\) to \(y = 164.29\), plateau from \(y = 164.29\) to \(y = 250\) — so both of the concrete branch points at \(\varepsilon = 0\) and \(\varepsilon = \varepsilon_{c2}\) lie inside the integration domain. The tension rebar yields (\(\varepsilon_s = -4.55\!\times\!10^{-3}\)); the compression rebar also yields (\(\varepsilon_s = +2.80\!\times\!10^{-3}\)).
Exact reference obtained by closed-form integration of each region of the concrete law and exact point evaluation of the rebar contributions (with embedded bulk subtraction):
\(N_x\) |
\(N_y\) |
\(n_{\mathrm{fib}}\) |
\(N\) [kN] |
\(M\) [kN·m] |
\(N\) rel. err. |
\(M\) rel. err. |
\(t\) [ms] |
|---|---|---|---|---|---|---|---|
30 |
50 |
1,500 |
541.956 |
232.957 |
\(6.80\!\times\!10^{-4}\) |
\(1.38\!\times\!10^{-5}\) |
0.12 |
60 |
100 |
6,000 |
541.681 |
232.960 |
\(1.73\!\times\!10^{-4}\) |
\(2.47\!\times\!10^{-6}\) |
0.15 |
120 |
200 |
24,000 |
541.611 |
232.960 |
\(4.31\!\times\!10^{-5}\) |
\(6.43\!\times\!10^{-7}\) |
4.00 |
240 |
400 |
96,000 |
541.593 |
232.960 |
\(1.07\!\times\!10^{-5}\) |
\(1.83\!\times\!10^{-7}\) |
9.63 |
480 |
800 |
384,000 |
541.589 |
232.960 |
\(2.68\!\times\!10^{-6}\) |
\(4.38\!\times\!10^{-8}\) |
26.27 |
960 |
1,600 |
1,536,000 |
541.588 |
232.960 |
\(6.71\!\times\!10^{-7}\) |
\(1.10\!\times\!10^{-8}\) |
133.7 |
The ratio between consecutive rows is \(\approx 4\) for both \(N\) and \(M\), confirming the \(\mathcal{O}(h^2)\) rate despite the two concrete branch points inside the integration domain.
A second observation worth noting: up to about 10,000 fibers the per-call time is dominated by Python and NumPy setup overhead rather than by the actual floating-point work. The step from 1,500 to 6,000 fibers costs essentially nothing in wall time but buys four times more accuracy; the step from 6,000 to 24,000 costs 25× more time for the next factor of 4. This explains why the default GenSec mesh sizes never go below 2–3 mm: there is no reason to.
Comparison with the EC2 rectangular stress-block approximation¶
The rectangular stress-block method — used in engineering hand calculations and in most national annexes to EC2 — replaces the full parabola–rectangle law by an equivalent uniform stress \(\eta f_{cd}\) acting over a reduced depth \(\lambda x\) from the compressed edge, with \(\lambda = 0.8\), \(\eta = 1.0\) for concrete classes up to C50/60. Applying it to the same RC section at the same strain state:
Two points.
The stress block is calibrated to match the bending moment accurately at ultimate — the reported \(0.4\%\) error on \(M\) is fully within engineering tolerance. The \(1.8\%\) discrepancy on \(N\), on the other hand, is the model-choice price paid for the simplification: the block has the correct resultant moment but not the correct resultant axial force. This is harmless in pure flexural verification (where \(N = 0\) is imposed by equilibrium and the neutral axis is recalibrated to match) but becomes significant in flexural–axial verification of columns, where a \(1.8\%\) error in \(N\) shifts the interaction curve visibly.
The fiber solver, in contrast, integrates the real PR law with relative error of \(10^{-4}\) at 6,000 fibers and \(10^{-5}\) at 24,000 fibers — several orders of magnitude tighter than the stress-block intrinsic error, and orders of magnitude tighter than any plausible modelling uncertainty on the constitutive parameters.
Engineering cost budget¶
A concrete example. A full 3-D N–Mx–My
interaction surface generated by GenSec (capacity.py)
typically samples a few hundred axial-force levels times a few
hundred angular directions, giving
\(n_{\mathrm{pts}} \sim 10^{4}\) points on the surface.
Each point costs one capacity search that internally runs a
Newton solver with \(n_{\mathrm{iter}} \lesssim 15\)
iterations, each requiring one fiber integration.
With 6,000 fibers (~5 mm mesh, error \(10^{-4}\)):
With 24,000 fibers (~2.5 mm mesh, error \(6\!\times\!10^{-7}\)):
The first is the regime GenSec targets by default; the second is unnecessary for any realistic engineering verification and is reserved for convergence studies or for sections where fiber-level accuracy matters (e.g. high-strength concrete with sharp softening branches, or sections near the boundary of the stress-block validity range).
Summary¶
Use Green’s theorem on the polygon whenever the integrand is polynomial and the domain is polygonal: it is exact, cheap, and stable. Use fiber summation — midpoint rule on a Shapely-clipped mesh — whenever the integrand depends on a non-polynomial constitutive law: clipping eliminates the boundary error, and the branch points of the constitutive law still allow \(\mathcal{O}(h^2)\) global convergence. Do not mix the two for the same quantity; the accuracy and consistency arguments above make the separation worth enforcing. Do not reach for higher-order quadrature rules: they provide no benefit on \(C^0\) integrands and can in fact be worse.
On a representative RC section with realistic rebars, 6,000 fibers give \(\sim 10^{-4}\) accuracy in 0.15 ms per integration — about three orders of magnitude more accurate than the EC2 rectangular stress-block hand calculation, and several orders of magnitude more accurate than the modelling uncertainty on the constitutive constants themselves.