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
00027
00028 #include "cellClassification.H"
00029 #include <meshTools/triSurfaceSearch.H>
00030 #include <meshTools/indexedOctree.H>
00031 #include <meshTools/treeDataFace.H>
00032 #include <meshTools/meshSearch.H>
00033 #include "cellInfo.H"
00034 #include <OpenFOAM/polyMesh.H>
00035 #include <OpenFOAM/MeshWave.H>
00036 #include <OpenFOAM/ListOps.H>
00037 #include <meshTools/meshTools.H>
00038
00039
00040
00041 namespace Foam
00042 {
00043 defineTypeNameAndDebug(cellClassification, 0);
00044 }
00045
00046
00047
00048
00049 Foam::label Foam::cellClassification::count
00050 (
00051 const labelList& elems,
00052 const label elem
00053 )
00054 {
00055 label cnt = 0;
00056
00057 forAll(elems, elemI)
00058 {
00059 if (elems[elemI] == elem)
00060 {
00061 cnt++;
00062 }
00063 }
00064 return cnt;
00065
00066 }
00067
00068
00069
00070
00071
00072
00073
00074 Foam::boolList Foam::cellClassification::markFaces
00075 (
00076 const triSurfaceSearch& search
00077 ) const
00078 {
00079 cpuTime timer;
00080
00081 boolList cutFace(mesh_.nFaces(), false);
00082
00083 label nCutFaces = 0;
00084
00085
00086
00087
00088 forAll(mesh_.edges(), edgeI)
00089 {
00090 if (debug && (edgeI % 10000 == 0))
00091 {
00092 Pout<< "Intersecting mesh edge " << edgeI << " with surface"
00093 << endl;
00094 }
00095
00096 const edge& e = mesh_.edges()[edgeI];
00097
00098 const point& p0 = mesh_.points()[e.start()];
00099 const point& p1 = mesh_.points()[e.end()];
00100
00101 pointIndexHit pHit(search.tree().findLineAny(p0, p1));
00102
00103 if (pHit.hit())
00104 {
00105 const labelList& myFaces = mesh_.edgeFaces()[edgeI];
00106
00107 forAll(myFaces, myFaceI)
00108 {
00109 label faceI = myFaces[myFaceI];
00110
00111 if (!cutFace[faceI])
00112 {
00113 cutFace[faceI] = true;
00114
00115 nCutFaces++;
00116 }
00117 }
00118 }
00119 }
00120
00121 if (debug)
00122 {
00123 Pout<< "Intersected edges of mesh with surface in = "
00124 << timer.cpuTimeIncrement() << " s\n" << endl << endl;
00125 }
00126
00127
00128
00129
00130
00131 labelList allFaces(mesh_.nFaces() - nCutFaces);
00132
00133 label allFaceI = 0;
00134
00135 forAll(cutFace, faceI)
00136 {
00137 if (!cutFace[faceI])
00138 {
00139 allFaces[allFaceI++] = faceI;
00140 }
00141 }
00142
00143 if (debug)
00144 {
00145 Pout<< "Testing " << allFaceI << " faces for piercing by surface"
00146 << endl;
00147 }
00148
00149 treeBoundBox allBb(mesh_.points());
00150
00151
00152 scalar tol = 1E-6 * allBb.avgDim();
00153
00154 point& bbMin = allBb.min();
00155 bbMin.x() -= tol;
00156 bbMin.y() -= tol;
00157 bbMin.z() -= tol;
00158
00159 point& bbMax = allBb.max();
00160 bbMax.x() += 2*tol;
00161 bbMax.y() += 2*tol;
00162 bbMax.z() += 2*tol;
00163
00164 indexedOctree<treeDataFace> faceTree
00165 (
00166 treeDataFace(false, mesh_, allFaces),
00167 allBb,
00168 8,
00169 10,
00170 3.0
00171 );
00172
00173 const triSurface& surf = search.surface();
00174 const edgeList& edges = surf.edges();
00175 const pointField& localPoints = surf.localPoints();
00176
00177 label nAddFaces = 0;
00178
00179 forAll(edges, edgeI)
00180 {
00181 if (debug && (edgeI % 10000 == 0))
00182 {
00183 Pout<< "Intersecting surface edge " << edgeI
00184 << " with mesh faces" << endl;
00185 }
00186 const edge& e = edges[edgeI];
00187
00188 const point& start = localPoints[e.start()];
00189 const point& end = localPoints[e.end()];
00190
00191 vector edgeNormal(end - start);
00192 const scalar edgeMag = mag(edgeNormal);
00193 const vector smallVec = 1E-9*edgeNormal;
00194
00195 edgeNormal /= edgeMag+VSMALL;
00196
00197
00198 point pt = start;
00199
00200 while (true)
00201 {
00202 pointIndexHit pHit(faceTree.findLine(pt, end));
00203
00204 if (!pHit.hit())
00205 {
00206 break;
00207 }
00208 else
00209 {
00210 label faceI = faceTree.shapes().faceLabels()[pHit.index()];
00211
00212 if (!cutFace[faceI])
00213 {
00214 cutFace[faceI] = true;
00215
00216 nAddFaces++;
00217 }
00218
00219
00220 pt = pHit.hitPoint() + smallVec;
00221
00222 if (((pt-start) & edgeNormal) >= edgeMag)
00223 {
00224 break;
00225 }
00226 }
00227 }
00228 }
00229
00230 if (debug)
00231 {
00232 Pout<< "Detected an additional " << nAddFaces << " faces cut"
00233 << endl;
00234
00235 Pout<< "Intersected edges of surface with mesh faces in = "
00236 << timer.cpuTimeIncrement() << " s\n" << endl << endl;
00237 }
00238
00239 return cutFace;
00240 }
00241
00242
00243
00244
00245
00246 void Foam::cellClassification::markCells
00247 (
00248 const meshSearch& queryMesh,
00249 const boolList& piercedFace,
00250 const pointField& outsidePts
00251 )
00252 {
00253
00254
00255
00256
00257 List<cellInfo> cellInfoList(mesh_.nCells());
00258
00259
00260 forAll(piercedFace, faceI)
00261 {
00262 if (piercedFace[faceI])
00263 {
00264 cellInfoList[mesh_.faceOwner()[faceI]] =
00265 cellInfo(cellClassification::CUT);
00266
00267 if (mesh_.isInternalFace(faceI))
00268 {
00269 cellInfoList[mesh_.faceNeighbour()[faceI]] =
00270 cellInfo(cellClassification::CUT);
00271 }
00272 }
00273 }
00274
00275
00276
00277
00278
00279
00280 labelHashSet outsideFacesMap(outsidePts.size() * 6 * 2);
00281
00282 forAll(outsidePts, outsidePtI)
00283 {
00284
00285 label cellI = queryMesh.findCell(outsidePts[outsidePtI], -1, false);
00286
00287 if (returnReduce(cellI, maxOp<label>()) == -1)
00288 {
00289 FatalErrorIn
00290 (
00291 "List<cellClassification::cType> markCells"
00292 "(const meshSearch&, const boolList&, const pointField&)"
00293 ) << "outsidePoint " << outsidePts[outsidePtI]
00294 << " is not inside any cell"
00295 << nl << "It might be on a face or outside the geometry"
00296 << exit(FatalError);
00297 }
00298
00299 if (cellI >= 0)
00300 {
00301 cellInfoList[cellI] = cellInfo(cellClassification::OUTSIDE);
00302
00303
00304 const labelList& myFaces = mesh_.cells()[cellI];
00305 forAll(myFaces, myFaceI)
00306 {
00307 outsideFacesMap.insert(myFaces[myFaceI]);
00308 }
00309 }
00310 }
00311
00312
00313
00314
00315
00316
00317 labelList changedFaces(outsideFacesMap.toc());
00318
00319 List<cellInfo> changedFacesInfo
00320 (
00321 changedFaces.size(),
00322 cellInfo(cellClassification::OUTSIDE)
00323 );
00324
00325 MeshWave<cellInfo> cellInfoCalc
00326 (
00327 mesh_,
00328 changedFaces,
00329 changedFacesInfo,
00330 cellInfoList,
00331 mesh_.globalData().nTotalCells()
00332 );
00333
00334
00335 const List<cellInfo>& allInfo = cellInfoCalc.allCellInfo();
00336
00337 forAll(allInfo, cellI)
00338 {
00339 label t = allInfo[cellI].type();
00340
00341 if (t == cellClassification::NOTSET)
00342 {
00343 t = cellClassification::INSIDE;
00344 }
00345 operator[](cellI) = t;
00346 }
00347 }
00348
00349
00350 void Foam::cellClassification::classifyPoints
00351 (
00352 const label meshType,
00353 const labelList& cellType,
00354 List<pointStatus>& pointSide
00355 ) const
00356 {
00357 pointSide.setSize(mesh_.nPoints());
00358
00359 forAll(mesh_.pointCells(), pointI)
00360 {
00361 const labelList& pCells = mesh_.pointCells()[pointI];
00362
00363 pointSide[pointI] = UNSET;
00364
00365 forAll(pCells, i)
00366 {
00367 label type = cellType[pCells[i]];
00368
00369 if (type == meshType)
00370 {
00371 if (pointSide[pointI] == UNSET)
00372 {
00373 pointSide[pointI] = MESH;
00374 }
00375 else if (pointSide[pointI] == NONMESH)
00376 {
00377 pointSide[pointI] = MIXED;
00378
00379 break;
00380 }
00381 }
00382 else
00383 {
00384 if (pointSide[pointI] == UNSET)
00385 {
00386 pointSide[pointI] = NONMESH;
00387 }
00388 else if (pointSide[pointI] == MESH)
00389 {
00390 pointSide[pointI] = MIXED;
00391
00392 break;
00393 }
00394 }
00395 }
00396 }
00397 }
00398
00399
00400 bool Foam::cellClassification::usesMixedPointsOnly
00401 (
00402 const List<pointStatus>& pointSide,
00403 const label cellI
00404 ) const
00405 {
00406 const faceList& faces = mesh_.faces();
00407
00408 const cell& cFaces = mesh_.cells()[cellI];
00409
00410 forAll(cFaces, cFaceI)
00411 {
00412 const face& f = faces[cFaces[cFaceI]];
00413
00414 forAll(f, fp)
00415 {
00416 if (pointSide[f[fp]] != MIXED)
00417 {
00418 return false;
00419 }
00420 }
00421 }
00422
00423
00424 return true;
00425 }
00426
00427
00428 void Foam::cellClassification::getMeshOutside
00429 (
00430 const label meshType,
00431 faceList& outsideFaces,
00432 labelList& outsideOwner
00433 ) const
00434 {
00435 const labelList& own = mesh_.faceOwner();
00436 const labelList& nbr = mesh_.faceNeighbour();
00437 const faceList& faces = mesh_.faces();
00438
00439 outsideFaces.setSize(mesh_.nFaces());
00440 outsideOwner.setSize(mesh_.nFaces());
00441 label outsideI = 0;
00442
00443
00444
00445 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
00446 {
00447 label ownType = operator[](own[faceI]);
00448 label nbrType = operator[](nbr[faceI]);
00449
00450 if (ownType == meshType && nbrType != meshType)
00451 {
00452 outsideFaces[outsideI] = faces[faceI];
00453 outsideOwner[outsideI] = own[faceI];
00454 outsideI++;
00455 }
00456 else if (ownType != meshType && nbrType == meshType)
00457 {
00458 outsideFaces[outsideI] = faces[faceI];
00459 outsideOwner[outsideI] = nbr[faceI];
00460 outsideI++;
00461 }
00462 }
00463
00464
00465
00466 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
00467 {
00468 if (operator[](own[faceI]) == meshType)
00469 {
00470 outsideFaces[outsideI] = faces[faceI];
00471 outsideOwner[outsideI] = own[faceI];
00472 outsideI++;
00473 }
00474 }
00475 outsideFaces.setSize(outsideI);
00476 outsideOwner.setSize(outsideI);
00477 }
00478
00479
00480
00481
00482
00483 Foam::cellClassification::cellClassification
00484 (
00485 const polyMesh& mesh,
00486 const meshSearch& meshQuery,
00487 const triSurfaceSearch& surfQuery,
00488 const pointField& outsidePoints
00489 )
00490 :
00491 labelList(mesh.nCells(), cellClassification::NOTSET),
00492 mesh_(mesh)
00493 {
00494 markCells
00495 (
00496 meshQuery,
00497 markFaces(surfQuery),
00498 outsidePoints
00499 );
00500 }
00501
00502
00503
00504 Foam::cellClassification::cellClassification
00505 (
00506 const polyMesh& mesh,
00507 const labelList& cellType
00508 )
00509 :
00510 labelList(cellType),
00511 mesh_(mesh)
00512 {
00513 if (mesh_.nCells() != size())
00514 {
00515 FatalErrorIn
00516 (
00517 "cellClassification::cellClassification"
00518 "(const polyMesh&, const labelList&)"
00519 ) << "Number of elements of cellType argument is not equal to the"
00520 << " number of cells" << abort(FatalError);
00521 }
00522 }
00523
00524
00525
00526 Foam::cellClassification::cellClassification(const cellClassification& cType)
00527 :
00528 labelList(cType),
00529 mesh_(cType.mesh())
00530 {}
00531
00532
00533
00534
00535
00536
00537
00538 Foam::label Foam::cellClassification::trimCutCells
00539 (
00540 const label nLayers,
00541 const label meshType,
00542 const label fillType
00543 )
00544 {
00545
00546 labelList newCellType(*this);
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563 newCellType = *this;
00564
00565
00566
00567 for (label iter = 0; iter < nLayers; iter++)
00568 {
00569
00570
00571
00572 List<pointStatus> pointSide(mesh_.nPoints());
00573 classifyPoints(meshType, newCellType, pointSide);
00574
00575
00576 forAll(pointSide, pointI)
00577 {
00578 if (pointSide[pointI] == MIXED)
00579 {
00580
00581 const labelList& pCells = mesh_.pointCells()[pointI];
00582
00583 forAll(pCells, i)
00584 {
00585 label type = newCellType[pCells[i]];
00586
00587 if (type == cellClassification::CUT)
00588 {
00589
00590
00591 newCellType[pCells[i]] = meshType;
00592 }
00593 }
00594 }
00595 }
00596 }
00597
00598
00599
00600
00601
00602
00603 label nChanged = 0;
00604
00605 forAll(newCellType, cellI)
00606 {
00607 if (operator[](cellI) == cellClassification::CUT)
00608 {
00609 if (newCellType[cellI] != meshType)
00610 {
00611
00612
00613 operator[](cellI) = fillType;
00614 nChanged++;
00615 }
00616 }
00617 }
00618
00619 return nChanged;
00620 }
00621
00622
00623
00624 Foam::label Foam::cellClassification::growSurface
00625 (
00626 const label meshType,
00627 const label fillType
00628 )
00629 {
00630 boolList hasMeshType(mesh_.nPoints(), false);
00631
00632
00633
00634 forAll(mesh_.pointCells(), pointI)
00635 {
00636 const labelList& myCells = mesh_.pointCells()[pointI];
00637
00638
00639 forAll(myCells, myCellI)
00640 {
00641 label type = operator[](myCells[myCellI]);
00642
00643 if (type == meshType)
00644 {
00645 hasMeshType[pointI] = true;
00646
00647 break;
00648 }
00649 }
00650 }
00651
00652
00653
00654 label nChanged = 0;
00655
00656 forAll(hasMeshType, pointI)
00657 {
00658 if (hasMeshType[pointI])
00659 {
00660 const labelList& myCells = mesh_.pointCells()[pointI];
00661
00662 forAll(myCells, myCellI)
00663 {
00664 if (operator[](myCells[myCellI]) != meshType)
00665 {
00666 operator[](myCells[myCellI]) = fillType;
00667
00668 nChanged++;
00669 }
00670 }
00671 }
00672 }
00673 return nChanged;
00674 }
00675
00676
00677
00678
00679
00680
00681 Foam::label Foam::cellClassification::fillHangingCells
00682 (
00683 const label meshType,
00684 const label fillType,
00685 const label maxIter
00686 )
00687 {
00688 label nTotChanged = 0;
00689
00690 for(label iter = 0; iter < maxIter; iter++)
00691 {
00692 label nChanged = 0;
00693
00694
00695 List<pointStatus> pointSide(mesh_.nPoints());
00696 classifyPoints(meshType, *this, pointSide);
00697
00698
00699
00700
00701 forAll(pointSide, pointI)
00702 {
00703 if (pointSide[pointI] == MIXED)
00704 {
00705 const labelList& pCells = mesh_.pointCells()[pointI];
00706
00707 forAll(pCells, i)
00708 {
00709 label cellI = pCells[i];
00710
00711 if (operator[](cellI) == meshType)
00712 {
00713 if (usesMixedPointsOnly(pointSide, cellI))
00714 {
00715 operator[](cellI) = fillType;
00716
00717 nChanged++;
00718 }
00719 }
00720 }
00721 }
00722 }
00723 nTotChanged += nChanged;
00724
00725 Pout<< "removeHangingCells : changed " << nChanged
00726 << " hanging cells" << endl;
00727
00728 if (nChanged == 0)
00729 {
00730 break;
00731 }
00732 }
00733
00734 return nTotChanged;
00735 }
00736
00737
00738 Foam::label Foam::cellClassification::fillRegionEdges
00739 (
00740 const label meshType,
00741 const label fillType,
00742 const label maxIter
00743 )
00744 {
00745 label nTotChanged = 0;
00746
00747 for(label iter = 0; iter < maxIter; iter++)
00748 {
00749
00750
00751 faceList outsideFaces;
00752 labelList outsideOwner;
00753
00754 getMeshOutside(meshType, outsideFaces, outsideOwner);
00755
00756
00757 primitiveFacePatch fp(outsideFaces, mesh_.points());
00758
00759 const labelListList& edgeFaces = fp.edgeFaces();
00760
00761 label nChanged = 0;
00762
00763
00764
00765 forAll(edgeFaces, edgeI)
00766 {
00767 const labelList& eFaces = edgeFaces[edgeI];
00768
00769 if (eFaces.size() > 2)
00770 {
00771
00772
00773
00774 forAll(eFaces, i)
00775 {
00776 label patchFaceI = eFaces[i];
00777
00778 label ownerCell = outsideOwner[patchFaceI];
00779
00780 if (operator[](ownerCell) == meshType)
00781 {
00782 operator[](ownerCell) = fillType;
00783
00784 nChanged++;
00785
00786 break;
00787 }
00788 }
00789 }
00790 }
00791
00792 nTotChanged += nChanged;
00793
00794 Pout<< "fillRegionEdges : changed " << nChanged
00795 << " cells using multiply connected edges" << endl;
00796
00797 if (nChanged == 0)
00798 {
00799 break;
00800 }
00801 }
00802
00803 return nTotChanged;
00804 }
00805
00806
00807 Foam::label Foam::cellClassification::fillRegionPoints
00808 (
00809 const label meshType,
00810 const label fillType,
00811 const label maxIter
00812 )
00813 {
00814 label nTotChanged = 0;
00815
00816 for(label iter = 0; iter < maxIter; iter++)
00817 {
00818
00819
00820 faceList outsideFaces;
00821 labelList outsideOwner;
00822
00823 getMeshOutside(meshType, outsideFaces, outsideOwner);
00824
00825
00826 primitiveFacePatch fp(outsideFaces, mesh_.points());
00827
00828 labelHashSet nonManifoldPoints;
00829
00830
00831 fp.checkPointManifold(false, &nonManifoldPoints);
00832
00833 const Map<label>& meshPointMap = fp.meshPointMap();
00834
00835 label nChanged = 0;
00836
00837 for
00838 (
00839 labelHashSet::const_iterator iter = nonManifoldPoints.begin();
00840 iter != nonManifoldPoints.end();
00841 ++iter
00842 )
00843 {
00844
00845 label patchPointI = meshPointMap[iter.key()];
00846
00847 const labelList& pFaces = fp.pointFaces()[patchPointI];
00848
00849
00850
00851
00852 forAll(pFaces, i)
00853 {
00854 label patchFaceI = pFaces[i];
00855
00856 label ownerCell = outsideOwner[patchFaceI];
00857
00858 if (operator[](ownerCell) == meshType)
00859 {
00860 operator[](ownerCell) = fillType;
00861
00862 nChanged++;
00863
00864 break;
00865 }
00866 }
00867 }
00868
00869 nTotChanged += nChanged;
00870
00871 Pout<< "fillRegionPoints : changed " << nChanged
00872 << " cells using multiply connected points" << endl;
00873
00874 if (nChanged == 0)
00875 {
00876 break;
00877 }
00878 }
00879
00880 return nTotChanged;
00881 }
00882
00883
00884 void Foam::cellClassification::writeStats(Ostream& os) const
00885 {
00886 os << "Cells:" << size() << endl
00887 << " notset : " << count(*this, NOTSET) << endl
00888 << " cut : " << count(*this, CUT) << endl
00889 << " inside : " << count(*this, INSIDE) << endl
00890 << " outside : " << count(*this, OUTSIDE) << endl;
00891 }
00892
00893
00894
00895
00896 void Foam::cellClassification::operator=(const Foam::cellClassification& rhs)
00897 {
00898 labelList::operator=(rhs);
00899 }
00900
00901
00902