src/triSurface/triSurface/triSurfaceAddressing.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Contains fix for PrimitivePatch addressing (which doesn't work if surface
  27     is non-manifold). Should be moved into PrimitivePatch.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "triSurface.H"
  32 #include "HashTable.H"
  33 #include "SortableList.H"
  34 #include "transform.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 void triSurface::calcSortedEdgeFaces() const
  44 {
  45     if (sortedEdgeFacesPtr_)
  46     {
  47         FatalErrorIn("triSurface::calcSortedEdgeFaces()")
  48             << "sortedEdgeFacesPtr_ already set"
  49             << abort(FatalError);
  50     }
  51
  52     const labelListList& eFaces = edgeFaces();
  53
  54     // create the lists for the various results. (resized on completion)
  55     sortedEdgeFacesPtr_ = new labelListList(eFaces.size());
  56     labelListList& sortedEdgeFaces = *sortedEdgeFacesPtr_;
  57
  58     forAll(eFaces, edgeI)
  59     {
  60         const labelList& myFaceNbs = eFaces[edgeI];
  61
  62         if (myFaceNbs.size() > 2)
  63         {
  64             // Get point on edge and normalized direction of edge (= e2 base
  65             // of our coordinate system)
  66             const edge& e = edges()[edgeI];
  67
  68             const point& edgePt = localPoints()[e.start()];
  69
  70             vector e2 = e.vec(localPoints());
  71             e2 /= mag(e2) + VSMALL;
  72
  73
  74             // Get opposite vertex for 0th face
  75             const labelledTri& f = localFaces()[myFaceNbs[0]];
  76             label fp0 = findIndex(f, e[0]);
  77             label fp1 = f.fcIndex(fp0);
  78             label vertI = (f[fp1] != e[1] ? f[fp1] : f.fcIndex(fp1));
  79
  80             // Get vector normal both to e2 and to edge from opposite vertex
  81             // to edge (will be x-axis of our coordinate system)
  82             vector e0 = e2 ^ (localPoints()[vertI] - edgePt);
  83             e0 /= mag(e0) + VSMALL;
  84
  85             // Get y-axis of coordinate system
  86             vector e1 = e2 ^ e0;
  87
  88
  89             SortableList<scalar> faceAngles(myFaceNbs.size());
  90
  91             // e0 is reference so angle is 0
  92             faceAngles[0] = 0;
  93
  94             for(label nbI = 1; nbI < myFaceNbs.size(); nbI++)
  95             {
  96                 // Get opposite vertex
  97                 const labelledTri& f = localFaces()[myFaceNbs[nbI]];
  98                 label fp0 = findIndex(f, e[0]);
  99                 label fp1 = f.fcIndex(fp0);
 100                 label vertI = (f[fp1] != e[1] ? f[fp1] : f.fcIndex(fp1));
 101
 102                 vector vec = e2 ^ (localPoints()[vertI] - edgePt);
 103                 vec /= mag(vec) + VSMALL;
 104
 105                 faceAngles[nbI] = pseudoAngle
 106                 (
 107                     e0,
 108                     e1,
 109                     vec
 110                 );
 111             }
 112
 113             faceAngles.sort();
 114
 115             sortedEdgeFaces[edgeI] = UIndirectList<label>
 116             (
 117                 myFaceNbs,
 118                 faceAngles.indices()
 119             );
 120         }
 121         else
 122         {
 123             // No need to sort. Just copy.
 124             sortedEdgeFaces[edgeI] = myFaceNbs;
 125         }
 126     }
 127 }
 128
 129
 130 void triSurface::calcEdgeOwner() const
 131 {
 132     if (edgeOwnerPtr_)
 133     {
 134         FatalErrorIn("triSurface::calcEdgeOwner()")
 135             << "edgeOwnerPtr_ already set"
 136             << abort(FatalError);
 137     }
 138
 139     // create the owner list
 140     edgeOwnerPtr_ = new labelList(nEdges());
 141     labelList& edgeOwner = *edgeOwnerPtr_;
 142
 143     forAll(edges(), edgeI)
 144     {
 145         const edge& e = edges()[edgeI];
 146
 147         const labelList& myFaces = edgeFaces()[edgeI];
 148
 149         if (myFaces.size() == 1)
 150         {
 151             edgeOwner[edgeI] = myFaces[0];
 152         }
 153         else
 154         {
 155             // Find the first face whose vertices are aligned with the edge.
 156             // (in case of multiply connected edge the best we can do)
 157             edgeOwner[edgeI] = -1;
 158
 159             forAll(myFaces, i)
 160             {
 161                 const labelledTri& f = localFaces()[myFaces[i]];
 162
 163                 if
 164                 (
 165                     ((f[0] == e.start()) && (f[1] == e.end()))
 166                  || ((f[1] == e.start()) && (f[2] == e.end()))
 167                  || ((f[2] == e.start()) && (f[0] == e.end()))
 168                 )
 169                 {
 170                     edgeOwner[edgeI] = myFaces[i];
 171
 172                     break;
 173                 }
 174             }
 175
 176             if (edgeOwner[edgeI] == -1)
 177             {
 178                 FatalErrorIn("triSurface::calcEdgeOwner()")
 179                     << "Edge " << edgeI << " vertices:" << e
 180                     << " is used by faces " << myFaces
 181                     << " vertices:"
 182                     << UIndirectList<labelledTri>(localFaces(), myFaces)()
 183                     << " none of which use the edge vertices in the same order"
 184                     << nl << "I give up" << abort(FatalError);
 185             }
 186         }
 187     }
 188 }
 189
 190
 191 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 192
 193 } // End namespace Foam
 194
 195 // ************************************************************************* //