| blob | 490d7f6b9e6f5240dfbed6f613ee10c6e4276c9d |
1 /*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | Copyright held by original author
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 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
19 for more details.
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 "Stack.H"
28 #include "objectMap.H"
29 #include "changeMap.H"
30 #include "multiThreader.H"
31 #include "dynamicTopoFvMesh.H"
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 namespace Foam
36 {
38 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
40 // Method for the collapse of a quad-face in 2D
41 // - Returns a changeMap with a type specifying:
42 // -1: Collapse failed since max number of topo-changes was reached.
43 // 0: Collapse could not be performed.
44 // 1: Collapsed to first node.
45 // 2: Collapsed to second node.
46 // - overRideCase is used to force a certain collapse configuration.
47 // -1: Use this value to let collapseQuadFace decide a case.
48 // 1: Force collapse to first node.
49 // 2: Force collapse to second node.
50 // - checkOnly performs a feasibility check and returns without modifications.
51 const changeMap dynamicTopoFvMesh::collapseQuadFace
52 (
53 const label fIndex,
54 label overRideCase,
55 bool checkOnly,
56 bool forceOp
57 )
58 {
59 // Figure out which thread this is...
60 label tIndex = self();
62 // Prepare the changeMap
63 changeMap map;
65 if
66 (
67 (statistics_[0] > maxModifications_)
68 && (maxModifications_ > -1)
69 )
70 {
71 // Reached the max allowable topo-changes.
72 stack(tIndex).clear();
74 return map;
75 }
77 // Check if edgeRefinements are to be avoided on patch.
78 if (lengthEstimator().checkRefinementPatch(whichPatch(fIndex)))
79 {
80 return map;
81 }
83 // Sanity check: Is the index legitimate?
84 if (fIndex < 0)
85 {
86 FatalErrorIn
87 (
88 "\n"
89 "const changeMap "
90 "dynamicTopoFvMesh::collapseQuadFace\n"
91 "(\n"
92 " const label fIndex,\n"
93 " label overRideCase,\n"
94 " bool checkOnly\n"
95 ")\n"
96 )
97 << " Invalid index: " << fIndex
98 << abort(FatalError);
99 }
101 // Define the edges on the face to be collapsed
102 FixedList<edge,4> checkEdge(edge(-1,-1));
103 FixedList<label,4> checkEdgeIndex(-1);
105 // Define checkEdges
106 checkEdgeIndex[0] = getTriBoundaryEdge(fIndex);
107 checkEdge[0] = edges_[checkEdgeIndex[0]];
109 const labelList& fEdges = faceEdges_[fIndex];
111 forAll(fEdges, edgeI)
112 {
113 if (checkEdgeIndex[0] != fEdges[edgeI])
114 {
115 const edge& thisEdge = edges_[fEdges[edgeI]];
117 if
118 (
119 checkEdge[0].start() == thisEdge[0] ||
120 checkEdge[0].start() == thisEdge[1]
121 )
122 {
123 checkEdgeIndex[1] = fEdges[edgeI];
124 checkEdge[1] = thisEdge;
126 // Update the map
127 map.firstEdge() = checkEdgeIndex[1];
128 }
129 else
130 if
131 (
132 checkEdge[0].end() == thisEdge[0] ||
133 checkEdge[0].end() == thisEdge[1]
134 )
135 {
136 checkEdgeIndex[2] = fEdges[edgeI];
137 checkEdge[2] = thisEdge;
139 // Update the map
140 map.secondEdge() = checkEdgeIndex[2];
141 }
142 else
143 {
144 checkEdgeIndex[3] = fEdges[edgeI];
145 checkEdge[3] = thisEdge;
146 }
147 }
148 }
150 // Build a hull of cells and tri-faces that are connected to each edge
151 FixedList<labelList, 2> hullCells;
152 FixedList<labelList, 2> hullTriFaces;
154 meshOps::constructPrismHull
155 (
156 checkEdgeIndex[1],
157 faces_,
158 cells_,
159 owner_,
160 neighbour_,
161 edgeFaces_,
162 hullTriFaces[0],
163 hullCells[0]
164 );
166 meshOps::constructPrismHull
167 (
168 checkEdgeIndex[2],
169 faces_,
170 cells_,
171 owner_,
172 neighbour_,
173 edgeFaces_,
174 hullTriFaces[1],
175 hullCells[1]
176 );
178 // Determine the neighbouring cells
179 label c0 = owner_[fIndex], c1 = neighbour_[fIndex];
181 // Define variables for the prism-face calculation
182 FixedList<face,2> c0BdyFace, c0IntFace, c1BdyFace, c1IntFace;
183 FixedList<label,2> c0BdyIndex, c0IntIndex, c1BdyIndex, c1IntIndex;
185 // Find the prism-faces
186 meshOps::findPrismFaces
187 (
188 fIndex,
189 c0,
190 faces_,
191 cells_,
192 neighbour_,
193 c0BdyFace,
194 c0BdyIndex,
195 c0IntFace,
196 c0IntIndex
197 );
199 if (c1 != -1)
200 {
201 meshOps::findPrismFaces
202 (
203 fIndex,
204 c1,
205 faces_,
206 cells_,
207 neighbour_,
208 c1BdyFace,
209 c1BdyIndex,
210 c1IntFace,
211 c1IntIndex
212 );
213 }
215 // Determine if either edge belongs to a boundary
216 FixedList<bool, 2> nBoundCurves(false), edgeBoundary(false);
218 edgeBoundary[0] = (whichEdgePatch(checkEdgeIndex[1]) > -1);
219 edgeBoundary[1] = (whichEdgePatch(checkEdgeIndex[2]) > -1);
221 // Decide on collapseCase
222 label collapseCase = -1;
224 if (edgeBoundary[0] && !edgeBoundary[1])
225 {
226 collapseCase = 1;
227 }
228 else
229 if (!edgeBoundary[0] && edgeBoundary[1])
230 {
231 collapseCase = 2;
232 }
233 else
234 if (edgeBoundary[0] && edgeBoundary[1])
235 {
236 if (c1 != -1)
237 {
238 if (debug > 2)
239 {
240 WarningIn
241 (
242 "\n"
243 "const changeMap "
244 "dynamicTopoFvMesh::collapseQuadFace\n"
245 "(\n"
246 " const label fIndex,\n"
247 " label overRideCase,\n"
248 " bool checkOnly\n"
249 ")\n"
250 ) << "Collapsing an internal face that "
251 << "lies on two boundary patches. "
252 << "Face: " << fIndex << ": " << faces_[fIndex]
253 << endl;
254 }
256 // Bail out for now. If proximity based refinement is
257 // switched on, mesh may be sliced at this point.
258 return map;
259 }
261 // Check if either edge lies on a bounding curve.
262 if (checkBoundingCurve(checkEdgeIndex[1]))
263 {
264 nBoundCurves[0] = true;
265 }
267 if (checkBoundingCurve(checkEdgeIndex[2]))
268 {
269 nBoundCurves[1] = true;
270 }
272 // Collapse towards a bounding curve
273 if (nBoundCurves[0] && !nBoundCurves[1])
274 {
275 collapseCase = 1;
276 }
277 else
278 if (!nBoundCurves[0] && nBoundCurves[1])
279 {
280 collapseCase = 2;
281 }
282 else
283 if (!nBoundCurves[0] && !nBoundCurves[1])
284 {
285 // No bounding curves. Collapse to mid-point.
286 collapseCase = 3;
287 }
288 else
289 {
290 // Two bounding curves? This might cause pinching.
291 // Bail out for now.
292 return map;
293 }
294 }
295 else
296 {
297 // Looks like this is an interior face.
298 // Collapse case [3] by default
299 collapseCase = 3;
300 }
302 // Perform an override if requested.
303 if (overRideCase != -1)
304 {
305 collapseCase = overRideCase;
306 }
308 // Configure the new point-positions
309 FixedList<bool, 2> check(false);
310 FixedList<FixedList<label, 2>, 2> checkPoints(FixedList<label, 2>(-1));
311 FixedList<point, 2> newPoint(vector::zero);
312 FixedList<point, 2> oldPoint(vector::zero);
314 // Determine the common vertices for the first and second edges
315 label cv0 = checkEdge[1].commonVertex(checkEdge[0]);
316 label cv1 = checkEdge[1].commonVertex(checkEdge[3]);
317 label cv2 = checkEdge[2].commonVertex(checkEdge[0]);
318 label cv3 = checkEdge[2].commonVertex(checkEdge[3]);
320 // Replacement check points
321 FixedList<label,2> original(-1), replacement(-1);
323 switch (collapseCase)
324 {
325 case 1:
326 {
327 // Collapse to the first node
328 original[0] = cv2; original[1] = cv3;
329 replacement[0] = cv0; replacement[1] = cv1;
331 newPoint[0] = points_[replacement[0]];
332 newPoint[1] = points_[replacement[1]];
334 oldPoint[0] = oldPoints_[replacement[0]];
335 oldPoint[1] = oldPoints_[replacement[1]];
337 // Define check-points
338 check[1] = true;
339 checkPoints[1][0] = original[0];
340 checkPoints[1][1] = original[1];
342 break;
343 }
345 case 2:
346 {
347 // Collapse to the second node
348 original[0] = cv0; original[1] = cv1;
349 replacement[0] = cv2; replacement[1] = cv3;
351 newPoint[0] = points_[replacement[0]];
352 newPoint[1] = points_[replacement[1]];
354 oldPoint[0] = oldPoints_[replacement[0]];
355 oldPoint[1] = oldPoints_[replacement[1]];
357 // Define check-points
358 check[0] = true;
359 checkPoints[0][0] = original[0];
360 checkPoints[0][1] = original[1];
362 break;
363 }
365 case 3:
366 {
367 // Collapse to the mid-point
368 original[0] = cv0; original[1] = cv1;
369 replacement[0] = cv2; replacement[1] = cv3;
371 // Define new point-positions
372 newPoint[0] =
373 (
374 0.5 *
375 (
376 points_[original[0]]
377 + points_[replacement[0]]
378 )
379 );
381 newPoint[1] =
382 (
383 0.5 *
384 (
385 points_[original[1]]
386 + points_[replacement[1]]
387 )
388 );
390 // Specify off-centering
391 scalar offCentre = (c1 == -1) ? 0.0 : 1.0;
393 FixedList<vector,2> te(vector::zero), xf(vector::zero);
394 FixedList<vector,2> ne(vector::zero), nf(vector::zero);
396 // Compute tangent-to-edge
397 te[0] = (oldPoints_[replacement[0]] - oldPoints_[original[0]]);
398 te[1] = (oldPoints_[replacement[1]] - oldPoints_[original[1]]);
400 // Compute face position / normal
401 if (c0BdyFace[0].which(original[0]) > -1)
402 {
403 xf[0] = c0BdyFace[0].centre(oldPoints_);
404 nf[0] = c0BdyFace[0].normal(oldPoints_);
406 xf[1] = c0BdyFace[1].centre(oldPoints_);
407 nf[1] = c0BdyFace[1].normal(oldPoints_);
408 }
409 else
410 if (c0BdyFace[1].which(original[0]) > -1)
411 {
412 xf[0] = c0BdyFace[1].centre(oldPoints_);
413 nf[0] = c0BdyFace[1].normal(oldPoints_);
415 xf[1] = c0BdyFace[0].centre(oldPoints_);
416 nf[1] = c0BdyFace[0].normal(oldPoints_);
417 }
418 else
419 {
420 FatalErrorIn
421 (
422 "\n"
423 "const changeMap "
424 "dynamicTopoFvMesh::collapseQuadFace\n"
425 "(\n"
426 " const label fIndex,\n"
427 " label overRideCase,\n"
428 " bool checkOnly\n"
429 ")\n"
430 ) << "Could not find point in face."
431 << endl;
432 }
434 // Compute edge-normals
435 ne[0] = (te[0] ^ nf[0]);
436 ne[1] = (te[1] ^ nf[1]);
438 ne[0] /= mag(ne[0]) + VSMALL;
439 ne[1] /= mag(ne[1]) + VSMALL;
441 // Reverse the vector, if necessary
442 if ((ne[0] & ne[1]) < 0.0)
443 {
444 ne[1] *= -1.0;
445 }
447 // Define modified old point-positions,
448 // with off-centering, if necessary
449 oldPoint[0] =
450 (
451 oldPoints_[original[0]]
452 + (0.5 * te[0])
453 + (((0.05 * mag(te[0])) * ne[0]) * offCentre)
454 );
456 oldPoint[1] =
457 (
458 oldPoints_[original[1]]
459 + (0.5 * te[1])
460 + (((0.05 * mag(te[1])) * ne[1]) * offCentre)
461 );
463 // Define check-points
464 check[0] = true;
465 checkPoints[0][0] = original[0];
466 checkPoints[0][1] = original[1];
468 check[1] = true;
469 checkPoints[1][0] = replacement[0];
470 checkPoints[1][1] = replacement[1];
472 break;
473 }
475 default:
476 {
477 // Don't think this will ever happen.
478 FatalErrorIn
479 (
480 "\n"
481 "const changeMap "
482 "dynamicTopoFvMesh::collapseQuadFace\n"
483 "(\n"
484 " const label fIndex,\n"
485 " label overRideCase,\n"
486 " bool checkOnly\n"
487 ")\n"
488 )
489 << "Edge: " << fIndex << ": " << faces_[fIndex]
490 << ". Couldn't decide on collapseCase."
491 << abort(FatalError);
493 break;
494 }
495 }
497 // Keep track of resulting cell quality,
498 // if collapse is indeed feasible
499 scalar collapseQuality(GREAT);
501 // Check collapsibility of cells around edges
502 // with the re-configured point
503 forAll(check, indexI)
504 {
505 if (!check[indexI])
506 {
507 continue;
508 }
510 // Check whether the collapse is possible.
511 if
512 (
513 checkCollapse
514 (
515 hullTriFaces[indexI],
516 c0BdyIndex,
517 c1BdyIndex,
518 checkPoints[indexI],
519 newPoint,
520 oldPoint,
521 collapseQuality,
522 (c1 != -1),
523 forceOp
524 )
525 )
526 {
527 map.type() = 0;
528 return map;
529 }
530 }
532 // Are we only performing checks?
533 if (checkOnly)
534 {
535 map.type() = collapseCase;
537 if (debug > 2)
538 {
539 Info << "Face: " << fIndex
540 << ":: " << faces_[fIndex] << nl
541 << " collapseCase determined to be: "
542 << collapseCase << nl
543 << " Resulting quality: "
544 << collapseQuality
545 << endl;
546 }
548 return map;
549 }
551 if (debug > 1)
552 {
553 const labelList& fE = faceEdges_[fIndex];
555 Info << nl << nl
556 << "Face: " << fIndex << ": " << faces_[fIndex] << nl
557 << "faceEdges: " << fE
558 << " is to be collapsed. " << endl;
560 label epIndex = whichPatch(fIndex);
562 Info << "Patch: ";
564 if (epIndex == -1)
565 {
566 Info << "Internal" << endl;
567 }
568 else
569 {
570 Info << boundaryMesh()[epIndex].name() << endl;
571 }
573 if (debug > 2)
574 {
575 Info << endl;
576 Info << "~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
577 Info << "Hulls before modification" << endl;
578 Info << "~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
580 Info << nl << "Cells belonging to first Edge Hull: "
581 << hullCells[0] << endl;
583 forAll(hullCells[0],cellI)
584 {
585 const cell& firstCurCell = cells_[hullCells[0][cellI]];
587 Info << "Cell: " << hullCells[0][cellI]
588 << ": " << firstCurCell << endl;
590 forAll(firstCurCell,faceI)
591 {
592 Info << firstCurCell[faceI]
593 << ": " << faces_[firstCurCell[faceI]] << endl;
594 }
595 }
597 const labelList& firstEdgeFaces = edgeFaces_[checkEdgeIndex[1]];
599 Info << nl << "First Edge Face Hull: "
600 << firstEdgeFaces << endl;
602 forAll(firstEdgeFaces,indexI)
603 {
604 Info << firstEdgeFaces[indexI]
605 << ": " << faces_[firstEdgeFaces[indexI]] << endl;
606 }
608 Info << nl << "Cells belonging to second Edge Hull: "
609 << hullCells[1] << endl;
611 forAll(hullCells[1], cellI)
612 {
613 const cell& secondCurCell = cells_[hullCells[1][cellI]];
615 Info << "Cell: " << hullCells[1][cellI]
616 << ": " << secondCurCell << endl;
618 forAll(secondCurCell, faceI)
619 {
620 Info << secondCurCell[faceI]
621 << ": " << faces_[secondCurCell[faceI]] << endl;
622 }
623 }
625 const labelList& secondEdgeFaces = edgeFaces_[checkEdgeIndex[2]];
627 Info << nl << "Second Edge Face Hull: "
628 << secondEdgeFaces << endl;
630 forAll(secondEdgeFaces, indexI)
631 {
632 Info << secondEdgeFaces[indexI]
633 << ": " << faces_[secondEdgeFaces[indexI]] << endl;
634 }
636 // Write out VTK files prior to change
637 if (debug > 3)
638 {
639 labelHashSet vtkCells;
641 forAll(hullCells[0], cellI)
642 {
643 if (!vtkCells.found(hullCells[0][cellI]))
644 {
645 vtkCells.insert(hullCells[0][cellI]);
646 }
647 }
649 forAll(hullCells[1], cellI)
650 {
651 if (!vtkCells.found(hullCells[1][cellI]))
652 {
653 vtkCells.insert(hullCells[1][cellI]);
654 }
655 }
657 writeVTK
658 (
659 Foam::name(fIndex)
660 + "_Collapse_0",
661 vtkCells.toc()
662 );
663 }
664 }
665 }
667 // Edges / Quad-faces to throw or keep during collapse
668 FixedList<label,2> ends(-1);
669 FixedList<label,2> faceToKeep(-1), faceToThrow(-1);
670 FixedList<label,4> edgeToKeep(-1), edgeToThrow(-1);
672 // Maintain a list of modified faces for mapping
673 labelHashSet modifiedFaces;
675 // Case 2 & 3 use identical connectivity,
676 // but different point locations
677 if (collapseCase == 2 || collapseCase == 3)
678 {
679 const labelList& firstEdgeFaces = edgeFaces_[checkEdgeIndex[1]];
681 // Collapse to the second node...
682 forAll(firstEdgeFaces,faceI)
683 {
684 // Replace point indices on faces.
685 meshOps::replaceLabel
686 (
687 cv0,
688 cv2,
689 faces_[firstEdgeFaces[faceI]]
690 );
692 meshOps::replaceLabel
693 (
694 cv1,
695 cv3,
696 faces_[firstEdgeFaces[faceI]]
697 );
699 // Add an entry for mapping
700 if (!modifiedFaces.found(firstEdgeFaces[faceI]))
701 {
702 modifiedFaces.insert(firstEdgeFaces[faceI]);
703 }
705 // Determine the quad-face in cell[0] & cell[1]
706 // that belongs to firstEdgeFaces
707 if (firstEdgeFaces[faceI] == c0IntIndex[0])
708 {
709 faceToKeep[0] = c0IntIndex[1];
710 faceToThrow[0] = c0IntIndex[0];
711 }
713 if (firstEdgeFaces[faceI] == c0IntIndex[1])
714 {
715 faceToKeep[0] = c0IntIndex[0];
716 faceToThrow[0] = c0IntIndex[1];
717 }
719 if (c1 != -1)
720 {
721 if (firstEdgeFaces[faceI] == c1IntIndex[0])
722 {
723 faceToKeep[1] = c1IntIndex[1];
724 faceToThrow[1] = c1IntIndex[0];
725 }
727 if (firstEdgeFaces[faceI] == c1IntIndex[1])
728 {
729 faceToKeep[1] = c1IntIndex[0];
730 faceToThrow[1] = c1IntIndex[1];
731 }
732 }
733 }
735 // Find common edges between quad / tri faces...
736 meshOps::findCommonEdge
737 (
738 c0BdyIndex[0],
739 faceToKeep[0],
740 faceEdges_,
741 edgeToKeep[0]
742 );
744 meshOps::findCommonEdge
745 (
746 c0BdyIndex[1],
747 faceToKeep[0],
748 faceEdges_,
749 edgeToKeep[1]
750 );
752 meshOps::findCommonEdge
753 (
754 c0BdyIndex[0],
755 faceToThrow[0],
756 faceEdges_,
757 edgeToThrow[0]
758 );
760 meshOps::findCommonEdge
761 (
762 c0BdyIndex[1],
763 faceToThrow[0],
764 faceEdges_,
765 edgeToThrow[1]
766 );
768 // Size down edgeFaces for the ends.
769 meshOps::findCommonEdge
770 (
771 faceToThrow[0],
772 faceToKeep[0],
773 faceEdges_,
774 ends[0]
775 );
777 meshOps::sizeDownList
778 (
779 faceToThrow[0],
780 edgeFaces_[ends[0]]
781 );
783 if (c1 != -1)
784 {
785 meshOps::findCommonEdge
786 (
787 c1BdyIndex[0],
788 faceToKeep[1],
789 faceEdges_,
790 edgeToKeep[2]
791 );
793 meshOps::findCommonEdge
794 (
795 c1BdyIndex[1],
796 faceToKeep[1],
797 faceEdges_,
798 edgeToKeep[3]
799 );
801 meshOps::findCommonEdge
802 (
803 c1BdyIndex[0],
804 faceToThrow[1],
805 faceEdges_,
806 edgeToThrow[2]
807 );
809 meshOps::findCommonEdge
810 (
811 c1BdyIndex[1],
812 faceToThrow[1],
813 faceEdges_,
814 edgeToThrow[3]
815 );
817 // Size down edgeFaces for the ends.
818 meshOps::findCommonEdge
819 (
820 faceToThrow[1],
821 faceToKeep[1],
822 faceEdges_,
823 ends[1]
824 );
826 meshOps::sizeDownList
827 (
828 faceToThrow[1],
829 edgeFaces_[ends[1]]
830 );
831 }
833 // Correct edgeFaces for triangular faces...
834 forAll(edgeToThrow, indexI)
835 {
836 if (edgeToThrow[indexI] == -1)
837 {
838 continue;
839 }
841 const labelList& eF = edgeFaces_[edgeToThrow[indexI]];
843 label origTriFace = -1, retTriFace = -1;
845 // Find the original triangular face index
846 forAll(eF, faceI)
847 {
848 if (eF[faceI] == c0BdyIndex[0])
849 {
850 origTriFace = c0BdyIndex[0];
851 break;
852 }
854 if (eF[faceI] == c0BdyIndex[1])
855 {
856 origTriFace = c0BdyIndex[1];
857 break;
858 }
860 if (eF[faceI] == c1BdyIndex[0])
861 {
862 origTriFace = c1BdyIndex[0];
863 break;
864 }
866 if (eF[faceI] == c1BdyIndex[1])
867 {
868 origTriFace = c1BdyIndex[1];
869 break;
870 }
871 }
873 // Find the retained triangular face index
874 forAll(eF, faceI)
875 {
876 if
877 (
878 (eF[faceI] != origTriFace) &&
879 (eF[faceI] != faceToThrow[0]) &&
880 (eF[faceI] != faceToThrow[1])
881 )
882 {
883 retTriFace = eF[faceI];
884 break;
885 }
886 }
888 // Finally replace the face index
889 if (retTriFace == -1)
890 {
891 // Couldn't find a retained face.
892 // This must be a boundary edge, so size-down instead.
893 meshOps::sizeDownList
894 (
895 origTriFace,
896 edgeFaces_[edgeToKeep[indexI]]
897 );
898 }
899 else
900 {
901 meshOps::replaceLabel
902 (
903 origTriFace,
904 retTriFace,
905 edgeFaces_[edgeToKeep[indexI]]
906 );
908 meshOps::replaceLabel
909 (
910 edgeToThrow[indexI],
911 edgeToKeep[indexI],
912 faceEdges_[retTriFace]
913 );
914 }
915 }
917 // Correct faceEdges / edgeFaces for quad-faces...
918 forAll(firstEdgeFaces,faceI)
919 {
920 meshOps::replaceLabel
921 (
922 checkEdgeIndex[1],
923 checkEdgeIndex[2],
924 faceEdges_[firstEdgeFaces[faceI]]
925 );
927 // Renumber the edges on this face
928 const labelList& fE = faceEdges_[firstEdgeFaces[faceI]];
930 forAll(fE, edgeI)
931 {
932 if (edges_[fE[edgeI]].start() == cv0)
933 {
934 edges_[fE[edgeI]].start() = cv2;
935 }
937 if (edges_[fE[edgeI]].end() == cv0)
938 {
939 edges_[fE[edgeI]].end() = cv2;
940 }
942 if (edges_[fE[edgeI]].start() == cv1)
943 {
944 edges_[fE[edgeI]].start() = cv3;
945 }
947 if (edges_[fE[edgeI]].end() == cv1)
948 {
949 edges_[fE[edgeI]].end() = cv3;
950 }
951 }
953 // Size-up edgeFaces for the replacement
954 if
955 (
956 (firstEdgeFaces[faceI] != faceToThrow[0]) &&
957 (firstEdgeFaces[faceI] != faceToThrow[1]) &&
958 (firstEdgeFaces[faceI] != fIndex)
959 )
960 {
961 meshOps::sizeUpList
962 (
963 firstEdgeFaces[faceI],
964 edgeFaces_[checkEdgeIndex[2]]
965 );
966 }
967 }
969 // Remove the current face from the replacement edge
970 meshOps::sizeDownList
971 (
972 fIndex,
973 edgeFaces_[checkEdgeIndex[2]]
974 );
976 // Replace point labels on all triangular boundary faces.
977 forAll(hullCells[0],cellI)
978 {
979 const cell& cellToCheck = cells_[hullCells[0][cellI]];
981 forAll(cellToCheck,faceI)
982 {
983 face& faceToCheck = faces_[cellToCheck[faceI]];
985 if (faceToCheck.size() == 3)
986 {
987 forAll(faceToCheck,pointI)
988 {
989 if (faceToCheck[pointI] == cv0)
990 {
991 faceToCheck[pointI] = cv2;
992 }
994 if (faceToCheck[pointI] == cv1)
995 {
996 faceToCheck[pointI] = cv3;
997 }
998 }
999 }
1000 }
1001 }
1003 // Now that we're done with all edges, remove them.
1004 removeEdge(checkEdgeIndex[0]);
1005 removeEdge(checkEdgeIndex[1]);
1006 removeEdge(checkEdgeIndex[3]);
1008 forAll(edgeToThrow, indexI)
1009 {
1010 if (edgeToThrow[indexI] == -1)
1011 {
1012 continue;
1013 }
1015 removeEdge(edgeToThrow[indexI]);
1016 }
1018 // Delete the two points...
1019 removePoint(cv0);
1020 removePoint(cv1);
1021 }
1022 else
1023 {
1024 const labelList& secondEdgeFaces = edgeFaces_[checkEdgeIndex[2]];
1026 // Collapse to the first node
1027 forAll(secondEdgeFaces,faceI)
1028 {
1029 // Replace point indices on faces.
1030 meshOps::replaceLabel
1031 (
1032 cv2,
1033 cv0,
1034 faces_[secondEdgeFaces[faceI]]
1035 );
1037 meshOps::replaceLabel
1038 (
1039 cv3,
1040 cv1,
1041 faces_[secondEdgeFaces[faceI]]
1042 );
1044 // Add an entry for mapping
1045 if (!modifiedFaces.found(secondEdgeFaces[faceI]))
1046 {
1047 modifiedFaces.insert(secondEdgeFaces[faceI]);
1048 }
1050 // Determine the quad-face(s) in cell[0] & cell[1]
1051 // that belongs to secondEdgeFaces
1052 if (secondEdgeFaces[faceI] == c0IntIndex[0])
1053 {
1054 faceToKeep[0] = c0IntIndex[1];
1055 faceToThrow[0] = c0IntIndex[0];
1056 }
1058 if (secondEdgeFaces[faceI] == c0IntIndex[1])
1059 {
1060 faceToKeep[0] = c0IntIndex[0];
1061 faceToThrow[0] = c0IntIndex[1];
1062 }
1064 if (c1 != -1)
1065 {
1066 if (secondEdgeFaces[faceI] == c1IntIndex[0])
1067 {
1068 faceToKeep[1] = c1IntIndex[1];
1069 faceToThrow[1] = c1IntIndex[0];
1070 }
1072 if (secondEdgeFaces[faceI] == c1IntIndex[1])
1073 {
1074 faceToKeep[1] = c1IntIndex[0];
1075 faceToThrow[1] = c1IntIndex[1];
1076 }
1077 }
1078 }
1080 // Find common edges between quad / tri faces...
1081 meshOps::findCommonEdge
1082 (
1083 c0BdyIndex[0],
1084 faceToKeep[0],
1085 faceEdges_,
1086 edgeToKeep[0]
1087 );
1089 meshOps::findCommonEdge
1090 (
1091 c0BdyIndex[1],
1092 faceToKeep[0],
1093 faceEdges_,
1094 edgeToKeep[1]
1095 );
1097 meshOps::findCommonEdge
1098 (
1099 c0BdyIndex[0],
1100 faceToThrow[0],
1101 faceEdges_,
1102 edgeToThrow[0]
1103 );
1105 meshOps::findCommonEdge
1106 (
1107 c0BdyIndex[1],
1108 faceToThrow[0],
1109 faceEdges_,
1110 edgeToThrow[1]
1111 );
1113 // Size down edgeFaces for the ends.
1114 meshOps::findCommonEdge
1115 (
1116 faceToThrow[0],
1117 faceToKeep[0],
1118 faceEdges_,
1119 ends[0]
1120 );
1122 meshOps::sizeDownList
1123 (
1124 faceToThrow[0],
1125 edgeFaces_[ends[0]]
1126 );
1128 if (c1 != -1)
1129 {
1130 meshOps::findCommonEdge
1131 (
1132 c1BdyIndex[0],
1133 faceToKeep[1],
1134 faceEdges_,
1135 edgeToKeep[2]
1136 );
1138 meshOps::findCommonEdge
1139 (
1140 c1BdyIndex[1],
1141 faceToKeep[1],
1142 faceEdges_,
1143 edgeToKeep[3]
1144 );
1146 meshOps::findCommonEdge
1147 (
1148 c1BdyIndex[0],
1149 faceToThrow[1],
1150 faceEdges_,
1151 edgeToThrow[2]
1152 );
1154 meshOps::findCommonEdge
1155 (
1156 c1BdyIndex[1],
1157 faceToThrow[1],
1158 faceEdges_,
1159 edgeToThrow[3]
1160 );
1162 // Size down edgeFaces for the ends.
1163 meshOps::findCommonEdge
1164 (
1165 faceToThrow[1],
1166 faceToKeep[1],
1167 faceEdges_,
1168 ends[1]
1169 );
1171 meshOps::sizeDownList
1172 (
1173 faceToThrow[1],
1174 edgeFaces_[ends[1]]
1175 );
1176 }
1178 // Correct edgeFaces for triangular faces...
1179 forAll(edgeToThrow, indexI)
1180 {
1181 if (edgeToThrow[indexI] == -1)
1182 {
1183 continue;
1184 }
1186 const labelList& eF = edgeFaces_[edgeToThrow[indexI]];
1188 label origTriFace = -1, retTriFace = -1;
1190 // Find the original triangular face index
1191 forAll(eF, faceI)
1192 {
1193 if (eF[faceI] == c0BdyIndex[0])
1194 {
1195 origTriFace = c0BdyIndex[0];
1196 break;
1197 }
1199 if (eF[faceI] == c0BdyIndex[1])
1200 {
1201 origTriFace = c0BdyIndex[1];
1202 break;
1203 }
1205 if (eF[faceI] == c1BdyIndex[0])
1206 {
1207 origTriFace = c1BdyIndex[0];
1208 break;
1209 }
1211 if (eF[faceI] == c1BdyIndex[1])
1212 {
1213 origTriFace = c1BdyIndex[1];
1214 break;
1215 }
1216 }
1218 // Find the retained triangular face index
1219 forAll(eF, faceI)
1220 {
1221 if
1222 (
1223 (eF[faceI] != origTriFace) &&
1224 (eF[faceI] != faceToThrow[0]) &&
1225 (eF[faceI] != faceToThrow[1])
1226 )
1227 {
1228 retTriFace = eF[faceI];
1229 break;
1230 }
1231 }
1233 // Finally replace the face/edge indices
1234 if (retTriFace == -1)
1235 {
1236 // Couldn't find a retained face.
1237 // This must be a boundary edge, so size-down instead.
1238 meshOps::sizeDownList
1239 (
1240 origTriFace,
1241 edgeFaces_[edgeToKeep[indexI]]
1242 );
1243 }
1244 else
1245 {
1246 meshOps::replaceLabel
1247 (
1248 origTriFace,
1249 retTriFace,
1250 edgeFaces_[edgeToKeep[indexI]]
1251 );
1253 meshOps::replaceLabel
1254 (
1255 edgeToThrow[indexI],
1256 edgeToKeep[indexI],
1257 faceEdges_[retTriFace]
1258 );
1259 }
1260 }
1262 // Correct faceEdges / edgeFaces for quad-faces...
1263 forAll(secondEdgeFaces,faceI)
1264 {
1265 meshOps::replaceLabel
1266 (
1267 checkEdgeIndex[2],
1268 checkEdgeIndex[1],
1269 faceEdges_[secondEdgeFaces[faceI]]
1270 );
1272 // Renumber the edges on this face
1273 const labelList& fE = faceEdges_[secondEdgeFaces[faceI]];
1275 forAll(fE, edgeI)
1276 {
1277 if (edges_[fE[edgeI]].start() == cv2)
1278 {
1279 edges_[fE[edgeI]].start() = cv0;
1280 }
1282 if (edges_[fE[edgeI]].end() == cv2)
1283 {
1284 edges_[fE[edgeI]].end() = cv0;
1285 }
1287 if (edges_[fE[edgeI]].start() == cv3)
1288 {
1289 edges_[fE[edgeI]].start() = cv1;
1290 }
1292 if (edges_[fE[edgeI]].end() == cv3)
1293 {
1294 edges_[fE[edgeI]].end() = cv1;
1295 }
1296 }
1298 // Size-up edgeFaces for the replacement
1299 if
1300 (
1301 (secondEdgeFaces[faceI] != faceToThrow[0]) &&
1302 (secondEdgeFaces[faceI] != faceToThrow[1]) &&
1303 (secondEdgeFaces[faceI] != fIndex)
1304 )
1305 {
1306 meshOps::sizeUpList
1307 (
1308 secondEdgeFaces[faceI],
1309 edgeFaces_[checkEdgeIndex[1]]
1310 );
1311 }
1312 }
1314 // Remove the current face from the replacement edge
1315 meshOps::sizeDownList
1316 (
1317 fIndex,
1318 edgeFaces_[checkEdgeIndex[1]]
1319 );
1321 // Replace point labels on all triangular boundary faces.
1322 forAll(hullCells[1], cellI)
1323 {
1324 const cell& cellToCheck = cells_[hullCells[1][cellI]];
1326 forAll(cellToCheck, faceI)
1327 {
1328 face& faceToCheck = faces_[cellToCheck[faceI]];
1330 if (faceToCheck.size() == 3)
1331 {
1332 forAll(faceToCheck, pointI)
1333 {
1334 if (faceToCheck[pointI] == cv2)
1335 {
1336 faceToCheck[pointI] = cv0;
1337 }
1339 if (faceToCheck[pointI] == cv3)
1340 {
1341 faceToCheck[pointI] = cv1;
1342 }
1343 }
1344 }
1345 }
1346 }
1348 // Now that we're done with all edges, remove them.
1349 removeEdge(checkEdgeIndex[0]);
1350 removeEdge(checkEdgeIndex[2]);
1351 removeEdge(checkEdgeIndex[3]);
1353 forAll(edgeToThrow, indexI)
1354 {
1355 if (edgeToThrow[indexI] == -1)
1356 {
1357 continue;
1358 }
1360 removeEdge(edgeToThrow[indexI]);
1361 }
1363 // Delete the two points...
1364 removePoint(cv2);
1365 removePoint(cv3);
1366 }
1368 if (debug > 2)
1369 {
1370 Info << endl;
1371 Info << "~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
1372 Info << "Hulls after modification" << endl;
1373 Info << "~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
1375 Info << nl << "Cells belonging to first Edge Hull: "
1376 << hullCells[0] << endl;
1378 forAll(hullCells[0], cellI)
1379 {
1380 const cell& firstCurCell = cells_[hullCells[0][cellI]];
1382 Info << "Cell: " << hullCells[0][cellI]
1383 << ": " << firstCurCell << endl;
1385 forAll(firstCurCell, faceI)
1386 {
1387 Info << firstCurCell[faceI]
1388 << ": " << faces_[firstCurCell[faceI]] << endl;
1389 }
1390 }
1392 const labelList& firstEdgeFaces = edgeFaces_[checkEdgeIndex[1]];
1394 Info << nl << "First Edge Face Hull: " << firstEdgeFaces << endl;
1396 forAll(firstEdgeFaces, indexI)
1397 {
1398 Info << firstEdgeFaces[indexI]
1399 << ": " << faces_[firstEdgeFaces[indexI]] << endl;
1400 }
1402 Info << nl << "Cells belonging to second Edge Hull: "
1403 << hullCells[1] << endl;
1405 forAll(hullCells[1], cellI)
1406 {
1407 const cell& secondCurCell = cells_[hullCells[1][cellI]];
1409 Info << "Cell: " << hullCells[1][cellI]
1410 << ": " << secondCurCell << endl;
1412 forAll(secondCurCell, faceI)
1413 {
1414 Info << secondCurCell[faceI]
1415 << ": " << faces_[secondCurCell[faceI]] << endl;
1416 }
1417 }
1419 const labelList& secondEdgeFaces = edgeFaces_[checkEdgeIndex[2]];
1421 Info << nl << "Second Edge Face Hull: " << secondEdgeFaces << endl;
1423 forAll(secondEdgeFaces, indexI)
1424 {
1425 Info << secondEdgeFaces[indexI]
1426 << ": " << faces_[secondEdgeFaces[indexI]] << endl;
1427 }
1429 Info << endl;
1431 Info << "Retained face: "
1432 << faceToKeep[0] << ": "
1433 << " owner: " << owner_[faceToKeep[0]]
1434 << " neighbour: " << neighbour_[faceToKeep[0]]
1435 << endl;
1437 Info << "Discarded face: "
1438 << faceToThrow[0] << ": "
1439 << " owner: " << owner_[faceToThrow[0]]
1440 << " neighbour: " << neighbour_[faceToThrow[0]]
1441 << endl;
1443 if (c1 != -1)
1444 {
1445 Info << "Retained face: "
1446 << faceToKeep[1] << ": "
1447 << " owner: " << owner_[faceToKeep[1]]
1448 << " neighbour: " << neighbour_[faceToKeep[1]]
1449 << endl;
1451 Info << "Discarded face: "
1452 << faceToThrow[1] << ": "
1453 << " owner: " << owner_[faceToThrow[1]]
1454 << " neighbour: " << neighbour_[faceToThrow[1]]
1455 << endl;
1456 }
1457 }
1459 // Ensure proper orientation for the two retained faces
1460 FixedList<label,2> cellCheck(0);
1462 if (owner_[faceToThrow[0]] == c0)
1463 {
1464 cellCheck[0] = neighbour_[faceToThrow[0]];
1466 if (owner_[faceToKeep[0]] == c0)
1467 {
1468 if
1469 (
1470 (neighbour_[faceToThrow[0]] > neighbour_[faceToKeep[0]])
1471 && (neighbour_[faceToKeep[0]] != -1)
1472 )
1473 {
1474 // This face is to be flipped
1475 faces_[faceToKeep[0]] =
1476 (
1477 faces_[faceToKeep[0]].reverseFace()
1478 );
1480 owner_[faceToKeep[0]] = neighbour_[faceToKeep[0]];
1481 neighbour_[faceToKeep[0]] = neighbour_[faceToThrow[0]];
1483 setFlip(faceToKeep[0]);
1484 }
1485 else
1486 {
1487 if (neighbour_[faceToThrow[0]] != -1)
1488 {
1489 // Keep orientation intact, and update the owner
1490 owner_[faceToKeep[0]] = neighbour_[faceToThrow[0]];
1491 }
1492 else
1493 {
1494 // This face will need to be flipped and converted
1495 // to a boundary face. Flip it now, so that conversion
1496 // happens later.
1497 faces_[faceToKeep[0]] =
1498 (
1499 faces_[faceToKeep[0]].reverseFace()
1500 );
1502 owner_[faceToKeep[0]] = neighbour_[faceToKeep[0]];
1503 neighbour_[faceToKeep[0]] = -1;
1505 setFlip(faceToKeep[0]);
1506 }
1507 }
1508 }
1509 else
1510 {
1511 // Keep orientation intact, and update the neighbour
1512 neighbour_[faceToKeep[0]] = neighbour_[faceToThrow[0]];
1513 }
1514 }
1515 else
1516 {
1517 cellCheck[0] = owner_[faceToThrow[0]];
1519 if (neighbour_[faceToKeep[0]] == c0)
1520 {
1521 if (owner_[faceToThrow[0]] < owner_[faceToKeep[0]])
1522 {
1523 // This face is to be flipped
1524 faces_[faceToKeep[0]] =
1525 (
1526 faces_[faceToKeep[0]].reverseFace()
1527 );
1529 neighbour_[faceToKeep[0]] = owner_[faceToKeep[0]];
1530 owner_[faceToKeep[0]] = owner_[faceToThrow[0]];
1532 setFlip(faceToKeep[0]);
1533 }
1534 else
1535 {
1536 // Keep orientation intact, and update the neighbour
1537 neighbour_[faceToKeep[0]] = owner_[faceToThrow[0]];
1538 }
1539 }
1540 else
1541 {
1542 // Keep orientation intact, and update the owner
1543 owner_[faceToKeep[0]] = owner_[faceToThrow[0]];
1544 }
1545 }
1547 if (c1 != -1)
1548 {
1549 if (owner_[faceToThrow[1]] == c1)
1550 {
1551 cellCheck[1] = neighbour_[faceToThrow[1]];
1553 if (owner_[faceToKeep[1]] == c1)
1554 {
1555 if
1556 (
1557 (neighbour_[faceToThrow[1]] > neighbour_[faceToKeep[1]])
1558 && (neighbour_[faceToKeep[1]] != -1)
1559 )
1560 {
1561 // This face is to be flipped
1562 faces_[faceToKeep[1]] =
1563 (
1564 faces_[faceToKeep[1]].reverseFace()
1565 );
1567 owner_[faceToKeep[1]] = neighbour_[faceToKeep[1]];
1568 neighbour_[faceToKeep[1]] = neighbour_[faceToThrow[1]];
1570 setFlip(faceToKeep[1]);
1571 }
1572 else
1573 {
1574 if (neighbour_[faceToThrow[1]] != -1)
1575 {
1576 // Keep orientation intact, and update the owner
1577 owner_[faceToKeep[1]] = neighbour_[faceToThrow[1]];
1578 }
1579 else
1580 {
1581 // This face will need to be flipped and converted
1582 // to a boundary face. Flip it now, so that conversion
1583 // happens later.
1584 faces_[faceToKeep[1]] =
1585 (
1586 faces_[faceToKeep[1]].reverseFace()
1587 );
1589 owner_[faceToKeep[1]] = neighbour_[faceToKeep[1]];
1590 neighbour_[faceToKeep[1]] = -1;
1592 setFlip(faceToKeep[1]);
1593 }
1594 }
1595 }
1596 else
1597 {
1598 // Keep orientation intact, and update the neighbour
1599 neighbour_[faceToKeep[1]] = neighbour_[faceToThrow[1]];
1600 }
1601 }
1602 else
1603 {
1604 cellCheck[1] = owner_[faceToThrow[1]];
1606 if (neighbour_[faceToKeep[1]] == c1)
1607 {
1608 if (owner_[faceToThrow[1]] < owner_[faceToKeep[1]])
1609 {
1610 // This face is to be flipped
1611 faces_[faceToKeep[1]] =
1612 (
1613 faces_[faceToKeep[1]].reverseFace()
1614 );
1616 neighbour_[faceToKeep[1]] = owner_[faceToKeep[1]];
1617 owner_[faceToKeep[1]] = owner_[faceToThrow[1]];
1619 setFlip(faceToKeep[1]);
1620 }
1621 else
1622 {
1623 // Keep orientation intact, and update the neighbour
1624 neighbour_[faceToKeep[1]] = owner_[faceToThrow[1]];
1625 }
1626 }
1627 else
1628 {
1629 // Keep orientation intact, and update the owner
1630 owner_[faceToKeep[1]] = owner_[faceToThrow[1]];
1631 }
1632 }
1633 }
1635 // Remove orphaned faces
1636 if (owner_[faceToKeep[0]] == -1)
1637 {
1638 removeFace(faceToKeep[0]);
1639 }
1640 else
1641 if
1642 (
1643 (neighbour_[faceToKeep[0]] == -1)
1644 && (whichPatch(faceToKeep[0]) < 0)
1645 )
1646 {
1647 // Obtain a copy before adding the new face,
1648 // since the reference might become invalid during list resizing.
1649 face newFace = faces_[faceToKeep[0]];
1650 label newOwn = owner_[faceToKeep[0]];
1651 labelList newFaceEdges = faceEdges_[faceToKeep[0]];
1653 // This face is being converted from interior to boundary. Remove
1654 // from the interior list and add as a boundary face to the end.
1655 label newFaceIndex =
1656 (
1657 insertFace
1658 (
1659 whichPatch(faceToThrow[0]),
1660 newFace,
1661 newOwn,
1662 -1
1663 )
1664 );
1666 // Add an entry for mapping
1667 if (!modifiedFaces.found(newFaceIndex))
1668 {
1669 modifiedFaces.insert(newFaceIndex);
1670 }
1672 // Add a faceEdges entry as well.
1673 // Edges don't have to change, since they're
1674 // all on the boundary anyway.
1675 faceEdges_.append(newFaceEdges);
1677 meshOps::replaceLabel
1678 (
1679 faceToKeep[0],
1680 newFaceIndex,
1681 cells_[newOwn]
1682 );
1684 // Correct edgeFaces with the new face label.
1685 forAll(newFaceEdges, edgeI)
1686 {
1687 meshOps::replaceLabel
1688 (
1689 faceToKeep[0],
1690 newFaceIndex,
1691 edgeFaces_[newFaceEdges[edgeI]]
1692 );
1693 }
1695 // Renumber the neighbour so that this face is removed correctly.
1696 neighbour_[faceToKeep[0]] = 0;
1697 removeFace(faceToKeep[0]);
1698 }
1700 // Remove the unwanted faces in the cell(s) adjacent to this face,
1701 // and correct the cells that contain discarded faces
1702 const cell& cell_0 = cells_[c0];
1704 forAll(cell_0,faceI)
1705 {
1706 if (cell_0[faceI] != fIndex && cell_0[faceI] != faceToKeep[0])
1707 {
1708 removeFace(cell_0[faceI]);
1709 }
1710 }
1712 if (cellCheck[0] != -1)
1713 {
1714 meshOps::replaceLabel
1715 (
1716 faceToThrow[0],
1717 faceToKeep[0],
1718 cells_[cellCheck[0]]
1719 );
1720 }
1722 // Remove the cell
1723 removeCell(c0);
1725 if (c1 == -1)
1726 {
1727 // Increment the surface-collapse counter
1728 statistics_[6]++;
1729 }
1730 else
1731 {
1732 // Remove orphaned faces
1733 if (owner_[faceToKeep[1]] == -1)
1734 {
1735 removeFace(faceToKeep[1]);
1736 }
1737 else
1738 if
1739 (
1740 (neighbour_[faceToKeep[1]] == -1)
1741 && (whichPatch(faceToKeep[1]) < 0)
1742 )
1743 {
1744 // Obtain a copy before adding the new face,
1745 // since the reference might become invalid during list resizing.
1746 face newFace = faces_[faceToKeep[1]];
1747 label newOwn = owner_[faceToKeep[1]];
1748 labelList newFaceEdges = faceEdges_[faceToKeep[1]];
1750 // This face is being converted from interior to boundary. Remove
1751 // from the interior list and add as a boundary face to the end.
1752 label newFaceIndex =
1753 (
1754 insertFace
1755 (
1756 whichPatch(faceToThrow[1]),
1757 newFace,
1758 newOwn,
1759 -1
1760 )
1761 );
1763 // Add an entry for mapping
1764 if (!modifiedFaces.found(newFaceIndex))
1765 {
1766 modifiedFaces.insert(newFaceIndex);
1767 }
1769 // Add a faceEdges entry as well.
1770 // Edges don't have to change, since they're
1771 // all on the boundary anyway.
1772 faceEdges_.append(newFaceEdges);
1774 meshOps::replaceLabel
1775 (
1776 faceToKeep[1],
1777 newFaceIndex,
1778 cells_[newOwn]
1779 );
1781 // Correct edgeFaces with the new face label.
1782 forAll(newFaceEdges, edgeI)
1783 {
1784 meshOps::replaceLabel
1785 (
1786 faceToKeep[1],
1787 newFaceIndex,
1788 edgeFaces_[newFaceEdges[edgeI]]
1789 );
1790 }
1792 // Renumber the neighbour so that this face is removed correctly.
1793 neighbour_[faceToKeep[1]] = 0;
1794 removeFace(faceToKeep[1]);
1795 }
1797 const cell& cell_1 = cells_[c1];
1799 forAll(cell_1, faceI)
1800 {
1801 if (cell_1[faceI] != fIndex && cell_1[faceI] != faceToKeep[1])
1802 {
1803 removeFace(cell_1[faceI]);
1804 }
1805 }
1807 if (cellCheck[1] != -1)
1808 {
1809 meshOps::replaceLabel
1810 (
1811 faceToThrow[1],
1812 faceToKeep[1],
1813 cells_[cellCheck[1]]
1814 );
1815 }
1817 // Remove the cell
1818 removeCell(c1);
1819 }
1821 // Move old / new points
1822 oldPoints_[replacement[0]] = oldPoint[0];
1823 oldPoints_[replacement[1]] = oldPoint[1];
1825 points_[replacement[0]] = newPoint[0];
1826 points_[replacement[1]] = newPoint[1];
1828 // Finally remove the face
1829 removeFace(fIndex);
1831 // Write out VTK files after change
1832 if (debug > 3)
1833 {
1834 labelHashSet vtkCells;
1836 forAll(hullCells[0], cellI)
1837 {
1838 if (hullCells[0][cellI] == c0 || hullCells[0][cellI] == c1)
1839 {
1840 continue;
1841 }
1843 if (!vtkCells.found(hullCells[0][cellI]))
1844 {
1845 vtkCells.insert(hullCells[0][cellI]);
1846 }
1847 }
1849 forAll(hullCells[1], cellI)
1850 {
1851 if (hullCells[1][cellI] == c0 || hullCells[1][cellI] == c1)
1852 {
1853 continue;
1854 }
1856 if (!vtkCells.found(hullCells[1][cellI]))
1857 {
1858 vtkCells.insert(hullCells[1][cellI]);
1859 }
1860 }
1862 writeVTK
1863 (
1864 Foam::name(fIndex)
1865 + "_Collapse_1",
1866 vtkCells.toc()
1867 );
1868 }
1870 // Specify that an old point-position
1871 // has been modified, if necessary
1872 if (collapseCase == 3 && c1 > -1)
1873 {
1874 labelList mP(2, -1);
1876 mP[0] = original[0];
1877 mP[1] = replacement[0];
1879 modPoints_.set(replacement[0], mP);
1881 mP[0] = original[1];
1882 mP[1] = replacement[1];
1884 modPoints_.set(replacement[1], mP);
1885 }
1887 // Fill-in candidate mapping information
1888 labelList mC(2, -1);
1889 mC[0] = c0, mC[1] = c1;
1891 // Now that all old / new cells possess correct connectivity,
1892 // compute mapping information.
1893 forAll(hullCells, indexI)
1894 {
1895 forAll(hullCells[indexI], cellI)
1896 {
1897 label mcIndex = hullCells[indexI][cellI];
1899 // Skip collapsed cells
1900 if (mcIndex == c0 || mcIndex == c1)
1901 {
1902 continue;
1903 }
1905 // Set the mapping for this cell
1906 setCellMapping(mcIndex, mC);
1907 }
1908 }
1910 // Set face mapping information for modified faces
1911 forAllConstIter(labelHashSet, modifiedFaces, fIter)
1912 {
1913 // Exclude deleted faces
1914 if (faces_[fIter.key()].empty())
1915 {
1916 continue;
1917 }
1919 // Decide between default / weighted mapping
1920 // based on boundary information
1921 label fPatch = whichPatch(fIter.key());
1923 if (fPatch == -1)
1924 {
1925 setFaceMapping(fIter.key());
1926 }
1927 else
1928 {
1929 // Fill-in candidate mapping information
1930 labelList faceCandidates;
1932 const labelList& fEdges = faceEdges_[fIter.key()];
1934 forAll(fEdges, edgeI)
1935 {
1936 if (whichEdgePatch(fEdges[edgeI]) == fPatch)
1937 {
1938 // Loop through associated edgeFaces
1939 const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
1941 forAll(eFaces, faceI)
1942 {
1943 if
1944 (
1945 (eFaces[faceI] != fIter.key()) &&
1946 (whichPatch(eFaces[faceI]) == fPatch)
1947 )
1948 {
1949 faceCandidates.setSize
1950 (
1951 faceCandidates.size() + 1,
1952 eFaces[faceI]
1953 );
1954 }
1955 }
1956 }
1957 }
1959 // Set the mapping for this face
1960 setFaceMapping(fIter.key(), faceCandidates);
1961 }
1962 }
1964 // Set the flag
1965 topoChangeFlag_ = true;
1967 // Increment the counter
1968 statistics_[4]++;
1970 // Increment the number of modifications
1971 statistics_[0]++;
1973 // Return a succesful collapse
1974 map.type() = collapseCase;
1976 return map;
1977 }
1980 // Method for the collapse of an edge in 3D
1981 // - Returns a changeMap with a type specifying:
1982 // -1: Collapse failed since max number of topo-changes was reached.
1983 // 0: Collapse could not be performed.
1984 // 1: Collapsed to first node.
1985 // 2: Collapsed to second node.
1986 // 3: Collapsed to mid-point (default)
1987 // - overRideCase is used to force a certain collapse configuration.
1988 // -1: Use this value to let collapseEdge decide a case.
1989 // 1: Force collapse to first node.
1990 // 2: Force collapse to second node.
1991 // 3: Force collapse to mid-point
1992 // - checkOnly performs a feasibility check and returns without modifications.
1993 // - forceOp to force the collapse.
1994 const changeMap dynamicTopoFvMesh::collapseEdge
1995 (
1996 const label eIndex,
1997 label overRideCase,
1998 bool checkOnly,
1999 bool forceOp
2000 )
2001 {
2002 // Edge collapse performs the following operations:
2003 // [1] Checks if either vertex of the edge is on a boundary
2004 // [2] Checks whether cells attached to deleted vertices will be valid
2005 // after the edge-collapse operation
2006 // [3] Deletes all cells surrounding this edge
2007 // [4] Deletes all faces surrounding this edge
2008 // [5] Deletes all faces surrounding the deleted vertex attached
2009 // to the cells in [3]
2010 // [6] Checks the orientation of faces connected to the retained
2011 // vertices
2012 // [7] Remove one of the vertices of the edge
2013 // Update faceEdges, edgeFaces and edgePoints information
2015 // For 2D meshes, perform face-collapse
2016 if (twoDMesh_)
2017 {
2018 return collapseQuadFace(eIndex, overRideCase, checkOnly);
2019 }
2021 // Figure out which thread this is...
2022 label tIndex = self();
2024 // Prepare the changeMaps
2025 changeMap map, slaveMap;
2027 if
2028 (
2029 (statistics_[0] > maxModifications_)
2030 && (maxModifications_ > -1)
2031 )
2032 {
2033 // Reached the max allowable topo-changes.
2034 stack(tIndex).clear();
2036 return map;
2037 }
2039 // Check if edgeRefinements are to be avoided on patch.
2040 if (lengthEstimator().checkRefinementPatch(whichEdgePatch(eIndex)))
2041 {
2042 return map;
2043 }
2045 // Sanity check: Is the index legitimate?
2046 if (eIndex < 0)
2047 {
2048 FatalErrorIn
2049 (
2050 "\n"
2051 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
2052 "(\n"
2053 " const label eIndex,\n"
2054 " label overRideCase,\n"
2055 " bool checkOnly,\n"
2056 " bool forceOp\n"
2057 ")\n"
2058 )
2059 << " Invalid index: " << eIndex
2060 << abort(FatalError);
2061 }
2063 // Hull variables
2064 bool found = false;
2065 label replaceIndex = -1, m = edgePoints_[eIndex].size();
2067 // Size up the hull lists
2068 labelList cellHull(m, -1);
2069 labelList faceHull(m, -1);
2070 labelList edgeHull(m, -1);
2071 labelListList ringEntities(4, labelList(m, -1));
2073 // Construct a hull around this edge
2074 meshOps::constructHull
2075 (
2076 eIndex,
2077 faces_,
2078 edges_,
2079 cells_,
2080 owner_,
2081 neighbour_,
2082 faceEdges_,
2083 edgeFaces_,
2084 edgePoints_,
2085 edgeHull,
2086 faceHull,
2087 cellHull,
2088 ringEntities
2089 );
2091 // Check whether points of the edge lies on a boundary
2092 FixedList<label, 2> nBoundCurves(0), checkPoints(-1);
2093 const FixedList<bool,2> edgeBoundary = checkEdgeBoundary(eIndex);
2095 // Decide on collapseCase
2096 label collapseCase = -1;
2098 if (edgeBoundary[0] && !edgeBoundary[1])
2099 {
2100 collapseCase = 1;
2101 }
2102 else
2103 if (!edgeBoundary[0] && edgeBoundary[1])
2104 {
2105 collapseCase = 2;
2106 }
2107 else
2108 if (edgeBoundary[0] && edgeBoundary[1])
2109 {
2110 // If this is an interior edge with two boundary points.
2111 // Bail out for now. If proximity based refinement is
2112 // switched on, mesh may be sliced at this point.
2113 if (whichEdgePatch(eIndex) == -1)
2114 {
2115 return map;
2116 }
2118 // Check if either point lies on a bounding curve
2119 // Used to ensure that collapses happen towards bounding curves.
2120 // If the edge itself is on a bounding curve, collapse is valid.
2121 forAll(edges_[eIndex], pointI)
2122 {
2123 const labelList& pEdges = pointEdges_[edges_[eIndex][pointI]];
2125 forAll(pEdges, edgeI)
2126 {
2127 if (checkBoundingCurve(pEdges[edgeI]))
2128 {
2129 nBoundCurves[pointI]++;
2130 }
2131 }
2132 }
2134 // Pick the point which is connected to more bounding curves
2135 if (nBoundCurves[0] > nBoundCurves[1])
2136 {
2137 collapseCase = 1;
2138 }
2139 else
2140 if (nBoundCurves[1] > nBoundCurves[0])
2141 {
2142 collapseCase = 2;
2143 }
2144 else
2145 {
2146 // Bounding edge: collapseEdge can collapse this edge
2147 collapseCase = 3;
2148 }
2149 }
2150 else
2151 {
2152 // Looks like this is an interior edge.
2153 // Collapse case [3] by default
2154 collapseCase = 3;
2155 }
2157 // Perform an override if requested.
2158 if (overRideCase != -1)
2159 {
2160 collapseCase = overRideCase;
2161 }
2163 // Configure the new point-position
2164 point newPoint = vector::zero;
2165 point oldPoint = vector::zero;
2167 label collapsePoint = -1, replacePoint = -1;
2169 switch (collapseCase)
2170 {
2171 case 1:
2172 {
2173 // Collapse to the first node
2174 replacePoint = edges_[eIndex][0];
2175 collapsePoint = edges_[eIndex][1];
2177 newPoint = points_[replacePoint];
2178 oldPoint = oldPoints_[replacePoint];
2180 checkPoints[0] = collapsePoint;
2182 break;
2183 }
2185 case 2:
2186 {
2187 // Collapse to the second node
2188 replacePoint = edges_[eIndex][1];
2189 collapsePoint = edges_[eIndex][0];
2191 newPoint = points_[replacePoint];
2192 oldPoint = oldPoints_[replacePoint];
2194 checkPoints[0] = collapsePoint;
2196 break;
2197 }
2199 case 3:
2200 {
2201 // Collapse to the mid-point
2202 replacePoint = edges_[eIndex][1];
2203 collapsePoint = edges_[eIndex][0];
2205 newPoint =
2206 (
2207 0.5 *
2208 (
2209 points_[replacePoint]
2210 + points_[collapsePoint]
2211 )
2212 );
2214 oldPoint =
2215 (
2216 0.5 *
2217 (
2218 oldPoints_[replacePoint]
2219 + oldPoints_[collapsePoint]
2220 )
2221 );
2223 checkPoints[0] = replacePoint;
2224 checkPoints[1] = collapsePoint;
2226 break;
2227 }
2229 default:
2230 {
2231 // Don't think this will ever happen.
2232 FatalErrorIn
2233 (
2234 "\n"
2235 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
2236 "(\n"
2237 " const label eIndex,\n"
2238 " label overRideCase,\n"
2239 " bool checkOnly,\n"
2240 " bool forceOp\n"
2241 ")\n"
2242 )
2243 << "Edge: " << eIndex << ": " << edges_[eIndex]
2244 << ". Couldn't decide on collapseCase."
2245 << abort(FatalError);
2247 break;
2248 }
2249 }
2251 // Loop through edges and check for feasibility of collapse
2252 // Also, keep track of resulting cell quality,
2253 // if collapse is indeed feasible
2254 scalar collapseQuality(GREAT);
2255 labelHashSet cellsChecked;
2257 // Add all hull cells as 'checked',
2258 // and therefore, feasible
2259 forAll(cellHull, cellI)
2260 {
2261 if (cellHull[cellI] == -1)
2262 {
2263 continue;
2264 }
2266 cellsChecked.insert(cellHull[cellI]);
2267 }
2269 // Check collapsibility of cells around edges
2270 // with the re-configured point
2271 forAll(checkPoints, pointI)
2272 {
2273 if (checkPoints[pointI] == -1)
2274 {
2275 continue;
2276 }
2278 const labelList& checkPointEdges = pointEdges_[checkPoints[pointI]];
2280 forAll(checkPointEdges, edgeI)
2281 {
2282 const labelList& eFaces = edgeFaces_[checkPointEdges[edgeI]];
2284 // Build a list of cells to check
2285 forAll(eFaces, faceI)
2286 {
2287 label own = owner_[eFaces[faceI]];
2288 label nei = neighbour_[eFaces[faceI]];
2290 // Check owner cell
2291 if (!cellsChecked.found(own))
2292 {
2293 // Check if a collapse is feasible
2294 if
2295 (
2296 checkCollapse
2297 (
2298 newPoint,
2299 oldPoint,
2300 checkPoints[pointI],
2301 own,
2302 cellsChecked,
2303 collapseQuality,
2304 forceOp
2305 )
2306 )
2307 {
2308 map.type() = 0;
2309 return map;
2310 }
2311 }
2313 // Check neighbour cell
2314 if (!cellsChecked.found(nei) && nei != -1)
2315 {
2316 // Check if a collapse is feasible
2317 if
2318 (
2319 checkCollapse
2320 (
2321 newPoint,
2322 oldPoint,
2323 checkPoints[pointI],
2324 nei,
2325 cellsChecked,
2326 collapseQuality,
2327 forceOp
2328 )
2329 )
2330 {
2331 map.type() = 0;
2332 return map;
2333 }
2334 }
2335 }
2336 }
2337 }
2339 // Are we only performing checks?
2340 if (checkOnly)
2341 {
2342 map.type() = collapseCase;
2344 if (debug > 2)
2345 {
2346 Info << "Edge: " << eIndex
2347 << ":: " << edges_[eIndex] << nl
2348 << " collapseCase determined to be: "
2349 << collapseCase << nl
2350 << " Resulting quality: "
2351 << collapseQuality
2352 << endl;
2353 }
2355 return map;
2356 }
2358 // Update number of surface collapses, if necessary.
2359 if (whichEdgePatch(eIndex) > -1)
2360 {
2361 statistics_[6]++;
2362 }
2364 // Define indices on the hull for removal / replacement
2365 label removeEdgeIndex = -1, replaceEdgeIndex = -1;
2366 label removeFaceIndex = -1, replaceFaceIndex = -1;
2368 if (replacePoint == edges_[eIndex][0])
2369 {
2370 replaceEdgeIndex = 0;
2371 replaceFaceIndex = 1;
2372 removeEdgeIndex = 2;
2373 removeFaceIndex = 3;
2374 }
2375 else
2376 if (replacePoint == edges_[eIndex][1])
2377 {
2378 removeEdgeIndex = 0;
2379 removeFaceIndex = 1;
2380 replaceEdgeIndex = 2;
2381 replaceFaceIndex = 3;
2382 }
2383 else
2384 {
2385 // Don't think this will ever happen.
2386 FatalErrorIn
2387 (
2388 "\n"
2389 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
2390 "(\n"
2391 " const label eIndex,\n"
2392 " label overRideCase,\n"
2393 " bool checkOnly,\n"
2394 " bool forceOp\n"
2395 ")\n"
2396 )
2397 << "Edge: " << eIndex << ": " << edges_[eIndex]
2398 << ". Couldn't decide on removal / replacement indices."
2399 << abort(FatalError);
2400 }
2402 if (debug > 1)
2403 {
2404 Info << nl << nl
2405 << "Edge: " << eIndex
2406 << ": " << edges_[eIndex]
2407 << " is to be collapsed. " << endl;
2409 label epIndex = whichEdgePatch(eIndex);
2411 Info << "Patch: ";
2413 if (epIndex == -1)
2414 {
2415 Info << "Internal" << endl;
2416 }
2417 else
2418 {
2419 Info << boundaryMesh()[epIndex].name() << endl;
2420 }
2422 Info << " nBoundCurves: " << nBoundCurves << endl;
2423 Info << " collapseCase: " << collapseCase << endl;
2425 Info << " Resulting quality: " << collapseQuality << endl;
2427 if (debug > 2)
2428 {
2429 Info << "Vertices: " << edgePoints_[eIndex] << endl;
2430 Info << "Edges: " << edgeHull << endl;
2431 Info << "Faces: " << faceHull << endl;
2432 Info << "Cells: " << cellHull << endl;
2433 Info << "replacePoint: " << replacePoint << endl;
2434 Info << "collapsePoint: " << collapsePoint << endl;
2435 Info << "checkPoints: " << checkPoints << endl;;
2436 Info << "ringEntities (removed faces): " << endl;
2438 forAll(ringEntities[removeFaceIndex], faceI)
2439 {
2440 label fIndex = ringEntities[removeFaceIndex][faceI];
2442 if (fIndex == -1)
2443 {
2444 continue;
2445 }
2447 Info << fIndex << ": " << faces_[fIndex] << endl;
2448 }
2450 Info << "ringEntities (removed edges): " << endl;
2451 forAll(ringEntities[removeEdgeIndex], edgeI)
2452 {
2453 label ieIndex = ringEntities[removeEdgeIndex][edgeI];
2455 if (ieIndex == -1)
2456 {
2457 continue;
2458 }
2460 Info << ieIndex << ": " << edges_[ieIndex] << endl;
2461 }
2463 Info << "ringEntities (replacement faces): " << endl;
2464 forAll(ringEntities[replaceFaceIndex], faceI)
2465 {
2466 label fIndex = ringEntities[replaceFaceIndex][faceI];
2468 if (fIndex == -1)
2469 {
2470 continue;
2471 }
2473 Info << fIndex << ": " << faces_[fIndex] << endl;
2474 }
2476 Info << "ringEntities (replacement edges): " << endl;
2477 forAll(ringEntities[replaceEdgeIndex], edgeI)
2478 {
2479 label ieIndex = ringEntities[replaceEdgeIndex][edgeI];
2481 if (ieIndex == -1)
2482 {
2483 continue;
2484 }
2486 Info << ieIndex << ": " << edges_[ieIndex] << endl;
2487 }
2489 labelList& collapsePointEdges = pointEdges_[collapsePoint];
2491 Info << "pointEdges (collapsePoint): ";
2493 forAll(collapsePointEdges, edgeI)
2494 {
2495 Info << collapsePointEdges[edgeI] << " ";
2496 }
2498 Info << endl;
2500 // Write out VTK files prior to change
2501 if (debug > 3)
2502 {
2503 labelList vtkCells = cellsChecked.toc();
2505 writeVTK
2506 (
2507 Foam::name(eIndex)
2508 + '(' + Foam::name(edges_[eIndex][0])
2509 + ',' + Foam::name(edges_[eIndex][1]) + ')'
2510 + "_Collapse_0",
2511 vtkCells
2512 );
2513 }
2514 }
2515 }
2517 // Maintain a list of modified faces for mapping
2518 labelHashSet modifiedFaces;
2520 // Renumber all hull faces and edges
2521 forAll(faceHull, indexI)
2522 {
2523 // Loop through all faces of the edge to be removed
2524 // and reassign them to the replacement edge
2525 label edgeToRemove = ringEntities[removeEdgeIndex][indexI];
2526 label faceToRemove = ringEntities[removeFaceIndex][indexI];
2527 label cellToRemove = cellHull[indexI];
2528 label replaceEdge = ringEntities[replaceEdgeIndex][indexI];
2529 label replaceFace = ringEntities[replaceFaceIndex][indexI];
2531 const labelList& rmvEdgeFaces = edgeFaces_[edgeToRemove];
2533 // Replace edgePoints for all edges emanating from hullVertices
2534 // except ring-edges; those are sized-down later
2535 const labelList& hullPointEdges =
2536 (
2537 pointEdges_[edgePoints_[eIndex][indexI]]
2538 );
2540 forAll(hullPointEdges, edgeI)
2541 {
2542 if
2543 (
2544 findIndex
2545 (
2546 edgePoints_[hullPointEdges[edgeI]],
2547 collapsePoint
2548 ) != -1
2549 && findIndex
2550 (
2551 edgePoints_[hullPointEdges[edgeI]],
2552 replacePoint
2553 ) == -1
2554 )
2555 {
2556 meshOps::replaceLabel
2557 (
2558 collapsePoint,
2559 replacePoint,
2560 edgePoints_[hullPointEdges[edgeI]]
2561 );
2562 }
2563 }
2565 forAll(rmvEdgeFaces, faceI)
2566 {
2567 // Replace edge labels for faces
2568 meshOps::replaceLabel
2569 (
2570 edgeToRemove,
2571 replaceEdge,
2572 faceEdges_[rmvEdgeFaces[faceI]]
2573 );
2575 // Loop through faces associated with this edge,
2576 // and renumber them as well.
2577 const face& faceToCheck = faces_[rmvEdgeFaces[faceI]];
2579 if ((replaceIndex = faceToCheck.which(collapsePoint)) > -1)
2580 {
2581 if (debug > 2)
2582 {
2583 Info << "Renumbering face: "
2584 << rmvEdgeFaces[faceI] << ": "
2585 << faceToCheck << endl;
2586 }
2588 // Renumber the face...
2589 faces_[rmvEdgeFaces[faceI]][replaceIndex] = replacePoint;
2591 // Add an entry for mapping
2592 if (!modifiedFaces.found(rmvEdgeFaces[faceI]))
2593 {
2594 modifiedFaces.insert(rmvEdgeFaces[faceI]);
2595 }
2596 }
2598 // Hull faces should be removed for the replacement edge
2599 if (rmvEdgeFaces[faceI] == faceHull[indexI])
2600 {
2601 meshOps::sizeDownList
2602 (
2603 faceHull[indexI],
2604 edgeFaces_[replaceEdge]
2605 );
2607 continue;
2608 }
2610 found = false;
2612 // Need to avoid ring faces as well.
2613 forAll(ringEntities[removeFaceIndex], faceII)
2614 {
2615 if
2616 (
2617 rmvEdgeFaces[faceI]
2618 == ringEntities[removeFaceIndex][faceII]
2619 )
2620 {
2621 found = true;
2622 break;
2623 }
2624 }
2626 // Size-up the replacement edge list if the face hasn't been found.
2627 // These faces are connected to the edge slated for
2628 // removal, but do not belong to the hull.
2629 if (!found)
2630 {
2631 meshOps::sizeUpList
2632 (
2633 rmvEdgeFaces[faceI],
2634 edgeFaces_[replaceEdge]
2635 );
2636 }
2637 }
2639 if (cellToRemove == -1)
2640 {
2641 continue;
2642 }
2644 // Size down edgeFaces for the ring edges
2645 meshOps::sizeDownList
2646 (
2647 faceToRemove,
2648 edgeFaces_[edgeHull[indexI]]
2649 );
2651 // Size down edgePoints for the ring edges
2652 meshOps::sizeDownList
2653 (
2654 collapsePoint,
2655 edgePoints_[edgeHull[indexI]]
2656 );
2658 // Ensure proper orientation of retained faces
2659 if (owner_[faceToRemove] == cellToRemove)
2660 {
2661 if (owner_[replaceFace] == cellToRemove)
2662 {
2663 if
2664 (
2665 (neighbour_[faceToRemove] > neighbour_[replaceFace])
2666 && (neighbour_[replaceFace] != -1)
2667 )
2668 {
2669 // This face is to be flipped
2670 faces_[replaceFace] = faces_[replaceFace].reverseFace();
2671 owner_[replaceFace] = neighbour_[replaceFace];
2672 neighbour_[replaceFace] = neighbour_[faceToRemove];
2674 setFlip(replaceFace);
2675 }
2676 else
2677 if
2678 (
2679 (neighbour_[faceToRemove] == -1) &&
2680 (neighbour_[replaceFace] != -1)
2681 )
2682 {
2683 // This interior face would need to be converted
2684 // to a boundary one, and flipped as well.
2685 face newFace = faces_[replaceFace].reverseFace();
2686 label newOwner = neighbour_[replaceFace];
2687 label newNeighbour = neighbour_[faceToRemove];
2688 labelList newFE = faceEdges_[replaceFace];
2690 label newFaceIndex =
2691 (
2692 insertFace
2693 (
2694 whichPatch(faceToRemove),
2695 newFace,
2696 newOwner,
2697 newNeighbour
2698 )
2699 );
2701 // Set this face aside for mapping
2702 if (!modifiedFaces.found(newFaceIndex))
2703 {
2704 modifiedFaces.insert(newFaceIndex);
2705 }
2707 // Ensure that all edges of this face are
2708 // on the boundary.
2709 forAll(newFE, edgeI)
2710 {
2711 if (whichEdgePatch(newFE[edgeI]) == -1)
2712 {
2713 edge newEdge = edges_[newFE[edgeI]];
2714 labelList newEF = edgeFaces_[newFE[edgeI]];
2715 labelList newEP = edgePoints_[newFE[edgeI]];
2717 // Need patch information for the new edge.
2718 // Find the corresponding edge in ringEntities.
2719 // Note that hullEdges doesn't need to be checked,
2720 // since they are common to both faces.
2721 label i =
2722 (
2723 findIndex
2724 (
2725 ringEntities[replaceEdgeIndex],
2726 newFE[edgeI]
2727 )
2728 );
2730 label repIndex =
2731 (
2732 whichEdgePatch
2733 (
2734 ringEntities[removeEdgeIndex][i]
2735 )
2736 );
2738 // Insert the new edge
2739 label newEdgeIndex =
2740 (
2741 insertEdge
2742 (
2743 repIndex,
2744 newEdge,
2745 newEF,
2746 newEP
2747 )
2748 );
2750 // Replace faceEdges for all
2751 // connected faces.
2752 forAll(newEF, faceI)
2753 {
2754 meshOps::replaceLabel
2755 (
2756 newFE[edgeI],
2757 newEdgeIndex,
2758 faceEdges_[newEF[faceI]]
2759 );
2760 }
2762 // Remove the edge
2763 removeEdge(newFE[edgeI]);
2765 // Replace faceEdges with new edge index
2766 newFE[edgeI] = newEdgeIndex;
2768 // Modify ringEntities
2769 ringEntities[replaceEdgeIndex][i] = newEdgeIndex;
2770 }
2771 }
2773 // Add the new faceEdges
2774 faceEdges_.append(newFE);
2776 // Replace edgeFaces with the new face index
2777 const labelList& newFEdges = faceEdges_[newFaceIndex];
2779 forAll(newFEdges, edgeI)
2780 {
2781 meshOps::replaceLabel
2782 (
2783 replaceFace,
2784 newFaceIndex,
2785 edgeFaces_[newFEdges[edgeI]]
2786 );
2787 }
2789 // Remove the face
2790 removeFace(replaceFace);
2792 // Replace label for the new owner
2793 meshOps::replaceLabel
2794 (
2795 replaceFace,
2796 newFaceIndex,
2797 cells_[newOwner]
2798 );
2800 // Modify ringEntities and replaceFace
2801 replaceFace = newFaceIndex;
2802 ringEntities[replaceFaceIndex][indexI] = newFaceIndex;
2803 }
2804 else
2805 if
2806 (
2807 (neighbour_[faceToRemove] == -1) &&
2808 (neighbour_[replaceFace] == -1)
2809 )
2810 {
2811 // Wierd overhanging cell. Since replaceFace
2812 // would be an orphan if this continued, remove
2813 // the face entirely.
2814 labelList rmFE = faceEdges_[replaceFace];
2816 forAll(rmFE, edgeI)
2817 {
2818 if
2819 (
2820 (edgeFaces_[rmFE[edgeI]].size() == 1) &&
2821 (edgeFaces_[rmFE[edgeI]][0] == replaceFace)
2822 )
2823 {
2824 // This edge has to be removed entirely.
2825 removeEdge(rmFE[edgeI]);
2827 label i =
2828 (
2829 findIndex
2830 (
2831 ringEntities[replaceEdgeIndex],
2832 rmFE[edgeI]
2833 )
2834 );
2836 // Modify ringEntities
2837 ringEntities[replaceEdgeIndex][i] = -1;
2838 }
2839 else
2840 {
2841 // Size-down edgeFaces
2842 meshOps::sizeDownList
2843 (
2844 replaceFace,
2845 edgeFaces_[rmFE[edgeI]]
2846 );
2847 }
2848 }
2850 removeFace(replaceFace);
2852 // Modify ringEntities and replaceFace
2853 replaceFace = -1;
2854 ringEntities[replaceFaceIndex][indexI] = -1;
2855 }
2856 else
2857 {
2858 // Keep orientation intact, and update the owner
2859 owner_[replaceFace] = neighbour_[faceToRemove];
2860 }
2861 }
2862 else
2863 if (neighbour_[faceToRemove] == -1)
2864 {
2865 // This interior face would need to be converted
2866 // to a boundary one, but with orientation intact.
2867 face newFace = faces_[replaceFace];
2868 label newOwner = owner_[replaceFace];
2869 label newNeighbour = neighbour_[faceToRemove];
2870 labelList newFE = faceEdges_[replaceFace];
2872 label newFaceIndex =
2873 (
2874 insertFace
2875 (
2876 whichPatch(faceToRemove),
2877 newFace,
2878 newOwner,
2879 newNeighbour
2880 )
2881 );
2883 // Set this face aside for mapping
2884 if (!modifiedFaces.found(newFaceIndex))
2885 {
2886 modifiedFaces.insert(newFaceIndex);
2887 }
2889 // Ensure that all edges of this face are
2890 // on the boundary.
2891 forAll(newFE, edgeI)
2892 {
2893 if (whichEdgePatch(newFE[edgeI]) == -1)
2894 {
2895 edge newEdge = edges_[newFE[edgeI]];
2896 labelList newEF = edgeFaces_[newFE[edgeI]];
2897 labelList newEP = edgePoints_[newFE[edgeI]];
2899 // Need patch information for the new edge.
2900 // Find the corresponding edge in ringEntities.
2901 // Note that hullEdges doesn't need to be checked,
2902 // since they are common to both faces.
2903 label i =
2904 (
2905 findIndex
2906 (
2907 ringEntities[replaceEdgeIndex],
2908 newFE[edgeI]
2909 )
2910 );
2912 label repIndex =
2913 (
2914 whichEdgePatch
2915 (
2916 ringEntities[removeEdgeIndex][i]
2917 )
2918 );
2920 // Insert the new edge
2921 label newEdgeIndex =
2922 (
2923 insertEdge
2924 (
2925 repIndex,
2926 newEdge,
2927 newEF,
2928 newEP
2929 )
2930 );
2932 // Replace faceEdges for all
2933 // connected faces.
2934 forAll(newEF, faceI)
2935 {
2936 meshOps::replaceLabel
2937 (
2938 newFE[edgeI],
2939 newEdgeIndex,
2940 faceEdges_[newEF[faceI]]
2941 );
2942 }
2944 // Remove the edge
2945 removeEdge(newFE[edgeI]);
2947 // Replace faceEdges with new edge index
2948 newFE[edgeI] = newEdgeIndex;
2950 // Modify ringEntities
2951 ringEntities[replaceEdgeIndex][i] = newEdgeIndex;
2952 }
2953 }
2955 // Add the new faceEdges
2956 faceEdges_.append(newFE);
2958 // Replace edgeFaces with the new face index
2959 const labelList& newFEdges = faceEdges_[newFaceIndex];
2961 forAll(newFEdges, edgeI)
2962 {
2963 meshOps::replaceLabel
2964 (
2965 replaceFace,
2966 newFaceIndex,
2967 edgeFaces_[newFEdges[edgeI]]
2968 );
2969 }
2971 // Remove the face
2972 removeFace(replaceFace);
2974 // Replace label for the new owner
2975 meshOps::replaceLabel
2976 (
2977 replaceFace,
2978 newFaceIndex,
2979 cells_[newOwner]
2980 );
2982 // Modify ringEntities and replaceFace
2983 replaceFace = newFaceIndex;
2984 ringEntities[replaceFaceIndex][indexI] = newFaceIndex;
2985 }
2986 else
2987 {
2988 // Keep orientation intact, and update the neighbour
2989 neighbour_[replaceFace] = neighbour_[faceToRemove];
2990 }
2992 // Update the cell
2993 if (neighbour_[faceToRemove] != -1)
2994 {
2995 meshOps::replaceLabel
2996 (
2997 faceToRemove,
2998 replaceFace,
2999 cells_[neighbour_[faceToRemove]]
3000 );
3001 }
3002 }
3003 else
3004 {
3005 if (neighbour_[replaceFace] == cellToRemove)
3006 {
3007 if (owner_[faceToRemove] < owner_[replaceFace])
3008 {
3009 // This face is to be flipped
3010 faces_[replaceFace] = faces_[replaceFace].reverseFace();
3011 neighbour_[replaceFace] = owner_[replaceFace];
3012 owner_[replaceFace] = owner_[faceToRemove];
3014 setFlip(replaceFace);
3015 }
3016 else
3017 {
3018 // Keep orientation intact, and update the neighbour
3019 neighbour_[replaceFace] = owner_[faceToRemove];
3020 }
3021 }
3022 else
3023 {
3024 // Keep orientation intact, and update the owner
3025 owner_[replaceFace] = owner_[faceToRemove];
3026 }
3028 // Update the cell
3029 meshOps::replaceLabel
3030 (
3031 faceToRemove,
3032 replaceFace,
3033 cells_[owner_[faceToRemove]]
3034 );
3035 }
3036 }
3038 // Remove all hull entities
3039 forAll(faceHull, indexI)
3040 {
3041 label edgeToRemove = ringEntities[removeEdgeIndex][indexI];
3042 label faceToRemove = ringEntities[removeFaceIndex][indexI];
3043 label cellToRemove = cellHull[indexI];
3045 if (cellToRemove != -1)
3046 {
3047 // Remove faceToRemove and associated faceEdges
3048 removeFace(faceToRemove);
3050 // Remove the hull cell
3051 removeCell(cellToRemove);
3052 }
3054 // Remove the hull edge and associated edgeFaces
3055 removeEdge(edgeToRemove);
3057 // Remove the hull face
3058 removeFace(faceHull[indexI]);
3059 }
3061 // Loop through pointEdges for the collapsePoint,
3062 // and replace all occurrences with replacePoint.
3063 // Size-up pointEdges for the replacePoint as well.
3064 const labelList& pEdges = pointEdges_[collapsePoint];
3066 forAll(pEdges, edgeI)
3067 {
3068 // Renumber edges
3069 label edgeIndex = pEdges[edgeI];
3071 if (edgeIndex != eIndex)
3072 {
3073 if (debug > 2)
3074 {
3075 Info << "Renumbering [edge]: "
3076 << edgeIndex << ": "
3077 << edges_[edgeIndex] << endl;
3078 }
3080 if (edges_[edgeIndex][0] == collapsePoint)
3081 {
3082 edges_[edgeIndex][0] = replacePoint;
3084 meshOps::sizeUpList
3085 (
3086 edgeIndex,
3087 pointEdges_[replacePoint]
3088 );
3089 }
3090 else
3091 if (edges_[edgeIndex][1] == collapsePoint)
3092 {
3093 edges_[edgeIndex][1] = replacePoint;
3095 meshOps::sizeUpList
3096 (
3097 edgeIndex,
3098 pointEdges_[replacePoint]
3099 );
3100 }
3101 else
3102 {
3103 // Looks like pointEdges is inconsistent
3104 FatalErrorIn
3105 (
3106 "\n"
3107 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
3108 "(\n"
3109 " const label eIndex,\n"
3110 " label overRideCase,\n"
3111 " bool checkOnly,\n"
3112 " bool forceOp\n"
3113 ")\n"
3114 )
3115 << "pointEdges is inconsistent." << nl
3116 << "Point: " << collapsePoint << nl
3117 << "pointEdges: " << pEdges << nl
3118 << abort(FatalError);
3119 }
3121 // Loop through faces associated with this edge,
3122 // and renumber them as well.
3123 const labelList& eFaces = edgeFaces_[edgeIndex];
3125 forAll(eFaces, faceI)
3126 {
3127 const face& faceToCheck = faces_[eFaces[faceI]];
3129 if ((replaceIndex = faceToCheck.which(collapsePoint)) > -1)
3130 {
3131 if (debug > 2)
3132 {
3133 Info << "Renumbering face: "
3134 << eFaces[faceI] << ": "
3135 << faceToCheck << endl;
3136 }
3138 // Renumber the face...
3139 faces_[eFaces[faceI]][replaceIndex] = replacePoint;
3141 // Set this face aside for mapping
3142 if (!modifiedFaces.found(eFaces[faceI]))
3143 {
3144 modifiedFaces.insert(eFaces[faceI]);
3145 }
3147 // Look for an edge on this face that doesn't
3148 // contain collapsePoint or replacePoint.
3149 label rplIndex = -1;
3150 const labelList& fEdges = faceEdges_[eFaces[faceI]];
3152 forAll(fEdges, edgeI)
3153 {
3154 const edge& eCheck = edges_[fEdges[edgeI]];
3156 if
3157 (
3158 eCheck[0] != collapsePoint
3159 && eCheck[1] != collapsePoint
3160 && eCheck[0] != replacePoint
3161 && eCheck[1] != replacePoint
3162 )
3163 {
3164 rplIndex = fEdges[edgeI];
3165 break;
3166 }
3167 }
3169 // Modify edgePoints for this edge
3170 meshOps::replaceLabel
3171 (
3172 collapsePoint,
3173 replacePoint,
3174 edgePoints_[rplIndex]
3175 );
3176 }
3177 }
3178 }
3179 }
3181 // At this point, edgePoints for the replacement edges are broken,
3182 // but edgeFaces are consistent. So use this information to re-build
3183 // edgePoints for all replacement edges.
3184 forAll(ringEntities[replaceEdgeIndex], edgeI)
3185 {
3186 // If the ring edge was removed, don't bother.
3187 if (ringEntities[replaceEdgeIndex][edgeI] == -1)
3188 {
3189 continue;
3190 }
3192 if (debug > 2)
3193 {
3194 Info << "Building edgePoints for edge: "
3195 << ringEntities[replaceEdgeIndex][edgeI] << ": "
3196 << edges_[ringEntities[replaceEdgeIndex][edgeI]]
3197 << endl;
3198 }
3200 buildEdgePoints(ringEntities[replaceEdgeIndex][edgeI]);
3201 }
3203 // Move old / new points
3204 oldPoints_[replacePoint] = oldPoint;
3205 points_[replacePoint] = newPoint;
3207 // Remove the collapse point
3208 removePoint(collapsePoint);
3210 // Remove the edge
3211 removeEdge(eIndex);
3213 // For cell-mapping, exclude all hull-cells
3214 forAll(cellHull, indexI)
3215 {
3216 if (cellsChecked.found(cellHull[indexI]))
3217 {
3218 cellsChecked.erase(cellHull[indexI]);
3219 }
3220 }
3222 labelList mapCells = cellsChecked.toc();
3224 // Write out VTK files after change
3225 if (debug > 3)
3226 {
3227 writeVTK
3228 (
3229 Foam::name(eIndex)
3230 + '(' + Foam::name(edges_[eIndex][0])
3231 + ',' + Foam::name(edges_[eIndex][1]) + ')'
3232 + "_Collapse_1",
3233 mapCells
3234 );
3235 }
3237 // Specify that an old point-position
3238 // has been modified, if necessary
3239 if (collapseCase == 3)
3240 {
3241 labelList mP(2, -1);
3243 mP[0] = collapsePoint;
3244 mP[1] = replacePoint;
3246 modPoints_.set(replacePoint, mP);
3247 }
3249 // Now that all old / new cells possess correct connectivity,
3250 // compute mapping information.
3251 forAll(mapCells, cellI)
3252 {
3253 // Fill-in candidate mapping information
3254 labelList mC(1, mapCells[cellI]);
3256 // Set the mapping for this cell
3257 setCellMapping(mapCells[cellI], mC);
3258 }
3260 // Set face mapping information for modified faces
3261 forAllConstIter(labelHashSet, modifiedFaces, fIter)
3262 {
3263 // Exclude deleted faces
3264 if (faces_[fIter.key()].empty())
3265 {
3266 continue;
3267 }
3269 // Decide between default / weighted mapping
3270 // based on boundary information
3271 label fPatch = whichPatch(fIter.key());
3273 if (fPatch == -1)
3274 {
3275 setFaceMapping(fIter.key());
3276 }
3277 else
3278 {
3279 // Fill-in candidate mapping information
3280 labelList faceCandidates;
3282 const labelList& fEdges = faceEdges_[fIter.key()];
3284 forAll(fEdges, edgeI)
3285 {
3286 if (whichEdgePatch(fEdges[edgeI]) == fPatch)
3287 {
3288 // Loop through associated edgeFaces
3289 const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
3291 forAll(eFaces, faceI)
3292 {
3293 if
3294 (
3295 (eFaces[faceI] != fIter.key()) &&
3296 (whichPatch(eFaces[faceI]) == fPatch)
3297 )
3298 {
3299 faceCandidates.setSize
3300 (
3301 faceCandidates.size() + 1,
3302 eFaces[faceI]
3303 );
3304 }
3305 }
3306 }
3307 }
3309 // Set the mapping for this face
3310 setFaceMapping(fIter.key(), faceCandidates);
3311 }
3312 }
3314 // Set the flag
3315 topoChangeFlag_ = true;
3317 // Increment the counter
3318 statistics_[4]++;
3320 // Increment the number of modifications
3321 statistics_[0]++;
3323 // Return a succesful collapse
3324 map.type() = collapseCase;
3326 return map;
3327 }
3330 // Remove the specified cells from the mesh,
3331 // and add internal faces/edges to the specified patch
3332 const changeMap dynamicTopoFvMesh::removeCells
3333 (
3334 const labelList& cList,
3335 const label patch
3336 )
3337 {
3338 changeMap map;
3340 labelHashSet pointsToRemove, edgesToRemove, facesToRemove;
3341 Map<label> facesToConvert, edgesToConvert;
3343 // First loop through all cells and accumulate
3344 // a set of faces to be removed/converted.
3345 forAll(cList, cellI)
3346 {
3347 const cell& cellToCheck = cells_[cList[cellI]];
3349 forAll(cellToCheck, faceI)
3350 {
3351 label own = owner_[cellToCheck[faceI]];
3352 label nei = neighbour_[cellToCheck[faceI]];
3354 if (nei == -1)
3355 {
3356 if (!facesToRemove.found(cellToCheck[faceI]))
3357 {
3358 facesToRemove.insert(cellToCheck[faceI]);
3359 }
3360 }
3361 else
3362 if
3363 (
3364 (findIndex(cList, own) != -1) &&
3365 (findIndex(cList, nei) != -1)
3366 )
3367 {
3368 if (!facesToRemove.found(cellToCheck[faceI]))
3369 {
3370 facesToRemove.insert(cellToCheck[faceI]);
3371 }
3372 }
3373 else
3374 {
3375 facesToConvert.set(cellToCheck[faceI], -1);
3376 }
3377 }
3378 }
3380 // Add all edges as candidates for conversion.
3381 // Some of these will be removed altogether.
3382 forAllConstIter(labelHashSet, facesToRemove, fIter)
3383 {
3384 const labelList& fEdges = faceEdges_[fIter.key()];
3386 forAll(fEdges, edgeI)
3387 {
3388 if (whichEdgePatch(fEdges[edgeI]) == patch)
3389 {
3390 // Make an identical map
3391 edgesToConvert.set(fEdges[edgeI], fEdges[edgeI]);
3392 }
3393 else
3394 {
3395 edgesToConvert.set(fEdges[edgeI], -1);
3396 }
3397 }
3398 }
3400 forAllConstIter(Map<label>, facesToConvert, fIter)
3401 {
3402 const labelList& fEdges = faceEdges_[fIter.key()];
3404 forAll(fEdges, edgeI)
3405 {
3406 if (whichEdgePatch(fEdges[edgeI]) == patch)
3407 {
3408 // Make an identical map
3409 edgesToConvert.set(fEdges[edgeI], fEdges[edgeI]);
3410 }
3411 else
3412 {
3413 edgesToConvert.set(fEdges[edgeI], -1);
3414 }
3415 }
3416 }
3418 // Build a list of edges to be removed.
3419 forAllConstIter(Map<label>, edgesToConvert, eIter)
3420 {
3421 const labelList& eFaces = edgeFaces_[eIter.key()];
3423 bool allRemove = true;
3425 forAll(eFaces, faceI)
3426 {
3427 if (facesToConvert.found(eFaces[faceI]))
3428 {
3429 allRemove = false;
3430 break;
3431 }
3432 }
3434 if (allRemove)
3435 {
3436 if (!edgesToRemove.found(eIter.key()))
3437 {
3438 edgesToRemove.insert(eIter.key());
3439 }
3440 }
3441 }
3443 // Weed-out the conversion list.
3444 forAllConstIter(labelHashSet, edgesToRemove, eIter)
3445 {
3446 edgesToConvert.erase(eIter.key());
3447 }
3449 // Build a set of points to be removed.
3450 if (!twoDMesh_)
3451 {
3452 forAllConstIter(labelHashSet, edgesToRemove, eIter)
3453 {
3454 const edge& edgeToCheck = edges_[eIter.key()];
3456 forAll(edgeToCheck, pointI)
3457 {
3458 const labelList& pEdges = pointEdges_[edgeToCheck[pointI]];
3460 bool allRemove = true;
3462 forAll(pEdges, edgeI)
3463 {
3464 if (!edgesToRemove.found(pEdges[edgeI]))
3465 {
3466 allRemove = false;
3467 break;
3468 }
3469 }
3471 if (allRemove)
3472 {
3473 if (!pointsToRemove.found(edgeToCheck[pointI]))
3474 {
3475 pointsToRemove.insert(edgeToCheck[pointI]);
3476 }
3477 }
3478 }
3479 }
3480 }
3482 forAllIter(Map<label>, edgesToConvert, eIter)
3483 {
3484 const labelList& eFaces = edgeFaces_[eIter.key()];
3486 label nConvFaces = 0;
3488 forAll(eFaces, faceI)
3489 {
3490 if (facesToConvert.found(eFaces[faceI]))
3491 {
3492 nConvFaces++;
3493 }
3494 }
3496 if (nConvFaces > 2)
3497 {
3498 Info << "Invalid conversion. Bailing out." << endl;
3499 return map;
3500 }
3501 }
3503 // Write out candidates for post-processing
3504 if (debug > 2)
3505 {
3506 writeVTK("pointsToRemove", pointsToRemove.toc(), 0);
3507 writeVTK("edgesToRemove", edgesToRemove.toc(), 1);
3508 writeVTK("facesToRemove", facesToRemove.toc(), 2);
3509 writeVTK("cellsToRemove", cList, 3);
3510 writeVTK("edgesToConvert", edgesToConvert.toc(), 1);
3511 writeVTK("facesToConvert", facesToConvert.toc(), 2);
3512 }
3514 // Loop through all faces for conversion, check orientation
3515 // and create new faces in their place.
3516 forAllIter(Map<label>, facesToConvert, fIter)
3517 {
3518 // Check if this internal face is oriented properly.
3519 face newFace;
3520 label newOwner = -1;
3521 labelList fEdges = faceEdges_[fIter.key()];
3523 if (findIndex(cList, neighbour_[fIter.key()]) != -1)
3524 {
3525 // Orientation is correct
3526 newFace = faces_[fIter.key()];
3527 newOwner = owner_[fIter.key()];
3528 }
3529 else
3530 if (findIndex(cList, owner_[fIter.key()]) != -1)
3531 {
3532 // Face is to be reversed.
3533 newFace = faces_[fIter.key()].reverseFace();
3534 newOwner = neighbour_[fIter.key()];
3536 setFlip(fIter.key());
3537 }
3538 else
3539 {
3540 // Something's terribly wrong
3541 FatalErrorIn
3542 (
3543 "\n"
3544 "const changeMap dynamicTopoFvMesh::removeCells\n"
3545 "(\n"
3546 " const labelList& cList,\n"
3547 " const label patch\n"
3548 ")\n"
3549 )
3550 << nl << " Invalid mesh. "
3551 << abort(FatalError);
3552 }
3554 // Insert the reconfigured face at the boundary.
3555 fIter() =
3556 (
3557 insertFace
3558 (
3559 patch,
3560 newFace,
3561 newOwner,
3562 -1
3563 )
3564 );
3566 // Add the faceEdges entry.
3567 // Edges will be corrected later.
3568 faceEdges_.append(fEdges);
3570 // Add this face to the map.
3571 map.addFace(fIter());
3573 // Replace cell with the new face label
3574 meshOps::replaceLabel
3575 (
3576 fIter.key(),
3577 fIter(),
3578 cells_[newOwner]
3579 );
3581 // Remove the internal face.
3582 removeFace(fIter.key());
3583 }
3585 // Create a new edge for each converted edge
3586 forAllIter(Map<label>, edgesToConvert, eIter)
3587 {
3588 if (eIter() == -1)
3589 {
3590 // Create copies before appending.
3591 edge newEdge = edges_[eIter.key()];
3592 labelList eFaces = edgeFaces_[eIter.key()];
3593 labelList ePoints = edgePoints_[eIter.key()];
3595 eIter() =
3596 (
3597 insertEdge
3598 (
3599 patch,
3600 newEdge,
3601 eFaces,
3602 ePoints
3603 )
3604 );
3606 // Add this edge to the map.
3607 map.addEdge(eIter());
3609 // Remove the edge
3610 removeEdge(eIter.key());
3611 }
3612 }
3614 // Loop through all faces for conversion, and replace edgeFaces.
3615 forAllConstIter(Map<label>, facesToConvert, fIter)
3616 {
3617 // Make a copy, because this list is going to
3618 // be modified within this loop.
3619 labelList fEdges = faceEdges_[fIter()];
3621 forAll(fEdges, edgeI)
3622 {
3623 if (edgesToConvert.found(fEdges[edgeI]))
3624 {
3625 meshOps::replaceLabel
3626 (
3627 fIter.key(),
3628 fIter(),
3629 edgeFaces_[edgesToConvert[fEdges[edgeI]]]
3630 );
3632 meshOps::replaceLabel
3633 (
3634 fEdges[edgeI],
3635 edgesToConvert[fEdges[edgeI]],
3636 faceEdges_[fIter()]
3637 );
3638 }
3639 }
3640 }
3642 // Loop through all edges for conversion, and size-down edgeFaces.
3643 forAllConstIter(Map<label>, edgesToConvert, eIter)
3644 {
3645 // Make a copy, because this list is going to
3646 // be modified within this loop.
3647 labelList eFaces = edgeFaces_[eIter()];
3649 forAll(eFaces, faceI)
3650 {
3651 if (facesToRemove.found(eFaces[faceI]))
3652 {
3653 meshOps::sizeDownList
3654 (
3655 eFaces[faceI],
3656 edgeFaces_[eIter()]
3657 );
3658 }
3660 // Replace old edges with new ones.
3661 labelList& fEdges = faceEdges_[eFaces[faceI]];
3663 forAll(fEdges, edgeI)
3664 {
3665 if (edgesToConvert.found(fEdges[edgeI]))
3666 {
3667 fEdges[edgeI] = edgesToConvert[fEdges[edgeI]];
3668 }
3669 }
3670 }
3671 }
3673 // At this point, edgeFaces is consistent.
3674 // Correct edge-points for all converted edges
3675 if (!twoDMesh_)
3676 {
3677 forAllConstIter(Map<label>, edgesToConvert, eIter)
3678 {
3679 buildEdgePoints(eIter());
3680 }
3681 }
3683 // Remove unwanted faces
3684 forAllConstIter(labelHashSet, facesToRemove, fIter)
3685 {
3686 removeFace(fIter.key());
3687 }
3689 // Remove unwanted edges
3690 forAllConstIter(labelHashSet, edgesToRemove, eIter)
3691 {
3692 removeEdge(eIter.key());
3693 }
3695 // Remove unwanted points
3696 forAllConstIter(labelHashSet, pointsToRemove, pIter)
3697 {
3698 removePoint(pIter.key());
3699 }
3701 // Remove all cells
3702 forAll(cList, cellI)
3703 {
3704 removeCell(cList[cellI]);
3705 }
3707 // Set the flag
3708 topoChangeFlag_ = true;
3710 return map;
3711 }
3714 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
3716 } // End namespace Foam
3718 // ************************************************************************* //