FreeFOAM The Cross-Platform CFD Toolkit
directionInfo.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 
25 \*---------------------------------------------------------------------------*/
26 
27 #include "directionInfo.H"
28 #include <OpenFOAM/hexMatcher.H>
29 #include <meshTools/meshTools.H>
30 #include <OpenFOAM/polyMesh.H>
31 
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 bool Foam::directionInfo::equal(const edge& e, const label v0, const label v1)
36 {
37  return
38  (e.start() == v0 && e.end() == v1)
39  || (e.start() == v1 && e.end() == v0);
40 }
41 
42 
43 Foam::point Foam::directionInfo::eMid
44 (
45  const primitiveMesh& mesh,
46  const label edgeI
47 )
48 {
49  const edge& e = mesh.edges()[edgeI];
50 
51  return
52  0.5
53  * (mesh.points()[e.start()] + mesh.points()[e.end()]);
54 }
55 
56 
57 // Find edge among edgeLabels that uses v0 and v1
59 (
60  const primitiveMesh& mesh,
61  const labelList& edgeLabels,
62  const label v1,
63  const label v0
64 )
65 {
66  forAll(edgeLabels, edgeLabelI)
67  {
68  label edgeI = edgeLabels[edgeLabelI];
69 
70  const edge& e = mesh.edges()[edgeI];
71 
72  if (equal(e, v0, v1))
73  {
74  return edgeI;
75  }
76  }
77 
78  FatalErrorIn("directionInfo::findEdge")
79  << "Cannot find an edge among " << edgeLabels << endl
80  << "that uses vertices " << v0
81  << " and " << v1
82  << abort(FatalError);
83 
84  return -1;
85 }
86 
87 
88 Foam::label Foam::directionInfo::lowest
89 (
90  const label size,
91  const label a,
92  const label b
93 )
94 {
95  // Get next point
96  label a1 = (a + 1) % size;
97 
98  if (a1 == b)
99  {
100  return a;
101  }
102  else
103  {
104  label b1 = (b + 1) % size;
105 
106  if (b1 != a)
107  {
108  FatalErrorIn("directionInfo::lowest")
109  << "Problem : a:" << a << " b:" << b << " size:" << size
110  << abort(FatalError);
111  }
112 
113  return b;
114  }
115 }
116 
117 
118 // Have edge on hex cell. Find corresponding edge on face. Return -1 if none
119 // found.
121 (
122  const primitiveMesh& mesh,
123  const label cellI,
124  const label faceI,
125  const label edgeI
126 )
127 {
128  if ((edgeI < 0) || (edgeI >= mesh.nEdges()))
129  {
130  FatalErrorIn("directionInfo::edgeToFaceIndex")
131  << "Illegal edge label:" << edgeI
132  << " when projecting cut edge from cell " << cellI
133  << " to face " << faceI
134  << abort(FatalError);
135  }
136 
137  const edge& e = mesh.edges()[edgeI];
138 
139  const face& f = mesh.faces()[faceI];
140 
141  // edgeI is either
142  // - in faceI. Convert into index in face.
143  // - connected (but not in) to face. Return -1.
144  // - in face opposite faceI. Convert into index in face.
145 
146  label fpA = findIndex(f, e.start());
147  label fpB = findIndex(f, e.end());
148 
149  if (fpA != -1)
150  {
151  if (fpB != -1)
152  {
153  return lowest(f.size(), fpA, fpB);
154  }
155  else
156  {
157  // e.start() in face, e.end() not
158  return -1;
159  }
160  }
161  else
162  {
163  if (fpB != -1)
164  {
165  // e.end() in face, e.start() not
166  return -1;
167  }
168  else
169  {
170  // Both not in face.
171  // e is on opposite face. Determine corresponding edge on this face:
172  // - determine two faces using edge (one is the opposite face,
173  // one is 'side' face
174  // - walk on both these faces to opposite edge
175  // - check if this opposite edge is on faceI
176 
177  label f0I, f1I;
178 
179  meshTools::getEdgeFaces(mesh, cellI, edgeI, f0I, f1I);
180 
181  // Walk to opposite edge on face f0
182  label edge0I =
183  meshTools::walkFace(mesh, f0I, edgeI, e.start(), 2);
184 
185  // Check if edge on faceI.
186 
187  const edge& e0 = mesh.edges()[edge0I];
188 
189  fpA = findIndex(f, e0.start());
190  fpB = findIndex(f, e0.end());
191 
192  if ((fpA != -1) && (fpB != -1))
193  {
194  return lowest(f.size(), fpA, fpB);
195  }
196 
197  // Face0 is doesn't have an edge on faceI (so must be the opposite
198  // face) so try face1.
199 
200  // Walk to opposite edge on face f1
201  label edge1I =
202  meshTools::walkFace(mesh, f1I, edgeI, e.start(), 2);
203 
204  // Check if edge on faceI.
205  const edge& e1 = mesh.edges()[edge1I];
206 
207  fpA = findIndex(f, e1.start());
208  fpB = findIndex(f, e1.end());
209 
210  if ((fpA != -1) && (fpB != -1))
211  {
212  return lowest(f.size(), fpA, fpB);
213  }
214 
215  FatalErrorIn("directionInfo::edgeToFaceIndex")
216  << "Found connected faces " << mesh.faces()[f0I] << " and "
217  << mesh.faces()[f1I] << " sharing edge " << edgeI << endl
218  << "But none seems to be connected to face " << faceI
219  << " vertices:" << f
220  << abort(FatalError);
221 
222  return -1;
223  }
224  }
225 }
226 
227 
228 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
229 
230 // Update this cell with neighbouring face information
232 (
233  const polyMesh& mesh,
234  const label thisCellI,
235  const label neighbourFaceI,
236  const directionInfo& neighbourInfo,
237  const scalar // tol
238 )
239 {
240  if (index_ >= -2)
241  {
242  // Already determined.
243  return false;
244  }
245 
246  if (hexMatcher().isA(mesh, thisCellI))
247  {
248  const face& f = mesh.faces()[neighbourFaceI];
249 
250  if (neighbourInfo.index() == -2)
251  {
252  // Geometric information from neighbour
253  index_ = -2;
254  }
255  else if (neighbourInfo.index() == -1)
256  {
257  // Cut tangential to face. Take any edge connected to face
258  // but not used in face.
259 
260  // Get first edge on face.
261  label edgeI = mesh.faceEdges()[neighbourFaceI][0];
262 
263  const edge& e = mesh.edges()[edgeI];
264 
265  // Find face connected to face through edgeI and on same cell.
266  label faceI =
268  (
269  mesh,
270  thisCellI,
271  neighbourFaceI,
272  edgeI
273  );
274 
275  // Find edge on faceI which is connected to e.start() but not edgeI.
276  index_ =
278  (
279  mesh,
280  mesh.faceEdges()[faceI],
281  edgeI,
282  e.start()
283  );
284  }
285  else
286  {
287  // Index is a vertex on the face. Convert to mesh edge.
288 
289  // Get mesh edge between f[index_] and f[index_+1]
290  label v0 = f[neighbourInfo.index()];
291  label v1 = f[(neighbourInfo.index() + 1) % f.size()];
292 
293  index_ = findEdge(mesh, mesh.faceEdges()[neighbourFaceI], v0, v1);
294  }
295  }
296  else
297  {
298  // Not a hex so mark this as geometric.
299  index_ = -2;
300  }
301 
302 
303  n_ = neighbourInfo.n();
304 
305  return true;
306 }
307 
308 
309 // Update this face with neighbouring cell information
311 (
312  const polyMesh& mesh,
313  const label thisFaceI,
314  const label neighbourCellI,
315  const directionInfo& neighbourInfo,
316  const scalar // tol
317 )
318 {
319  // Handle special cases first
320 
321  if (index_ >= -2)
322  {
323  // Already determined
324  return false;
325  }
326 
327  // Handle normal cases where topological or geometrical info comes from
328  // neighbouring cell
329 
330  if (neighbourInfo.index() >= 0)
331  {
332  // Neighbour has topological direction (and hence is hex). Find cut
333  // edge on face.
334  index_ =
335  edgeToFaceIndex
336  (
337  mesh,
338  neighbourCellI,
339  thisFaceI,
340  neighbourInfo.index()
341  );
342  }
343  else
344  {
345  // Neighbour has geometric information. Use.
346  index_ = -2;
347  }
348 
349 
350  n_ = neighbourInfo.n();
351 
352  return true;
353 }
354 
355 
356 // Merge this with information on same face
358 (
359  const polyMesh& mesh,
360  const label, // thisFaceI
361  const directionInfo& neighbourInfo,
362  const scalar // tol
363 )
364 {
365  if (index_ >= -2)
366  {
367  // Already visited.
368  return false;
369  }
370  else
371  {
372  index_ = neighbourInfo.index();
373 
374  n_ = neighbourInfo.n();
375 
376  return true;
377  }
378 }
379 
380 
381 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
382 
383 Foam::Ostream& Foam::operator<<
384 (
385  Foam::Ostream& os,
386  const Foam::directionInfo& wDist
387 )
388 {
389  return os << wDist.index_ << wDist.n_;
390 }
391 
392 
394 {
395  return is >> wDist.index_ >> wDist.n_;
396 }
397 
398 // ************************ vim: set sw=4 sts=4 et: ************************ //