FreeFOAM The Cross-Platform CFD Toolkit
fieldviewTopology.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 \*---------------------------------------------------------------------------*/
27 
28 #include "fieldviewTopology.H"
29 #include <OpenFOAM/polyMesh.H>
30 #include <OpenFOAM/cellShape.H>
31 #include <OpenFOAM/cellModeller.H>
32 #include <OpenFOAM/wallPolyPatch.H>
34 
35 
36 #include "fv_reader_tags.h"
37 
38 extern "C"
39 {
40  unsigned int fv_encode_elem_header(int elem_type, int wall_info[]);
41 }
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
46 Foam::labelList Foam::fieldviewTopology::calcFaceAddressing
47 (
48  const faceList& allFaces, // faces given faceLabels
49  const cellShape& shape,
50  const labelList& faces, // faceLabels for given cell
51  const label cellI
52 )
53 {
54  // return value.
55  labelList shapeToMesh(shape.nFaces(), -1);
56 
57  const faceList modelFaces(shape.faces());
58 
59  // Loop over all faces of cellShape
60  forAll(modelFaces, cellFaceI)
61  {
62  // face (vertex list)
63  const face& modelFace = modelFaces[cellFaceI];
64 
65  // Loop over all face labels
66  forAll(faces, faceI)
67  {
68  const face& vertLabels = allFaces[faces[faceI]];
69 
70  if (vertLabels == modelFace)
71  {
72  shapeToMesh[cellFaceI] = faces[faceI];
73  break;
74  }
75  }
76 
77  if (shapeToMesh[cellFaceI] == -1)
78  {
79  FatalErrorIn("foamToFieldview : calcFaceAddressing")
80  << "calcFaceAddressing : can't match face to shape.\n"
81  << " shape face:" << modelFace << endl
82  << " face labels:" << faces << endl
83  << " cellI:" << cellI << endl;
84 
85  FatalError << "Faces consist of vertices:" << endl;
86  forAll(faces, faceI)
87  {
89  << " face:" << faces[faceI]
90  << allFaces[faces[faceI]] << endl;
91  }
93  }
94  }
95  return shapeToMesh;
96 }
97 
98 
99 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
100 
101 // Construct from components
102 Foam::fieldviewTopology::fieldviewTopology
103 (
104  const polyMesh& mesh,
105  const bool setWallInfo
106 )
107 :
108  hexLabels_((1+8)*mesh.nCells()),
109  prismLabels_((1+6)*mesh.nCells()),
110  pyrLabels_((1+5)*mesh.nCells()),
111  tetLabels_((1+4)*mesh.nCells()),
112  nPoly_(0),
113  quadFaceLabels_(mesh.boundaryMesh().size()),
114  nPolyFaces_(mesh.boundaryMesh().size())
115 {
116  // Mark all faces that are to be seen as wall for particle
117  // tracking and all cells that use one or more of these walls
118 
119  labelList wallFace(mesh.nFaces(), NOT_A_WALL);
120  boolList wallCell(mesh.nCells(), false);
121 
122  if (setWallInfo)
123  {
124  forAll (mesh.boundaryMesh(), patchI)
125  {
126  const polyPatch& currPatch = mesh.boundaryMesh()[patchI];
127  if
128  (
129  isA<wallPolyPatch>(currPatch)
130  || isA<symmetryPolyPatch>(currPatch)
131  )
132  {
133  forAll(currPatch, patchFaceI)
134  {
135  label meshFaceI = currPatch.start() + patchFaceI;
136 
137  wallFace[meshFaceI] = A_WALL;
138  wallCell[mesh.faceOwner()[meshFaceI]] = true;
139  }
140  }
141  }
142  }
143 
144 
145 
146  const cellModel& tet = *(cellModeller::lookup("tet"));
147  const cellModel& pyr = *(cellModeller::lookup("pyr"));
148  const cellModel& prism = *(cellModeller::lookup("prism"));
149  const cellModel& wedge = *(cellModeller::lookup("wedge"));
150  const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
151  const cellModel& hex = *(cellModeller::lookup("hex"));
152 
153  // Pre calculate headers for cells not on walls
154  labelList notWallFlags(6, NOT_A_WALL);
155  label tetNotWall = fv_encode_elem_header
156  (
157  FV_TET_ELEM_ID, notWallFlags.begin()
158  );
159  label pyrNotWall = fv_encode_elem_header
160  (
161  FV_PYRA_ELEM_ID, notWallFlags.begin()
162  );
163  label prismNotWall = fv_encode_elem_header
164  (
165  FV_PRISM_ELEM_ID, notWallFlags.begin()
166  );
167  label hexNotWall = fv_encode_elem_header
168  (
169  FV_HEX_ELEM_ID, notWallFlags.begin()
170  );
171 
172  // Some aliases
173  const cellList& cellFaces = mesh.cells();
174  const cellShapeList& cellShapes = mesh.cellShapes();
175 
176 
177  label hexi = 0;
178  label prismi = 0;
179  label pyri = 0;
180  label teti = 0;
181 
182  const faceList& allFaces = mesh.faces();
183 
184  labelList wallFlags(6);
185  forAll(cellShapes, celli)
186  {
187  const cellShape& cellShape = cellShapes[celli];
188  const cellModel& cellModel = cellShape.model();
189 
190  if (cellModel == tet)
191  {
192  if (!wallCell[celli])
193  {
194  tetLabels_[teti++] = tetNotWall;
195  }
196  else
197  {
198  labelList modelToMesh = calcFaceAddressing
199  (
200  allFaces, cellShape, cellFaces[celli], celli
201  );
202 
203  wallFlags[0] = wallFace[modelToMesh[0]];
204  wallFlags[1] = wallFace[modelToMesh[1]];
205  wallFlags[2] = wallFace[modelToMesh[2]];
206  wallFlags[3] = wallFace[modelToMesh[3]];
207 
208  tetLabels_[teti++] = fv_encode_elem_header
209  (
210  FV_TET_ELEM_ID, wallFlags.begin()
211  );
212  }
213 
214  tetLabels_[teti++] = cellShape[0] + 1;
215  tetLabels_[teti++] = cellShape[1] + 1;
216  tetLabels_[teti++] = cellShape[2] + 1;
217  tetLabels_[teti++] = cellShape[3] + 1;
218  }
219  else if (cellModel == pyr)
220  {
221  if (!wallCell[celli])
222  {
223  pyrLabels_[pyri++] = pyrNotWall;
224  }
225  else
226  {
227  labelList modelToMesh = calcFaceAddressing
228  (
229  allFaces, cellShape, cellFaces[celli], celli
230  );
231 
232  wallFlags[0] = wallFace[modelToMesh[0]];
233  wallFlags[1] = wallFace[modelToMesh[3]];
234  wallFlags[2] = wallFace[modelToMesh[2]];
235  wallFlags[3] = wallFace[modelToMesh[1]];
236  wallFlags[4] = wallFace[modelToMesh[4]];
237 
238  pyrLabels_[pyri++] = fv_encode_elem_header
239  (
240  FV_PYRA_ELEM_ID, wallFlags.begin()
241  );
242  }
243 
244  pyrLabels_[pyri++] = cellShape[0] + 1;
245  pyrLabels_[pyri++] = cellShape[1] + 1;
246  pyrLabels_[pyri++] = cellShape[2] + 1;
247  pyrLabels_[pyri++] = cellShape[3] + 1;
248  pyrLabels_[pyri++] = cellShape[4] + 1;
249  }
250  else if (cellModel == prism)
251  {
252  if (!wallCell[celli])
253  {
254  prismLabels_[prismi++] = prismNotWall;
255  }
256  else
257  {
258  labelList modelToMesh = calcFaceAddressing
259  (
260  allFaces, cellShape, cellFaces[celli], celli
261  );
262 
263  wallFlags[0] = wallFace[modelToMesh[4]];
264  wallFlags[1] = wallFace[modelToMesh[2]];
265  wallFlags[2] = wallFace[modelToMesh[3]];
266  wallFlags[3] = wallFace[modelToMesh[0]];
267  wallFlags[4] = wallFace[modelToMesh[1]];
268 
269  prismLabels_[prismi++] = fv_encode_elem_header
270  (
271  FV_PRISM_ELEM_ID, wallFlags.begin()
272  );
273  }
274 
275  prismLabels_[prismi++] = cellShape[0] + 1;
276  prismLabels_[prismi++] = cellShape[3] + 1;
277  prismLabels_[prismi++] = cellShape[4] + 1;
278  prismLabels_[prismi++] = cellShape[1] + 1;
279  prismLabels_[prismi++] = cellShape[5] + 1;
280  prismLabels_[prismi++] = cellShape[2] + 1;
281  }
282  else if (cellModel == tetWedge)
283  {
284  // Treat as prism with collapsed edge
285  if (!wallCell[celli])
286  {
287  prismLabels_[prismi++] = prismNotWall;
288  }
289  else
290  {
291  labelList modelToMesh = calcFaceAddressing
292  (
293  allFaces, cellShape, cellFaces[celli], celli
294  );
295 
296  wallFlags[0] = wallFace[modelToMesh[1]];
297  wallFlags[1] = wallFace[modelToMesh[2]];
298  wallFlags[2] = wallFace[modelToMesh[3]];
299  wallFlags[3] = wallFace[modelToMesh[0]];
300  wallFlags[4] = wallFace[modelToMesh[3]];
301 
302  prismLabels_[prismi++] = fv_encode_elem_header
303  (
304  FV_PRISM_ELEM_ID, wallFlags.begin()
305  );
306  }
307 
308  prismLabels_[prismi++] = cellShape[0] + 1;
309  prismLabels_[prismi++] = cellShape[3] + 1;
310  prismLabels_[prismi++] = cellShape[4] + 1;
311  prismLabels_[prismi++] = cellShape[1] + 1;
312  prismLabels_[prismi++] = cellShape[4] + 1;
313  prismLabels_[prismi++] = cellShape[2] + 1;
314  }
315  else if (cellModel == wedge)
316  {
317  if (!wallCell[celli])
318  {
319  hexLabels_[hexi++] = hexNotWall;
320  }
321  else
322  {
323  labelList modelToMesh = calcFaceAddressing
324  (
325  allFaces, cellShape, cellFaces[celli], celli
326  );
327 
328  wallFlags[0] = wallFace[modelToMesh[2]];
329  wallFlags[1] = wallFace[modelToMesh[3]];
330  wallFlags[2] = wallFace[modelToMesh[0]];
331  wallFlags[3] = wallFace[modelToMesh[1]];
332  wallFlags[4] = wallFace[modelToMesh[4]];
333  wallFlags[5] = wallFace[modelToMesh[5]];
334 
335  hexLabels_[hexi++] = fv_encode_elem_header
336  (
337  FV_HEX_ELEM_ID, wallFlags.begin()
338  );
339  }
340  hexLabels_[hexi++] = cellShape[0] + 1;
341  hexLabels_[hexi++] = cellShape[1] + 1;
342  hexLabels_[hexi++] = cellShape[0] + 1;
343  hexLabels_[hexi++] = cellShape[2] + 1;
344  hexLabels_[hexi++] = cellShape[3] + 1;
345  hexLabels_[hexi++] = cellShape[4] + 1;
346  hexLabels_[hexi++] = cellShape[6] + 1;
347  hexLabels_[hexi++] = cellShape[5] + 1;
348  }
349  else if (cellModel == hex)
350  {
351  if (!wallCell[celli])
352  {
353  hexLabels_[hexi++] = hexNotWall;
354  }
355  else
356  {
357  labelList modelToMesh = calcFaceAddressing
358  (
359  allFaces, cellShape, cellFaces[celli], celli
360  );
361 
362  wallFlags[0] = wallFace[modelToMesh[0]];
363  wallFlags[1] = wallFace[modelToMesh[1]];
364  wallFlags[2] = wallFace[modelToMesh[4]];
365  wallFlags[3] = wallFace[modelToMesh[5]];
366  wallFlags[4] = wallFace[modelToMesh[2]];
367  wallFlags[5] = wallFace[modelToMesh[3]];
368 
369  hexLabels_[hexi++] = fv_encode_elem_header
370  (
371  FV_HEX_ELEM_ID, wallFlags.begin()
372  );
373  }
374  hexLabels_[hexi++] = cellShape[0] + 1;
375  hexLabels_[hexi++] = cellShape[1] + 1;
376  hexLabels_[hexi++] = cellShape[3] + 1;
377  hexLabels_[hexi++] = cellShape[2] + 1;
378  hexLabels_[hexi++] = cellShape[4] + 1;
379  hexLabels_[hexi++] = cellShape[5] + 1;
380  hexLabels_[hexi++] = cellShape[7] + 1;
381  hexLabels_[hexi++] = cellShape[6] + 1;
382  }
383  else
384  {
385  nPoly_++;
386  }
387  }
388 
389  hexLabels_.setSize(hexi);
390  prismLabels_.setSize(prismi);
391  pyrLabels_.setSize(pyri);
392  tetLabels_.setSize(teti);
393 
394 
395  //
396  // Patches
397  //
398  forAll(mesh.boundaryMesh(), patchI)
399  {
400  const polyPatch& patchFaces = mesh.boundaryMesh()[patchI];
401 
402  labelList& faceLabels = quadFaceLabels_[patchI];
403 
404  // Faces, each 4 labels. Size big enough
405  faceLabels.setSize(patchFaces.size()*4);
406 
407  label labelI = 0;
408 
409  forAll(patchFaces, faceI)
410  {
411  const face& patchFace = patchFaces[faceI];
412 
413  if (patchFace.size() == 3)
414  {
415  faceLabels[labelI++] = patchFace[0] + 1;
416  faceLabels[labelI++] = patchFace[1] + 1;
417  faceLabels[labelI++] = patchFace[2] + 1;
418  faceLabels[labelI++] = 0; // Fieldview:triangle definition
419  }
420  else if (patchFace.size() == 4)
421  {
422  faceLabels[labelI++] = patchFace[0] + 1;
423  faceLabels[labelI++] = patchFace[1] + 1;
424  faceLabels[labelI++] = patchFace[2] + 1;
425  faceLabels[labelI++] = patchFace[3] + 1;
426  }
427  }
428 
429  faceLabels.setSize(labelI);
430 
431  label nFaces = labelI/4;
432 
433  nPolyFaces_[patchI] = patchFaces.size() - nFaces;
434  }
435 }
436 
437 
438 // ************************ vim: set sw=4 sts=4 et: ************************ //