Extract continuous-time system matrices

[jjalora]

Hello! I was wondering how one could go about extracting the system matrices (i.e., M, C, K, and H) at a given equilibrium configuration. Currently, I'm extracting the global matrix and then algebraically computing the jacobians. This, unfortunately, introduces a dependency on the time discretization, which for my purpose is undesirable.

I'm aware there is some ongoing work on making these matrices available (ref), but rather than wait, I'd like to try to see if there is some way to extract this on my own. Currently, I'm taking a look at computeStiffnessMatrix in TetrahedronFEMForceField here and it looks like I should perhaps extract SR (which looks to be a rotated JKJt, but I'm not 100% sure). I'm sure I'm not the only one who has tried to do this - I'm wondering if anyone else has figured this out and if there is an easy hack to make this matrix accessible through the python interface? (I'm interfacing with SOFA exclusively through SofaPython)

I'm new to SOFA so even the simplest details/suggestions would be very instructive. Thank you so much for the help!

[alxbilger]

Hi @jjalora

you are right that there is an ongoing work on this topic. It should be ready for the next release in June (I hope).

The matrix K is computed in the addKToMatrix functions of all force fields. Unfortunately, the contributions are all summed up into a global matrix. You can extract the global matrix, but cannot retrieve K from it.

If you can't wait for the next release, you'll need to hack the C++ code to get the K matrix. Here is what I suggest:

• In MatrixLinearSolver, create another member dedicated to the K matrix (Matrix kMatrix;).
• Go to MatrixLinearSolver::setSystemMBKMatrix.
• After mops.addMBK_ToMatrix, call it again but with mFactor set to 0, bFactor set to 0 and kFactor set to 1. The matrix accessor must reference your new matrix kMatrix. kMatrix will be populated from the contributions of all the force fields. It is also important to note that this matrix will be "global", in the sense that it will have the dimension of the problem solved by your solver and not the dimension of the associated mechanical object. Both numbers are different if you solve multiple objects with the same solver.
• Make kMatrix accessible to Python.

Note that the hack I provided has not been tested. You'll need to setup the kMatrix properly, along with the matrix accessor.

Good luck,

Alex

[jjalora]

@alxbilger - Thank you for the detailed response! I'll try this out and will report back here with any findings. Cheers!

-John

[artXuRo]

Hi,
I was trying to extract M, D, K, and H matrix from sofa, it seems to work when I set 'savingMatrixToFile=True', and I can get a TXT file that contains the 3*n(n is the number of nodes) matrix. As explained in the sofa, this file contains M, D, and K matrix. I was wondering how to get the H matrix which distributes inputs to all the elements which contain the cable nodes.

Best
Art

[hugtalbot]

Great to read that you succeeded @artXuRo !!!!

You can find the H matrix in a specific format which I recently described in another post. Is this what you are looking for?

[artXuRo]

Thanks for the response!
As of now, I can extract the constraint matrix through the following code:
constraint = self.cable1.dof1.constraint.value and the data type is char.
whether I can gain the matrix in a vec3d type?

And I have another question is, when I gain the system matrix from SparseLDLSolver, I get a large TXT file that contains 1132 matrices.
As the node is 377. And I can't distinguish which are Mass matrices, Dumping matrices, and Stiffness matrices.

Best
Art

[alxbilger]

Yes, the SparseLDLSolver Data savingMatrixToFile saves the global matrix which is an accumulation of the mass and stiffness matrices. And it is impossible to extract them individually from it.
With the pull requested I mentioned, you can assemble the stiffness matrix or the mass matrix.
Note that the Data savingMatrixToFile in SparseLDLSolver will be removed in v22.06, in favor of a more generic matrix exporter. It is already available in v21.12. The component name is GlobalSystemMatrixExporter, and can work with any linear solver, not just SparseLDLSolver.

[alxbilger]

@jjalora and @artXuRo

Accessing the matrices is still a work in progress, but you can find an open pull request of SofaPython3 where you can do what you want: sofa-framework/SofaPython3#262

Any feedback is welcome.

[artXuRo]

It would be great if I can extract the matrix I want directly.
I can't waiting for the new vesion~

[jjalora]

Indeed, this would be extremely useful. For now we've resorted to inferring the modes of the system using PCA but since, in principal, we should have access to the system matrices, we'd like to extract the modes exactly.

@alxbilger, thanks for doing the work to make this happen and thanks for pointing out the pull request. I'll play around with it and let you know!

[alxbilger]

It would be great if I can extract the matrix I want directly. I can't waiting for the new vesion~

@artXuRo if you compile SOFA and SofaPython3 you can do it right now ;)