[mitchell-tesch]

While testing the geometrical non-linear solver, I found that solve(geometrical_non_linear=True, discretize_kwargs=dict(n=20)), as shown in the docs, fails for a horizontal member with:

ValueError: operands could not be broadcast together with shapes (20,20) (60,60)

That came from two separate bugs.

1. Solving a model a second time can give wrong results. This affects plain linear analysis too:

ss = SystemElements(EI=1810, EA=460_000)
ss.add_multiple_elements([[0, 0], [6, 0]], 4)
ss.add_support_hinged(1)
ss.add_support_roll(ss.id_last_node)
ss.point_load(ss.id_last_node, Fx=-10)
ss.solve()

ss.remove_loads()
ss.point_load(3, Fy=-10)
ss.solve()
ss.get_node_displacements(3)["uy"]  # 0.0, but a fresh model gives -0.0249

process_conditions() treats any DOF that's exactly 0 in system_displacement_vector as a support. After a solve, that vector
holds the results, so any free DOF that happened to come out as 0 (here, every vertical displacement and rotation) is locked on the next solve. The non-linear solver always solves more than once, so it runs into this every time.
The fix records the restrained DOFs once in process_supports(), and process_conditions() uses that record.

2. discretize() loses distributed loads. It re-applies each q-load with its two components swapped, and only to the first sub-element. It now splits each element's load over all of its sub-elements, the same way insert_node() already does.

Docs (second commit). The docs recommend discretize_kwargs as if it refines the whole second-order calculation. It only meshes a copy of the model that's used for the buckling factor (the solve() docstring says this), so the moments within a single-element member stay first order. I caught myself out with this. The docs now suggest ss.discretize(n) on
the model itself, and describe discretize_kwargs as affecting the buckling factor only.

Tests

tests/test_resolve_and_discretize.py checks that a re-solved model matches a fresh one. It also checks that the reactions after discretize() match the original model for uniform, varying and inclined loads. All four tests fail on master. black, mypy, pylint and pytest pass locally.

This PR doesn't depend on #490, and they merge cleanly together. With both applied, discretize_kwargs gives the Euler load for the column in examples/ex_18_discretize.py (493.48 with n=15, against 600 without).

A question

Would you like discretize_kwargs to also apply to the second-order results, or to be deprecated in favour of ss.discretize()? I've left the behaviour as it is and only corrected the docs, but I'm happy tofollow up either way.

Note: with this PR alone, solve(geometrical_non_linear=True, discretize_kwargs=...) gets past the shape mismatch but still fails in det_linear_buckling() (np.linalg.inv(kg): singular matrix). #490 fixes that eigen-solve, so the call needs both PRs. They are independent and merge cleanly in either order.

Known limitation (unchanged): discretize() passes the element's spring to every sub-element, so a release ends up at every internal sub-node and the model becomes a mechanism. #482 already reworks this, so I've left it out here.