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