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