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: overlapping domain decomposition     *
 8  *                utility routines                     *
 9  *******************************************************/
 
11 assert(mpisize == 2);
12 load "hpddm"
13 macro dimension()2// EOM
14 include "macro_ddm.idp"
 
16 int[int] l = [2, 2, 2, 2];
17 mesh T = square(5, 6, label = l, [5.0/6.0 * x, y]);
18 mesh Ttop = trunc(T, x < 1.0/3.0, label = 2);
19 Ttop = trunc(Ttop, y > 0.7, label = 10);
20 Ttop = movemesh(Ttop, [1.0/3.0 - x, y]);
21 mesh Tt = trunc(T, !(y > 0.7 || (x > 1.0/3.0 && y < 0.3)), label = 2);
 
23 macro grad(u)[dx(u), dy(u)]// EOM
24 meshN Th = Ttop + Tt;
25 fespace Vh(Th, P1);
26 int[int][int] intersection;
27 real[int] D;
28 {
29     fespace Ph(Th, P0);
30     Ph part = x < 1.0/3.0 ? 1 : 0;
31     buildWithPartitioning(Th, part[], 1, intersection, D, P1, mpiCommWorld);
32 }
 
34 varf vPb(u, v) = intN(Th)(grad(u)' * grad(v)) + intN(Th)(v) + on(2, u = 1.0);
35 matrix Loc = vPb(Vh, Vh, tgv = -1);
36 real[int] b = vPb(0, Vh, tgv = -1);
 
38 schwarz A(Loc, intersection, D);
39 Vh x;
40 x[] = A^-1 * b;
41 macro def(i)i// EOM
42 plotMPI(Th, x, P1, def, real, cmm = "Solution");
43 statistics(A);
 
45 Vh u;
46 u[] = mpirank;
47 plotMPI(Th, u, P1, def, real, cmm = "Not consistent function");
48 exchange(A, u[], scaled = true);
49 plotMPI(Th, u, P1, def, real, cmm = "Consistent function");
50 u[] = D;
51 // exchange(A, u[]); // equivalent to the line below
52 exchange(A, u[], scaled = false);
53 plotMPI(Th, u, P1, def, real, cmm = "Partition of unity equal to 1?");