FreeFOAM The Cross-Platform CFD Toolkit
ensightPartCells.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 \*----------------------------------------------------------------------------*/
25 
26 #include "ensightPartCells.H"
28 #include <OpenFOAM/IOstream.H>
29 #include <OpenFOAM/IStringStream.H>
30 #include <OpenFOAM/dictionary.H>
31 #include <OpenFOAM/cellModeller.H>
32 
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(ensightPartCells, 0);
39  addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream);
40 }
41 
43 (
44  IStringStream
45  (
46  "(tetra4 pyramid5 penta6 hexa8 nfaced)"
47  )()
48 );
49 
50 
51 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
52 
53 void Foam::ensightPartCells::classify(const labelList& idList)
54 {
55  // References to cell shape models
56  const cellModel& tet = *(cellModeller::lookup("tet"));
57  const cellModel& pyr = *(cellModeller::lookup("pyr"));
58  const cellModel& prism = *(cellModeller::lookup("prism"));
59  const cellModel& hex = *(cellModeller::lookup("hex"));
60 
61  const polyMesh& mesh = *meshPtr_;
62  const cellShapeList& cellShapes = mesh.cellShapes();
63 
64  offset_ = 0;
65  size_ = mesh.nCells();
66 
67  bool limited = false;
68  if (&idList)
69  {
70  limited = true;
71  size_ = idList.size();
72  }
73 
74  // count the shapes
75  label nTet = 0;
76  label nPyr = 0;
77  label nPrism = 0;
78  label nHex = 0;
79  label nPoly = 0;
80 
81 
82  // TODO: allow tet-decomposition of polyhedral cells
83 #if 0
84  label nTetDecomp = 0;
85  label nPyrDecomp = 0;
86 #endif
87 
88  for (label listI = 0; listI < size_; ++listI)
89  {
90  label cellId = listI;
91  if (limited)
92  {
93  cellId = idList[listI];
94  }
95 
96  const cellShape& cellShape = cellShapes[cellId];
97  const cellModel& cellModel = cellShape.model();
98 
99  if (cellModel == tet)
100  {
101  nTet++;
102  }
103  else if (cellModel == pyr)
104  {
105  nPyr++;
106  }
107  else if (cellModel == prism)
108  {
109  nPrism++;
110  }
111  else if (cellModel == hex)
112  {
113  nHex++;
114  }
115  else
116  {
117  nPoly++;
118 
119  // TODO: allow tet-decomposition of polyhedral cells
120 #if 0
121  const cell& cFaces = mesh.cells()[cellI];
122 
123  forAll(cFaces, cFaceI)
124  {
125  const face& f = mesh.faces()[cFaces[cFaceI]];
126 
127  label nQuads = 0;
128  label nTris = 0;
129  f.nTrianglesQuads(mesh.points(), nTris, nQuads);
130 
131  nTetDecomp += nTris;
132  nPyrDecomp += nQuads;
133  }
134 
135  nAddCells--;
136  nAddPoints++;
137 #endif
138  }
139  }
140 
141 
142  // we can avoid double looping, but at the cost of allocation
143  labelList tetCells(nTet);
144  labelList pyramidCells(nPyr);
145  labelList prismCells(nPrism);
146  labelList hexCells(nHex);
147  labelList polyCells(nPoly);
148 
149  nTet = 0,
150  nPyr = 0;
151  nPrism = 0;
152  nHex = 0;
153  nPoly = 0;
154 
155  // classify the shapes
156  for (label listI = 0; listI < size_; ++listI)
157  {
158  label cellId = listI;
159  if (limited)
160  {
161  cellId = idList[listI];
162  }
163 
164  const cellShape& cellShape = cellShapes[cellId];
165  const cellModel& cellModel = cellShape.model();
166 
167  if (cellModel == tet)
168  {
169  tetCells[nTet++] = cellId;
170  }
171  else if (cellModel == pyr)
172  {
173  pyramidCells[nPyr++] = cellId;
174  }
175  else if (cellModel == prism)
176  {
177  prismCells[nPrism++] = cellId;
178  }
179  else if (cellModel == hex)
180  {
181  hexCells[nHex++] = cellId;
182  }
183  else
184  {
185  polyCells[nPoly++] = cellId;
186 
187  // TODO: allow tet-decomposition of polyhedral cells
188 #if 0
189  // Mapping from additional point to cell
190  addPointCellLabels_[api] = cellId;
191 
192  const cell& cFaces = mesh.cells()[cellId];
193 
194  forAll(cFaces, cFaceI)
195  {
196  const face& f = mesh.faces()[cFaces[cFaceI]];
197 
198  label nQuads = 0;
199  label nTris = 0;
200  f.nTrianglesQuads(mesh.points(), nTris, nQuads);
201 
202  nTetDecomp += nTris;
203  nPyrDecomp += nQuads;
204  }
205 
206  nAddCells--;
207  nAddPoints++;
208 #endif
209  }
210  }
211 
212 
213  // MUST match with elementTypes
215 
216  elemLists_[tetra4Elements].transfer( tetCells );
217  elemLists_[pyramid5Elements].transfer( pyramidCells );
218  elemLists_[penta6Elements].transfer( prismCells );
219  elemLists_[hexa8Elements].transfer( hexCells );
220  elemLists_[nfacedElements].transfer( polyCells );
221 }
222 
223 
224 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
225 
227 (
228  label partNumber,
229  const string& partDescription
230 )
231 :
232  ensightPart(partNumber, partDescription)
233 {}
234 
235 
237 (
238  label partNumber,
239  const polyMesh& pMesh
240 )
241 :
242  ensightPart(partNumber, "cells", pMesh)
243 {
244  classify();
245 }
246 
247 
249 (
250  label partNumber,
251  const polyMesh& pMesh,
252  const labelList& idList
253 )
254 :
255  ensightPart(partNumber, "cells", pMesh)
256 {
257  classify(idList);
258 }
259 
260 
262 (
263  label partNumber,
264  const polyMesh& pMesh,
265  const cellZone& cZone
266 )
267 :
268  ensightPart(partNumber, cZone.name(), pMesh)
269 {
270  classify(cZone);
271 }
272 
273 
275 :
276  ensightPart(part)
277 {}
278 
279 
281 :
282  ensightPart()
283 {
284  reconstruct(is);
285 }
286 
287 
288 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
289 
291 {}
292 
293 
294 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
295 
296 Foam::ensightPart::localPoints Foam::ensightPartCells::calcLocalPoints() const
297 {
298  const polyMesh& mesh = *meshPtr_;
299 
300  localPoints ptList(mesh);
301  labelList& usedPoints = ptList.list;
302  label nPoints = 0;
303 
304  forAll(elemLists_, typeI)
305  {
306  const labelList& idList = elemLists_[typeI];
307 
308  // add all points from cells
309  forAll(idList, i)
310  {
311  label id = idList[i] + offset_;
312  const labelList& cFaces = mesh.cells()[id];
313 
314  forAll(cFaces, cFaceI)
315  {
316  const face& f = mesh.faces()[cFaces[cFaceI]];
317 
318  forAll(f, fp)
319  {
320  if (usedPoints[f[fp]] == -1)
321  {
322  usedPoints[f[fp]] = nPoints++;
323  }
324  }
325  }
326  }
327  }
328 
329  // this is not absolutely necessary, but renumber anyhow
330  nPoints = 0;
331  forAll(usedPoints, ptI)
332  {
333  if (usedPoints[ptI] > -1)
334  {
335  usedPoints[ptI] = nPoints++;
336  }
337  }
338 
339  ptList.nPoints = nPoints;
340  return ptList;
341 }
342 
343 
344 void Foam::ensightPartCells::writeConnectivity
345 (
346  ensightGeoFile& os,
347  const string& key,
348  const labelList& idList,
349  const labelList& pointMap
350 ) const
351 {
352  os.writeKeyword(key);
353  os.write(idList.size());
354  os.newline();
355 
356  const polyMesh& mesh = *meshPtr_;
357 
358  // write polyhedral
359  if (word(key) == "nfaced")
360  {
361  const faceList& meshFaces = mesh.faces();
362 
363  // write the number of faces per element
364  forAll(idList, i)
365  {
366  label id = idList[i] + offset_;
367  const labelList& cFace = mesh.cells()[id];
368 
369  os.write( cFace.size() );
370  os.newline();
371  }
372 
373  // write the number of points per element face
374  forAll(idList, i)
375  {
376  label id = idList[i] + offset_;
377  const labelList& cFace = mesh.cells()[id];
378 
379  forAll(cFace, faceI)
380  {
381  const face& cf = meshFaces[cFace[faceI]];
382 
383  os.write( cf.size() );
384  os.newline();
385  }
386  }
387 
388  // write the points describing each element face
389  forAll(idList, i)
390  {
391  label id = idList[i] + offset_;
392  const labelList& cFace = mesh.cells()[id];
393 
394  forAll(cFace, faceI)
395  {
396  const face& cf = meshFaces[cFace[faceI]];
397 
398  forAll(cf, ptI)
399  {
400  // convert global -> local index
401  // (note: Ensight indices start with 1)
402  os.write( pointMap[cf[ptI]] + 1);
403  }
404  os.newline();
405  }
406  }
407  }
408  else
409  {
410  // write primitive
411  const cellShapeList& cellShapes = mesh.cellShapes();
412 
413  forAll(idList, i)
414  {
415  label id = idList[i] + offset_;
416  const cellShape& cellPoints = cellShapes[id];
417 
418  // convert global -> local index
419  // (note: Ensight indices start with 1)
420  forAll(cellPoints, ptI)
421  {
422  os.write( pointMap[cellPoints[ptI]] + 1 );
423  }
424  os.newline();
425  }
426  }
427 }
428 
429 
430 // ************************ vim: set sw=4 sts=4 et: ************************ //