FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

isoSurfaceCell.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "isoSurfaceCell.H"
00027 #include <OpenFOAM/dictionary.H>
00028 #include <OpenFOAM/polyMesh.H>
00029 #include <OpenFOAM/mergePoints.H>
00030 #include <OpenFOAM/tetMatcher.H>
00031 #include <OpenFOAM/syncTools.H>
00032 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00033 #include <triSurface/faceTriangulation.H>
00034 
00035 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00036 
00037 namespace Foam
00038 {
00039     defineTypeNameAndDebug(isoSurfaceCell, 0);
00040 }
00041 
00042 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00043 
00044 Foam::scalar Foam::isoSurfaceCell::isoFraction
00045 (
00046     const scalar s0,
00047     const scalar s1
00048 ) const
00049 {
00050     scalar d = s1-s0;
00051 
00052     if (mag(d) > VSMALL)
00053     {
00054         return (iso_-s0)/d;
00055     }
00056     else
00057     {
00058         return -1.0;
00059     }
00060 }
00061 
00062 
00063 //Foam::List<Foam::triFace> Foam::isoSurfaceCell::triangulate(const face& f)
00064 // const
00065 //{
00066 //    faceList triFaces(f.nTriangles(mesh_.points()));
00067 //    label triFaceI = 0;
00068 //    f.triangles(mesh_.points(), triFaceI, triFaces);
00069 //
00070 //    List<triFace> tris(triFaces.size());
00071 //    forAll(triFaces, i)
00072 //    {
00073 //        tris[i][0] = triFaces[i][0];
00074 //        tris[i][1] = triFaces[i][1];
00075 //        tris[i][2] = triFaces[i][2];
00076 //    }
00077 //    return tris;
00078 //}
00079 
00080 
00081 Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
00082 (
00083     const PackedBoolList& isTet,
00084     const scalarField& cellValues,
00085     const scalarField& pointValues,
00086     const label cellI
00087 ) const
00088 {
00089     const cell& cFaces = mesh_.cells()[cellI];
00090 
00091     if (isTet.get(cellI) == 1)
00092     {
00093         forAll(cFaces, cFaceI)
00094         {
00095             const face& f = mesh_.faces()[cFaces[cFaceI]];
00096 
00097             for (label fp = 1; fp < f.size() - 1; fp++)
00098             {
00099                 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
00100 
00101                 bool aLower = (pointValues[tri[0]] < iso_);
00102                 bool bLower = (pointValues[tri[1]] < iso_);
00103                 bool cLower = (pointValues[tri[2]] < iso_);
00104 
00105                 if (aLower == bLower && aLower == cLower)
00106                 {}
00107                 else
00108                 {
00109                     return CUT;
00110                 }
00111             }
00112         }
00113         return NOTCUT;
00114     }
00115     else
00116     {
00117         bool cellLower = (cellValues[cellI] < iso_);
00118 
00119         // First check if there is any cut in cell
00120         bool edgeCut = false;
00121 
00122         forAll(cFaces, cFaceI)
00123         {
00124             const face& f = mesh_.faces()[cFaces[cFaceI]];
00125 
00126             // Check pyramids cut
00127             forAll(f, fp)
00128             {
00129                 if ((pointValues[f[fp]] < iso_) != cellLower)
00130                 {
00131                     edgeCut = true;
00132                     break;
00133                 }
00134             }
00135 
00136             if (edgeCut)
00137             {
00138                 break;
00139             }
00140 
00141             for (label fp = 1; fp < f.size() - 1; fp++)
00142             {
00143                 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
00144             //List<triFace> tris(triangulate(f));
00145             //forAll(tris, i)
00146             //{
00147             //    const triFace& tri = tris[i];
00148 
00149                 bool aLower = (pointValues[tri[0]] < iso_);
00150                 bool bLower = (pointValues[tri[1]] < iso_);
00151                 bool cLower = (pointValues[tri[2]] < iso_);
00152 
00153                 if (aLower == bLower && aLower == cLower)
00154                 {}
00155                 else
00156                 {
00157                     edgeCut = true;
00158                     break;
00159                 }
00160             }
00161 
00162             if (edgeCut)
00163             {
00164                 break;
00165             }
00166         }
00167 
00168         if (edgeCut)
00169         {
00170             // Count actual cuts (expensive since addressing needed)
00171             // Note: not needed if you don't want to preserve maxima/minima
00172             // centred around cellcentre. In that case just always return CUT
00173 
00174             const labelList& cPoints = mesh_.cellPoints(cellI);
00175 
00176             label nPyrCuts = 0;
00177 
00178             forAll(cPoints, i)
00179             {
00180                 if ((pointValues[cPoints[i]] < iso_) != cellLower)
00181                 {
00182                     nPyrCuts++;
00183                 }
00184             }
00185 
00186             if (nPyrCuts == cPoints.size())
00187             {
00188                 return SPHERE;
00189             }
00190             else
00191             {
00192                 return CUT;
00193             }
00194         }
00195         else
00196         {
00197             return NOTCUT;
00198         }
00199     }
00200 }
00201 
00202 
00203 void Foam::isoSurfaceCell::calcCutTypes
00204 (
00205     const PackedBoolList& isTet,
00206     const scalarField& cVals,
00207     const scalarField& pVals
00208 )
00209 {
00210     cellCutType_.setSize(mesh_.nCells());
00211     nCutCells_ = 0;
00212     forAll(mesh_.cells(), cellI)
00213     {
00214         cellCutType_[cellI] = calcCutType(isTet, cVals, pVals, cellI);
00215 
00216         if (cellCutType_[cellI] == CUT)
00217         {
00218             nCutCells_++;
00219         }
00220     }
00221 
00222     if (debug)
00223     {
00224         Pout<< "isoSurfaceCell : detected " << nCutCells_
00225             << " candidate cut cells." << endl;
00226     }
00227 }
00228 
00229 
00230 
00231 // Return the two common points between two triangles
00232 Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
00233 (
00234     const labelledTri& tri0,
00235     const labelledTri& tri1
00236 )
00237 {
00238     labelPair common(-1, -1);
00239 
00240     label fp0 = 0;
00241     label fp1 = findIndex(tri1, tri0[fp0]);
00242 
00243     if (fp1 == -1)
00244     {
00245         fp0 = 1;
00246         fp1 = findIndex(tri1, tri0[fp0]);
00247     }
00248 
00249     if (fp1 != -1)
00250     {
00251         // So tri0[fp0] is tri1[fp1]
00252 
00253         // Find next common point
00254         label fp0p1 = tri0.fcIndex(fp0);
00255         label fp1p1 = tri1.fcIndex(fp1);
00256         label fp1m1 = tri1.rcIndex(fp1);
00257 
00258         if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
00259         {
00260             common[0] = tri0[fp0];
00261             common[1] = tri0[fp0p1];
00262         }
00263     }
00264     return common;
00265 }
00266 
00267 
00268 // Caculate centre of surface.
00269 Foam::point Foam::isoSurfaceCell::calcCentre(const triSurface& s)
00270 {
00271     vector sum = vector::zero;
00272 
00273     forAll(s, i)
00274     {
00275         sum += s[i].centre(s.points());
00276     }
00277     return sum/s.size();
00278 }
00279 
00280 
00281 // Replace surface (localPoints, localTris) with single point. Returns
00282 // point. Destructs arguments.
00283 Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
00284 (
00285     pointField& localPoints,
00286     DynamicList<labelledTri, 64>& localTris
00287 )
00288 {
00289     pointIndexHit info(false, vector::zero, localTris.size());
00290 
00291     if (localTris.size() == 1)
00292     {
00293         const labelledTri& tri = localTris[0];
00294         info.setPoint(tri.centre(localPoints));
00295         info.setHit();
00296     }
00297     else if (localTris.size() == 2)
00298     {
00299         // Check if the two triangles share an edge.
00300         const labelledTri& tri0 = localTris[0];
00301         const labelledTri& tri1 = localTris[0];
00302 
00303         labelPair shared = findCommonPoints(tri0, tri1);
00304 
00305         if (shared[0] != -1)
00306         {
00307             info.setPoint
00308             (
00309                 0.5
00310               * (
00311                     tri0.centre(localPoints)
00312                   + tri1.centre(localPoints)
00313                 )
00314             );
00315             info.setHit();
00316         }
00317     }
00318     else if (localTris.size())
00319     {
00320         // Check if single region. Rare situation.
00321         triSurface surf
00322         (
00323             localTris,
00324             geometricSurfacePatchList(0),
00325             localPoints,
00326             true
00327         );
00328         localTris.clearStorage();
00329 
00330         labelList faceZone;
00331         label nZones = surf.markZones
00332         (
00333             boolList(surf.nEdges(), false),
00334             faceZone
00335         );
00336 
00337         if (nZones == 1)
00338         {
00339             info.setPoint(calcCentre(surf));
00340             info.setHit();
00341         }
00342     }
00343 
00344     return info;
00345 }
00346 
00347 
00348 void Foam::isoSurfaceCell::calcSnappedCc
00349 (
00350     const PackedBoolList& isTet,
00351     const scalarField& cVals,
00352     const scalarField& pVals,
00353 
00354     DynamicList<point>& snappedPoints,
00355     labelList& snappedCc
00356 ) const
00357 {
00358     const pointField& cc = mesh_.cellCentres();
00359     const pointField& pts = mesh_.points();
00360 
00361     snappedCc.setSize(mesh_.nCells());
00362     snappedCc = -1;
00363 
00364     // Work arrays
00365     DynamicList<point, 64> localPoints(64);
00366     DynamicList<labelledTri, 64> localTris(64);
00367     Map<label> pointToLocal(64);
00368 
00369     forAll(mesh_.cells(), cellI)
00370     {
00371         if (cellCutType_[cellI] == CUT && isTet.get(cellI) == 0)
00372         {
00373             scalar cVal = cVals[cellI];
00374 
00375             const cell& cFaces = mesh_.cells()[cellI];
00376 
00377             localPoints.clear();
00378             localTris.clear();
00379             pointToLocal.clear();
00380 
00381             // Create points for all intersections close to cell centre
00382             // (i.e. from pyramid edges)
00383 
00384             forAll(cFaces, cFaceI)
00385             {
00386                 const face& f = mesh_.faces()[cFaces[cFaceI]];
00387 
00388                 forAll(f, fp)
00389                 {
00390                     label pointI = f[fp];
00391 
00392                     scalar s = isoFraction(cVal, pVals[pointI]);
00393 
00394                     if (s >= 0.0 && s <= 0.5)
00395                     {
00396                         if (pointToLocal.insert(pointI, localPoints.size()))
00397                         {
00398                             localPoints.append((1.0-s)*cc[cellI]+s*pts[pointI]);
00399                         }
00400                     }
00401                 }
00402             }
00403 
00404             if (localPoints.size() == 1)
00405             {
00406                 // No need for any analysis.
00407                 snappedCc[cellI] = snappedPoints.size();
00408                 snappedPoints.append(localPoints[0]);
00409 
00410                 //Pout<< "cell:" << cellI
00411                 //    << " at " << mesh_.cellCentres()[cellI]
00412                 //    << " collapsing " << localPoints
00413                 //    << " intersections down to "
00414                 //    << snappedPoints[snappedCc[cellI]] << endl;
00415             }
00416             else if (localPoints.size() == 2)
00417             {
00418                 //? No need for any analysis.???
00419                 snappedCc[cellI] = snappedPoints.size();
00420                 snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
00421 
00422                 //Pout<< "cell:" << cellI
00423                 //    << " at " << mesh_.cellCentres()[cellI]
00424                 //    << " collapsing " << localPoints
00425                 //    << " intersections down to "
00426                 //    << snappedPoints[snappedCc[cellI]] << endl;
00427             }
00428             else if (localPoints.size())
00429             {
00430                 // Need to analyse
00431                 forAll(cFaces, cFaceI)
00432                 {
00433                     const face& f = mesh_.faces()[cFaces[cFaceI]];
00434 
00435                     // Do a tetrahedrisation. Each face to cc becomes pyr.
00436                     // Each pyr gets split into tets by diagonalisation
00437                     // of face.
00438 
00439                     for (label fp = 1; fp < f.size() - 1; fp++)
00440                     {
00441                         triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
00442                     //List<triFace> tris(triangulate(f));
00443                     //forAll(tris, i)
00444                     //{
00445                     //    const triFace& tri = tris[i];
00446 
00447                         // Get fractions for the three edges to cell centre
00448                         FixedList<scalar, 3> s(3);
00449                         s[0] = isoFraction(cVal, pVals[tri[0]]);
00450                         s[1] = isoFraction(cVal, pVals[tri[1]]);
00451                         s[2] = isoFraction(cVal, pVals[tri[2]]);
00452 
00453                         if
00454                         (
00455                             (s[0] >= 0.0 && s[0] <= 0.5)
00456                          && (s[1] >= 0.0 && s[1] <= 0.5)
00457                          && (s[2] >= 0.0 && s[2] <= 0.5)
00458                         )
00459                         {
00460                             localTris.append
00461                             (
00462                                 labelledTri
00463                                 (
00464                                     pointToLocal[tri[0]],
00465                                     pointToLocal[tri[1]],
00466                                     pointToLocal[tri[2]],
00467                                     0
00468                                 )
00469                             );
00470                         }
00471                     }
00472                 }
00473 
00474                 pointField surfPoints;
00475                 surfPoints.transfer(localPoints);
00476                 pointIndexHit info = collapseSurface(surfPoints, localTris);
00477 
00478                 if (info.hit())
00479                 {
00480                     snappedCc[cellI] = snappedPoints.size();
00481                     snappedPoints.append(info.hitPoint());
00482 
00483                     //Pout<< "cell:" << cellI
00484                     //    << " at " << mesh_.cellCentres()[cellI]
00485                     //    << " collapsing " << surfPoints
00486                     //    << " intersections down to "
00487                     //    << snappedPoints[snappedCc[cellI]] << endl;
00488                 }
00489             }
00490         }
00491     }
00492 }
00493 
00494 
00495 // Generate triangles for face connected to pointI
00496 void Foam::isoSurfaceCell::genPointTris
00497 (
00498     const scalarField& cellValues,
00499     const scalarField& pointValues,
00500     const label pointI,
00501     const face& f,
00502     const label cellI,
00503     DynamicList<point, 64>& localTriPoints
00504 ) const
00505 {
00506     const pointField& cc = mesh_.cellCentres();
00507     const pointField& pts = mesh_.points();
00508 
00509     for (label fp = 1; fp < f.size() - 1; fp++)
00510     {
00511         triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
00512     //List<triFace> tris(triangulate(f));
00513     //forAll(tris, i)
00514     //{
00515     //    const triFace& tri = tris[i];
00516 
00517         label index = findIndex(tri, pointI);
00518 
00519         if (index == -1)
00520         {
00521             continue;
00522         }
00523 
00524         // Tet between index..index-1, index..index+1, index..cc
00525         label b = tri[tri.fcIndex(index)];
00526         label c = tri[tri.rcIndex(index)];
00527 
00528         // Get fractions for the three edges emanating from point
00529         FixedList<scalar, 3> s(3);
00530         s[0] = isoFraction(pointValues[pointI], pointValues[b]);
00531         s[1] = isoFraction(pointValues[pointI], pointValues[c]);
00532         s[2] = isoFraction(pointValues[pointI], cellValues[cellI]);
00533 
00534         if
00535         (
00536             (s[0] >= 0.0 && s[0] <= 0.5)
00537          && (s[1] >= 0.0 && s[1] <= 0.5)
00538          && (s[2] >= 0.0 && s[2] <= 0.5)
00539         )
00540         {
00541             localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[b]);
00542             localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[c]);
00543             localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*cc[cellI]);
00544         }
00545     }
00546 }
00547 
00548 
00549 // Generate triangle for tet connected to pointI
00550 void Foam::isoSurfaceCell::genPointTris
00551 (
00552     const scalarField& pointValues,
00553     const label pointI,
00554     const label cellI,
00555     DynamicList<point, 64>& localTriPoints
00556 ) const
00557 {
00558     const pointField& pts = mesh_.points();
00559     const cell& cFaces = mesh_.cells()[cellI];
00560 
00561     FixedList<label, 4> tet;
00562 
00563     label face0 = cFaces[0];
00564     const face& f0 = mesh_.faces()[face0];
00565 
00566     if (mesh_.faceOwner()[face0] != cellI)
00567     {
00568         tet[0] = f0[0];
00569         tet[1] = f0[1];
00570         tet[2] = f0[2];
00571     }
00572     else
00573     {
00574         tet[0] = f0[0];
00575         tet[1] = f0[2];
00576         tet[2] = f0[1];
00577     }
00578 
00579     // Find the point on the next face that is not on f0
00580     const face& f1 = mesh_.faces()[cFaces[1]];
00581 
00582     forAll(f1, fp)
00583     {
00584         label p1 = f1[fp];
00585 
00586         if (p1 != tet[0] && p1 != tet[1] && p1 != tet[2])
00587         {
00588             tet[3] = p1;
00589             break;
00590         }
00591     }
00592 
00593 
00594     // Get the index of pointI
00595 
00596     label i0 = findIndex(tet, pointI);
00597     label i1 = tet.fcIndex(i0);
00598     label i2 = tet.fcIndex(i1);
00599     label i3 = tet.fcIndex(i2);
00600 
00601     // Get fractions for the three edges emanating from point
00602     FixedList<scalar, 3> s(3);
00603     s[0] = isoFraction(pointValues[pointI], pointValues[tet[i1]]);
00604     s[1] = isoFraction(pointValues[pointI], pointValues[tet[i2]]);
00605     s[2] = isoFraction(pointValues[pointI], pointValues[tet[i3]]);
00606 
00607     if
00608     (
00609         (s[0] >= 0.0 && s[0] <= 0.5)
00610      && (s[1] >= 0.0 && s[1] <= 0.5)
00611      && (s[2] >= 0.0 && s[2] <= 0.5)
00612     )
00613     {
00614         localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[tet[i1]]);
00615         localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[tet[i2]]);
00616         localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*pts[tet[i3]]);
00617     }
00618 }
00619 
00620 
00621 void Foam::isoSurfaceCell::calcSnappedPoint
00622 (
00623     const PackedBoolList& isBoundaryPoint,
00624     const PackedBoolList& isTet,
00625     const scalarField& cVals,
00626     const scalarField& pVals,
00627 
00628     DynamicList<point>& snappedPoints,
00629     labelList& snappedPoint
00630 ) const
00631 {
00632     const point greatPoint(VGREAT, VGREAT, VGREAT);
00633     pointField collapsedPoint(mesh_.nPoints(), greatPoint);
00634 
00635 
00636     // Work arrays
00637     DynamicList<point, 64> localTriPoints(100);
00638     labelHashSet localPointCells(100);
00639 
00640     forAll(mesh_.pointFaces(), pointI)
00641     {
00642         if (isBoundaryPoint.get(pointI) == 1)
00643         {
00644             continue;
00645         }
00646 
00647         const labelList& pFaces = mesh_.pointFaces()[pointI];
00648 
00649         bool anyCut = false;
00650 
00651         forAll(pFaces, i)
00652         {
00653             label faceI = pFaces[i];
00654 
00655             if
00656             (
00657                 cellCutType_[mesh_.faceOwner()[faceI]] == CUT
00658              || (
00659                     mesh_.isInternalFace(faceI)
00660                  && cellCutType_[mesh_.faceNeighbour()[faceI]] == CUT
00661                 )
00662             )
00663             {
00664                 anyCut = true;
00665                 break;
00666             }
00667         }
00668 
00669         if (!anyCut)
00670         {
00671             continue;
00672         }
00673 
00674 
00675         localPointCells.clear();
00676         localTriPoints.clear();
00677 
00678         forAll(pFaces, pFaceI)
00679         {
00680             label faceI = pFaces[pFaceI];
00681             const face& f = mesh_.faces()[faceI];
00682             label own = mesh_.faceOwner()[faceI];
00683 
00684             // Triangulate around f[0] on owner side
00685             if (isTet.get(own) == 1)
00686             {
00687                 if (localPointCells.insert(own))
00688                 {
00689                     genPointTris(pVals, pointI, own, localTriPoints);
00690                 }
00691             }
00692             else
00693             {
00694                 genPointTris(cVals, pVals, pointI, f, own, localTriPoints);
00695             }
00696 
00697             if (mesh_.isInternalFace(faceI))
00698             {
00699                 label nei = mesh_.faceNeighbour()[faceI];
00700 
00701                 if (isTet.get(nei) == 1)
00702                 {
00703                     if (localPointCells.insert(nei))
00704                     {
00705                         genPointTris(pVals, pointI, nei, localTriPoints);
00706                     }
00707                 }
00708                 else
00709                 {
00710                     genPointTris(cVals, pVals, pointI, f, nei, localTriPoints);
00711                 }
00712             }
00713         }
00714 
00715         if (localTriPoints.size() == 3)
00716         {
00717             // Single triangle. No need for any analysis. Average points.
00718             pointField points;
00719             points.transfer(localTriPoints);
00720             collapsedPoint[pointI] = sum(points)/points.size();
00721 
00722             //Pout<< "    point:" << pointI
00723             //    << " replacing coord:" << mesh_.points()[pointI]
00724             //    << " by average:" << collapsedPoint[pointI] << endl;
00725         }
00726         else if (localTriPoints.size())
00727         {
00728             // Convert points into triSurface.
00729 
00730             // Merge points and compact out non-valid triangles
00731             labelList triMap;               // merged to unmerged triangle
00732             labelList triPointReverseMap;   // unmerged to merged point
00733             triSurface surf
00734             (
00735                 stitchTriPoints
00736                 (
00737                     false,                  // do not check for duplicate tris
00738                     localTriPoints,
00739                     triPointReverseMap,
00740                     triMap
00741                 )
00742             );
00743 
00744             labelList faceZone;
00745             label nZones = surf.markZones
00746             (
00747                 boolList(surf.nEdges(), false),
00748                 faceZone
00749             );
00750 
00751             if (nZones == 1)
00752             {
00753                 collapsedPoint[pointI] = calcCentre(surf);
00754                 //Pout<< "    point:" << pointI << " nZones:" << nZones
00755                 //    << " replacing coord:" << mesh_.points()[pointI]
00756                 //    << " by average:" << collapsedPoint[pointI] << endl;
00757             }
00758         }
00759     }
00760 
00761     syncTools::syncPointList
00762     (
00763         mesh_,
00764         collapsedPoint,
00765         minEqOp<point>(),
00766         greatPoint,
00767         true                // are coordinates so separate
00768     );
00769 
00770     snappedPoint.setSize(mesh_.nPoints());
00771     snappedPoint = -1;
00772 
00773     forAll(collapsedPoint, pointI)
00774     {
00775         if (collapsedPoint[pointI] != greatPoint)
00776         {
00777             snappedPoint[pointI] = snappedPoints.size();
00778             snappedPoints.append(collapsedPoint[pointI]);
00779         }
00780     }
00781 }
00782 
00783 
00784 
00785 
00786 Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
00787 (
00788     const bool checkDuplicates,
00789     const List<point>& triPoints,
00790     labelList& triPointReverseMap,  // unmerged to merged point
00791     labelList& triMap               // merged to unmerged triangle
00792 ) const
00793 {
00794     label nTris = triPoints.size()/3;
00795 
00796     if ((triPoints.size() % 3) != 0)
00797     {
00798         FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
00799             << "Problem: number of points " << triPoints.size()
00800             << " not a multiple of 3." << abort(FatalError);
00801     }
00802 
00803     pointField newPoints;
00804     mergePoints
00805     (
00806         triPoints,
00807         mergeDistance_,
00808         false,
00809         triPointReverseMap,
00810         newPoints
00811     );
00812 
00813     // Check that enough merged.
00814     if (debug)
00815     {
00816         pointField newNewPoints;
00817         labelList oldToNew;
00818         bool hasMerged = mergePoints
00819         (
00820             newPoints,
00821             mergeDistance_,
00822             true,
00823             oldToNew,
00824             newNewPoints
00825         );
00826 
00827         if (hasMerged)
00828         {
00829             FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
00830                 << "Merged points contain duplicates"
00831                 << " when merging with distance " << mergeDistance_ << endl
00832                 << "merged:" << newPoints.size() << " re-merged:"
00833                 << newNewPoints.size()
00834                 << abort(FatalError);
00835         }
00836     }
00837 
00838 
00839     List<labelledTri> tris;
00840     {
00841         DynamicList<labelledTri> dynTris(nTris);
00842         label rawPointI = 0;
00843         DynamicList<label> newToOldTri(nTris);
00844 
00845         for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
00846         {
00847             labelledTri tri
00848             (
00849                 triPointReverseMap[rawPointI],
00850                 triPointReverseMap[rawPointI+1],
00851                 triPointReverseMap[rawPointI+2],
00852                 0
00853             );
00854             rawPointI += 3;
00855 
00856             if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
00857             {
00858                 newToOldTri.append(oldTriI);
00859                 dynTris.append(tri);
00860             }
00861         }
00862 
00863         triMap.transfer(newToOldTri);
00864         tris.transfer(dynTris);
00865     }
00866 
00867 
00868     // Use face centres to determine 'flat hole' situation (see RMT paper).
00869     // Two unconnected triangles get connected because (some of) the edges
00870     // separating them get collapsed. Below only checks for duplicate triangles,
00871     // not non-manifold edge connectivity.
00872     if (checkDuplicates)
00873     {
00874         if (debug)
00875         {
00876             Pout<< "isoSurfaceCell : merged from " << nTris
00877                 << " down to " << tris.size() << " triangles." << endl;
00878         }
00879 
00880         pointField centres(tris.size());
00881         forAll(tris, triI)
00882         {
00883             centres[triI] = tris[triI].centre(newPoints);
00884         }
00885 
00886         pointField mergedCentres;
00887         labelList oldToMerged;
00888         bool hasMerged = mergePoints
00889         (
00890             centres,
00891             mergeDistance_,
00892             false,
00893             oldToMerged,
00894             mergedCentres
00895         );
00896 
00897         if (debug)
00898         {
00899             Pout<< "isoSurfaceCell : detected "
00900                 << centres.size()-mergedCentres.size()
00901                 << " duplicate triangles." << endl;
00902         }
00903 
00904         if (hasMerged)
00905         {
00906             // Filter out duplicates.
00907             label newTriI = 0;
00908             DynamicList<label> newToOldTri(tris.size());
00909             labelList newToMaster(mergedCentres.size(), -1);
00910             forAll(tris, triI)
00911             {
00912                 label mergedI = oldToMerged[triI];
00913 
00914                 if (newToMaster[mergedI] == -1)
00915                 {
00916                     newToMaster[mergedI] = triI;
00917                     newToOldTri.append(triMap[triI]);
00918                     tris[newTriI++] = tris[triI];
00919                 }
00920             }
00921 
00922             triMap.transfer(newToOldTri);
00923             tris.setSize(newTriI);
00924         }
00925     }
00926 
00927     return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
00928 }
00929 
00930 
00931 // Does face use valid vertices?
00932 bool Foam::isoSurfaceCell::validTri(const triSurface& surf, const label faceI)
00933 {
00934     // Simple check on indices ok.
00935 
00936     const labelledTri& f = surf[faceI];
00937 
00938     if
00939     (
00940         (f[0] < 0) || (f[0] >= surf.points().size())
00941      || (f[1] < 0) || (f[1] >= surf.points().size())
00942      || (f[2] < 0) || (f[2] >= surf.points().size())
00943     )
00944     {
00945         WarningIn("validTri(const triSurface&, const label)")
00946             << "triangle " << faceI << " vertices " << f
00947             << " uses point indices outside point range 0.."
00948             << surf.points().size()-1 << endl;
00949 
00950         return false;
00951     }
00952 
00953     if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
00954     {
00955         WarningIn("validTri(const triSurface&, const label)")
00956             << "triangle " << faceI
00957             << " uses non-unique vertices " << f
00958             << endl;
00959         return false;
00960     }
00961 
00962     // duplicate triangle check
00963 
00964     const labelList& fFaces = surf.faceFaces()[faceI];
00965 
00966     // Check if faceNeighbours use same points as this face.
00967     // Note: discards normal information - sides of baffle are merged.
00968     forAll(fFaces, i)
00969     {
00970         label nbrFaceI = fFaces[i];
00971 
00972         if (nbrFaceI <= faceI)
00973         {
00974             // lower numbered faces already checked
00975             continue;
00976         }
00977 
00978         const labelledTri& nbrF = surf[nbrFaceI];
00979 
00980         if
00981         (
00982             ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
00983          && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
00984          && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
00985         )
00986         {
00987             WarningIn("validTri(const triSurface&, const label)")
00988                 << "triangle " << faceI << " vertices " << f
00989                 << " has the same vertices as triangle " << nbrFaceI
00990                 << " vertices " << nbrF
00991                 << endl;
00992 
00993             return false;
00994         }
00995     }
00996     return true;
00997 }
00998 
00999 
01000 void Foam::isoSurfaceCell::calcAddressing
01001 (
01002     const triSurface& surf,
01003     List<FixedList<label, 3> >& faceEdges,
01004     labelList& edgeFace0,
01005     labelList& edgeFace1,
01006     Map<labelList>& edgeFacesRest
01007 ) const
01008 {
01009     const pointField& points = surf.points();
01010 
01011     pointField edgeCentres(3*surf.size());
01012     label edgeI = 0;
01013     forAll(surf, triI)
01014     {
01015         const labelledTri& tri = surf[triI];
01016         edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
01017         edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
01018         edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
01019     }
01020 
01021     pointField mergedCentres;
01022     labelList oldToMerged;
01023     bool hasMerged = mergePoints
01024     (
01025         edgeCentres,
01026         mergeDistance_,
01027         false,
01028         oldToMerged,
01029         mergedCentres
01030     );
01031 
01032     if (debug)
01033     {
01034         Pout<< "isoSurfaceCell : detected "
01035             << mergedCentres.size()
01036             << " edges on " << surf.size() << " triangles." << endl;
01037     }
01038 
01039     if (!hasMerged)
01040     {
01041         if (surf.size() == 1)
01042         {
01043             faceEdges.setSize(1);
01044             faceEdges[0][0] = 0;
01045             faceEdges[0][1] = 1;
01046             faceEdges[0][2] = 2;
01047             edgeFace0.setSize(1);
01048             edgeFace0[0] = 0;
01049             edgeFace1.setSize(1);
01050             edgeFace1[0] = -1;
01051             edgeFacesRest.clear();
01052         }
01053         return;
01054     }
01055 
01056 
01057     // Determine faceEdges
01058     faceEdges.setSize(surf.size());
01059     edgeI = 0;
01060     forAll(surf, triI)
01061     {
01062         faceEdges[triI][0] = oldToMerged[edgeI++];
01063         faceEdges[triI][1] = oldToMerged[edgeI++];
01064         faceEdges[triI][2] = oldToMerged[edgeI++];
01065     }
01066 
01067 
01068     // Determine edgeFaces
01069     edgeFace0.setSize(mergedCentres.size());
01070     edgeFace0 = -1;
01071     edgeFace1.setSize(mergedCentres.size());
01072     edgeFace1 = -1;
01073     edgeFacesRest.clear();
01074 
01075     forAll(oldToMerged, oldEdgeI)
01076     {
01077         label triI = oldEdgeI / 3;
01078         label edgeI = oldToMerged[oldEdgeI];
01079 
01080         if (edgeFace0[edgeI] == -1)
01081         {
01082             edgeFace0[edgeI] = triI;
01083         }
01084         else if (edgeFace1[edgeI] == -1)
01085         {
01086             edgeFace1[edgeI] = triI;
01087         }
01088         else
01089         {
01090             //WarningIn("orientSurface(triSurface&)")
01091             //    << "Edge " << edgeI << " with centre " << mergedCentres[edgeI]
01092             //    << " used by more than two triangles: " << edgeFace0[edgeI]
01093             //    << ", "
01094             //    << edgeFace1[edgeI] << " and " << triI << endl;
01095             Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
01096 
01097             if (iter != edgeFacesRest.end())
01098             {
01099                 labelList& eFaces = iter();
01100                 label sz = eFaces.size();
01101                 eFaces.setSize(sz+1);
01102                 eFaces[sz] = triI;
01103             }
01104             else
01105             {
01106                 edgeFacesRest.insert(edgeI, labelList(1, triI));
01107             }
01108         }
01109     }
01110 }
01111 
01112 
01113 void Foam::isoSurfaceCell::walkOrientation
01114 (
01115     const triSurface& surf,
01116     const List<FixedList<label, 3> >& faceEdges,
01117     const labelList& edgeFace0,
01118     const labelList& edgeFace1,
01119     const label seedTriI,
01120     labelList& flipState
01121 )
01122 {
01123     // Do walk for consistent orientation.
01124     DynamicList<label> changedFaces(surf.size());
01125 
01126     changedFaces.append(seedTriI);
01127 
01128     while (changedFaces.size())
01129     {
01130         DynamicList<label> newChangedFaces(changedFaces.size());
01131 
01132         forAll(changedFaces, i)
01133         {
01134             label triI = changedFaces[i];
01135             const labelledTri& tri = surf[triI];
01136             const FixedList<label, 3>& fEdges = faceEdges[triI];
01137 
01138             forAll(fEdges, fp)
01139             {
01140                 label edgeI = fEdges[fp];
01141 
01142                 // my points:
01143                 label p0 = tri[fp];
01144                 label p1 = tri[tri.fcIndex(fp)];
01145 
01146                 label nbrI =
01147                 (
01148                     edgeFace0[edgeI] != triI
01149                   ? edgeFace0[edgeI]
01150                   : edgeFace1[edgeI]
01151                 );
01152 
01153                 if (nbrI != -1 && flipState[nbrI] == -1)
01154                 {
01155                     const labelledTri& nbrTri = surf[nbrI];
01156 
01157                     // nbr points
01158                     label nbrFp = findIndex(nbrTri, p0);
01159                     label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
01160 
01161                     bool sameOrientation = (p1 == nbrP1);
01162 
01163                     if (flipState[triI] == 0)
01164                     {
01165                         flipState[nbrI] = (sameOrientation ? 0 : 1);
01166                     }
01167                     else
01168                     {
01169                         flipState[nbrI] = (sameOrientation ? 1 : 0);
01170                     }
01171                     newChangedFaces.append(nbrI);
01172                 }
01173             }
01174         }
01175 
01176         changedFaces.transfer(newChangedFaces);
01177     }
01178 }
01179 
01180 
01181 void Foam::isoSurfaceCell::orientSurface
01182 (
01183     triSurface& surf,
01184     const List<FixedList<label, 3> >& faceEdges,
01185     const labelList& edgeFace0,
01186     const labelList& edgeFace1,
01187     const Map<labelList>& edgeFacesRest
01188 )
01189 {
01190     // -1 : unvisited
01191     //  0 : leave as is
01192     //  1 : flip
01193     labelList flipState(surf.size(), -1);
01194 
01195     label seedTriI = 0;
01196 
01197     while (true)
01198     {
01199         // Find first unvisited triangle
01200         for
01201         (
01202             ;
01203             seedTriI < surf.size() && flipState[seedTriI] != -1;
01204             seedTriI++
01205         )
01206         {}
01207 
01208         if (seedTriI == surf.size())
01209         {
01210             break;
01211         }
01212 
01213         // Note: Determine orientation of seedTriI?
01214         // for now assume it is ok
01215         flipState[seedTriI] = 0;
01216 
01217         walkOrientation
01218         (
01219             surf,
01220             faceEdges,
01221             edgeFace0,
01222             edgeFace1,
01223             seedTriI,
01224             flipState
01225         );
01226     }
01227 
01228     // Do actual flipping
01229     surf.clearOut();
01230     forAll(surf, triI)
01231     {
01232         if (flipState[triI] == 1)
01233         {
01234             labelledTri tri(surf[triI]);
01235 
01236             surf[triI][0] = tri[0];
01237             surf[triI][1] = tri[2];
01238             surf[triI][2] = tri[1];
01239         }
01240         else if (flipState[triI] == -1)
01241         {
01242             FatalErrorIn
01243             (
01244                 "isoSurfaceCell::orientSurface(triSurface&, const label)"
01245             )   << "problem" << abort(FatalError);
01246         }
01247     }
01248 }
01249 
01250 
01251 // Checks if triangle is connected through edgeI only.
01252 bool Foam::isoSurfaceCell::danglingTriangle
01253 (
01254     const FixedList<label, 3>& fEdges,
01255     const labelList& edgeFace1
01256 )
01257 {
01258     label nOpen = 0;
01259     forAll(fEdges, i)
01260     {
01261         if (edgeFace1[fEdges[i]] == -1)
01262         {
01263             nOpen++;
01264         }
01265     }
01266 
01267     if (nOpen == 1 || nOpen == 2 || nOpen == 3)
01268     {
01269         return true;
01270     }
01271     else
01272     {
01273         return false;
01274     }
01275 }
01276 
01277 
01278 // Mark triangles to keep. Returns number of dangling triangles.
01279 Foam::label Foam::isoSurfaceCell::markDanglingTriangles
01280 (
01281     const List<FixedList<label, 3> >& faceEdges,
01282     const labelList& edgeFace0,
01283     const labelList& edgeFace1,
01284     const Map<labelList>& edgeFacesRest,
01285     boolList& keepTriangles
01286 )
01287 {
01288     keepTriangles.setSize(faceEdges.size());
01289     keepTriangles = true;
01290 
01291     label nDangling = 0;
01292 
01293     // Remove any dangling triangles
01294     forAllConstIter(Map<labelList>,  edgeFacesRest, iter)
01295     {
01296         // These are all the non-manifold edges. Filter out all triangles
01297         // with only one connected edge (= this edge)
01298 
01299         label edgeI = iter.key();
01300         const labelList& otherEdgeFaces = iter();
01301 
01302         // Remove all dangling triangles
01303         if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
01304         {
01305             keepTriangles[edgeFace0[edgeI]] = false;
01306             nDangling++;
01307         }
01308         if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
01309         {
01310             keepTriangles[edgeFace1[edgeI]] = false;
01311             nDangling++;
01312         }
01313         forAll(otherEdgeFaces, i)
01314         {
01315             label triI = otherEdgeFaces[i];
01316             if (danglingTriangle(faceEdges[triI], edgeFace1))
01317             {
01318                 keepTriangles[triI] = false;
01319                 nDangling++;
01320             }
01321         }
01322     }
01323     return nDangling;
01324 }
01325 
01326 
01327 Foam::triSurface Foam::isoSurfaceCell::subsetMesh
01328 (
01329     const triSurface& s,
01330     const labelList& newToOldFaces,
01331     labelList& oldToNewPoints,
01332     labelList& newToOldPoints
01333 )
01334 {
01335     const boolList include
01336     (
01337         createWithValues<boolList>
01338         (
01339             s.size(),
01340             false,
01341             newToOldFaces,
01342             true
01343         )
01344     );
01345 
01346     newToOldPoints.setSize(s.points().size());
01347     oldToNewPoints.setSize(s.points().size());
01348     oldToNewPoints = -1;
01349     {
01350         label pointI = 0;
01351 
01352         forAll(include, oldFacei)
01353         {
01354             if (include[oldFacei])
01355             {
01356                 // Renumber labels for triangle
01357                 const labelledTri& tri = s[oldFacei];
01358 
01359                 forAll(tri, fp)
01360                 {
01361                     label oldPointI = tri[fp];
01362 
01363                     if (oldToNewPoints[oldPointI] == -1)
01364                     {
01365                         oldToNewPoints[oldPointI] = pointI;
01366                         newToOldPoints[pointI++] = oldPointI;
01367                     }
01368                 }
01369             }
01370         }
01371         newToOldPoints.setSize(pointI);
01372     }
01373 
01374     // Extract points
01375     pointField newPoints(newToOldPoints.size());
01376     forAll(newToOldPoints, i)
01377     {
01378         newPoints[i] = s.points()[newToOldPoints[i]];
01379     }
01380     // Extract faces
01381     List<labelledTri> newTriangles(newToOldFaces.size());
01382 
01383     forAll(newToOldFaces, i)
01384     {
01385         // Get old vertex labels
01386         const labelledTri& tri = s[newToOldFaces[i]];
01387 
01388         newTriangles[i][0] = oldToNewPoints[tri[0]];
01389         newTriangles[i][1] = oldToNewPoints[tri[1]];
01390         newTriangles[i][2] = oldToNewPoints[tri[2]];
01391         newTriangles[i].region() = tri.region();
01392     }
01393 
01394     // Reuse storage.
01395     return triSurface(newTriangles, s.patches(), newPoints, true);
01396 }
01397 
01398 
01399 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
01400 
01401 Foam::isoSurfaceCell::isoSurfaceCell
01402 (
01403     const polyMesh& mesh,
01404     const scalarField& cVals,
01405     const scalarField& pVals,
01406     const scalar iso,
01407     const bool regularise,
01408     const scalar mergeTol
01409 )
01410 :
01411     mesh_(mesh),
01412     iso_(iso),
01413     mergeDistance_(mergeTol*mesh.bounds().mag())
01414 {
01415     // Determine if cell is tet
01416     PackedBoolList isTet(mesh_.nCells());
01417     {
01418         tetMatcher tet;
01419 
01420         forAll(isTet, cellI)
01421         {
01422             if (tet.isA(mesh_, cellI))
01423             {
01424                 isTet.set(cellI, 1);
01425             }
01426         }
01427     }
01428 
01429     // Determine if point is on boundary. Points on boundaries are never
01430     // snapped. Coupled boundaries are handled explicitly so not marked here.
01431     PackedBoolList isBoundaryPoint(mesh_.nPoints());
01432     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
01433     forAll(patches, patchI)
01434     {
01435         const polyPatch& pp = patches[patchI];
01436 
01437         if (!pp.coupled())
01438         {
01439             label faceI = pp.start();
01440             forAll(pp, i)
01441             {
01442                 const face& f = mesh_.faces()[faceI++];
01443 
01444                 forAll(f, fp)
01445                 {
01446                     isBoundaryPoint.set(f[fp], 1);
01447                 }
01448             }
01449         }
01450     }
01451 
01452 
01453     // Determine if any cut through cell
01454     calcCutTypes(isTet, cVals, pVals);
01455 
01456     DynamicList<point> snappedPoints(nCutCells_);
01457 
01458     // Per cc -1 or a point inside snappedPoints.
01459     labelList snappedCc;
01460     if (regularise)
01461     {
01462         calcSnappedCc
01463         (
01464             isTet,
01465             cVals,
01466             pVals,
01467             snappedPoints,
01468             snappedCc
01469         );
01470     }
01471     else
01472     {
01473         snappedCc.setSize(mesh_.nCells());
01474         snappedCc = -1;
01475     }
01476 
01477     if (debug)
01478     {
01479         Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()
01480             << " cell centres to intersection." << endl;
01481     }
01482 
01483     snappedPoints.shrink();
01484     label nCellSnaps = snappedPoints.size();
01485 
01486     // Per point -1 or a point inside snappedPoints.
01487     labelList snappedPoint;
01488     if (regularise)
01489     {
01490         calcSnappedPoint
01491         (
01492             isBoundaryPoint,
01493             isTet,
01494             cVals,
01495             pVals,
01496             snappedPoints,
01497             snappedPoint
01498         );
01499     }
01500     else
01501     {
01502         snappedPoint.setSize(mesh_.nPoints());
01503         snappedPoint = -1;
01504     }
01505 
01506     if (debug)
01507     {
01508         Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps
01509             << " vertices to intersection." << endl;
01510     }
01511 
01512 
01513 
01514     DynamicList<point> triPoints(nCutCells_);
01515     DynamicList<label> triMeshCells(nCutCells_);
01516 
01517     generateTriPoints
01518     (
01519         cVals,
01520         pVals,
01521 
01522         mesh_.cellCentres(),
01523         mesh_.points(),
01524 
01525         snappedPoints,
01526         snappedCc,
01527         snappedPoint,
01528 
01529         triPoints,
01530         triMeshCells
01531     );
01532 
01533     if (debug)
01534     {
01535         Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
01536             << " unmerged triangles." << endl;
01537     }
01538 
01539     // Merge points and compact out non-valid triangles
01540     labelList triMap;           // merged to unmerged triangle
01541     triSurface::operator=
01542     (
01543         stitchTriPoints
01544         (
01545             true,               // check for duplicate tris
01546             triPoints,
01547             triPointMergeMap_,  // unmerged to merged point
01548             triMap
01549         )
01550     );
01551 
01552     if (debug)
01553     {
01554         Pout<< "isoSurfaceCell : generated " << triMap.size()
01555             << " merged triangles." << endl;
01556     }
01557 
01558     meshCells_.setSize(triMap.size());
01559     forAll(triMap, i)
01560     {
01561         meshCells_[i] = triMeshCells[triMap[i]];
01562     }
01563 
01564     if (debug)
01565     {
01566         Pout<< "isoSurfaceCell : checking " << size()
01567             << " triangles for validity." << endl;
01568 
01569         forAll(*this, triI)
01570         {
01571             // Copied from surfaceCheck
01572             validTri(*this, triI);
01573         }
01574     }
01575 
01576 
01577 //if (false)
01578 {
01579     List<FixedList<label, 3> > faceEdges;
01580     labelList edgeFace0, edgeFace1;
01581     Map<labelList> edgeFacesRest;
01582 
01583 
01584     while (true)
01585     {
01586         // Calculate addressing
01587         calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
01588 
01589         // See if any dangling triangles
01590         boolList keepTriangles;
01591         label nDangling = markDanglingTriangles
01592         (
01593             faceEdges,
01594             edgeFace0,
01595             edgeFace1,
01596             edgeFacesRest,
01597             keepTriangles
01598         );
01599 
01600         if (debug)
01601         {
01602             Pout<< "isoSurfaceCell : detected " << nDangling
01603                 << " dangling triangles." << endl;
01604         }
01605 
01606         if (nDangling == 0)
01607         {
01608             break;
01609         }
01610 
01611         // Create face map (new to old)
01612         labelList subsetTriMap(findIndices(keepTriangles, true));
01613 
01614         labelList subsetPointMap;
01615         labelList reversePointMap;
01616         triSurface::operator=
01617         (
01618             subsetMesh
01619             (
01620                 *this,
01621                 subsetTriMap,
01622                 reversePointMap,
01623                 subsetPointMap
01624             )
01625         );
01626         meshCells_ = labelField(meshCells_, subsetTriMap);
01627         inplaceRenumber(reversePointMap, triPointMergeMap_);
01628     }
01629 
01630     orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
01631 }
01632 
01633     //combineCellTriangles();
01634 }
01635 
01636 
01641 //void Foam::isoSurfaceCell::combineCellTriangles()
01642 //{
01643 //    if (size())
01644 //    {
01645 //        DynamicList<labelledTri> newTris(size());
01646 //        DynamicList<label> newTriToCell(size());
01647 //
01648 //        label startTriI = 0;
01649 //
01650 //        DynamicList<labelledTri> tris;
01651 //
01652 //        for (label triI = 1; triI <= meshCells_.size(); triI++)
01653 //        {
01654 //            if
01655 //            (
01656 //                triI == meshCells_.size()
01657 //             || meshCells_[triI] != meshCells_[startTriI]
01658 //            )
01659 //            {
01660 //                label nTris = triI-startTriI;
01661 //
01662 //                if (nTris == 1)
01663 //                {
01664 //                    newTris.append(operator[](startTriI));
01665 //                    newTriToCell.append(meshCells_[startTriI]);
01666 //                }
01667 //                else
01668 //                {
01669 //                    // Collect from startTriI to triI in a triSurface
01670 //                    tris.clear();
01671 //                    for (label i = startTriI; i < triI; i++)
01672 //                    {
01673 //                        tris.append(operator[](i));
01674 //                    }
01675 //                    triSurface cellTris(tris, patches(), points());
01676 //                    tris.clear();
01677 //
01678 //                    // Get outside
01679 //                    const labelListList& loops = cellTris.edgeLoops();
01680 //
01681 //                    forAll(loops, i)
01682 //                    {
01683 //                        // Do proper triangulation of loop
01684 //                        face loop(renumber(cellTris.meshPoints(), loops[i]));
01685 //
01686 //                        faceTriangulation faceTris
01687 //                        (
01688 //                            points(),
01689 //                            loop,
01690 //                            true
01691 //                        );
01692 //
01693 //                        // Copy into newTris
01694 //                        forAll(faceTris, faceTriI)
01695 //                        {
01696 //                            const triFace& tri = faceTris[faceTriI];
01697 //
01698 //                            newTris.append
01699 //                            (
01700 //                                labelledTri
01701 //                                (
01702 //                                    tri[0],
01703 //                                    tri[1],
01704 //                                    tri[2],
01705 //                                    operator[](startTriI).region()
01706 //                                )
01707 //                            );
01708 //                            newTriToCell.append(meshCells_[startTriI]);
01709 //                        }
01710 //                    }
01711 //                }
01712 //
01713 //                startTriI = triI;
01714 //            }
01715 //        }
01716 //        newTris.shrink();
01717 //        newTriToCell.shrink();
01718 //
01719 //        // Compact
01720 //        pointField newPoints(points().size());
01721 //        label newPointI = 0;
01722 //        labelList oldToNewPoint(points().size(), -1);
01723 //
01724 //        forAll(newTris, i)
01725 //        {
01726 //            labelledTri& tri = newTris[i];
01727 //            forAll(tri, j)
01728 //            {
01729 //                label pointI = tri[j];
01730 //
01731 //                if (oldToNewPoint[pointI] == -1)
01732 //                {
01733 //                    oldToNewPoint[pointI] = newPointI;
01734 //                    newPoints[newPointI++] = points()[pointI];
01735 //                }
01736 //                tri[j] = oldToNewPoint[pointI];
01737 //            }
01738 //        }
01739 //        newPoints.setSize(newPointI);
01740 //
01741 //        triSurface::operator=
01742 //        (
01743 //            triSurface
01744 //            (
01745 //                newTris,
01746 //                patches(),
01747 //                newPoints,
01748 //                true
01749 //            )
01750 //        );
01751 //        meshCells_.transfer(newTriToCell);
01752 //    }
01753 //}
01755 
01756 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines