FreeFOAM The Cross-Platform CFD Toolkit
triSurfaceAddressing.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 Description
25  Contains fix for PrimitivePatch addressing (which doesn't work if surface
26  is non-manifold). Should be moved into PrimitivePatch.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include "triSurface.H"
31 #include <OpenFOAM/HashTable.H>
32 #include <OpenFOAM/SortableList.H>
33 #include <OpenFOAM/transform.H>
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void triSurface::calcSortedEdgeFaces() const
43 {
44  if (sortedEdgeFacesPtr_)
45  {
46  FatalErrorIn("triSurface::calcSortedEdgeFaces()")
47  << "sortedEdgeFacesPtr_ already set"
48  << abort(FatalError);
49  }
50 
51  const labelListList& eFaces = edgeFaces();
52 
53  // create the lists for the various results. (resized on completion)
54  sortedEdgeFacesPtr_ = new labelListList(eFaces.size());
55  labelListList& sortedEdgeFaces = *sortedEdgeFacesPtr_;
56 
57  forAll(eFaces, edgeI)
58  {
59  const labelList& myFaceNbs = eFaces[edgeI];
60 
61  if (myFaceNbs.size() > 2)
62  {
63  // Get point on edge and normalized direction of edge (= e2 base
64  // of our coordinate system)
65  const edge& e = edges()[edgeI];
66 
67  const point& edgePt = localPoints()[e.start()];
68 
69  vector e2 = e.vec(localPoints());
70  e2 /= mag(e2) + VSMALL;
71 
72 
73  // Get opposite vertex for 0th face
74  const labelledTri& f = localFaces()[myFaceNbs[0]];
75  label fp0 = findIndex(f, e[0]);
76  label fp1 = f.fcIndex(fp0);
77  label vertI = (f[fp1] != e[1] ? f[fp1] : f.fcIndex(fp1));
78 
79  // Get vector normal both to e2 and to edge from opposite vertex
80  // to edge (will be x-axis of our coordinate system)
81  vector e0 = e2 ^ (localPoints()[vertI] - edgePt);
82  e0 /= mag(e0) + VSMALL;
83 
84  // Get y-axis of coordinate system
85  vector e1 = e2 ^ e0;
86 
87 
88  SortableList<scalar> faceAngles(myFaceNbs.size());
89 
90  // e0 is reference so angle is 0
91  faceAngles[0] = 0;
92 
93  for(label nbI = 1; nbI < myFaceNbs.size(); nbI++)
94  {
95  // Get opposite vertex
96  const labelledTri& f = localFaces()[myFaceNbs[nbI]];
97  label fp0 = findIndex(f, e[0]);
98  label fp1 = f.fcIndex(fp0);
99  label vertI = (f[fp1] != e[1] ? f[fp1] : f.fcIndex(fp1));
100 
101  vector vec = e2 ^ (localPoints()[vertI] - edgePt);
102  vec /= mag(vec) + VSMALL;
103 
104  faceAngles[nbI] = pseudoAngle
105  (
106  e0,
107  e1,
108  vec
109  );
110  }
111 
112  faceAngles.sort();
113 
114  sortedEdgeFaces[edgeI] = UIndirectList<label>
115  (
116  myFaceNbs,
117  faceAngles.indices()
118  );
119  }
120  else
121  {
122  // No need to sort. Just copy.
123  sortedEdgeFaces[edgeI] = myFaceNbs;
124  }
125  }
126 }
127 
128 
129 void triSurface::calcEdgeOwner() const
130 {
131  if (edgeOwnerPtr_)
132  {
133  FatalErrorIn("triSurface::calcEdgeOwner()")
134  << "edgeOwnerPtr_ already set"
135  << abort(FatalError);
136  }
137 
138  // create the owner list
139  edgeOwnerPtr_ = new labelList(nEdges());
140  labelList& edgeOwner = *edgeOwnerPtr_;
141 
142  forAll(edges(), edgeI)
143  {
144  const edge& e = edges()[edgeI];
145 
146  const labelList& myFaces = edgeFaces()[edgeI];
147 
148  if (myFaces.size() == 1)
149  {
150  edgeOwner[edgeI] = myFaces[0];
151  }
152  else
153  {
154  // Find the first face whose vertices are aligned with the edge.
155  // (in case of multiply connected edge the best we can do)
156  edgeOwner[edgeI] = -1;
157 
158  forAll(myFaces, i)
159  {
160  const labelledTri& f = localFaces()[myFaces[i]];
161 
162  if
163  (
164  ((f[0] == e.start()) && (f[1] == e.end()))
165  || ((f[1] == e.start()) && (f[2] == e.end()))
166  || ((f[2] == e.start()) && (f[0] == e.end()))
167  )
168  {
169  edgeOwner[edgeI] = myFaces[i];
170 
171  break;
172  }
173  }
174 
175  if (edgeOwner[edgeI] == -1)
176  {
177  FatalErrorIn("triSurface::calcEdgeOwner()")
178  << "Edge " << edgeI << " vertices:" << e
179  << " is used by faces " << myFaces
180  << " vertices:"
181  << UIndirectList<labelledTri>(localFaces(), myFaces)()
182  << " none of which use the edge vertices in the same order"
183  << nl << "I give up" << abort(FatalError);
184  }
185  }
186  }
187 }
188 
189 
190 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
191 
192 } // End namespace Foam
193 
194 // ************************ vim: set sw=4 sts=4 et: ************************ //