Skip to content

[Spring] SpringForceField now handles configuration singularity #5822

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
remibessard wants to merge 1 commit into
base: master
Choose a base branch
from
Open

[Spring] SpringForceField now handles configuration singularity #5822

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 8 additions & 1 deletion ...onent/SolidMechanics/Spring/src/sofa/component/solidmechanics/spring/SpringForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -542,11 +542,18 @@ auto SpringForceField<DataTypes>::computeSpringForce(
// dF = k_s.U.U^T.dX + f/l.(I-U.U^T).dX = ((k_s-f/l).U.U^T + f/l.I).dX
auto& m = springForce->dForce_dX;
Real tgt = forceIntensity * inverseLength;
Real tol = 1e-8;
Copy link
Contributor

@hugtalbot hugtalbot Dec 18, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should we not use std::limit ?

I agree that a user-defined data would not be appropriate, but I do not like hard-coded value. Not sure there is an ideal solution..

Real regParam = 1e-4;
bool isSingular = std::abs(elongation) < tol && std::abs(elongationVelocity) < tol;
if (isSingular)
msg_warning(this) << "!!! We detected a degenerated Spring configuration, we've to add a "
"regularization term to make the matrix invertible !!!";
for(sofa::Index j=0; j<N; ++j )
{
for(sofa::Index k=0; k<N; ++k )
{
m(j,k) = ((Real)spring.ks-tgt) * u[j] * u[k];
m[j][k] = ((Real)spring.ks-tgt) * u[j] * u[k];
m[j][k] += (isSingular) ? regParam * spring.ks * (((j==k) ? 1.0: 0.0) - u[j] * u[k]) : 0.0;
Copy link
Contributor

@alxbilger alxbilger Dec 11, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
m[j][k] += (isSingular) ? regParam * spring.ks * (((j==k) ? 1.0: 0.0) - u[j] * u[k]) : 0.0;
if (isSingular && (j==k))
m[j][k] += regParam * spring.ks * (1 - u[j] * u[k]);

This is more readable.

Comment on lines +555 to +556
Copy link
Contributor

@bakpaul bakpaul Dec 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
m[j][k] = ((Real)spring.ks-tgt) * u[j] * u[k];
m[j][k] += (isSingular) ? regParam * spring.ks * (((j==k) ? 1.0: 0.0) - u[j] * u[k]) : 0.0;
m(j,k) = ((Real)spring.ks-tgt) * u[j] * u[k];
m(j,k) += (isSingular) ? regParam * spring.ks * (((j==k) ? 1.0: 0.0) - u[j] * u[k]) : 0.0;

}
m(j,j) += tgt;
}
Expand Down
Loading