[LinearSystem] MatrixProjectionMethod: Fix out-of-bounds matrix writes

[fredroy]

MatrixProjectionMethod::addMappedMatrixToGlobalMatrixEigen iterated over the cartesian product of the union of both mechanical states' top-most ancestors, instead of inputs1 x inputs2.
For the resulting spurious pairs no Jacobian exists, and computeProjection silently falls back to an unprojected block whose dimensions belong to the mapped state, written at the global offsets of an unrelated state.

[with-all-tests]

[ci-depends-on https://github.com/sofa-framework/Regression/pull/126]

---

By submitting this pull request, I acknowledge that
I have read, understand, and agree SOFA Developer Certificate of Origin (DCO).

---

Reviewers will merge this pull-request only if

• it builds with SUCCESS for all platforms on the CI.
• it does not generate new warnings.
• it does not generate new unit test failures.
• it does not generate new scene test failures.
• it does not break API compatibility.
• it is more than 1 week old (or has fast-merge label).

[fredroy]

for windows....
[ci-build][with-all-tests]

[sofabot]

[ci-depends-on] detected during build #4.

To unlock the merge button, you must

• Merge or close Update regression due to PR6256 (MatrixProjectionMethod fix) Regression#126
  For this build, the following CMake flags will be set
  -DREGRESSION_GIT_REPOSITORY='https://github.com/sofa-framework/Regression'
  -DREGRESSION_GIT_TAG='06148ea71569a494d56a7e1b9e6975a54500065e'

[alxbilger]

Strongly related to #6222

[alxbilger]

I am afraid Claude is wrong!

Consider a spring that has 2 input mapped states resulting into 2 mapping jacobian matrices: $J_1$ and $J_2$. The cartesian product suggested by Claude would lead to the pairs $(J_1, J_2)$, while the previous code would lead to $(J_1,J_1)$, $(J_1,J_2)$, $(J_2,J_1)$, $(J_2,J_2)$. The latter is the correct solution. It would lead to $J_1^T K J_1$, $J_1^T K J_2$, $J_2^T K J_1$, $J_2^T K J_2$, i.e. the totality of the coupling terms.

I suggest to close this PR.

[fredroy]

His/Its answer 👀:

---

Thanks for pushing back on this - the requirement you state is right: all four coupling terms J1^T K J1, J1^T K J2, J2^T K J1, J2^T K J2 must end up in the global matrix. Where I think we differ is on where those four terms come from.

They do not come from expanding one K over the union of both ancestor sets. They come from four separate mapped matrices, created upstream in MatrixLinearSystem.inl:253:

inline auto generatePairs(const vector<BaseMechanicalState*>& mstates)
{
for (auto* a : mstates)
for (auto* b : mstates)
pairs.emplace_back({a, b}); // all ordered pairs
}

A spring on states {M1, M2} therefore gets four mapped matrices K11, K12, K21, K22, each with its own MatrixProjectionMethod, and each projected with the jacobians that match its own row and column spaces.

I built exactly your scenario to check - a spring between two mapped states with different parents, global matrix P1(3) + P2(3) = 6:

... same, with M2 ...

With the PR, logging every block handed to addToGlobalMatrix:

K=(M1,M1) a=P1 b=P1 J0=1 J1=1 pos=(0,0) block=3x3 nnz=9
K=(M1,M2) a=P1 b=P2 J0=1 J1=1 pos=(0,3) block=3x3 nnz=9
K=(M2,M1) a=P2 b=P1 J0=1 J1=1 pos=(3,0) block=3x3 nnz=9
K=(M2,M2) a=P2 b=P2 J0=1 J1=1 pos=(3,3) block=3x3 nnz=9

That is precisely J1^T K J1, J1^T K J2, J2^T K J1, J2^T K J2 - all four terms, at the four expected offsets, each with both jacobians present.

On the same scene before the PR, ten blocks are produced. The interesting ones:

K=(M1,M2) a=P1 b=P1 J0=1 J1=0 pos=(0,0) <- spurious
K=(M1,M2) a=P1 b=P2 J0=1 J1=1 pos=(0,3) <- the correct term
K=(M1,M2) a=P2 b=P1 J0=0 J1=0 pos=(3,0) <- spurious
K=(M1,M2) a=P2 b=P2 J0=0 J1=1 pos=(3,3) <- spurious

J0=0 / J1=0 means getJacobianFrom() returned null for that pair - no such jacobian exists - and computeProjection() then silently falls back to treating the missing one as identity. Those blocks are not J^T K J at all. The net effect is that K12 and K21 get smeared onto the (0,0) and (3,3) diagonal blocks, which already receive the correct J1^T K11 J1 and J2^T K22 J2, so the spring is counted three times on the diagonal.

On the case where a single state is mapped from two parents (inputs1 == inputs2 == {P1,P2}), union and cartesian product are identical: 2x2 = 4 pairs either way. Measured on MatrixLinearSystem.scn, grouping blocks by ancestor-set size:

n1,n2	blocks	before	after
2,2	16	4 per call	4 per call - unchanged
1,2 and 2,1	4	9 per call	2 per call

So nothing with a genuine two-jacobian mapping changes. The PR only affects calls where the row state and the column state have different ancestor sets.

The structural reason the union cannot be right: J[0] has to be a jacobian of the row state (|mstatePair[0]| x |a|) and J[1] of the column state. Drawing a from inputs2 asks for a jacobian that does not exist, which is why getJacobianFrom returns null on exactly the union-only pairs. In the crash this PR fixes, those blocks also had dimensions that did not match the states they were written against (a 120x960 block placed where the state is 120), which is how they ran off the end of a 1200x1200 matrix.

One case would change my mind: if you know a mapping configuration where a single K legitimately needs a jacobian from the other state's ancestor set. If so, getJacobianFrom returning null for it is itself the bug, and the identity fallback in computeProjection is what should be fixed rather than the pair enumeration. Happy to dig into that if you have one in mind - and I can share the test scene so you can run it both ways.