Parameters for macro_ddm.idp

General Discussion
1

Hi FreeFem developers!

I use FreeFem to implement shape optimization for nonlinear mechanics and I am trying make the entire solver parallel. To begin with, I started with the example elasticity-3d_PETSc.edp from FreeFem examples. The first few lines read

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

…
Mat A;
buildMat(Th, getARGV("-split", 1), A, Pk, mpiCommWorld, 3)

This works perfectly fine and I have the correct displacement solution. However, I want the store the global solution and in one of your previous posts, I read that I should call the restrict operator, for which I must use
buildDmesh(Th)
createMat(Th, A, Pk);
instead of buildmat.
When I do that, I get an error

build(Th, 1, intersection, privateDmeshThkhi[0], P1, ThComm)Error in macro expansion build, we wait for , and we get  )
 number of macro parameter in definition is 7

 	0	 in elasticity-3d-PETSc.edp line : 21
 	2	 in  macro: buildDmesh in /usr/local/lib/ff++/4.7-1/idp/macro_ddm.idp line : 470
  current line = 469 mpirank 0 / 4
Compile error :  Wrong number of parameter in  macro call
	line number :469, )
error Compile error :  Wrong number of parameter in  macro call

However, if I declare
macro vectorialfe()P1// EOM
after including the macro_ddm.idp instead of declaring it before, buildDmesh and createmat work fine. But when I plot the displacement u, it seems erroneous.

Kindly help me troubleshoot this issue and to understand macro_ddm.idp.
If this query has been addressed earlier somewhere, I apologize.

2

You don’t really need vectorialfe(), it’s from the old API (I know realize that elasticity-3d-PETSc.edp uses the old API, my apologizes…).
Here is a fixed script that should do what you want.

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

macro def(i)[i, i#B, i#C]// EOM     // vector field definition
macro init(i)[i, i, i]// EOM        // vector field initialization
real Sqrt = sqrt(2.0);
macro epsilon(u)[dx(u), dy(u#B), dz(u#C), (dz(u#B) + dy(u#C)) / Sqrt, (dz(u) + dx(u#C)) / Sqrt, (dy(u) + dx(u#B)) / Sqrt]// EOM
macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM
func Pk = [P2, P2, P2]; // finite element space


int[int] LL = [2,3, 2,1, 2,2];
mesh3 Th = cube(10 * getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5), [10 * x, y, z], label = LL);
Mat A;
int[int] n2o;
macro ThN2O()n2o// EOM
buildDmesh(Th);
createMat(Th, A, Pk)

real f = -9000.0;
real strain = 100.0;
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));
varf vPb(def(u), def(v)) = int3d(Th)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + int3d(Th)(f * vC) + on(1, u = 0.0, uB = 0.0, uC = 0.0); // '
fespace Wh(Th, 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);
Wh<real> def(u);

u[] = A^-1 * rhs;
if(!NoGraphicWindow) {
    real[int] sol;
    ChangeNumbering(A, u[], sol);
    ChangeNumbering(A, u[], sol, inverse = true);
    mesh3 ThPlt = cube(10 * getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5), [10 * x, y, z], label = LL);
    fespace WhPlt(ThPlt, Pk);
    WhPlt [pltx, plty, pltz];
    WhPlt [reducex, reducey, reducez];
    int[int] rest = restrict(Wh, WhPlt, n2o);
    for[i, v : rest] pltx[][v] = u[][i];
    mpiReduce(pltx[], reducex[], processor(0, mpiCommWorld), mpiSUM);
    if(mpirank == 0)
        medit("Global solution", ThPlt, [real(reducex), real(reducey), real(reducez)]);
}
3

This script works perfectly fine for me.
Thanks a lot for your help!!