1include "getARGV.idp"
2IFMACRO(!partitioner)
3macro partitioner()metis// EOM
4ENDIFMACRO
5IFMACRO(partitioner,metis)
6load "metis"
7macro partitionerSeq(part, Th, size){ if(size <= 1) part = 0; else metisdual(part, Th, size); }// EOM
8macro partitionerPar(part, Th, comm, size)broadcast(processor(0, comm), part)// EOM
9ENDIFMACRO
10IFMACRO(partitioner,scotch)
11load "scotch"
12macro partitionerSeq(part, Th, size)scotch(part, Th, size)// EOM
13macro partitionerPar(part, Th, comm, size)broadcast(processor(0, comm), part)// EOM
14ENDIFMACRO
15IFMACRO(partitioner,parmetis)
16load "parmetis"
17macro partitionerSeq(part, Th, size)// EOM
18macro partitionerPar(part, Th, comm, size)parmetis(part, Th, size, communicator = comm, worker = getARGV("-parmetis_worker", 1))// EOM
19ENDIFMACRO
20IFMACRO(!partitionerSeq)
21cout << "The macro 'partitioner' must be set to 'metis', 'scotch', or 'parmetis'" << endl;
22exit(1);
23ENDIFMACRO
24IFMACRO(!dimension)
25cout << "The macro 'dimension' must be defined" << endl;
26exit(1);
27ENDIFMACRO
28IFMACRO(dimension,2)
29macro meshN()mesh// EOM             // two-dimensional problem
30macro intN()int2d// EOM             // two-dimensional integral
31macro intN1()int1d// EOM            // one-dimensional integral
32macro readmeshN()readmesh// EOM     // two-dimensional problem
33ENDIFMACRO
34IFMACRO(dimension,3)
35load "msh3"
36macro meshN()mesh3// EOM            // three-dimensional problem
37macro intN()int3d// EOM             // three-dimensional integral
38macro intN1()int2d// EOM            // two-dimensional integral
39macro readmeshN()readmesh3// EOM    // three-dimensional problem
40ENDIFMACRO
41IFMACRO(dimension,3S)
42load "msh3"
43macro meshN()meshS// EOM            // three-dimensional surface problem
44macro intN()int2d// EOM             // two-dimensional integral
45macro intN1()int1d// EOM            // one-dimensional integral
46macro intNxN()int2dx2d// EOM        // two-dimensional integral for BEM
47ENDIFMACRO
48IFMACRO(dimension,3L)
49load "msh3"
50macro meshN()meshL// EOM            // three-dimensional line problem
51macro intN()int1d// EOM             // one-dimensional integral
52macro intN1()int0d// EOM            // zero-dimensional integral
53macro intNxN()int1dx1d// EOM        // one-dimensional integral for BEM
54ENDIFMACRO
55IFMACRO(!meshN)
56cout << "The macro 'dimension' must be set to '2', '3', '3S' or '3L'" << endl;
57exit(1);
58ENDIFMACRO
59
60macro plotDmesh(Th, params)
61if(!NoGraphicWindow || usedARGV("-fglut") != -1) {
62    fespace PhPlot(Th, P0);
63    PhPlot plt;
64    if(Th.nt)
65        plt[] = mpirank;
66NewMacro defPlt#Th(u)u EndMacro
67    plotMPI(Th, plt, P0, defPlt#Th, real, params)
68}//
69
70macro plotD(Th, u, params)
71if(!NoGraphicWindow || usedARGV("-fglut") != -1) {
72    fespace VhPlot(Th, P1);
73    VhPlot plt;
74    if(Th.nt)
75        plt = u;
76NewMacro defPlt#Th(v)v EndMacro
77    plotMPI(Th, plt, P1, defPlt#Th, real, params)
78}//
79
80macro plotMPI(Th, u, Pk, def, K, params)
81if(!NoGraphicWindow || usedARGV("-fglut") != -1) {
82    meshN ThCurrent = Th;
83    fespace Xh(ThCurrent, Pk);
84    Xh<K> def(uSend);
85    if(ThCurrent.nt)
86        def(uSend) = u;
87    if(mpirank == 0) {
88        meshN[int] meshTab(mpisize);
89        Xh<K>[int] def(uTab)(mpisize);
90        if(ThCurrent.nt)
91            uTab[0][] = uSend[];
92        meshTab[0] = ThCurrent;
93        mpiRequest[int] rq(mpisize - 1);
94        for(int i = 1; i < mpisize; ++i)
95            Irecv(processor(i, mpiCommWorld, rq[i - 1]), meshTab[i]);
96        mpiWaitAll(rq);
97        for(int i = 1; i < mpisize; ++i) {
98            ThCurrent = meshTab[i];
99            if(ThCurrent.nt)
100                Irecv(processor(i, mpiCommWorld, rq[i - 1]), uTab[i][]);
101        }
102        mpiWaitAll(rq);
103        plot(def(uTab), params);
104    }
105    else {
106        mpiRequest[int] rq(2);
107        Isend(processor(0, rq[0]), ThCurrent);
108        if(ThCurrent.nt)
109            Isend(processor(0, rq[1]), uSend[]);
110        mpiWaitAll(rq);
111    }
112}// EOM
113
114macro partition(meshName, borderName, globalName, PhGlobal, VhGlobal, part, rank, size, s, overlap, level, prolongation, D, P, intersection, comm, fakeInterface, PkPart, defPart, initPart, bs) {
115    int backupSM = searchMethod;
116    searchMethod = 1;
117    assert(level >= 1);
118IFMACRO(!privateCreatePartition)
119IFMACRO(!privateCreateMat)
120    meshName[level - 1] = trunc(globalName, abs(part - rank) < 0.1, label = fakeInterface);
121    intersection.resize(1);
122    intersection[0].resize(0);
123    PhGlobal supp = abs(part - rank) < 0.1;
124    VhGlobal suppSmooth;
125    AddLayers(globalName, supp[], 2 * overlap, suppSmooth[]);
126    Unique(part[], intersection[0], remove = rank);
127    fespace VhLocal(meshName[level - 1], P1);
128    VhLocal[int] partitionIntersection(intersection[0].n);
129    if(s > 1 && level <= 1) {
130        globalName = trunc(globalName, suppSmooth > 0.001, split = s);
131        supp = abs(part - rank) < 0.1;
132        suppSmooth = 0;
133        AddLayers(globalName, supp[], 2 * overlap, suppSmooth[]);
134    }
135IFMACRO(!privateDmesh#N2O)
136    globalName = trunc(globalName, suppSmooth > 0.001, label = 9999);
137ENDIFMACRO
138    real eps = globalName.measure;
139    real[int] epsTab(intersection[0].n);
140    mpiRequest[int] rq(2 * intersection[0].n);
141    if(mpiSize(comm) == size) {
142        for(int j = 0; j < intersection[0].n; ++j)
143            Irecv(processor(intersection[0][j], comm, rq[j]), epsTab[j]);
144        for(int j = 0; j < intersection[0].n; ++j)
145            Isend(processor(intersection[0][j], comm, rq[intersection[0].n + j]), eps);
146    }
147    else
148        epsTab = 1.0e+30;
149IFMACRO(!privateDmesh#N2O)
150    suppSmooth = suppSmooth;
151    meshName[level - 1] = trunc(globalName, suppSmooth > 0.501, label = fakeInterface);
152ENDIFMACRO
153IFMACRO(privateDmesh#N2O)
154    meshName[level - 1] = trunc(globalName, suppSmooth > 0.501, label = fakeInterface, new2old = privateDmesh#N2O);
155    globalName = trunc(globalName, suppSmooth > 0.001, label = 9999);
156    suppSmooth = suppSmooth;
157ENDIFMACRO
158    if(level > 1) {
159        prolongation.resize(level - 1);
160        if(s > 1) {
161            meshN globalNameRefined = globalName;
162            for(int i = level - 1; i > 0; --i) {
163                globalNameRefined = trunc(globalNameRefined, 1, split = s);
164                meshName[i - 1] = trunc(globalNameRefined, suppSmooth > 0.501, label = fakeInterface);
165                fespace WhLocalRefined(meshName[i - 1], P);
166                fespace WhLocalCoarse(meshName[i], P);
167                prolongation[i - 1] = interpolate(WhLocalRefined, WhLocalCoarse);
168            }
169        }
170        else
171            for(int i = level - 1; i > 0; --i)
172                meshName[i - 1] = meshName[i];
173    }
174    if(!removeZeros && (fakeInterface != -111111 || overlap != 1)) {
175        if(suppSmooth[].min < 0.501) {
176            supp = supp;
177            borderName[level - 1] = trunc(globalName, (suppSmooth > (overlap - 0.999) / real(2 * overlap)) && (suppSmooth < 0.501), label = (abs(fakeInterface) + 1) * 100);
178            if(s > 1)
179                for(int i = level - 2; i >= 0; --i) {
180                    borderName[i] = trunc(borderName[i + 1], 1, split = s, label = (abs(fakeInterface) + 1) * 100);
181                    meshN tempRefined = meshName[i] + borderName[i];
182                    fespace PhRefined(tempRefined, P0);
183                    PhRefined suppRefined = supp;
184                    fespace VhBorderRefined(borderName[i], P1);
185                    VhBorderRefined suppBorder = suppRefined;
186                    borderName[i] = trunc(borderName[i], suppBorder > 0.01);
187                }
188            else
189                for(int i = level - 2; i >= 0; --i)
190                    borderName[i] = borderName[i + 1];
191        }
192    }
193    VhLocal khi = max(2 * suppSmooth - 1.0, 0.0);
194    VhLocal sum = khi;
195    VhGlobal phi;
196    part = part;
197    int numberIntersection = 0;
198    mpiWaitAll(rq);
199    for(int i = 0; i < intersection[0].n; ++i) {
200        PhGlobal suppPartition = abs(part - intersection[0][i]) < 0.1;
201        AddLayers(globalName, suppPartition[], overlap, phi[]);
202        if(min(eps, epsTab[i]) > 0.0) {
203            if(intN(globalName)(phi) / min(eps, epsTab[i]) > 1.0e-10) {
204                partitionIntersection[numberIntersection] = phi;
205                if(!trueRestrict)
206                    sum[] += partitionIntersection[numberIntersection][];
207                intersection[0][numberIntersection++] = intersection[0][i];
208            }
209        }
210    }
211    if(numberIntersection != intersection[0].n)
212        intersection[0].resize(numberIntersection);
213    intersection.resize(1 + level * numberIntersection);
214ENDIFMACRO
215IFMACRO(privateCreateMat)
216    assert(level == 1);
217    int numberIntersection = privateDmesh#meshName#intersectionDef.n - 1;
218    intersection.resize(1 + level * numberIntersection);
219    intersection[0].resize(numberIntersection);
220    fespace VhLocal(meshName[level - 1], P1);
221    VhLocal[int] partitionIntersection(numberIntersection);
222    for(int i = 0; i < numberIntersection; ++i) {
223        intersection[0][i] = privateDmesh#meshName#intersectionDef[0][i];
224        partitionIntersection[i][] = privateDmesh#meshName#intersectionDef[1 + i];
225    }
226ENDIFMACRO
227IFMACRO(privateBuildDmesh)
228    privateDmesh#meshName#intersectionDef.resize(1 + numberIntersection);
229    privateDmesh#meshName#intersectionDef[0].resize(numberIntersection);
230    for(int i = 0; i < numberIntersection; ++i) {
231        privateDmesh#meshName#intersectionDef[0][i] = intersection[0][i];
232        privateDmesh#meshName#intersectionDef[1 + i].resize(VhLocal.ndof);
233        privateDmesh#meshName#intersectionDef[1 + i] = partitionIntersection[i][];
234    }
235ENDIFMACRO
236    meshN[int] meshIntersection(numberIntersection);
237    for(int j = 0; j < (s == 1 ? 1 : level); ++j) {
238        for(int i = 0; i < numberIntersection; ++i) {
239            int[int] n2o;
240            meshIntersection[i] = trunc(meshName[j], partitionIntersection[i] > 1.0e-6, new2old = n2o, label = 9999);
241IFMACRO(!privateCreateMat)
242            if(!removeZeros)
243ENDIFMACRO
244            {
245IFMACRO(vectorialfe)
246                fespace singleComponentWh(meshName[j], vectorialfe);
247                fespace WhIntersection(meshIntersection[i], vectorialfe);
248ENDIFMACRO
249IFMACRO(!vectorialfe)
250                fespace singleComponentWh(meshName[j], P);
251                fespace WhIntersection(meshIntersection[i], P);
252ENDIFMACRO
253                intersection[1 + i + j * numberIntersection] = restrict(WhIntersection, singleComponentWh, n2o);
254            }
255        }
256    }
257IFMACRO(!privateCreateMat)
258    if(s == 1 && level > 1 && !removeZeros)
259        for(int j = 1; j < level; ++j)
260            for(int i = 0; i < numberIntersection; ++i) {
261                intersection[1 + i + j * numberIntersection].resize(intersection[1 + i].n);
262                intersection[1 + i + j * numberIntersection] = intersection[1 + i];
263            }
264    partitionIntersection.resize(0);
265    for(int i = 0; i < (trueRestrict ? level : level - 1); ++i) {
266        fespace VhRefined(meshName[i], P1);
267        fespace PhRefined(meshName[i], P0);
268        PhRefined partRefined = part;
269        PhRefined supp = abs(partRefined - rank) < 0.1;
270        varf vSupp(u, v) = intN(meshName[i], qforder = 1)(supp * v);
271        VhRefined khiL;
272        khiL[] = vSupp(0, VhRefined);
273        khiL = khiL > 0.0;
274        VhRefined sum = khiL;
275        for(int j = 0; j < numberIntersection; ++j) {
276            supp = abs(partRefined - intersection[0][j]) < 0.1;
277            VhRefined phiL;
278            phiL[] = vSupp(0, VhRefined);
279            phiL = phiL > 0.0;
280            sum[] += phiL[];
281        }
282        khiL[] ./= sum[];
283        if(i < level - 1) {
284            fespace WhRefined(meshName[i], PkPart);
285            WhRefined defPart(func2vec);
286            defPart(func2vec) = initPart(khiL);
287            D[i].resize(WhRefined.ndof);
288            D[i] = func2vec[];
289        }
290        else
291            khi[] = khiL[];
292    }
293    if(!trueRestrict)
294        khi[] = khi[] ./= sum[];
295    if(trueRestrict && mpiSize(comm) == size && removeZeros) {
296        assert(level == 1);
297        meshN ThIntersection;
298        fespace PhIntersection(ThIntersection, P0);
299        PhIntersection[int] recv(numberIntersection);
300        PhIntersection[int] send(numberIntersection);
301        mpiRequest[int] rq(2 * numberIntersection);
302        for(int i = 0; i < numberIntersection; ++i) {
303            ThIntersection = meshIntersection[i];
304            Irecv(processor(intersection[0][i], comm, rq[i]), recv[i][]);
305            send[i] = khi;
306            Isend(processor(intersection[0][i], comm, rq[numberIntersection + i]), send[i][]);
307        }
308        meshName[0] = trunc(meshName[0], khi > 1.0e-6, label = 9999);
309        khi = khi;
310        int[int] skip(0);
311        for(int k = 0; k < 2 * numberIntersection; ++k) {
312            int i = mpiWaitAny(rq);
313            if(i < numberIntersection) {
314                ThIntersection = meshIntersection[i];
315                PhIntersection intersection = send[i] > 1.0e-6 && recv[i] > 1.0e-6;
316                if(intersection[].l2 > 1.0e-6)
317                    meshIntersection[i] = trunc(meshIntersection[i], intersection  > 1.0e-6, label = 9999);
318                else {
319                    skip.resize(skip.n + 1);
320                    skip[skip.n - 1] = i;
321                }
322            }
323        }
324        skip.sort;
325        intersection.resize(1 + numberIntersection - skip.n);
326        int j = 0;
327        for(int i = 0; i < numberIntersection; ++i) {
328            bool skipped = false;
329            if(j < skip.n) {
330                if(skip[j] == i) {
331                    ++j;
332                    skipped = true;
333                }
334            }
335            if(!skipped) {
336IFMACRO(vectorialfe)
337                fespace singleComponentWh(meshName[0], vectorialfe);
338                fespace WhIntersection(meshIntersection[i], vectorialfe);
339ENDIFMACRO
340IFMACRO(!vectorialfe)
341                fespace singleComponentWh(meshName[0], P);
342                fespace WhIntersection(meshIntersection[i], P);
343ENDIFMACRO
344                matrix meshName#R = interpolate(WhIntersection, singleComponentWh);
345                meshName#R.thresholding(1.0e-10);
346                real[int] meshName#C;
347                int[int] meshName#I;
348                [meshName#I, intersection[1 + i - j], meshName#C] = meshName#R;
349                intersection[1 + i - j].resize(meshName#R.nbcoef);
350                intersection[0][i - j] = intersection[0][i];
351            }
352        }
353        numberIntersection -= skip.n;
354        intersection[0].resize(numberIntersection);
355        if(fakeInterface != -111111 || overlap != 1) {
356            PhGlobal suppPartition = khi > 0.1;
357            AddLayers(globalName, suppPartition[], 1, phi[]);
358            borderName[0] = trunc(globalName, phi > 0.001 && phi < 0.501, label = (abs(fakeInterface) + 1) * 100);
359        }
360    }
361ENDIFMACRO
362IFMACRO(vectorialfe)
363    if(bs > 1)
364        for(int i = 0; i < intersection.n - 1; ++i) {
365            int n = intersection[1 + i].n;
366            intersection[1 + i].resize(n * bs);
367            for(int j = n - 1; j != -1; --j)
368                for(int k = bs - 1; k != -1; --k)
369                    intersection[1 + i][j * bs + k] = intersection[1 + i][j] * bs + k;
370        }
371ENDIFMACRO
372ENDIFMACRO
373IFMACRO(privateCreatePartition)
374    fespace VhLocal(meshName[level - 1], P1);
375IFMACRO(!privateCreateMat)
376    VhLocal khi;
377ENDIFMACRO
378ENDIFMACRO
379IFMACRO(privateCreateMat)
380    VhLocal khi;
381    khi[] = privateDmesh#meshName#khiDef[0];
382ENDIFMACRO
383    fespace WhPart(meshName[level - 1], PkPart);
384    WhPart defPart(func2vec);
385    D[level - 1].resize(WhPart.ndof);
386    if((WhPart.ndof % meshName[level - 1].nt) == 0) {
387IFMACRO(privateCreateMat)
388        fespace PhLocal(meshName[level - 1], P0);
389        PhLocal partLocal;
390        partLocal[] = privateDmesh#meshName#khiDef[1];
391        defPart(func2vec) = initPart(abs(partLocal - rank) < 0.1);
392ENDIFMACRO
393IFMACRO(!privateCreateMat)
394        defPart(func2vec) = initPart(abs(part - rank) < 0.1);
395ENDIFMACRO
396    }
397    else if(WhPart.ndof == meshName[level - 1].nv) {
398        func2vec[] = khi[];
399    }
400    else {
401        defPart(func2vec) = initPart(khi);
402    }
403    D[level - 1] = func2vec[];
404IFMACRO(!privateCreatePartition)
405IFMACRO(!privateCreateMat)
406IFMACRO(privateBuildDmesh)
407    fespace PhLocal(meshName[level - 1], P0);
408    PhLocal partLocal;
409    partLocal = part;
410    privateDmesh#meshName#khiDef[1].resize(partLocal[].n);
411    privateDmesh#meshName#khiDef[1] = partLocal[];
412ENDIFMACRO
413ENDIFMACRO
414ENDIFMACRO
415    searchMethod = backupSM;
416}// EOM
417
418macro saveDmesh(ThName, name)
419IFMACRO(privateDmesh#ThName)
420{
421IFMACRO(!ThName#Comm)
422NewMacro ThName#Comm()mpiCommWorld EndMacro
423ENDIFMACRO
424IFMACRO(dimension,3)
425savemesh(ThName, name + "_" + mpiRank(ThName#Comm) + "_" + mpiSize(ThName#Comm) + ".meshb");
426ENDIFMACRO
427IFMACRO(dimension,2)
428savemesh(ThName, name + "_" + mpiRank(ThName#Comm) + "_" + mpiSize(ThName#Comm) + ".msh");
429ENDIFMACRO
430ofstream khi(name + "_" + mpiRank(ThName#Comm) + "_" + mpiSize(ThName#Comm) + ".khi");
431khi << privateDmesh#ThName#khi << endl;
432khi << privateDmesh#ThName#intersection << endl;
433IFMACRO(ThName#N2O)
434khi << ThName#N2O << endl;
435ENDIFMACRO
436}
437ENDIFMACRO
438IFMACRO(!privateDmesh#ThName)
439assert(0);
440ENDIFMACRO
441EndMacro
442
443macro loadDmesh(ThName, name)
444IFMACRO(!privateDmesh#ThName)
445NewMacro privateDmesh#ThName()privateDmesh#ThName EndMacro
446NewMacro privateDmesh#ThName#khi()privateDmesh#ThName#khiDef EndMacro
447NewMacro privateDmesh#ThName#intersection()privateDmesh#ThName#intersectionDef EndMacro
448real[int][int] privateDmesh#ThName#khi(2);
449real[int][int] privateDmesh#ThName#intersection;
450ENDIFMACRO
451{
452IFMACRO(!ThName#Comm)
453NewMacro ThName#Comm()mpiCommWorld EndMacro
454ENDIFMACRO
455IFMACRO(dimension,3)
456ThName = readmesh3(name + "_" + mpiRank(ThName#Comm) + "_" + mpiSize(ThName#Comm) + ".meshb");
457ENDIFMACRO
458IFMACRO(dimension,2)
459ThName = readmesh(name + "_" + mpiRank(ThName#Comm) + "_" + mpiSize(ThName#Comm) + ".msh");
460ENDIFMACRO
461privateDmesh#ThName#khi.resize(2);
462privateDmesh#ThName#khi[0].resize(ThName.nv);
463privateDmesh#ThName#khi[1].resize(ThName.nt);
464if(mpiSize(ThName#Comm) > 1) {
465    ifstream khi(name + "_" + mpiRank(ThName#Comm) + "_" + mpiSize(ThName#Comm) + ".khi");
466    int m;
467    khi >> m;
468    assert(m == 2);
469    khi >> privateDmesh#ThName#khi[0];
470    khi >> privateDmesh#ThName#khi[1];
471    khi >> m;
472    privateDmesh#ThName#intersection.resize(m);
473    for(int j = 0; j < m; ++j) {
474        int n;
475        khi >> n;
476        real[int] test;
477        privateDmesh#ThName#intersection[j].resize(n);
478        for[i, v : privateDmesh#ThName#intersection[j]]
479            khi >> v;
480    }
481IFMACRO(ThName#N2O)
482    ThName#N2O.resize(ThName.nt);
483    khi >> ThName#N2O;
484ENDIFMACRO
485}
486else {
487    privateDmesh#ThName#khi[0] = 1.0;
488    privateDmesh#ThName#khi[1] = 1.0;
489}
490}
491EndMacro
492
493macro buildDmesh(ThName)
494IFMACRO(!privateDmesh#ThName)
495NewMacro privateDmesh#ThName()privateDmesh#ThName EndMacro
496NewMacro privateDmesh#ThName#khi()privateDmesh#ThName#khiDef EndMacro
497NewMacro privateDmesh#ThName#intersection()privateDmesh#ThName#intersectionDef EndMacro
498real[int][int] privateDmesh#ThName#khi(2);
499real[int][int] privateDmesh#ThName#intersection;
500ENDIFMACRO
501{
502IFMACRO(!ThName#Comm)
503NewMacro ThName#Comm()mpiCommWorld EndMacro
504ENDIFMACRO
505NewMacro privateBuildDmesh()1 EndMacro
506int[int][int] intersection;
507NewMacro privateDmesh#ThTab()privateDmesh#ThName EndMacro
508NewMacro privateDmesh#ThTab#khi()privateDmesh#ThName#khiDef EndMacro
509NewMacro privateDmesh#ThTab#intersection()privateDmesh#ThName#intersectionDef EndMacro
510IFMACRO(ThName#N2O)
511NewMacro privateDmesh#N2O()ThName#N2O EndMacro
512ENDIFMACRO
513build(ThName, 1, intersection, privateDmesh#ThName#khi[0], P1, ThName#Comm);
514}
515EndMacro
516macro reconstructDmesh(ThName)
517IFMACRO(!privateDmesh#ThName)
518NewMacro privateDmesh#ThName()privateDmesh#ThName EndMacro
519NewMacro privateDmesh#ThName#khi()privateDmesh#ThName#khiDef EndMacro
520NewMacro privateDmesh#ThName#intersection()privateDmesh#ThName#intersectionDef EndMacro
521real[int][int] privateDmesh#ThName#khi(2);
522real[int][int] privateDmesh#ThName#intersection;
523ENDIFMACRO
524IFMACRO(!ThName#Comm)
525NewMacro ThName#Comm()mpiCommWorld EndMacro
526ENDIFMACRO
527if(ThName#Comm) {
528    int[int] neighbors;
529    {
530        real[int] bb(2 * dimension);
531        boundingbox(ThName, bb);
532        real[int] bbAll(2 * dimension * mpiSize(ThName#Comm));
533        mpiAllgather(bb, bbAll, ThName#Comm);
534        real hmax;
535        {
536            real tmp = ThName.hmax;
537            mpiAllReduce(tmp, hmax, ThName#Comm, mpiMAX);
538        }
539        int between = 0;
540        for(int i = 0; i < mpiSize(ThName#Comm); ++i) {
541            if(i != mpiRank(ThName#Comm) &&
542IFMACRO(dimension, 2)
543            !(bbAll[1 + 4 * i] < bb[0] - hmax
544              || bbAll[0 + 4 * i] > bb[1] + hmax
545              || bbAll[3 + 4 * i] < bb[2] - hmax
546              || bbAll[2 + 4 * i] > bb[3] + hmax)
547ENDIFMACRO
548IFMACRO(dimension, 3)
549            !(bbAll[1 + 6 * i] < bb[0] - hmax
550              || bbAll[0 + 6 * i] > bb[1] + hmax
551              || bbAll[3 + 6 * i] < bb[2] - hmax
552              || bbAll[2 + 6 * i] > bb[3] + hmax
553              || bbAll[5 + 6 * i] < bb[4] - hmax
554              || bbAll[4 + 6 * i] > bb[5] + hmax)
555ENDIFMACRO
556                                                 ) {
557                neighbors.resize(neighbors.n + 1);
558                neighbors[neighbors.n - 1] = i;
559            }
560        }
561    }
562    reconstructDmeshWithNeighbors(ThName, neighbors)
563}
564EndMacro
565macro reconstructDmeshWithNeighbors(ThName, neighborsName)
566IFMACRO(!privateDmesh#ThName)
567NewMacro privateDmesh#ThName()privateDmesh#ThName EndMacro
568NewMacro privateDmesh#ThName#khi()privateDmesh#ThName#khiDef EndMacro
569NewMacro privateDmesh#ThName#intersection()privateDmesh#ThName#intersectionDef EndMacro
570real[int][int] privateDmesh#ThName#khi(2);
571real[int][int] privateDmesh#ThName#intersection;
572ENDIFMACRO
573{
574    real[int] part;
575    {
576        if(verbosity > 0)
577            mpiBarrier(ThName#Comm);
578        real timerReconstruction = mpiWtime();
579        varf vG(u, v) = on(labels(ThName), u = 1.0);
580        fespace VhGamma(ThName, P1);
581        fespace PhGamma(ThName, P0);
582        VhGamma gamma;
583        gamma[] = vG(0, VhGamma, tgv = -1);
584        PhGamma gammaElt = gamma > 0.1;
585        meshN ThLocalInit = trunc(ThName, gammaElt > 0.1, label = -111112);
586        meshN ThLocalInitInterior = trunc(ThName, gammaElt < 0.1, label = -111112);
587        neighborsName.sort;
588        int between = 0;
589        for(int i = 0; i < neighborsName.n; ++i)
590            if(neighborsName[i] > mpiRank(ThName#Comm)) {
591                between = i;
592                break;
593            }
594        if(neighborsName.n)
595            if(neighborsName[neighborsName.n - 1] < mpiRank(ThName#Comm))
596                between = neighborsName.n;
597        meshN[int] ThTab(neighborsName.n + 1);
598        ThTab[between] = ThLocalInit;
599        mpiRequest[int] rqRecv(neighborsName.n);
600        mpiRequest[int] rqSend(neighborsName.n);
601        for(int i = 0; i < neighborsName.n; ++i)
602            Isend(processor(neighborsName[i], ThName#Comm, rqSend[i]), ThLocalInit);
603        for(int i = 0; i < between; ++i)
604            Irecv(processor(neighborsName[i], ThName#Comm, rqRecv[i]), ThTab[i]);
605        for(int i = between; i < neighborsName.n; ++i)
606            Irecv(processor(neighborsName[i], ThName#Comm, rqRecv[i]), ThTab[i + 1]);
607        mpiWaitAll(rqRecv);
608        meshN ThLocalNew = gluemesh(ThTab);
609        int m = 0;
610        for(int i = 0; i < between; ++i)
611            m += ThTab[i].nt;
612        ThTab[between] = trunc(ThLocalNew, nuTriangle >= m && nuTriangle < m + ThTab[between].nt, label = -111111);
613        mpiWaitAll(rqSend);
614        for(int i = 0; i < neighborsName.n; ++i)
615            Isend(processor(neighborsName[i], ThName#Comm, rqSend[i]), ThTab[between]);
616        for(int i = 0; i < between; ++i)
617            Irecv(processor(neighborsName[i], ThName#Comm, rqRecv[i]), ThTab[i]);
618        for(int i = between; i < neighborsName.n; ++i)
619            Irecv(processor(neighborsName[i], ThName#Comm, rqRecv[i]), ThTab[i + 1]);
620        mpiWaitAll(rqRecv);
621        ThTab.resize(neighborsName.n + 2);
622        ThTab[neighborsName.n + 1] = ThLocalInitInterior;
623        ThName = gluemesh(ThTab);
624IFMACRO(dimension, 3)
625        ThName = change(ThName, rmlfaces = -111112);
626ENDIFMACRO
627IFMACRO(dimension, 2)
628        ThName = change(ThName, rmledges = -111112);
629ENDIFMACRO
630        part.resize(ThName.nt);
631        m = 0;
632        for(int i = 0; i < between; ++i) {
633            part(m:m + ThTab[i].nt - 1) = neighborsName[i];
634            m += ThTab[i].nt;
635        }
636        part(m:m + ThTab[between].nt - 1) = mpiRank(ThName#Comm);
637        m += ThTab[between].nt;
638        for(int i = between; i < neighborsName.n; ++i) {
639            part(m:m + ThTab[i + 1].nt - 1) = neighborsName[i];
640            m += ThTab[i + 1].nt;
641        }
642        part(part.n - ThLocalInitInterior.nt:ThName.nt - 1) = mpiRank(ThName#Comm);
643        mpiWaitAll(rqSend);
644        if(verbosity > 0) {
645            mpiBarrier(ThName#Comm);
646            if(mpiRank(ThName#Comm) == 0)
647                cout.scientific << " --- distributed mesh reconstructed (in " << mpiWtime() - timerReconstruction << ")" << endl;
648        }
649    }
650    NewMacro privateBuildDmesh()1 EndMacro
651    NewMacro privateReconstructDmesh()1 EndMacro
652    int[int][int] intersection;
653    NewMacro privateDmesh#ThTab()privateDmesh#ThName EndMacro
654    NewMacro privateDmesh#ThTab#khi()privateDmesh#ThName#khiDef EndMacro
655    NewMacro privateDmesh#ThTab#intersection()privateDmesh#ThName#intersectionDef EndMacro
656IFMACRO(ThName#N2O)
657    NewMacro privateDmesh#N2O()ThName#N2O EndMacro
658ENDIFMACRO
659    buildWithPartitioning(ThName, part, 1, intersection, privateDmesh#ThName#khi[0], P1, ThName#Comm)
660}
661EndMacro
662macro copyDmesh(OldName, NewName)
663IFMACRO(!privateDmesh#NewName)
664NewMacro privateDmesh#NewName()privateDmesh#NewName EndMacro
665NewMacro privateDmesh#NewName#khi()privateDmesh#NewName#khiDef EndMacro
666NewMacro privateDmesh#NewName#intersection()privateDmesh#NewName#intersectionDef EndMacro
667real[int][int] privateDmesh#NewName#khi(2);
668real[int][int] privateDmesh#NewName#intersection;
669ENDIFMACRO
670IFMACRO(privateDmesh#OldName)
671NewName = OldName;
672privateDmesh#NewName#khi[0].resize(privateDmesh#OldName#khi[0].n);
673privateDmesh#NewName#khi[0] = privateDmesh#OldName#khi[0];
674privateDmesh#NewName#khi[1].resize(privateDmesh#OldName#khi[1].n);
675privateDmesh#NewName#khi[1] = privateDmesh#OldName#khi[1];
676privateDmesh#NewName#intersection.resize(privateDmesh#OldName#intersection.n);
677for(int i = 0; i < privateDmesh#NewName#intersection.n; ++i) {
678    privateDmesh#NewName#intersection[i].resize(privateDmesh#OldName#intersection[i].n);
679    privateDmesh#NewName#intersection[i] = privateDmesh#OldName#intersection[i];
680}
681ENDIFMACRO
682EndMacro
683macro createMat(ThName, MatName, PkName)
684IFMACRO(privateDmesh#ThName)
685{
686IFMACRO(!ThName#Comm)
687NewMacro ThName#Comm()mpiCommWorld EndMacro
688ENDIFMACRO
689IFMACRO(!privateCreateMatCheckDmesh)
690if(ThName.nv != privateDmesh#ThName#khi[0].n || (privateDmesh#ThName#khi[1].n && ThName.nt != privateDmesh#ThName#khi[1].n)) {
691    buildDmesh(ThName)
692}
693ENDIFMACRO
694NewMacro privateCreateMat()1 EndMacro
695int[int][int] intersection;
696real[int][int] DTab(1);
697meshN[int] ThTab(1);
698ThTab[0] = ThName;
699NewMacro privateDmesh#ThTab()privateDmesh#ThName EndMacro
700NewMacro privateDmesh#ThTab#khi()privateDmesh#ThName#khiDef EndMacro
701NewMacro privateDmesh#ThTab#intersection()privateDmesh#ThName#intersectionDef EndMacro
702IFMACRO(!def)
703NewMacro def(i)i EndMacro
704ENDIFMACRO
705IFMACRO(!init)
706NewMacro init(i)i EndMacro
707ENDIFMACRO
708if(mpiSize(ThName#Comm) > 1) {
709    partition(ThTab, privateCreateMat, privateCreateMat, privateCreateMat, privateCreateMat, privateCreateMat, mpiRank(ThName#Comm), mpiSize(ThName#Comm), 1, 1, 1, privateCreateMat, DTab, PkName, intersection, ThName#Comm, -111111, PkName, def, init, 1)
710}
711else {
712    fespace WhGlobal(ThName, PkName);
713    DTab[0].resize(WhGlobal.ndof);
714    DTab[0] = 1;
715    intersection.resize(0);
716}
717IFMACRO(!privateCreatePartition)
718constructor(MatName, DTab[0].n, intersection, DTab[0], communicator = ThName#Comm);
719ENDIFMACRO
720IFMACRO(privateCreatePartition)
721privateCreatePartition.resize(DTab[0].n);
722privateCreatePartition = DTab[0];
723ENDIFMACRO
724}
725ENDIFMACRO
726IFMACRO(!privateDmesh#ThName)
727buildDmesh(ThName)
728{
729    NewMacro privateCreateMatCheckDmesh()1 EndMacro
730    createMat(ThName, MatName, PkName)
731}
732ENDIFMACRO
733EndMacro
734
735macro createPartition(ThName, PartName, PkName)
736IFMACRO(!privateDmesh#ThName)
737buildDmesh(ThName)
738ENDIFMACRO
739{
740    NewMacro privateCreateMatCheckDmesh()1 EndMacro
741    NewMacro privateCreatePartition()PartName EndMacro
742    createMat(ThName, privateCreatePartition, PkName)
743}
744EndMacro
745
746macro buildOverlapEdgePeriodicRecursive(Th, ThBorder, fakeInterface, s, overlap, level, prolongation, intersection, DTab, P, comm, excluded, PkPart, defPart, initPart, labPeriodic, userPartitioning, bs) {
747IFMACRO(!def)
748    NewMacro def(i)i EndMacro
749ENDIFMACRO
750IFMACRO(!init)
751    NewMacro init(i)i EndMacro
752ENDIFMACRO
753    Th.resize(level);
754    ThBorder.resize(level);
755    prolongation.resize(level - 1);
756    real timerPartition = mpiWtime();
757    if(mpiSize(comm) > 1 && !excluded) {
758        meshN ThGlobal = Th[level - 1];
759        fespace PhGlobal(ThGlobal, P0);
760        fespace VhGlobal(ThGlobal, P1);
761        PhGlobal partGlobal;
762IFMACRO(!privateReconstructDmesh)
763        if(userPartitioning.n != PhGlobal.ndof || labPeriodic.n > 0) {
764            timerPartition = mpiWtime();
765            meshN ThGlobalPeriodic;
766            if(labPeriodic.n > 0) {
767                VhGlobal marker;
768                for(int i = 0; i < labPeriodic.n; ++i) {
769                    varf vMarker(u, v) = on(labPeriodic[i], u = 1.0);
770                    marker[] += vMarker(0, VhGlobal, tgv = -1);
771                }
772                PhGlobal partPeriodic = marker > 0.1;
773                while(1) {
774                    AddLayers(ThGlobal, partPeriodic[], 1 + overlap, marker[]);
775                    partPeriodic = marker > 0.001;
776                    ThGlobalPeriodic = trunc(ThGlobal, partPeriodic < 0.999);
777                    if(ThGlobal.nt / real(ThGlobalPeriodic.nt) > mpisize / real(mpisize - 1))
778                        break;
779                }
780            }
781            if(mpiRank(comm) == 0) {
782                if(verbosity > 0)
783                    cout.scientific << " --- global mesh of " << ThGlobal.nt << " elements (prior to refinement) partitioned with " << Stringification(partitioner);
784                if(labPeriodic.n > 0) {
785                    fespace PhPeriodic(ThGlobalPeriodic, P0);
786                    PhPeriodic partPeriodic;
787                    if(mpiSize(comm) > 2) {
788                        partitionerSeq(partPeriodic[], ThGlobalPeriodic, mpiSize(comm) - 1);
789                        partPeriodic[] += 1.0;
790                    }
791                    else
792                        partPeriodic[] = 1.0;
793                    partGlobal = partPeriodic;
794                }
795                else {
796                    partitionerSeq(partGlobal[], ThGlobal, mpiSize(comm));
797                }
798            }
799            if(labPeriodic.n > 0 && Stringification(partitioner) != "metis" && Stringification(partitioner) != "scotch") {
800                fespace PhPeriodic(ThGlobalPeriodic, P0);
801                PhPeriodic partPeriodic;
802                if(mpiSize(comm) > 2) {
803                    partitionerPar(partPeriodic[], ThGlobalPeriodic, comm, mpiSize(comm) - 1);
804                    partPeriodic[] += 1.0;
805                }
806                else
807                    partPeriodic[] = 1.0;
808                partGlobal = partPeriodic;
809            }
810            else
811                partitionerPar(partGlobal[], ThGlobal, comm, mpiSize(comm));
812            if(mpiRank(comm) == 0 && verbosity > 0)
813                cout.scientific << " (in " << mpiWtime() - timerPartition << ")" << endl;
814            timerPartition = mpiWtime();
815        }
816        else {
817            partGlobal[] = userPartitioning;
818        }
819ENDIFMACRO
820IFMACRO(privateReconstructDmesh)
821        partGlobal[] = userPartitioning;
822ENDIFMACRO
823IFMACRO(!trueRestrict)
824        bool trueRestrict = usedARGV("-true_restrict") != -1;
825ENDIFMACRO
826IFMACRO(!removeZeros)
827        bool removeZeros = trueRestrict && overlap == 1 && usedARGV("-remove_zeros") != -1;
828ENDIFMACRO
829        if(verbosity > 0) {
830            mpiBarrier(comm);
831            timerPartition = mpiWtime();
832        }
833IFMACRO(privateBuildDmesh)
834        NewMacro defP1(i)i EndMacro
835        NewMacro initP1(i)i EndMacro
836        partition(Th, ThBorder, ThGlobal, PhGlobal, VhGlobal, partGlobal, mpiRank(comm), mpiSize(comm), s, overlap, level, prolongation, DTab, P, intersection, comm, fakeInterface, PkPart, defP1, initP1, bs)
837ENDIFMACRO
838IFMACRO(!privateBuildDmesh)
839        partition(Th, ThBorder, ThGlobal, PhGlobal, VhGlobal, partGlobal, mpiRank(comm), mpiSize(comm), s, overlap, level, prolongation, DTab, P, intersection, comm, fakeInterface, PkPart, defPart, initPart, bs)
840ENDIFMACRO
841    }
842    else if(mpiSize(comm) == 1) {
843        for(int i = level - 1; i > 0; --i) {
844            Th[i - 1] = trunc(Th[i], 1, split = s);
845            fespace WhLocalRefined(Th[i - 1], P);
846            fespace WhLocalCoarse(Th[i], P);
847            prolongation[i - 1] = interpolate(WhLocalRefined, WhLocalCoarse);
848            DTab[i].resize(WhLocalCoarse.ndof);
849            DTab[i] = 1.0;
850        }
851        if(level == 1) {
852IFMACRO(privateBuildDmesh)
853IFMACRO(privateDmesh#N2O)
854            if(s > 1)
855                Th[0] = trunc(Th[0], 1, split = s, new2old = privateDmesh#N2O);
856            else {
857                privateDmesh#N2O.resize(Th[0].nt);
858                privateDmesh#N2O = 0:Th[0].nt-1;
859            }
860ENDIFMACRO
861IFMACRO(!privateDmesh#N2O)
862            if(s > 1)
863                Th[0] = trunc(Th[0], 1, split = s);
864ENDIFMACRO
865ENDIFMACRO
866IFMACRO(!privateBuildDmesh)
867            if(s > 1)
868                Th[0] = trunc(Th[0], 1, split = s);
869ENDIFMACRO
870        }
871        fespace WhLocal(Th[0], P);
872        DTab[0].resize(WhLocal.ndof);
873        DTab[0] = 1.0;
874    }
875    if(verbosity > 0) {
876        mpiBarrier(comm);
877        if(mpiRank(comm) == 0)
878            cout.scientific << " --- partition of unity built (in " << mpiWtime() - timerPartition << ")" << endl;
879    }
880}// EOM
881
882macro buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart, labPeriodic, userPartitioning, bs) {
883    meshN[int] ThTab(1);
884    meshN[int] ThBorderTab(1);
885    real[int][int] DTab(1);
886    ThTab[0] = Th;
887    matrix[int] prolongation(0);
888    buildOverlapEdgePeriodicRecursive(ThTab, ThBorderTab, fakeInterface, s, overlap, 1, prolongation, intersection, DTab, P, comm, excluded, PkPart, defPart, initPart, labPeriodic, userPartitioning, bs)
889    Th = ThTab[0];
890    ThBorder = ThBorderTab[0];
891    D.resize(DTab[0].n);
892    D = DTab[0];
893}// EOM
894
895IFMACRO(vectorialfe)
896macro buildOverlapEdgeRecursive(Th, ThBorder, fakeInterface, s, overlap, level, prolongation, intersection, D, P, comm, excluded, PkPart, defPart, initPart, bs) {
897    int[int] emptyArray(0);
898    real[int] emptyRealArray(0);
899    buildOverlapEdgePeriodicRecursive(Th, ThBorder, fakeInterface, s, overlap, level, prolongation, intersection, D, P, comm, excluded, PkPart, defPart, initPart, emptyArray, emptyRealArray, bs)
900}// EOM
901macro buildOverlapEdge(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart, bs) {
902    int[int] emptyArray(0);
903    real[int] emptyRealArray(0);
904    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart, emptyArray, emptyRealArray, bs)
905}// EOM
906macro buildOverlapEdgeWithPartitioning(Th, ThBorder, part, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart, bs) {
907    int[int] emptyArray(0);
908    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart, emptyArray, part, bs)
909}// EOM
910macro buildOverlapWithPartitioning(Th, ThBorder, part, fakeInterface, s, overlap, intersection, D, P, comm, excluded, bs) {
911    int[int] emptyArray(0);
912    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, emptyArray, part, bs)
913}// EOM
914macro buildOverlap(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, bs) {
915    int[int] emptyArray(0);
916    real[int] emptyRealArray(0);
917    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, emptyArray, emptyRealArray, bs)
918}// EOM
919macro buildOverlapPeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, labPeriodic, bs) {
920    real[int] emptyArray(0);
921    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, labPeriodic, emptyArray, bs)
922}// EOM
923macro buildEdgeWithPartitioning(Th, part, s, intersection, D, P, comm, PkPart, defPart, initPart, bs) {
924    int[int] emptyArray(0);
925    meshN ThBorder;
926    int fakeInterface = -111111;
927    int overlap = 1;
928    bool excluded = false;
929    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, intersection, D, P, comm, excluded, PkPart, defPart, initPart, emptyArray, part, bs)
930}// EOM
931macro buildWithPartitioning(Th, part, s, intersection, D, P, comm, bs) {
932    int[int] emptyArray(0);
933    meshN ThBorder;
934    int fakeInterface = -111111;
935    int overlap = 1;
936    bool excluded = false;
937    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, emptyArray, part, bs)
938}// EOM
939macro build(Th, s, intersection, D, P, comm, bs) {
940    int[int] emptyArray(0);
941    real[int] emptyRealArray(0);
942    meshN ThBorder;
943    int fakeInterface = -111111;
944    int overlap = 1;
945    bool excluded = false;
946    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, emptyArray, emptyRealArray, bs)
947}// EOM
948macro buildPeriodic(Th, s, intersection, D, P, comm, labPeriodic, bs) {
949    int[int] emptyArray(0);
950    real[int] emptyRealArray(0);
951    meshN ThBorder;
952    int fakeInterface = -111111;
953    int overlap = 1;
954    bool excluded = false;
955    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, labPeriodic, emptyRealArray, bs)
956}// EOM
957macro buildMinimalist(Th, intersection, D, P, bs) {
958    int[int] emptyArray(0);
959    real[int] emptyRealArray(0);
960    meshN ThBorder;
961    int fakeInterface = -111111;
962    int overlap = 1;
963    bool excluded = false;
964    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, 1, overlap, intersection, D, P, mpiCommWorld, excluded, P, def, init, emptyArray, emptyRealArray, bs)
965}// EOM
966macro buildRecursive(Th, s, level, prolongation, intersectionMat, DTab, P, comm, bsMat) {
967    int[int] emptyArray(0);
968    real[int] emptyRealArray(0);
969    meshN[int] ThBorderTab(level);
970    DTab.resize(level);
971    buildOverlapEdgePeriodicRecursive(Th, ThBorderTab, -111111, s, 1, level, prolongation, intersectionMat, DTab, P, comm, false, P, def, init, emptyArray, emptyRealArray, bsMat)
972}// EOM
973macro buildMatRecursive(Th, s, level, prolongation, A, P, comm, bsMat) {
974    int[int] emptyArray(0);
975    real[int] emptyRealArray(0);
976    meshN[int] ThBorderTab(level);
977    int[int][int] intersectionMat;
978    real[int][int] DTab(level);
979    buildOverlapEdgePeriodicRecursive(Th, ThBorderTab, -111111, s, 1, level, prolongation, intersectionMat, DTab, P, comm, false, P, def, init, emptyArray, emptyRealArray, bsMat)
980    for(int i = 0; i < level; ++i)
981        constructor(A[i], DTab[i].n, intersectionMat, DTab[i], bs = bsMat, communicator = comm, level = i);
982}// EOM
983macro buildMatEdgeRecursive(Th, s, level, prolongation, A, P, comm, PkPart, defPart, initPart, bsMat) {
984    int[int] emptyArray(0);
985    real[int] emptyRealArray(0);
986    meshN[int] ThBorderTab(level);
987    int[int][int] intersectionMat;
988    real[int][int] DTab(level);
989    buildOverlapEdgePeriodicRecursive(Th, ThBorderTab, -111111, s, 1, level, prolongation, intersectionMat, DTab, P, comm, false, PkPart, defPart, initPart, emptyArray, emptyRealArray, bsMat)
990    for(int i = 0; i < level; ++i)
991        constructor(A[i], DTab[i].n, intersectionMat, DTab[i], bs = bsMat, communicator = comm, level = i);
992}// EOM
993macro buildMatEdgeWithPartitioning(Th, part, s, A, P, comm, PkPart, defPart, initPart, bsMat) {
994    real[int] DMat;
995    int[int][int] intersectionMat;
996    buildEdgeWithPartitioning(Th, part, s, intersectionMat, DMat, P, comm, PkPart, defPart, initPart, bsMat)
997    constructor(A, DMat.n, intersectionMat, DMat, bs = bsMat, communicator = comm);
998}// EOM
999macro buildMatWithPartitioning(Th, part, s, A, P, comm, bsMat) {
1000    real[int] DMat;
1001    int[int][int] intersectionMat;
1002    buildWithPartitioning(Th, part, s, intersectionMat, DMat, P, comm, bsMat)
1003    constructor(A, DMat.n, intersectionMat, DMat, bs = bsMat, communicator = comm);
1004}// EOM
1005macro buildMat(Th, s, A, P, comm, bsMat) {
1006    real[int] DMat;
1007    int[int][int] intersectionMat;
1008    build(Th, s, intersectionMat, DMat, P, comm, bsMat)
1009    constructor(A, DMat.n, intersectionMat, DMat, bs = bsMat, communicator = comm);
1010}// EOM
1011macro buildMatPeriodic(Th, s, A, P, comm, labPeriodic, bsMat) {
1012    real[int] DMat;
1013    int[int][int] intersectionMat;
1014    buildPeriodic(Th, s, intersectionMat, DMat, P, comm, labPeriodic, bsMat)
1015    constructor(A, DMat.n, intersectionMat, DMat, bs = bsMat, communicator = comm);
1016}// EOM
1017macro buildMatMinimalist(Th, A, P, bsMat) {
1018    real[int] DMat;
1019    int[int][int] intersectionMat;
1020    buildMinimalist(Th, intersectionMat, DMat, P, bsMat)
1021    constructor(A, DMat.n, intersectionMat, DMat, bs = bsMat, communicator = comm);
1022}// EOM
1023ENDIFMACRO
1024IFMACRO(!vectorialfe)
1025macro buildOverlapEdgeRecursive(Th, ThBorder, fakeInterface, s, overlap, level, prolongation, intersection, D, P, comm, excluded, PkPart, defPart, initPart) {
1026    int[int] emptyArray(0);
1027    real[int] emptyRealArray(0);
1028    buildOverlapEdgePeriodicRecursive(Th, ThBorder, fakeInterface, s, overlap, level, prolongation, intersection, D, P, comm, excluded, PkPart, defPart, initPart, emptyArray, emptyRealArray, 1)
1029}// EOM
1030macro buildOverlapEdge(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart) {
1031    int[int] emptyArray(0);
1032    real[int] emptyRealArray(0);
1033    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart, emptyArray, emptyRealArray, 1)
1034}// EOM
1035macro buildOverlapEdgeWithPartitioning(Th, ThBorder, part, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart) {
1036    int[int] emptyArray(0);
1037    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart, emptyArray, part, 1)
1038}// EOM
1039macro buildOverlapWithPartitioning(Th, ThBorder, part, fakeInterface, s, overlap, intersection, D, P, comm, excluded) {
1040    int[int] emptyArray(0);
1041    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, emptyArray, part, 1)
1042}// EOM
1043macro buildOverlap(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded) {
1044    int[int] emptyArray(0);
1045    real[int] emptyRealArray(0);
1046    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, emptyArray, emptyRealArray, 1)
1047}// EOM
1048macro buildOverlapPeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, labPeriodic) {
1049    real[int] emptyArray(0);
1050    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, labPeriodic, emptyArray, 1)
1051}// EOM
1052macro buildEdgeWithPartitioning(Th, part, s, intersection, D, P, comm, PkPart, defPart, initPart) {
1053    int[int] emptyArray(0);
1054    meshN ThBorder;
1055    int fakeInterface = -111111;
1056    int overlap = 1;
1057    bool excluded = false;
1058    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, PkPart, defPart, initPart, emptyArray, part, 1)
1059}// EOM
1060macro buildWithPartitioning(Th, part, s, intersection, D, P, comm) {
1061    int[int] emptyArray(0);
1062    meshN ThBorder;
1063    int fakeInterface = -111111;
1064    int overlap = 1;
1065    bool excluded = false;
1066    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, emptyArray, part, 1)
1067}// EOM
1068macro build(Th, s, intersection, D, P, comm) {
1069    int[int] emptyArray(0);
1070    real[int] emptyRealArray(0);
1071    meshN ThBorder;
1072    int fakeInterface = -111111;
1073    int overlap = 1;
1074    bool excluded = false;
1075    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, emptyArray, emptyRealArray, 1)
1076}// EOM
1077macro buildPeriodic(Th, s, intersection, D, P, comm, labPeriodic) {
1078    int[int] emptyArray(0);
1079    real[int] emptyRealArray(0);
1080    meshN ThBorder;
1081    int fakeInterface = -111111;
1082    int overlap = 1;
1083    bool excluded = false;
1084    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, s, overlap, intersection, D, P, comm, excluded, P, def, init, labPeriodic, emptyRealArray, 1)
1085}// EOM
1086macro buildMinimalist(Th, intersection, D, P) {
1087    int[int] emptyArray(0);
1088    real[int] emptyRealArray(0);
1089    meshN ThBorder;
1090    int fakeInterface = -111111;
1091    int overlap = 1;
1092    bool excluded = false;
1093    buildOverlapEdgePeriodic(Th, ThBorder, fakeInterface, 1, overlap, intersection, D, P, mpiCommWorld, excluded, P, def, init, emptyArray, emptyRealArray, 1)
1094}// EOM
1095macro buildRecursive(Th, s, level, prolongation, intersectionMat, DTab, P, comm) {
1096    int[int] emptyArray(0);
1097    real[int] emptyRealArray(0);
1098    meshN[int] ThBorderTab(level);
1099    DTab.resize(level);
1100    buildOverlapEdgePeriodicRecursive(Th, ThBorderTab, -111111, s, 1, level, prolongation, intersectionMat, DTab, P, comm, false, P, def, init, emptyArray, emptyRealArray, 1)
1101}// EOM
1102macro buildMatRecursive(Th, s, level, prolongation, A, P, comm) {
1103    int[int] emptyArray(0);
1104    real[int] emptyRealArray(0);
1105    meshN[int] ThBorderTab(level);
1106    int[int][int] intersectionMat;
1107    real[int][int] DTab(level);
1108    buildOverlapEdgePeriodicRecursive(Th, ThBorderTab, -111111, s, 1, level, prolongation, intersectionMat, DTab, P, comm, false, P, def, init, emptyArray, emptyRealArray, 1)
1109    for(int i = 0; i < level; ++i)
1110        constructor(A[i], DTab[i].n, intersectionMat, DTab[i], communicator = comm, level = i);
1111}// EOM
1112macro buildMatEdgeRecursive(Th, s, level, prolongation, A, P, comm, PkPart, defPart, initPart) {
1113    int[int] emptyArray(0);
1114    real[int] emptyRealArray(0);
1115    meshN[int] ThBorderTab(level);
1116    int[int][int] intersectionMat;
1117    real[int][int] DTab(level);
1118    buildOverlapEdgePeriodicRecursive(Th, ThBorderTab, -111111, s, 1, level, prolongation, intersectionMat, DTab, P, comm, false, PkPart, defPart, initPart, emptyArray, emptyRealArray, 1)
1119    for(int i = 0; i < level; ++i)
1120        constructor(A[i], DTab[i].n, intersectionMat, DTab[i], communicator = comm, level = i);
1121}// EOM
1122macro buildMatEdgeWithPartitioning(Th, part, s, A, P, comm, PkPart, defPart, initPart) {
1123    real[int] DMat;
1124    int[int][int] intersectionMat;
1125    buildEdgeWithPartitioning(Th, part, s, intersectionMat, DMat, P, comm, PkPart, defPart, initPart)
1126    constructor(A, DMat.n, intersectionMat, DMat, communicator = comm);
1127}// EOM
1128macro buildMatWithPartitioning(Th, part, s, A, P, comm) {
1129    real[int] DMat;
1130    int[int][int] intersectionMat;
1131    buildWithPartitioning(Th, part, s, intersectionMat, DMat, P, comm)
1132    constructor(A, DMat.n, intersectionMat, DMat, communicator = comm);
1133}// EOM
1134macro buildMat(Th, s, A, P, comm) {
1135    real[int] DMat;
1136    int[int][int] intersectionMat;
1137    build(Th, s, intersectionMat, DMat, P, comm)
1138    constructor(A, DMat.n, intersectionMat, DMat, communicator = comm);
1139}// EOM
1140macro buildMatPeriodic(Th, s, A, P, comm, labPeriodic) {
1141    real[int] DMat;
1142    int[int][int] intersectionMat;
1143    buildPeriodic(Th, s, intersectionMat, DMat, P, comm, labPeriodic)
1144    constructor(A, DMat.n, intersectionMat, DMat, communicator = comm);
1145}// EOM
1146macro buildMatMinimalist(Th, A, P) {
1147    real[int] DMat;
1148    int[int][int] intersectionMat;
1149    buildMinimalist(Th, intersectionMat, DMat, P)
1150    constructor(A, DMat.n, intersectionMat, DMat);
1151}// EOM
1152ENDIFMACRO
1153
1154macro transfer(ThName, Pk, u, ThNew, PkNew, uNew)
1155IFMACRO(privateDmesh#ThName)
1156{
1157IFMACRO(!ThName#Comm)
1158NewMacro ThName#Comm()mpiCommWorld EndMacro
1159ENDIFMACRO
1160if(verbosity > 0)
1161    mpiBarrier(ThName#Comm);
1162real timerTransfer = mpiWtime();
1163IFMACRO(!def)
1164NewMacro def(i)i EndMacro
1165ENDIFMACRO
1166if(mpiSize(ThName#Comm) == 1) {
1167    def(uNew) = def(u);
1168}
1169else {
1170    int backupSM = searchMethod;
1171    searchMethod = 0;
1172    fespace VhLocalOld(ThName, Pk);
1173    assert(u[].n == VhLocalOld.ndof);
1174    fespace VhLocalNew(ThNew, PkNew);
1175    assert(uNew[].n == VhLocalNew.ndof);
1176    real[int] bb(4 * dimension);
1177    {
1178        real[int] tmp(2 * dimension);
1179        boundingbox(ThName, tmp);
1180        bb(0:2 * dimension - 1) = tmp;
1181        boundingbox(ThNew, tmp);
1182        bb(2 * dimension:4 * dimension - 1) = tmp;
1183    }
1184    real[int] bbAll(4 * dimension * mpiSize(ThName#Comm));
1185    mpiAllgather(bb, bbAll, ThName#Comm);
1186    int[int] rankSend(0);
1187    int[int] rankRecv(0);
1188    for(int i = 0; i < mpiSize(ThName#Comm); ++i) {
1189IFMACRO(dimension, 2)
1190        if(!(bbAll[1 + 8 * i] < bb[4]
1191          || bbAll[0 + 8 * i] > bb[5]
1192          || bbAll[3 + 8 * i] < bb[6]
1193          || bbAll[2 + 8 * i] > bb[7]))
1194ENDIFMACRO
1195IFMACRO(dimension, 3)
1196        if(!(bbAll[1 + 12 * i] < bb[6]
1197          || bbAll[0 + 12 * i] > bb[7]
1198          || bbAll[3 + 12 * i] < bb[8]
1199          || bbAll[2 + 12 * i] > bb[9]
1200          || bbAll[5 + 12 * i] < bb[10]
1201          || bbAll[4 + 12 * i] > bb[11]))
1202ENDIFMACRO
1203                                                    {
1204            rankRecv.resize(rankRecv.n + 1);
1205            rankRecv[rankRecv.n - 1] = i;
1206        }
1207IFMACRO(dimension, 2)
1208        if(!(bbAll[5 + 8 * i] < bb[0]
1209          || bbAll[4 + 8 * i] > bb[1]
1210          || bbAll[7 + 8 * i] < bb[2]
1211          || bbAll[6 + 8 * i] > bb[3]))
1212ENDIFMACRO
1213IFMACRO(dimension, 3)
1214        if(!(bbAll[7 + 12 * i] < bb[0]
1215          || bbAll[6 + 12 * i] > bb[1]
1216          || bbAll[9 + 12 * i] < bb[2]
1217          || bbAll[8 + 12 * i] > bb[3]
1218          || bbAll[11 + 12 * i] < bb[4]
1219          || bbAll[10 + 12 * i] > bb[5]))
1220ENDIFMACRO
1221                                                    {
1222            rankSend.resize(rankSend.n + 1);
1223            rankSend[rankSend.n - 1] = i;
1224        }
1225    }
1226    real[int] D, backupRegion(ThName.nt);
1227    createPartition(ThName, D, Pk)
1228    VhLocalOld def(scaledU);
1229    scaledU[] = u[];
1230    scaledU[] .*= D;
1231    fespace PhPart(ThName, P0);
1232    {
1233        PhPart backup = region;
1234        backupRegion = backup[];
1235        int[int] newRegion(ThName.nt);
1236        for[i, v : privateDmesh#ThName#khiDef[1]] newRegion[i] = abs(v - mpiRank(ThName#Comm)) < 1.0e-2 ? 1 : 0;
1237        ThName = change(ThName, fregion = newRegion[nuTriangle]);
1238    }
1239    meshN[int] recvTh(rankRecv.n);
1240    meshN[int] sendTh(rankSend.n);
1241    real[int][int] exchangeU(rankSend.n + rankRecv.n);
1242    mpiRequest[int] rqSendTh(rankSend.n);
1243    mpiRequest[int] rqSendU(rankSend.n);
1244    mpiRequest[int] rqRecvTh(rankRecv.n);
1245    mpiRequest[int] rqRecvU(rankRecv.n);
1246    for[i, v : rankRecv]
1247        Irecv(processor(v, rqRecvTh[i]), recvTh[i]);
1248    for[i, v : rankSend] {
1249        PhPart part;
1250        part[] = privateDmesh#ThName#khiDef[1];
1251IFMACRO(dimension, 2)
1252        part = (bbAll[4 + 8 * v] < x
1253             && bbAll[5 + 8 * v] > x
1254             && bbAll[6 + 8 * v] < y
1255             && bbAll[7 + 8 * v] > y) ? 1.0 : 0.0;
1256ENDIFMACRO
1257IFMACRO(dimension, 3)
1258        part = (bbAll[6 + 12 * v] < x
1259             && bbAll[7 + 12 * v] > x
1260             && bbAll[8 + 12 * v] < y
1261             && bbAll[9 + 12 * v] > y
1262             && bbAll[10 + 12 * v] < z
1263             && bbAll[11 + 12 * v] > z) ? 1.0 : 0.0;
1264ENDIFMACRO
1265        if(part[].linfty > 1.0e-2) {
1266            int[int] n2o;
1267            sendTh[i] = trunc(ThName, part > 1.0e-2, new2old = n2o);
1268            fespace VhRestriction(sendTh[i], Pk);
1269            int[int] map;
1270            map = restrict(VhRestriction, VhLocalOld, n2o);
1271            exchangeU[rankRecv.n + i].resize(VhRestriction.ndof);
1272            for[j, w : map] exchangeU[rankRecv.n + i][j] = scaledU[][w];
1273            Isend(processor(v, rqSendTh[i]), sendTh[i]);
1274            Isend(processor(v, rqSendU[i]), exchangeU[rankRecv.n + i]);
1275        }
1276        else
1277            Isend(processor(v, rqSendTh[i]), sendTh[i]);
1278    }
1279    meshN gluedExchange;
1280    {
1281        meshN[int] toGlue(rankRecv.n);
1282        int j = 0;
1283        for[i, v : rankRecv] {
1284            int index = mpiWaitAny(rqRecvTh);
1285            if(recvTh[index].nt) {
1286                fespace VhRestriction(recvTh[index], Pk);
1287                exchangeU[index].resize(VhRestriction.ndof);
1288                Irecv(processor(rankRecv[index], rqRecvU[index]), exchangeU[index]);
1289                fespace PhRestriction(recvTh[index], P0);
1290                PhRestriction ind = region;
1291                if(abs(ind[].max - 1.0) < 1.0e-2) {
1292                    toGlue[j] = trunc(recvTh[index], ind > 1.0e-2);
1293                    ++j;
1294                }
1295            }
1296        }
1297        toGlue.resize(j);
1298        gluedExchange = gluemesh(toGlue);
1299    }
1300    meshN interpolateExchange;
1301    fespace VhExchange(gluedExchange, Pk);
1302    VhExchange def(uExchange);
1303    for[i, v : rankRecv] {
1304        int index = mpiWaitAny(rqRecvU);
1305        if(index != mpiUndefined) {
1306            if(recvTh[index].nt) {
1307                fespace VhRestriction(recvTh[index], Pk);
1308                matrix R = interpolate(VhRestriction, VhExchange);
1309                uExchange[] += R' * exchangeU[index];
1310            }
1311        }
1312    }
1313    searchMethod = backupSM;
1314    def(uNew) = def(uExchange);
1315    ThName = change(ThName, fregion = backupRegion[nuTriangle]);
1316    mpiWaitAll(rqSendTh);
1317    mpiWaitAll(rqSendU);
1318}
1319if(verbosity > 0) {
1320    mpiBarrier(ThName#Comm);
1321    if(mpiRank(ThName#Comm) == 0)
1322        cout.scientific << " --- distributed solution transferred (in " << mpiWtime() - timerTransfer << ")" << endl;
1323}
1324}
1325ENDIFMACRO
1326// EOM
1327
1328macro createParMmgCommunicators(ThName, ThParMmgName, ThN2O, ThCommunicators) {
1329IFMACRO(!privateDmesh#ThName)
1330assert(0);
1331ENDIFMACRO
1332    Mat A;
1333    createMat(ThName, A, P1);
1334    real[int] D(ThName.nt);
1335    createPartition(ThName, D, P0);
1336    fespace Ph(ThName, P0);
1337    Ph d;
1338    d[] = D;
1339    ThParMmgName = trunc(ThName, abs(d) > 1.0e-2, label = -111111, new2old = ThN2O);
1340    fespace VhWithoutOverlap(ThParMmgName, P1);
1341    varf vG(u, v) = on(-111111, u = 1.0);
1342    real[int] gamma(ThParMmgName.nv);
1343    gamma = vG(0, VhWithoutOverlap, tgv = -1);
1344    fespace VhWithOverlap(ThName, P1);
1345    int[int] rest = restrict(VhWithoutOverlap, VhWithOverlap, ThN2O);
1346    ParMmgCommunicators(A, gamma, rest, ThCommunicators);
1347}// EOM
1348
1349macro gatherDmesh(ThName, comm, ThGatherName) {
1350IFMACRO(!privateDmesh#ThName)
1351assert(0);
1352ENDIFMACRO
1353IFMACRO(!ThName#Comm)
1354NewMacro ThName#Comm()mpiCommWorld EndMacro
1355ENDIFMACRO
1356IFMACRO(!ThGatherName#Comm)
1357NewMacro ThGatherName#Comm()mpiCommWorld EndMacro
1358ENDIFMACRO
1359    if(verbosity > 0 && ThName#Comm)
1360        mpiBarrier(ThName#Comm);
1361    real timerGather = mpiWtime();
1362    int size;
1363    if(ThGatherName#Comm)
1364        size = mpiSize(comm);
1365    else
1366        size = 0;
1367    int reduce;
1368    mpiAllReduce(size, reduce, ThName#Comm, mpiSUM);
1369    assert(reduce == mpiSize(ThName#Comm));
1370    fespace PhPart(ThName, P0);
1371    meshN ThNoOverlap = trunc(ThName, abs(privateDmesh#ThName#khiDef[1][nuTriangle] - mpiRank(ThName#Comm)) < 1.0e-2, label = -111112);
1372    if(ThGatherName#Comm) {
1373        meshN[int] recvTh(size);
1374        mpiRequest[int] rqRecv(size - 1);
1375        for(int i = 1; i < size; ++i)
1376            Irecv(processor(i, comm, rqRecv[i - 1]), recvTh[i]);
1377        recvTh[0] = ThNoOverlap;
1378        mpiWaitAll(rqRecv);
1379        ThGatherName = gluemesh(recvTh);
1380    }
1381    else {
1382        mpiRequest rqSend;
1383        Isend(processor(0, comm, rqSend), ThNoOverlap);
1384        mpiWait(rqSend);
1385    }
1386    if(verbosity > 0 && ThName#Comm) {
1387        mpiBarrier(ThName#Comm);
1388        if(mpiRank(ThName#Comm) == 0)
1389            cout.scientific << " --- distributed mesh gathered (in " << mpiWtime() - timerGather << ")" << endl;
1390    }
1391}
1392reconstructDmesh(ThGatherName)
1393// EOM
1394
1395macro scatterDmesh(ThName, comm, ThScatterName) {
1396IFMACRO(!privateDmesh#ThName)
1397assert(0);
1398ENDIFMACRO
1399IFMACRO(!ThName#Comm)
1400NewMacro ThName#Comm()mpiCommWorld EndMacro
1401ENDIFMACRO
1402IFMACRO(!ThScatterName#Comm)
1403NewMacro ThScatterName#Comm()mpiCommWorld EndMacro
1404ENDIFMACRO
1405    if(verbosity > 0 && ThScatterName#Comm)
1406        mpiBarrier(ThScatterName#Comm);
1407    real timerScatter = mpiWtime();
1408    int size;
1409    if(ThName#Comm) {
1410        size = mpiSize(comm);
1411    }
1412    else
1413        size = 0;
1414    int reduce;
1415    mpiAllReduce(size, reduce, ThScatterName#Comm, mpiSUM);
1416    assert(reduce == mpiSize(ThScatterName#Comm));
1417    if(ThName#Comm) {
1418        meshN ThNoOverlap;
1419        if(mpiSize(ThName#Comm) == 1)
1420            ThNoOverlap = ThName;
1421        else
1422            ThNoOverlap = trunc(ThName, abs(privateDmesh#ThName#khiDef[1][nuTriangle] - mpiRank(ThName#Comm)) < 1.0e-2, label = -111112);
1423        fespace PhPart(ThNoOverlap, P0);
1424        PhPart part;
1425        partitionerSeq(part[], ThNoOverlap, mpiSize(comm));
1426        partitionerPar(part[], ThNoOverlap, mpiCommSelf, mpiSize(comm));
1427        meshN[int] sendTh(mpiSize(comm) - 1);
1428        mpiRequest[int] rqSend(mpiSize(comm) - 1);
1429        for(int i = 1; i < mpiSize(comm); ++i) {
1430            sendTh[i - 1] = trunc(ThNoOverlap, abs(part - i) < 1.0e-2, label = -111112);
1431            Isend(processor(i, comm, rqSend[i - 1]), sendTh[i - 1]);
1432        }
1433        ThScatterName = trunc(ThNoOverlap, abs(part) < 1.0e-2, label = -111112);
1434        mpiWaitAll(rqSend);
1435    }
1436    else if(ThScatterName#Comm) {
1437        mpiRequest rqRecv;
1438        Irecv(processor(0, comm, rqRecv), ThScatterName);
1439        mpiWait(rqRecv);
1440    }
1441    if(verbosity > 0 && ThScatterName#Comm) {
1442        mpiBarrier(ThScatterName#Comm);
1443        if(mpiRank(ThScatterName#Comm) == 0)
1444            cout.scientific << " --- distributed mesh scattered (in " << mpiWtime() - timerScatter << ")" << endl;
1445    }
1446}
1447reconstructDmesh(ThScatterName)
1448// EOM
1449
1450macro gatherSolution(ThName, comm, ThGatherName, Pk, u, uNew) {
1451IFMACRO(!privateDmesh#ThName)
1452assert(0);
1453ENDIFMACRO
1454IFMACRO(!privateDmesh#ThGatherName)
1455assert(0);
1456ENDIFMACRO
1457IFMACRO(!ThName#Comm)
1458NewMacro ThName#Comm()mpiCommWorld EndMacro
1459ENDIFMACRO
1460IFMACRO(!ThGatherName#Comm)
1461NewMacro ThGatherName#Comm()mpiCommWorld EndMacro
1462ENDIFMACRO
1463    if(verbosity > 0 && ThName#Comm)
1464        mpiBarrier(ThName#Comm);
1465    real timerGather = mpiWtime();
1466    if(ThGatherName#Comm) {
1467        meshN[int] recvTh(mpiSize(comm) - 1);
1468        real[int][int] recvU(mpiSize(comm) - 1);
1469        mpiRequest[int] rqRecvTh(mpiSize(comm) - 1);
1470        mpiRequest[int] rqRecvU(mpiSize(comm) - 1);
1471        for(int i = 0; i < mpiSize(comm) - 1; ++i)
1472            Irecv(processor(i + 1, comm, rqRecvTh[i]), recvTh[i]);
1473        for(int i = 0; i < mpiSize(comm) - 1; ++i) {
1474            int index = mpiWaitAny(rqRecvTh);
1475            fespace VhRecv(recvTh[index], Pk);
1476            recvU[index].resize(VhRecv.ndof);
1477            Irecv(processor(index + 1, comm, rqRecvU[index]), recvU[index]);
1478        }
1479        fespace VhGlobalGather(ThGatherName, Pk);
1480        real[int] visited(VhGlobalGather.ndof);
1481        visited = 1.0;
1482        {
1483            fespace VhRestriction(ThName, Pk);
1484            matrix R = interpolate(VhRestriction, VhGlobalGather);
1485            real[int] buffer = R' * u[];
1486            buffer .*= visited;
1487            real[int] ones(VhRestriction.ndof);
1488            ones = -1.0;
1489            visited += R' * ones;
1490            for[j, v : visited] v = max(v, 0.0);
1491            uNew[] += buffer;
1492        }
1493        for(int i = 0; i < mpiSize(comm) - 1; ++i) {
1494            int index = mpiWaitAny(rqRecvU);
1495            fespace VhRestriction(recvTh[index], Pk);
1496            matrix R = interpolate(VhRestriction, VhGlobalGather);
1497            real[int] buffer = R' * recvU[index];
1498            buffer .*= visited;
1499            real[int] ones(VhRestriction.ndof);
1500            ones = -1.0;
1501            visited += R' * ones;
1502            for[j, v : visited] v = max(v, 0.0);
1503            uNew[] += buffer;
1504        }
1505    }
1506    else {
1507        mpiRequest[int] rqSend(2);
1508        Isend(processor(0, comm, rqSend[0]), ThName);
1509        fespace VhLocalGather(ThName, Pk);
1510        assert(u[].n == VhLocalGather.ndof);
1511        Isend(processor(0, comm, rqSend[1]), u[]);
1512        mpiWaitAll(rqSend);
1513    }
1514    if(verbosity > 0 && ThName#Comm) {
1515        mpiBarrier(ThName#Comm);
1516        if(mpiRank(ThName#Comm) == 0)
1517            cout.scientific << " --- distributed solution gathered (in " << mpiWtime() - timerGather << ")" << endl;
1518    }
1519} // EOM
1520
1521macro scatterSolution(ThName, comm, ThScatterName, Pk, u, uNew) {
1522IFMACRO(!privateDmesh#ThName)
1523assert(0);
1524ENDIFMACRO
1525IFMACRO(!privateDmesh#ThScatterName)
1526assert(0);
1527ENDIFMACRO
1528IFMACRO(!def)
1529NewMacro def(i)i EndMacro
1530ENDIFMACRO
1531IFMACRO(!ThName#Comm)
1532NewMacro ThName#Comm()mpiCommWorld EndMacro
1533ENDIFMACRO
1534IFMACRO(!ThScatterName#Comm)
1535NewMacro ThScatterName#Comm()mpiCommWorld EndMacro
1536ENDIFMACRO
1537    if(verbosity > 0 && ThScatterName#Comm)
1538        mpiBarrier(ThScatterName#Comm);
1539    real timerScatter = mpiWtime();
1540    if(mpiRank(comm) == 0) {
1541        broadcast(processor(0, comm), ThName);
1542        broadcast(processor(0, comm), u[]);
1543        def(uNew) = def(u);
1544    }
1545    else {
1546        meshN ThGlobalScatter;
1547        broadcast(processor(0, comm), ThGlobalScatter);
1548        fespace VhGlobalScatter(ThGlobalScatter, Pk);
1549        VhGlobalScatter uGlobalScatter;
1550        broadcast(processor(0, comm), uGlobalScatter[]);
1551        def(uNew) = def(uGlobalScatter);
1552    }
1553    if(verbosity > 0 && ThScatterName#Comm) {
1554        mpiBarrier(ThScatterName#Comm);
1555        if(mpiRank(ThScatterName#Comm) == 0)
1556            cout.scientific << " --- distributed solution scattered (in " << mpiWtime() - timerScatter << ")" << endl;
1557    }
1558}
1559// EOM
1560