meshTriangulation.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-2021 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 "polygonTriangulate.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();
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  // Create a triangulation engine
356  polygonTriangulate triEngine;
357 
358  // Triangulate internal faces
359  forAll(faceIsCut, facei)
360  {
361  if (faceIsCut[facei] && isInternalFace(mesh, includedCell, facei))
362  {
363  triEngine.triangulate
364  (
365  UIndirectList<point>(points, faces[facei])
366  );
367 
368  insertTriangles
369  (
370  triEngine.triPoints(faces[facei]),
371  facei,
372  internalFacesPatch,
373  false, // no reverse
374  triangles,
375  triI
376  );
377  }
378  }
379  nInternalFaces_ = triI;
380 
381 
382  // Triangulate external faces
383  forAll(faceIsCut, facei)
384  {
385  if (faceIsCut[facei] && !isInternalFace(mesh, includedCell, facei))
386  {
387  // Face will become outside of the surface.
388 
389  label patchi = -1;
390  bool reverse = false;
391 
392  if (mesh.isInternalFace(facei))
393  {
394  patchi = internalFacesPatch;
395 
396  // Check orientation. Check which side of the face gets
397  // included (note: only one side is).
398  if (includedCell[mesh.faceOwner()[facei]])
399  {
400  reverse = false;
401  }
402  else
403  {
404  reverse = true;
405  }
406  }
407  else
408  {
409  // Face was already outside so orientation ok.
410 
411  patchi = patches.whichPatch(facei);
412 
413  reverse = false;
414  }
415 
416  triEngine.triangulate
417  (
418  UIndirectList<point>(points, faces[facei])
419  );
420 
421  insertTriangles
422  (
423  triEngine.triPoints(faces[facei]),
424  facei,
425  patchi,
426  reverse, // whether to reverse
427  triangles,
428  triI
429  );
430  }
431  }
432  }
433 
434  // Shrink if necessary (because of invalid triangulations)
435  triangles.setSize(triI);
436  faceMap_.setSize(triI);
437 
438  Pout<< "nInternalFaces_:" << nInternalFaces_ << endl;
439  Pout<< "triangles:" << triangles.size() << endl;
440 
441 
442  geometricSurfacePatchList surfPatches(patches.size());
443 
444  forAll(patches, patchi)
445  {
446  surfPatches[patchi] =
448  (
449  patches[patchi].physicalType(),
450  patches[patchi].name(),
451  patchi
452  );
453  }
454 
455  // Create globally numbered tri surface
456  if (faceCentreDecomposition)
457  {
458  // Use newPoints (mesh points + face centres)
459  triSurface globalSurf(triangles, surfPatches, newPoints);
460 
461  // Create locally numbered tri surface
463  (
464  triSurface
465  (
466  globalSurf.localFaces(),
467  surfPatches,
468  globalSurf.localPoints()
469  )
470  );
471  }
472  else
473  {
474  // Use mesh points
475  triSurface globalSurf(triangles, surfPatches, mesh.points());
476 
477  // Create locally numbered tri surface
479  (
480  triSurface
481  (
482  globalSurf.localFaces(),
483  surfPatches,
484  globalSurf.localPoints()
485  )
486  );
487  }
488 }
489 
490 
491 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:100
const Field< PointType > & localPoints() const
Return pointField of points in patch.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
List< triFace > triFaceList
list of triFaces
Definition: triFaceList.H:42
face triFace(3)
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
const Field< PointType > & points() const
Return reference to global points.
List< label > labelList
A List of labels.
Definition: labelList.H:56
Triangle with additional region number.
Definition: labelledTri.H:57
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
Foam::polyBoundaryMesh.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
void reverse(UList< T > &, const label n)
Definition: UListI.H:334
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
meshTriangulation()
Construct null.
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
const vectorField & faceCentres() const
void setSize(const label)
Reset size of List.
Definition: List.C:281
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:322
triSurface()
Construct null.
Definition: triSurface.C:534
label patchi
Triangulation of three-dimensional polygons.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
A List with indirect addressing.
Definition: fvMatrix.H:106
The geometricSurfacePatch is like patchIdentifier but for surfaces. Holds type, name and index...
label nPoints() const
faceList faces() const
Return the list of triangles as a faceList.
Definition: triSurface.C:1131
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Triangulated surface description with patch information.
Definition: triSurface.H:66
label nInternalFaces() const
Number of triangles in *this which are internal to the surface.
const UList< triFace > & triPoints() const
Get the triangles&#39; points.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.