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

topoCellLooper.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 "topoCellLooper.H"
00027 #include <meshTools/cellFeatures.H>
00028 #include <OpenFOAM/polyMesh.H>
00029 #include <OpenFOAM/mathematicalConstants.H>
00030 #include <OpenFOAM/DynamicList.H>
00031 #include <OpenFOAM/ListOps.H>
00032 #include <meshTools/meshTools.H>
00033 #include <OpenFOAM/hexMatcher.H>
00034 
00035 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00036 
00037 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00038 
00039 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00040 namespace Foam
00041 {
00042    defineTypeNameAndDebug(topoCellLooper, 0);
00043    addToRunTimeSelectionTable(cellLooper, topoCellLooper, word);
00044 }
00045 
00046 // Angle for polys to be considered splitHexes.
00047 const Foam::scalar Foam::topoCellLooper::featureCos =
00048     Foam::cos(10.0 * mathematicalConstant::pi/180.0);
00049 
00050 
00051 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00052 
00053 // In-memory truncate a list
00054 template <class T>
00055 void Foam::topoCellLooper::subsetList
00056 (
00057     const label startI,
00058     const label freeI,
00059     DynamicList<T>& lst
00060 )
00061 {
00062     if (startI == 0)
00063     {
00064         // Truncate (setSize decides itself not to do anything if nothing
00065         // changed)
00066         if (freeI < 0)
00067         {
00068             FatalErrorIn("topoCellLooper::subsetList")
00069                 << "startI:" << startI << "  freeI:" << freeI
00070                 << "  lst:" << lst << abort(FatalError);
00071         }
00072         lst.setCapacity(freeI);
00073     }
00074     else
00075     {
00076         // Shift elements down
00077         label newI = 0;
00078         for (label elemI = startI; elemI < freeI; elemI++)
00079         {
00080             lst[newI++] = lst[elemI];
00081         }
00082 
00083         if ((freeI - startI) < 0)
00084         {
00085             FatalErrorIn("topoCellLooper::subsetList")
00086                 << "startI:" << startI << "  freeI:" << freeI
00087                 << "  lst:" << lst << abort(FatalError);
00088         }
00089 
00090         lst.setCapacity(freeI - startI);
00091     }
00092 }
00093 
00094 
00095 // Returns label of edge nFeaturePts away from startEdge (in the direction of
00096 // startVertI) and not counting non-featurePoints.
00097 //
00098 // When stepping to this face it can happen in 3 different ways:
00099 //
00100 //  --|------
00101 //    |
00102 //  1 |   0
00103 //    |A
00104 //    |
00105 //    |
00106 //  --|------
00107 //
00108 //  A: jumping from face0 to face1 across edge A.
00109 //      startEdge != -1
00110 //      startVert == -1
00111 //
00112 //  --|------
00113 //    |
00114 //  1 |   0
00115 //    +B
00116 //    |
00117 //    |
00118 //  --|------
00119 //
00120 //  B: jumping from face0 to face1 across (non-feature) point B
00121 //      startEdge == -1
00122 //      startVert != -1
00123 //
00124 //  --|------
00125 //  0 |   1
00126 //    |C
00127 //  --+
00128 //    |
00129 //    |
00130 //  --|------
00131 //
00132 //  C: jumping from face0 to face1 across (feature) point C on edge.
00133 //      startEdge != -1
00134 //      startVert != -1
00135 //
00136 void Foam::topoCellLooper::walkFace
00137 (
00138     const cellFeatures& features,
00139     const label faceI,
00140     const label startEdgeI,
00141     const label startVertI,
00142     const label nFeaturePts,
00143 
00144     label& edgeI,
00145     label& vertI
00146 ) const
00147 {
00148     const labelList& fEdges = mesh().faceEdges()[faceI];
00149 
00150     edgeI = startEdgeI;
00151 
00152     vertI = startVertI;
00153 
00154     // Number of feature points crossed so far
00155     label nVisited = 0;
00156 
00157     if (vertI == -1)
00158     {
00159         // Started on edge. Go to one of its endpoints.
00160         vertI = mesh().edges()[edgeI].start();
00161 
00162         if (features.isFeatureVertex(faceI, vertI))
00163         {
00164             nVisited++;
00165         }
00166     }
00167 
00168     if ((edgeI == -1) || !meshTools::edgeOnFace(mesh(), faceI, edgeI))
00169     {
00170         // Either edge is not set or not on current face.  Just take one of
00171         // the edges on this face as starting edge.
00172         edgeI = getFirstVertEdge(faceI, vertI);
00173     }
00174 
00175     // Now we should have starting edge on face and a vertex on that edge.
00176 
00177     do
00178     {
00179 
00180         edgeI = meshTools::otherEdge(mesh(), fEdges, edgeI, vertI);
00181 
00182         if (nVisited == nFeaturePts)
00183         {
00184             break;
00185         }
00186 
00187         vertI = mesh().edges()[edgeI].otherVertex(vertI);
00188 
00189         if (features.isFeatureVertex(faceI, vertI))
00190         {
00191             nVisited++;
00192         }
00193     }
00194     while (true);
00195 }
00196 
00197 
00198 // Returns list of vertices on 'superEdge' i.e. list of edges connected by
00199 // non-feature points. First and last are feature points, ones inbetween are
00200 // not.
00201 Foam::labelList Foam::topoCellLooper::getSuperEdge
00202 (
00203     const cellFeatures& features,
00204     const label faceI,
00205     const label startEdgeI,
00206     const label startVertI
00207 ) const
00208 {
00209     const labelList& fEdges = mesh().faceEdges()[faceI];
00210 
00211     labelList superVerts(fEdges.size());
00212     label superVertI = 0;
00213 
00214 
00215     label edgeI = startEdgeI;
00216 
00217     label vertI = startVertI;
00218 
00219     superVerts[superVertI++] = vertI;
00220 
00221     label prevEdgeI = -1;
00222 
00223     do
00224     {
00225         vertI = mesh().edges()[edgeI].otherVertex(vertI);
00226 
00227         superVerts[superVertI++] = vertI;
00228 
00229         prevEdgeI = edgeI;
00230 
00231         edgeI = meshTools::otherEdge(mesh(), fEdges, edgeI, vertI);
00232     }
00233     while (!features.isFeaturePoint(prevEdgeI, edgeI));
00234 
00235     superVerts.setSize(superVertI);
00236 
00237     return superVerts;
00238 }
00239 
00240 
00241 // Return non-feature edge from cells' edges that is most perpendicular
00242 // to refinement direction.
00243 Foam::label Foam::topoCellLooper::getAlignedNonFeatureEdge
00244 (
00245     const vector& refDir,
00246     const label cellI,
00247     const cellFeatures& features
00248 ) const
00249 {
00250     const labelList& cEdges = mesh().cellEdges()[cellI];
00251 
00252     const point& ctr = mesh().cellCentres()[cellI];
00253 
00254     label cutEdgeI = -1;
00255     scalar maxCos = -GREAT;
00256 
00257     forAll(cEdges, cEdgeI)
00258     {
00259         label edgeI = cEdges[cEdgeI];
00260 
00261         if (!features.isFeatureEdge(edgeI))
00262         {
00263             const edge& e = mesh().edges()[edgeI];
00264 
00265             // Get plane spanned by e.start, e.end and cell centre.
00266             vector e0 = mesh().points()[e.start()] - ctr;
00267             vector e1 = mesh().points()[e.end()] - ctr;
00268 
00269             vector n = e0 ^ e1;
00270             n /= mag(n);
00271 
00272             scalar cosAngle = mag(refDir & n);
00273 
00274             if (cosAngle > maxCos)
00275             {
00276                 maxCos = cosAngle;
00277 
00278                 cutEdgeI = edgeI;
00279             }
00280         }
00281     }
00282 
00283     return cutEdgeI;
00284 }
00285 
00286 
00287 // Starts from edge and vertex on edge on face (or neighbouring face)
00288 // and steps either to existing vertex (vertI != -1) or to edge (vertI == -1)
00289 // by walking point-edge and crossing nFeats featurePoints.
00290 void Foam::topoCellLooper::walkAcrossFace
00291 (
00292     const cellFeatures& features,
00293     const label faceI,
00294     const label startEdgeI,
00295     const label startVertI,
00296     const label nFeats,
00297 
00298     label& edgeI,
00299     label& vertI
00300 ) const
00301 {
00302     label oppositeVertI = -1;
00303     label oppositeEdgeI = -1;
00304 
00305     // Go to oppositeEdge and a vertex on that.
00306     walkFace
00307     (
00308         features,
00309         faceI,
00310         startEdgeI,
00311         startVertI,
00312         nFeats,
00313 
00314         oppositeEdgeI,
00315         oppositeVertI
00316     );
00317 
00318     // Loop over super edge to find internal points if there are any.
00319 
00320     labelList superEdge =
00321         getSuperEdge
00322         (
00323             features,
00324             faceI,
00325             oppositeEdgeI,
00326             oppositeVertI
00327         );
00328 
00329     label sz = superEdge.size();
00330 
00331     if (sz == 2)
00332     {
00333         // No non-feature point inbetween feature points.
00334         // Mark edge.
00335 
00336         vertI = -1;
00337         edgeI = oppositeEdgeI;
00338     }
00339     else if (sz == 3)
00340     {
00341         vertI = superEdge[1];
00342         edgeI = -1;
00343     }
00344     else
00345     {
00346         //Should choose acc. to geometry!
00347         label index = sz/2;
00348 
00349         if (debug)
00350         {
00351             Pout<< "    Don't know what to do. Stepped to non-feature point "
00352                 << "at index " << index << " in superEdge:" << superEdge
00353                 << endl;
00354         }
00355 
00356         vertI = superEdge[index];
00357         edgeI = -1;
00358     }
00359 }
00360 
00361 
00362 // Walks cell circumference. Updates face, edge, vertex.
00363 //
00364 // Position on face is given by:
00365 //
00366 //  vertI == -1, faceI != -1, edgeI != -1
00367 //      on edge of face. Cross edge to neighbouring face.
00368 //
00369 //  vertI != -1, edgeI != -1, faceI == -1
00370 //      coming from edge onto vertex vertI. Need to step to one
00371 //      of the faces not using edgeI.
00372 //
00373 //  vertI != -1, edgeI == -1, faceI != -1
00374 //      coming from vertex on side of face. Step to one of the faces
00375 //      using vertI but not faceI
00376 void Foam::topoCellLooper::walkSplitHex
00377 (
00378     const label cellI,
00379     const cellFeatures& features,
00380     const label fromFaceI,
00381     const label fromEdgeI,
00382     const label fromVertI,
00383 
00384     DynamicList<label>& loop,
00385     DynamicList<scalar>& loopWeights
00386 ) const
00387 {
00388     // Work vars giving position on cell
00389     label faceI = fromFaceI;
00390     label edgeI = fromEdgeI;
00391     label vertI = fromVertI;
00392 
00393     do
00394     {
00395         if (debug)
00396         {
00397             Pout<< "Entering walk with : cell:" << cellI << " face:" << faceI;
00398             if (faceI != -1)
00399             {
00400                 Pout<< " verts:" << mesh().faces()[faceI];
00401             }
00402             Pout<< " edge:" << edgeI;
00403             if (edgeI != -1)
00404             {
00405                 Pout<< " verts:" << mesh().edges()[edgeI];
00406             }
00407             Pout<< " vert:" << vertI << endl;
00408         }
00409 
00410         label startLoop = -1;
00411 
00412         if
00413         (
00414             (vertI != -1)
00415          && (
00416                 (startLoop =
00417                     findIndex
00418                     (
00419                         loop,
00420                         vertToEVert(vertI)
00421                     )
00422                 )
00423             != -1
00424             )
00425         )
00426         {
00427             // Breaking walk since vertI already cut
00428             label firstFree = loop.size();
00429 
00430             subsetList(startLoop, firstFree, loop);
00431             subsetList(startLoop, firstFree, loopWeights);
00432 
00433             break;
00434         }
00435         if
00436         (
00437             (edgeI != -1)
00438          && (
00439                 (startLoop =
00440                     findIndex
00441                     (
00442                         loop,
00443                         edgeToEVert(edgeI)
00444                     )
00445                 )
00446              != -1
00447             )
00448         )
00449         {
00450             // Breaking walk since edgeI already cut
00451             label firstFree = loop.size();
00452 
00453             subsetList(startLoop, firstFree, loop);
00454             subsetList(startLoop, firstFree, loopWeights);
00455 
00456             break;
00457         }
00458 
00459 
00460         if (vertI == -1)
00461         {
00462             // On edge
00463             if (edgeI == -1)
00464             {
00465                 FatalErrorIn("walkSplitHex") << "Illegal edge and vert"
00466                     << abort(FatalError);
00467             }
00468 
00469             loop.append(edgeToEVert(edgeI));
00470             loopWeights.append(0.5);
00471 
00472             // Cross edge to next face
00473             faceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);
00474 
00475             if (debug)
00476             {
00477                 Pout<< "    stepped across edge " << mesh().edges()[edgeI]
00478                     << " to face " << faceI << " verts:"
00479                     << mesh().faces()[faceI] << endl;
00480             }
00481 
00482             label nextEdgeI = -1;
00483             label nextVertI = -1;
00484 
00485             // Walk face along its edges
00486             walkAcrossFace
00487             (
00488                 features,
00489                 faceI,
00490                 edgeI,
00491                 vertI,
00492                 2,
00493 
00494                 nextEdgeI,
00495                 nextVertI
00496             );
00497 
00498             edgeI = nextEdgeI;
00499             vertI = nextVertI;
00500         }
00501         else
00502         {
00503             // On vertex.
00504 
00505             loop.append(vertToEVert(vertI));
00506             loopWeights.append(-GREAT);
00507 
00508             if (edgeI == -1)
00509             {
00510                 // Normal vertex on edge of face. Get edges connected to it
00511                 // which are not on faceI.
00512                 labelList nextEdges = getVertEdgesNonFace
00513                 (
00514                     cellI,
00515                     faceI,
00516                     vertI
00517                 );
00518 
00519                 if (nextEdges.empty())
00520                 {
00521                     // Cross to other face (there is only one since no edges)
00522                     const labelList& pFaces = mesh().pointFaces()[vertI];
00523 
00524                     forAll(pFaces, pFaceI)
00525                     {
00526                         label thisFaceI = pFaces[pFaceI];
00527 
00528                         if
00529                         (
00530                             (thisFaceI != faceI)
00531                          && meshTools::faceOnCell(mesh(), cellI, thisFaceI)
00532                         )
00533                         {
00534                             faceI = thisFaceI;
00535                             break;
00536                         }
00537                     }
00538 
00539                     if (debug)
00540                     {
00541                         Pout<< "    stepped from non-edge vertex " << vertI
00542                             << " to face " << faceI << " verts:"
00543                             << mesh().faces()[faceI]
00544                             << " since candidate edges:" << nextEdges << endl;
00545                     }
00546 
00547                     label nextEdgeI = -1;
00548                     label nextVertI = -1;
00549 
00550                     walkAcrossFace
00551                     (
00552                         features,
00553                         faceI,
00554                         edgeI,
00555                         vertI,
00556                         2,          // 2 vertices to cross
00557 
00558                         nextEdgeI,
00559                         nextVertI
00560                     );
00561 
00562                     edgeI = nextEdgeI;
00563                     vertI = nextVertI;
00564                 }
00565                 else if (nextEdges.size() == 1)
00566                 {
00567                     // One edge. Go along this one.
00568                     edgeI = nextEdges[0];
00569 
00570                     if (debug)
00571                     {
00572                         Pout<< "    stepped from non-edge vertex " << vertI
00573                             << " along edge " << edgeI << " verts:"
00574                             << mesh().edges()[edgeI]
00575                             << " out of candidate edges:"
00576                             << nextEdges << endl;
00577                     }
00578 
00579                     vertI = mesh().edges()[edgeI].otherVertex(vertI);
00580 
00581                     faceI = -1;
00582                 }
00583                 else
00584                 {
00585                     // Multiple edges to choose from. Pick any one.
00586                     // (ideally should be geometric)
00587 
00588                     label index = nextEdges.size()/2;
00589 
00590                     edgeI = nextEdges[index];
00591 
00592                     if (debug)
00593                     {
00594                         Pout<< "    stepped from non-edge vertex " << vertI
00595                             << " along edge " << edgeI << " verts:"
00596                             << mesh().edges()[edgeI]
00597                             << " out of candidate edges:" << nextEdges << endl;
00598                     }
00599 
00600                     vertI = mesh().edges()[edgeI].otherVertex(vertI);
00601 
00602                     faceI = -1;
00603                 }
00604             }
00605             else
00606             {
00607                 // Get faces connected to startVertI but not startEdgeI
00608                 labelList nextFaces =
00609                     getVertFacesNonEdge
00610                     (
00611                         cellI,
00612                         edgeI,
00613                         vertI
00614                     );
00615 
00616                 if (nextFaces.size() == 1)
00617                 {
00618                     // Only one face to cross.
00619                     faceI = nextFaces[0];
00620 
00621                     label nextEdgeI = -1;
00622                     label nextVertI = -1;
00623 
00624                     walkAcrossFace
00625                     (
00626                         features,
00627                         faceI,
00628                         edgeI,
00629                         vertI,
00630                         2,          // 2 vertices to cross
00631 
00632                         nextEdgeI,
00633                         nextVertI
00634                     );
00635 
00636                     edgeI = nextEdgeI;
00637                     vertI = nextVertI;
00638                 }
00639                 else if (nextFaces.size() == 2)
00640                 {
00641                     // Split face. Get edge inbetween.
00642                     faceI = -1;
00643 
00644                     edgeI =
00645                         meshTools::getSharedEdge
00646                         (
00647                             mesh(),
00648                             nextFaces[0],
00649                             nextFaces[1]
00650                         );
00651 
00652                     vertI = mesh().edges()[edgeI].otherVertex(vertI);
00653                 }
00654                 else
00655                 {
00656                     FatalErrorIn("walkFromVert") << "Not yet implemented"
00657                         << "Choosing from more than "
00658                         << "two candidates:" << nextFaces
00659                         << " when coming from vertex " << vertI << " on cell "
00660                         << cellI << abort(FatalError);
00661                 }
00662             }
00663         }
00664 
00665         if (debug)
00666         {
00667             Pout<< "Walked to : face:" << faceI;
00668             if (faceI != -1)
00669             {
00670                 Pout<< " verts:" << mesh().faces()[faceI];
00671             }
00672             Pout<< " edge:" << edgeI;
00673             if (edgeI != -1)
00674             {
00675                 Pout<< " verts:" << mesh().edges()[edgeI];
00676             }
00677             Pout<< " vert:" << vertI << endl;
00678         }
00679     }
00680     while (true);
00681 }
00682 
00683 
00684 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00685 
00686 // Construct from components
00687 Foam::topoCellLooper::topoCellLooper(const polyMesh& mesh)
00688 :
00689     hexCellLooper(mesh)
00690 {}
00691 
00692 
00693 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00694 
00695 Foam::topoCellLooper::~topoCellLooper()
00696 {}
00697 
00698 
00699 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00700 
00701 bool Foam::topoCellLooper::cut
00702 (
00703     const vector& refDir,
00704     const label cellI,
00705     const boolList& vertIsCut,
00706     const boolList& edgeIsCut,
00707     const scalarField& edgeWeight,
00708 
00709     labelList& loop,
00710     scalarField& loopWeights
00711 ) const
00712 {
00713     if (mesh().cellShapes()[cellI].model() == hex_)
00714     {
00715         // Let parent handle hex case.
00716         return
00717             hexCellLooper::cut
00718             (
00719                 refDir,
00720                 cellI,
00721                 vertIsCut,
00722                 edgeIsCut,
00723                 edgeWeight,
00724                 loop,
00725                 loopWeights
00726             );
00727     }
00728     else
00729     {
00730         cellFeatures superCell(mesh(), featureCos, cellI);
00731 
00732         if (hexMatcher().isA(superCell.faces()))
00733         {
00734             label edgeI =
00735                 getAlignedNonFeatureEdge
00736                 (
00737                     refDir,
00738                     cellI,
00739                     superCell
00740                 );
00741 
00742             label vertI = -1;
00743 
00744             label faceI = -1;
00745 
00746             if (edgeI != -1)
00747             {
00748                 // Found non-feature edge. Start walking from vertex on edge.
00749                 vertI = mesh().edges()[edgeI].start();
00750             }
00751             else
00752             {
00753                 // No 'matching' non-feature edge found on cell. Get starting
00754                 // normal i.e. feature edge.
00755                 edgeI = getMisAlignedEdge(refDir, cellI);
00756 
00757                 // Get any face using edge
00758                 label face0;
00759                 label face1;
00760                 meshTools::getEdgeFaces(mesh(), cellI, edgeI, face0, face1);
00761 
00762                 faceI = face0;
00763             }
00764 
00765             label nEstCuts = 2*mesh().cells()[cellI].size();
00766 
00767             DynamicList<label> localLoop(nEstCuts);
00768             DynamicList<scalar> localLoopWeights(nEstCuts);
00769 
00770             walkSplitHex
00771             (
00772                 cellI,
00773                 superCell,
00774                 faceI,
00775                 edgeI,
00776                 vertI,
00777 
00778                 localLoop,
00779                 localLoopWeights
00780             );
00781 
00782             if (localLoop.size() <=2)
00783             {
00784                 return false;
00785             }
00786             else
00787             {
00788                 loop.transfer(localLoop);
00789                 loopWeights.transfer(localLoopWeights);
00790 
00791                 return true;
00792             }
00793         }
00794         else
00795         {
00796             // Let parent handle poly case.
00797             return hexCellLooper::cut
00798             (
00799                 refDir,
00800                 cellI,
00801                 vertIsCut,
00802                 edgeIsCut,
00803                 edgeWeight,
00804                 loop,
00805                 loopWeights
00806             );
00807         }
00808     }
00809 }
00810 
00811 
00812 bool Foam::topoCellLooper::cut
00813 (
00814     const plane& cutPlane,
00815     const label cellI,
00816     const boolList& vertIsCut,
00817     const boolList& edgeIsCut,
00818     const scalarField& edgeWeight,
00819 
00820     labelList& loop,
00821     scalarField& loopWeights
00822 ) const
00823 {
00824     // Let parent handle cut with plane.
00825     return
00826         hexCellLooper::cut
00827         (
00828             cutPlane,
00829             cellI,
00830             vertIsCut,
00831             edgeIsCut,
00832             edgeWeight,
00833 
00834             loop,
00835             loopWeights
00836         );
00837 }
00838 
00839 
00840 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines