FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 
26  Note: bug in vtk displaying wedges? Seems to display ok if we decompose
27  them. Should be thoroughly tested!
28  (they appear rarely in polyhedral meshes, do appear in some cut meshes)
29 
30 \*---------------------------------------------------------------------------*/
31 
32 #include "vtkTopo.H"
33 #include <OpenFOAM/polyMesh.H>
34 #include <OpenFOAM/cellShape.H>
35 #include <OpenFOAM/cellModeller.H>
36 
37 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
38 
39 // Construct from components
40 Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
41 :
42  mesh_(mesh),
43  vertLabels_(),
44  cellTypes_(),
45  addPointCellLabels_(),
46  superCells_()
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& tetWedge = *(cellModeller::lookup("tetWedge"));
52  const cellModel& hex = *(cellModeller::lookup("hex"));
53 
54 
55  const cellShapeList& cellShapes = mesh_.cellShapes();
56 
57 
58  // Number of additional points needed by the decomposition of polyhedra
59  label nAddPoints = 0;
60 
61  // Number of additional cells generated by the decomposition of polyhedra
62  label nAddCells = 0;
63 
64  // Scan for cells which need to be decomposed and count additional points
65  // and cells
66 
67  forAll(cellShapes, cellI)
68  {
69  const cellModel& model = cellShapes[cellI].model();
70 
71  if
72  (
73  model != hex
74 // && model != wedge // See above.
75  && model != prism
76  && model != pyr
77  && model != tet
78  && model != tetWedge
79  )
80  {
81  const cell& cFaces = mesh_.cells()[cellI];
82 
83  forAll(cFaces, cFaceI)
84  {
85  const face& f = mesh_.faces()[cFaces[cFaceI]];
86 
87  label nQuads = 0;
88  label nTris = 0;
89  f.nTrianglesQuads(mesh_.points(), nTris, nQuads);
90 
91  nAddCells += nQuads + nTris;
92  }
93 
94  nAddCells--;
95  nAddPoints++;
96  }
97  }
98 
99  // Set size of additional point addressing array
100  // (from added point to original cell)
101  addPointCellLabels_.setSize(nAddPoints);
102 
103  // Set size of additional cells mapping array
104  // (from added cell to original cell)
105  superCells_.setSize(nAddCells);
106 
107  // List of vertex labels in VTK ordering
108  vertLabels_.setSize(cellShapes.size() + nAddCells);
109 
110  // Label of vtk type
111  cellTypes_.setSize(cellShapes.size() + nAddCells);
112 
113  // Set counters for additional points and additional cells
114  label api = 0, aci = 0;
115 
116  forAll(cellShapes, cellI)
117  {
118  const cellShape& cellShape = cellShapes[cellI];
119  const cellModel& cellModel = cellShape.model();
120 
121  labelList& vtkVerts = vertLabels_[cellI];
122 
123  if (cellModel == tet)
124  {
125  vtkVerts = cellShape;
126 
127  cellTypes_[cellI] = VTK_TETRA;
128  }
129  else if (cellModel == pyr)
130  {
131  vtkVerts = cellShape;
132 
133  cellTypes_[cellI] = VTK_PYRAMID;
134  }
135  else if (cellModel == prism)
136  {
137  vtkVerts.setSize(6);
138  vtkVerts[0] = cellShape[0];
139  vtkVerts[1] = cellShape[2];
140  vtkVerts[2] = cellShape[1];
141  vtkVerts[3] = cellShape[3];
142  vtkVerts[4] = cellShape[5];
143  vtkVerts[5] = cellShape[4];
144 
145  // VTK calls this a wedge.
146  cellTypes_[cellI] = VTK_WEDGE;
147  }
148  else if (cellModel == tetWedge)
149  {
150  // Treat as squeezed prism
151  vtkVerts.setSize(6);
152  vtkVerts[0] = cellShape[0];
153  vtkVerts[1] = cellShape[2];
154  vtkVerts[2] = cellShape[1];
155  vtkVerts[3] = cellShape[3];
156  vtkVerts[4] = cellShape[4];
157  vtkVerts[5] = cellShape[4];
158 
159  cellTypes_[cellI] = VTK_WEDGE;
160  }
161 // else if (cellModel == wedge)
162 // {
163 // // Treat as squeezed hex
164 // vtkVerts.setSize(8);
165 // vtkVerts[0] = cellShape[0];
166 // vtkVerts[1] = cellShape[1];
167 // vtkVerts[2] = cellShape[2];
168 // vtkVerts[3] = cellShape[0];
169 // vtkVerts[4] = cellShape[3];
170 // vtkVerts[5] = cellShape[4];
171 // vtkVerts[6] = cellShape[5];
172 // vtkVerts[7] = cellShape[6];
173 //
174 // cellTypes_[cellI] = VTK_HEXAHEDRON;
175 // }
176  else if (cellModel == hex)
177  {
178  vtkVerts.setSize(8);
179  vtkVerts[0] = cellShape[0];
180  vtkVerts[1] = cellShape[1];
181  vtkVerts[2] = cellShape[2];
182  vtkVerts[3] = cellShape[3];
183  vtkVerts[4] = cellShape[4];
184  vtkVerts[5] = cellShape[5];
185  vtkVerts[6] = cellShape[6];
186  vtkVerts[7] = cellShape[7];
187 
188  cellTypes_[cellI] = VTK_HEXAHEDRON;
189  }
190  else
191  {
192  // Polyhedral cell. Decompose into tets + prisms.
193  // (see dxFoamExec/createDxConnections.C)
194 
195  // Mapping from additional point to cell
196  addPointCellLabels_[api] = cellI;
197 
198  // Whether to insert cell in place of original or not.
199  bool substituteCell = true;
200 
201  const labelList& cFaces = mesh_.cells()[cellI];
202 
203  forAll(cFaces, cFaceI)
204  {
205  const face& f = mesh_.faces()[cFaces[cFaceI]];
206 
207  // Number of triangles and quads in decomposition
208  label nTris = 0;
209  label nQuads = 0;
210  f.nTrianglesQuads(mesh_.points(), nTris, nQuads);
211 
212  // Do actual decomposition into triFcs and quadFcs.
213  faceList triFcs(nTris);
214  faceList quadFcs(nQuads);
215  label trii = 0;
216  label quadi = 0;
217  f.trianglesQuads(mesh_.points(), trii, quadi, triFcs, quadFcs);
218 
219  forAll(quadFcs, quadi)
220  {
221  label thisCellI = -1;
222 
223  if (substituteCell)
224  {
225  thisCellI = cellI;
226 
227  substituteCell = false;
228  }
229  else
230  {
231  thisCellI = mesh_.nCells() + aci;
232 
233  superCells_[aci] = cellI;
234 
235  aci++;
236  }
237 
238  labelList& addVtkVerts = vertLabels_[thisCellI];
239 
240  addVtkVerts.setSize(5);
241 
242  const face& quad = quadFcs[quadi];
243 
244  addVtkVerts[0] = quad[0];
245  addVtkVerts[1] = quad[1];
246  addVtkVerts[2] = quad[2];
247  addVtkVerts[3] = quad[3];
248  addVtkVerts[4] = mesh_.nPoints() + api;
249 
250  cellTypes_[thisCellI] = VTK_PYRAMID;
251  }
252 
253  forAll(triFcs, trii)
254  {
255  label thisCellI = -1;
256 
257  if (substituteCell)
258  {
259  thisCellI = cellI;
260 
261  substituteCell = false;
262  }
263  else
264  {
265  thisCellI = mesh_.nCells() + aci;
266 
267  superCells_[aci] = cellI;
268 
269  aci++;
270  }
271 
272 
273  labelList& addVtkVerts = vertLabels_[thisCellI];
274 
275  const face& tri = triFcs[trii];
276 
277  addVtkVerts.setSize(4);
278  addVtkVerts[0] = tri[0];
279  addVtkVerts[1] = tri[1];
280  addVtkVerts[2] = tri[2];
281  addVtkVerts[3] = mesh_.nPoints() + api;
282 
283  cellTypes_[thisCellI] = VTK_TETRA;
284  }
285  }
286 
287  api++;
288  }
289  }
290 
291  Pout<< " Original cells:" << mesh_.nCells()
292  << " points:" << mesh_.nPoints()
293  << " Additional cells:" << superCells_.size()
294  << " additional points:" << addPointCellLabels_.size()
295  << nl << endl;
296 }
297 
298 // ************************ vim: set sw=4 sts=4 et: ************************ //