prop2 returns Cw = 0 for every open section (wrong column + wrong array in the unit-warping loop)

[miami3pl]

prop2 returns Cw = 0 for every open section

Version: pycufsm 0.2.0 · Python: 3.14 · numpy: 2.4.6 (see issue 2)

Summary

pycufsm.pre.cutwp.prop2 computes a correct area, Ixx, Iyy, J and shear
centre for an open section, and then returns Cw = 0. Every other property
is right, so the zero is easy to miss — and examples/example_1_new.py passes
"Cw": 0 in its hardcoded sect_props, which suggests it has been.

Reproducer

repro_cw_zero.py (attached below) models an ArcelorMittal UPE 240 channel
at its mid-line, so the expected answer is independently published:
Iw = 26,400 cm⁶ = 2.64e10 mm⁴ (fillet-aware; a mid-line model gives ≈2.78e10).

Output on 0.2.0:

A   = 3755.0 mm^2  Ixx = 3.4849e+07  J = 1.3864e+05  x0 = -34.73 mm
Cw  = 0.0

A, Ixx, J and x0 all match the published values (J is ~8% low only
because a mid-line model omits the root fillet). Only Cw is wrong.

Cause

The "Compute unit warping" loop (pre/cutwp.py, ~line 269) is a copy of the
sectorial-area loop above it, intended to run about the shear centre and
fill wo_vals. Two things did not get changed in the copy:

start_node = int(ends[i, 0]) - 1
end_node   = int(ends[i, 0]) - 1      # (1) column 0 — should be column 1
...
if w_vals[start_node, 0] == 0:        # (2) w_vals — should be wo_vals
    w_vals[start_node, 0] = start_node + 1
    w_vals[start_node, 1] = w_vals[end_node, 1] - po_vals * lengths[i]
elif w_vals[end_node, 0] == 0:
    ...
w_no = w_no + 1 / (2 * area) * (wo_vals[start_node, 1] + wo_vals[end_node, 1]) * ...

1. With start_node == end_node, po_vals is the cross product of a point with
   itself, i.e. identically 0, so no sectorial coordinate accumulates.
2. w_vals is already fully populated by the preceding loop, so neither branch
   of the if/elif ever fires and wo_vals stays zero.

The w_no accumulation on the last line already reads wo_vals, which is the
strongest hint that wo_vals is what the loop was meant to fill.

Downstream: wn_vals = w_no - wo_vals[:, 1] is all zeros, so C_warping sums
zeros.

Fix

             start_node = int(ends[i, 0]) - 1
-            end_node = int(ends[i, 0]) - 1
+            end_node = int(ends[i, 1]) - 1
             po_vals = (...)
-            if w_vals[start_node, 0] == 0:
-                w_vals[start_node, 0] = start_node + 1
-                w_vals[start_node, 1] = w_vals[end_node, 1] - po_vals * lengths[i]
-            elif w_vals[end_node, 0] == 0:
-                w_vals[end_node, 0] = end_node + 1
-                w_vals[end_node, 1] = w_vals[start_node, 1] + po_vals * lengths[i]
+            if wo_vals[start_node, 0] == 0:
+                wo_vals[start_node, 0] = start_node + 1
+                wo_vals[start_node, 1] = wo_vals[end_node, 1] - po_vals * lengths[i]
+            elif wo_vals[end_node, 0] == 0:
+                wo_vals[end_node, 0] = end_node + 1
+                wo_vals[end_node, 1] = wo_vals[start_node, 1] + po_vals * lengths[i]

The traversal guard at the top of the same loop also tests w_vals and must
test the set being built, or the while runs to the end of ends every time:

-                (np.any(w_vals[:, 0] == ends[i, 0]) and np.any(w_vals[:, 0] == ends[i, 1]))
-                or (not (np.any(w_vals[:, 0] == ends[i, 0])) and (not np.any(w_vals[:, 0] == ends[i, 1])))
+                (np.any(wo_vals[:, 0] == ends[i, 0]) and np.any(wo_vals[:, 0] == ends[i, 1]))
+                or (not (np.any(wo_vals[:, 0] == ends[i, 0])) and (not np.any(wo_vals[:, 0] == ends[i, 1])))

Validation of the fix

Three plain hot-rolled channels, against published Iw and against an
independent mid-line sectorial integration (a separate implementation
normalising on ∮ω·t·ds = 0, written without reference to CUTWP):

section	prop2 after fix	independent integration	published Iw
UPE 200	1.157e10 mm⁶	1.156e10	1.100e10
UPE 240	2.776e10	2.775e10	2.640e10
UPE 300	7.546e10	7.543e10	7.270e10

The two implementations agree to 0.032%. Both sit ~4% above the published
value, which is the expected mid-line-vs-fillet-aware difference and is
consistent with the same ~6–8% gap seen in J.

A suggested regression test: assert Cw > 0 and within a few percent of
2.78e10 mm⁶ for the UPE 240 geometry in the reproducer.

Reproducer (standalone)

"""Minimal reproducer: pycufsm returns Cw = 0 for any open section.

Run with the ORIGINAL (unpatched) pycufsm 0.2.0.
A plain hot-rolled channel (ArcelorMittal UPE 240) is used because its warping
constant is published, so the expected answer is independently known:
    Iw = 26,400 cm^6 = 2.64e10 mm^6   (fillet-aware; a mid-line model gives ~2.78e10)
"""
import numpy as np
from pycufsm.pre.cutwp import prop2_new

# UPE 240 mid-line: h=240, b=90, tw=7, tf=12.5  ->  hm = 227.5, bm = 86.5
hm, bm, tw, tf = 227.5, 86.5, 7.0, 12.5
n_web, n_fl = 16, 6
nodes = [[bm * (1 - i / n_fl), 0.0] for i in range(n_fl + 1)]
nodes += [[0.0, hm * i / n_web] for i in range(1, n_web + 1)]
nodes += [[bm * i / n_fl, hm] for i in range(1, n_fl + 1)]

k = 0
elements = [{"nodes": list(range(k, k + n_fl + 1)), "t": tf, "mat": "S"}]
k += n_fl
elements.append({"nodes": list(range(k, k + n_web + 1)), "t": tw, "mat": "S"})
k += n_web
elements.append({"nodes": list(range(k, k + n_fl + 1)), "t": tf, "mat": "S"})

sp = prop2_new(np.array(nodes), elements)
print("A   = %.1f mm^2      (published 3850)" % float(sp["A"]))
print("Ixx = %.4e mm^4   (published 3.60e7)" % float(sp["Ixx"]))
print("J   = %.4e mm^4   (published 1.51e5; no fillet here, so ~8%% low)" % float(sp["J"]))
print("x0  = %.2f mm       (shear centre; expected -34.7)" % float(sp["x0"]))
print()
print("Cw  = %.4e mm^6   <-- EXPECTED ~2.78e10 (mid-line), published 2.64e10"
      % float(sp["Cw"]))
assert float(sp["Cw"]) > 0, "BUG: Cw is zero for an open section"