# PETSc and Mesh interpolation
**URL:** https://community.freefem.org/t/petsc-and-mesh-interpolation/4109
**Category:** General Discussion
**Created:** [October 21, 2025, 10:16am UTC](https://community.freefem.org/t/petsc-and-mesh-interpolation/4109 "2025-10-21T10:16:25Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@HMasson](https://community.freefem.org/u/HMasson)
#### Post date: [October 21, 2025, 10:16am UTC](https://community.freefem.org/t/petsc-and-mesh-interpolation/4109/1 "2025-10-21T10:16:25Z")
Hello FreeFem community,
I am trying to solve a 3D elastic problem using parallel computing with PETSc and multiple meshes.
In my problem, I use 2 different meshes, one is an optimization domain (Inner mesh), and the second is a bigger mesh (Full Mesh) encompassing the Inner Mesh. Parallel FEA should be performed on the Full Mesh only. The reason behind the use of 2 meshes is that the shape of the mesh is complex, and only a small portion inside is subject to Topology Optimization. The 2 meshes were constructed such that the location of nodes in the Full and Inner mesh correspond.
Following the examples I managed to compute displacements on the Full Mesh using PETSc, and I would like to project those displacements onto the Inner Mesh (to compute sensibilities relative to element densities).
When no parallel computing is used, I was able to make it work simply by using Fsens=Isens;
When PETSc is used, I am unable to retrieve the displacements in the Inner Mesh.
Below is the code I tried to run :
```auto
include "Parameters.idp"
load "PETSc"
macro dimension()3 // EOM
include "macro_ddm.idp"
macro def(i) [i, i#2, i#3] // EOM
macro init(i)[i, i, i] // EOM
/////////// Optimization domain ///////////
// Boundary conditions
int free = 0, loaded = 1, fixed = 2, out = 4;
mesh3 ThF = readmesh3("Meshes/CantileverFull.mesh");
mesh3 ThI = readmesh3("Meshes/BeamInner.mesh");
plot(ThF);
// Finite element function
func Pk = [P1, P1, P1];
fespace VhF (ThF, Pk);
VhF [ux, uy, uz], [vx, vy, vz];
fespace VhF0 (ThF, P0);
VhF0 E, theta, Fsens;
fespace VhI0 (ThI, P0);
VhI0 Isens;
// Strain and stiffness tensor
macro e (x1,x2,x3) [dx(x1),dy(x2),dz(x3),(dz(x2)+dy(x3)),(dz(x1)+dx(x3)),(dy(x1)+dx(x2))] // EOM
macro D(Ei) (Ei*[[2*mu + lambda, lambda, lambda, 0, 0, 0], [lambda, 2*mu + lambda, lambda, 0, 0, 0], [lambda, lambda, 2*mu + lambda, 0, 0, 0], [0, 0, 0, mu, 0, 0], [0, 0, 0, 0, mu, 0], [0, 0, 0, 0, 0, mu]]) // EOM
// Density and elastic modulus field
theta = 1;
E = theta^p*(E1-E0) + E0;
// Variational form of Elasticity problem
varf Elasticity([ux, uy, uz], [vx, vy, vz])
= int3d(ThF) (E*((e(ux,uy,uz)'*D(E1)*e(vx,vy,vz))))
+ int3d(ThF) (theta*omega^2*Radius*vz)
+ int2d(ThF, loaded)(fx*vx + fy*vy + fz*vz)
+ on(fixed, ux=0, uy=0, uz=0);
Mat K;
createMat(ThF, K, Pk);
set(K, sparams="-pc_type lu");
real[int] F(VhF.ndof);
K = Elasticity(VhF, VhF);
F = Elasticity(0, VhF);
ux[] = K^-1 * F;
mpiBarrier(mpiCommWorld);
// Sensibility field on ThF
Fsens = p*theta^(p-1)*(E1-E0)*(e(ux, uy, uz)'*D(E1)*e(ux, uy, uz));
// Sensibility field on ThI
Isens = Fsens;
savevtk("Outputs/Isens_ThI.vtu", ThI, Isens);
```
Parameters.idp only contains scalar values relative to material properties and is attached to this post.
[Parameters.idp](https://community.freefem.org/uploads/short-url/rIjvCw8inUqB59KpgPaY33qJEEB.idp) (1.6 KB)
Meshes are attached too big to fit in the message, but you can find them here :
> **[SwissTransfer - Send large files securely and free of charge](https://www.swisstransfer.com/d/c3c61893-d6e8-4f29-a5d7-cb00407adc13)**
>
> Send up to 50 GB - Free and without registration - Keep your transfers for up to 30 days.
I am fairly new to parallel computing using freefem, any advice on what could cause an issue ?
PS : I also tried using the Transfer function or the VecInterpolate to not sucess (error in macro\_dmm.idp) (Following advice from [Copy variables between different meshes with PETSc](https://community.freefem.org/t/copy-variables-between-different-meshes-with-petsc/390))
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [October 21, 2025, 12:47pm UTC](https://community.freefem.org/t/petsc-and-mesh-interpolation/4109/2 "2025-10-21T12:47:40Z")
You likely want to rebuild a centralized solution `Fsens` using the `N2O` functionality of the `MatCreate()` macro. And then your `Isens = Fsens` will work.
---
### Author:  [@HMasson](https://community.freefem.org/u/HMasson)
#### Post date: [October 21, 2025, 1:33pm UTC](https://community.freefem.org/t/petsc-and-mesh-interpolation/4109/3 "2025-10-21T13:33:54Z")
Thank you for your quick response.
Do you have in mind an example in which this is implemented ? (There is no proper documentation on MatCreate or n2o in the online doc)
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [October 21, 2025, 1:35pm UTC](https://community.freefem.org/t/petsc-and-mesh-interpolation/4109/4 "2025-10-21T13:35:35Z")
There are numerous examples using it, e.g., [FreeFem-sources/examples/hpddm/laplace-adapt-3d-PETSc.edp at b1e524c8ca6a9b44cb9eff780459bacfc12a3ad7 · FreeFem/FreeFem-sources · GitHub](https://github.com/FreeFem/FreeFem-sources/blob/b1e524c8ca6a9b44cb9eff780459bacfc12a3ad7/examples/hpddm/laplace-adapt-3d-PETSc.edp#L22).