search

sofa-framework / sofa

Real-time multi-physics simulation with an emphasis on medical simulation.
https://www.sofa-framework.org
GNU Lesser General Public License v2.1
931 stars 312 forks source link

Why the mass and stiffness matrices extracted from SOFA do not match the right matrices extracted from COMSOL ? #4954

Open SustechShenCong opened 2 months ago

SustechShenCong commented 2 months ago

In my work I need to extract the dynamical model from SOFA, i.e. extracting mass and stiffness matrices from SOFA. I follow the SOFA instructions which use "mass_matrix = root.mass.assembleMMatrix()" & "stiffness_matrix = root.force_field.assembleKMatrix()" to obtain the mass and stiffness matrices. My test example is a cube with 8 Square unit, which results 27 nodes, i.e. 81dofs. see mesh The COMSOL result of mass matrix is below: Comsol_M The SOFA result of mass matrix is below: SOFA_M

SustechShenCong commented 2 months ago

In this issue I only consider the shape of the matrix, not the detailed values!!!

SustechShenCong commented 2 months ago

Solution refers to https://github.com/sofa-framework/sofa/discussions/3134

alxbilger commented 2 months ago

I assume that you use MeshMatrixMass.

In MeshMatrixMass, the mass integration is distributed on vertices and edges. In a hexahedron, some vertices are not connected by edges. Therefore, this type of interaction is not taken into account in the mass matrix. And that is why, you see a lot more zeros in the matrix in SOFA.

I guess you will have more accurate results with tetrahedra, where this situation does not happen.

SustechShenCong commented 2 months ago

Thank you very much for your prompt reply! Your suggestions will be very helpful to me!

SustechShenCong commented 2 months ago

Another question, If there are some ways to extract the Constraints Matrix and "Internal Force Vector and its Jacobian". Some simple examples would be greatly appreciated!

alxbilger commented 2 months ago

There is a work in progress to extract the constraint matrix. The internal force vector is available in the MechanicalObject, in the force Data field (but it is merged with all the other forces).

SustechShenCong commented 2 months ago

Thank you very much for your reply!