polyMeshZipUpCells.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 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 "polyMeshZipUpCells.H"
27 #include "polyMesh.H"
28 #include "Time.H"
29 
30 // #define DEBUG_ZIPUP 1
31 // #define DEBUG_CHAIN 1
32 // #define DEBUG_ORDER 1
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
37 {
38  if (polyMesh::debug)
39  {
40  Info<< "bool polyMeshZipUpCells(polyMesh& mesh) const: "
41  << "zipping up topologically open cells" << endl;
42  }
43 
44  // Algorithm:
45  // Take the original mesh and visit all cells. For every cell
46  // calculate the edges of all faces on the cells. A cell is
47  // correctly topologically closed when all the edges are referenced
48  // by exactly two faces. If the edges are referenced only by a
49  // single face, additional vertices need to be inserted into some
50  // of the faces (topological closedness). If an edge is
51  // referenced by more that two faces, there is an error in
52  // topological closedness.
53  // Point insertion into the faces is done by attempting to create
54  // closed loops and inserting the intermediate points into the
55  // defining edge
56  // Note:
57  // The algorithm is recursive and changes the mesh faces in each
58  // pass. It is therefore essential to discard the addressing
59  // after every pass. The algorithm is completed when the mesh
60  // stops changing.
61 
62  label nChangedFacesInMesh = 0;
63  label nCycles = 0;
64 
65  labelHashSet problemCells;
66 
67  do
68  {
69  nChangedFacesInMesh = 0;
70 
71  const cellList& Cells = mesh.cells();
72  const pointField& Points = mesh.points();
73 
74  faceList newFaces = mesh.faces();
75 
76  const faceList& oldFaces = mesh.faces();
77  const labelListList& pFaces = mesh.pointFaces();
78 
79  forAll(Cells, celli)
80  {
81  const labelList& curFaces = Cells[celli];
82  const edgeList cellEdges = Cells[celli].edges(oldFaces);
83  const labelList cellPoints = Cells[celli].labels(oldFaces);
84 
85  // Find the edges used only once in the cell
86 
87  labelList edgeUsage(cellEdges.size(), 0);
88 
89  forAll(curFaces, facei)
90  {
91  edgeList curFaceEdges = oldFaces[curFaces[facei]].edges();
92 
93  forAll(curFaceEdges, faceEdgeI)
94  {
95  const edge& curEdge = curFaceEdges[faceEdgeI];
96 
97  forAll(cellEdges, cellEdgeI)
98  {
99  if (cellEdges[cellEdgeI] == curEdge)
100  {
101  edgeUsage[cellEdgeI]++;
102  break;
103  }
104  }
105  }
106  }
107 
108  edgeList singleEdges(cellEdges.size());
109  label nSingleEdges = 0;
110 
111  forAll(edgeUsage, edgeI)
112  {
113  if (edgeUsage[edgeI] == 1)
114  {
115  singleEdges[nSingleEdges] = cellEdges[edgeI];
116  nSingleEdges++;
117  }
118  else if (edgeUsage[edgeI] != 2)
119  {
121  << "edge " << cellEdges[edgeI] << " in cell " << celli
122  << " used " << edgeUsage[edgeI] << " times. " << nl
123  << "Should be 1 or 2 - serious error "
124  << "in mesh structure. " << endl;
125 
126  #ifdef DEBUG_ZIPUP
127  forAll(curFaces, facei)
128  {
129  Info<< "face: " << oldFaces[curFaces[facei]]
130  << endl;
131  }
132 
133  Info<< "Cell edges: " << cellEdges << nl
134  << "Edge usage: " << edgeUsage << nl
135  << "Cell points: " << cellPoints << endl;
136 
137  forAll(cellPoints, cpI)
138  {
139  Info<< "vertex create \"" << cellPoints[cpI]
140  << "\" coordinates "
141  << Points[cellPoints[cpI]] << endl;
142  }
143  #endif
144 
145  // Gather the problem cell
146  problemCells.insert(celli);
147  }
148  }
149 
150  // Check if the cell is already zipped up
151  if (nSingleEdges == 0) continue;
152 
153  singleEdges.setSize(nSingleEdges);
154 
155  #ifdef DEBUG_ZIPUP
156  Info<< "Cell " << celli << endl;
157 
158  forAll(curFaces, facei)
159  {
160  Info<< "face: " << oldFaces[curFaces[facei]] << endl;
161  }
162 
163  Info<< "Cell edges: " << cellEdges << nl
164  << "Edge usage: " << edgeUsage << nl
165  << "Single edges: " << singleEdges << nl
166  << "Cell points: " << cellPoints << endl;
167 
168  forAll(cellPoints, cpI)
169  {
170  Info<< "vertex create \"" << cellPoints[cpI]
171  << "\" coordinates "
172  << points()[cellPoints[cpI]] << endl;
173  }
174  #endif
175 
176  // Loop through all single edges and mark the points they use
177  // points marked twice are internal to edge; those marked more than
178  // twice are corners
179 
180  labelList pointUsage(cellPoints.size(), 0);
181 
182  forAll(singleEdges, edgeI)
183  {
184  const edge& curEdge = singleEdges[edgeI];
185 
186  forAll(cellPoints, pointi)
187  {
188  if
189  (
190  cellPoints[pointi] == curEdge.start()
191  || cellPoints[pointi] == curEdge.end()
192  )
193  {
194  pointUsage[pointi]++;
195  }
196  }
197  }
198 
199  boolList singleEdgeUsage(singleEdges.size(), false);
200 
201  // loop through all edges and eliminate the ones that are
202  // blocked out
203  forAll(singleEdges, edgeI)
204  {
205  bool blockedHead = false;
206  bool blockedTail = false;
207 
208  label newEdgeStart = singleEdges[edgeI].start();
209  label newEdgeEnd = singleEdges[edgeI].end();
210 
211  // check that the edge has not got all ends blocked
212  forAll(cellPoints, pointi)
213  {
214  if (cellPoints[pointi] == newEdgeStart)
215  {
216  if (pointUsage[pointi] > 2)
217  {
218  blockedHead = true;
219  }
220  }
221  else if (cellPoints[pointi] == newEdgeEnd)
222  {
223  if (pointUsage[pointi] > 2)
224  {
225  blockedTail = true;
226  }
227  }
228  }
229 
230  if (blockedHead && blockedTail)
231  {
232  // Eliminating edge singleEdges[edgeI] as blocked
233  singleEdgeUsage[edgeI] = true;
234  }
235  }
236 
237  // Go through the points and start from the point used twice
238  // check all the edges to find the edges starting from this point
239  // add the
240 
241  labelListList edgesToInsert(singleEdges.size());
242  label nEdgesToInsert = 0;
243 
244  // Find a good edge
245  forAll(singleEdges, edgeI)
246  {
247  SLList<label> pointChain;
248 
249  bool blockHead = false;
250  bool blockTail = false;
251 
252  if (!singleEdgeUsage[edgeI])
253  {
254  // found a new edge
255  singleEdgeUsage[edgeI] = true;
256 
257  label newEdgeStart = singleEdges[edgeI].start();
258  label newEdgeEnd = singleEdges[edgeI].end();
259 
260  pointChain.insert(newEdgeStart);
261  pointChain.append(newEdgeEnd);
262 
263  #ifdef DEBUG_CHAIN
264  Info<< "found edge to start with: "
265  << singleEdges[edgeI] << endl;
266  #endif
267 
268  // Check if head or tail are blocked
269  forAll(cellPoints, pointi)
270  {
271  if (cellPoints[pointi] == newEdgeStart)
272  {
273  if (pointUsage[pointi] > 2)
274  {
275  #ifdef DEBUG_CHAIN
276  Info<< "start head blocked" << endl;
277  #endif
278 
279  blockHead = true;
280  }
281  }
282  else if (cellPoints[pointi] == newEdgeEnd)
283  {
284  if (pointUsage[pointi] > 2)
285  {
286  #ifdef DEBUG_CHAIN
287  Info<< "start tail blocked" << endl;
288  #endif
289 
290  blockTail = true;
291  }
292  }
293  }
294 
295  bool stopSearching = false;
296 
297  // Go through the unused edges and try to chain them up
298  do
299  {
300  stopSearching = false;
301 
302  forAll(singleEdges, addEdgeI)
303  {
304  if (!singleEdgeUsage[addEdgeI])
305  {
306  // Grab start and end of the candidate
307  label addStart =
308  singleEdges[addEdgeI].start();
309 
310  label addEnd =
311  singleEdges[addEdgeI].end();
312 
313  #ifdef DEBUG_CHAIN
314  Info<< "Trying candidate "
315  << singleEdges[addEdgeI] << endl;
316  #endif
317 
318  // Try to add the edge onto the head
319  if (!blockHead)
320  {
321  if (pointChain.first() == addStart)
322  {
323  // Added at start mark as used
324  pointChain.insert(addEnd);
325 
326  singleEdgeUsage[addEdgeI] = true;
327  }
328  else if (pointChain.first() == addEnd)
329  {
330  pointChain.insert(addStart);
331 
332  singleEdgeUsage[addEdgeI] = true;
333  }
334  }
335 
336  // Try the other end only if the first end
337  // did not add it
338  if (!blockTail && !singleEdgeUsage[addEdgeI])
339  {
340  if (pointChain.last() == addStart)
341  {
342  // Added at start mark as used
343  pointChain.append(addEnd);
344 
345  singleEdgeUsage[addEdgeI] = true;
346  }
347  else if (pointChain.last() == addEnd)
348  {
349  pointChain.append(addStart);
350 
351  singleEdgeUsage[addEdgeI] = true;
352  }
353  }
354 
355  // check if the new head or tail are blocked
356  label curEdgeStart = pointChain.first();
357  label curEdgeEnd = pointChain.last();
358 
359  #ifdef DEBUG_CHAIN
360  Info<< "curEdgeStart: " << curEdgeStart
361  << " curEdgeEnd: " << curEdgeEnd << endl;
362  #endif
363 
364  forAll(cellPoints, pointi)
365  {
366  if (cellPoints[pointi] == curEdgeStart)
367  {
368  if (pointUsage[pointi] > 2)
369  {
370  #ifdef DEBUG_CHAIN
371  Info<< "head blocked" << endl;
372  #endif
373 
374  blockHead = true;
375  }
376  }
377  else if (cellPoints[pointi] == curEdgeEnd)
378  {
379  if (pointUsage[pointi] > 2)
380  {
381  #ifdef DEBUG_CHAIN
382  Info<< "tail blocked" << endl;
383  #endif
384 
385  blockTail = true;
386  }
387  }
388  }
389 
390  // Check if the loop is closed
391  if (curEdgeStart == curEdgeEnd)
392  {
393  #ifdef DEBUG_CHAIN
394  Info<< "closed loop" << endl;
395  #endif
396 
397  pointChain.removeHead();
398 
399  blockHead = true;
400  blockTail = true;
401 
402  stopSearching = true;
403  }
404 
405  #ifdef DEBUG_CHAIN
406  Info<< "current pointChain: " << pointChain
407  << endl;
408  #endif
409 
410  if (stopSearching) break;
411  }
412  }
413  } while (stopSearching);
414  }
415 
416  #ifdef DEBUG_CHAIN
417  Info<< "completed patch chain: " << pointChain << endl;
418  #endif
419 
420  if (pointChain.size() > 2)
421  {
422  edgesToInsert[nEdgesToInsert] = pointChain;
423  nEdgesToInsert++;
424  }
425  }
426 
427  edgesToInsert.setSize(nEdgesToInsert);
428 
429  #ifdef DEBUG_ZIPUP
430  Info<< "edgesToInsert: " << edgesToInsert << endl;
431  #endif
432 
433  // Insert the edges into a list of faces
434  forAll(edgesToInsert, edgeToInsertI)
435  {
436  // Order the points of the edge
437  // Warning: the ordering must be parametric, because in
438  // the case of multiple point insertion onto the same edge
439  // it is possible to get non-cyclic loops
440  //
441 
442  const labelList& unorderedEdge = edgesToInsert[edgeToInsertI];
443 
444  scalarField dist(unorderedEdge.size());
445 
446  // Calculate distance
447  point startPoint = Points[unorderedEdge[0]];
448  dist[0] = 0;
449 
450  vector dir = Points[unorderedEdge.last()] - startPoint;
451 
452  for (label i = 1; i < dist.size(); i++)
453  {
454  dist[i] = (Points[unorderedEdge[i]] - startPoint) & dir;
455  }
456 
457  // Sort points
458  labelList orderedEdge(unorderedEdge.size(), -1);
459  boolList used(unorderedEdge.size(), false);
460 
461  forAll(orderedEdge, epI)
462  {
463  label nextPoint = -1;
464  scalar minDist = great;
465 
466  forAll(dist, i)
467  {
468  if (!used[i] && dist[i] < minDist)
469  {
470  minDist = dist[i];
471  nextPoint = i;
472  }
473  }
474 
475  // Insert the next point
476  orderedEdge[epI] = unorderedEdge[nextPoint];
477  used[nextPoint] = true;
478  }
479 
480  #ifdef DEBUG_ORDER
481  Info<< "unorderedEdge: " << unorderedEdge << nl
482  << "orderedEdge: " << orderedEdge << endl;
483  #endif
484 
485  // check for duplicate points in the ordered edge
486  forAll(orderedEdge, checkI)
487  {
488  for
489  (
490  label checkJ = checkI + 1;
491  checkJ < orderedEdge.size();
492  checkJ++
493  )
494  {
495  if (orderedEdge[checkI] == orderedEdge[checkJ])
496  {
498  << "Duplicate point found in edge to insert. "
499  << nl << "Point: " << orderedEdge[checkI]
500  << " edge: " << orderedEdge << endl;
501 
502  problemCells.insert(celli);
503  }
504  }
505  }
506 
507  edge testEdge
508  (
509  orderedEdge[0],
510  orderedEdge.last()
511  );
512 
513  // In order to avoid edge-to-edge comparison, get faces using
514  // point-face addressing in two goes.
515  const labelList& startPF = pFaces[testEdge.start()];
516  const labelList& endPF = pFaces[testEdge.start()];
517 
518  labelList facesSharingEdge(startPF.size() + endPF.size());
519  label nfse = 0;
520 
521  forAll(startPF, pfI)
522  {
523  facesSharingEdge[nfse++] = startPF[pfI];
524  }
525  forAll(endPF, pfI)
526  {
527  facesSharingEdge[nfse++] = endPF[pfI];
528  }
529 
530  forAll(facesSharingEdge, facei)
531  {
532  bool faceChanges = false;
533 
534  // Label of the face being analysed
535  const label currentFaceIndex = facesSharingEdge[facei];
536 
537  const edgeList curFaceEdges =
538  oldFaces[currentFaceIndex].edges();
539 
540  forAll(curFaceEdges, cfeI)
541  {
542  if (curFaceEdges[cfeI] == testEdge)
543  {
544  faceChanges = true;
545  break;
546  }
547  }
548 
549  if (faceChanges)
550  {
551  nChangedFacesInMesh++;
552  // In order to avoid losing point from multiple
553  // insertions into the same face, the new face
554  // will be change incrementally.
555  // 1) Check if all the internal points of the edge
556  // to add already exist in the face. If so, the
557  // edge has already been included 2) Check if the
558  // point insertion occurs on an edge which is
559  // still untouched. If so, simply insert
560  // additional points into the face. 3) If not,
561  // the edge insertion occurs on an already
562  // modified edge. ???
563 
564  face& newFace = newFaces[currentFaceIndex];
565 
566  bool allPointsPresent = true;
567 
568  forAll(orderedEdge, oeI)
569  {
570  bool curPointFound = false;
571 
572  forAll(newFace, nfI)
573  {
574  if (newFace[nfI] == orderedEdge[oeI])
575  {
576  curPointFound = true;
577  break;
578  }
579  }
580 
581  allPointsPresent =
582  allPointsPresent && curPointFound;
583  }
584 
585  #ifdef DEBUG_ZIPUP
586  if (allPointsPresent)
587  {
588  Info<< "All points present" << endl;
589  }
590  #endif
591 
592  if (!allPointsPresent)
593  {
594  // Not all points are already present. The
595  // new edge will need to be inserted into the
596  // face.
597 
598  // Check to see if a new edge fits onto an
599  // untouched edge of the face. Make sure the
600  // edges are grabbed before the face is
601  // resized.
602  edgeList newFaceEdges = newFace.edges();
603 
604  #ifdef DEBUG_ZIPUP
605  Info<< "Not all points present." << endl;
606  #endif
607 
608  label nNewFacePoints = 0;
609 
610  bool edgeAdded = false;
611 
612  forAll(newFaceEdges, curFacEdgI)
613  {
614  // Does the current edge change?
615  if (newFaceEdges[curFacEdgI] == testEdge)
616  {
617  // Found an edge match
618  edgeAdded = true;
619 
620  // Resize the face to accept additional
621  // points
622  newFace.setSize
623  (
624  newFace.size()
625  + orderedEdge.size() - 2
626  );
627 
628  if
629  (
630  newFaceEdges[curFacEdgI].start()
631  == testEdge.start()
632  )
633  {
634  // insertion in ascending order
635  for
636  (
637  label i = 0;
638  i < orderedEdge.size() - 1;
639  i++
640  )
641  {
642  newFace[nNewFacePoints] =
643  orderedEdge[i];
644  nNewFacePoints++;
645  }
646  }
647  else
648  {
649  // insertion in reverse order
650  for
651  (
652  label i = orderedEdge.size() - 1;
653  i > 0;
654  i--
655  )
656  {
657  newFace[nNewFacePoints] =
658  orderedEdge[i];
659  nNewFacePoints++;
660  }
661  }
662  }
663  else
664  {
665  // Does not fit onto this edge.
666  // Copy the next point into the face
667  newFace[nNewFacePoints] =
668  newFaceEdges[curFacEdgI].start();
669  nNewFacePoints++;
670  }
671  }
672 
673  #ifdef DEBUG_ZIPUP
674  Info<< "oldFace: "
675  << oldFaces[currentFaceIndex] << nl
676  << "newFace: " << newFace << endl;
677  #endif
678 
679  // Check for duplicate points in the new face
680  forAll(newFace, checkI)
681  {
682  for
683  (
684  label checkJ = checkI + 1;
685  checkJ < newFace.size();
686  checkJ++
687  )
688  {
689  if (newFace[checkI] == newFace[checkJ])
690  {
692  << "Duplicate point found "
693  << "in the new face. " << nl
694  << "Point: "
695  << orderedEdge[checkI]
696  << " face: "
697  << newFace << endl;
698 
699  problemCells.insert(celli);
700  }
701  }
702  }
703 
704  // Check if the edge is added.
705  // If not, then it comes on top of an already
706  // modified edge and they need to be
707  // merged in together.
708  if (!edgeAdded)
709  {
710  Info<< "This edge modifies an already modified "
711  << "edge. Point insertions skipped."
712  << endl;
713  }
714  }
715  }
716  }
717  }
718  }
719 
720  if (problemCells.size())
721  {
722  // This cycle has failed. Print out the problem cells
723  labelList toc(problemCells.toc());
724  sort(toc);
725 
727  << "Found " << problemCells.size() << " problem cells." << nl
728  << "Cells: " << toc
729  << abort(FatalError);
730  }
731 
732  Info<< "Cycle " << ++nCycles
733  << " changed " << nChangedFacesInMesh << " faces." << endl;
734 
735 
736  const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
737 
738  // Reset the polyMesh. Number of points/faces/cells/patches stays the
739  // same, only the faces themselves have changed so clear all derived
740  // (edge, point) addressing.
741 
742  // Collect the patch sizes
743  labelList patchSizes(bMesh.size(), 0);
744  labelList patchStarts(bMesh.size(), 0);
745 
747  {
748  patchSizes[patchi] = bMesh[patchi].size();
749  patchStarts[patchi] = bMesh[patchi].start();
750  }
751 
752  // Reset the mesh. Number of active faces is one beyond the last patch
753  // (patches guaranteed to be in increasing order)
754  mesh.resetPrimitives
755  (
756  NullObjectMove<pointField>(),
757  move(newFaces),
758  NullObjectMove<labelList>(),
759  NullObjectMove<labelList>(),
760  patchSizes,
761  patchStarts,
762  true // boundary forms valid boundary mesh.
763  );
764 
765  // Reset any addressing on face zones.
766  mesh.faceZones().clearAddressing();
767 
768  // Clear the addressing
769  mesh.clearOut();
770 
771  } while (nChangedFacesInMesh > 0 || nCycles > 100);
772 
773  // Flags the mesh files as being changed
774  mesh.setInstance(mesh.time().name());
775 
776  if (nChangedFacesInMesh > 0)
777  {
779  << "with the original mesh"
780  << abort(FatalError);
781  }
782 
783  return nCycles != 1;
784 }
785 
786 
787 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
List< Key > toc() const
Return the table of contents.
Definition: HashTable.C:202
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
Template class for non-intrusive linked lists.
Definition: LList.H:76
void insert(const T &a)
Add at head of list.
Definition: LList.H:166
T & first()
Return the first entry added.
Definition: LList.H:139
T & last()
Return the last entry added.
Definition: LList.H:151
T removeHead()
Remove and return head.
Definition: LList.H:178
void append(const T &a)
Add at tail of list.
Definition: LList.H:172
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void setSize(const label)
Reset size of List.
Definition: List.C:281
void clearAddressing()
Clear addressing.
Definition: MeshZones.C:387
iterator end()
Return an iterator to end traversing the UList.
Definition: UListI.H:224
T & last()
Return the last element of the list.
Definition: UListI.H:128
const word & name() const
Return const reference to name.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:61
label end() const
Return end vertex label.
Definition: edgeI.H:92
label start() const
Return start vertex label.
Definition: edgeI.H:81
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
edgeList edges() const
Return edges in face point ordering,.
Definition: face.C:406
const Time & time() const
Return time.
Foam::polyBoundaryMesh.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:445
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1374
void resetPrimitives(pointField &&points, faceList &&faces, labelList &&owner, labelList &&neighbour, const labelList &patchSizes, const labelList &patchStarts, const bool validBoundary=true)
Reset mesh primitive data. Assumes all patch info correct.
Definition: polyMesh.C:701
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:403
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1361
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:78
void clearOut()
Clear all geometry and addressing unnecessary for CFD.
const labelListList & pointFaces() const
const cellList & cells() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
const polyBoundaryMesh & bMesh
label patchi
const pointField & points
#define WarningInFunction
Report a warning using Foam::Warning.
scalar minDist(const List< pointIndexHit > &hitList)
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
errorManip< error > abort(error &err)
Definition: errorManip.H:131
messageStream Info
void sort(UList< T > &)
Definition: UList.C:115
bool polyMeshZipUpCells(polyMesh &mesh)
error FatalError
static const char nl
Definition: Ostream.H:260
Cell zip-up tool. This function modifies the list of faces such that all the cells are topologically ...
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< small) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &mergedCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:229