1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright held by original author
7 -------------------------------------------------------------------------------
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "polyMeshZipUpCells.H"
31 // #define DEBUG_ZIPUP 1
32 // #define DEBUG_CHAIN 1
33 // #define DEBUG_ORDER 1
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 bool Foam::polyMeshZipUpCells(polyMesh& mesh)
41 Info<< "bool polyMeshZipUpCells(polyMesh& mesh) const: "
42 << "zipping up topologically open cells" << endl;
46 // Take the original mesh and visit all cells. For every cell
47 // calculate the edges of all faces on the cells. A cell is
48 // correctly topologically closed when all the edges are referenced
49 // by exactly two faces. If the edges are referenced only by a
50 // single face, additional vertices need to be inserted into some
51 // of the faces (topological closedness). If an edge is
52 // referenced by more that two faces, there is an error in
53 // topological closedness.
54 // Point insertion into the faces is done by attempting to create
55 // closed loops and inserting the intermediate points into the
58 // The algorithm is recursive and changes the mesh faces in each
59 // pass. It is therefore essential to discard the addressing
60 // after every pass. The algorithm is completed when the mesh
63 label nChangedFacesInMesh = 0;
66 labelHashSet problemCells;
70 nChangedFacesInMesh = 0;
72 const cellList& Cells = mesh.cells();
73 const pointField& Points = mesh.allPoints();
75 faceList newFaces = mesh.allFaces();
77 const faceList& oldFaces = mesh.allFaces();
78 const labelListList& pFaces = mesh.pointFaces();
82 const labelList& curFaces = Cells[cellI];
83 const edgeList cellEdges = Cells[cellI].edges(oldFaces);
84 const labelList cellPoints = Cells[cellI].labels(oldFaces);
86 // Find the edges used only once in the cell
88 labelList edgeUsage(cellEdges.size(), 0);
90 forAll (curFaces, faceI)
92 edgeList curFaceEdges = oldFaces[curFaces[faceI]].edges();
94 forAll (curFaceEdges, faceEdgeI)
96 const edge& curEdge = curFaceEdges[faceEdgeI];
98 forAll (cellEdges, cellEdgeI)
100 if (cellEdges[cellEdgeI] == curEdge)
102 edgeUsage[cellEdgeI]++;
109 edgeList singleEdges(cellEdges.size());
110 label nSingleEdges = 0;
112 forAll (edgeUsage, edgeI)
114 if (edgeUsage[edgeI] == 1)
116 singleEdges[nSingleEdges] = cellEdges[edgeI];
119 else if (edgeUsage[edgeI] != 2)
121 WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
122 << "edge " << cellEdges[edgeI] << " in cell " << cellI
123 << " used " << edgeUsage[edgeI] << " times. " << nl
124 << "Should be 1 or 2 - serious error "
125 << "in mesh structure. " << endl;
128 forAll (curFaces, faceI)
130 Info<< "face: " << oldFaces[curFaces[faceI]]
134 Info<< "Cell edges: " << cellEdges << nl
135 << "Edge usage: " << edgeUsage << nl
136 << "Cell points: " << cellPoints << endl;
138 forAll (cellPoints, cpI)
140 Info<< "vertex create \"" << cellPoints[cpI]
142 << Points[cellPoints[cpI]] << endl;
146 // Gather the problem cell
147 problemCells.insert(cellI);
151 // Check if the cell is already zipped up
152 if (nSingleEdges == 0) continue;
154 singleEdges.setSize(nSingleEdges);
157 Info << "Cell " << cellI << endl;
159 forAll (curFaces, faceI)
161 Info<< "face: " << oldFaces[curFaces[faceI]] << endl;
164 Info<< "Cell edges: " << cellEdges << nl
165 << "Edge usage: " << edgeUsage << nl
166 << "Single edges: " << singleEdges << nl
167 << "Cell points: " << cellPoints << endl;
169 forAll (cellPoints, cpI)
171 Info<< "vertex create \"" << cellPoints[cpI]
173 << points()[cellPoints[cpI]] << endl;
177 // Loop through all single edges and mark the points they use
178 // points marked twice are internal to edge; those marked more than
181 labelList pointUsage(cellPoints.size(), 0);
183 forAll (singleEdges, edgeI)
185 const edge& curEdge = singleEdges[edgeI];
187 forAll (cellPoints, pointI)
191 cellPoints[pointI] == curEdge.start()
192 || cellPoints[pointI] == curEdge.end()
195 pointUsage[pointI]++;
200 boolList singleEdgeUsage(singleEdges.size(), false);
202 // loop through all edges and eliminate the ones that are
204 forAll (singleEdges, edgeI)
206 bool blockedHead = false;
207 bool blockedTail = false;
209 label newEdgeStart = singleEdges[edgeI].start();
210 label newEdgeEnd = singleEdges[edgeI].end();
212 // check that the edge has not got all ends blocked
213 forAll (cellPoints, pointI)
215 if (cellPoints[pointI] == newEdgeStart)
217 if (pointUsage[pointI] > 2)
222 else if (cellPoints[pointI] == newEdgeEnd)
224 if (pointUsage[pointI] > 2)
231 if (blockedHead && blockedTail)
233 // Eliminating edge singleEdges[edgeI] as blocked
234 singleEdgeUsage[edgeI] = true;
238 // Go through the points and start from the point used twice
239 // check all the edges to find the edges starting from this point
242 labelListList edgesToInsert(singleEdges.size());
243 label nEdgesToInsert = 0;
246 forAll (singleEdges, edgeI)
248 SLList<label> pointChain;
250 bool blockHead = false;
251 bool blockTail = false;
253 if (!singleEdgeUsage[edgeI])
256 singleEdgeUsage[edgeI] = true;
258 label newEdgeStart = singleEdges[edgeI].start();
259 label newEdgeEnd = singleEdges[edgeI].end();
261 pointChain.insert(newEdgeStart);
262 pointChain.append(newEdgeEnd);
265 Info<< "found edge to start with: "
266 << singleEdges[edgeI] << endl;
269 // Check if head or tail are blocked
270 forAll (cellPoints, pointI)
272 if (cellPoints[pointI] == newEdgeStart)
274 if (pointUsage[pointI] > 2)
277 Info << "start head blocked" << endl;
283 else if(cellPoints[pointI] == newEdgeEnd)
285 if (pointUsage[pointI] > 2)
288 Info << "start tail blocked" << endl;
296 bool stopSearching = false;
298 // Go through the unused edges and try to chain them up
301 stopSearching = false;
303 forAll (singleEdges, addEdgeI)
305 if (!singleEdgeUsage[addEdgeI])
307 // Grab start and end of the candidate
309 singleEdges[addEdgeI].start();
312 singleEdges[addEdgeI].end();
315 Info<< "Trying candidate "
316 << singleEdges[addEdgeI] << endl;
319 // Try to add the edge onto the head
322 if (pointChain.first() == addStart)
324 // Added at start mark as used
325 pointChain.insert(addEnd);
327 singleEdgeUsage[addEdgeI] = true;
329 else if (pointChain.first() == addEnd)
331 pointChain.insert(addStart);
333 singleEdgeUsage[addEdgeI] = true;
337 // Try the other end only if the first end
339 if (!blockTail && !singleEdgeUsage[addEdgeI])
341 if (pointChain.last() == addStart)
343 // Added at start mark as used
344 pointChain.append(addEnd);
346 singleEdgeUsage[addEdgeI] = true;
348 else if (pointChain.last() == addEnd)
350 pointChain.append(addStart);
352 singleEdgeUsage[addEdgeI] = true;
356 // check if the new head or tail are blocked
357 label curEdgeStart = pointChain.first();
358 label curEdgeEnd = pointChain.last();
361 Info<< "curEdgeStart: " << curEdgeStart
362 << " curEdgeEnd: " << curEdgeEnd << endl;
365 forAll (cellPoints, pointI)
367 if (cellPoints[pointI] == curEdgeStart)
369 if (pointUsage[pointI] > 2)
372 Info << "head blocked" << endl;
378 else if(cellPoints[pointI] == curEdgeEnd)
380 if (pointUsage[pointI] > 2)
383 Info << "tail blocked" << endl;
391 // Check if the loop is closed
392 if (curEdgeStart == curEdgeEnd)
395 Info << "closed loop" << endl;
398 pointChain.removeHead();
403 stopSearching = true;
407 Info<< "current pointChain: " << pointChain
411 if (stopSearching) break;
414 } while (stopSearching);
418 Info << "completed patch chain: " << pointChain << endl;
421 if (pointChain.size() > 2)
423 edgesToInsert[nEdgesToInsert] = pointChain;
428 edgesToInsert.setSize(nEdgesToInsert);
431 Info << "edgesToInsert: " << edgesToInsert << endl;
434 // Insert the edges into a list of faces
435 forAll (edgesToInsert, edgeToInsertI)
437 // Order the points of the edge
438 // Warning: the ordering must be parametric, because in
439 // the case of multiple point insertion onto the same edge
440 // it is possible to get non-cyclic loops
443 const labelList& unorderedEdge = edgesToInsert[edgeToInsertI];
445 scalarField dist(unorderedEdge.size());
447 // Calculate distance
448 point startPoint = Points[unorderedEdge[0]];
452 Points[unorderedEdge[unorderedEdge.size() - 1]]
455 for (label i = 1; i < dist.size(); i++)
457 dist[i] = (Points[unorderedEdge[i]] - startPoint) & dir;
461 labelList orderedEdge(unorderedEdge.size(), -1);
462 boolList used(unorderedEdge.size(), false);
464 forAll (orderedEdge, epI)
466 label nextPoint = -1;
467 scalar minDist = GREAT;
471 if (!used[i] && dist[i] < minDist)
478 // Insert the next point
479 orderedEdge[epI] = unorderedEdge[nextPoint];
480 used[nextPoint] = true;
484 Info<< "unorderedEdge: " << unorderedEdge << nl
485 << "orderedEdge: " << orderedEdge << endl;
488 // check for duplicate points in the ordered edge
489 forAll (orderedEdge, checkI)
493 label checkJ = checkI + 1;
494 checkJ < orderedEdge.size();
498 if (orderedEdge[checkI] == orderedEdge[checkJ])
500 WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
501 << "Duplicate point found in edge to insert. "
502 << nl << "Point: " << orderedEdge[checkI]
503 << " edge: " << orderedEdge << endl;
505 problemCells.insert(cellI);
513 orderedEdge[orderedEdge.size() - 1]
516 // In order to avoid edge-to-edge comparison, get faces using
517 // point-face addressing in two goes.
518 const labelList& startPF = pFaces[testEdge.start()];
519 const labelList& endPF = pFaces[testEdge.start()];
521 labelList facesSharingEdge(startPF.size() + endPF.size());
524 forAll (startPF, pfI)
526 facesSharingEdge[nfse++] = startPF[pfI];
530 facesSharingEdge[nfse++] = endPF[pfI];
533 forAll (facesSharingEdge, faceI)
535 bool faceChanges = false;
537 // Label of the face being analysed
538 const label currentFaceIndex = facesSharingEdge[faceI];
540 const edgeList curFaceEdges =
541 oldFaces[currentFaceIndex].edges();
543 forAll (curFaceEdges, cfeI)
545 if (curFaceEdges[cfeI] == testEdge)
554 nChangedFacesInMesh++;
555 // In order to avoid loosing point from multiple
556 // insertions into the same face, the new face
557 // will be change incrementally.
558 // 1) Check if all the internal points of the edge
559 // to add already exist in the face. If so, the
560 // edge has already been included 2) Check if the
561 // point insertion occurs on an edge which is
562 // still untouched. If so, simply insert
563 // additional points into the face. 3) If not,
564 // the edge insertion occurs on an already
565 // modified edge. ???
567 face& newFace = newFaces[currentFaceIndex];
569 bool allPointsPresent = true;
571 forAll (orderedEdge, oeI)
573 bool curPointFound = false;
575 forAll (newFace, nfI)
577 if (newFace[nfI] == orderedEdge[oeI])
579 curPointFound = true;
585 allPointsPresent && curPointFound;
589 if (allPointsPresent)
591 Info << "All points present" << endl;
595 if (!allPointsPresent)
597 // Not all points are already present. The
598 // new edge will need to be inserted into the
601 // Check to see if a new edge fits onto an
602 // untouched edge of the face. Make sure the
603 // edges are grabbed before the face is
605 edgeList newFaceEdges = newFace.edges();
608 Info << "Not all points present." << endl;
611 label nNewFacePoints = 0;
613 bool edgeAdded = false;
615 forAll (newFaceEdges, curFacEdgI)
617 // Does the current edge change?
618 if (newFaceEdges[curFacEdgI] == testEdge)
620 // Found an edge match
623 // Resize the face to accept additional
628 + orderedEdge.size() - 2
633 newFaceEdges[curFacEdgI].start()
637 // insertion in ascending order
641 i < orderedEdge.size() - 1;
645 newFace[nNewFacePoints] =
652 // insertion in reverse order
655 label i = orderedEdge.size() - 1;
660 newFace[nNewFacePoints] =
668 // Does not fit onto this edge.
669 // Copy the next point into the face
670 newFace[nNewFacePoints] =
671 newFaceEdges[curFacEdgI].start();
678 << oldFaces[currentFaceIndex] << nl
679 << "newFace: " << newFace << endl;
682 // Check for duplicate points in the new face
683 forAll (newFace, checkI)
687 label checkJ = checkI + 1;
688 checkJ < newFace.size();
692 if (newFace[checkI] == newFace[checkJ])
694 WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
695 << "Duplicate point found "
696 << "in the new face. " << nl
698 << orderedEdge[checkI]
702 problemCells.insert(cellI);
707 // Check if the edge is added.
708 // If not, then it comes on top of an already
709 // modified edge and they need to be
710 // merged in together.
713 Info<< "This edge modifies an already modified "
714 << "edge. Point insertions skipped."
723 if (problemCells.size())
725 // This cycle has failed. Print out the problem cells
726 labelList toc(problemCells.toc());
729 FatalErrorIn("void polyMeshZipUpCells(polyMesh& mesh)")
730 << "Found " << problemCells.size() << " problem cells." << nl
732 << abort(FatalError);
735 Info<< "Cycle " << ++nCycles
736 << " changed " << nChangedFacesInMesh << " faces." << endl;
739 const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
741 // Reset the polyMesh. Number of points/faces/cells/patches stays the
742 // same, only the faces themselves have changed so clear all derived
743 // (edge, point) addressing.
745 // Collect the patch sizes
746 labelList patchSizes(bMesh.size(), 0);
747 labelList patchStarts(bMesh.size(), 0);
749 forAll (bMesh, patchI)
751 patchSizes[patchI] = bMesh[patchI].size();
752 patchStarts[patchI] = bMesh[patchI].start();
755 // Reset the mesh. Number of active faces is one beyond the last patch
756 // (patches guaranteed to be in increasing order)
757 // 1.6.x merge change. Reconsider. HJ, 22/Aug/2010
760 Xfer<pointField>::null(),
762 Xfer<labelList>::null(),
763 Xfer<labelList>::null(),
766 true // boundary forms valid boundary mesh.
769 // Clear the addressing
772 } while (nChangedFacesInMesh > 0 || nCycles > 100);
774 // Flags the mesh files as being changed
775 mesh.setInstance(mesh.time().timeName());
777 if (nChangedFacesInMesh > 0)
779 FatalErrorIn("void polyMeshZipUpCells(polyMesh& mesh)")
780 << "cell zip-up failed after 100 cycles. Probable problem "
781 << "with the original mesh"
782 << abort(FatalError);
789 // ************************************************************************* //