vtkPVFoamMeshVolume.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-2019 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 "vtkPVFoam.H"
27 #include "vtkPVFoamReader.h"
28 
29 // OpenFOAM includes
30 #include "fvMesh.H"
31 #include "cellModeller.H"
32 #include "vtkOpenFOAMPoints.H"
33 #include "Swap.H"
34 
35 // VTK includes
36 #include "vtkCellArray.h"
37 #include "vtkIdTypeArray.h"
38 #include "vtkUnstructuredGrid.h"
39 
40 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
41 
42 vtkUnstructuredGrid* Foam::vtkPVFoam::volumeVTKMesh
43 (
44  const fvMesh& mesh,
45  polyDecomp& decompInfo
46 )
47 {
48  const cellModel& tet = *(cellModeller::lookup("tet"));
49  const cellModel& pyr = *(cellModeller::lookup("pyr"));
50  const cellModel& prism = *(cellModeller::lookup("prism"));
51  const cellModel& wedge = *(cellModeller::lookup("wedge"));
52  const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
53  const cellModel& hex = *(cellModeller::lookup("hex"));
54 
55  vtkUnstructuredGrid* vtkmesh = vtkUnstructuredGrid::New();
56 
57  if (debug)
58  {
60  printMemory();
61  }
62 
63  const cellShapeList& cellShapes = mesh.cellShapes();
64 
65  // Number of additional points needed by the decomposition of polyhedra
66  label nAddPoints = 0;
67 
68  // Number of additional cells generated by the decomposition of polyhedra
69  label nAddCells = 0;
70 
71  // face owner is needed to determine the face orientation
72  const labelList& owner = mesh.faceOwner();
73 
74  labelList& superCells = decompInfo.superCells();
75  labelList& addPointCellLabels = decompInfo.addPointCellLabels();
76 
77  // Scan for cells which need to be decomposed and count additional points
78  // and cells
79  if (!reader_->GetUseVTKPolyhedron())
80  {
81  if (debug)
82  {
84  << "... scanning for polyhedra" << endl;
85  }
86 
87  forAll(cellShapes, celli)
88  {
89  const cellModel& model = cellShapes[celli].model();
90 
91  if
92  (
93  model != hex
94  && model != wedge
95  && model != prism
96  && model != pyr
97  && model != tet
98  && model != tetWedge
99  )
100  {
101  const cell& cFaces = mesh.cells()[celli];
102 
103  forAll(cFaces, cFacei)
104  {
105  const face& f = mesh.faces()[cFaces[cFacei]];
106 
107  label nQuads = 0;
108  label nTris = 0;
109  f.nTrianglesQuads(mesh.points(), nTris, nQuads);
110 
111  nAddCells += nQuads + nTris;
112  }
113 
114  nAddCells--;
115  nAddPoints++;
116  }
117  }
118  }
119 
120  // Set size of additional point addressing array
121  // (from added point to original cell)
122  addPointCellLabels.setSize(nAddPoints);
123 
124  // Set size of additional cells mapping array
125  // (from added cell to original cell)
126 
127  if (debug)
128  {
129  Info<< " mesh nCells = " << mesh.nCells() << nl
130  << " nPoints = " << mesh.nPoints() << nl
131  << " nAddCells = " << nAddCells << nl
132  << " nAddPoints = " << nAddPoints << endl;
133  }
134 
135  superCells.setSize(mesh.nCells() + nAddCells);
136 
137  if (debug)
138  {
139  Info<< " ... converting points" << endl;
140  }
141 
142  // Convert OpenFOAM mesh vertices to VTK
143  vtkPoints* vtkpoints = vtkPoints::New();
144  vtkpoints->Allocate(mesh.nPoints() + nAddPoints);
145 
146  const Foam::pointField& points = mesh.points();
147 
148  forAll(points, i)
149  {
150  vtkInsertNextOpenFOAMPoint(vtkpoints, points[i]);
151  }
152 
153 
154  if (debug)
155  {
156  Info<< " ... converting cells" << endl;
157  }
158 
159  vtkmesh->Allocate(mesh.nCells() + nAddCells);
160 
161  // Set counters for additional points and additional cells
162  label addPointi = 0, addCelli = 0;
163 
164  // Create storage for points - needed for mapping from OpenFOAM to VTK
165  // data types - max 'order' = hex = 8 points
166  vtkIdType nodeIds[8];
167 
168  // face-stream for a polyhedral cell
169  // [numFace0Pts, id1, id2, id3, numFace1Pts, id1, id2, id3, ...]
170  DynamicList<vtkIdType> faceStream(256);
171 
172  forAll(cellShapes, celli)
173  {
174  const cellShape& cellShape = cellShapes[celli];
175  const cellModel& cellModel = cellShape.model();
176 
177  superCells[addCelli++] = celli;
178 
179  if (cellModel == tet)
180  {
181  for (int j = 0; j < 4; j++)
182  {
183  nodeIds[j] = cellShape[j];
184  }
185  vtkmesh->InsertNextCell
186  (
187  VTK_TETRA,
188  4,
189  nodeIds
190  );
191  }
192  else if (cellModel == pyr)
193  {
194  for (int j = 0; j < 5; j++)
195  {
196  nodeIds[j] = cellShape[j];
197  }
198  vtkmesh->InsertNextCell
199  (
200  VTK_PYRAMID,
201  5,
202  nodeIds
203  );
204  }
205  else if (cellModel == prism)
206  {
207  // VTK has a different node order for VTK_WEDGE
208  // their triangles point outwards!
209  nodeIds[0] = cellShape[0];
210  nodeIds[1] = cellShape[2];
211  nodeIds[2] = cellShape[1];
212  nodeIds[3] = cellShape[3];
213  nodeIds[4] = cellShape[5];
214  nodeIds[5] = cellShape[4];
215 
216  vtkmesh->InsertNextCell
217  (
218  VTK_WEDGE,
219  6,
220  nodeIds
221  );
222  }
223  else if (cellModel == tetWedge && !reader_->GetUseVTKPolyhedron())
224  {
225  // Treat as squeezed prism (VTK_WEDGE)
226 
227  nodeIds[0] = cellShape[0];
228  nodeIds[1] = cellShape[2];
229  nodeIds[2] = cellShape[1];
230  nodeIds[3] = cellShape[3];
231  nodeIds[4] = cellShape[4];
232  nodeIds[5] = cellShape[3];
233 
234  vtkmesh->InsertNextCell
235  (
236  VTK_WEDGE,
237  6,
238  nodeIds
239  );
240  }
241  else if (cellModel == wedge)
242  {
243  // Treat as squeezed hex
244 
245  nodeIds[0] = cellShape[0];
246  nodeIds[1] = cellShape[1];
247  nodeIds[2] = cellShape[2];
248  nodeIds[3] = cellShape[2];
249  nodeIds[4] = cellShape[3];
250  nodeIds[5] = cellShape[4];
251  nodeIds[6] = cellShape[5];
252  nodeIds[7] = cellShape[6];
253 
254  vtkmesh->InsertNextCell
255  (
256  VTK_HEXAHEDRON,
257  8,
258  nodeIds
259  );
260  }
261  else if (cellModel == hex)
262  {
263  for (int j = 0; j < 8; j++)
264  {
265  nodeIds[j] = cellShape[j];
266  }
267  vtkmesh->InsertNextCell
268  (
269  VTK_HEXAHEDRON,
270  8,
271  nodeIds
272  );
273  }
274  else if (reader_->GetUseVTKPolyhedron())
275  {
276  // Polyhedral cell - use VTK_POLYHEDRON
277  const labelList& cFaces = mesh.cells()[celli];
278 
279 #ifdef HAS_VTK_POLYHEDRON
280  vtkIdType nFaces = cFaces.size();
281  vtkIdType nLabels = nFaces;
282 
283  // count size for face stream
284  forAll(cFaces, cFacei)
285  {
286  const face& f = mesh.faces()[cFaces[cFacei]];
287  nLabels += f.size();
288  }
289 
290  // build face-stream
291  // [numFace0Pts, id1, id2, id3, numFace1Pts, id1, id2, id3, ...]
292  // point Ids are global
293  faceStream.clear();
294  faceStream.reserve(nLabels + nFaces);
295 
296  forAll(cFaces, cFacei)
297  {
298  const face& f = mesh.faces()[cFaces[cFacei]];
299  const bool isOwner = (owner[cFaces[cFacei]] == celli);
300  const label nFacePoints = f.size();
301 
302  // number of labels for this face
303  faceStream.append(nFacePoints);
304 
305  if (isOwner)
306  {
307  forAll(f, fp)
308  {
309  faceStream.append(f[fp]);
310  }
311  }
312  else
313  {
314  // Fairly immaterial if we reverse the list
315  // or use face::reverseFace()
316  forAllReverse(f, fp)
317  {
318  faceStream.append(f[fp]);
319  }
320  }
321  }
322 
323  vtkmesh->InsertNextCell(VTK_POLYHEDRON, nFaces, faceStream.data());
324 #else
325  // This is a horrible substitute
326  // but avoids crashes when there is no vtkPolyhedron support
327 
328  // Establish unique node ids used
329  HashSet<vtkIdType, Hash<label>> hashUniqId(2*256);
330 
331  forAll(cFaces, cFacei)
332  {
333  const face& f = mesh.faces()[cFaces[cFacei]];
334 
335  forAll(f, fp)
336  {
337  hashUniqId.insert(f[fp]);
338  }
339  }
340 
341  // Use face stream to store unique node ids:
342  faceStream = hashUniqId.sortedToc();
343 
344  vtkmesh->InsertNextCell
345  (
346  VTK_CONVEX_POINT_SET,
347  vtkIdType(faceStream.size()),
348  faceStream.data()
349  );
350 #endif
351  }
352  else
353  {
354  // Polyhedral cell. Decompose into tets + prisms.
355 
356  // Mapping from additional point to cell
357  addPointCellLabels[addPointi] = celli;
358 
359  // The new vertex from the cell-centre
360  const label newVertexLabel = mesh.nPoints() + addPointi;
361  vtkInsertNextOpenFOAMPoint(vtkpoints, mesh.C()[celli]);
362 
363  // Whether to insert cell in place of original or not.
364  bool substituteCell = true;
365 
366  const labelList& cFaces = mesh.cells()[celli];
367  forAll(cFaces, cFacei)
368  {
369  const face& f = mesh.faces()[cFaces[cFacei]];
370  const bool isOwner = (owner[cFaces[cFacei]] == celli);
371 
372  // Number of triangles and quads in decomposition
373  label nTris = 0;
374  label nQuads = 0;
375  f.nTrianglesQuads(mesh.points(), nTris, nQuads);
376 
377  // Do actual decomposition into triFcs and quadFcs.
378  faceList triFcs(nTris);
379  faceList quadFcs(nQuads);
380  label trii = 0;
381  label quadi = 0;
382  f.trianglesQuads(mesh.points(), trii, quadi, triFcs, quadFcs);
383 
384  forAll(quadFcs, quadI)
385  {
386  if (substituteCell)
387  {
388  substituteCell = false;
389  }
390  else
391  {
392  superCells[addCelli++] = celli;
393  }
394 
395  const face& quad = quadFcs[quadI];
396 
397  // Ensure we have the correct orientation for the
398  // base of the primitive cell shape.
399  // If the cell is face owner, the orientation needs to be
400  // flipped.
401  // At the moment, VTK doesn't actually seem to care if
402  // negative cells are defined, but we'll do it anyhow
403  // (for safety).
404  if (isOwner)
405  {
406  nodeIds[0] = quad[3];
407  nodeIds[1] = quad[2];
408  nodeIds[2] = quad[1];
409  nodeIds[3] = quad[0];
410  }
411  else
412  {
413  nodeIds[0] = quad[0];
414  nodeIds[1] = quad[1];
415  nodeIds[2] = quad[2];
416  nodeIds[3] = quad[3];
417  }
418  nodeIds[4] = newVertexLabel;
419  vtkmesh->InsertNextCell
420  (
421  VTK_PYRAMID,
422  5,
423  nodeIds
424  );
425  }
426 
427  forAll(triFcs, triI)
428  {
429  if (substituteCell)
430  {
431  substituteCell = false;
432  }
433  else
434  {
435  superCells[addCelli++] = celli;
436  }
437 
438  const face& tri = triFcs[triI];
439 
440  // See note above about the orientation.
441  if (isOwner)
442  {
443  nodeIds[0] = tri[2];
444  nodeIds[1] = tri[1];
445  nodeIds[2] = tri[0];
446  }
447  else
448  {
449  nodeIds[0] = tri[0];
450  nodeIds[1] = tri[1];
451  nodeIds[2] = tri[2];
452  }
453  nodeIds[3] = newVertexLabel;
454 
455  vtkmesh->InsertNextCell
456  (
457  VTK_TETRA,
458  4,
459  nodeIds
460  );
461  }
462  }
463 
464  addPointi++;
465  }
466  }
467 
468  vtkmesh->SetPoints(vtkpoints);
469  vtkpoints->Delete();
470 
471  if (debug)
472  {
473  printMemory();
474  }
475 
476  return vtkmesh;
477 }
478 
479 
480 // ************************************************************************* //
#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
static void printMemory()
Simple memory used debugging information.
void vtkInsertNextOpenFOAMPoint(vtkPoints *points, const Foam::point &p)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
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
Pre-declare SubField and related Field type.
Definition: Field.H:56
autoPtr< BasicCompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleMomentumTransportModel::transportModel &transport)
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const char nl
Definition: Ostream.H:260
Swap its arguments.
messageStream Info
#define InfoInFunction
Report an information message using Foam::Info.