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