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

directionInfo.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 
00027 #include "directionInfo.H"
00028 #include <OpenFOAM/hexMatcher.H>
00029 #include <meshTools/meshTools.H>
00030 #include <OpenFOAM/polyMesh.H>
00031 
00032 
00033 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00034 
00035 bool Foam::directionInfo::equal(const edge& e, const label v0, const label v1)
00036 {
00037     return
00038         (e.start() == v0 && e.end() == v1)
00039      || (e.start() == v1 && e.end() == v0);
00040 }
00041 
00042 
00043 Foam::point Foam::directionInfo::eMid
00044 (
00045     const primitiveMesh& mesh,
00046     const label edgeI
00047 )
00048 {
00049     const edge& e = mesh.edges()[edgeI];
00050 
00051     return
00052         0.5
00053       * (mesh.points()[e.start()] + mesh.points()[e.end()]);
00054 }
00055 
00056 
00057 // Find edge among edgeLabels that uses v0 and v1
00058 Foam::label Foam::directionInfo::findEdge
00059 (
00060     const primitiveMesh& mesh,
00061     const labelList& edgeLabels,
00062     const label v1,
00063     const label v0
00064 )
00065 {
00066     forAll(edgeLabels, edgeLabelI)
00067     {
00068         label edgeI = edgeLabels[edgeLabelI];
00069 
00070         const edge& e = mesh.edges()[edgeI];
00071 
00072         if (equal(e, v0, v1))
00073         {
00074             return edgeI;
00075         }
00076     }
00077 
00078     FatalErrorIn("directionInfo::findEdge")
00079         << "Cannot find an edge among " << edgeLabels << endl
00080         << "that uses vertices " << v0
00081         << " and " << v1
00082         << abort(FatalError);
00083 
00084     return -1;
00085 }
00086 
00087 
00088 Foam::label Foam::directionInfo::lowest
00089 (
00090     const label size,
00091     const label a,
00092     const label b
00093 )
00094 {
00095     // Get next point
00096     label a1 = (a + 1) % size;
00097 
00098     if (a1 == b)
00099     {
00100         return a;
00101     }
00102     else
00103     {
00104         label b1 = (b + 1) % size;
00105 
00106         if (b1 != a)
00107         {
00108             FatalErrorIn("directionInfo::lowest")
00109                 << "Problem : a:" << a << " b:" << b << " size:" << size
00110                 << abort(FatalError);
00111         }
00112 
00113         return b;
00114     }
00115 }
00116 
00117 
00118 // Have edge on hex cell. Find corresponding edge on face. Return -1 if none
00119 // found.
00120 Foam::label Foam::directionInfo::edgeToFaceIndex
00121 (
00122     const primitiveMesh& mesh,
00123     const label cellI,
00124     const label faceI,
00125     const label edgeI
00126 )
00127 {
00128     if ((edgeI < 0) || (edgeI >= mesh.nEdges()))
00129     {
00130         FatalErrorIn("directionInfo::edgeToFaceIndex")
00131             << "Illegal edge label:" << edgeI
00132             << " when projecting cut edge from cell " << cellI
00133             << " to face " << faceI
00134             << abort(FatalError);
00135     }
00136 
00137     const edge& e = mesh.edges()[edgeI];
00138 
00139     const face& f = mesh.faces()[faceI];
00140 
00141     // edgeI is either
00142     // - in faceI. Convert into index in face.
00143     // - connected (but not in) to face. Return -1.
00144     // - in face opposite faceI. Convert into index in face.
00145 
00146     label fpA = findIndex(f, e.start());
00147     label fpB = findIndex(f, e.end());
00148 
00149     if (fpA != -1)
00150     {
00151         if (fpB != -1)
00152         {
00153             return lowest(f.size(), fpA, fpB);
00154         }
00155         else
00156         {
00157             // e.start() in face, e.end() not
00158             return -1;
00159         }
00160     }
00161     else
00162     {
00163         if (fpB != -1)
00164         {
00165             // e.end() in face, e.start() not
00166             return -1;
00167         }
00168         else
00169         {
00170             // Both not in face.
00171             // e is on opposite face. Determine corresponding edge on this face:
00172             // - determine two faces using edge (one is the opposite face, 
00173             //   one is 'side' face
00174             // - walk on both these faces to opposite edge
00175             // - check if this opposite edge is on faceI
00176 
00177             label f0I, f1I;
00178 
00179             meshTools::getEdgeFaces(mesh, cellI, edgeI, f0I, f1I);
00180 
00181             // Walk to opposite edge on face f0
00182             label edge0I =
00183                 meshTools::walkFace(mesh, f0I, edgeI, e.start(), 2);
00184 
00185             // Check if edge on faceI.
00186 
00187             const edge& e0 = mesh.edges()[edge0I];
00188 
00189             fpA = findIndex(f, e0.start());
00190             fpB = findIndex(f, e0.end());
00191 
00192             if ((fpA != -1) && (fpB != -1))
00193             {
00194                 return lowest(f.size(), fpA, fpB);
00195             }
00196 
00197             // Face0 is doesn't have an edge on faceI (so must be the opposite
00198             // face) so try face1.
00199 
00200             // Walk to opposite edge on face f1
00201             label edge1I =
00202                 meshTools::walkFace(mesh, f1I, edgeI, e.start(), 2);
00203 
00204             // Check if edge on faceI.
00205             const edge& e1 = mesh.edges()[edge1I];
00206 
00207             fpA = findIndex(f, e1.start());
00208             fpB = findIndex(f, e1.end());
00209 
00210             if ((fpA != -1) && (fpB != -1))
00211             {
00212                 return lowest(f.size(), fpA, fpB);
00213             }
00214 
00215             FatalErrorIn("directionInfo::edgeToFaceIndex")
00216                 << "Found connected faces " << mesh.faces()[f0I] << " and "
00217                 << mesh.faces()[f1I] << " sharing edge " << edgeI << endl
00218                 << "But none seems to be connected to face " << faceI
00219                 << " vertices:" << f
00220                 << abort(FatalError);
00221 
00222             return -1;
00223         }
00224     }
00225 }
00226 
00227 
00228 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00229 
00230 // Update this cell with neighbouring face information
00231 bool Foam::directionInfo::updateCell
00232 (
00233     const polyMesh& mesh,
00234     const label thisCellI,
00235     const label neighbourFaceI,
00236     const directionInfo& neighbourInfo,
00237     const scalar        // tol
00238 )
00239 {
00240     if (index_ >= -2)
00241     {
00242         // Already determined.
00243         return false;
00244     }
00245 
00246     if (hexMatcher().isA(mesh, thisCellI))
00247     {
00248         const face& f = mesh.faces()[neighbourFaceI];
00249 
00250         if (neighbourInfo.index() == -2)
00251         {
00252             // Geometric information from neighbour
00253             index_ = -2;
00254         }
00255         else if (neighbourInfo.index() == -1)
00256         {
00257             // Cut tangential to face. Take any edge connected to face
00258             // but not used in face.
00259 
00260             // Get first edge on face.
00261             label edgeI = mesh.faceEdges()[neighbourFaceI][0];
00262 
00263             const edge& e = mesh.edges()[edgeI];
00264 
00265             // Find face connected to face through edgeI and on same cell.
00266             label faceI = 
00267                 meshTools::otherFace
00268                 (
00269                     mesh,
00270                     thisCellI,
00271                     neighbourFaceI,
00272                     edgeI
00273                 );
00274 
00275             // Find edge on faceI which is connected to e.start() but not edgeI.
00276             index_ =
00277                 meshTools::otherEdge
00278                 (
00279                     mesh,
00280                     mesh.faceEdges()[faceI],
00281                     edgeI,
00282                     e.start()
00283                 );
00284         }
00285         else
00286         {
00287             // Index is a vertex on the face. Convert to mesh edge.
00288 
00289             // Get mesh edge between f[index_] and f[index_+1]
00290             label v0 = f[neighbourInfo.index()];
00291             label v1 = f[(neighbourInfo.index() + 1) % f.size()];
00292 
00293             index_ = findEdge(mesh, mesh.faceEdges()[neighbourFaceI], v0, v1);
00294         }
00295     }
00296     else
00297     {
00298         // Not a hex so mark this as geometric.
00299         index_ = -2;
00300     }
00301 
00302 
00303     n_ = neighbourInfo.n();
00304 
00305     return true;
00306 }    
00307 
00308 
00309 // Update this face with neighbouring cell information
00310 bool Foam::directionInfo::updateFace
00311 (
00312     const polyMesh& mesh,
00313     const label thisFaceI,
00314     const label neighbourCellI,
00315     const directionInfo& neighbourInfo,
00316     const scalar    // tol
00317 )
00318 {
00319     // Handle special cases first
00320 
00321     if (index_ >= -2)
00322     {
00323         // Already determined
00324         return false;
00325     }
00326 
00327     // Handle normal cases where topological or geometrical info comes from
00328     // neighbouring cell
00329 
00330     if (neighbourInfo.index() >= 0)
00331     {
00332         // Neighbour has topological direction (and hence is hex). Find cut
00333         // edge on face.
00334         index_ =
00335             edgeToFaceIndex
00336             (
00337                 mesh,
00338                 neighbourCellI,
00339                 thisFaceI,
00340                 neighbourInfo.index()
00341             );
00342     }
00343     else
00344     {
00345         // Neighbour has geometric information. Use.
00346         index_ = -2;
00347     }
00348 
00349 
00350     n_ = neighbourInfo.n();
00351 
00352     return true;
00353 }    
00354 
00355 
00356 // Merge this with information on same face
00357 bool Foam::directionInfo::updateFace
00358 (
00359     const polyMesh& mesh,
00360     const label,    // thisFaceI
00361     const directionInfo& neighbourInfo,
00362     const scalar    // tol
00363 )
00364 {
00365     if (index_ >= -2)
00366     {
00367         // Already visited.
00368         return false;
00369     }
00370     else
00371     {
00372         index_ = neighbourInfo.index();
00373 
00374         n_ = neighbourInfo.n();
00375 
00376         return true;
00377     }
00378 }    
00379 
00380 
00381 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
00382 
00383 Foam::Ostream& Foam::operator<<
00384 (
00385     Foam::Ostream& os,
00386     const Foam::directionInfo& wDist
00387 )
00388 {
00389     return os << wDist.index_ << wDist.n_;
00390 }
00391 
00392 
00393 Foam::Istream& Foam::operator>>(Foam::Istream& is, Foam::directionInfo& wDist)
00394 {
00395     return is >> wDist.index_ >> wDist.n_;
00396 }
00397 
00398 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines