meshTriangulation.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-2016 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 "meshTriangulation.H"
27 #include "polyMesh.H"
28 #include "faceTriangulation.H"
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 bool Foam::meshTriangulation::isInternalFace
33 (
34  const primitiveMesh& mesh,
35  const boolList& includedCell,
36  const label facei
37 )
38 {
39  if (mesh.isInternalFace(facei))
40  {
41  label own = mesh.faceOwner()[facei];
42  label nei = mesh.faceNeighbour()[facei];
43 
44  if (includedCell[own] && includedCell[nei])
45  {
46  // Neighbouring cell will get included in subset
47  // as well so face is internal.
48  return true;
49  }
50  else
51  {
52  return false;
53  }
54  }
55  else
56  {
57  return false;
58  }
59 }
60 
61 
62 void Foam::meshTriangulation::getFaces
63 (
64  const primitiveMesh& mesh,
65  const boolList& includedCell,
66  boolList& faceIsCut,
67  label& nFaces,
68  label& nInternalFaces
69 )
70 {
71  // All faces to be triangulated.
72  faceIsCut.setSize(mesh.nFaces());
73  faceIsCut = false;
74 
75  nFaces = 0;
76  nInternalFaces = 0;
77 
78  forAll(includedCell, celli)
79  {
80  // Include faces of cut cells only.
81  if (includedCell[celli])
82  {
83  const labelList& cFaces = mesh.cells()[celli];
84 
85  forAll(cFaces, i)
86  {
87  label facei = cFaces[i];
88 
89  if (!faceIsCut[facei])
90  {
91  // First visit of face.
92  nFaces++;
93  faceIsCut[facei] = true;
94 
95  // See if would become internal or external face
96  if (isInternalFace(mesh, includedCell, facei))
97  {
98  nInternalFaces++;
99  }
100  }
101  }
102  }
103  }
104 
105  Pout<< "Subset consists of " << nFaces << " faces out of " << mesh.nFaces()
106  << " of which " << nInternalFaces << " are internal" << endl;
107 }
108 
109 
110 void Foam::meshTriangulation::insertTriangles
111 (
112  const triFaceList& faceTris,
113  const label facei,
114  const label regionI,
115  const bool reverse,
116 
117  List<labelledTri>& triangles,
118  label& triI
119 )
120 {
121  // Copy triangles. Optionally reverse them
122  forAll(faceTris, i)
123  {
124  const triFace& f = faceTris[i];
125 
126  labelledTri& tri = triangles[triI];
127 
128  if (reverse)
129  {
130  tri[0] = f[0];
131  tri[2] = f[1];
132  tri[1] = f[2];
133  }
134  else
135  {
136  tri[0] = f[0];
137  tri[1] = f[1];
138  tri[2] = f[2];
139  }
140 
141  tri.region() = regionI;
142 
143  faceMap_[triI] = facei;
144 
145  triI++;
146  }
147 }
148 
149 
150 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
151 
152 // Null constructor
154 :
155  triSurface(),
156  nInternalFaces_(0),
157  faceMap_()
158 {}
159 
160 
161 // Construct from faces of cells
163 (
164  const polyMesh& mesh,
165  const label internalFacesPatch,
166  const boolList& includedCell,
167  const bool faceCentreDecomposition
168 )
169 :
170  triSurface(),
171  nInternalFaces_(0),
172  faceMap_()
173 {
174  const faceList& faces = mesh.faces();
175  const pointField& points = mesh.points();
176  const polyBoundaryMesh& patches = mesh.boundaryMesh();
177 
178  // All faces to be triangulated.
179  boolList faceIsCut;
180  label nFaces, nInternalFaces;
181 
182  getFaces
183  (
184  mesh,
185  includedCell,
186  faceIsCut,
187  nFaces,
188  nInternalFaces
189  );
190 
191 
192  // Find upper limit for number of triangles
193  // (can be less if triangulation fails)
194  label nTotTri = 0;
195 
196  if (faceCentreDecomposition)
197  {
198  forAll(faceIsCut, facei)
199  {
200  if (faceIsCut[facei])
201  {
202  nTotTri += faces[facei].size();
203  }
204  }
205  }
206  else
207  {
208  forAll(faceIsCut, facei)
209  {
210  if (faceIsCut[facei])
211  {
212  nTotTri += faces[facei].nTriangles(points);
213  }
214  }
215  }
216  Pout<< "nTotTri : " << nTotTri << endl;
217 
218 
219  // Storage for new and old points (only for faceCentre decomposition;
220  // for triangulation uses only existing points)
221  pointField newPoints;
222 
223  if (faceCentreDecomposition)
224  {
225  newPoints.setSize(mesh.nPoints() + faces.size());
226  forAll(mesh.points(), pointi)
227  {
228  newPoints[pointi] = mesh.points()[pointi];
229  }
230  // Face centres
231  forAll(faces, facei)
232  {
233  newPoints[mesh.nPoints() + facei] = mesh.faceCentres()[facei];
234  }
235  }
236 
237  // Storage for all triangles
238  List<labelledTri> triangles(nTotTri);
239  faceMap_.setSize(nTotTri);
240  label triI = 0;
241 
242 
243  if (faceCentreDecomposition)
244  {
245  // Decomposition around face centre
246  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
247 
248  // Triangulate internal faces
249  forAll(faceIsCut, facei)
250  {
251  if (faceIsCut[facei] && isInternalFace(mesh, includedCell, facei))
252  {
253  // Face was internal to the mesh and will be 'internal' to
254  // the surface.
255 
256  // Triangulate face
257  const face& f = faces[facei];
258 
259  forAll(f, fp)
260  {
261  faceMap_[triI] = facei;
262 
263  triangles[triI++] =
265  (
266  f[fp],
267  f.nextLabel(fp),
268  mesh.nPoints() + facei, // face centre
269  internalFacesPatch
270  );
271  }
272  }
273  }
274  nInternalFaces_ = triI;
275 
276 
277  // Triangulate external faces
278  forAll(faceIsCut, facei)
279  {
280  if (faceIsCut[facei] && !isInternalFace(mesh, includedCell, facei))
281  {
282  // Face will become outside of the surface.
283 
284  label patchi = -1;
285  bool reverse = false;
286 
287  if (mesh.isInternalFace(facei))
288  {
289  patchi = internalFacesPatch;
290 
291  // Check orientation. Check which side of the face gets
292  // included (note: only one side is).
293  if (includedCell[mesh.faceOwner()[facei]])
294  {
295  reverse = false;
296  }
297  else
298  {
299  reverse = true;
300  }
301  }
302  else
303  {
304  // Face was already outside so orientation ok.
305 
306  patchi = patches.whichPatch(facei);
307 
308  reverse = false;
309  }
310 
311 
312  // Triangulate face
313  const face& f = faces[facei];
314 
315  if (reverse)
316  {
317  forAll(f, fp)
318  {
319  faceMap_[triI] = facei;
320 
321  triangles[triI++] =
323  (
324  f.nextLabel(fp),
325  f[fp],
326  mesh.nPoints() + facei, // face centre
327  patchi
328  );
329  }
330  }
331  else
332  {
333  forAll(f, fp)
334  {
335  faceMap_[triI] = facei;
336 
337  triangles[triI++] =
339  (
340  f[fp],
341  f.nextLabel(fp),
342  mesh.nPoints() + facei, // face centre
343  patchi
344  );
345  }
346  }
347  }
348  }
349  }
350  else
351  {
352  // Triangulation using existing vertices
353  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
354 
355  // Triangulate internal faces
356  forAll(faceIsCut, facei)
357  {
358  if (faceIsCut[facei] && isInternalFace(mesh, includedCell, facei))
359  {
360  // Face was internal to the mesh and will be 'internal' to
361  // the surface.
362 
363  // Triangulate face. Fall back to naive triangulation if failed.
364  faceTriangulation faceTris(points, faces[facei], true);
365 
366  if (faceTris.empty())
367  {
369  << "Could not find triangulation for face " << facei
370  << " vertices " << faces[facei] << " coords "
371  << IndirectList<point>(points, faces[facei])() << endl;
372  }
373  else
374  {
375  // Copy triangles. Make them internalFacesPatch
376  insertTriangles
377  (
378  faceTris,
379  facei,
380  internalFacesPatch,
381  false, // no reverse
382 
383  triangles,
384  triI
385  );
386  }
387  }
388  }
389  nInternalFaces_ = triI;
390 
391 
392  // Triangulate external faces
393  forAll(faceIsCut, facei)
394  {
395  if (faceIsCut[facei] && !isInternalFace(mesh, includedCell, facei))
396  {
397  // Face will become outside of the surface.
398 
399  label patchi = -1;
400  bool reverse = false;
401 
402  if (mesh.isInternalFace(facei))
403  {
404  patchi = internalFacesPatch;
405 
406  // Check orientation. Check which side of the face gets
407  // included (note: only one side is).
408  if (includedCell[mesh.faceOwner()[facei]])
409  {
410  reverse = false;
411  }
412  else
413  {
414  reverse = true;
415  }
416  }
417  else
418  {
419  // Face was already outside so orientation ok.
420 
421  patchi = patches.whichPatch(facei);
422 
423  reverse = false;
424  }
425 
426  // Triangulate face
427  faceTriangulation faceTris(points, faces[facei], true);
428 
429  if (faceTris.empty())
430  {
432  << "Could not find triangulation for face " << facei
433  << " vertices " << faces[facei] << " coords "
434  << IndirectList<point>(points, faces[facei])() << endl;
435  }
436  else
437  {
438  // Copy triangles. Optionally reverse them
439  insertTriangles
440  (
441  faceTris,
442  facei,
443  patchi,
444  reverse, // whether to reverse
445 
446  triangles,
447  triI
448  );
449  }
450  }
451  }
452  }
453 
454  // Shrink if necessary (because of invalid triangulations)
455  triangles.setSize(triI);
456  faceMap_.setSize(triI);
457 
458  Pout<< "nInternalFaces_:" << nInternalFaces_ << endl;
459  Pout<< "triangles:" << triangles.size() << endl;
460 
461 
462  geometricSurfacePatchList surfPatches(patches.size());
463 
464  forAll(patches, patchi)
465  {
466  surfPatches[patchi] =
468  (
469  patches[patchi].physicalType(),
470  patches[patchi].name(),
471  patchi
472  );
473  }
474 
475  // Create globally numbered tri surface
476  if (faceCentreDecomposition)
477  {
478  // Use newPoints (mesh points + face centres)
479  triSurface globalSurf(triangles, surfPatches, newPoints);
480 
481  // Create locally numbered tri surface
483  (
484  triSurface
485  (
486  globalSurf.localFaces(),
487  surfPatches,
488  globalSurf.localPoints()
489  )
490  );
491  }
492  else
493  {
494  // Use mesh points
495  triSurface globalSurf(triangles, surfPatches, mesh.points());
496 
497  // Create locally numbered tri surface
499  (
500  triSurface
501  (
502  globalSurf.localFaces(),
503  surfPatches,
504  globalSurf.localPoints()
505  )
506  );
507  }
508 }
509 
510 
511 // ************************************************************************* //
Triangulation of faces. Handles concave polygons as well (inefficiently)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
friend Ostream & operator(Ostream &, const UList< T > &)
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:313
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const Field< point > & points() const
Return reference to global points.
const vectorField & faceCentres() const
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
List< triFace > triFaceList
list of triFaces
Definition: triFaceList.H:42
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:314
face triFace(3)
label nInternalFaces() const
Number of triangles in *this which are internal to the surface.
List< label > labelList
A List of labels.
Definition: labelList.H:56
Triangle with additional region number.
Definition: labelledTri.H:57
Foam::polyBoundaryMesh.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
meshTriangulation()
Construct null.
void setSize(const label)
Reset size of List.
Definition: List.C:295
triSurface()
Construct null.
Definition: triSurface.C:594
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
The geometricSurfacePatch is like patchIdentifier but for surfaces. Holds type, name and index...
label nPoints() const
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:117
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
Triangulated surface description with patch information.
Definition: triSurface.H:65
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
A List with indirect addressing.
Definition: IndirectList.H:102
const Field< PointType > & localPoints() const
Return pointField of points in patch.
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29