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: DmeshCreate with N2O and restrict     *
 8  *******************************************************/
  
10 load "PETSc"
11 macro dimension()2// EOM
12 include "macro_ddm.idp"
 
14 border a(t=0, 2*pi) { x=cos(t); y=sin(t); label=1; }
15 mesh ThGlobal;
16 if(!mpirank)
17     ThGlobal = buildmesh(a(80));
18 broadcast(processor(0), ThGlobal);
 
20 int[int] n2o;
21 mesh Th = ThGlobal;
22 macro ThN2O()n2o// EOM
23 DmeshCreate(Th);
 
25 func Pk = P2; // consistent finite element space
26 fespace Vh(Th, Pk); // full global mesh
27 fespace VhGlobal(ThGlobal, Pk); // distributed local mesh
 
29 Vh u = x^2 + y^2;
30 int[int] R = restrict(Vh, VhGlobal, n2o);
31 VhGlobal uGlobal, uReduce;
32 uGlobal[](R) = u[];
33 mpiReduce(uGlobal[], uReduce[], processor(0, mpiCommWorld), mpiSUM);
34 plot(uReduce, wait = 1, cmm = "Wrong u on the global mesh (redundant degrees of freedom)");
 
36 real[int] D;
37 PartitionCreate(Th, D, Pk);
38 u[] .*= D;
39 uGlobal[](R) = u[];
40 mpiReduce(uGlobal[], uReduce[], processor(0, mpiCommWorld), mpiSUM);
41 plot(uReduce, wait = 1, cmm = "u on the global mesh (weighted before summed)");
 
43 Mat A;
44 MatCreate(Th, A, Pk);
45 u = x^2 + y^2;
46 u[] .*= A.D;
47 uGlobal[](R) = u[];
48 mpiReduce(uGlobal[], uReduce[], processor(0, mpiCommWorld), mpiSUM);
49 plot(uReduce, wait = 1, cmm = "u on the global mesh (weighted before summed by the partition of unity attached to the Mat)");
 
51 u = x^2 + y^2;
52 real[int] tmp;
53 ChangeNumbering(A, u[], tmp);
54 ChangeNumbering(A, u[], tmp, inverse = true);
55 uGlobal[](R) = u[];
56 mpiReduce(uGlobal[], uReduce[], processor(0, mpiCommWorld), mpiSUM);
57 plot(uReduce, wait = 1, cmm = "u on the global mesh (using PETSc unique numbering)");
 
59 uGlobal = x^2 + y^2;
60 u[] = uGlobal[](R);
61 plotD(Th, u, cmm = "u on each local mesh");