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

meshTools.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 "meshTools.H"
00027 #include <OpenFOAM/polyMesh.H>
00028 #include <OpenFOAM/hexMatcher.H>
00029 
00030 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00031 
00032 // Check if n is in same direction as normals of all faceLabels
00033 bool Foam::meshTools::visNormal
00034 (
00035     const vector& n,
00036     const vectorField& faceNormals,
00037     const labelList& faceLabels
00038 )
00039 {
00040     forAll(faceLabels, i)
00041     {
00042         if ((faceNormals[faceLabels[i]] & n) < SMALL)
00043         {
00044             // Found normal in different direction from n.
00045             return false;
00046         }
00047     }
00048 
00049     return true;
00050 }
00051 
00052 
00053 Foam::vectorField Foam::meshTools::calcBoxPointNormals(const primitivePatch& pp)
00054 {
00055     vectorField octantNormal(8);
00056     octantNormal[mXmYmZ] = vector(-1, -1, -1);
00057     octantNormal[pXmYmZ] = vector(1, -1, -1);
00058     octantNormal[mXpYmZ] = vector(-1, 1, -1);
00059     octantNormal[pXpYmZ] = vector(1, 1, -1);
00060 
00061     octantNormal[mXmYpZ] = vector(-1, -1, 1);
00062     octantNormal[pXmYpZ] = vector(1, -1, 1);
00063     octantNormal[mXpYpZ] = vector(-1, 1, 1);
00064     octantNormal[pXpYpZ] = vector(1, 1, 1);
00065 
00066     octantNormal /= mag(octantNormal);
00067 
00068 
00069     vectorField pn(pp.nPoints());
00070 
00071     const vectorField& faceNormals = pp.faceNormals();
00072     const vectorField& pointNormals = pp.pointNormals();
00073     const labelListList& pointFaces = pp.pointFaces();
00074 
00075     forAll(pointFaces, pointI)
00076     {
00077         const labelList& pFaces = pointFaces[pointI];
00078 
00079         if (visNormal(pointNormals[pointI], faceNormals, pFaces))
00080         {
00081             pn[pointI] = pointNormals[pointI];
00082         }
00083         else
00084         {
00085             WarningIn
00086             (
00087                 "Foam::meshTools::calcBoxPointNormals(const primitivePatch& pp)"
00088             )   << "Average point normal not visible for point:"
00089                 << pp.meshPoints()[pointI] << endl;
00090 
00091             label visOctant =
00092                 mXmYmZMask
00093               | pXmYmZMask
00094               | mXpYmZMask
00095               | pXpYmZMask
00096               | mXmYpZMask
00097               | pXmYpZMask
00098               | mXpYpZMask
00099               | pXpYpZMask;
00100 
00101             forAll(pFaces, i)
00102             {
00103                 const vector& n = faceNormals[pFaces[i]];
00104 
00105                 if (n.x() > SMALL)
00106                 {
00107                     // All -x octants become invisible
00108                     visOctant &= ~mXmYmZMask;
00109                     visOctant &= ~mXmYpZMask;
00110                     visOctant &= ~mXpYmZMask;
00111                     visOctant &= ~mXpYpZMask;
00112                 }
00113                 else if (n.x() < -SMALL)
00114                 {
00115                     // All +x octants become invisible
00116                     visOctant &= ~pXmYmZMask;
00117                     visOctant &= ~pXmYpZMask;
00118                     visOctant &= ~pXpYmZMask;
00119                     visOctant &= ~pXpYpZMask;
00120                 }
00121 
00122                 if (n.y() > SMALL)
00123                 {
00124                     visOctant &= ~mXmYmZMask;
00125                     visOctant &= ~mXmYpZMask;
00126                     visOctant &= ~pXmYmZMask;
00127                     visOctant &= ~pXmYpZMask;
00128                 }
00129                 else if (n.y() < -SMALL)
00130                 {
00131                     visOctant &= ~mXpYmZMask;
00132                     visOctant &= ~mXpYpZMask;
00133                     visOctant &= ~pXpYmZMask;
00134                     visOctant &= ~pXpYpZMask;
00135                 }
00136 
00137                 if (n.z() > SMALL)
00138                 {
00139                     visOctant &= ~mXmYmZMask;
00140                     visOctant &= ~mXpYmZMask;
00141                     visOctant &= ~pXmYmZMask;
00142                     visOctant &= ~pXpYmZMask;
00143                 }
00144                 else if (n.z() < -SMALL)
00145                 {
00146                     visOctant &= ~mXmYpZMask;
00147                     visOctant &= ~mXpYpZMask;
00148                     visOctant &= ~pXmYpZMask;
00149                     visOctant &= ~pXpYpZMask;
00150                 }
00151             }
00152 
00153             label visI = -1;
00154 
00155             label mask = 1;
00156 
00157             for (label octant = 0; octant < 8; octant++)
00158             {
00159                 if (visOctant & mask)
00160                 {
00161                     // Take first visible octant found. Note: could use
00162                     // combination of visible octants here.
00163                     visI = octant;
00164 
00165                     break;
00166                 }
00167                 mask <<= 1;
00168             }
00169 
00170             if (visI != -1)
00171             {
00172                 // Take a vector in this octant.
00173                 pn[pointI] = octantNormal[visI];
00174             }
00175             else
00176             {
00177                 pn[pointI] = vector::zero;
00178 
00179                 WarningIn
00180                 (
00181                     "Foam::meshTools::calcBoxPointNormals"
00182                     "(const primitivePatch& pp)"
00183                 )   << "No visible octant for point:" << pp.meshPoints()[pointI]
00184                     << " cooord:" << pp.points()[pp.meshPoints()[pointI]] << nl
00185                     << "Normal set to " << pn[pointI] << endl;
00186             }
00187         }
00188     }
00189 
00190     return pn;
00191 }
00192 
00193 
00194 Foam::vector Foam::meshTools::normEdgeVec
00195 (
00196     const primitiveMesh& mesh,
00197     const label edgeI
00198 )
00199 {
00200     vector eVec = mesh.edges()[edgeI].vec(mesh.points());
00201 
00202     eVec /= mag(eVec);
00203 
00204     return eVec;
00205 }
00206 
00207 
00208 void Foam::meshTools::writeOBJ
00209 (
00210     Ostream& os,
00211     const point& pt
00212 )
00213 {
00214     os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
00215 }
00216 
00217 
00218 void Foam::meshTools::writeOBJ
00219 (
00220     Ostream& os,
00221     const faceList& faces,
00222     const pointField& points,
00223     const labelList& faceLabels
00224 )
00225 {
00226     Map<label> foamToObj(4*faceLabels.size());
00227 
00228     label vertI = 0;
00229 
00230     forAll(faceLabels, i)
00231     {
00232         const face& f = faces[faceLabels[i]];
00233 
00234         forAll(f, fp)
00235         {
00236             if (foamToObj.insert(f[fp], vertI))
00237             {
00238                 writeOBJ(os, points[f[fp]]);
00239                 vertI++;
00240             }
00241         }
00242 
00243         os << 'l';
00244         forAll(f, fp)
00245         {
00246             os << ' ' << foamToObj[f[fp]]+1;
00247         }
00248         os << ' ' << foamToObj[f[0]]+1 << endl;
00249     }
00250 }
00251 
00252 
00253 void Foam::meshTools::writeOBJ
00254 (
00255     Ostream& os,
00256     const faceList& faces,
00257     const pointField& points
00258 )
00259 {
00260     labelList allFaces(faces.size());
00261     forAll(allFaces, i)
00262     {
00263         allFaces[i] = i;
00264     }
00265     writeOBJ(os, faces, points, allFaces);
00266 }
00267 
00268 
00269 void Foam::meshTools::writeOBJ
00270 (
00271     Ostream& os,
00272     const cellList& cells,
00273     const faceList& faces,
00274     const pointField& points,
00275     const labelList& cellLabels
00276 )
00277 {
00278     labelHashSet usedFaces(4*cellLabels.size());
00279 
00280     forAll(cellLabels, i)
00281     {
00282         const cell& cFaces = cells[cellLabels[i]];
00283 
00284         forAll(cFaces, j)
00285         {
00286             usedFaces.insert(cFaces[j]);
00287         }
00288     }
00289 
00290     writeOBJ(os, faces, points, usedFaces.toc());
00291 }
00292 
00293 
00294 bool Foam::meshTools::edgeOnCell
00295 (
00296     const primitiveMesh& mesh,
00297     const label cellI,
00298     const label edgeI
00299 )
00300 {
00301     return findIndex(mesh.edgeCells(edgeI), cellI) != -1;
00302 }
00303 
00304 
00305 bool Foam::meshTools::edgeOnFace
00306 (
00307     const primitiveMesh& mesh,
00308     const label faceI,
00309     const label edgeI
00310 )
00311 {
00312     return findIndex(mesh.faceEdges(faceI), edgeI) != -1;
00313 }
00314 
00315 
00316 // Return true if faceI part of cellI
00317 bool Foam::meshTools::faceOnCell
00318 (
00319     const primitiveMesh& mesh,
00320     const label cellI,
00321     const label faceI
00322 )
00323 {
00324     if (mesh.isInternalFace(faceI))
00325     {
00326         if
00327         (
00328             (mesh.faceOwner()[faceI] == cellI)
00329          || (mesh.faceNeighbour()[faceI] == cellI)
00330         )
00331         {
00332             return true;
00333         }
00334     }
00335     else
00336     {
00337         if (mesh.faceOwner()[faceI] == cellI)
00338         {
00339             return true;
00340         }
00341     }
00342     return false;
00343 }
00344 
00345 
00346 Foam::label Foam::meshTools::findEdge
00347 (
00348     const edgeList& edges,
00349     const labelList& candidates,
00350     const label v0,
00351     const label v1
00352 )
00353 {
00354     forAll(candidates, i)
00355     {
00356         label edgeI = candidates[i];
00357 
00358         const edge& e = edges[edgeI];
00359 
00360         if ((e[0] == v0 && e[1] == v1) || (e[0] == v1 && e[1] == v0))
00361         {
00362             return edgeI;
00363         }
00364     }
00365     return -1;
00366 }
00367 
00368 
00369 Foam::label Foam::meshTools::findEdge
00370 (
00371     const primitiveMesh& mesh,
00372     const label v0,
00373     const label v1
00374 )
00375 {
00376     const edgeList& edges = mesh.edges();
00377 
00378     const labelList& v0Edges = mesh.pointEdges()[v0];
00379 
00380     forAll(v0Edges, i)
00381     {
00382         label edgeI = v0Edges[i];
00383 
00384         const edge& e = edges[edgeI];
00385 
00386         if ((e.start() == v1) || (e.end() == v1))
00387         {
00388             return edgeI;
00389         }
00390     }
00391     return -1;
00392 }
00393 
00394 
00395 Foam::label Foam::meshTools::getSharedEdge
00396 (
00397     const primitiveMesh& mesh,
00398     const label f0,
00399     const label f1
00400 )
00401 {
00402     const labelList& f0Edges = mesh.faceEdges()[f0];
00403     const labelList& f1Edges = mesh.faceEdges()[f1];
00404 
00405     forAll(f0Edges, f0EdgeI)
00406     {
00407         label edge0 = f0Edges[f0EdgeI];
00408 
00409         forAll(f1Edges, f1EdgeI)
00410         {
00411             label edge1 = f1Edges[f1EdgeI];
00412 
00413             if (edge0 == edge1)
00414             {
00415                 return edge0;
00416             }
00417         }
00418     }
00419     FatalErrorIn
00420     (
00421         "meshTools::getSharedEdge(const primitiveMesh&, const label"
00422         ", const label)"
00423     )   << "Faces " << f0 << " and " << f1 << " do not share an edge"
00424         << abort(FatalError);
00425 
00426     return -1;
00427 
00428 }
00429 
00430 
00431 Foam::label Foam::meshTools::getSharedFace
00432 (
00433     const primitiveMesh& mesh,
00434     const label cell0I,
00435     const label cell1I
00436 )
00437 {
00438     const cell& cFaces = mesh.cells()[cell0I];
00439 
00440     forAll(cFaces, cFaceI)
00441     {
00442         label faceI = cFaces[cFaceI];
00443 
00444         if
00445         (
00446             mesh.isInternalFace(faceI)
00447          && (
00448                 mesh.faceOwner()[faceI] == cell1I
00449              || mesh.faceNeighbour()[faceI] == cell1I
00450             )
00451         )
00452         {
00453             return faceI;
00454         }
00455     }
00456 
00457 
00458     FatalErrorIn
00459     (
00460         "meshTools::getSharedFace(const primitiveMesh&, const label"
00461         ", const label)"
00462     )   << "No common face for"
00463         << "  cell0I:" << cell0I << "  faces:" << cFaces
00464         << "  cell1I:" << cell1I << "  faces:"
00465         << mesh.cells()[cell1I]
00466         << abort(FatalError);
00467 
00468     return -1;
00469 }
00470 
00471 
00472 // Get the two faces on cellI using edgeI.
00473 void Foam::meshTools::getEdgeFaces
00474 (
00475     const primitiveMesh& mesh,
00476     const label cellI,
00477     const label edgeI,
00478     label& face0,
00479     label& face1
00480 )
00481 {
00482     const labelList& eFaces = mesh.edgeFaces(edgeI);
00483 
00484     face0 = -1;
00485     face1 = -1;
00486 
00487     forAll(eFaces, eFaceI)
00488     {
00489         label faceI = eFaces[eFaceI];
00490 
00491         if (faceOnCell(mesh, cellI, faceI))
00492         {
00493             if (face0 == -1)
00494             {
00495                 face0 = faceI;
00496             }
00497             else
00498             {
00499                 face1 = faceI;
00500 
00501                 return;
00502             }
00503         }
00504     }
00505 
00506     if ((face0 == -1) || (face1 == -1))
00507     {
00508         FatalErrorIn
00509         (
00510             "meshTools::getEdgeFaces(const primitiveMesh&, const label"
00511             ", const label, label&, label&"
00512         )   << "Can not find faces using edge " << mesh.edges()[edgeI]
00513             << " on cell " << cellI << abort(FatalError);
00514     }
00515 }
00516 
00517 
00518 // Return label of other edge connected to vertex
00519 Foam::label Foam::meshTools::otherEdge
00520 (
00521     const primitiveMesh& mesh,
00522     const labelList& edgeLabels,
00523     const label thisEdgeI,
00524     const label thisVertI
00525 )
00526 {
00527     forAll(edgeLabels, edgeLabelI)
00528     {
00529         label edgeI = edgeLabels[edgeLabelI];
00530 
00531         if (edgeI != thisEdgeI)
00532         {
00533             const edge& e = mesh.edges()[edgeI];
00534 
00535             if ((e.start() == thisVertI) || (e.end() == thisVertI))
00536             {
00537                 return edgeI;
00538             }
00539         }
00540     }
00541 
00542     FatalErrorIn
00543     (
00544         "meshTools::otherEdge(const primitiveMesh&, const labelList&"
00545         ", const label, const label)"
00546     )   << "Can not find edge in "
00547         << UIndirectList<edge>(mesh.edges(), edgeLabels)()
00548         << " connected to edge "
00549         << thisEdgeI << " with vertices " << mesh.edges()[thisEdgeI]
00550         << " on side " << thisVertI << abort(FatalError);
00551 
00552     return -1;
00553 }
00554 
00555 
00556 // Return face on other side of edgeI
00557 Foam::label Foam::meshTools::otherFace
00558 (
00559     const primitiveMesh& mesh,
00560     const label cellI,
00561     const label faceI,
00562     const label edgeI
00563 )
00564 {
00565     label face0;
00566     label face1;
00567 
00568     getEdgeFaces(mesh, cellI, edgeI, face0, face1);
00569 
00570     if (face0 == faceI)
00571     {
00572         return face1;
00573     }
00574     else
00575     {
00576         return face0;
00577     }
00578 }
00579 
00580 
00581 // Return face on other side of edgeI
00582 Foam::label Foam::meshTools::otherCell
00583 (
00584     const primitiveMesh& mesh,
00585     const label otherCellI,
00586     const label faceI
00587 )
00588 {
00589     if (!mesh.isInternalFace(faceI))
00590     {
00591         FatalErrorIn
00592         (
00593             "meshTools::otherCell(const primitiveMesh&, const label"
00594             ", const label)"
00595         )   << "Face " << faceI << " is not internal"
00596             << abort(FatalError);
00597     }
00598 
00599     label newCellI = mesh.faceOwner()[faceI];
00600 
00601     if (newCellI == otherCellI)
00602     {
00603         newCellI = mesh.faceNeighbour()[faceI];
00604     }
00605     return newCellI;
00606 }
00607 
00608 
00609 // Returns label of edge nEdges away from startEdge (in the direction of
00610 // startVertI)
00611 Foam::label Foam::meshTools::walkFace
00612 (
00613     const primitiveMesh& mesh,
00614     const label faceI,
00615     const label startEdgeI,
00616     const label startVertI,
00617     const label nEdges
00618 )
00619 {
00620     const labelList& fEdges = mesh.faceEdges(faceI);
00621 
00622     label edgeI = startEdgeI;
00623 
00624     label vertI = startVertI;
00625 
00626     for (label iter = 0; iter < nEdges; iter++)
00627     {
00628         edgeI = otherEdge(mesh, fEdges, edgeI, vertI);
00629 
00630         vertI = mesh.edges()[edgeI].otherVertex(vertI);
00631     }
00632 
00633     return edgeI;
00634 }
00635 
00636 
00637 void Foam::meshTools::constrainToMeshCentre
00638 (
00639     const polyMesh& mesh,
00640     point& pt
00641 )
00642 {
00643     const Vector<label>& dirs = mesh.geometricD();
00644 
00645     const point& min = mesh.bounds().min();
00646     const point& max = mesh.bounds().max();
00647 
00648     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00649     {
00650         if (dirs[cmpt] == -1)
00651         {
00652             pt[cmpt] = 0.5*(min[cmpt] + max[cmpt]);
00653         }
00654     }
00655 }
00656 
00657 
00658 void Foam::meshTools::constrainToMeshCentre
00659 (
00660     const polyMesh& mesh,
00661     pointField& pts
00662 )
00663 {
00664     const Vector<label>& dirs = mesh.geometricD();
00665 
00666     const point& min = mesh.bounds().min();
00667     const point& max = mesh.bounds().max();
00668 
00669     bool isConstrained = false;
00670     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00671     {
00672         if (dirs[cmpt] == -1)
00673         {
00674             isConstrained = true;
00675             break;
00676         }
00677     }
00678 
00679     if (isConstrained)
00680     {
00681         forAll(pts, i)
00682         {
00683             for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00684             {
00685                 if (dirs[cmpt] == -1)
00686                 {
00687                     pts[i][cmpt] = 0.5*(min[cmpt] + max[cmpt]);
00688                 }
00689             }
00690         }
00691     }
00692 }
00693 
00694 
00695 //- Set the constrained components of directions/velocity to zero
00696 void Foam::meshTools::constrainDirection
00697 (
00698     const polyMesh& mesh,
00699     const Vector<label>& dirs,
00700     vector& d
00701 )
00702 {
00703     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00704     {
00705         if (dirs[cmpt] == -1)
00706         {
00707             d[cmpt] = 0.0;
00708         }
00709     }
00710 }
00711 
00712 
00713 void Foam::meshTools::constrainDirection
00714 (
00715     const polyMesh& mesh,
00716     const Vector<label>& dirs,
00717     vectorField& d
00718 )
00719 {
00720     bool isConstrained = false;
00721     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00722     {
00723         if (dirs[cmpt] == -1)
00724         {
00725             isConstrained = true;
00726             break;
00727         }
00728     }
00729 
00730     if (isConstrained)
00731     {
00732         forAll(d, i)
00733         {
00734             for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
00735             {
00736                 if (dirs[cmpt] == -1)
00737                 {
00738                     d[i][cmpt] = 0.0;
00739                 }
00740             }
00741         }
00742     }
00743 }
00744 
00745 
00746 void Foam::meshTools::getParallelEdges
00747 (
00748     const primitiveMesh& mesh,
00749     const label cellI,
00750     const label e0,
00751     label& e1,
00752     label& e2,
00753     label& e3
00754 )
00755 {
00756     // Go to any face using e0
00757     label faceI = meshTools::otherFace(mesh, cellI, -1, e0);
00758 
00759     // Opposite edge on face
00760     e1 = meshTools::walkFace(mesh, faceI, e0, mesh.edges()[e0].end(), 2);
00761 
00762     faceI = meshTools::otherFace(mesh, cellI, faceI, e1);
00763 
00764     e2 = meshTools::walkFace(mesh, faceI, e1, mesh.edges()[e1].end(), 2);
00765 
00766     faceI = meshTools::otherFace(mesh, cellI, faceI, e2);
00767 
00768     e3 = meshTools::walkFace(mesh, faceI, e2, mesh.edges()[e2].end(), 2);
00769 }
00770 
00771 
00772 Foam::vector Foam::meshTools::edgeToCutDir
00773 (
00774     const primitiveMesh& mesh,
00775     const label cellI,
00776     const label startEdgeI
00777 )
00778 {
00779     if (!hexMatcher().isA(mesh, cellI))
00780     {
00781         FatalErrorIn
00782         (
00783             "Foam::meshTools::getCutDir(const label, const label)"
00784         )   << "Not a hex : cell:" << cellI << abort(FatalError);
00785     }
00786 
00787 
00788     vector avgVec(normEdgeVec(mesh, startEdgeI));
00789 
00790     label edgeI = startEdgeI;
00791 
00792     label faceI = -1;
00793 
00794     for (label i = 0; i < 3; i++)
00795     {
00796         // Step to next face, next edge
00797         faceI = meshTools::otherFace(mesh, cellI, faceI, edgeI);
00798 
00799         vector eVec(normEdgeVec(mesh, edgeI));
00800 
00801         if ((eVec & avgVec) > 0)
00802         {
00803             avgVec += eVec;
00804         }
00805         else
00806         {
00807             avgVec -= eVec;
00808         }
00809 
00810         label vertI = mesh.edges()[edgeI].end();
00811 
00812         edgeI = meshTools::walkFace(mesh, faceI, edgeI, vertI, 2);
00813     }
00814 
00815     avgVec /= mag(avgVec) + VSMALL;
00816 
00817     return avgVec;
00818 }
00819 
00820 
00821 // Find edges most aligned with cutDir
00822 Foam::label Foam::meshTools::cutDirToEdge
00823 (
00824     const primitiveMesh& mesh,
00825     const label cellI,
00826     const vector& cutDir
00827 )
00828 {
00829     if (!hexMatcher().isA(mesh, cellI))
00830     {
00831         FatalErrorIn
00832         (
00833             "Foam::meshTools::getCutDir(const label, const vector&)"
00834         )   << "Not a hex : cell:" << cellI << abort(FatalError);
00835     }
00836 
00837     const labelList& cEdges = mesh.cellEdges()[cellI];
00838 
00839     labelHashSet doneEdges(2*cEdges.size());
00840 
00841     scalar maxCos = -GREAT;
00842     label maxEdgeI = -1;
00843 
00844     for (label i = 0; i < 4; i++)
00845     {
00846         forAll(cEdges, cEdgeI)
00847         {
00848             label e0 = cEdges[cEdgeI];
00849 
00850             if (!doneEdges.found(e0))
00851             {
00852                 vector avgDir(edgeToCutDir(mesh, cellI, e0));
00853 
00854                 scalar cosAngle = mag(avgDir & cutDir);
00855 
00856                 if (cosAngle > maxCos)
00857                 {
00858                     maxCos = cosAngle;
00859                     maxEdgeI = e0;
00860                 }
00861 
00862                 // Mark off edges in cEdges.
00863                 label e1, e2, e3;
00864                 getParallelEdges(mesh, cellI, e0, e1, e2, e3);
00865 
00866                 doneEdges.insert(e0);
00867                 doneEdges.insert(e1);
00868                 doneEdges.insert(e2);
00869                 doneEdges.insert(e3);
00870             }
00871         }
00872     }
00873 
00874     forAll(cEdges, cEdgeI)
00875     {
00876         if (!doneEdges.found(cEdges[cEdgeI]))
00877         {
00878             FatalErrorIn
00879             (
00880                 "meshTools::cutDirToEdge(const label, const vector&)"
00881             )   << "Cell:" << cellI << " edges:" << cEdges << endl
00882                 << "Edge:" << cEdges[cEdgeI] << " not yet handled"
00883                 << abort(FatalError);
00884         }
00885     }
00886 
00887     if (maxEdgeI == -1)
00888     {
00889         FatalErrorIn
00890         (
00891             "meshTools::cutDirToEdge(const label, const vector&)"
00892         )   << "Problem : did not find edge aligned with " << cutDir
00893             << " on cell " << cellI << abort(FatalError);
00894     }
00895 
00896     return maxEdgeI;
00897 }
00898 
00899 
00900 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines