# Iterative parallel computing with PETSc

Hi FreeFem++ developers!

Thank you very much for using FF.
I’d like to repeat the parallel computation using PETSc for implementing topology optimization.
So, I begin with the example elasticity-3d_PETSc.edp from FreeFem examples.

At first, I repeated the same analysis with two different meshes, run with the following code;

//Begin Program

macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

mesh3		ShPlt=Sh;
mesh3 		ShPltMoved;
plot(Sh,cmm="FirstMesh");

// Define Boundary Label
int[int] wall(3);     		wall    	= [1,7,13];			// 	fixed displacement
int[int] traction(1);     	traction    = [9];				//	traction

//
int[int] n2o;
macro ShN2O()n2o// EOM
func Pk = [P2, P2, P2]; // finite element space
fespace Wh(Sh, Pk);
fespace WhPlt(ShPlt, Pk);

// Parameters
real Young = 1.0e8;
real poisson = 0.45;
real tmp = 1.0 + poisson;
real mu = Young  / (2.0 * tmp);
real lambda = Young * poisson / (tmp * (1.0 - 2.0 * poisson));

// Define Solid Mechanics
macro def(i)[i, i#B, i#C]// EOM     // vector field definition
macro init(i)[i, i, i]// EOM        // vector field initialization
macro gu 	[gu1,gu2,gu3] 				// Glabal	displacement vector
macro epsilon(u) 	[dx(u),dy(u#B),dz(u#C),(dz(u#B)+dy(u#C)), (dz(u)+dx(u#C)), (dy(u)+dx(u#B))] // strain tensor
macro D [
[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,				2.*mu+lambda,	0,				0			],
[0,				0,				0,				0,				2.*mu+lambda,	0			],
[0,				0,				0,				0,				0,				2.*mu+lambda]
] //elastic tensor
//macro g 	[0,0,-1000000.0] 				// traction vector
real gZ;	gZ = -1000000.0;
macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM

WhPlt [pltx, plty, pltz];
WhPlt [reducex, reducey, reducez];
WhPlt [gu1, gu2, gu3];
Wh def(u);
real[int] sol;
int[int] rest;

varf vPb(def(u), def(v)) = int3d(Sh)((D*epsilon(u))'*epsilon(v))
+int2d(Sh,traction)(gZ*vC)
+ on(wall, u = 0.0, uB = 0.0, uC = 0.0); // '
Mat A;

//---------------------------------1stAnalysis---------------------------------

buildDmesh(Sh);
createMat(Sh, A, Pk)

A = vPb(Wh, Wh);
real[int] rhs = vPb(0, Wh);
set(A, sparams = "-pc_type cholesky -ksp_view -ksp_max_it 100", bs = 3);
u[] = A^-1 * rhs;

if(!NoGraphicWindow) {
ChangeNumbering(A, u[], sol);
ChangeNumbering(A, u[], sol, inverse = true);
rest = restrict(Wh, WhPlt, n2o);
for[i, v : rest] pltx[][v] = u[][i];
mpiReduce(pltx[], reducex[], processor(0, mpiCommWorld), mpiSUM);
}
gu1[]=reducex[];
gu2[]=reducey[];
gu3[]=reducez[];
ShPltMoved = movemesh3(ShPlt,transfo=[x+gu1,y+gu2,z+gu3]);
plot(ShPltMoved,cmm="Deformation1");

//-----------------------------------------------------------------------------

ShPlt=Sh;
plot(Sh,cmm="SecondMesh");

//---------------------------------2ndAnalysis----------------------------------

buildDmesh(Sh);
createMat(Sh, A, Pk)

A = vPb(Wh, Wh);
rhs = vPb(0, Wh);
set(A, sparams = "-pc_type cholesky -ksp_view -ksp_max_it 100", bs = 3);
u[] = A^-1 * rhs;

if(!NoGraphicWindow) {
ChangeNumbering(A, u[], sol);
ChangeNumbering(A, u[], sol, inverse = true);
rest = restrict(Wh, WhPlt, n2o);
for[i, v : rest] pltx[][v] = u[][i];
mpiReduce(pltx[], reducex[], processor(0, mpiCommWorld), mpiSUM);
}
gu1[]=reducex[];
gu2[]=reducey[];
gu3[]=reducez[];
ShPltMoved = movemesh3(ShPlt,transfo=[x+gu1,y+gu2,z+gu3]);
plot(ShPltMoved,cmm="Deformation2");

//End Program

When I do that, I get an error.

current line = 96 mpirank 3 / 4
current line = 96 mpirank 0 / 4
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** current line = 96 mpirank 2 / 4**
** current line = 96 mpirank 1 / 4**
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** err code 8 , mpirank 0**
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** err code 8 , mpirank 3**
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** err code 8 , mpirank 1**
Exec error : AssembleLinearForm size rhs and nb of DF of Vh
** – number :1**
** err code 8 , mpirank 2**

I confirmed that multiple analysis runs can be performed on the same mesh.
I don’t know what procedures are required in case of a mesh change.
I’m very sorry if a similar question has already been asked and this is a rudimentary question.
I am sure you are very busy, but I would appreciate your guidance.

We cannot copy/paste your code. Make sure you use preformatted text.