1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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
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
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 <OpenFOAM/polyMesh.H>
28 #include <OpenFOAM/Time.H>
30 // #define DEBUG_ZIPUP 1
31 // #define DEBUG_CHAIN 1
32 // #define DEBUG_ORDER 1
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 bool Foam::polyMeshZipUpCells(polyMesh& mesh)
40 Info<< "bool polyMeshZipUpCells(polyMesh& mesh) const: "
41 << "zipping up topologically open cells" << endl;
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
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
62 label nChangedFacesInMesh = 0;
65 labelHashSet problemCells;
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();
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)
91 edgeList curFaceEdges = oldFaces[curFaces[faceI]].edges();
93 forAll (curFaceEdges, faceEdgeI)
95 const edge& curEdge = curFaceEdges[faceEdgeI];
97 forAll (cellEdges, cellEdgeI)
99 if (cellEdges[cellEdgeI] == curEdge)
101 edgeUsage[cellEdgeI]++;
108 edgeList singleEdges(cellEdges.size());
109 label nSingleEdges = 0;
111 forAll (edgeUsage, edgeI)
113 if (edgeUsage[edgeI] == 1)
115 singleEdges[nSingleEdges] = cellEdges[edgeI];
118 else if (edgeUsage[edgeI] != 2)
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;
127 forAll (curFaces, faceI)
129 Info<< "face: " << oldFaces[curFaces[faceI]]
133 Info<< "Cell edges: " << cellEdges << nl
134 << "Edge usage: " << edgeUsage << nl
135 << "Cell points: " << cellPoints << endl;
137 forAll (cellPoints, cpI)
139 Info<< "vertex create \"" << cellPoints[cpI]
141 << Points[cellPoints[cpI]] << endl;
145 // Gather the problem cell
146 problemCells.insert(cellI);
150 // Check if the cell is already zipped up
151 if (nSingleEdges == 0) continue;
153 singleEdges.setSize(nSingleEdges);
156 Info << "Cell " << cellI << endl;
158 forAll (curFaces, faceI)
160 Info<< "face: " << oldFaces[curFaces[faceI]] << endl;
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)
170 Info<< "vertex create \"" << cellPoints[cpI]
172 << points()[cellPoints[cpI]] << endl;
176 // Loop through all single edges and mark the points they use
177 // points marked twice are internal to edge; those marked more than
180 labelList pointUsage(cellPoints.size(), 0);
182 forAll (singleEdges, edgeI)
184 const edge& curEdge = singleEdges[edgeI];
186 forAll (cellPoints, pointI)
190 cellPoints[pointI] == curEdge.start()
191 || cellPoints[pointI] == curEdge.end()
194 pointUsage[pointI]++;
199 boolList singleEdgeUsage(singleEdges.size(), false);
201 // loop through all edges and eliminate the ones that are
203 forAll (singleEdges, edgeI)
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)
214 if (cellPoints[pointI] == newEdgeStart)
216 if (pointUsage[pointI] > 2)
221 else if (cellPoints[pointI] == newEdgeEnd)
223 if (pointUsage[pointI] > 2)
230 if (blockedHead && blockedTail)
232 // Eliminating edge singleEdges[edgeI] as blocked
233 singleEdgeUsage[edgeI] = true;
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
241 labelListList edgesToInsert(singleEdges.size());
242 label nEdgesToInsert = 0;
245 forAll (singleEdges, edgeI)
247 SLList<label> pointChain;
249 bool blockHead = false;
250 bool blockTail = false;
252 if (!singleEdgeUsage[edgeI])
255 singleEdgeUsage[edgeI] = true;
257 label newEdgeStart = singleEdges[edgeI].start();
258 label newEdgeEnd = singleEdges[edgeI].end();
260 pointChain.insert(newEdgeStart);
261 pointChain.append(newEdgeEnd);
264 Info<< "found edge to start with: "
265 << singleEdges[edgeI] << endl;
268 // Check if head or tail are blocked
269 forAll (cellPoints, pointI)
271 if (cellPoints[pointI] == newEdgeStart)
273 if (pointUsage[pointI] > 2)
276 Info << "start head blocked" << endl;
282 else if(cellPoints[pointI] == newEdgeEnd)
284 if (pointUsage[pointI] > 2)
287 Info << "start tail blocked" << endl;
295 bool stopSearching = false;
297 // Go through the unused edges and try to chain them up
300 stopSearching = false;
302 forAll (singleEdges, addEdgeI)
304 if (!singleEdgeUsage[addEdgeI])
306 // Grab start and end of the candidate
308 singleEdges[addEdgeI].start();
311 singleEdges[addEdgeI].end();
314 Info<< "Trying candidate "
315 << singleEdges[addEdgeI] << endl;
318 // Try to add the edge onto the head
321 if (pointChain.first() == addStart)
323 // Added at start mark as used
324 pointChain.insert(addEnd);
326 singleEdgeUsage[addEdgeI] = true;
328 else if (pointChain.first() == addEnd)
330 pointChain.insert(addStart);
332 singleEdgeUsage[addEdgeI] = true;
336 // Try the other end only if the first end
338 if (!blockTail && !singleEdgeUsage[addEdgeI])
340 if (pointChain.last() == addStart)
342 // Added at start mark as used
343 pointChain.append(addEnd);
345 singleEdgeUsage[addEdgeI] = true;
347 else if (pointChain.last() == addEnd)
349 pointChain.append(addStart);
351 singleEdgeUsage[addEdgeI] = true;
355 // check if the new head or tail are blocked
356 label curEdgeStart = pointChain.first();
357 label curEdgeEnd = pointChain.last();
360 Info<< "curEdgeStart: " << curEdgeStart
361 << " curEdgeEnd: " << curEdgeEnd << endl;
364 forAll (cellPoints, pointI)
366 if (cellPoints[pointI] == curEdgeStart)
368 if (pointUsage[pointI] > 2)
371 Info << "head blocked" << endl;
377 else if(cellPoints[pointI] == curEdgeEnd)
379 if (pointUsage[pointI] > 2)
382 Info << "tail blocked" << endl;
390 // Check if the loop is closed
391 if (curEdgeStart == curEdgeEnd)
394 Info << "closed loop" << endl;
397 pointChain.removeHead();
402 stopSearching = true;
406 Info<< "current pointChain: " << pointChain
410 if (stopSearching) break;
413 } while (stopSearching);
417 Info << "completed patch chain: " << pointChain << endl;
420 if (pointChain.size() > 2)
422 edgesToInsert[nEdgesToInsert] = pointChain;
427 edgesToInsert.setSize(nEdgesToInsert);
430 Info << "edgesToInsert: " << edgesToInsert << endl;
433 // Insert the edges into a list of faces
434 forAll (edgesToInsert, edgeToInsertI)
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
442 const labelList& unorderedEdge = edgesToInsert[edgeToInsertI];
444 scalarField dist(unorderedEdge.size());
446 // Calculate distance
447 point startPoint = Points[unorderedEdge[0]];
451 Points[unorderedEdge[unorderedEdge.size() - 1]]
454 for (label i = 1; i < dist.size(); i++)
456 dist[i] = (Points[unorderedEdge[i]] - startPoint) & dir;
460 labelList orderedEdge(unorderedEdge.size(), -1);
461 boolList used(unorderedEdge.size(), false);
463 forAll (orderedEdge, epI)
465 label nextPoint = -1;
466 scalar minDist = GREAT;
470 if (!used[i] && dist[i] < minDist)
477 // Insert the next point
478 orderedEdge[epI] = unorderedEdge[nextPoint];
479 used[nextPoint] = true;
483 Info<< "unorderedEdge: " << unorderedEdge << nl
484 << "orderedEdge: " << orderedEdge << endl;
487 // check for duplicate points in the ordered edge
488 forAll (orderedEdge, checkI)
492 label checkJ = checkI + 1;
493 checkJ < orderedEdge.size();
497 if (orderedEdge[checkI] == orderedEdge[checkJ])
499 WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
500 << "Duplicate point found in edge to insert. "
501 << nl << "Point: " << orderedEdge[checkI]
502 << " edge: " << orderedEdge << endl;
504 problemCells.insert(cellI);
512 orderedEdge[orderedEdge.size() - 1]
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()];
520 labelList facesSharingEdge(startPF.size() + endPF.size());
523 forAll (startPF, pfI)
525 facesSharingEdge[nfse++] = startPF[pfI];
529 facesSharingEdge[nfse++] = endPF[pfI];
532 forAll (facesSharingEdge, faceI)
534 bool faceChanges = false;
536 // Label of the face being analysed
537 const label currentFaceIndex = facesSharingEdge[faceI];
539 const edgeList curFaceEdges =
540 oldFaces[currentFaceIndex].edges();
542 forAll (curFaceEdges, cfeI)
544 if (curFaceEdges[cfeI] == testEdge)
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. ???
566 face& newFace = newFaces[currentFaceIndex];
568 bool allPointsPresent = true;
570 forAll (orderedEdge, oeI)
572 bool curPointFound = false;
574 forAll (newFace, nfI)
576 if (newFace[nfI] == orderedEdge[oeI])
578 curPointFound = true;
584 allPointsPresent && curPointFound;
588 if (allPointsPresent)
590 Info << "All points present" << endl;
594 if (!allPointsPresent)
596 // Not all points are already present. The
597 // new edge will need to be inserted into the
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
604 edgeList newFaceEdges = newFace.edges();
607 Info << "Not all points present." << endl;
610 label nNewFacePoints = 0;
612 bool edgeAdded = false;
614 forAll (newFaceEdges, curFacEdgI)
616 // Does the current edge change?
617 if (newFaceEdges[curFacEdgI] == testEdge)
619 // Found an edge match
622 // Resize the face to accept additional
627 + orderedEdge.size() - 2
632 newFaceEdges[curFacEdgI].start()
636 // insertion in ascending order
640 i < orderedEdge.size() - 1;
644 newFace[nNewFacePoints] =
651 // insertion in reverse order
654 label i = orderedEdge.size() - 1;
659 newFace[nNewFacePoints] =
667 // Does not fit onto this edge.
668 // Copy the next point into the face
669 newFace[nNewFacePoints] =
670 newFaceEdges[curFacEdgI].start();
677 << oldFaces[currentFaceIndex] << nl
678 << "newFace: " << newFace << endl;
681 // Check for duplicate points in the new face
682 forAll (newFace, checkI)
686 label checkJ = checkI + 1;
687 checkJ < newFace.size();
691 if (newFace[checkI] == newFace[checkJ])
693 WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
694 << "Duplicate point found "
695 << "in the new face. " << nl
697 << orderedEdge[checkI]
701 problemCells.insert(cellI);
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.
712 Info<< "This edge modifies an already modified "
713 << "edge. Point insertions skipped."
722 if (problemCells.size())
724 // This cycle has failed. Print out the problem cells
725 labelList toc(problemCells.toc());
728 FatalErrorIn("void polyMeshZipUpCells(polyMesh& mesh)")
729 << "Found " << problemCells.size() << " problem cells." << nl
731 << abort(FatalError);
734 Info<< "Cycle " << ++nCycles
735 << " changed " << nChangedFacesInMesh << " faces." << endl;
738 const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
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.
744 // Collect the patch sizes
745 labelList patchSizes(bMesh.size(), 0);
746 labelList patchStarts(bMesh.size(), 0);
748 forAll (bMesh, patchI)
750 patchSizes[patchI] = bMesh[patchI].size();
751 patchStarts[patchI] = bMesh[patchI].start();
754 // Reset the mesh. Number of active faces is one beyond the last patch
755 // (patches guaranteed to be in increasing order)
758 Xfer<pointField>::null(),
760 Xfer<labelList>::null(),
761 Xfer<labelList>::null(),
764 true // boundary forms valid boundary mesh.
767 // Reset any addressing on face zones.
768 mesh.faceZones().clearAddressing();
770 // Clear the addressing
773 } while (nChangedFacesInMesh > 0 || nCycles > 100);
775 // Flags the mesh files as being changed
776 mesh.setInstance(mesh.time().timeName());
778 if (nChangedFacesInMesh > 0)
780 FatalErrorIn("void polyMeshZipUpCells(polyMesh& mesh)")
781 << "cell zip-up failed after 100 cycles. Probable problem "
782 << "with the original mesh"
783 << abort(FatalError);
790 // ************************ vim: set sw=4 sts=4 et: ************************ //