Hi FreeFEM developers!

I thank you for having integrated PETSc with FreeFEM!

This has resulted in an easy implementation of large scale topology optimization.

I am attempting to implement a contact mechanics solver in parallel (for topology optimization).

Implemented via fixed point and penalization, the contact problem reads

(K + M(U_n))U_{n+1} = F

where U_n is the solution iterate, K is the linear elasticity matrix, F is the external force and M(U_n) is the contact penalization term.

To simplify the construction of M(U_n), I work with co-incident contact surfaces (in 2D and in 3D).

In this case, it suffices to determine indices of the coincident vertices, and the corresponding pairing matrix.

Using the pairing matrix, I construct the matrix M(U_n) by specifying values at certain indices (ones corresponding to the contact surface).

For PETSc, I need to create subdomains with ghost elements. This complicates the creation of the pairing matrix and the construction of the matrix M(U_n).

Is there a possibility to simply create the global matrix M(U_n) in sequential and then feed it to a macro like createMat?

If not, is there any other solution?