FreeFOAM The Cross-Platform CFD Toolkit
polyMeshFromShapeMesh.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  Create polyMesh from cell and patch shapes
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "polyMesh.H"
30 #include <OpenFOAM/Time.H>
31 #include <OpenFOAM/primitiveMesh.H>
32 #include <OpenFOAM/DynamicList.H>
33 
34 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
35 
36 Foam::labelListList Foam::polyMesh::cellShapePointCells
37 (
38  const cellShapeList& c
39 ) const
40 {
41  List<DynamicList<label, primitiveMesh::cellsPerPoint_> >
42  pc(points().size());
43 
44  // For each cell
45  forAll(c, i)
46  {
47  // For each vertex
48  const labelList& labels = c[i];
49 
50  forAll(labels, j)
51  {
52  // Set working point label
53  label curPoint = labels[j];
54  DynamicList<label, primitiveMesh::cellsPerPoint_>& curPointCells =
55  pc[curPoint];
56 
57  // Enter the cell label in the point's cell list
58  curPointCells.append(i);
59  }
60  }
61 
62  labelListList pointCellAddr(pc.size());
63 
64  forAll (pc, pointI)
65  {
66  pointCellAddr[pointI].transfer(pc[pointI]);
67  }
68 
69  return pointCellAddr;
70 }
71 
72 
73 Foam::labelList Foam::polyMesh::facePatchFaceCells
74 (
75  const faceList& patchFaces,
76  const labelListList& pointCells,
77  const faceListList& cellsFaceShapes,
78  const label patchID
79 ) const
80 {
81  register bool found;
82 
83  labelList FaceCells(patchFaces.size());
84 
85  forAll(patchFaces, fI)
86  {
87  found = false;
88 
89  const face& curFace = patchFaces[fI];
90  const labelList& facePoints = patchFaces[fI];
91 
92  forAll(facePoints, pointI)
93  {
94  const labelList& facePointCells = pointCells[facePoints[pointI]];
95 
96  forAll(facePointCells, cellI)
97  {
98  faceList cellFaces = cellsFaceShapes[facePointCells[cellI]];
99 
100  forAll(cellFaces, cellFace)
101  {
102  if (cellFaces[cellFace] == curFace)
103  {
104  // Found the cell corresponding to this face
105  FaceCells[fI] = facePointCells[cellI];
106 
107  found = true;
108  }
109  if (found) break;
110  }
111  if (found) break;
112  }
113  if (found) break;
114  }
115 
116  if (!found)
117  {
119  (
120  "polyMesh::facePatchFaceCells(const faceList& patchFaces,"
121  "const labelListList& pointCells,"
122  "const faceListList& cellsFaceShapes,"
123  "const label patchID)"
124  ) << "face " << fI << " in patch " << patchID
125  << " does not have neighbour cell"
126  << " face: " << patchFaces[fI]
127  << abort(FatalError);
128  }
129  }
130 
131  return FaceCells;
132 }
133 
134 
135 Foam::polyMesh::polyMesh
136 (
137  const IOobject& io,
138  const Xfer<pointField>& points,
139  const cellShapeList& cellsAsShapes,
140  const faceListList& boundaryFaces,
141  const wordList& boundaryPatchNames,
142  const wordList& boundaryPatchTypes,
143  const word& defaultBoundaryPatchName,
144  const word& defaultBoundaryPatchType,
145  const wordList& boundaryPatchPhysicalTypes,
146  const bool syncPar
147 )
148 :
149  objectRegistry(io),
150  primitiveMesh(),
151  points_
152  (
153  IOobject
154  (
155  "points",
156  instance(),
157  meshSubDir,
158  *this,
161  ),
162  points
163  ),
164  faces_
165  (
166  IOobject
167  (
168  "faces",
169  instance(),
170  meshSubDir,
171  *this,
174  ),
175  0
176  ),
177  owner_
178  (
179  IOobject
180  (
181  "owner",
182  instance(),
183  meshSubDir,
184  *this,
187  ),
188  0
189  ),
190  neighbour_
191  (
192  IOobject
193  (
194  "neighbour",
195  instance(),
196  meshSubDir,
197  *this,
200  ),
201  0
202  ),
203  clearedPrimitives_(false),
204  boundary_
205  (
206  IOobject
207  (
208  "boundary",
209  instance(),
210  meshSubDir,
211  *this,
214  ),
215  *this,
216  boundaryFaces.size() + 1 // add room for a default patch
217  ),
218  bounds_(points_, syncPar),
219  geometricD_(Vector<label>::zero),
220  solutionD_(Vector<label>::zero),
221  pointZones_
222  (
223  IOobject
224  (
225  "pointZones",
226  instance(),
227  meshSubDir,
228  *this,
231  ),
232  *this,
233  0
234  ),
235  faceZones_
236  (
237  IOobject
238  (
239  "faceZones",
240  instance(),
241  meshSubDir,
242  *this,
245  ),
246  *this,
247  0
248  ),
249  cellZones_
250  (
251  IOobject
252  (
253  "cellZones",
254  instance(),
255  meshSubDir,
256  *this,
259  ),
260  *this,
261  0
262  ),
263  globalMeshDataPtr_(NULL),
264  moving_(false),
265  curMotionTimeIndex_(time().timeIndex()),
266  oldPointsPtr_(NULL)
267 {
268  if (debug)
269  {
270  Info<<"Constructing polyMesh from cell and boundary shapes." << endl;
271  }
272 
273  // Remove all of the old mesh files if they exist
274  removeFiles(instance());
275 
276  // Calculate the faces of all cells
277  // Initialise maximum possible numer of mesh faces to 0
278  label maxFaces = 0;
279 
280  // Set up a list of face shapes for each cell
281  faceListList cellsFaceShapes(cellsAsShapes.size());
282  cellList cells(cellsAsShapes.size());
283 
284  forAll(cellsFaceShapes, cellI)
285  {
286  cellsFaceShapes[cellI] = cellsAsShapes[cellI].faces();
287 
288  cells[cellI].setSize(cellsFaceShapes[cellI].size());
289 
290  // Initialise cells to -1 to flag undefined faces
291  static_cast<labelList&>(cells[cellI]) = -1;
292 
293  // Count maximum possible numer of mesh faces
294  maxFaces += cellsFaceShapes[cellI].size();
295  }
296 
297  // Set size of faces array to maximum possible number of mesh faces
298  faces_.setSize(maxFaces);
299 
300  // Initialise number of faces to 0
301  label nFaces = 0;
302 
303  // set reference to point-cell addressing
304  labelListList PointCells = cellShapePointCells(cellsAsShapes);
305 
306  bool found = false;
307 
308  forAll(cells, cellI)
309  {
310  // Note:
311  // Insertion cannot be done in one go as the faces need to be
312  // added into the list in the increasing order of neighbour
313  // cells. Therefore, all neighbours will be detected first
314  // and then added in the correct order.
315 
316  const faceList& curFaces = cellsFaceShapes[cellI];
317 
318  // Record the neighbour cell
319  labelList neiCells(curFaces.size(), -1);
320 
321  // Record the face of neighbour cell
322  labelList faceOfNeiCell(curFaces.size(), -1);
323 
324  label nNeighbours = 0;
325 
326  // For all faces ...
327  forAll(curFaces, faceI)
328  {
329  // Skip faces that have already been matched
330  if (cells[cellI][faceI] >= 0) continue;
331 
332  found = false;
333 
334  const face& curFace = curFaces[faceI];
335 
336  // Get the list of labels
337  const labelList& curPoints = curFace;
338 
339  // For all points
340  forAll(curPoints, pointI)
341  {
342  // dGget the list of cells sharing this point
343  const labelList& curNeighbours =
344  PointCells[curPoints[pointI]];
345 
346  // For all neighbours
347  forAll(curNeighbours, neiI)
348  {
349  label curNei = curNeighbours[neiI];
350 
351  // Reject neighbours with the lower label
352  if (curNei > cellI)
353  {
354  // Get the list of search faces
355  const faceList& searchFaces = cellsFaceShapes[curNei];
356 
357  forAll(searchFaces, neiFaceI)
358  {
359  if (searchFaces[neiFaceI] == curFace)
360  {
361  // Match!!
362  found = true;
363 
364  // Record the neighbour cell and face
365  neiCells[faceI] = curNei;
366  faceOfNeiCell[faceI] = neiFaceI;
367  nNeighbours++;
368 
369  break;
370  }
371  }
372  if (found) break;
373  }
374  if (found) break;
375  }
376  if (found) break;
377  } // End of current points
378  } // End of current faces
379 
380  // Add the faces in the increasing order of neighbours
381  for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
382  {
383  // Find the lowest neighbour which is still valid
384  label nextNei = -1;
385  label minNei = cells.size();
386 
387  forAll (neiCells, ncI)
388  {
389  if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
390  {
391  nextNei = ncI;
392  minNei = neiCells[ncI];
393  }
394  }
395 
396  if (nextNei > -1)
397  {
398  // Add the face to the list of faces
399  faces_[nFaces] = curFaces[nextNei];
400 
401  // Set cell-face and cell-neighbour-face to current face label
402  cells[cellI][nextNei] = nFaces;
403  cells[neiCells[nextNei]][faceOfNeiCell[nextNei]] = nFaces;
404 
405  // Stop the neighbour from being used again
406  neiCells[nextNei] = -1;
407 
408  // Increment number of faces counter
409  nFaces++;
410  }
411  else
412  {
414  (
415  "polyMesh::polyMesh\n"
416  "(\n"
417  " const IOobject&,\n"
418  " const Xfer<pointField>&,\n"
419  " const cellShapeList& cellsAsShapes,\n"
420  " const faceListList& boundaryFaces,\n"
421  " const wordList& boundaryPatchTypes,\n"
422  " const wordList& boundaryPatchNames,\n"
423  " const word& defaultBoundaryPatchType\n"
424  ")"
425  ) << "Error in internal face insertion"
426  << abort(FatalError);
427  }
428  }
429  }
430 
431  // Do boundary faces
432 
433  labelList patchSizes(boundaryFaces.size(), -1);
434  labelList patchStarts(boundaryFaces.size(), -1);
435 
436  forAll (boundaryFaces, patchI)
437  {
438  const faceList& patchFaces = boundaryFaces[patchI];
439 
440  labelList curPatchFaceCells =
441  facePatchFaceCells
442  (
443  patchFaces,
444  PointCells,
445  cellsFaceShapes,
446  patchI
447  );
448 
449  // Grab the start label
450  label curPatchStart = nFaces;
451 
452  forAll (patchFaces, faceI)
453  {
454  const face& curFace = patchFaces[faceI];
455 
456  const label cellInside = curPatchFaceCells[faceI];
457 
458  faces_[nFaces] = curFace;
459 
460  // get faces of the cell inside
461  const faceList& facesOfCellInside = cellsFaceShapes[cellInside];
462 
463  bool found = false;
464 
465  forAll (facesOfCellInside, cellFaceI)
466  {
467  if (facesOfCellInside[cellFaceI] == curFace)
468  {
469  if (cells[cellInside][cellFaceI] >= 0)
470  {
472  (
473  "polyMesh::polyMesh\n"
474  "(\n"
475  " const IOobject&,\n"
476  " const Xfer<pointField>&,\n"
477  " const cellShapeList& cellsAsShapes,\n"
478  " const faceListList& boundaryFaces,\n"
479  " const wordList& boundaryPatchTypes,\n"
480  " const wordList& boundaryPatchNames,\n"
481  " const word& defaultBoundaryPatchType\n"
482  ")"
483  ) << "Trying to specify a boundary face " << curFace
484  << " on the face on cell " << cellInside
485  << " which is either an internal face or already "
486  << "belongs to some other patch. This is face "
487  << faceI << " of patch "
488  << patchI << " named "
489  << boundaryPatchNames[patchI] << "."
490  << abort(FatalError);
491  }
492 
493  found = true;
494 
495  cells[cellInside][cellFaceI] = nFaces;
496 
497  break;
498  }
499  }
500 
501  if (!found)
502  {
503  FatalErrorIn("polyMesh::polyMesh(... construct from shapes...)")
504  << "face " << faceI << " of patch " << patchI
505  << " does not seem to belong to cell " << cellInside
506  << " which, according to the addressing, "
507  << "should be next to it."
508  << abort(FatalError);
509  }
510 
511  // increment the counter of faces
512  nFaces++;
513  }
514 
515  patchSizes[patchI] = nFaces - curPatchStart;
516  patchStarts[patchI] = curPatchStart;
517  }
518 
519  // Grab "non-existing" faces and put them into a default patch
520 
521  label defaultPatchStart = nFaces;
522 
523  forAll(cells, cellI)
524  {
525  labelList& curCellFaces = cells[cellI];
526 
527  forAll(curCellFaces, faceI)
528  {
529  if (curCellFaces[faceI] == -1) // "non-existent" face
530  {
531  curCellFaces[faceI] = nFaces;
532  faces_[nFaces] = cellsFaceShapes[cellI][faceI];
533 
534  nFaces++;
535  }
536  }
537  }
538 
539  // Reset the size of the face list
540  faces_.setSize(nFaces);
541 
542  // Warning: Patches can only be added once the face list is
543  // completed, as they hold a subList of the face list
544  forAll (boundaryFaces, patchI)
545  {
546  // add the patch to the list
547  boundary_.set
548  (
549  patchI,
551  (
552  boundaryPatchTypes[patchI],
553  boundaryPatchNames[patchI],
554  patchSizes[patchI],
555  patchStarts[patchI],
556  patchI,
557  boundary_
558  )
559  );
560 
561  if
562  (
563  boundaryPatchPhysicalTypes.size()
564  && boundaryPatchPhysicalTypes[patchI].size()
565  )
566  {
567  boundary_[patchI].physicalType() =
568  boundaryPatchPhysicalTypes[patchI];
569  }
570  }
571 
572  label nAllPatches = boundaryFaces.size();
573 
574  if (nFaces > defaultPatchStart)
575  {
576  WarningIn("polyMesh::polyMesh(... construct from shapes...)")
577  << "Found " << nFaces - defaultPatchStart
578  << " undefined faces in mesh; adding to default patch." << endl;
579 
580  boundary_.set
581  (
582  nAllPatches,
584  (
585  defaultBoundaryPatchType,
586  defaultBoundaryPatchName,
587  nFaces - defaultPatchStart,
588  defaultPatchStart,
589  boundary_.size() - 1,
590  boundary_
591  )
592  );
593 
594  nAllPatches++;
595  }
596 
597  // Reset the size of the boundary
598  boundary_.setSize(nAllPatches);
599 
600  // Set the primitive mesh
601  initMesh(cells);
602 
603  if (syncPar)
604  {
605  // Calculate topology for the patches (processor-processor comms etc.)
606  boundary_.updateMesh();
607 
608  // Calculate the geometry for the patches (transformation tensors etc.)
609  boundary_.calcGeometry();
610  }
611 
612  if (debug)
613  {
614  if (checkMesh())
615  {
616  Info << "Mesh OK" << endl;
617  }
618  }
619 }
620 
621 
622 // ************************ vim: set sw=4 sts=4 et: ************************ //