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