Add quadratic manufactured solution for 1D bar validation#10
Add quadratic manufactured solution for 1D bar validation#10Fimache wants to merge 3 commits intomain-freefem-sofafrom
Conversation
…gence script and parameters
th-skam
left a comment
th-skam
left a comment
There was a problem hiding this comment.
A few things to be careful of. We'll merge this PR after we make sure the results are correct.
Quick Points:
- This is an equilibrium problem, so we should use the static solver, not a time integrator.
- For SOFA, we must use components from the Elasticity plugin. The
SpringForceFieldin the scene should not be there. - The mesh is finer in FreeFEM than in SOFA.
- The errors are not consistently computed between FreeFEM and SOFA. To see this, look at the following images, generated by these scripts:
What is obvious from the 1st image, is that although the solutions are very close, the error is orders of magnitude off for SOFA. Which indicates that the way we compute errors is wrong.
The second image also reveals something is wrong. In a previous PR, we simulated the distributed load over a bar. We saw super-convergence of the nodal error with quadratic polynomial solutions. For the continuous L2 norm, we saw the expected convergence. However, what I see here is:
- For FreeFEM: Super convergence of the nodal error is observed. Also, we see the expected convergence of the L2 norm (integrated over the 1-dimensional domain). That's also correct.
- For SOFA:
- We are computing the RMS in the code which should also exhibit super-convergence. It doesn't because the meshes are not identical and require interpolation within elements. This converges by refining but at a wrong rate, caused by using the wrong
hto compute it.
- We are computing the RMS in the code which should also exhibit super-convergence. It doesn't because the meshes are not identical and require interpolation within elements. This converges by refining but at a wrong rate, caused by using the wrong
Let me know if anything is not clear.
|
|
||
| rootNode.gravity = [0, 0, 0] | ||
|
|
||
| rootNode.dt = 1e6 |
There was a problem hiding this comment.
| rootNode.dt = 1e6 | |
| rootNode.dt = 1 |
Time shouldn't matter for solving this equilibrium problem.
| rootNode.addObject("EulerImplicitSolver", | ||
| name="odeSolver", | ||
| rayleighStiffness=0.0, | ||
| rayleighMass=0.0) |
There was a problem hiding this comment.
We should not use time marching to solve a static equilibrium problem.
Look at https://github.com/InfinyTech3D/Elasticity/blob/main-freefem-sofa/examples/Freefem/solution_man/cas_1D/test_sofa/scenario2.py to see how to set up the static solver. You've already done it.
| rootNode.addObject("StiffSpringForceField", name="springs", | ||
| spring=springs_str) |
There was a problem hiding this comment.
All of this is testing the wrong thing. You should use the LinearStrainFEMForceField with Vec1d from the Elasticity plugin. The component you used only mimics what linear FEM does using a spring model.
| L = length | ||
| f_const = 2.0 * young_modulus / (L * L) |
There was a problem hiding this comment.
| L = length | |
| f_const = 2.0 * young_modulus / (L * L) | |
| f_const = 2.0 * young_modulus / (length * length) |
No reason to redefine this, just use length. Substitute this further down.
| raise ValueError( | ||
| f"Dimensions incompatibles : u_num={len(u_num)}, u_ref={len(u_ref)}." | ||
| ) | ||
| return np.linalg.norm(u_num - u_ref) / np.sqrt(len(u_ref)) |
There was a problem hiding this comment.
This is RMS not just the L2 norm, right? Or something mixed. It's not resulting from integration over the domain as FreeFEM does. So the comparison is meaningless.
| def convergence_rate(h_arr, err_arr, floor=1e-13): | ||
| h_arr = np.asarray(h_arr) | ||
| err_arr = np.asarray(err_arr) | ||
| mask = err_arr > floor | ||
| if mask.sum() < 2: | ||
| return np.nan | ||
| return np.polyfit(np.log(h_arr[mask]), np.log(err_arr[mask]), 1)[0] |
There was a problem hiding this comment.
Is this a fitting line on the log-log scatter plot to derive a slope value?
| real E = $youngModulus; | ||
| int N = $nx; | ||
|
|
||
| meshL Th = segment(N, [x*L, 0, 0]); |
There was a problem hiding this comment.
| meshL Th = segment(N, [x*L, 0, 0]); | |
| meshL Th = segment(N-1, [x*L, 0, 0]); |
Be careful here. Look at this:
https://github.com/InfinyTech3D/Elasticity/pull/10/changes#diff-0597ab6b63a9f899e584c61ab4c1aaa3b98264f9ecbce0eae4c7ec166d2bc0c2R98-R119
You are passing the same nx to SOFA and FreeFEM.
- SOFA creates nx nodes
- FreeFEM creates nx elements -> nx + 1 nodes.
So the meshes used are not the same and they don't have the same h. However, the same h is used in the convergence rate analysis.
Let's set it now: N will always refer to the number of nodes (never the elements).
| u_sofa_interp = np.interp(x_ff, x_sofa, u_sofa) | ||
| u_ex_ff = u_exact(x_ff, L) | ||
| err_sofa = l2_error_discrete(u_sofa_interp, u_ex_ff) |
There was a problem hiding this comment.
Why interpolate on a simple 1D example ? If you look at the work we did in previous PRs, we don't need to interpolate.
2 participants
This PR adds a quadratic manufactured solution for 1D bar elasticity to validate SOFA.