00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "regionSplit.H"
00027 #include <OpenFOAM/cyclicPolyPatch.H>
00028 #include <OpenFOAM/processorPolyPatch.H>
00029 #include <OpenFOAM/globalIndex.H>
00030 #include <OpenFOAM/syncTools.H>
00031
00032
00033
00034 namespace Foam
00035 {
00036
00037 defineTypeNameAndDebug(regionSplit, 0);
00038
00039 }
00040
00041
00042
00043
00044 void Foam::regionSplit::transferCoupledFaceRegion
00045 (
00046 const label faceI,
00047 const label otherFaceI,
00048
00049 labelList& faceRegion,
00050 DynamicList<label>& newChangedFaces
00051 ) const
00052 {
00053 if (faceRegion[faceI] >= 0)
00054 {
00055 if (faceRegion[otherFaceI] == -1)
00056 {
00057 faceRegion[otherFaceI] = faceRegion[faceI];
00058 newChangedFaces.append(otherFaceI);
00059 }
00060 else if (faceRegion[otherFaceI] == -2)
00061 {
00062
00063
00064 }
00065 else if (faceRegion[otherFaceI] != faceRegion[faceI])
00066 {
00067 FatalErrorIn
00068 (
00069 "regionSplit::transferCoupledFaceRegion"
00070 "(const label, const label, labelList&, labelList&) const"
00071 ) << "Problem : coupled face " << faceI
00072 << " on patch " << mesh_.boundaryMesh().whichPatch(faceI)
00073 << " has region " << faceRegion[faceI]
00074 << " but coupled face " << otherFaceI
00075 << " has region " << faceRegion[otherFaceI]
00076 << endl
00077 << "Is your blocked faces specification"
00078 << " synchronized across coupled boundaries?"
00079 << abort(FatalError);
00080 }
00081 }
00082 else if (faceRegion[faceI] == -1)
00083 {
00084 if (faceRegion[otherFaceI] >= 0)
00085 {
00086 faceRegion[faceI] = faceRegion[otherFaceI];
00087 newChangedFaces.append(faceI);
00088 }
00089 else if (faceRegion[otherFaceI] == -2)
00090 {
00091
00092
00093 }
00094 }
00095 }
00096
00097
00098 void Foam::regionSplit::fillSeedMask
00099 (
00100 const List<labelPair>& explicitConnections,
00101 labelList& cellRegion,
00102 labelList& faceRegion,
00103 const label seedCellID,
00104 const label markValue
00105 ) const
00106 {
00107
00108 cellRegion[seedCellID] = markValue;
00109
00110
00111
00112 const cell& cFaces = mesh_.cells()[seedCellID];
00113
00114 label nFaces = 0;
00115
00116 labelList changedFaces(cFaces.size());
00117
00118 forAll(cFaces, i)
00119 {
00120 label faceI = cFaces[i];
00121
00122 if (faceRegion[faceI] == -1)
00123 {
00124 faceRegion[faceI] = markValue;
00125 changedFaces[nFaces++] = faceI;
00126 }
00127 }
00128 changedFaces.setSize(nFaces);
00129
00130
00131
00132
00133 while (changedFaces.size())
00134 {
00135
00136
00137
00138
00139
00140
00141 DynamicList<label> changedCells(changedFaces.size());
00142
00143 forAll(changedFaces, i)
00144 {
00145 label faceI = changedFaces[i];
00146
00147 label own = mesh_.faceOwner()[faceI];
00148
00149 if (cellRegion[own] == -1)
00150 {
00151 cellRegion[own] = markValue;
00152 changedCells.append(own);
00153 }
00154
00155 if (mesh_.isInternalFace(faceI))
00156 {
00157 label nei = mesh_.faceNeighbour()[faceI];
00158
00159 if (cellRegion[nei] == -1)
00160 {
00161 cellRegion[nei] = markValue;
00162 changedCells.append(nei);
00163 }
00164 }
00165 }
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175 DynamicList<label> newChangedFaces(changedCells.size());
00176
00177 forAll(changedCells, i)
00178 {
00179 label cellI = changedCells[i];
00180
00181 const cell& cFaces = mesh_.cells()[cellI];
00182
00183 forAll(cFaces, cFaceI)
00184 {
00185 label faceI = cFaces[cFaceI];
00186
00187 if (faceRegion[faceI] == -1)
00188 {
00189 faceRegion[faceI] = markValue;
00190 newChangedFaces.append(faceI);
00191 }
00192 }
00193 }
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00207
00208 forAll(patches, patchI)
00209 {
00210 const polyPatch& pp = patches[patchI];
00211
00212 if (isA<cyclicPolyPatch>(pp))
00213 {
00214 label faceI = pp.start();
00215
00216 label halfSz = pp.size()/2;
00217
00218 for (label i = 0; i < halfSz; i++)
00219 {
00220 label otherFaceI = refCast<const cyclicPolyPatch>(pp)
00221 .transformGlobalFace(faceI);
00222
00223 transferCoupledFaceRegion
00224 (
00225 faceI,
00226 otherFaceI,
00227 faceRegion,
00228 newChangedFaces
00229 );
00230
00231 faceI++;
00232 }
00233 }
00234 }
00235 forAll(explicitConnections, i)
00236 {
00237 transferCoupledFaceRegion
00238 (
00239 explicitConnections[i][0],
00240 explicitConnections[i][1],
00241 faceRegion,
00242 newChangedFaces
00243 );
00244 }
00245
00246
00247
00248
00249
00250
00251
00252 changedFaces.transfer(newChangedFaces);
00253 }
00254 }
00255
00256
00257 Foam::label Foam::regionSplit::calcRegionSplit
00258 (
00259 const boolList& blockedFace,
00260 const List<labelPair>& explicitConnections,
00261
00262 labelList& cellRegion
00263 ) const
00264 {
00265 if (debug)
00266 {
00267 if (blockedFace.size())
00268 {
00269
00270 boolList syncBlockedFace(blockedFace);
00271 syncTools::swapFaceList(mesh_, syncBlockedFace, false);
00272
00273 forAll(syncBlockedFace, faceI)
00274 {
00275 if (syncBlockedFace[faceI] != blockedFace[faceI])
00276 {
00277 FatalErrorIn
00278 (
00279 "regionSplit::calcRegionSplit(..)"
00280 ) << "Face " << faceI << " not synchronised. My value:"
00281 << blockedFace[faceI] << " coupled value:"
00282 << syncBlockedFace[faceI]
00283 << abort(FatalError);
00284 }
00285 }
00286 }
00287 }
00288
00289
00290
00291
00292 labelList faceRegion(mesh_.nFaces(), -1);
00293
00294 if (blockedFace.size())
00295 {
00296 forAll(blockedFace, faceI)
00297 {
00298 if (blockedFace[faceI])
00299 {
00300 faceRegion[faceI] = -2;
00301 }
00302 }
00303 }
00304
00305
00306
00307
00308
00309
00310 label nRegions = 0;
00311
00312 label unsetCellI = 0;
00313
00314 do
00315 {
00316
00317
00318 for (; unsetCellI < mesh_.nCells(); unsetCellI++)
00319 {
00320 if (cellRegion[unsetCellI] == -1)
00321 {
00322 break;
00323 }
00324 }
00325
00326 if (unsetCellI >= mesh_.nCells())
00327 {
00328 break;
00329 }
00330
00331 fillSeedMask
00332 (
00333 explicitConnections,
00334 cellRegion,
00335 faceRegion,
00336 unsetCellI,
00337 nRegions
00338 );
00339
00340
00341 nRegions++;
00342 unsetCellI++;
00343 }
00344 while(true);
00345
00346
00347 if (debug)
00348 {
00349 forAll(cellRegion, cellI)
00350 {
00351 if (cellRegion[cellI] < 0)
00352 {
00353 FatalErrorIn("regionSplit::calcRegionSplit(..)")
00354 << "cell:" << cellI << " region:" << cellRegion[cellI]
00355 << abort(FatalError);
00356 }
00357 }
00358
00359 forAll(faceRegion, faceI)
00360 {
00361 if (faceRegion[faceI] == -1)
00362 {
00363 FatalErrorIn("regionSplit::calcRegionSplit(..)")
00364 << "face:" << faceI << " region:" << faceRegion[faceI]
00365 << abort(FatalError);
00366 }
00367 }
00368 }
00369
00370
00371
00372
00373
00374
00375
00376 globalIndex globalRegions(nRegions);
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387 labelList mergedGlobal(identity(globalRegions.size()));
00388
00389
00390
00391 while (Pstream::parRun())
00392 {
00393 if (debug)
00394 {
00395 Pout<< nl << "-- Starting Iteration --" << endl;
00396 }
00397
00398 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00399
00400
00401 forAll(patches, patchI)
00402 {
00403 const polyPatch& pp = patches[patchI];
00404
00405 if (isA<processorPolyPatch>(pp))
00406 {
00407 labelList myGlobalRegions(pp.size());
00408
00409 label faceI = pp.start();
00410
00411 forAll(pp, i)
00412 {
00413 if (faceRegion[faceI] < 0)
00414 {
00415 myGlobalRegions[i] = faceRegion[faceI];
00416 }
00417 else
00418 {
00419 myGlobalRegions[i] = mergedGlobal
00420 [globalRegions.toGlobal(faceRegion[faceI])];
00421 }
00422
00423 faceI++;
00424 }
00425
00426 OPstream toProcNbr
00427 (
00428 Pstream::blocking,
00429 refCast<const processorPolyPatch>(pp).neighbProcNo()
00430 );
00431
00432 toProcNbr << myGlobalRegions;
00433 }
00434 }
00435
00436
00437
00438
00439 label nMerged = 0;
00440
00441 forAll(patches, patchI)
00442 {
00443 const polyPatch& pp = patches[patchI];
00444
00445 if (isA<processorPolyPatch>(pp))
00446 {
00447 const processorPolyPatch& procPp =
00448 refCast<const processorPolyPatch>(pp);
00449
00450 IPstream fromProcNbr(Pstream::blocking, procPp.neighbProcNo());
00451
00452 labelList nbrRegions(fromProcNbr);
00453
00454
00455
00456
00457 label faceI = pp.start();
00458
00459 forAll(pp, i)
00460 {
00461 if
00462 (
00463 faceRegion[faceI] < 0
00464 || nbrRegions[i] < 0
00465 )
00466 {
00467 if (faceRegion[faceI] != nbrRegions[i])
00468 {
00469 FatalErrorIn("regionSplit::calcRegionSplit(..)")
00470 << "On patch:" << pp.name()
00471 << " face:" << faceI
00472 << " my local region:" << faceRegion[faceI]
00473 << " neighbouring region:"
00474 << nbrRegions[i] << nl
00475 << "Maybe your blockedFaces are not"
00476 << " synchronized across coupled faces?"
00477 << abort(FatalError);
00478 }
00479 }
00480 else
00481 {
00482 label uncompactGlobal =
00483 globalRegions.toGlobal(faceRegion[faceI]);
00484
00485 label myGlobal = mergedGlobal[uncompactGlobal];
00486
00487 if (myGlobal != nbrRegions[i])
00488 {
00489 label minRegion = min(myGlobal, nbrRegions[i]);
00490
00491 if (debug)
00492 {
00493 Pout<< "Merging region " << myGlobal
00494 << " (on proc " << Pstream::myProcNo()
00495 << ") and region " << nbrRegions[i]
00496 << " (on proc " << procPp.neighbProcNo()
00497 << ") into region " << minRegion << endl;
00498 }
00499
00500 mergedGlobal[uncompactGlobal] = minRegion;
00501 mergedGlobal[myGlobal] = minRegion;
00502 mergedGlobal[nbrRegions[i]] = minRegion;
00503
00504 nMerged++;
00505 }
00506 }
00507
00508 faceI++;
00509 }
00510 }
00511 }
00512
00513
00514 reduce(nMerged, sumOp<label>());
00515
00516 if (debug)
00517 {
00518 Pout<< "nMerged:" << nMerged << endl;
00519 }
00520
00521 if (nMerged == 0)
00522 {
00523 break;
00524 }
00525
00526
00527 Pstream::listCombineGather(mergedGlobal, minEqOp<label>());
00528 Pstream::listCombineScatter(mergedGlobal);
00529 }
00530
00531
00532
00533
00534
00535
00536
00537
00538 labelList mergedToCompacted(globalRegions.size(), -1);
00539
00540 label compactI = 0;
00541
00542 forAll(mergedGlobal, i)
00543 {
00544 label merged = mergedGlobal[i];
00545
00546 if (mergedToCompacted[merged] == -1)
00547 {
00548 mergedToCompacted[merged] = compactI++;
00549 }
00550 }
00551
00552 if (debug)
00553 {
00554 Pout<< "Compacted down to " << compactI << " regions." << endl;
00555 }
00556
00557
00558 forAll(cellRegion, cellI)
00559 {
00560 label region = cellRegion[cellI];
00561
00562 if (region >= 0)
00563 {
00564 label merged = mergedGlobal[globalRegions.toGlobal(region)];
00565
00566 cellRegion[cellI] = mergedToCompacted[merged];
00567 }
00568 }
00569
00570 return compactI;
00571 }
00572
00573
00574
00575
00576 Foam::regionSplit::regionSplit(const polyMesh& mesh)
00577 :
00578 labelList(mesh.nCells(), -1),
00579 mesh_(mesh),
00580 nRegions_(calcRegionSplit(boolList(0, false), List<labelPair>(0), *this))
00581 {}
00582
00583
00584 Foam::regionSplit::regionSplit
00585 (
00586 const polyMesh& mesh,
00587 const boolList& blockedFace
00588 )
00589 :
00590 labelList(mesh.nCells(), -1),
00591 mesh_(mesh),
00592 nRegions_(calcRegionSplit(blockedFace, List<labelPair>(0), *this))
00593 {}
00594
00595
00596 Foam::regionSplit::regionSplit
00597 (
00598 const polyMesh& mesh,
00599 const boolList& blockedFace,
00600 const List<labelPair>& explicitConnections
00601 )
00602 :
00603 labelList(mesh.nCells(), -1),
00604 mesh_(mesh),
00605 nRegions_(calcRegionSplit(blockedFace, explicitConnections, *this))
00606 {}
00607
00608
00609
00610
00611
00612