1 /*******************************************************
 2  *  This file is part of the FreeFEM tutorial made by  *
 3  *      Pierre Jolivet <pierre@joliv.et>               *
 4  *                                                     *
 5  *  See https://joliv.et/FreeFem-tutorial for more     *
 6  *                                                     *
 7  *   Description: MUMPS distributed-memory             *
 8  *                parallel solver                      *
 9  *******************************************************/
 
11 load "MUMPS"
12 load "hpddm"
13 load "metis"
14 include "cube.idp"
15 macro grad(u)[dx(u), dy(u), dz(u)]// EOM
 
17 int[int] LL = [1,2, 1,1, 1,1];
18 mesh3 cb = cube(10, 10, 10, [x, y, z], label = LL);
19 fespace Ph(cb, P0);
20 Ph part;
21 if(mpisize > 1)
22     metisdual(part[], cb, mpisize);
23 else
24     part[] = 0.0;
25 cb = change(cb, fregion = part[][nuTriangle]);
26 fespace Vh(cb, P2);
27 real dt = 0.001;
28 varf vPb(u, v) = int3d(cb, mpirank)(u * v + dt * (grad(u)' * grad(v))) + int3d(cb, mpirank)(dt * cos(x) * v) + on(1, u = 1.0);
29 varf vPbM(u, v) = int3d(cb, mpirank)(u * v);
30 matrix A = vPb(Vh, Vh);
31 matrix M = vPbM(Vh, Vh);
32 Vh b, u;
33 b[] = vPb(0, Vh);
34 set(A, solver = sparsesolver, master = -1);
35 u[] = 0.0;
36 int iMax = 20;
37 for(int i = 0; i < iMax; ++i) {
38     real[int] inc = b[];
39     inc += M * u[];
40     u[] = A^-1 * inc;
41     plot(u, wait = 1);
42     int[int] fforder(1);
43     fforder = 1;
44     mesh3 cbLocal = trunc(cb, (abs(part - mpirank) < 1.0e-12 ? 1 : 0));
45     fespace VhLocal(cbLocal, P2);
46     VhLocal exportU = u;
47     savevtk("Temperature.vtu", cbLocal, exportU, dataname = "Temperature", order = fforder, append = i ? true : false);
48 }