Fix product of inertia (Ixy) weighting in cutwp.prop2 (part of #128)

[gaoflow]

Problem

In cutwp.prop2, the product of inertia Ixy is computed with a misplaced closing parenthesis:

I_xy = np.sum((x_diffs * y_diffs / 12 + (x_means - x_centroid) * (y_means - y_centroid) * lengths * thicknesses))

Compare with how Ixx and Iyy are formed just above:

I_xx = np.sum((y_diffs**2 / 12 + (y_means - y_centroid) ** 2) * lengths * thicknesses)
I_yy = np.sum((x_diffs**2 / 12 + (x_means - x_centroid) ** 2) * lengths * thicknesses)

For Ixx/Iyy the * lengths * thicknesses factor multiplies the whole bracket (the element's own second moment d²/12 plus the parallel-axis term). For Ixy the parenthesis closes too late, so * lengths * thicknesses only multiplies the parallel-axis term and the element's own contribution x_diffs * y_diffs / 12 is left unweighted by the element area.

The result is an incorrect Ixy (and therefore an incorrect principal angle phi and principal moments I11/I22) for any section containing inclined elements. Sections whose elements are all axis-aligned are unaffected, because then x_diffs * y_diffs is zero for every element — which is likely why this went unnoticed.

This is one of the two issues reported in #128.

Fix

Move the closing parenthesis so the full bracket is weighted by lengths * thicknesses, consistent with Ixx/Iyy:

I_xy = np.sum((x_diffs * y_diffs / 12 + (x_means - x_centroid) * (y_means - y_centroid)) * lengths * thicknesses)

Verification

For a small section with an inclined element (so the x_diffs * y_diffs / 12 term is non-zero):

coord = np.array([[0.0, 0.0], [4.0, 0.0], [6.0, 3.0]])
ends  = np.array([[0, 1, 0.1], [1, 2, 0.1]], dtype=float)
cutwp.prop2(coord, ends)["Ixy"]
# before: 1.3533230611135743
# after : 1.0336006248867737

The fixed value matches both an independent thin-walled product-of-inertia formula and direct numerical integration of (x - xc)(y - yc) along the elements (agreement to ~1e-9). The rotation invariant Ixx + Iyy == I11 + I22 continues to hold, and a section symmetric about an axis still yields Ixy == 0.

Added tests/test_cutwp.py with these checks; the inclined-element test fails on master and passes with this change. The rest of the suite is unchanged (the one pre-existing it_results_in_correct_signature_curve failure fails identically before and after this change).

Note

#128 mentions two problems in prop2. This PR addresses the product-of-inertia one, which is self-contained and straightforward to verify. The other appears to be in the shear-center/unit-warping loop (the warping constant path) and needs separate, carefully validated treatment, so I've kept it out of this change rather than ship a partial fix there.

[brookshsmith]

Thank you so much, @gaoflow ! And thank you for adding the test file for cutwp as well. Simple fix, but you've added a lot value in your PR - merged.