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