[SofaKernel] ADD: add polynomial springs force fields

[sergeiNikolaev]

This components allow to simulate springs with Polynomial stress strain behavior. The general force computation looks like on image below:

Despite the fact I am using only polynomial springs, in general it is possible to set any first order smooth function as stress-strain relation. The polynomial parameters are set as two arrays. d_polynomialDegree describe set of polynomial degrees for every spring. d_polynomialStiffness describe set of polynomial coefficients sequentially combined in one vector. The coefficients are put from smaller degree to bigger one, and the free coefficient is always zero (since for no strain we have no stress). For examples the coefficients for polynomials [3,2,4] will be put as [a1,a2,a3,b1,b2,c1,c2,c3,c4]. The dedication of jacobian matrix for PolynomialSpringForceField is given below I guess the dedication has to be added as documentation, but I do not know how to add it.

For RestShapeSpring I added an exponetial addition to the denominator to avoid Nan problems for cases when spring has a zero length. As a result, the stress simulation is shifted compared with polynomial values, but it keeps its nonlinearity. Here is the dedication of a derivative

---

This PR:

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

Reviewers will merge only if all these checks are true.

[hugtalbot]

Thanks @sergeiNikolaev for this contribution and making this nice effort of bringing back this work.

• Is any test available with these new classes?
• is there any example scene?

[hugtalbot]

Sorry wrong button.

[sergeiNikolaev]

Thanks @sergeiNikolaev for this contribution and making this nice effort of bringing back this work.

• Is any test available with these new classes?
• is there any example scene?

Thank you for your answer. Yes, there is an example scene, which I also use as a test. But to add it to Sofa I have to remove plugin dependencies from it.

[fredroy]

Amazing work @sergeiNikolaev ! However, as it seems you want into SofaCommon, we have to be a bit strict on the guidelines 👍 So could you:

• remove all useless commented code?
• set prefix for Data<> (d) and member (m) ?
• remove serr with msg_error
• maybe uncomment all the msg_info (as it is not print if you dont enable printLog)
• would be cool to use pragma once instead of the preprocessor guards (ifdef) Thanks !

[adagolodjo]

Hi @sergeiNikolaev, thanks for this contribution. I just stared reading your equations and I will probably have somes question in order to understand more about it.
1) First the reference (paper where you got this force equation from) of the equation 1 2) In the first image why you have absolut value on the current lenght (∆l) if this is a length it probably positif rigth? 3) Why current length ∆l = (∆x^2 + ∆y^2 +∆z^2 )^(1/2) it looks more like the formula of the displacement right ?

[sergeiNikolaev]

Hello. Thank you for your comments, guys.

@Younesssss , answering to your questions:

1. First the reference (paper where you got this force equation from) of the equation 1

It is a general spring equation. You could find it for example here: Nealen A. et al "Physically Based Deformable Models in Computer Graphics" https://matthias-research.github.io/pages/publications/egstar2005.pdf or in StiffSpringForceField SOFA component. A spring force contains two parts: value and direction. Compared with linear spring the value now is computed like a result of polynomial. The computation of direction remains the same. The S (cross-section area) technically has to describe the characteristics of material, but since it depends on sample shape we are trying to simulate, finally we could say that we compute force intensity instead of object force and, thus, have to recompute it based on sample shape and amount of springs we are using for simulation.

2. In the first image why you have absolute value on the current length (∆l) if this is a length it probably positif rigth?

Yes, the length is always positive. I added absolute value since square root, in general, might also mean negative value.

3. Why current length ∆l = (∆x^2 + ∆y^2 +∆z^2 )^(1/2) it looks more like the formula of the displacement right ?

Yes, it is a displacement, but this displacement describes a distance between spring endpoints (actually, p1 - p2), not the displacement between the object and its rest shape. Therefore, I marked these displacements with notation ∆. Spring length is a norm of this distance.

[sergeiNikolaev]

Amazing work @sergeiNikolaev ! However, as it seems you want into SofaCommon, we have to be a bit strict on the guidelines So could you:

• remove all useless commented code?
• set prefix for Data<> (d) and member (m) ?
• remove serr with msg_error
• maybe uncomment all the msg_info (as it is not print if you dont enable printLog)
• would be cool to use pragma once instead of the preprocessor guards (ifdef) Thanks !

@fredroy I have submitted the correction to the code.

[sergeiNikolaev]

Thanks @sergeiNikolaev for this contribution and making this nice effort of bringing back this work.

• Is any test available with these new classes?
• is there any example scene?

@hugtalbot Here is the test to verify polynomial springs (Python2 scripts).

polynomial_spring_test.zip

The test "spring_verification.py" depicts two points (p1, p2) connected with a spring. One point (p1) is fixed, another (p2) is moving in space. The external impact (f_external) to the system could be applied using LinearForceField (you need to make target_force as a Data<> object in order to collect its information, but I have this modification) or ControlPoint (third point (p_control), attached with a linear rest shape spring). In case of a second option, the external force is computed as (p_control - p2)*k_rest_shape_spring. During the simulation process, the position "p2" and forces "f2" for point p2, and also the external force or the position of the control point are saved to data files. The script "spring_verification_draw_result.py" loads all data. The strain for spring is computed as (norm(p2-p1) - l0) / l0 for PolynomialSpring or as norm(p2-p1) / l0 for RestShapePolynomialSpring. I use static solver with many iterations to make spring force equal to external force, but I guess it is also possible to compute the spring force as difference f2 - f_external. Based on this data, a force-strain relation is plotted and compared with polynomial values. And this is the end of the test.

NOTE: Please, keep in mind, that configuration parameters have to be the same for both scripts, in order to get a valid result.

[sergeiNikolaev]

Ah, yes. The comparison of a force-strain relation is done on the last figure (figure 5). All other figures just show the tracking of different parameters through the process.

[damienmarchal]

@guparan, @Younesssss, @hugtalbot any volunteer to integrate the test in the normal test framework so we can merge this asap ?

[hugtalbot]

Hi @sergeiNikolaev thank you for sharing your zip test. I think we will be able to work it out from here. TODO for me:

• Need a test in c++ for checking new functionality, to be implemented from py example.
• Add an example in SOFA.

[sergeiNikolaev]

Hello.

Hi @sergeiNikolaev thank you for sharing your zip test. I think we will be able to work it out from here. TODO for me:

• Need a test in c++ for checking new functionality, to be implemented from py example.
• Add an example in SOFA.

My question is what test it is possible to write for Polynomial springs? Personally me, I verified them by generating some data using the scene from .zip and then compare it visually with the polynomial result. But, I am not sure it is possible to perform the same stuff as a functional test.

[hugtalbot]

Hi @sergeiNikolaev Sorry for the latency on this topic, tests can be made so that we create a simple scene in c++ imposing rest shape and an imposed displacement. Thus, the resulting force should be analytically defined. By the way do you have any latex file for the polynomial equations?

[sergeiNikolaev]

Hi @sergeiNikolaev Sorry for the latency on this topic, tests can be made so that we create a simple scene in c++ imposing rest shape and an imposed displacement. Thus, the resulting force should be analytically defined. By the way do you have any latex file for the polynomial equations?

Thank you for your answer @hugtalbot. Regarding your question, I do not, sorry. I wrote these equations since for me it is easier to process the stuff like this with a "writing memory". And then I decided that it will be enough just to scan them and send to you.

[hugtalbot]

@sergeiNikolaev would you also have a simple example scene ?

[hugtalbot]

[ci-build][with-all-tests]

[hugtalbot]

Test and scene examples are added. Doc is ready here: https://github.com/sofa-framework/doc/pull/38/files

Looks now fine to me. Anyone else I merge by the end of the week.

[hugtalbot]

changes taken into account, only the test to go!