Arostica Benchmark

[javijv4]

I have been running the idealized left-ventricle benchmark (case 1A) from Aróstica et al. (2025) https://doi.org/10.1016/j.cma.2024.117485. The simulation as it stands in tests/cases/struct/LV_HolzapfelOgden_active does not match the published displacements well (see the black trace in the Figure). I narrowed down the discrepancy to two issues:

1. Displacement order. Switching to quadratic (P2) elements improves the result substantially (red dashed curve). I think this is mainly volumetric locking of the P1 tetrahedra.
2. Fiber interpolation. svMultiPhysics currently reads fibers only as cell data (a P0 field). That represents poorly how the fibers change under higher-order displacements. I asked Codex to add nodal fibers on this branch. Keeping the displacements linear and moving the fibers onto the nodes (P1) also improves the result (blue dashed curve), though by less than the quadratic mesh.

Neither change alone puts the displacements inside the published range. Together they do. Quadratic displacements with nodal fibers of the same order (P2) fall inside that range (purple curve). The gray band is the minimum-to-maximum envelope of the groups in the benchmark dataset.

A few points for discussion:

1. Do we want to add support for nodal fibers? The design I have in mind is a fiber-loader class that, while reading the file, interpolates the directions to the quadrature points, orthonormalizes them while keeping the primary fiber direction, and stores that array in ComMod (where the current cell fibers are stored). The physics that uses fibers then stays almost untouched (just needs to point to this new array). That is what I asked Codex to implement. I have not reviewed that code in depth, so it probably needs a decent amount of refinement.
2. I have not tried ustruct. A mixed formulation would address the locking, and the displacements stay linear, so the fiber-field order may not matter as much. Most groups in the benchmark used a compressible P2 formulation, though, so I lean towards getting that formulation working as expected.

Very small note: The Python code to generate the pressure/active stress curve used scipy solve_ivp with no restriction on the integrator step size. This made the active stress curve a bit off. We should update that file in the tests, but I will wait until we know whether a subsequent PR will fix it. It won't affect the reference solutions, since this discrepancy appeared after about 180 timesteps.