FreeFOAM The Cross-Platform CFD Toolkit
polyMeshZipUpCells.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 "polyMeshZipUpCells.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include <OpenFOAM/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  {
120  WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
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 =
451  Points[unorderedEdge[unorderedEdge.size() - 1]]
452  - startPoint;
453 
454  for (label i = 1; i < dist.size(); i++)
455  {
456  dist[i] = (Points[unorderedEdge[i]] - startPoint) & dir;
457  }
458 
459  // Sort points
460  labelList orderedEdge(unorderedEdge.size(), -1);
461  boolList used(unorderedEdge.size(), false);
462 
463  forAll (orderedEdge, epI)
464  {
465  label nextPoint = -1;
466  scalar minDist = GREAT;
467 
468  forAll (dist, i)
469  {
470  if (!used[i] && dist[i] < minDist)
471  {
472  minDist = dist[i];
473  nextPoint = i;
474  }
475  }
476 
477  // Insert the next point
478  orderedEdge[epI] = unorderedEdge[nextPoint];
479  used[nextPoint] = true;
480  }
481 
482 # ifdef DEBUG_ORDER
483  Info<< "unorderedEdge: " << unorderedEdge << nl
484  << "orderedEdge: " << orderedEdge << endl;
485 # endif
486 
487  // check for duplicate points in the ordered edge
488  forAll (orderedEdge, checkI)
489  {
490  for
491  (
492  label checkJ = checkI + 1;
493  checkJ < orderedEdge.size();
494  checkJ++
495  )
496  {
497  if (orderedEdge[checkI] == orderedEdge[checkJ])
498  {
499  WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
500  << "Duplicate point found in edge to insert. "
501  << nl << "Point: " << orderedEdge[checkI]
502  << " edge: " << orderedEdge << endl;
503 
504  problemCells.insert(cellI);
505  }
506  }
507  }
508 
509  edge testEdge
510  (
511  orderedEdge[0],
512  orderedEdge[orderedEdge.size() - 1]
513  );
514 
515  // In order to avoid edge-to-edge comparison, get faces using
516  // point-face addressing in two goes.
517  const labelList& startPF = pFaces[testEdge.start()];
518  const labelList& endPF = pFaces[testEdge.start()];
519 
520  labelList facesSharingEdge(startPF.size() + endPF.size());
521  label nfse = 0;
522 
523  forAll (startPF, pfI)
524  {
525  facesSharingEdge[nfse++] = startPF[pfI];
526  }
527  forAll (endPF, pfI)
528  {
529  facesSharingEdge[nfse++] = endPF[pfI];
530  }
531 
532  forAll (facesSharingEdge, faceI)
533  {
534  bool faceChanges = false;
535 
536  // Label of the face being analysed
537  const label currentFaceIndex = facesSharingEdge[faceI];
538 
539  const edgeList curFaceEdges =
540  oldFaces[currentFaceIndex].edges();
541 
542  forAll (curFaceEdges, cfeI)
543  {
544  if (curFaceEdges[cfeI] == testEdge)
545  {
546  faceChanges = true;
547  break;
548  }
549  }
550 
551  if (faceChanges)
552  {
553  nChangedFacesInMesh++;
554  // In order to avoid loosing point from multiple
555  // insertions into the same face, the new face
556  // will be change incrementally.
557  // 1) Check if all the internal points of the edge
558  // to add already exist in the face. If so, the
559  // edge has already been included 2) Check if the
560  // point insertion occurs on an edge which is
561  // still untouched. If so, simply insert
562  // additional points into the face. 3) If not,
563  // the edge insertion occurs on an already
564  // modified edge. ???
565 
566  face& newFace = newFaces[currentFaceIndex];
567 
568  bool allPointsPresent = true;
569 
570  forAll (orderedEdge, oeI)
571  {
572  bool curPointFound = false;
573 
574  forAll (newFace, nfI)
575  {
576  if (newFace[nfI] == orderedEdge[oeI])
577  {
578  curPointFound = true;
579  break;
580  }
581  }
582 
583  allPointsPresent =
584  allPointsPresent && curPointFound;
585  }
586 
587 # ifdef DEBUG_ZIPUP
588  if (allPointsPresent)
589  {
590  Info << "All points present" << endl;
591  }
592 # endif
593 
594  if (!allPointsPresent)
595  {
596  // Not all points are already present. The
597  // new edge will need to be inserted into the
598  // face.
599 
600  // Check to see if a new edge fits onto an
601  // untouched edge of the face. Make sure the
602  // edges are grabbed before the face is
603  // resized.
604  edgeList newFaceEdges = newFace.edges();
605 
606 # ifdef DEBUG_ZIPUP
607  Info << "Not all points present." << endl;
608 # endif
609 
610  label nNewFacePoints = 0;
611 
612  bool edgeAdded = false;
613 
614  forAll (newFaceEdges, curFacEdgI)
615  {
616  // Does the current edge change?
617  if (newFaceEdges[curFacEdgI] == testEdge)
618  {
619  // Found an edge match
620  edgeAdded = true;
621 
622  // Resize the face to accept additional
623  // points
624  newFace.setSize
625  (
626  newFace.size()
627  + orderedEdge.size() - 2
628  );
629 
630  if
631  (
632  newFaceEdges[curFacEdgI].start()
633  == testEdge.start()
634  )
635  {
636  // insertion in ascending order
637  for
638  (
639  label i = 0;
640  i < orderedEdge.size() - 1;
641  i++
642  )
643  {
644  newFace[nNewFacePoints] =
645  orderedEdge[i];
646  nNewFacePoints++;
647  }
648  }
649  else
650  {
651  // insertion in reverse order
652  for
653  (
654  label i = orderedEdge.size() - 1;
655  i > 0;
656  i--
657  )
658  {
659  newFace[nNewFacePoints] =
660  orderedEdge[i];
661  nNewFacePoints++;
662  }
663  }
664  }
665  else
666  {
667  // Does not fit onto this edge.
668  // Copy the next point into the face
669  newFace[nNewFacePoints] =
670  newFaceEdges[curFacEdgI].start();
671  nNewFacePoints++;
672  }
673  }
674 
675 # ifdef DEBUG_ZIPUP
676  Info<< "oldFace: "
677  << oldFaces[currentFaceIndex] << nl
678  << "newFace: " << newFace << endl;
679 # endif
680 
681  // Check for duplicate points in the new face
682  forAll (newFace, checkI)
683  {
684  for
685  (
686  label checkJ = checkI + 1;
687  checkJ < newFace.size();
688  checkJ++
689  )
690  {
691  if (newFace[checkI] == newFace[checkJ])
692  {
693  WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
694  << "Duplicate point found "
695  << "in the new face. " << nl
696  << "Point: "
697  << orderedEdge[checkI]
698  << " face: "
699  << newFace << endl;
700 
701  problemCells.insert(cellI);
702  }
703  }
704  }
705 
706  // Check if the edge is added.
707  // If not, then it comes on top of an already
708  // modified edge and they need to be
709  // merged in together.
710  if (!edgeAdded)
711  {
712  Info<< "This edge modifies an already modified "
713  << "edge. Point insertions skipped."
714  << endl;
715  }
716  }
717  }
718  }
719  }
720  }
721 
722  if (problemCells.size())
723  {
724  // This cycle has failed. Print out the problem cells
725  labelList toc(problemCells.toc());
726  sort(toc);
727 
728  FatalErrorIn("void polyMeshZipUpCells(polyMesh& mesh)")
729  << "Found " << problemCells.size() << " problem cells." << nl
730  << "Cells: " << toc
731  << abort(FatalError);
732  }
733 
734  Info<< "Cycle " << ++nCycles
735  << " changed " << nChangedFacesInMesh << " faces." << endl;
736 
737 
738  const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
739 
740  // Reset the polyMesh. Number of points/faces/cells/patches stays the
741  // same, only the faces themselves have changed so clear all derived
742  // (edge, point) addressing.
743 
744  // Collect the patch sizes
745  labelList patchSizes(bMesh.size(), 0);
746  labelList patchStarts(bMesh.size(), 0);
747 
748  forAll (bMesh, patchI)
749  {
750  patchSizes[patchI] = bMesh[patchI].size();
751  patchStarts[patchI] = bMesh[patchI].start();
752  }
753 
754  // Reset the mesh. Number of active faces is one beyond the last patch
755  // (patches guaranteed to be in increasing order)
756  mesh.resetPrimitives
757  (
759  xferMove(newFaces),
762  patchSizes,
763  patchStarts,
764  true // boundary forms valid boundary mesh.
765  );
766 
767  // Reset any addressing on face zones.
768  mesh.faceZones().clearAddressing();
769 
770  // Clear the addressing
771  mesh.clearOut();
772 
773  } while (nChangedFacesInMesh > 0 || nCycles > 100);
774 
775  // Flags the mesh files as being changed
776  mesh.setInstance(mesh.time().timeName());
777 
778  if (nChangedFacesInMesh > 0)
779  {
780  FatalErrorIn("void polyMeshZipUpCells(polyMesh& mesh)")
781  << "cell zip-up failed after 100 cycles. Probable problem "
782  << "with the original mesh"
783  << abort(FatalError);
784  }
785 
786  return nCycles != 1;
787 }
788 
789 
790 // ************************ vim: set sw=4 sts=4 et: ************************ //