corotational MMS 3D#27
Conversation
Fimache
commented
Jun 19, 2026
Fimache
commented
- SinusCorotational: u_ex = RBM (45 deg) + sinusoidal u_d
- HexahedronCorotationalFEMForceField with rotationMethod=large
- Issue: mesh explodes at init due to large rigid rotation in BCs
- Tried: shift trick on all nodes, boundary-only shift, rest_position sync
- None converge: Newton residual ~ 1e11 after 5 iterations
th-skam
left a comment
th-skam
left a comment
There was a problem hiding this comment.
You seem to be on an old branch. You'll have to rebase to get a clean diff.
There are some conflicts now that are not resolved automatically.
| nodes_3d = get_nodes_3d(L, nx, ny, nz) | ||
|
|
||
| Grid = rootNode.addChild("Grid") | ||
| Grid.addObject("RegularGridTopology", name="grid", | ||
| nx=nx, ny=ny, nz=nz, | ||
| min=[0.0, 0.0, 0.0], max=[L, L, L]) | ||
|
|
||
| Solid = rootNode.addChild("Solid") | ||
| Solid.addObject("StaticSolver", name="staticSolver", printLog=False) | ||
| Solid.addObject("NewtonRaphsonSolver", name="newtonSolver", | ||
| maxNbIterationsNewton=5, | ||
| absoluteResidualStoppingThreshold=1e-10, | ||
| printLog=False) | ||
| Solid.addObject("SparseLDLSolver", name="linearSolver", | ||
| template="CompressedRowSparseMatrixd") | ||
|
|
||
| dofs = Solid.addObject("MechanicalObject", name="dofs", template="Vec3d", | ||
| position=nodes_3d.tolist(), | ||
| showObject=with_visual, showObjectScale=0.005 * L) |
There was a problem hiding this comment.
It's better to just create the RegularGridTopology once and use the positions from here moving forward.
Instead, what you do now, is to use the get_nodes_3d to create nodes, you also create the mesh with RegularGridTopology, then you copy the nodes_3d into the MechObject. Just use the RegularGridTopology.
There was a problem hiding this comment.
Ties into the init-order bug I found. With src="@../Grid/grid", positions resolve during init(), so apply_bcs has to run after init() now (writing to dofs.position directly), not before, or SOFA overwrites it.
Plan: drop get_nodes_3d, use src="@../Grid/grid" on the MechanicalObject, move apply_bcs into solve_solid post-init. For case_scene (GUI), I'll need a controller on onSimulationInitDoneEvent to apply BCs at the right time.
There was a problem hiding this comment.
Yes exactly, best would be a controller. You cannot manually rotate your dofs' positions at the createScene time. I said it in previous comment but you can take a look at the AffineMovementProjectiveConstraint I mentioned in previous comment.
| # ── Shift trick sur TOUS les nœuds ──────────────────────────────── | ||
| # → Newton part d'un état cohérent, proche de la solution | ||
| with dofs.position.writeable() as pos: | ||
| for i, (x, y, z) in enumerate(xyz): | ||
| ux, uy, uz = self.u_ex(x, y, z, L) | ||
| pos[i] = [x + ux, y + uy, z + uz] | ||
|
|
||
|
|
||
| boundary_indices = [ | ||
| i for i, (x, y, z) in enumerate(xyz) | ||
| if x < eps or x > L - eps | ||
| or y < eps or y > L - eps | ||
| or z < eps or z > L - eps | ||
| ] | ||
| Solid.addObject("FixedProjectiveConstraint", | ||
| name="fix_boundary", | ||
| indices=boundary_indices) |
There was a problem hiding this comment.
What is the purpose of this shift trick? Could you explain what you are trying to achieve?
Are you changing the FixedProjectiveConstraint on the fly?
There was a problem hiding this comment.
-
The manufactured solution u_ex contains a 45° rigid-body rotation component. If Newton starts from position = X (reference configuration), the corotational polar decomposition sees identity rotation on every element at iteration 0; the linearised tangent is completely wrong for a 45° rotation, and Newton diverges. The shift trick sets position <== X + u_ex(X) as the initial guess, so the polar decomposition immediately extracts the correct rotation. Newton then only has to converge the small deformational residual, not "discover" a 45° rotation from scratch.
-
No; that was a bug in the previous version. apply_bcs was calling dofs.position.writeable() during scene construction (before Sofa.Simulation.init()), which has no effect since the MechanicalObject is not yet initialised. In the corrected version, the shifted positions X + u_ex(X) are passed directly to the MechanicalObject constructor in solidCoro.py, so SOFA sets both position and rest_position to the shifted values at init time. apply_bcs now only adds the FixedProjectiveConstraint once; it is never modified after init.
There was a problem hiding this comment.
- Thanks, this is clear now.
- Not so much of a bug. The strategy with the FixedProjectiveConstraint is not ideal, though I get why you tried to use it. You would need to manually orchestrate the simulation through python or, better, use a controller to do this properly. Another possibility is to use the
AffineMovementProjectiveConstraintwhich applies the rotations from within SOFA incrementally. I am doing some tests with it. I suggest we review that.
Also, I saw somewhere you are additionally changing the rest_positions. Problem with that is that the force field does not make a correct comparison to determine the affine rotation and subtract it from the force computation due to strain.
| # ── AJOUT 2 : rest_position = position shiftée après init ───────── | ||
| # apply_bcs a placé dofs.position à X + u_ex(X) | ||
| # On synchronise rest_position pour que SOFA parte de là | ||
| with dofs.rest_position.writeable() as rpos: | ||
| rpos[:] = dofs.position.array().copy() | ||
|
|
||
| conn = elem.read_connectivity(topology) | ||
| Sofa.Simulation.animate(root, root.dt.value) | ||
| pos1 = dofs.position.array().copy() | ||
| Sofa.Simulation.unload(root) | ||
|
|
||
| # ── AJOUT 3 : déplacement = pos_finale - X_original ─────────────── | ||
| # (pas pos_finale - pos0_shiftée comme avant) | ||
| ux = pos1[:, 0] - nodes_3d_original[:, 0] | ||
| uy = pos1[:, 1] - nodes_3d_original[:, 1] | ||
| uz = pos1[:, 2] - nodes_3d_original[:, 2] | ||
|
|
||
| return SolidSolution3D( | ||
| nodes=nodes_3d_original, conn=conn, ux=ux, uy=uy, uz=uz | ||
| ) |
There was a problem hiding this comment.
Let's revisit this logic, I don't think you should play with the reference position yourself. Let's discuss this
|
One thing that is wrong here is the amplitude of the sinusoids. 
you can't use an amplitude that makes your mesh invert like that. We've seen this already in 2D. |
2 participants
