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 "meshRefinement.H"
00027 #include <finiteVolume/fvMesh.H>
00028 #include <OpenFOAM/syncTools.H>
00029 #include <OpenFOAM/Time.H>
00030 #include <autoMesh/refinementSurfaces.H>
00031 #include <meshTools/pointSet.H>
00032 #include <meshTools/faceSet.H>
00033 #include <OpenFOAM/indirectPrimitivePatch.H>
00034 #include <dynamicMesh/polyTopoChange.H>
00035 #include <meshTools/meshTools.H>
00036 #include <dynamicMesh/polyModifyFace.H>
00037 #include <dynamicMesh/polyModifyCell.H>
00038 #include <dynamicMesh/polyAddFace.H>
00039 #include <dynamicMesh/polyRemoveFace.H>
00040 #include <dynamicMesh/polyAddPoint.H>
00041 #include <dynamicMesh/localPointRegion.H>
00042 #include <dynamicMesh/duplicatePoints.H>
00043 #include <OpenFOAM/OFstream.H>
00044 #include <meshTools/regionSplit.H>
00045 #include <dynamicMesh/removeCells.H>
00046
00047
00048
00049
00050
00051
00052 Foam::label Foam::meshRefinement::createBaffle
00053 (
00054 const label faceI,
00055 const label ownPatch,
00056 const label neiPatch,
00057 polyTopoChange& meshMod
00058 ) const
00059 {
00060 const face& f = mesh_.faces()[faceI];
00061 label zoneID = mesh_.faceZones().whichZone(faceI);
00062 bool zoneFlip = false;
00063
00064 if (zoneID >= 0)
00065 {
00066 const faceZone& fZone = mesh_.faceZones()[zoneID];
00067 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
00068 }
00069
00070 meshMod.setAction
00071 (
00072 polyModifyFace
00073 (
00074 f,
00075 faceI,
00076 mesh_.faceOwner()[faceI],
00077 -1,
00078 false,
00079 ownPatch,
00080 false,
00081 zoneID,
00082 zoneFlip
00083 )
00084 );
00085
00086
00087 label dupFaceI = -1;
00088
00089 if (mesh_.isInternalFace(faceI))
00090 {
00091 if (neiPatch == -1)
00092 {
00093 FatalErrorIn
00094 (
00095 "meshRefinement::createBaffle"
00096 "(const label, const label, const label, polyTopoChange&)"
00097 ) << "No neighbour patch for internal face " << faceI
00098 << " fc:" << mesh_.faceCentres()[faceI]
00099 << " ownPatch:" << ownPatch << abort(FatalError);
00100 }
00101
00102 bool reverseFlip = false;
00103 if (zoneID >= 0)
00104 {
00105 reverseFlip = !zoneFlip;
00106 }
00107
00108 dupFaceI = meshMod.setAction
00109 (
00110 polyAddFace
00111 (
00112 f.reverseFace(),
00113 mesh_.faceNeighbour()[faceI],
00114 -1,
00115 -1,
00116 -1,
00117 faceI,
00118 true,
00119 neiPatch,
00120 zoneID,
00121 reverseFlip
00122 )
00123 );
00124 }
00125 return dupFaceI;
00126 }
00127
00128
00129
00130 Foam::label Foam::meshRefinement::getBafflePatch
00131 (
00132 const labelList& facePatch,
00133 const label faceI
00134 ) const
00135 {
00136 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00137
00138
00139
00140
00141
00142
00143
00144
00145 forAll(mesh_.faces()[faceI], fp)
00146 {
00147 label pointI = mesh_.faces()[faceI][fp];
00148
00149 forAll(mesh_.pointFaces()[pointI], pf)
00150 {
00151 label pFaceI = mesh_.pointFaces()[pointI][pf];
00152
00153 label patchI = patches.whichPatch(pFaceI);
00154
00155 if (patchI != -1 && !patches[patchI].coupled())
00156 {
00157 return patchI;
00158 }
00159 else if (facePatch[pFaceI] != -1)
00160 {
00161 return facePatch[pFaceI];
00162 }
00163 }
00164 }
00165
00166
00167
00168 const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
00169
00170 forAll(ownFaces, i)
00171 {
00172 label cFaceI = ownFaces[i];
00173
00174 label patchI = patches.whichPatch(cFaceI);
00175
00176 if (patchI != -1 && !patches[patchI].coupled())
00177 {
00178 return patchI;
00179 }
00180 else if (facePatch[cFaceI] != -1)
00181 {
00182 return facePatch[cFaceI];
00183 }
00184 }
00185
00186 if (mesh_.isInternalFace(faceI))
00187 {
00188 const cell& neiFaces = mesh_.cells()[mesh_.faceNeighbour()[faceI]];
00189
00190 forAll(neiFaces, i)
00191 {
00192 label cFaceI = neiFaces[i];
00193
00194 label patchI = patches.whichPatch(cFaceI);
00195
00196 if (patchI != -1 && !patches[patchI].coupled())
00197 {
00198 return patchI;
00199 }
00200 else if (facePatch[cFaceI] != -1)
00201 {
00202 return facePatch[cFaceI];
00203 }
00204 }
00205 }
00206
00207 WarningIn
00208 (
00209 "meshRefinement::getBafflePatch(const labelList&, const label)"
00210 ) << "Could not find boundary face neighbouring internal face "
00211 << faceI << " with face centre " << mesh_.faceCentres()[faceI]
00212 << nl
00213 << "Using arbitrary patch " << 0 << " instead." << endl;
00214
00215 return 0;
00216 }
00217
00218
00219
00220 void Foam::meshRefinement::getBafflePatches
00221 (
00222 const labelList& globalToPatch,
00223 const labelList& neiLevel,
00224 const pointField& neiCc,
00225
00226 labelList& ownPatch,
00227 labelList& neiPatch
00228 ) const
00229 {
00230 autoPtr<OFstream> str;
00231 label vertI = 0;
00232 if (debug&OBJINTERSECTIONS)
00233 {
00234 str.reset
00235 (
00236 new OFstream
00237 (
00238 mesh_.time().path()/timeName()/"intersections.obj"
00239 )
00240 );
00241
00242 Pout<< "getBafflePatches : Writing surface intersections to file "
00243 << str().name() << nl << endl;
00244 }
00245
00246 const pointField& cellCentres = mesh_.cellCentres();
00247
00248
00249 const wordList& cellZoneNames = surfaces_.cellZoneNames();
00250
00251 labelList surfacesToBaffle(cellZoneNames.size());
00252 label baffleI = 0;
00253 forAll(cellZoneNames, surfI)
00254 {
00255 if (cellZoneNames[surfI].size())
00256 {
00257 if (debug)
00258 {
00259 Pout<< "getBafflePatches : Not baffling surface "
00260 << surfaces_.names()[surfI] << endl;
00261 }
00262 }
00263 else
00264 {
00265 surfacesToBaffle[baffleI++] = surfI;
00266 }
00267 }
00268 surfacesToBaffle.setSize(baffleI);
00269
00270 ownPatch.setSize(mesh_.nFaces());
00271 ownPatch = -1;
00272 neiPatch.setSize(mesh_.nFaces());
00273 neiPatch = -1;
00274
00275
00276
00277
00278
00279 labelList testFaces(intersectedFaces());
00280
00281
00282
00283
00284 pointField start(testFaces.size());
00285 pointField end(testFaces.size());
00286
00287 forAll(testFaces, i)
00288 {
00289 label faceI = testFaces[i];
00290
00291 label own = mesh_.faceOwner()[faceI];
00292
00293 if (mesh_.isInternalFace(faceI))
00294 {
00295 start[i] = cellCentres[own];
00296 end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
00297 }
00298 else
00299 {
00300 start[i] = cellCentres[own];
00301 end[i] = neiCc[faceI-mesh_.nInternalFaces()];
00302 }
00303 }
00304
00305
00306
00307
00308 labelList surface1;
00309 List<pointIndexHit> hit1;
00310 labelList region1;
00311 labelList surface2;
00312 List<pointIndexHit> hit2;
00313 labelList region2;
00314 surfaces_.findNearestIntersection
00315 (
00316 surfacesToBaffle,
00317 start,
00318 end,
00319
00320 surface1,
00321 hit1,
00322 region1,
00323
00324 surface2,
00325 hit2,
00326 region2
00327 );
00328
00329 forAll(testFaces, i)
00330 {
00331 label faceI = testFaces[i];
00332
00333 if (hit1[i].hit() && hit2[i].hit())
00334 {
00335 if (str.valid())
00336 {
00337 meshTools::writeOBJ(str(), start[i]);
00338 vertI++;
00339 meshTools::writeOBJ(str(), hit1[i].rawPoint());
00340 vertI++;
00341 meshTools::writeOBJ(str(), hit2[i].rawPoint());
00342 vertI++;
00343 meshTools::writeOBJ(str(), end[i]);
00344 vertI++;
00345 str()<< "l " << vertI-3 << ' ' << vertI-2 << nl;
00346 str()<< "l " << vertI-2 << ' ' << vertI-1 << nl;
00347 str()<< "l " << vertI-1 << ' ' << vertI << nl;
00348 }
00349
00350
00351 ownPatch[faceI] = globalToPatch
00352 [
00353 surfaces_.globalRegion(surface1[i], region1[i])
00354 ];
00355 neiPatch[faceI] = globalToPatch
00356 [
00357 surfaces_.globalRegion(surface2[i], region2[i])
00358 ];
00359
00360 if (ownPatch[faceI] == -1 || neiPatch[faceI] == -1)
00361 {
00362 FatalErrorIn("getBafflePatches(..)")
00363 << "problem." << abort(FatalError);
00364 }
00365 }
00366 }
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
00377 syncTools::syncFaceList(mesh_, neiPatch, maxEqOp<label>(), false);
00378 }
00379
00380
00381
00382 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
00383 (
00384 const labelList& ownPatch,
00385 const labelList& neiPatch
00386 )
00387 {
00388 if
00389 (
00390 ownPatch.size() != mesh_.nFaces()
00391 || neiPatch.size() != mesh_.nFaces()
00392 )
00393 {
00394 FatalErrorIn
00395 (
00396 "meshRefinement::createBaffles"
00397 "(const labelList&, const labelList&)"
00398 ) << "Illegal size :"
00399 << " ownPatch:" << ownPatch.size()
00400 << " neiPatch:" << neiPatch.size()
00401 << ". Should be number of faces:" << mesh_.nFaces()
00402 << abort(FatalError);
00403 }
00404
00405 if (debug)
00406 {
00407 labelList syncedOwnPatch(ownPatch);
00408 syncTools::syncFaceList(mesh_, syncedOwnPatch, maxEqOp<label>(), false);
00409 labelList syncedNeiPatch(neiPatch);
00410 syncTools::syncFaceList(mesh_, syncedNeiPatch, maxEqOp<label>(), false);
00411
00412 forAll(syncedOwnPatch, faceI)
00413 {
00414 if
00415 (
00416 (ownPatch[faceI] == -1 && syncedOwnPatch[faceI] != -1)
00417 || (neiPatch[faceI] == -1 && syncedNeiPatch[faceI] != -1)
00418 )
00419 {
00420 FatalErrorIn
00421 (
00422 "meshRefinement::createBaffles"
00423 "(const labelList&, const labelList&)"
00424 ) << "Non synchronised at face:" << faceI
00425 << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
00426 << " fc:" << mesh_.faceCentres()[faceI] << endl
00427 << "ownPatch:" << ownPatch[faceI]
00428 << " syncedOwnPatch:" << syncedOwnPatch[faceI]
00429 << " neiPatch:" << neiPatch[faceI]
00430 << " syncedNeiPatch:" << syncedNeiPatch[faceI]
00431 << abort(FatalError);
00432 }
00433 }
00434 }
00435
00436
00437 polyTopoChange meshMod(mesh_);
00438
00439 label nBaffles = 0;
00440
00441 forAll(ownPatch, faceI)
00442 {
00443 if (ownPatch[faceI] != -1)
00444 {
00445
00446
00447 createBaffle
00448 (
00449 faceI,
00450 ownPatch[faceI],
00451 neiPatch[faceI],
00452 meshMod
00453 );
00454 nBaffles++;
00455 }
00456 }
00457
00458
00459 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
00460
00461
00462 mesh_.updateMesh(map);
00463
00464
00465 if (map().hasMotionPoints())
00466 {
00467 mesh_.movePoints(map().preMotionPoints());
00468 }
00469 else
00470 {
00471
00472 mesh_.clearOut();
00473 }
00474
00475 if (overwrite())
00476 {
00477 mesh_.setInstance(oldInstance());
00478 }
00479
00480
00481
00482 faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
00483
00484 const labelList& reverseFaceMap = map().reverseFaceMap();
00485 const labelList& faceMap = map().faceMap();
00486
00487
00488 forAll(ownPatch, oldFaceI)
00489 {
00490 label faceI = reverseFaceMap[oldFaceI];
00491
00492 if (ownPatch[oldFaceI] != -1 && faceI >= 0)
00493 {
00494 const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
00495
00496 forAll(ownFaces, i)
00497 {
00498 baffledFacesSet.insert(ownFaces[i]);
00499 }
00500 }
00501 }
00502
00503 forAll(faceMap, faceI)
00504 {
00505 label oldFaceI = faceMap[faceI];
00506
00507 if (oldFaceI >= 0 && reverseFaceMap[oldFaceI] != faceI)
00508 {
00509 const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
00510
00511 forAll(ownFaces, i)
00512 {
00513 baffledFacesSet.insert(ownFaces[i]);
00514 }
00515 }
00516 }
00517 baffledFacesSet.sync(mesh_);
00518
00519 updateMesh(map, baffledFacesSet.toc());
00520
00521 return map;
00522 }
00523
00524
00525
00526
00527
00528 Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
00529 (
00530 const labelList& testFaces
00531 ) const
00532 {
00533 labelList duplicateFace
00534 (
00535 localPointRegion::findDuplicateFaces
00536 (
00537 mesh_,
00538 testFaces
00539 )
00540 );
00541
00542 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00543
00544
00545 List<labelPair> duplicateFaces(testFaces.size());
00546 label dupI = 0;
00547
00548 forAll(duplicateFace, i)
00549 {
00550 label otherFaceI = duplicateFace[i];
00551
00552 if (otherFaceI != -1 && i < otherFaceI)
00553 {
00554 label meshFace0 = testFaces[i];
00555 label patch0 = patches.whichPatch(meshFace0);
00556 label meshFace1 = testFaces[otherFaceI];
00557 label patch1 = patches.whichPatch(meshFace1);
00558
00559 if
00560 (
00561 (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
00562 || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
00563 )
00564 {
00565 FatalErrorIn
00566 (
00567 "meshRefinement::getDuplicateFaces"
00568 "(const bool, const labelList&)"
00569 ) << "One of two duplicate faces is on"
00570 << " processorPolyPatch."
00571 << "This is not allowed." << nl
00572 << "Face:" << meshFace0
00573 << " is on patch:" << patches[patch0].name()
00574 << nl
00575 << "Face:" << meshFace1
00576 << " is on patch:" << patches[patch1].name()
00577 << abort(FatalError);
00578 }
00579
00580 duplicateFaces[dupI++] = labelPair(meshFace0, meshFace1);
00581 }
00582 }
00583 duplicateFaces.setSize(dupI);
00584
00585 Info<< "getDuplicateFaces : found " << returnReduce(dupI, sumOp<label>())
00586 << " pairs of duplicate faces." << nl << endl;
00587
00588
00589 if (debug)
00590 {
00591 faceSet duplicateFaceSet(mesh_, "duplicateFaces", 2*dupI);
00592
00593 forAll(duplicateFaces, i)
00594 {
00595 duplicateFaceSet.insert(duplicateFaces[i][0]);
00596 duplicateFaceSet.insert(duplicateFaces[i][1]);
00597 }
00598 Pout<< "Writing duplicate faces (baffles) to faceSet "
00599 << duplicateFaceSet.name() << nl << endl;
00600 duplicateFaceSet.write();
00601 }
00602
00603 return duplicateFaces;
00604 }
00605
00606
00607
00608
00609
00610
00611
00612 Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
00613 (
00614 const List<labelPair>& couples
00615 ) const
00616 {
00617
00618
00619
00620
00621
00622 labelList duplicateFaces(couples.size());
00623
00624
00625
00626
00627
00628 labelList nBafflesPerEdge(mesh_.nEdges(), 0);
00629
00630
00631
00632
00633
00634
00635 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00636
00637 forAll(patches, patchI)
00638 {
00639 const polyPatch& pp = patches[patchI];
00640
00641
00642 if (!pp.coupled())
00643 {
00644 label faceI = pp.start();
00645
00646 forAll(pp, i)
00647 {
00648 const labelList& fEdges = mesh_.faceEdges(faceI);
00649
00650 forAll(fEdges, fEdgeI)
00651 {
00652 nBafflesPerEdge[fEdges[fEdgeI]]++;
00653 }
00654 faceI++;
00655 }
00656 }
00657 }
00658
00659
00660 DynamicList<label> fe0;
00661 DynamicList<label> fe1;
00662
00663
00664
00665
00666
00667 forAll(couples, i)
00668 {
00669 const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0);
00670
00671 forAll(fEdges0, fEdgeI)
00672 {
00673 nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000;
00674 }
00675
00676 const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1);
00677
00678 forAll(fEdges1, fEdgeI)
00679 {
00680 nBafflesPerEdge[fEdges1[fEdgeI]] += 1000000;
00681 }
00682 }
00683
00684
00685 syncTools::syncEdgeList
00686 (
00687 mesh_,
00688 nBafflesPerEdge,
00689 plusEqOp<label>(),
00690 0,
00691 false
00692 );
00693
00694
00695
00696
00697
00698 List<labelPair> filteredCouples(couples.size());
00699 label filterI = 0;
00700
00701 forAll(couples, i)
00702 {
00703 const labelPair& couple = couples[i];
00704
00705 if
00706 (
00707 patches.whichPatch(couple.first())
00708 == patches.whichPatch(couple.second())
00709 )
00710 {
00711 const labelList& fEdges = mesh_.faceEdges(couples[i].first());
00712
00713 forAll(fEdges, fEdgeI)
00714 {
00715 label edgeI = fEdges[fEdgeI];
00716
00717 if (nBafflesPerEdge[edgeI] == 2*1000000+2*1)
00718 {
00719 filteredCouples[filterI++] = couples[i];
00720 break;
00721 }
00722 }
00723 }
00724 }
00725 filteredCouples.setSize(filterI);
00726
00727
00728
00729
00730
00731
00732
00733 return filteredCouples;
00734 }
00735
00736
00737
00738 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
00739 (
00740 const List<labelPair>& couples
00741 )
00742 {
00743
00744 polyTopoChange meshMod(mesh_);
00745
00746 const faceList& faces = mesh_.faces();
00747 const labelList& faceOwner = mesh_.faceOwner();
00748 const faceZoneMesh& faceZones = mesh_.faceZones();
00749
00750 forAll(couples, i)
00751 {
00752 label face0 = couples[i].first();
00753 label face1 = couples[i].second();
00754
00755
00756
00757 label own0 = faceOwner[face0];
00758 label own1 = faceOwner[face1];
00759
00760 if (face1 < 0 || own0 < own1)
00761 {
00762
00763 label zoneID = faceZones.whichZone(face0);
00764 bool zoneFlip = false;
00765
00766 if (zoneID >= 0)
00767 {
00768 const faceZone& fZone = faceZones[zoneID];
00769 zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
00770 }
00771
00772 label nei = (face1 < 0 ? -1 : own1);
00773
00774 meshMod.setAction(polyRemoveFace(face1));
00775 meshMod.setAction
00776 (
00777 polyModifyFace
00778 (
00779 faces[face0],
00780 face0,
00781 own0,
00782 nei,
00783 false,
00784 -1,
00785 false,
00786 zoneID,
00787 zoneFlip
00788 )
00789 );
00790 }
00791 else
00792 {
00793
00794 label zoneID = faceZones.whichZone(face1);
00795 bool zoneFlip = false;
00796
00797 if (zoneID >= 0)
00798 {
00799 const faceZone& fZone = faceZones[zoneID];
00800 zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
00801 }
00802
00803 meshMod.setAction(polyRemoveFace(face0));
00804 meshMod.setAction
00805 (
00806 polyModifyFace
00807 (
00808 faces[face1],
00809 face1,
00810 own1,
00811 own0,
00812 false,
00813 -1,
00814 false,
00815 zoneID,
00816 zoneFlip
00817 )
00818 );
00819 }
00820 }
00821
00822
00823 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
00824
00825
00826 mesh_.updateMesh(map);
00827
00828
00829 if (map().hasMotionPoints())
00830 {
00831 mesh_.movePoints(map().preMotionPoints());
00832 }
00833 else
00834 {
00835
00836 mesh_.clearOut();
00837 }
00838
00839 if (overwrite())
00840 {
00841 mesh_.setInstance(oldInstance());
00842 }
00843
00844
00845
00846
00847 labelList newExposedFaces(2*couples.size());
00848 label newI = 0;
00849
00850 forAll(couples, i)
00851 {
00852 label newFace0 = map().reverseFaceMap()[couples[i].first()];
00853 if (newFace0 != -1)
00854 {
00855 newExposedFaces[newI++] = newFace0;
00856 }
00857
00858 label newFace1 = map().reverseFaceMap()[couples[i].second()];
00859 if (newFace1 != -1)
00860 {
00861 newExposedFaces[newI++] = newFace1;
00862 }
00863 }
00864 newExposedFaces.setSize(newI);
00865 updateMesh(map, newExposedFaces);
00866
00867 return map;
00868 }
00869
00870
00871
00872 void Foam::meshRefinement::findCellZoneGeometric
00873 (
00874 const labelList& closedNamedSurfaces,
00875 labelList& namedSurfaceIndex,
00876 const labelList& surfaceToCellZone,
00877
00878 labelList& cellToZone
00879 ) const
00880 {
00881 const pointField& cellCentres = mesh_.cellCentres();
00882 const labelList& faceOwner = mesh_.faceOwner();
00883 const labelList& faceNeighbour = mesh_.faceNeighbour();
00884
00885
00886 labelList insideSurfaces;
00887 surfaces_.findInside
00888 (
00889 closedNamedSurfaces,
00890 cellCentres,
00891 insideSurfaces
00892 );
00893
00894 forAll(insideSurfaces, cellI)
00895 {
00896 if (cellToZone[cellI] == -2)
00897 {
00898 label surfI = insideSurfaces[cellI];
00899
00900 if (surfI != -1)
00901 {
00902 cellToZone[cellI] = surfaceToCellZone[surfI];
00903 }
00904 }
00905 }
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915 label nCandidates = 0;
00916 forAll(namedSurfaceIndex, faceI)
00917 {
00918 label surfI = namedSurfaceIndex[faceI];
00919
00920 if (surfI != -1)
00921 {
00922 if (mesh_.isInternalFace(faceI))
00923 {
00924 nCandidates += 2;
00925 }
00926 else
00927 {
00928 nCandidates += 1;
00929 }
00930 }
00931 }
00932
00933
00934 pointField candidatePoints(nCandidates);
00935 nCandidates = 0;
00936 forAll(namedSurfaceIndex, faceI)
00937 {
00938 label surfI = namedSurfaceIndex[faceI];
00939
00940 if (surfI != -1)
00941 {
00942 label own = faceOwner[faceI];
00943 const point& ownCc = cellCentres[own];
00944
00945 if (mesh_.isInternalFace(faceI))
00946 {
00947 label nei = faceNeighbour[faceI];
00948 const point& neiCc = cellCentres[nei];
00949
00950 const vector d = 1E-4*(neiCc - ownCc);
00951 candidatePoints[nCandidates++] = ownCc-d;
00952 candidatePoints[nCandidates++] = neiCc+d;
00953 }
00954 else
00955 {
00956 const point& neiFc = mesh_.faceCentres()[faceI];
00957
00958 const vector d = 1E-4*(neiFc - ownCc);
00959 candidatePoints[nCandidates++] = ownCc-d;
00960 }
00961 }
00962 }
00963
00964
00965
00966
00967 surfaces_.findInside
00968 (
00969 closedNamedSurfaces,
00970 candidatePoints,
00971 insideSurfaces
00972 );
00973
00974
00975
00976
00977 nCandidates = 0;
00978 forAll(namedSurfaceIndex, faceI)
00979 {
00980 label surfI = namedSurfaceIndex[faceI];
00981
00982 if (surfI != -1)
00983 {
00984 label own = faceOwner[faceI];
00985
00986 if (mesh_.isInternalFace(faceI))
00987 {
00988 label ownSurfI = insideSurfaces[nCandidates++];
00989 if (ownSurfI != -1)
00990 {
00991 cellToZone[own] = surfaceToCellZone[ownSurfI];
00992 }
00993
00994 label neiSurfI = insideSurfaces[nCandidates++];
00995 if (neiSurfI != -1)
00996 {
00997 label nei = faceNeighbour[faceI];
00998
00999 cellToZone[nei] = surfaceToCellZone[neiSurfI];
01000 }
01001 }
01002 else
01003 {
01004 label ownSurfI = insideSurfaces[nCandidates++];
01005 if (ownSurfI != -1)
01006 {
01007 cellToZone[own] = surfaceToCellZone[ownSurfI];
01008 }
01009 }
01010 }
01011 }
01012
01013
01014
01015
01016
01017
01018 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
01019 {
01020 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
01021 label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
01022
01023 if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
01024 {
01025
01026 namedSurfaceIndex[faceI] = findIndex
01027 (
01028 surfaceToCellZone,
01029 max(ownZone, neiZone)
01030 );
01031 }
01032 }
01033
01034 labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
01035 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
01036
01037 forAll(patches, patchI)
01038 {
01039 const polyPatch& pp = patches[patchI];
01040
01041 if (pp.coupled())
01042 {
01043 forAll(pp, i)
01044 {
01045 label faceI = pp.start()+i;
01046 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
01047 neiCellZone[faceI-mesh_.nInternalFaces()] = ownZone;
01048 }
01049 }
01050 }
01051 syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
01052
01053 forAll(patches, patchI)
01054 {
01055 const polyPatch& pp = patches[patchI];
01056
01057 if (pp.coupled())
01058 {
01059 forAll(pp, i)
01060 {
01061 label faceI = pp.start()+i;
01062 label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
01063 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
01064
01065 if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
01066 {
01067
01068 namedSurfaceIndex[faceI] = findIndex
01069 (
01070 surfaceToCellZone,
01071 max(ownZone, neiZone)
01072 );
01073 }
01074 }
01075 }
01076 }
01077
01078
01079 syncTools::syncFaceList
01080 (
01081 mesh_,
01082 namedSurfaceIndex,
01083 maxEqOp<label>(),
01084 false
01085 );
01086 }
01087
01088
01089 bool Foam::meshRefinement::calcRegionToZone
01090 (
01091 const label surfZoneI,
01092 const label ownRegion,
01093 const label neiRegion,
01094
01095 labelList& regionToCellZone
01096 ) const
01097 {
01098 bool changed = false;
01099
01100
01101 if (ownRegion != neiRegion)
01102 {
01103
01104
01105
01106
01107
01108 if (regionToCellZone[ownRegion] == -2)
01109 {
01110 if (regionToCellZone[neiRegion] == surfZoneI)
01111 {
01112
01113
01114 regionToCellZone[ownRegion] = -1;
01115 changed = true;
01116 }
01117 else if (regionToCellZone[neiRegion] != -2)
01118 {
01119
01120
01121 regionToCellZone[ownRegion] = surfZoneI;
01122 changed = true;
01123 }
01124 }
01125 else if (regionToCellZone[neiRegion] == -2)
01126 {
01127 if (regionToCellZone[ownRegion] == surfZoneI)
01128 {
01129
01130
01131 regionToCellZone[neiRegion] = -1;
01132 changed = true;
01133 }
01134 else if (regionToCellZone[ownRegion] != -2)
01135 {
01136
01137
01138 regionToCellZone[neiRegion] = surfZoneI;
01139 changed = true;
01140 }
01141 }
01142 }
01143 return changed;
01144 }
01145
01146
01147
01148
01149
01150
01151 void Foam::meshRefinement::findCellZoneTopo
01152 (
01153 const point& keepPoint,
01154 const labelList& namedSurfaceIndex,
01155 const labelList& surfaceToCellZone,
01156 labelList& cellToZone
01157 ) const
01158 {
01159
01160 boolList blockedFace(mesh_.nFaces());
01161
01162 forAll(namedSurfaceIndex, faceI)
01163 {
01164 if (namedSurfaceIndex[faceI] == -1)
01165 {
01166 blockedFace[faceI] = false;
01167 }
01168 else
01169 {
01170 blockedFace[faceI] = true;
01171 }
01172 }
01173
01174
01175
01176 regionSplit cellRegion(mesh_, blockedFace);
01177 blockedFace.clear();
01178
01179
01180
01181
01182
01183 labelList regionToCellZone(cellRegion.nRegions(), -2);
01184
01185
01186
01187
01188 forAll(cellToZone, cellI)
01189 {
01190 if (cellToZone[cellI] != -2)
01191 {
01192 regionToCellZone[cellRegion[cellI]] = cellToZone[cellI];
01193 }
01194 }
01195
01196
01197
01198
01199 label keepRegionI = -1;
01200
01201 label cellI = mesh_.findCell(keepPoint);
01202
01203 if (cellI != -1)
01204 {
01205 keepRegionI = cellRegion[cellI];
01206 }
01207 reduce(keepRegionI, maxOp<label>());
01208
01209 Info<< "Found point " << keepPoint << " in cell " << cellI
01210 << " in global region " << keepRegionI
01211 << " out of " << cellRegion.nRegions() << " regions." << endl;
01212
01213 if (keepRegionI == -1)
01214 {
01215 FatalErrorIn
01216 (
01217 "meshRefinement::findCellZoneTopo"
01218 "(const point&, const labelList&, const labelList&, labelList&)"
01219 ) << "Point " << keepPoint
01220 << " is not inside the mesh." << nl
01221 << "Bounding box of the mesh:" << mesh_.globalData().bb()
01222 << exit(FatalError);
01223 }
01224
01225
01226 if (regionToCellZone[keepRegionI] == -2)
01227 {
01228 regionToCellZone[keepRegionI] = -1;
01229 }
01230
01231
01232
01233
01234 while (true)
01235 {
01236
01237
01238
01239
01240
01241
01242
01243 Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
01244 Pstream::listCombineScatter(regionToCellZone);
01245
01246
01247 bool changed = false;
01248
01249
01250
01251 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
01252 {
01253 label surfI = namedSurfaceIndex[faceI];
01254
01255 if (surfI != -1)
01256 {
01257
01258
01259 bool changedCell = calcRegionToZone
01260 (
01261 surfaceToCellZone[surfI],
01262 cellRegion[mesh_.faceOwner()[faceI]],
01263 cellRegion[mesh_.faceNeighbour()[faceI]],
01264 regionToCellZone
01265 );
01266
01267 changed = changed | changedCell;
01268 }
01269 }
01270
01271
01272
01273 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
01274
01275
01276 labelList neiCellRegion(mesh_.nFaces()-mesh_.nInternalFaces());
01277 forAll(patches, patchI)
01278 {
01279 const polyPatch& pp = patches[patchI];
01280
01281 if (pp.coupled())
01282 {
01283 forAll(pp, i)
01284 {
01285 label faceI = pp.start()+i;
01286 neiCellRegion[faceI-mesh_.nInternalFaces()] =
01287 cellRegion[mesh_.faceOwner()[faceI]];
01288 }
01289 }
01290 }
01291 syncTools::swapBoundaryFaceList(mesh_, neiCellRegion, false);
01292
01293
01294
01295 forAll(patches, patchI)
01296 {
01297 const polyPatch& pp = patches[patchI];
01298
01299 if (pp.coupled())
01300 {
01301 forAll(pp, i)
01302 {
01303 label faceI = pp.start()+i;
01304
01305 label surfI = namedSurfaceIndex[faceI];
01306
01307 if (surfI != -1)
01308 {
01309 bool changedCell = calcRegionToZone
01310 (
01311 surfaceToCellZone[surfI],
01312 cellRegion[mesh_.faceOwner()[faceI]],
01313 neiCellRegion[faceI-mesh_.nInternalFaces()],
01314 regionToCellZone
01315 );
01316
01317 changed = changed | changedCell;
01318 }
01319 }
01320 }
01321 }
01322
01323 if (!returnReduce(changed, orOp<bool>()))
01324 {
01325 break;
01326 }
01327 }
01328
01329
01330 forAll(regionToCellZone, regionI)
01331 {
01332 label zoneI = regionToCellZone[regionI];
01333
01334 if (zoneI == -2)
01335 {
01336 FatalErrorIn
01337 (
01338 "meshRefinement::findCellZoneTopo"
01339 "(const point&, const labelList&, const labelList&, labelList&)"
01340 ) << "For region " << regionI << " haven't set cell zone."
01341 << exit(FatalError);
01342 }
01343 }
01344
01345 if (debug)
01346 {
01347 forAll(regionToCellZone, regionI)
01348 {
01349 Pout<< "Region " << regionI
01350 << " becomes cellZone:" << regionToCellZone[regionI]
01351 << endl;
01352 }
01353 }
01354
01355
01356 forAll(cellToZone, cellI)
01357 {
01358 cellToZone[cellI] = regionToCellZone[cellRegion[cellI]];
01359 }
01360 }
01361
01362
01363
01364
01365 void Foam::meshRefinement::makeConsistentFaceIndex
01366 (
01367 const labelList& cellToZone,
01368 labelList& namedSurfaceIndex
01369 ) const
01370 {
01371 const labelList& faceOwner = mesh_.faceOwner();
01372 const labelList& faceNeighbour = mesh_.faceNeighbour();
01373
01374 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
01375 {
01376 label ownZone = cellToZone[faceOwner[faceI]];
01377 label neiZone = cellToZone[faceNeighbour[faceI]];
01378
01379 if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
01380 {
01381 namedSurfaceIndex[faceI] = -1;
01382 }
01383 else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
01384 {
01385 FatalErrorIn("meshRefinement::zonify()")
01386 << "Different cell zones on either side of face " << faceI
01387 << " at " << mesh_.faceCentres()[faceI]
01388 << " but face not marked with a surface."
01389 << abort(FatalError);
01390 }
01391 }
01392
01393 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
01394
01395
01396 labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
01397 forAll(patches, patchI)
01398 {
01399 const polyPatch& pp = patches[patchI];
01400
01401 if (pp.coupled())
01402 {
01403 forAll(pp, i)
01404 {
01405 label faceI = pp.start()+i;
01406 neiCellZone[faceI-mesh_.nInternalFaces()] =
01407 cellToZone[mesh_.faceOwner()[faceI]];
01408 }
01409 }
01410 }
01411 syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
01412
01413
01414 forAll(patches, patchI)
01415 {
01416 const polyPatch& pp = patches[patchI];
01417
01418 if (pp.coupled())
01419 {
01420 forAll(pp, i)
01421 {
01422 label faceI = pp.start()+i;
01423
01424 label ownZone = cellToZone[faceOwner[faceI]];
01425 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
01426
01427 if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
01428 {
01429 namedSurfaceIndex[faceI] = -1;
01430 }
01431 else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
01432 {
01433 FatalErrorIn("meshRefinement::zonify()")
01434 << "Different cell zones on either side of face "
01435 << faceI << " at " << mesh_.faceCentres()[faceI]
01436 << " but face not marked with a surface."
01437 << abort(FatalError);
01438 }
01439 }
01440 }
01441 }
01442 }
01443
01444
01445
01446
01447
01448 void Foam::meshRefinement::baffleAndSplitMesh
01449 (
01450 const bool handleSnapProblems,
01451 const bool removeEdgeConnectedCells,
01452 const scalarField& perpendicularAngle,
01453 const bool mergeFreeStanding,
01454 const dictionary& motionDict,
01455 Time& runTime,
01456 const labelList& globalToPatch,
01457 const point& keepPoint
01458 )
01459 {
01460
01461
01462
01463
01464
01465
01466 Info<< "Introducing baffles for "
01467 << returnReduce(countHits(), sumOp<label>())
01468 << " faces that are intersected by the surface." << nl << endl;
01469
01470
01471 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
01472 pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
01473 calcNeighbourData(neiLevel, neiCc);
01474
01475 labelList ownPatch, neiPatch;
01476 getBafflePatches
01477 (
01478 globalToPatch,
01479 neiLevel,
01480 neiCc,
01481
01482 ownPatch,
01483 neiPatch
01484 );
01485
01486 if (debug)
01487 {
01488 runTime++;
01489 }
01490
01491 createBaffles(ownPatch, neiPatch);
01492
01493 if (debug)
01494 {
01495
01496 checkData();
01497 }
01498
01499 Info<< "Created baffles in = "
01500 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01501
01502 printMeshInfo(debug, "After introducing baffles");
01503
01504 if (debug)
01505 {
01506 Pout<< "Writing baffled mesh to time " << timeName()
01507 << endl;
01508 write(debug, runTime.path()/"baffles");
01509 Pout<< "Dumped debug data in = "
01510 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01511 }
01512
01513
01514
01515
01516
01517
01518
01519 if (handleSnapProblems)
01520 {
01521 Info<< nl
01522 << "Introducing baffles to block off problem cells" << nl
01523 << "----------------------------------------------" << nl
01524 << endl;
01525
01526 labelList facePatch
01527 (
01528 markFacesOnProblemCells
01529 (
01530 motionDict,
01531 removeEdgeConnectedCells,
01532 perpendicularAngle,
01533 globalToPatch
01534 )
01535
01536 );
01537 Info<< "Analyzed problem cells in = "
01538 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01539
01540 if (debug)
01541 {
01542 faceSet problemTopo(mesh_, "problemFacesTopo", 100);
01543
01544 const labelList facePatchTopo
01545 (
01546 markFacesOnProblemCells
01547 (
01548 motionDict,
01549 removeEdgeConnectedCells,
01550 perpendicularAngle,
01551 globalToPatch
01552 )
01553 );
01554 forAll(facePatchTopo, faceI)
01555 {
01556 if (facePatchTopo[faceI] != -1)
01557 {
01558 problemTopo.insert(faceI);
01559 }
01560 }
01561 Pout<< "Dumping " << problemTopo.size()
01562 << " problem faces to " << problemTopo.objectPath() << endl;
01563 problemTopo.write();
01564 }
01565
01566 Info<< "Introducing baffles to delete problem cells." << nl << endl;
01567
01568 if (debug)
01569 {
01570 runTime++;
01571 }
01572
01573
01574 createBaffles(facePatch, facePatch);
01575
01576 if (debug)
01577 {
01578
01579 checkData();
01580 }
01581 Info<< "Created baffles in = "
01582 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01583
01584 printMeshInfo(debug, "After introducing baffles");
01585
01586 if (debug)
01587 {
01588 Pout<< "Writing extra baffled mesh to time "
01589 << timeName() << endl;
01590 write(debug, runTime.path()/"extraBaffles");
01591 Pout<< "Dumped debug data in = "
01592 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01593 }
01594 }
01595
01596
01597
01598
01599
01600 Info<< nl
01601 << "Remove unreachable sections of mesh" << nl
01602 << "-----------------------------------" << nl
01603 << endl;
01604
01605 if (debug)
01606 {
01607 runTime++;
01608 }
01609
01610 splitMeshRegions(keepPoint);
01611
01612 if (debug)
01613 {
01614
01615 checkData();
01616 }
01617 Info<< "Split mesh in = "
01618 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01619
01620 printMeshInfo(debug, "After subsetting");
01621
01622 if (debug)
01623 {
01624 Pout<< "Writing subsetted mesh to time " << timeName()
01625 << endl;
01626 write(debug, runTime.path()/timeName());
01627 Pout<< "Dumped debug data in = "
01628 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01629 }
01630
01631
01632
01633
01634
01635 if (mergeFreeStanding)
01636 {
01637 Info<< nl
01638 << "Merge free-standing baffles" << nl
01639 << "---------------------------" << nl
01640 << endl;
01641
01642 if (debug)
01643 {
01644 runTime++;
01645 }
01646
01647
01648 List<labelPair> couples
01649 (
01650 filterDuplicateFaces
01651 (
01652 getDuplicateFaces
01653 (
01654 identity(mesh_.nFaces()-mesh_.nInternalFaces())
01655 +mesh_.nInternalFaces()
01656 )
01657 )
01658 );
01659
01660 label nCouples = couples.size();
01661 reduce(nCouples, sumOp<label>());
01662
01663 Info<< "Detected free-standing baffles : " << nCouples << endl;
01664
01665 if (nCouples > 0)
01666 {
01667
01668
01669
01670 mergeBaffles(couples);
01671
01672 if (debug)
01673 {
01674
01675 checkData();
01676 }
01677 }
01678 Info<< "Merged free-standing baffles in = "
01679 << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
01680 }
01681 }
01682
01683
01684
01685 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
01686 (
01687 const label nBufferLayers,
01688 const labelList& globalToPatch,
01689 const point& keepPoint
01690 )
01691 {
01692
01693
01694
01695
01696 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
01697 pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
01698 calcNeighbourData(neiLevel, neiCc);
01699
01700 labelList ownPatch, neiPatch;
01701 getBafflePatches
01702 (
01703 globalToPatch,
01704 neiLevel,
01705 neiCc,
01706
01707 ownPatch,
01708 neiPatch
01709 );
01710
01711
01712 boolList blockedFace(mesh_.nFaces(), false);
01713
01714 forAll(ownPatch, faceI)
01715 {
01716 if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
01717 {
01718 blockedFace[faceI] = true;
01719 }
01720 }
01721 syncTools::syncFaceList(mesh_, blockedFace, orEqOp<bool>(), false);
01722
01723
01724 regionSplit cellRegion(mesh_, blockedFace);
01725 blockedFace.clear();
01726
01727
01728 label keepRegionI = -1;
01729
01730 label cellI = mesh_.findCell(keepPoint);
01731
01732 if (cellI != -1)
01733 {
01734 keepRegionI = cellRegion[cellI];
01735 }
01736 reduce(keepRegionI, maxOp<label>());
01737
01738 Info<< "Found point " << keepPoint << " in cell " << cellI
01739 << " in global region " << keepRegionI
01740 << " out of " << cellRegion.nRegions() << " regions." << endl;
01741
01742 if (keepRegionI == -1)
01743 {
01744 FatalErrorIn
01745 (
01746 "meshRefinement::splitMesh"
01747 "(const label, const labelList&, const point&)"
01748 ) << "Point " << keepPoint
01749 << " is not inside the mesh." << nl
01750 << "Bounding box of the mesh:" << mesh_.globalData().bb()
01751 << exit(FatalError);
01752 }
01753
01754
01755
01756
01757
01758
01759
01760
01761 const labelList& faceOwner = mesh_.faceOwner();
01762 const labelList& faceNeighbour = mesh_.faceNeighbour();
01763
01764
01765 label defaultPatch = 0;
01766 if (globalToPatch.size())
01767 {
01768 defaultPatch = globalToPatch[0];
01769 }
01770
01771 for (label i = 0; i < nBufferLayers; i++)
01772 {
01773
01774
01775 labelList pointBaffle(mesh_.nPoints(), -1);
01776
01777 forAll(faceNeighbour, faceI)
01778 {
01779 const face& f = mesh_.faces()[faceI];
01780
01781 label ownRegion = cellRegion[faceOwner[faceI]];
01782 label neiRegion = cellRegion[faceNeighbour[faceI]];
01783
01784 if (ownRegion == keepRegionI && neiRegion != keepRegionI)
01785 {
01786
01787
01788
01789 forAll(f, fp)
01790 {
01791 pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
01792 }
01793 }
01794 else if (ownRegion != keepRegionI && neiRegion == keepRegionI)
01795 {
01796 label newPatchI = neiPatch[faceI];
01797 if (newPatchI == -1)
01798 {
01799 newPatchI = max(defaultPatch, ownPatch[faceI]);
01800 }
01801 forAll(f, fp)
01802 {
01803 pointBaffle[f[fp]] = newPatchI;
01804 }
01805 }
01806 }
01807 for
01808 (
01809 label faceI = mesh_.nInternalFaces();
01810 faceI < mesh_.nFaces();
01811 faceI++
01812 )
01813 {
01814 const face& f = mesh_.faces()[faceI];
01815
01816 label ownRegion = cellRegion[faceOwner[faceI]];
01817
01818 if (ownRegion == keepRegionI)
01819 {
01820 forAll(f, fp)
01821 {
01822 pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
01823 }
01824 }
01825 }
01826
01827
01828 syncTools::syncPointList
01829 (
01830 mesh_,
01831 pointBaffle,
01832 maxEqOp<label>(),
01833 -1,
01834 false
01835 );
01836
01837
01838
01839
01840 const labelListList& pointFaces = mesh_.pointFaces();
01841
01842 forAll(pointFaces, pointI)
01843 {
01844 if (pointBaffle[pointI] != -1)
01845 {
01846 const labelList& pFaces = pointFaces[pointI];
01847
01848 forAll(pFaces, pFaceI)
01849 {
01850 label faceI = pFaces[pFaceI];
01851
01852 if (ownPatch[faceI] == -1)
01853 {
01854 ownPatch[faceI] = pointBaffle[pointI];
01855 }
01856 }
01857 }
01858 }
01859 syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
01860
01861
01862
01863
01864 labelList newOwnPatch(ownPatch);
01865
01866 forAll(ownPatch, faceI)
01867 {
01868 if (ownPatch[faceI] != -1)
01869 {
01870 label own = faceOwner[faceI];
01871
01872 if (cellRegion[own] != keepRegionI)
01873 {
01874 cellRegion[own] = keepRegionI;
01875
01876 const cell& ownFaces = mesh_.cells()[own];
01877 forAll(ownFaces, j)
01878 {
01879 if (ownPatch[ownFaces[j]] == -1)
01880 {
01881 newOwnPatch[ownFaces[j]] = ownPatch[faceI];
01882 }
01883 }
01884 }
01885 if (mesh_.isInternalFace(faceI))
01886 {
01887 label nei = faceNeighbour[faceI];
01888
01889 if (cellRegion[nei] != keepRegionI)
01890 {
01891 cellRegion[nei] = keepRegionI;
01892
01893 const cell& neiFaces = mesh_.cells()[nei];
01894 forAll(neiFaces, j)
01895 {
01896 if (ownPatch[neiFaces[j]] == -1)
01897 {
01898 newOwnPatch[neiFaces[j]] = ownPatch[faceI];
01899 }
01900 }
01901 }
01902 }
01903 }
01904 }
01905
01906 ownPatch.transfer(newOwnPatch);
01907
01908 syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
01909 }
01910
01911
01912
01913
01914
01915
01916
01917 DynamicList<label> cellsToRemove(mesh_.nCells());
01918 forAll(cellRegion, cellI)
01919 {
01920 if (cellRegion[cellI] != keepRegionI)
01921 {
01922 cellsToRemove.append(cellI);
01923 }
01924 }
01925 cellsToRemove.shrink();
01926
01927 label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
01928 reduce(nCellsToKeep, sumOp<label>());
01929
01930 Info<< "Keeping all cells in region " << keepRegionI
01931 << " containing point " << keepPoint << endl
01932 << "Selected for keeping : " << nCellsToKeep
01933 << " cells." << endl;
01934
01935
01936
01937 removeCells cellRemover(mesh_);
01938
01939
01940 labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
01941 labelList exposedPatches(exposedFaces.size());
01942
01943 forAll(exposedFaces, i)
01944 {
01945 label faceI = exposedFaces[i];
01946
01947 if (ownPatch[faceI] != -1)
01948 {
01949 exposedPatches[i] = ownPatch[faceI];
01950 }
01951 else
01952 {
01953 WarningIn("meshRefinement::splitMesh(..)")
01954 << "For exposed face " << faceI
01955 << " fc:" << mesh_.faceCentres()[faceI]
01956 << " found no patch." << endl
01957 << " Taking patch " << defaultPatch
01958 << " instead." << endl;
01959 exposedPatches[i] = defaultPatch;
01960 }
01961 }
01962
01963 return doRemoveCells
01964 (
01965 cellsToRemove,
01966 exposedFaces,
01967 exposedPatches,
01968 cellRemover
01969 );
01970 }
01971
01972
01973
01974
01975 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
01976 {
01977
01978 polyTopoChange meshMod(mesh_);
01979
01980
01981
01982 localPointRegion regionSide(mesh_);
01983
01984 label nNonManifPoints = returnReduce
01985 (
01986 regionSide.meshPointMap().size(),
01987 sumOp<label>()
01988 );
01989
01990 Info<< "dupNonManifoldPoints : Found : " << nNonManifPoints
01991 << " non-manifold points (out of "
01992 << mesh_.globalData().nTotalPoints()
01993 << ')' << endl;
01994
01995
01996 duplicatePoints pointDuplicator(mesh_);
01997
01998
01999 pointDuplicator.setRefinement(regionSide, meshMod);
02000
02001
02002 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
02003
02004
02005 mesh_.updateMesh(map);
02006
02007
02008 if (map().hasMotionPoints())
02009 {
02010 mesh_.movePoints(map().preMotionPoints());
02011 }
02012 else
02013 {
02014
02015 mesh_.clearOut();
02016 }
02017
02018 if (overwrite())
02019 {
02020 mesh_.setInstance(oldInstance());
02021 }
02022
02023
02024
02025 updateMesh(map, labelList(0));
02026
02027 return map;
02028 }
02029
02030
02031
02032 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
02033 (
02034 const point& keepPoint,
02035 const bool allowFreeStandingZoneFaces
02036 )
02037 {
02038 const wordList& cellZoneNames = surfaces_.cellZoneNames();
02039 const wordList& faceZoneNames = surfaces_.faceZoneNames();
02040
02041 labelList namedSurfaces(surfaces_.getNamedSurfaces());
02042
02043 boolList isNamedSurface(cellZoneNames.size(), false);
02044
02045 forAll(namedSurfaces, i)
02046 {
02047 label surfI = namedSurfaces[i];
02048
02049 isNamedSurface[surfI] = true;
02050
02051 Info<< "Surface : " << surfaces_.names()[surfI] << nl
02052 << " faceZone : " << faceZoneNames[surfI] << nl
02053 << " cellZone : " << cellZoneNames[surfI] << endl;
02054 }
02055
02056
02057
02058
02059 labelList surfaceToFaceZone(faceZoneNames.size(), -1);
02060 {
02061 faceZoneMesh& faceZones = mesh_.faceZones();
02062
02063 forAll(namedSurfaces, i)
02064 {
02065 label surfI = namedSurfaces[i];
02066
02067 label zoneI = faceZones.findZoneID(faceZoneNames[surfI]);
02068
02069 if (zoneI == -1)
02070 {
02071 zoneI = faceZones.size();
02072 faceZones.setSize(zoneI+1);
02073 faceZones.set
02074 (
02075 zoneI,
02076 new faceZone
02077 (
02078 faceZoneNames[surfI],
02079 labelList(0),
02080 boolList(0),
02081 zoneI,
02082 faceZones
02083 )
02084 );
02085 }
02086
02087 if (debug)
02088 {
02089 Pout<< "Faces on " << surfaces_.names()[surfI]
02090 << " will go into faceZone " << zoneI << endl;
02091 }
02092 surfaceToFaceZone[surfI] = zoneI;
02093 }
02094
02095
02096 List<wordList> allFaceZones(Pstream::nProcs());
02097 allFaceZones[Pstream::myProcNo()] = faceZones.names();
02098 Pstream::gatherList(allFaceZones);
02099 Pstream::scatterList(allFaceZones);
02100
02101 for (label procI = 1; procI < allFaceZones.size(); procI++)
02102 {
02103 if (allFaceZones[procI] != allFaceZones[0])
02104 {
02105 FatalErrorIn
02106 (
02107 "meshRefinement::zonify"
02108 "(const label, const point&)"
02109 ) << "Zones not synchronised among processors." << nl
02110 << " Processor0 has faceZones:" << allFaceZones[0]
02111 << " , processor" << procI
02112 << " has faceZones:" << allFaceZones[procI]
02113 << exit(FatalError);
02114 }
02115 }
02116 }
02117
02118 labelList surfaceToCellZone(cellZoneNames.size(), -1);
02119 {
02120 cellZoneMesh& cellZones = mesh_.cellZones();
02121
02122 forAll(namedSurfaces, i)
02123 {
02124 label surfI = namedSurfaces[i];
02125
02126 label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
02127
02128 if (zoneI == -1)
02129 {
02130 zoneI = cellZones.size();
02131 cellZones.setSize(zoneI+1);
02132 cellZones.set
02133 (
02134 zoneI,
02135 new cellZone
02136 (
02137 cellZoneNames[surfI],
02138 labelList(0),
02139 zoneI,
02140 cellZones
02141 )
02142 );
02143 }
02144
02145 if (debug)
02146 {
02147 Pout<< "Cells inside " << surfaces_.names()[surfI]
02148 << " will go into cellZone " << zoneI << endl;
02149 }
02150 surfaceToCellZone[surfI] = zoneI;
02151 }
02152
02153
02154 List<wordList> allCellZones(Pstream::nProcs());
02155 allCellZones[Pstream::myProcNo()] = cellZones.names();
02156 Pstream::gatherList(allCellZones);
02157 Pstream::scatterList(allCellZones);
02158
02159 for (label procI = 1; procI < allCellZones.size(); procI++)
02160 {
02161 if (allCellZones[procI] != allCellZones[0])
02162 {
02163 FatalErrorIn
02164 (
02165 "meshRefinement::zonify"
02166 "(const label, const point&)"
02167 ) << "Zones not synchronised among processors." << nl
02168 << " Processor0 has cellZones:" << allCellZones[0]
02169 << " , processor" << procI
02170 << " has cellZones:" << allCellZones[procI]
02171 << exit(FatalError);
02172 }
02173 }
02174 }
02175
02176
02177
02178 const pointField& cellCentres = mesh_.cellCentres();
02179 const labelList& faceOwner = mesh_.faceOwner();
02180 const labelList& faceNeighbour = mesh_.faceNeighbour();
02181 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
02182
02183
02184
02185
02186
02187
02188
02189 labelList namedSurfaceIndex(mesh_.nFaces(), -1);
02190
02191 {
02192
02193 labelList nSurfFaces(faceZoneNames.size(), 0);
02194
02195
02196 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
02197 pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
02198 calcNeighbourData(neiLevel, neiCc);
02199
02200
02201
02202
02203
02204
02205
02206
02207
02208 labelList testFaces(intersectedFaces());
02209
02210
02211
02212
02213 pointField start(testFaces.size());
02214 pointField end(testFaces.size());
02215
02216 forAll(testFaces, i)
02217 {
02218 label faceI = testFaces[i];
02219
02220 if (mesh_.isInternalFace(faceI))
02221 {
02222 start[i] = cellCentres[faceOwner[faceI]];
02223 end[i] = cellCentres[faceNeighbour[faceI]];
02224 }
02225 else
02226 {
02227 start[i] = cellCentres[faceOwner[faceI]];
02228 end[i] = neiCc[faceI-mesh_.nInternalFaces()];
02229 }
02230 }
02231
02232
02233
02234
02235
02236
02237
02238 labelList surface1;
02239 labelList surface2;
02240 {
02241 List<pointIndexHit> hit1;
02242 labelList region1;
02243 List<pointIndexHit> hit2;
02244 labelList region2;
02245 surfaces_.findNearestIntersection
02246 (
02247 namedSurfaces,
02248 start,
02249 end,
02250
02251 surface1,
02252 hit1,
02253 region1,
02254 surface2,
02255 hit2,
02256 region2
02257 );
02258 }
02259
02260 forAll(testFaces, i)
02261 {
02262 label faceI = testFaces[i];
02263
02264 if (surface1[i] != -1)
02265 {
02266
02267 namedSurfaceIndex[faceI] = surface1[i];
02268 nSurfFaces[surface1[i]]++;
02269 }
02270 else if (surface2[i] != -1)
02271 {
02272 namedSurfaceIndex[faceI] = surface2[i];
02273 nSurfFaces[surface2[i]]++;
02274 }
02275 }
02276
02277
02278
02279
02280
02281
02282 syncTools::syncFaceList
02283 (
02284 mesh_,
02285 namedSurfaceIndex,
02286 maxEqOp<label>(),
02287 false
02288 );
02289
02290
02291 if (debug)
02292 {
02293 forAll(nSurfFaces, surfI)
02294 {
02295 Pout<< "Surface:"
02296 << surfaces_.names()[surfI]
02297 << " nZoneFaces:" << nSurfFaces[surfI] << nl;
02298 }
02299 Pout<< endl;
02300 }
02301 }
02302
02303
02304
02305
02306
02307
02308 labelList closedNamedSurfaces(surfaces_.getClosedNamedSurfaces());
02309
02310
02311
02312
02313
02314 labelList cellToZone(mesh_.nCells(), -2);
02315
02316
02317
02318
02319
02320 if (closedNamedSurfaces.size())
02321 {
02322 findCellZoneGeometric
02323 (
02324 closedNamedSurfaces,
02325 namedSurfaceIndex,
02326 surfaceToCellZone,
02327 cellToZone
02328 );
02329 }
02330
02331
02332
02333
02334
02335 {
02336
02337 findCellZoneTopo
02338 (
02339 keepPoint,
02340 namedSurfaceIndex,
02341 surfaceToCellZone,
02342 cellToZone
02343 );
02344 }
02345
02346
02347
02348 if (!allowFreeStandingZoneFaces)
02349 {
02350 makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
02351 }
02352
02353
02354 polyTopoChange meshMod(mesh_);
02355
02356
02357
02358
02359
02360 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
02361 {
02362 label surfI = namedSurfaceIndex[faceI];
02363
02364 if (surfI != -1)
02365 {
02366
02367 label ownZone = cellToZone[faceOwner[faceI]];
02368 label neiZone = cellToZone[faceNeighbour[faceI]];
02369
02370 bool flip;
02371 if (ownZone == max(ownZone, neiZone))
02372 {
02373 flip = false;
02374 }
02375 else
02376 {
02377 flip = true;
02378 }
02379
02380 meshMod.setAction
02381 (
02382 polyModifyFace
02383 (
02384 mesh_.faces()[faceI],
02385 faceI,
02386 faceOwner[faceI],
02387 faceNeighbour[faceI],
02388 false,
02389 -1,
02390 false,
02391 surfaceToFaceZone[surfI],
02392 flip
02393 )
02394 );
02395 }
02396 }
02397
02398
02399 labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
02400 forAll(patches, patchI)
02401 {
02402 const polyPatch& pp = patches[patchI];
02403
02404 if (pp.coupled())
02405 {
02406 forAll(pp, i)
02407 {
02408 label faceI = pp.start()+i;
02409 neiCellZone[faceI-mesh_.nInternalFaces()] =
02410 cellToZone[mesh_.faceOwner()[faceI]];
02411 }
02412 }
02413 }
02414 syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
02415
02416
02417 PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
02418
02419
02420 forAll(patches, patchI)
02421 {
02422 const polyPatch& pp = patches[patchI];
02423
02424 label faceI = pp.start();
02425
02426 forAll(pp, i)
02427 {
02428 label surfI = namedSurfaceIndex[faceI];
02429
02430 if (surfI != -1)
02431 {
02432 label ownZone = cellToZone[faceOwner[faceI]];
02433 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
02434
02435 bool flip;
02436
02437 label maxZone = max(ownZone, neiZone);
02438
02439 if (maxZone == -1)
02440 {
02441 flip = false;
02442 }
02443 else if (ownZone == neiZone)
02444 {
02445
02446
02447 flip = !isMasterFace[faceI];
02448 }
02449 else if (neiZone == maxZone)
02450 {
02451 flip = true;
02452 }
02453 else
02454 {
02455 flip = false;
02456 }
02457
02458 meshMod.setAction
02459 (
02460 polyModifyFace
02461 (
02462 mesh_.faces()[faceI],
02463 faceI,
02464 faceOwner[faceI],
02465 -1,
02466 false,
02467 patchI,
02468 false,
02469 surfaceToFaceZone[surfI],
02470 flip
02471 )
02472 );
02473 }
02474 faceI++;
02475 }
02476 }
02477
02478
02479
02480
02481
02482 forAll(cellToZone, cellI)
02483 {
02484 label zoneI = cellToZone[cellI];
02485
02486 if (zoneI >= 0)
02487 {
02488 meshMod.setAction
02489 (
02490 polyModifyCell
02491 (
02492 cellI,
02493 false,
02494 zoneI
02495 )
02496 );
02497 }
02498 }
02499
02500
02501 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
02502
02503
02504 mesh_.updateMesh(map);
02505
02506
02507 if (map().hasMotionPoints())
02508 {
02509 mesh_.movePoints(map().preMotionPoints());
02510 }
02511 else
02512 {
02513
02514 mesh_.clearOut();
02515 }
02516
02517 if (overwrite())
02518 {
02519 mesh_.setInstance(oldInstance());
02520 }
02521
02522
02523 updateMesh(map, labelList());
02524
02525 return map;
02526 }
02527
02528
02529