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

polyDualMesh.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 InClass
00025     polyDualMesh
00026 
00027 \*---------------------------------------------------------------------------*/
00028 
00029 #include "polyDualMesh.H"
00030 #include <meshTools/meshTools.H>
00031 #include <OpenFOAM/OFstream.H>
00032 #include <OpenFOAM/Time.H>
00033 #include <OpenFOAM/SortableList.H>
00034 #include <meshTools/pointSet.H>
00035 
00036 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00037 
00038 namespace Foam
00039 {
00040     defineTypeNameAndDebug(polyDualMesh, 0);
00041 }
00042 
00043 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00044 
00045 // Determine order for faces:
00046 // - upper-triangular order for internal faces
00047 // - external faces after internal faces (were already so)
00048 Foam::labelList Foam::polyDualMesh::getFaceOrder
00049 (
00050     const labelList& faceOwner,
00051     const labelList& faceNeighbour,
00052     const cellList& cells,
00053     label& nInternalFaces
00054 )
00055 {
00056     labelList oldToNew(faceOwner.size(), -1);
00057 
00058     // First unassigned face
00059     label newFaceI = 0;
00060 
00061     forAll(cells, cellI)
00062     {
00063         const labelList& cFaces = cells[cellI];
00064 
00065         SortableList<label> nbr(cFaces.size());
00066 
00067         forAll(cFaces, i)
00068         {
00069             label faceI = cFaces[i];
00070 
00071             label nbrCellI = faceNeighbour[faceI];
00072 
00073             if (nbrCellI != -1)
00074             {
00075                 // Internal face. Get cell on other side.
00076                 if (nbrCellI == cellI)
00077                 {
00078                     nbrCellI = faceOwner[faceI];
00079                 }
00080 
00081                 if (cellI < nbrCellI)
00082                 {
00083                     // CellI is master
00084                     nbr[i] = nbrCellI;
00085                 }
00086                 else
00087                 {
00088                     // nbrCell is master. Let it handle this face.
00089                     nbr[i] = -1;
00090                 }
00091             }
00092             else
00093             {
00094                 // External face. Do later.
00095                 nbr[i] = -1;
00096             }
00097         }
00098 
00099         nbr.sort();
00100 
00101         forAll(nbr, i)
00102         {
00103             if (nbr[i] != -1)
00104             {
00105                 oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
00106             }
00107         }
00108     }
00109 
00110     nInternalFaces = newFaceI;
00111 
00112     Pout<< "nInternalFaces:" << nInternalFaces << endl;
00113     Pout<< "nFaces:" << faceOwner.size() << endl;
00114     Pout<< "nCells:" << cells.size() << endl;
00115 
00116     // Leave patch faces intact.
00117     for (label faceI = newFaceI; faceI < faceOwner.size(); faceI++)
00118     {
00119         oldToNew[faceI] = faceI;
00120     }
00121 
00122 
00123     // Check done all faces.
00124     forAll(oldToNew, faceI)
00125     {
00126         if (oldToNew[faceI] == -1)
00127         {
00128             FatalErrorIn
00129             (
00130                 "polyDualMesh::getFaceOrder"
00131                 "(const labelList&, const labelList&, const label) const"
00132             )   << "Did not determine new position"
00133                 << " for face " << faceI
00134                 << abort(FatalError);
00135         }
00136     }
00137 
00138     return oldToNew;
00139 }
00140 
00141 
00142 // Get the two edges on faceI using pointI. Returns them such that the order
00143 // - otherVertex of e0
00144 // - pointI
00145 // - otherVertex(pointI) of e1
00146 // is in face order
00147 void Foam::polyDualMesh::getPointEdges
00148 (
00149     const primitivePatch& patch,
00150     const label faceI,
00151     const label pointI,
00152     label& e0,
00153     label& e1
00154 )
00155 {
00156     const labelList& fEdges = patch.faceEdges()[faceI];
00157     const face& f = patch.localFaces()[faceI];
00158 
00159     e0 = -1;
00160     e1 = -1;
00161 
00162     forAll(fEdges, i)
00163     {
00164         label edgeI = fEdges[i];
00165 
00166         const edge& e = patch.edges()[edgeI];
00167 
00168         if (e[0] == pointI)
00169         {
00170             // One of the edges using pointI. Check which one.
00171             label index = findIndex(f, pointI);
00172 
00173             if (f.nextLabel(index) == e[1])
00174             {
00175                 e1 = edgeI;
00176             }
00177             else
00178             {
00179                 e0 = edgeI;
00180             }
00181 
00182             if (e0 != -1 && e1 != -1)
00183             {
00184                 return;
00185             }
00186         }
00187         else if (e[1] == pointI)
00188         {
00189             // One of the edges using pointI. Check which one.
00190             label index = findIndex(f, pointI);
00191 
00192             if (f.nextLabel(index) == e[0])
00193             {
00194                 e1 = edgeI;
00195             }
00196             else
00197             {
00198                 e0 = edgeI;
00199             }
00200 
00201             if (e0 != -1 && e1 != -1)
00202             {
00203                 return;
00204             }
00205         }
00206     }
00207 
00208     FatalErrorIn("getPointEdges") << "Cannot find two edges on face:" << faceI
00209         << " vertices:" << patch.localFaces()[faceI]
00210         << " that uses point:" << pointI
00211         << abort(FatalError);
00212 }
00213 
00214 
00215 // Collect the face on an external point of the patch.
00216 Foam::labelList Foam::polyDualMesh::collectPatchSideFace
00217 (
00218     const polyPatch& patch,
00219     const label patchToDualOffset,
00220     const labelList& edgeToDualPoint,
00221     const labelList& pointToDualPoint,
00222     const label pointI,
00223 
00224     label& edgeI
00225 )
00226 {
00227     // Construct face by walking around e[eI] starting from
00228     // patchEdgeI
00229 
00230     label meshPointI = patch.meshPoints()[pointI];
00231     const labelList& pFaces = patch.pointFaces()[pointI];
00232 
00233     DynamicList<label> dualFace;
00234 
00235     if (pointToDualPoint[meshPointI] >= 0)
00236     {
00237         // Number of pFaces + 2 boundary edge + feature point
00238         dualFace.setCapacity(pFaces.size()+2+1);
00239         // Store dualVertex for feature edge
00240         dualFace.append(pointToDualPoint[meshPointI]);
00241     }
00242     else
00243     {
00244         dualFace.setCapacity(pFaces.size()+2);
00245     }
00246 
00247     // Store dual vertex for starting edge.
00248     if (edgeToDualPoint[patch.meshEdges()[edgeI]] < 0)
00249     {
00250         FatalErrorIn("polyDualMesh::collectPatchSideFace") << edgeI
00251             << abort(FatalError);
00252     }
00253 
00254     dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
00255 
00256     label faceI = patch.edgeFaces()[edgeI][0];
00257 
00258     // Check order of vertices of edgeI in face to see if we need to reverse.
00259     bool reverseFace;
00260 
00261     label e0, e1;
00262     getPointEdges(patch, faceI, pointI, e0, e1);
00263 
00264     if (e0 == edgeI)
00265     {
00266         reverseFace = true;
00267     }
00268     else
00269     {
00270         reverseFace = false;
00271     }
00272 
00273     while (true)
00274     {
00275         // Store dual vertex for faceI.
00276         dualFace.append(faceI + patchToDualOffset);
00277 
00278         // Cross face to other edge on pointI
00279         label e0, e1;
00280         getPointEdges(patch, faceI, pointI, e0, e1);
00281 
00282         if (e0 == edgeI)
00283         {
00284             edgeI = e1;
00285         }
00286         else
00287         {
00288             edgeI = e0;
00289         }
00290 
00291         if (edgeToDualPoint[patch.meshEdges()[edgeI]] >= 0)
00292         {
00293             // Feature edge. Insert dual vertex for edge.
00294             dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
00295         }
00296 
00297         const labelList& eFaces = patch.edgeFaces()[edgeI];
00298 
00299         if (eFaces.size() == 1)
00300         {
00301             // Reached other edge of patch
00302             break;
00303         }
00304 
00305         // Cross edge to other face.
00306         if (eFaces[0] == faceI)
00307         {
00308             faceI = eFaces[1];
00309         }
00310         else
00311         {
00312             faceI = eFaces[0];
00313         }
00314     }
00315 
00316     dualFace.shrink();
00317 
00318     if (reverseFace)
00319     {
00320         reverse(dualFace);
00321     }
00322 
00323     return dualFace;
00324 }
00325 
00326 
00327 // Collect face around pointI which is not on the outside of the patch.
00328 // Returns the vertices of the face and the indices in these vertices of
00329 // any points which are on feature edges.
00330 void Foam::polyDualMesh::collectPatchInternalFace
00331 (
00332     const polyPatch& patch,
00333     const label patchToDualOffset,
00334     const labelList& edgeToDualPoint,
00335 
00336     const label pointI,
00337     const label startEdgeI,
00338 
00339     labelList& dualFace2,
00340     labelList& featEdgeIndices2
00341 )
00342 {
00343     // Construct face by walking around pointI starting from startEdgeI
00344     const labelList& meshEdges = patch.meshEdges();
00345     const labelList& pFaces = patch.pointFaces()[pointI];
00346 
00347     // Vertices of dualFace
00348     DynamicList<label> dualFace(pFaces.size());
00349     // Indices in dualFace of vertices added because of feature edge
00350     DynamicList<label> featEdgeIndices(dualFace.size());
00351 
00352 
00353     label edgeI = startEdgeI;
00354     label faceI = patch.edgeFaces()[edgeI][0];
00355 
00356     // Check order of vertices of edgeI in face to see if we need to reverse.
00357     bool reverseFace;
00358 
00359     label e0, e1;
00360     getPointEdges(patch, faceI, pointI, e0, e1);
00361 
00362     if (e0 == edgeI)
00363     {
00364         reverseFace = true;
00365     }
00366     else
00367     {
00368         reverseFace = false;
00369     }
00370 
00371     while (true)
00372     {
00373         // Insert dual vertex for face
00374         dualFace.append(faceI + patchToDualOffset);
00375 
00376         // Cross face to other edge on pointI
00377         label e0, e1;
00378         getPointEdges(patch, faceI, pointI, e0, e1);
00379 
00380         if (e0 == edgeI)
00381         {
00382             edgeI = e1;
00383         }
00384         else
00385         {
00386             edgeI = e0;
00387         }
00388 
00389         if (edgeToDualPoint[meshEdges[edgeI]] >= 0)
00390         {
00391             // Feature edge. Insert dual vertex for edge.
00392             dualFace.append(edgeToDualPoint[meshEdges[edgeI]]);
00393             featEdgeIndices.append(dualFace.size()-1);
00394         }
00395 
00396         if (edgeI == startEdgeI)
00397         {
00398             break;
00399         }
00400 
00401         // Cross edge to other face.
00402         const labelList& eFaces = patch.edgeFaces()[edgeI];
00403 
00404         if (eFaces[0] == faceI)
00405         {
00406             faceI = eFaces[1];
00407         }
00408         else
00409         {
00410             faceI = eFaces[0];
00411         }
00412     }
00413 
00414     dualFace2.transfer(dualFace);
00415 
00416     featEdgeIndices2.transfer(featEdgeIndices);
00417 
00418     if (reverseFace)
00419     {
00420         reverse(dualFace2);
00421 
00422         // Correct featEdgeIndices for change in dualFace2
00423         forAll(featEdgeIndices2, i)
00424         {
00425             featEdgeIndices2[i] = dualFace2.size() -1 - featEdgeIndices2[i];
00426         }
00427         // Reverse indices (might not be nessecary but do anyway)
00428         reverse(featEdgeIndices2);
00429     }
00430 }
00431 
00432 
00433 void Foam::polyDualMesh::splitFace
00434 (
00435     const polyPatch& patch,
00436     const labelList& pointToDualPoint,
00437 
00438     const label pointI,
00439     const labelList& dualFace,
00440     const labelList& featEdgeIndices,
00441 
00442     DynamicList<face>& dualFaces,
00443     DynamicList<label>& dualOwner,
00444     DynamicList<label>& dualNeighbour,
00445     DynamicList<label>& dualRegion
00446 )
00447 {
00448 
00449     // Split face because of feature edges/points
00450     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00451     label meshPointI = patch.meshPoints()[pointI];
00452 
00453     if (pointToDualPoint[meshPointI] >= 0)
00454     {
00455         // Feature point. Do face-centre decomposition.
00456 
00457         if (featEdgeIndices.size() < 2)
00458         {
00459             // Feature point but no feature edges. Not handled at the moment
00460             dualFaces.append(face(dualFace));
00461             dualOwner.append(meshPointI);
00462             dualNeighbour.append(-1);
00463             dualRegion.append(patch.index());
00464         }
00465         else
00466         {
00467             // Do 'face-centre' decomposition. Start from first feature
00468             // edge create face up until next feature edge.
00469 
00470             forAll(featEdgeIndices, i)
00471             {
00472                 label startFp = featEdgeIndices[i];
00473 
00474                 label endFp = featEdgeIndices[(i+1) % featEdgeIndices.size()];
00475 
00476                 // Collect face from startFp to endFp
00477                 label sz = 0;
00478 
00479                 if (endFp > startFp)
00480                 {
00481                     sz = endFp - startFp + 2;
00482                 }
00483                 else
00484                 {
00485                     sz = endFp + dualFace.size() - startFp + 2;
00486                 }
00487                 face subFace(sz);
00488 
00489                 // feature point becomes face centre.
00490                 subFace[0] = pointToDualPoint[patch.meshPoints()[pointI]];
00491 
00492                 // Copy from startFp up to endFp.
00493                 for (label subFp = 1; subFp < subFace.size(); subFp++)
00494                 {
00495                     subFace[subFp] = dualFace[startFp];
00496 
00497                     startFp = (startFp+1) % dualFace.size();
00498                 }
00499 
00500                 dualFaces.append(face(subFace));
00501                 dualOwner.append(meshPointI);
00502                 dualNeighbour.append(-1);
00503                 dualRegion.append(patch.index());
00504             }
00505         }
00506     }
00507     else
00508     {
00509         // No feature point. Check if feature edges.
00510         if (featEdgeIndices.size() < 2)
00511         {
00512             // Not enough feature edges. No split.
00513             dualFaces.append(face(dualFace));
00514             dualOwner.append(meshPointI);
00515             dualNeighbour.append(-1);
00516             dualRegion.append(patch.index());
00517         }
00518         else
00519         {
00520             // Multiple feature edges but no feature point.
00521             // Split face along feature edges. Gives weird result if
00522             // number of feature edges > 2.
00523 
00524             // Storage for new face
00525             DynamicList<label> subFace(dualFace.size());
00526 
00527             forAll(featEdgeIndices, featI)
00528             {
00529                 label startFp = featEdgeIndices[featI];
00530                 label endFp = featEdgeIndices[featEdgeIndices.fcIndex(featI)];
00531 
00532                 label fp = startFp;
00533 
00534                 while (true)
00535                 {
00536                     label vertI = dualFace[fp];
00537 
00538                     subFace.append(vertI);
00539 
00540                     if (fp == endFp)
00541                     {
00542                         break;
00543                     }
00544 
00545                     fp = dualFace.fcIndex(fp);
00546                 }
00547 
00548                 if (subFace.size() > 2)
00549                 {
00550                     // Enough vertices to create a face from.
00551                     subFace.shrink();
00552 
00553                     dualFaces.append(face(subFace));
00554                     dualOwner.append(meshPointI);
00555                     dualNeighbour.append(-1);
00556                     dualRegion.append(patch.index());
00557 
00558                     subFace.clear();
00559                 }
00560             }
00561             // Check final face.
00562             if (subFace.size() > 2)
00563             {
00564                 // Enough vertices to create a face from.
00565                 subFace.shrink();
00566 
00567                 dualFaces.append(face(subFace));
00568                 dualOwner.append(meshPointI);
00569                 dualNeighbour.append(-1);
00570                 dualRegion.append(patch.index());
00571 
00572                 subFace.clear();
00573             }
00574         }
00575     }
00576 }
00577 
00578 
00579 // Create boundary face for every point in patch
00580 void Foam::polyDualMesh::dualPatch
00581 (
00582     const polyPatch& patch,
00583     const label patchToDualOffset,
00584     const labelList& edgeToDualPoint,
00585     const labelList& pointToDualPoint,
00586 
00587     const pointField& dualPoints,
00588 
00589     DynamicList<face>& dualFaces,
00590     DynamicList<label>& dualOwner,
00591     DynamicList<label>& dualNeighbour,
00592     DynamicList<label>& dualRegion
00593 )
00594 {
00595     const labelList& meshEdges = patch.meshEdges();
00596 
00597     // Whether edge has been done.
00598     // 0 : not
00599     // 1 : done e.start()
00600     // 2 : done e.end()
00601     // 3 : done both
00602     labelList doneEdgeSide(meshEdges.size(), 0);
00603 
00604     boolList donePoint(patch.nPoints(), false);
00605 
00606 
00607     // Do points on edge of patch
00608     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
00609 
00610     forAll(doneEdgeSide, patchEdgeI)
00611     {
00612         const labelList& eFaces = patch.edgeFaces()[patchEdgeI];
00613 
00614         if (eFaces.size() == 1)
00615         {
00616             const edge& e = patch.edges()[patchEdgeI];
00617 
00618             forAll(e, eI)
00619             {
00620                 label bitMask = 1<<eI;
00621 
00622                 if ((doneEdgeSide[patchEdgeI] & bitMask) == 0)
00623                 {
00624                     // Construct face by walking around e[eI] starting from
00625                     // patchEdgeI
00626 
00627                     label pointI = e[eI];
00628 
00629                     label edgeI = patchEdgeI;
00630                     labelList dualFace
00631                     (
00632                         collectPatchSideFace
00633                         (
00634                             patch,
00635                             patchToDualOffset,
00636                             edgeToDualPoint,
00637                             pointToDualPoint,
00638 
00639                             pointI,
00640                             edgeI
00641                         )
00642                     );
00643 
00644                     // Now edgeI is end edge. Mark as visited
00645                     if (patch.edges()[edgeI][0] == pointI)
00646                     {
00647                         doneEdgeSide[edgeI] |= 1;
00648                     }
00649                     else
00650                     {
00651                         doneEdgeSide[edgeI] |= 2;
00652                     }
00653 
00654                     dualFaces.append(face(dualFace));
00655                     dualOwner.append(patch.meshPoints()[pointI]);
00656                     dualNeighbour.append(-1);
00657                     dualRegion.append(patch.index());
00658 
00659                     doneEdgeSide[patchEdgeI] |= bitMask;
00660                     donePoint[pointI] = true;
00661                 }
00662             }
00663         }
00664     }
00665 
00666 
00667 
00668     // Do patch-internal points
00669     // ~~~~~~~~~~~~~~~~~~~~~~~~
00670 
00671     forAll(donePoint, pointI)
00672     {
00673         if (!donePoint[pointI])
00674         {
00675             labelList dualFace, featEdgeIndices;
00676 
00677             collectPatchInternalFace
00678             (
00679                 patch,
00680                 patchToDualOffset,
00681                 edgeToDualPoint,
00682                 pointI,
00683                 patch.pointEdges()[pointI][0],  // Arbitrary starting edge
00684 
00685                 dualFace,
00686                 featEdgeIndices
00687             );
00688 
00689             //- Either keep in one piece or split face according to feature.
00690 
00692             //dualFaces.append(face(dualFace));
00693             //dualOwner.append(patch.meshPoints()[pointI]);
00694             //dualNeighbour.append(-1);
00695             //dualRegion.append(patch.index());
00696 
00697             splitFace
00698             (
00699                 patch,
00700                 pointToDualPoint,
00701                 pointI,
00702                 dualFace,
00703                 featEdgeIndices,
00704 
00705                 dualFaces,
00706                 dualOwner,
00707                 dualNeighbour,
00708                 dualRegion
00709             );
00710 
00711             donePoint[pointI] = true;
00712         }
00713     }
00714 }
00715 
00716 
00717 void Foam::polyDualMesh::calcDual
00718 (
00719     const polyMesh& mesh,
00720     const labelList& featureEdges,
00721     const labelList& featurePoints
00722 )
00723 {
00724     const polyBoundaryMesh& patches = mesh.boundaryMesh();
00725     const label nIntFaces = mesh.nInternalFaces();
00726 
00727 
00728     // Get patch for all of outside
00729     primitivePatch allBoundary
00730     (
00731         SubList<face>
00732         (
00733             mesh.faces(),
00734             mesh.nFaces() - nIntFaces,
00735             nIntFaces
00736         ),
00737         mesh.points()
00738     );
00739 
00740     // Correspondence from patch edge to mesh edge.
00741     labelList meshEdges
00742     (
00743         allBoundary.meshEdges
00744         (
00745             mesh.edges(),
00746             mesh.pointEdges()
00747         )
00748     );
00749 
00750     {
00751         pointSet nonManifoldPoints
00752         (
00753             mesh,
00754             "nonManifoldPoints",
00755             meshEdges.size() / 100
00756         );
00757 
00758         allBoundary.checkPointManifold(true, &nonManifoldPoints);
00759 
00760         if (nonManifoldPoints.size())
00761         {
00762             nonManifoldPoints.write();
00763 
00764             FatalErrorIn
00765             (
00766                 "polyDualMesh::calcDual(const polyMesh&, const labelList&"
00767                 ", const labelList&)"
00768             )   << "There are " << nonManifoldPoints.size() << " points where"
00769                 << " the outside of the mesh is non-manifold." << nl
00770                 << "Such a mesh cannot be converted to a dual." << nl
00771                 << "Writing points at non-manifold sites to pointSet "
00772                 << nonManifoldPoints.name()
00773                 << exit(FatalError);
00774         }
00775     }
00776 
00777 
00778     // Assign points
00779     // ~~~~~~~~~~~~~
00780 
00781     // mesh label   dualMesh vertex
00782     // ----------   ---------------
00783     // cellI        cellI
00784     // faceI        nCells+faceI-nIntFaces
00785     // featEdgeI    nCells+nFaces-nIntFaces+featEdgeI
00786     // featPointI   nCells+nFaces-nIntFaces+nFeatEdges+featPointI
00787 
00788     pointField dualPoints
00789     (
00790         mesh.nCells()                           // cell centres
00791       + mesh.nFaces() - nIntFaces               // boundary face centres
00792       + featureEdges.size()                     // additional boundary edges
00793       + featurePoints.size()                    // additional boundary points
00794     );
00795 
00796     label dualPointI = 0;
00797 
00798 
00799     // Cell centres.
00800     const pointField& cellCentres = mesh.cellCentres();
00801 
00802     cellPoint_.setSize(cellCentres.size());
00803 
00804     forAll(cellCentres, cellI)
00805     {
00806         cellPoint_[cellI] = dualPointI;
00807         dualPoints[dualPointI++] = cellCentres[cellI];
00808     }
00809 
00810     // Boundary faces centres
00811     const pointField& faceCentres = mesh.faceCentres();
00812 
00813     boundaryFacePoint_.setSize(mesh.nFaces() - nIntFaces);
00814 
00815     for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
00816     {
00817         boundaryFacePoint_[faceI - nIntFaces] = dualPointI;
00818         dualPoints[dualPointI++] = faceCentres[faceI];
00819     }
00820 
00821     // Edge status:
00822     //  >0 : dual point label (edge is feature edge)
00823     //  -1 : is boundary edge.
00824     //  -2 : is internal edge.
00825     labelList edgeToDualPoint(mesh.nEdges(), -2);
00826 
00827     forAll(meshEdges, patchEdgeI)
00828     {
00829         label edgeI = meshEdges[patchEdgeI];
00830 
00831         edgeToDualPoint[edgeI] = -1;
00832     }
00833 
00834     forAll(featureEdges, i)
00835     {
00836         label edgeI = featureEdges[i];
00837 
00838         const edge& e = mesh.edges()[edgeI];
00839 
00840         edgeToDualPoint[edgeI] = dualPointI;
00841         dualPoints[dualPointI++] = e.centre(mesh.points());
00842     }
00843 
00844 
00845 
00846     // Point status:
00847     //  >0 : dual point label (point is feature point)
00848     //  -1 : is point on edge of patch
00849     //  -2 : is point on patch (but not on edge)
00850     //  -3 : is internal point.
00851     labelList pointToDualPoint(mesh.nPoints(), -3);
00852 
00853     forAll(patches, patchI)
00854     {
00855         const labelList& meshPoints = patches[patchI].meshPoints();
00856 
00857         forAll(meshPoints, i)
00858         {
00859             pointToDualPoint[meshPoints[i]] = -2;
00860         }
00861 
00862         const labelListList& loops = patches[patchI].edgeLoops();
00863 
00864         forAll(loops, i)
00865         {
00866             const labelList& loop = loops[i];
00867 
00868             forAll(loop, j)
00869             {
00870                  pointToDualPoint[meshPoints[loop[j]]] = -1;
00871             }
00872         }
00873     }
00874 
00875     forAll(featurePoints, i)
00876     {
00877         label pointI = featurePoints[i];
00878 
00879         pointToDualPoint[pointI] = dualPointI;
00880         dualPoints[dualPointI++] = mesh.points()[pointI];
00881     }
00882 
00883 
00884     // Storage for new faces.
00885     // Dynamic sized since we don't know size.
00886 
00887     DynamicList<face> dynDualFaces(mesh.nEdges());
00888     DynamicList<label> dynDualOwner(mesh.nEdges());
00889     DynamicList<label> dynDualNeighbour(mesh.nEdges());
00890     DynamicList<label> dynDualRegion(mesh.nEdges());
00891 
00892 
00893     // Generate faces from edges on the boundary
00894     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00895 
00896     forAll(meshEdges, patchEdgeI)
00897     {
00898         label edgeI = meshEdges[patchEdgeI];
00899 
00900         const edge& e = mesh.edges()[edgeI];
00901 
00902         label owner = -1;
00903         label neighbour = -1;
00904 
00905         if (e[0] < e[1])
00906         {
00907             owner = e[0];
00908             neighbour = e[1];
00909         }
00910         else
00911         {
00912             owner = e[1];
00913             neighbour = e[0];
00914         }
00915 
00916         // Find the boundary faces using the edge.
00917         const labelList& patchFaces = allBoundary.edgeFaces()[patchEdgeI];
00918 
00919         if (patchFaces.size() != 2)
00920         {
00921             FatalErrorIn("polyDualMesh::calcDual")
00922                 << "Cannot handle edges with " << patchFaces.size()
00923                 << " connected boundary faces."
00924                 << abort(FatalError);
00925         }
00926 
00927         label face0 = patchFaces[0] + nIntFaces;
00928         const face& f0 = mesh.faces()[face0];
00929 
00930         label face1 = patchFaces[1] + nIntFaces;
00931 
00932         // We want to start walking from patchFaces[0] or patchFaces[1],
00933         // depending on which one uses owner,neighbour in the right order.
00934 
00935         label startFaceI = -1;
00936         label endFaceI = -1;
00937 
00938         label index = findIndex(f0, neighbour);
00939 
00940         if (f0.nextLabel(index) == owner)
00941         {
00942             startFaceI = face0;
00943             endFaceI = face1;
00944         }
00945         else
00946         {
00947             startFaceI = face1;
00948             endFaceI = face0;
00949         }
00950 
00951         // Now walk from startFaceI to cell to face to cell etc. until endFaceI
00952 
00953         DynamicList<label> dualFace;
00954 
00955         if (edgeToDualPoint[edgeI] >= 0)
00956         {
00957             // Number of cells + 2 boundary faces + feature edge point
00958             dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2+1);
00959             // Store dualVertex for feature edge
00960             dualFace.append(edgeToDualPoint[edgeI]);
00961         }
00962         else
00963         {
00964             dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2);
00965         }
00966 
00967         // Store dual vertex for starting face.
00968         dualFace.append(mesh.nCells() + startFaceI - nIntFaces);
00969 
00970         label cellI = mesh.faceOwner()[startFaceI];
00971         label faceI = startFaceI;
00972 
00973         while (true)
00974         {
00975             // Store dual vertex from cellI.
00976             dualFace.append(cellI);
00977 
00978             // Cross cell to other face on edge.
00979             label f0, f1;
00980             meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
00981 
00982             if (f0 == faceI)
00983             {
00984                 faceI = f1;
00985             }
00986             else
00987             {
00988                 faceI = f0;
00989             }
00990 
00991             // Cross face to other cell.
00992             if (faceI == endFaceI)
00993             {
00994                 break;
00995             }
00996 
00997             if (mesh.faceOwner()[faceI] == cellI)
00998             {
00999                 cellI = mesh.faceNeighbour()[faceI];
01000             }
01001             else
01002             {
01003                 cellI = mesh.faceOwner()[faceI];
01004             }
01005         }
01006 
01007         // Store dual vertex for endFace.
01008         dualFace.append(mesh.nCells() + endFaceI - nIntFaces);
01009 
01010         dynDualFaces.append(face(dualFace.shrink()));
01011         dynDualOwner.append(owner);
01012         dynDualNeighbour.append(neighbour);
01013         dynDualRegion.append(-1);
01014 
01015         {
01016             // Check orientation.
01017             const face& f = dynDualFaces[dynDualFaces.size()-1];
01018             vector n = f.normal(dualPoints);
01019             if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
01020             {
01021                 WarningIn("calcDual") << "Incorrect orientation"
01022                     << " on boundary edge:" << edgeI
01023                     << mesh.points()[mesh.edges()[edgeI][0]]
01024                     << mesh.points()[mesh.edges()[edgeI][1]]
01025                     << endl;
01026             }
01027         }
01028     }
01029 
01030 
01031     // Generate faces from internal edges
01032     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01033 
01034     forAll(edgeToDualPoint, edgeI)
01035     {
01036         if (edgeToDualPoint[edgeI] == -2)
01037         {
01038             // Internal edge.
01039 
01040             // Find dual owner, neighbour
01041 
01042             const edge& e = mesh.edges()[edgeI];
01043 
01044             label owner = -1;
01045             label neighbour = -1;
01046 
01047             if (e[0] < e[1])
01048             {
01049                 owner = e[0];
01050                 neighbour = e[1];
01051             }
01052             else
01053             {
01054                 owner = e[1];
01055                 neighbour = e[0];
01056             }
01057 
01058             // Get a starting cell
01059             const labelList& eCells = mesh.edgeCells()[edgeI];
01060 
01061             label cellI = eCells[0];
01062 
01063             // Get the two faces on the cell and edge.
01064             label face0, face1;
01065             meshTools::getEdgeFaces(mesh, cellI, edgeI, face0, face1);
01066 
01067             // Find the starting face by looking at the order in which
01068             // the face uses the owner, neighbour
01069             const face& f0 = mesh.faces()[face0];
01070 
01071             label index = findIndex(f0, neighbour);
01072 
01073             bool f0OrderOk = (f0.nextLabel(index) == owner);
01074 
01075             label startFaceI = -1;
01076 
01077             if (f0OrderOk == (mesh.faceOwner()[face0] == cellI))
01078             {
01079                 startFaceI = face0;
01080             }
01081             else
01082             {
01083                 startFaceI = face1;
01084             }
01085 
01086             // Walk face-cell-face until starting face reached.
01087             DynamicList<label> dualFace(mesh.edgeCells()[edgeI].size());
01088 
01089             label faceI = startFaceI;
01090 
01091             while (true)
01092             {
01093                 // Store dual vertex from cellI.
01094                 dualFace.append(cellI);
01095 
01096                 // Cross cell to other face on edge.
01097                 label f0, f1;
01098                 meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
01099 
01100                 if (f0 == faceI)
01101                 {
01102                     faceI = f1;
01103                 }
01104                 else
01105                 {
01106                     faceI = f0;
01107                 }
01108 
01109                 // Cross face to other cell.
01110                 if (faceI == startFaceI)
01111                 {
01112                     break;
01113                 }
01114 
01115                 if (mesh.faceOwner()[faceI] == cellI)
01116                 {
01117                     cellI = mesh.faceNeighbour()[faceI];
01118                 }
01119                 else
01120                 {
01121                     cellI = mesh.faceOwner()[faceI];
01122                 }
01123             }
01124 
01125             dynDualFaces.append(face(dualFace.shrink()));
01126             dynDualOwner.append(owner);
01127             dynDualNeighbour.append(neighbour);
01128             dynDualRegion.append(-1);
01129 
01130             {
01131                 // Check orientation.
01132                 const face& f = dynDualFaces[dynDualFaces.size()-1];
01133                 vector n = f.normal(dualPoints);
01134                 if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
01135                 {
01136                     WarningIn("calcDual") << "Incorrect orientation"
01137                         << " on internal edge:" << edgeI
01138                         << mesh.points()[mesh.edges()[edgeI][0]]
01139                         << mesh.points()[mesh.edges()[edgeI][1]]
01140                         << endl;
01141                 }
01142             }
01143         }
01144     }
01145 
01146     // Dump faces.
01147     if (debug)
01148     {
01149         dynDualFaces.shrink();
01150         dynDualOwner.shrink();
01151         dynDualNeighbour.shrink();
01152         dynDualRegion.shrink();
01153 
01154         OFstream str("dualInternalFaces.obj");
01155         Pout<< "polyDualMesh::calcDual : dumping internal faces to "
01156             << str.name() << endl;
01157 
01158         forAll(dualPoints, dualPointI)
01159         {
01160             meshTools::writeOBJ(str, dualPoints[dualPointI]);
01161         }
01162         forAll(dynDualFaces, dualFaceI)
01163         {
01164             const face& f = dynDualFaces[dualFaceI];
01165 
01166             str<< 'f';
01167             forAll(f, fp)
01168             {
01169                 str<< ' ' << f[fp]+1;
01170             }
01171             str<< nl;
01172         }
01173     }
01174 
01175     const label nInternalFaces = dynDualFaces.size();
01176 
01177     // Outside faces
01178     // ~~~~~~~~~~~~~
01179 
01180     forAll(patches, patchI)
01181     {
01182         const polyPatch& pp = patches[patchI];
01183 
01184         dualPatch
01185         (
01186             pp,
01187             mesh.nCells() + pp.start() - nIntFaces,
01188             edgeToDualPoint,
01189             pointToDualPoint,
01190 
01191             dualPoints,
01192 
01193             dynDualFaces,
01194             dynDualOwner,
01195             dynDualNeighbour,
01196             dynDualRegion
01197         );
01198     }
01199 
01200 
01201     // Transfer face info to straight lists
01202     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01203     faceList dualFaces(dynDualFaces.shrink(), true);
01204     dynDualFaces.clear();
01205 
01206     labelList dualOwner(dynDualOwner.shrink(), true);
01207     dynDualOwner.clear();
01208 
01209     labelList dualNeighbour(dynDualNeighbour.shrink(), true);
01210     dynDualNeighbour.clear();
01211 
01212     labelList dualRegion(dynDualRegion.shrink(), true);
01213     dynDualRegion.clear();
01214 
01215 
01216 
01217     // Dump faces.
01218     if (debug)
01219     {
01220         OFstream str("dualFaces.obj");
01221         Pout<< "polyDualMesh::calcDual : dumping all faces to "
01222             << str.name() << endl;
01223 
01224         forAll(dualPoints, dualPointI)
01225         {
01226             meshTools::writeOBJ(str, dualPoints[dualPointI]);
01227         }
01228         forAll(dualFaces, dualFaceI)
01229         {
01230             const face& f = dualFaces[dualFaceI];
01231 
01232             str<< 'f';
01233             forAll(f, fp)
01234             {
01235                 str<< ' ' << f[fp]+1;
01236             }
01237             str<< nl;
01238         }
01239     }
01240 
01241     // Create cells.
01242     cellList dualCells(mesh.nPoints());
01243 
01244     forAll(dualCells, cellI)
01245     {
01246         dualCells[cellI].setSize(0);
01247     }
01248 
01249     forAll(dualOwner, faceI)
01250     {
01251         label cellI = dualOwner[faceI];
01252 
01253         labelList& cFaces = dualCells[cellI];
01254 
01255         label sz = cFaces.size();
01256         cFaces.setSize(sz+1);
01257         cFaces[sz] = faceI;
01258     }
01259     forAll(dualNeighbour, faceI)
01260     {
01261         label cellI = dualNeighbour[faceI];
01262 
01263         if (cellI != -1)
01264         {
01265             labelList& cFaces = dualCells[cellI];
01266 
01267             label sz = cFaces.size();
01268             cFaces.setSize(sz+1);
01269             cFaces[sz] = faceI;
01270         }
01271     }
01272 
01273 
01274     // Do upper-triangular face ordering. Determines face reordering map and
01275     // number of internal faces.
01276     label dummy;
01277 
01278     labelList oldToNew
01279     (
01280         getFaceOrder
01281         (
01282             dualOwner,
01283             dualNeighbour,
01284             dualCells,
01285             dummy
01286         )
01287     );
01288 
01289     // Reorder faces.
01290     inplaceReorder(oldToNew, dualFaces);
01291     inplaceReorder(oldToNew, dualOwner);
01292     inplaceReorder(oldToNew, dualNeighbour);
01293     inplaceReorder(oldToNew, dualRegion);
01294     forAll(dualCells, cellI)
01295     {
01296         inplaceRenumber(oldToNew, dualCells[cellI]);
01297     }
01298 
01299 
01300     // Create patches
01301     labelList patchSizes(patches.size(), 0);
01302 
01303     forAll(dualRegion, faceI)
01304     {
01305         if (dualRegion[faceI] >= 0)
01306         {
01307             patchSizes[dualRegion[faceI]]++;
01308         }
01309     }
01310 
01311     labelList patchStarts(patches.size(), 0);
01312 
01313     label faceI = nInternalFaces;
01314 
01315     forAll(patches, patchI)
01316     {
01317         patchStarts[patchI] = faceI;
01318 
01319         faceI += patchSizes[patchI];
01320     }
01321 
01322 
01323     Pout<< "nFaces:" << dualFaces.size()
01324         << " patchSizes:" << patchSizes
01325         << " patchStarts:" << patchStarts
01326         << endl;
01327 
01328 
01329     // Add patches. First add zero sized (since mesh still 0 size)
01330     List<polyPatch*> dualPatches(patches.size());
01331 
01332     forAll(patches, patchI)
01333     {
01334         const polyPatch& pp = patches[patchI];
01335 
01336         dualPatches[patchI] = pp.clone
01337         (
01338             boundaryMesh(),
01339             patchI,
01340             0, //patchSizes[patchI],
01341             0  //patchStarts[patchI]
01342         ).ptr();
01343     }
01344     addPatches(dualPatches);
01345 
01346     // Assign to mesh.
01347     resetPrimitives
01348     (
01349         xferMove(dualPoints),
01350         xferMove(dualFaces),
01351         xferMove(dualOwner),
01352         xferMove(dualNeighbour),
01353         patchSizes,
01354         patchStarts
01355     );
01356 }
01357 
01358 
01359 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
01360 
01361 // Construct from components
01362 Foam::polyDualMesh::polyDualMesh(const IOobject& io)
01363 :
01364     polyMesh(io),
01365     cellPoint_
01366     (
01367         IOobject
01368         (
01369             "cellPoint",
01370             time().findInstance(meshDir(), "cellPoint"),
01371             meshSubDir,
01372             *this,
01373             IOobject::MUST_READ,
01374             IOobject::NO_WRITE
01375         )
01376     ),
01377     boundaryFacePoint_
01378     (
01379         IOobject
01380         (
01381             "boundaryFacePoint",
01382             time().findInstance(meshDir(), "boundaryFacePoint"),
01383             meshSubDir,
01384             *this,
01385             IOobject::MUST_READ,
01386             IOobject::NO_WRITE
01387         )
01388     )
01389 {}
01390 
01391 
01392 // Construct from polyMesh
01393 Foam::polyDualMesh::polyDualMesh
01394 (
01395     const polyMesh& mesh,
01396     const labelList& featureEdges,
01397     const labelList& featurePoints
01398 )
01399 :
01400     polyMesh
01401     (
01402         mesh,
01403         xferCopy(pointField()),   // to prevent any warnings "points not allocated"
01404         xferCopy(faceList()),     // to prevent any warnings "faces  not allocated"
01405         xferCopy(cellList())
01406     ),
01407     cellPoint_
01408     (
01409         IOobject
01410         (
01411             "cellPoint",
01412             time().findInstance(meshDir(), "faces"),
01413             meshSubDir,
01414             *this,
01415             IOobject::NO_READ,
01416             IOobject::AUTO_WRITE
01417         ),
01418         labelList(mesh.nCells(), -1)
01419     ),
01420     boundaryFacePoint_
01421     (
01422         IOobject
01423         (
01424             "boundaryFacePoint",
01425             time().findInstance(meshDir(), "faces"),
01426             meshSubDir,
01427             *this,
01428             IOobject::NO_READ,
01429             IOobject::AUTO_WRITE
01430         ),
01431         labelList(mesh.nFaces() - mesh.nInternalFaces())
01432     )
01433 {
01434     calcDual(mesh, featureEdges, featurePoints);
01435 }
01436 
01437 
01438 // Construct from polyMesh and feature angle
01439 Foam::polyDualMesh::polyDualMesh
01440 (
01441     const polyMesh& mesh,
01442     const scalar featureCos
01443 )
01444 :
01445     polyMesh
01446     (
01447         mesh,
01448         xferCopy(pointField()),   // to prevent any warnings "points not allocated"
01449         xferCopy(faceList()),     // to prevent any warnings "faces  not allocated"
01450         xferCopy(cellList())
01451     ),
01452     cellPoint_
01453     (
01454         IOobject
01455         (
01456             "cellPoint",
01457             time().findInstance(meshDir(), "faces"),
01458             meshSubDir,
01459             *this,
01460             IOobject::NO_READ,
01461             IOobject::AUTO_WRITE
01462         ),
01463         labelList(mesh.nCells(), -1)
01464     ),
01465     boundaryFacePoint_
01466     (
01467         IOobject
01468         (
01469             "boundaryFacePoint",
01470             time().findInstance(meshDir(), "faces"),
01471             meshSubDir,
01472             *this,
01473             IOobject::NO_READ,
01474             IOobject::AUTO_WRITE
01475         ),
01476         labelList(mesh.nFaces() - mesh.nInternalFaces(), -1)
01477     )
01478 {
01479     labelList featureEdges, featurePoints;
01480 
01481     calcFeatures(mesh, featureCos, featureEdges, featurePoints);
01482     calcDual(mesh, featureEdges, featurePoints);
01483 }
01484 
01485 
01486 void Foam::polyDualMesh::calcFeatures
01487 (
01488     const polyMesh& mesh,
01489     const scalar featureCos,
01490     labelList& featureEdges,
01491     labelList& featurePoints
01492 )
01493 {
01494     // Create big primitivePatch for all outside.
01495     primitivePatch allBoundary
01496     (
01497         SubList<face>
01498         (
01499             mesh.faces(),
01500             mesh.nFaces() - mesh.nInternalFaces(),
01501             mesh.nInternalFaces()
01502         ),
01503         mesh.points()
01504     );
01505 
01506     // For ease of use store patch number per face in allBoundary.
01507     labelList allRegion(allBoundary.size());
01508 
01509     const polyBoundaryMesh& patches = mesh.boundaryMesh();
01510 
01511     forAll(patches, patchI)
01512     {
01513         const polyPatch& pp = patches[patchI];
01514 
01515         forAll(pp, i)
01516         {
01517             allRegion[i + pp.start() - mesh.nInternalFaces()] = patchI;
01518         }
01519     }
01520 
01521 
01522     // Calculate patch/feature edges
01523     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01524 
01525     const labelListList& edgeFaces = allBoundary.edgeFaces();
01526     const vectorField& faceNormals = allBoundary.faceNormals();
01527     const labelList& meshPoints = allBoundary.meshPoints();
01528 
01529     boolList isFeatureEdge(edgeFaces.size(), false);
01530 
01531     forAll(edgeFaces, edgeI)
01532     {
01533         const labelList& eFaces = edgeFaces[edgeI];
01534 
01535         if (eFaces.size() != 2)
01536         {
01537             // Non-manifold. Problem?
01538             const edge& e = allBoundary.edges()[edgeI];
01539 
01540             WarningIn("polyDualMesh::calcFeatures") << "Edge "
01541                 << meshPoints[e[0]] << ' ' << meshPoints[e[1]]
01542                 << "  coords:" << mesh.points()[meshPoints[e[0]]]
01543                 << mesh.points()[meshPoints[e[1]]]
01544                 << " has more than 2 faces connected to it:"
01545                 << eFaces.size() << endl;
01546 
01547             isFeatureEdge[edgeI] = true;
01548         }
01549         else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
01550         {
01551             isFeatureEdge[edgeI] = true;
01552         }
01553         else if
01554         (
01555             (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
01556           < featureCos
01557         )
01558         {
01559             isFeatureEdge[edgeI] = true;
01560         }
01561     }
01562 
01563 
01564     // Calculate feature points
01565     // ~~~~~~~~~~~~~~~~~~~~~~~~
01566 
01567     const labelListList& pointEdges = allBoundary.pointEdges();
01568 
01569     DynamicList<label> allFeaturePoints(pointEdges.size());
01570 
01571     forAll(pointEdges, pointI)
01572     {
01573         const labelList& pEdges = pointEdges[pointI];
01574 
01575         label nFeatEdges = 0;
01576 
01577         forAll(pEdges, i)
01578         {
01579             if (isFeatureEdge[pEdges[i]])
01580             {
01581                 nFeatEdges++;
01582             }
01583         }
01584         if (nFeatEdges > 2)
01585         {
01586             // Store in mesh vertex label
01587             allFeaturePoints.append(allBoundary.meshPoints()[pointI]);
01588         }
01589     }
01590     featurePoints.transfer(allFeaturePoints);
01591 
01592     if (debug)
01593     {
01594         OFstream str("featurePoints.obj");
01595 
01596         Pout<< "polyDualMesh::calcFeatures : dumping feature points to "
01597             << str.name() << endl;
01598 
01599         forAll(featurePoints, i)
01600         {
01601             label pointI = featurePoints[i];
01602             meshTools::writeOBJ(str, mesh.points()[pointI]);
01603         }
01604     }
01605 
01606 
01607     // Get all feature edges.
01608     labelList meshEdges
01609     (
01610         allBoundary.meshEdges
01611         (
01612             mesh.edges(),
01613             mesh.cellEdges(),
01614             SubList<label>
01615             (
01616                 mesh.faceOwner(),
01617                 allBoundary.size(),
01618                 mesh.nInternalFaces()
01619             )
01620         )
01621     );
01622 
01623     DynamicList<label> allFeatureEdges(isFeatureEdge.size());
01624     forAll(isFeatureEdge, edgeI)
01625     {
01626         if (isFeatureEdge[edgeI])
01627         {
01628             // Store in mesh edge label.
01629             allFeatureEdges.append(meshEdges[edgeI]);
01630         }
01631     }
01632     featureEdges.transfer(allFeatureEdges);
01633 }
01634 
01635 
01636 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
01637 
01638 Foam::polyDualMesh::~polyDualMesh()
01639 {}
01640 
01641 
01642 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines