FreeFOAM The Cross-Platform CFD Toolkit
cellFeatures.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 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "cellFeatures.H"
29 #include <OpenFOAM/primitiveMesh.H>
30 #include <OpenFOAM/HashSet.H>
31 #include <OpenFOAM/Map.H>
33 #include <OpenFOAM/ListOps.H>
34 #include <meshTools/meshTools.H>
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 
39 // Return true if edge start and end are on increasing face vertices. (edge is
40 // guaranteed to be on face)
41 bool Foam::cellFeatures::faceAlignedEdge(const label faceI, const label edgeI)
42  const
43 {
44  const edge& e = mesh_.edges()[edgeI];
45 
46  const face& f = mesh_.faces()[faceI];
47 
48  forAll(f, fp)
49  {
50  if (f[fp] == e.start())
51  {
52  label fp1 = f.fcIndex(fp);
53 
54  return f[fp1] == e.end();
55  }
56  }
57 
59  (
60  "cellFeatures::faceAlignedEdge(const label, const label)"
61  ) << "Can not find edge " << mesh_.edges()[edgeI]
62  << " on face " << faceI << abort(FatalError);
63 
64  return false;
65 }
66 
67 
68 // Return edge in featureEdge that uses vertI and is on same superface
69 // but is not edgeI
70 Foam::label Foam::cellFeatures::nextEdge
71 (
72  const Map<label>& toSuperFace,
73  const label superFaceI,
74  const label thisEdgeI,
75  const label thisVertI
76 ) const
77 {
78  const labelList& pEdges = mesh_.pointEdges()[thisVertI];
79 
80  forAll(pEdges, pEdgeI)
81  {
82  label edgeI = pEdges[pEdgeI];
83 
84  if ((edgeI != thisEdgeI) && featureEdge_.found(edgeI))
85  {
86  // Check that edge is used by a face on same superFace
87 
88  const labelList& eFaces = mesh_.edgeFaces()[edgeI];
89 
90  forAll(eFaces, eFaceI)
91  {
92  label faceI = eFaces[eFaceI];
93 
94  if
95  (
96  meshTools::faceOnCell(mesh_, cellI_, faceI)
97  && (toSuperFace[faceI] == superFaceI)
98  )
99  {
100  return edgeI;
101  }
102  }
103  }
104  }
105 
107  (
108  "cellFeatures::nextEdge(const label, const Map<label>"
109  ", const labelHashSet&, const label, const label, const label)"
110  ) << "Can not find edge in " << featureEdge_ << " connected to edge "
111  << thisEdgeI << " at vertex " << thisVertI << endl
112  << "This might mean that the externalEdges do not form a closed loop"
113  << abort(FatalError);
114 
115  return -1;
116 }
117 
118 
119 // Return true if angle between faces using it is larger than certain value.
120 bool Foam::cellFeatures::isCellFeatureEdge
121 (
122  const scalar minCos,
123  const label edgeI
124 ) const
125 {
126  // Get the two faces using this edge
127 
128  label face0;
129  label face1;
130  meshTools::getEdgeFaces(mesh_, cellI_, edgeI, face0, face1);
131 
132  // Check the angle between them by comparing the face normals.
133 
134  vector n0 = mesh_.faceAreas()[face0];
135  n0 /= mag(n0);
136 
137  vector n1 = mesh_.faceAreas()[face1];
138  n1 /= mag(n1);
139 
140  scalar cosAngle = n0 & n1;
141 
142 
143  const edge& e = mesh_.edges()[edgeI];
144 
145  const face& f0 = mesh_.faces()[face0];
146 
147  label face0Start = findIndex(f0, e.start());
148  label face0End = f0.fcIndex(face0Start);
149 
150  const face& f1 = mesh_.faces()[face1];
151 
152  label face1Start = findIndex(f1, e.start());
153  label face1End = f1.fcIndex(face1Start);
154 
155  if
156  (
157  (
158  (f0[face0End] == e.end())
159  && (f1[face1End] != e.end())
160  )
161  || (
162  (f0[face0End] != e.end())
163  && (f1[face1End] == e.end())
164  )
165  )
166  {
167  }
168  else
169  {
170  cosAngle = -cosAngle;
171  }
172 
173  if (cosAngle < minCos)
174  {
175  return true;
176  }
177  else
178  {
179  return false;
180  }
181 }
182 
183 
184 // Recursively mark (on toSuperFace) all face reachable across non-feature
185 // edges.
186 void Foam::cellFeatures::walkSuperFace
187 (
188  const label faceI,
189  const label superFaceI,
190  Map<label>& toSuperFace
191 ) const
192 {
193  if (!toSuperFace.found(faceI))
194  {
195  toSuperFace.insert(faceI, superFaceI);
196 
197  const labelList& fEdges = mesh_.faceEdges()[faceI];
198 
199  forAll(fEdges, fEdgeI)
200  {
201  label edgeI = fEdges[fEdgeI];
202 
203  if (!featureEdge_.found(edgeI))
204  {
205  label face0;
206  label face1;
207  meshTools::getEdgeFaces(mesh_, cellI_, edgeI, face0, face1);
208 
209  if (face0 == faceI)
210  {
211  face0 = face1;
212  }
213 
214  walkSuperFace
215  (
216  face0,
217  superFaceI,
218  toSuperFace
219  );
220  }
221  }
222  }
223 }
224 
225 
226 void Foam::cellFeatures::calcSuperFaces() const
227 {
228  // Determine superfaces by edge walking across non-feature edges
229 
230  const labelList& cFaces = mesh_.cells()[cellI_];
231 
232  // Mapping from old to super face:
233  // <not found> : not visited
234  // >=0 : superFace
235  Map<label> toSuperFace(10*cFaces.size());
236 
237  label superFaceI = 0;
238 
239  forAll(cFaces, cFaceI)
240  {
241  label faceI = cFaces[cFaceI];
242 
243  if (!toSuperFace.found(faceI))
244  {
245  walkSuperFace
246  (
247  faceI,
248  superFaceI,
249  toSuperFace
250  );
251  superFaceI++;
252  }
253  }
254 
255  // Construct superFace-to-oldface mapping.
256 
257  faceMap_.setSize(superFaceI);
258 
259  forAll(cFaces, cFaceI)
260  {
261  label faceI = cFaces[cFaceI];
262 
263  faceMap_[toSuperFace[faceI]].append(faceI);
264  }
265 
266  forAll(faceMap_, superI)
267  {
268  faceMap_[superI].shrink();
269  }
270 
271 
272  // Construct superFaces
273 
274  facesPtr_ = new faceList(superFaceI);
275 
276  faceList& faces = *facesPtr_;
277 
278  forAll(cFaces, cFaceI)
279  {
280  label faceI = cFaces[cFaceI];
281 
282  label superFaceI = toSuperFace[faceI];
283 
284  if (faces[superFaceI].empty())
285  {
286  // Superface not yet constructed.
287 
288  // Find starting feature edge on face.
289  label startEdgeI = -1;
290 
291  const labelList& fEdges = mesh_.faceEdges()[faceI];
292 
293  forAll(fEdges, fEdgeI)
294  {
295  label edgeI = fEdges[fEdgeI];
296 
297  if (featureEdge_.found(edgeI))
298  {
299  startEdgeI = edgeI;
300 
301  break;
302  }
303  }
304 
305 
306  if (startEdgeI != -1)
307  {
308  // Walk point-edge-point along feature edges
309 
310  DynamicList<label> superFace(10*mesh_.faces()[faceI].size());
311 
312  const edge& e = mesh_.edges()[startEdgeI];
313 
314  // Walk either start-end or end-start depending on orientation
315  // of face. SuperFace will have cellI as owner.
316  bool flipOrientation =
317  (mesh_.faceOwner()[faceI] == cellI_)
318  ^ (faceAlignedEdge(faceI, startEdgeI));
319 
320  label startVertI = -1;
321 
322  if (flipOrientation)
323  {
324  startVertI = e.end();
325  }
326  else
327  {
328  startVertI = e.start();
329  }
330 
331  label edgeI = startEdgeI;
332 
333  label vertI = e.otherVertex(startVertI);
334 
335  do
336  {
337  label newEdgeI = nextEdge
338  (
339  toSuperFace,
340  superFaceI,
341  edgeI,
342  vertI
343  );
344 
345  // Determine angle between edges.
346  if (isFeaturePoint(edgeI, newEdgeI))
347  {
348  superFace.append(vertI);
349  }
350 
351  edgeI = newEdgeI;
352 
353  if (vertI == startVertI)
354  {
355  break;
356  }
357 
358  vertI = mesh_.edges()[edgeI].otherVertex(vertI);
359  }
360  while (true);
361 
362  if (superFace.size() <= 2)
363  {
364  WarningIn("cellFeatures::calcSuperFaces")
365  << " Can not collapse faces " << faceMap_[superFaceI]
366  << " into one big face on cell " << cellI_ << endl
367  << "Try decreasing minCos:" << minCos_ << endl;
368  }
369  else
370  {
371  faces[superFaceI].transfer(superFace);
372  }
373  }
374  }
375  }
376 }
377 
378 
379 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
380 
381 // Construct from components
382 Foam::cellFeatures::cellFeatures
383 (
384  const primitiveMesh& mesh,
385  const scalar minCos,
386  const label cellI
387 )
388 :
389  mesh_(mesh),
390  minCos_(minCos),
391  cellI_(cellI),
392  featureEdge_(10*mesh.cellEdges()[cellI].size()),
393  facesPtr_(NULL),
394  faceMap_(0)
395 {
396  const labelList& cEdges = mesh_.cellEdges()[cellI_];
397 
398  forAll(cEdges, cEdgeI)
399  {
400  label edgeI = cEdges[cEdgeI];
401 
402  if (isCellFeatureEdge(minCos_, edgeI))
403  {
404  featureEdge_.insert(edgeI);
405  }
406  }
407 
408 }
409 
410 
411 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
412 
414 {
415  deleteDemandDrivenData(facesPtr_);
416 }
417 
418 
419 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
420 
421 bool Foam::cellFeatures::isFeaturePoint(const label edge0, const label edge1)
422  const
423 {
424  if
425  (
426  (edge0 < 0)
427  || (edge0 >= mesh_.nEdges())
428  || (edge1 < 0)
429  || (edge1 >= mesh_.nEdges())
430  )
431  {
433  (
434  "cellFeatures::isFeatureVertex(const label, const label)"
435  ) << "Illegal edge labels : edge0:" << edge0 << " edge1:" << edge1
436  << abort(FatalError);
437  }
438 
439  const edge& e0 = mesh_.edges()[edge0];
440 
441  vector e0Vec = e0.vec(mesh_.points());
442  e0Vec /= mag(e0Vec);
443 
444  const edge& e1 = mesh_.edges()[edge1];
445 
446  vector e1Vec = e1.vec(mesh_.points());
447  e1Vec /= mag(e1Vec);
448 
449  scalar cosAngle;
450 
451  if
452  (
453  (e0.start() == e1.end())
454  || (e0.end() == e1.start())
455  )
456  {
457  // Same direction
458  cosAngle = e0Vec & e1Vec;
459  }
460  else if
461  (
462  (e0.start() == e1.start())
463  || (e0.end() == e1.end())
464  )
465  {
466  // back on back
467  cosAngle = - e0Vec & e1Vec;
468  }
469  else
470  {
471  cosAngle = GREAT; // satisfy compiler
472 
474  (
475  "cellFeatures::isFeaturePoint(const label, const label"
476  ", const label)"
477  ) << "Edges do not share common vertex. e0:" << e0
478  << " e1:" << e1 << abort(FatalError);
479  }
480 
481  if (cosAngle < minCos_)
482  {
483  // Angle larger than criterium
484  return true;
485  }
486  else
487  {
488  return false;
489  }
490 }
491 
492 
493 bool Foam::cellFeatures::isFeatureVertex(const label faceI, const label vertI)
494  const
495 {
496  if
497  (
498  (faceI < 0)
499  || (faceI >= mesh_.nFaces())
500  || (vertI < 0)
501  || (vertI >= mesh_.nPoints())
502  )
503  {
505  (
506  "cellFeatures::isFeatureVertex(const label, const label)"
507  ) << "Illegal face " << faceI << " or vertex " << vertI
508  << abort(FatalError);
509  }
510 
511  const labelList& pEdges = mesh_.pointEdges()[vertI];
512 
513  label edge0 = -1;
514  label edge1 = -1;
515 
516  forAll(pEdges, pEdgeI)
517  {
518  label edgeI = pEdges[pEdgeI];
519 
520  if (meshTools::edgeOnFace(mesh_, faceI, edgeI))
521  {
522  if (edge0 == -1)
523  {
524  edge0 = edgeI;
525  }
526  else
527  {
528  edge1 = edgeI;
529 
530  // Found the two edges.
531  break;
532  }
533  }
534  }
535 
536  if (edge1 == -1)
537  {
539  (
540  "cellFeatures::isFeatureVertex(const label, const label)"
541  ) << "Did not find two edges sharing vertex " << vertI
542  << " on face " << faceI << " vertices:" << mesh_.faces()[faceI]
543  << abort(FatalError);
544  }
545 
546  return isFeaturePoint(edge0, edge1);
547 }
548 
549 
550 // ************************ vim: set sw=4 sts=4 et: ************************ //