ENH: autoLayerDriver: better layering information message
[OpenFOAM-2.0.x.git] / src / meshTools / polyMeshZipUpCells / polyMeshZipUpCells.C
bloba1256b0d6f66e10873fe81d6fbd2ea495219f751
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
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.
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.
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/>.
24 \*---------------------------------------------------------------------------*/
26 #include "polyMeshZipUpCells.H"
27 #include "polyMesh.H"
28 #include "Time.H"
30 // #define DEBUG_ZIPUP 1
31 // #define DEBUG_CHAIN 1
32 // #define DEBUG_ORDER 1
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 bool Foam::polyMeshZipUpCells(polyMesh& mesh)
38     if (polyMesh::debug)
39     {
40         Info<< "bool polyMeshZipUpCells(polyMesh& mesh) const: "
41             << "zipping up topologically open cells" << endl;
42     }
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.
62     label nChangedFacesInMesh = 0;
63     label nCycles = 0;
65     labelHashSet problemCells;
67     do
68     {
69         nChangedFacesInMesh = 0;
71         const cellList& Cells = mesh.cells();
72         const pointField& Points = mesh.points();
74         faceList newFaces = mesh.faces();
76         const faceList& oldFaces = mesh.faces();
77         const labelListList& pFaces = mesh.pointFaces();
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);
85             // Find the edges used only once in the cell
87             labelList edgeUsage(cellEdges.size(), 0);
89             forAll(curFaces, faceI)
90             {
91                 edgeList curFaceEdges = oldFaces[curFaces[faceI]].edges();
93                 forAll(curFaceEdges, faceEdgeI)
94                 {
95                     const edge& curEdge = curFaceEdges[faceEdgeI];
97                     forAll(cellEdges, cellEdgeI)
98                     {
99                         if (cellEdges[cellEdgeI] == curEdge)
100                         {
101                             edgeUsage[cellEdgeI]++;
102                             break;
103                         }
104                     }
105                 }
106             }
108             edgeList singleEdges(cellEdges.size());
109             label nSingleEdges = 0;
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;
126 #                   ifdef DEBUG_ZIPUP
127                     forAll(curFaces, faceI)
128                     {
129                         Info<< "face: " << oldFaces[curFaces[faceI]]
130                             << endl;
131                     }
133                     Info<< "Cell edges: " << cellEdges << nl
134                         << "Edge usage: " << edgeUsage << nl
135                         << "Cell points: " << cellPoints << endl;
137                     forAll(cellPoints, cpI)
138                     {
139                         Info<< "vertex create \"" << cellPoints[cpI]
140                             << "\" coordinates "
141                             << Points[cellPoints[cpI]] << endl;
142                     }
143 #                   endif
145                     // Gather the problem cell
146                     problemCells.insert(cellI);
147                 }
148             }
150             // Check if the cell is already zipped up
151             if (nSingleEdges == 0) continue;
153             singleEdges.setSize(nSingleEdges);
155 #           ifdef DEBUG_ZIPUP
156             Info<< "Cell " << cellI << endl;
158             forAll(curFaces, faceI)
159             {
160                 Info<< "face: " << oldFaces[curFaces[faceI]] << endl;
161             }
163             Info<< "Cell edges: " << cellEdges << nl
164                 << "Edge usage: " << edgeUsage << nl
165                 << "Single edges: " << singleEdges << nl
166                 << "Cell points: " << cellPoints << endl;
168             forAll(cellPoints, cpI)
169             {
170                 Info<< "vertex create \"" << cellPoints[cpI]
171                     << "\" coordinates "
172                     << points()[cellPoints[cpI]] << endl;
173             }
174 #           endif
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
180             labelList pointUsage(cellPoints.size(), 0);
182             forAll(singleEdges, edgeI)
183             {
184                 const edge& curEdge = singleEdges[edgeI];
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             }
199             boolList singleEdgeUsage(singleEdges.size(), false);
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;
208                 label newEdgeStart = singleEdges[edgeI].start();
209                 label newEdgeEnd = singleEdges[edgeI].end();
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                 }
230                 if (blockedHead && blockedTail)
231                 {
232                     // Eliminating edge singleEdges[edgeI] as blocked
233                     singleEdgeUsage[edgeI] = true;
234                 }
235             }
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
241             labelListList edgesToInsert(singleEdges.size());
242             label nEdgesToInsert = 0;
244             // Find a good edge
245             forAll(singleEdges, edgeI)
246             {
247                 SLList<label> pointChain;
249                 bool blockHead = false;
250                 bool blockTail = false;
252                 if (!singleEdgeUsage[edgeI])
253                 {
254                     // found a new edge
255                     singleEdgeUsage[edgeI] = true;
257                     label newEdgeStart = singleEdges[edgeI].start();
258                     label newEdgeEnd = singleEdges[edgeI].end();
260                     pointChain.insert(newEdgeStart);
261                     pointChain.append(newEdgeEnd);
263 #                   ifdef DEBUG_CHAIN
264                     Info<< "found edge to start with: "
265                         << singleEdges[edgeI] << endl;
266 #                   endif
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
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
290                                 blockTail = true;
291                             }
292                         }
293                     }
295                     bool stopSearching = false;
297                     // Go through the unused edges and try to chain them up
298                     do
299                     {
300                         stopSearching = false;
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();
310                                 label addEnd =
311                                     singleEdges[addEdgeI].end();
313 #                               ifdef DEBUG_CHAIN
314                                 Info<< "Trying candidate "
315                                     << singleEdges[addEdgeI] << endl;
316 #                               endif
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);
326                                         singleEdgeUsage[addEdgeI] = true;
327                                     }
328                                     else if (pointChain.first() == addEnd)
329                                     {
330                                         pointChain.insert(addStart);
332                                         singleEdgeUsage[addEdgeI] = true;
333                                     }
334                                 }
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);
345                                         singleEdgeUsage[addEdgeI] = true;
346                                     }
347                                     else if (pointChain.last() == addEnd)
348                                     {
349                                         pointChain.append(addStart);
351                                         singleEdgeUsage[addEdgeI] = true;
352                                     }
353                                 }
355                                 // check if the new head or tail are blocked
356                                 label curEdgeStart = pointChain.first();
357                                 label curEdgeEnd = pointChain.last();
359 #                               ifdef DEBUG_CHAIN
360                                 Info<< "curEdgeStart: " << curEdgeStart
361                                     << " curEdgeEnd: " << curEdgeEnd << endl;
362 #                               endif
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
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
385                                             blockTail = true;
386                                         }
387                                     }
388                                 }
390                                 // Check if the loop is closed
391                                 if (curEdgeStart == curEdgeEnd)
392                                 {
393 #                                   ifdef DEBUG_CHAIN
394                                     Info<< "closed loop" << endl;
395 #                                   endif
397                                     pointChain.removeHead();
399                                     blockHead = true;
400                                     blockTail = true;
402                                     stopSearching = true;
403                                 }
405 #                               ifdef DEBUG_CHAIN
406                                 Info<< "current pointChain: " << pointChain
407                                     << endl;
408 #                               endif
410                                 if (stopSearching) break;
411                             }
412                         }
413                     } while (stopSearching);
414                 }
416 #               ifdef DEBUG_CHAIN
417                 Info<< "completed patch chain: " << pointChain << endl;
418 #               endif
420                 if (pointChain.size() > 2)
421                 {
422                     edgesToInsert[nEdgesToInsert] = pointChain;
423                     nEdgesToInsert++;
424                 }
425             }
427             edgesToInsert.setSize(nEdgesToInsert);
429 #           ifdef DEBUG_ZIPUP
430             Info<< "edgesToInsert: " << edgesToInsert << endl;
431 #           endif
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                 //
442                 const labelList& unorderedEdge = edgesToInsert[edgeToInsertI];
444                 scalarField dist(unorderedEdge.size());
446                 // Calculate distance
447                 point startPoint = Points[unorderedEdge[0]];
448                 dist[0] = 0;
450                 vector dir = Points[unorderedEdge.last()] - startPoint;
452                 for (label i = 1; i < dist.size(); i++)
453                 {
454                     dist[i] = (Points[unorderedEdge[i]] - startPoint) & dir;
455                 }
457                 // Sort points
458                 labelList orderedEdge(unorderedEdge.size(), -1);
459                 boolList used(unorderedEdge.size(), false);
461                 forAll(orderedEdge, epI)
462                 {
463                     label nextPoint = -1;
464                     scalar minDist = GREAT;
466                     forAll(dist, i)
467                     {
468                         if (!used[i] && dist[i] < minDist)
469                         {
470                             minDist = dist[i];
471                             nextPoint = i;
472                         }
473                     }
475                     // Insert the next point
476                     orderedEdge[epI] = unorderedEdge[nextPoint];
477                     used[nextPoint] = true;
478                 }
480 #               ifdef DEBUG_ORDER
481                 Info<< "unorderedEdge: " << unorderedEdge << nl
482                     << "orderedEdge: " << orderedEdge << endl;
483 #               endif
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                         {
497                             WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
498                                 << "Duplicate point found in edge to insert. "
499                                 << nl << "Point: " << orderedEdge[checkI]
500                                 << " edge: " << orderedEdge << endl;
502                             problemCells.insert(cellI);
503                         }
504                     }
505                 }
507                 edge testEdge
508                 (
509                     orderedEdge[0],
510                     orderedEdge.last()
511                 );
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()];
518                 labelList facesSharingEdge(startPF.size() + endPF.size());
519                 label nfse = 0;
521                 forAll(startPF, pfI)
522                 {
523                     facesSharingEdge[nfse++] = startPF[pfI];
524                 }
525                 forAll(endPF, pfI)
526                 {
527                     facesSharingEdge[nfse++] = endPF[pfI];
528                 }
530                 forAll(facesSharingEdge, faceI)
531                 {
532                     bool faceChanges = false;
534                     // Label of the face being analysed
535                     const label currentFaceIndex = facesSharingEdge[faceI];
537                     const edgeList curFaceEdges =
538                         oldFaces[currentFaceIndex].edges();
540                     forAll(curFaceEdges, cfeI)
541                     {
542                         if (curFaceEdges[cfeI] == testEdge)
543                         {
544                             faceChanges = true;
545                             break;
546                         }
547                     }
549                     if (faceChanges)
550                     {
551                         nChangedFacesInMesh++;
552                         // In order to avoid loosing 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. ???
564                         face& newFace = newFaces[currentFaceIndex];
566                         bool allPointsPresent = true;
568                         forAll(orderedEdge, oeI)
569                         {
570                             bool curPointFound = false;
572                             forAll(newFace, nfI)
573                             {
574                                 if (newFace[nfI] == orderedEdge[oeI])
575                                 {
576                                     curPointFound = true;
577                                     break;
578                                 }
579                             }
581                             allPointsPresent =
582                                 allPointsPresent && curPointFound;
583                         }
585 #                       ifdef DEBUG_ZIPUP
586                         if (allPointsPresent)
587                         {
588                             Info<< "All points present" << endl;
589                         }
590 #                       endif
592                         if (!allPointsPresent)
593                         {
594                             // Not all points are already present.  The
595                             // new edge will need to be inserted into the
596                             // face.
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();
604 #                           ifdef DEBUG_ZIPUP
605                             Info<< "Not all points present." << endl;
606 #                           endif
608                             label nNewFacePoints = 0;
610                             bool edgeAdded = false;
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;
620                                     // Resize the face to accept additional
621                                     // points
622                                     newFace.setSize
623                                     (
624                                         newFace.size()
625                                       + orderedEdge.size() - 2
626                                     );
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                             }
673 #                           ifdef DEBUG_ZIPUP
674                             Info<< "oldFace: "
675                                 << oldFaces[currentFaceIndex] << nl
676                                 << "newFace: " << newFace << endl;
677 #                           endif
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                                     {
691                                         WarningIn
692                                         (
693                                             "void polyMeshZipUpCells"
694                                             "("
695                                                 "polyMesh& mesh"
696                                             ")"
697                                         )
698                                             << "Duplicate point found "
699                                             << "in the new face. " << nl
700                                             << "Point: "
701                                             << orderedEdge[checkI]
702                                             << " face: "
703                                             << newFace << endl;
705                                         problemCells.insert(cellI);
706                                     }
707                                 }
708                             }
710                             // Check if the edge is added.
711                             // If not, then it comes on top of an already
712                             // modified edge and they need to be
713                             // merged in together.
714                             if (!edgeAdded)
715                             {
716                                 Info<< "This edge modifies an already modified "
717                                     << "edge.  Point insertions skipped."
718                                     << endl;
719                             }
720                         }
721                     }
722                 }
723             }
724         }
726         if (problemCells.size())
727         {
728             // This cycle has failed.  Print out the problem cells
729             labelList toc(problemCells.toc());
730             sort(toc);
732             FatalErrorIn("void polyMeshZipUpCells(polyMesh& mesh)")
733                 << "Found " << problemCells.size() << " problem cells." << nl
734                 << "Cells: " << toc
735                 << abort(FatalError);
736         }
738         Info<< "Cycle " << ++nCycles
739             << " changed " << nChangedFacesInMesh << " faces." << endl;
742         const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
744         // Reset the polyMesh. Number of points/faces/cells/patches stays the
745         // same, only the faces themselves have changed so clear all derived
746         // (edge, point) addressing.
748         // Collect the patch sizes
749         labelList patchSizes(bMesh.size(), 0);
750         labelList patchStarts(bMesh.size(), 0);
752         forAll(bMesh, patchI)
753         {
754             patchSizes[patchI] = bMesh[patchI].size();
755             patchStarts[patchI] = bMesh[patchI].start();
756         }
758         // Reset the mesh. Number of active faces is one beyond the last patch
759         // (patches guaranteed to be in increasing order)
760         mesh.resetPrimitives
761         (
762             Xfer<pointField>::null(),
763             xferMove(newFaces),
764             Xfer<labelList>::null(),
765             Xfer<labelList>::null(),
766             patchSizes,
767             patchStarts,
768             true                // boundary forms valid boundary mesh.
769         );
771         // Reset any addressing on face zones.
772         mesh.faceZones().clearAddressing();
774         // Clear the addressing
775         mesh.clearOut();
777     } while (nChangedFacesInMesh > 0 || nCycles > 100);
779     // Flags the mesh files as being changed
780     mesh.setInstance(mesh.time().timeName());
782     if (nChangedFacesInMesh > 0)
783     {
784         FatalErrorIn("void polyMeshZipUpCells(polyMesh& mesh)")
785             << "cell zip-up failed after 100 cycles.  Probable problem "
786             << "with the original mesh"
787             << abort(FatalError);
788     }
790     return nCycles != 1;
794 // ************************************************************************* //