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