All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
vtkTopo.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 Description
25  Note: bug in vtk displaying wedges? Seems to display ok if we decompose
26  them. Should be thoroughly tested!
27  (they appear rarely in polyhedral meshes, do appear in some cut meshes)
28 
29 \*---------------------------------------------------------------------------*/
30 
31 #include "vtkTopo.H"
32 #include "polyMesh.H"
33 #include "cellShape.H"
34 #include "cellModeller.H"
35 #include "Swap.H"
36 #include "polygonTriangulate.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
41 
42 
43 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
44 
45 Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
46 :
47  mesh_(mesh),
48  vertLabels_(),
49  cellTypes_(),
50  addPointCellLabels_(),
51  superCells_()
52 {
53  const cellModel& tet = *(cellModeller::lookup("tet"));
54  const cellModel& pyr = *(cellModeller::lookup("pyr"));
55  const cellModel& prism = *(cellModeller::lookup("prism"));
56  const cellModel& wedge = *(cellModeller::lookup("wedge"));
57  const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
58  const cellModel& hex = *(cellModeller::lookup("hex"));
59 
60  const cellShapeList& cellShapes = mesh_.cellShapes();
61 
62  // Number of additional points needed by the decomposition of polyhedra
63  label nAddPoints = 0;
64 
65  // Number of additional cells generated by the decomposition of polyhedra
66  label nAddCells = 0;
67 
68  // face owner is needed to determine the face orientation
69  const labelList& owner = mesh.faceOwner();
70 
71  // Scan for cells which need to be decomposed and count additional points
72  // and cells
73  if (decomposePoly)
74  {
75  forAll(cellShapes, celli)
76  {
77  const cellModel& model = cellShapes[celli].model();
78 
79  if
80  (
81  model != hex
82  && model != wedge // See above.
83  && model != prism
84  && model != pyr
85  && model != tet
86  && model != tetWedge
87  )
88  {
89  const cell& cFaces = mesh_.cells()[celli];
90 
91  forAll(cFaces, cFacei)
92  {
93  const face& f = mesh_.faces()[cFaces[cFacei]];
94 
95  if (f.size() != 4)
96  {
97  nAddCells += f.nTriangles();
98  }
99  else
100  {
101  nAddCells ++;
102  }
103  }
104 
105  nAddCells--;
106  nAddPoints++;
107  }
108  }
109  }
110 
111 
112  // Set size of additional point addressing array
113  // (from added point to original cell)
114  addPointCellLabels_.setSize(nAddPoints);
115 
116  // Set size of additional cells mapping array
117  // (from added cell to original cell)
118  superCells_.setSize(nAddCells);
119 
120  // List of vertex labels in VTK ordering
121  vertLabels_.setSize(cellShapes.size() + nAddCells);
122 
123  // Label of vtk type
124  cellTypes_.setSize(cellShapes.size() + nAddCells);
125 
126  // Create a triangulation engine
127  polygonTriangulate triEngine;
128 
129  // Set counters for additional points and additional cells
130  label addPointi = 0, addCelli = 0;
131 
132  forAll(cellShapes, celli)
133  {
134  const cellShape& cellShape = cellShapes[celli];
135  const cellModel& cellModel = cellShape.model();
136 
137  labelList& vtkVerts = vertLabels_[celli];
138 
139  if (cellModel == tet)
140  {
141  vtkVerts = cellShape;
142 
143  cellTypes_[celli] = VTK_TETRA;
144  }
145  else if (cellModel == pyr)
146  {
147  vtkVerts = cellShape;
148 
149  cellTypes_[celli] = VTK_PYRAMID;
150  }
151  else if (cellModel == prism)
152  {
153  // VTK has a different node order for VTK_WEDGE
154  // their triangles point outwards!
155  vtkVerts = cellShape;
156 
157  Foam::Swap(vtkVerts[1], vtkVerts[2]);
158  Foam::Swap(vtkVerts[4], vtkVerts[5]);
159 
160  cellTypes_[celli] = VTK_WEDGE;
161  }
162  else if (cellModel == tetWedge && decomposePoly)
163  {
164  // Treat as squeezed prism (VTK_WEDGE)
165  vtkVerts.setSize(6);
166  vtkVerts[0] = cellShape[0];
167  vtkVerts[1] = cellShape[2];
168  vtkVerts[2] = cellShape[1];
169  vtkVerts[3] = cellShape[3];
170  vtkVerts[4] = cellShape[4];
171  vtkVerts[5] = cellShape[3];
172 
173  cellTypes_[celli] = VTK_WEDGE;
174  }
175  else if (cellModel == wedge)
176  {
177  // Treat as squeezed hex
178  vtkVerts.setSize(8);
179  vtkVerts[0] = cellShape[0];
180  vtkVerts[1] = cellShape[1];
181  vtkVerts[2] = cellShape[2];
182  vtkVerts[3] = cellShape[2];
183  vtkVerts[4] = cellShape[3];
184  vtkVerts[5] = cellShape[4];
185  vtkVerts[6] = cellShape[5];
186  vtkVerts[7] = cellShape[6];
187 
188  cellTypes_[celli] = VTK_HEXAHEDRON;
189  }
190  else if (cellModel == hex)
191  {
192  vtkVerts = cellShape;
193 
194  cellTypes_[celli] = VTK_HEXAHEDRON;
195  }
196  else if (decomposePoly)
197  {
198  // Polyhedral cell. Decompose into tets + pyramids.
199 
200  // Mapping from additional point to cell
201  addPointCellLabels_[addPointi] = celli;
202 
203  // The new vertex from the cell-centre
204  const label newVertexLabel = mesh_.nPoints() + addPointi;
205 
206  // Whether to insert cell in place of original or not.
207  bool substituteCell = true;
208 
209  const labelList& cFaces = mesh_.cells()[celli];
210  forAll(cFaces, cFacei)
211  {
212  const face& f = mesh_.faces()[cFaces[cFacei]];
213  const bool isOwner = (owner[cFaces[cFacei]] == celli);
214 
215  // Do decomposition into triFcs and quadFcs.
216  faceList triFcs(f.size() == 4 ? 0 : f.nTriangles());
217  faceList quadFcs(f.size() == 4 ? 1 : 0);
218  if (f.size() != 4)
219  {
220  triEngine.triangulate
221  (
222  UIndirectList<point>(mesh.points(), f)
223  );
224  forAll(triEngine.triPoints(), trii)
225  {
226  triFcs[trii] = triEngine.triPoints(trii, f);
227  }
228  }
229  else
230  {
231  quadFcs[0] = f;
232  }
233 
234  forAll(quadFcs, quadI)
235  {
236  label thisCelli;
237 
238  if (substituteCell)
239  {
240  thisCelli = celli;
241  substituteCell = false;
242  }
243  else
244  {
245  thisCelli = mesh_.nCells() + addCelli;
246  superCells_[addCelli++] = celli;
247  }
248 
249  labelList& addVtkVerts = vertLabels_[thisCelli];
250 
251  addVtkVerts.setSize(5);
252 
253  const face& quad = quadFcs[quadI];
254 
255  // Ensure we have the correct orientation for the
256  // base of the primitive cell shape.
257  // If the cell is face owner, the orientation needs to be
258  // flipped.
259  // At the moment, VTK doesn't actually seem to care if
260  // negative cells are defined, but we'll do it anyhow
261  // (for safety).
262  if (isOwner)
263  {
264  addVtkVerts[0] = quad[3];
265  addVtkVerts[1] = quad[2];
266  addVtkVerts[2] = quad[1];
267  addVtkVerts[3] = quad[0];
268  }
269  else
270  {
271  addVtkVerts[0] = quad[0];
272  addVtkVerts[1] = quad[1];
273  addVtkVerts[2] = quad[2];
274  addVtkVerts[3] = quad[3];
275  }
276  addVtkVerts[4] = newVertexLabel;
277 
278  cellTypes_[thisCelli] = VTK_PYRAMID;
279  }
280 
281  forAll(triFcs, triI)
282  {
283  label thisCelli;
284 
285  if (substituteCell)
286  {
287  thisCelli = celli;
288  substituteCell = false;
289  }
290  else
291  {
292  thisCelli = mesh_.nCells() + addCelli;
293  superCells_[addCelli++] = celli;
294  }
295 
296 
297  labelList& addVtkVerts = vertLabels_[thisCelli];
298 
299  const face& tri = triFcs[triI];
300 
301  addVtkVerts.setSize(4);
302 
303  // See note above about the orientation.
304  if (isOwner)
305  {
306  addVtkVerts[0] = tri[2];
307  addVtkVerts[1] = tri[1];
308  addVtkVerts[2] = tri[0];
309  }
310  else
311  {
312  addVtkVerts[0] = tri[0];
313  addVtkVerts[1] = tri[1];
314  addVtkVerts[2] = tri[2];
315  }
316  addVtkVerts[3] = newVertexLabel;
317 
318  cellTypes_[thisCelli] = VTK_TETRA;
319  }
320  }
321 
322  addPointi++;
323  }
324  else
325  {
326  // Polyhedral cell - not decomposed
327  cellTypes_[celli] = VTK_POLYHEDRON;
328 
329  const labelList& cFaces = mesh_.cells()[celli];
330 
331  // space for the number of faces and size of each face
332  label nData = 1 + cFaces.size();
333 
334  // count total number of face points
335  forAll(cFaces, cFacei)
336  {
337  const face& f = mesh.faces()[cFaces[cFacei]];
338  nData += f.size(); // space for the face labels
339  }
340 
341  vtkVerts.setSize(nData);
342 
343  nData = 0;
344  vtkVerts[nData++] = cFaces.size();
345 
346  // build face stream
347  forAll(cFaces, cFacei)
348  {
349  const face& f = mesh.faces()[cFaces[cFacei]];
350  const bool isOwner = (owner[cFaces[cFacei]] == celli);
351 
352  // number of labels for this face
353  vtkVerts[nData++] = f.size();
354 
355  if (isOwner)
356  {
357  forAll(f, fp)
358  {
359  vtkVerts[nData++] = f[fp];
360  }
361  }
362  else
363  {
364  // fairly immaterial if we reverse the list
365  // or use face::reverseFace()
366  forAllReverse(f, fp)
367  {
368  vtkVerts[nData++] = f[fp];
369  }
370  }
371  }
372  }
373  }
374 
375  if (decomposePoly)
376  {
377  Pout<< " Original cells:" << mesh_.nCells()
378  << " points:" << mesh_.nPoints()
379  << " Additional cells:" << superCells_.size()
380  << " additional points:" << addPointCellLabels_.size()
381  << nl << endl;
382  }
383 
384 }
385 
386 
387 // ************************************************************************* //
#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
static const cellModel * lookup(const word &)
Look up a model by name and return a pointer to the model or nullptr.
Definition: cellModeller.C:100
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:446
List< cellShape > cellShapeList
List of cellShapes and PtrList of List of cellShape.
Definition: cellShapeList.H:43
void Swap(T &a, T &b)
Definition: Swap.H:43
static bool decomposePoly
Enable/disable polyhedron decomposition. Default = true.
Definition: vtkTopo.H:87
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const char nl
Definition: Ostream.H:260
labelList f(nPoints)
vtkTopo(const polyMesh &)
Construct from components.
void setSize(const label)
Reset size of List.
Definition: List.C:281
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Swap its arguments.