PETSc and Mesh interpolation

General Discussion
1

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 :

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 (1.6 KB)

Meshes are attached too big to fit in the message, but you can find them here :

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)

2

You likely want to rebuild a centralized solution Fsens using the N2O functionality of the MatCreate() macro. And then your Isens = Fsens will work.

3

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)

4

There are numerous examples using it, e.g., FreeFem-sources/examples/hpddm/laplace-adapt-3d-PETSc.edp at b1e524c8ca6a9b44cb9eff780459bacfc12a3ad7 · FreeFem/FreeFem-sources · GitHub.

1 Like