| blob | b4b6a63857bde860e47aba76eae4a69f1dc5077c |
1 /*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | foam-extend: Open Source CFD
4 \\ / O peration | Version: 3.2
5 \\ / A nd | Web: http://www.foam-extend.org
6 \\/ M anipulation | For copyright notice see file Copyright
7 -------------------------------------------------------------------------------
8 License
9 This file is part of foam-extend.
11 foam-extend 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 3 of the License, or (at your
14 option) any later version.
16 foam-extend is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 General Public License for more details.
21 You should have received a copy of the GNU General Public License
22 along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
24 \*---------------------------------------------------------------------------*/
26 #include "Stack.H"
27 #include "triFace.H"
28 #include "objectMap.H"
29 #include "changeMap.H"
30 #include "coupledInfo.H"
31 #include "multiThreader.H"
32 #include "dynamicTopoFvMesh.H"
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 namespace Foam
37 {
39 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
41 // Method to collapse a quad-face in 2D
42 // - Returns a changeMap with a type specifying:
43 // -3: Collapse failed since face was on a noRefinement patch.
44 // -1: Collapse failed since max number of topo-changes was reached.
45 // 0: Collapse could not be performed.
46 // 1: Collapsed to first node.
47 // 2: Collapsed to second node.
48 // 3: Collapse to mid-point.
49 // - overRideCase is used to force a certain collapse configuration.
50 // -1: Use this value to let collapseQuadFace decide a case.
51 // 1: Force collapse to first node.
52 // 2: Force collapse to second node.
53 // 3: Force collapse to mid-point.
54 // - checkOnly performs a feasibility check and returns without modifications.
55 const changeMap dynamicTopoFvMesh::collapseQuadFace
56 (
57 const label fIndex,
58 label overRideCase,
59 bool checkOnly,
60 bool forceOp
61 )
62 {
63 // Figure out which thread this is...
64 label tIndex = self();
66 // Prepare the changeMaps
67 changeMap map;
68 List<changeMap> slaveMaps;
69 bool collapsingSlave = false;
71 if
72 (
73 (status(TOTAL_MODIFICATIONS) > maxModifications_)
74 && (maxModifications_ > -1)
75 )
76 {
77 // Reached the max allowable topo-changes.
78 stack(tIndex).clear();
80 return map;
81 }
83 // Check if edgeRefinements are to be avoided on patch.
84 if (baseMesh_.lengthEstimator().checkRefinementPatch(whichPatch(fIndex)))
85 {
86 map.type() = -3;
88 return map;
89 }
91 // Sanity check: Is the index legitimate?
92 if (fIndex < 0)
93 {
94 FatalErrorIn
95 (
96 "\n"
97 "const changeMap "
98 "dynamicTopoFvMesh::collapseQuadFace\n"
99 "(\n"
100 " const label fIndex,\n"
101 " label overRideCase,\n"
102 " bool checkOnly,\n"
103 " bool forceOp\n"
104 ")\n"
105 )
106 << " Invalid index: " << fIndex
107 << abort(FatalError);
108 }
110 // Define the edges on the face to be collapsed
111 FixedList<edge,4> checkEdge(edge(-1, -1));
112 FixedList<label,4> checkEdgeIndex(-1);
114 // Define checkEdges
115 getCheckEdges(fIndex, (*this), map, checkEdge, checkEdgeIndex);
117 // Determine the common vertices for the first and second edges
118 label cv0 = checkEdge[1].commonVertex(checkEdge[0]);
119 label cv1 = checkEdge[1].commonVertex(checkEdge[3]);
120 label cv2 = checkEdge[2].commonVertex(checkEdge[0]);
121 label cv3 = checkEdge[2].commonVertex(checkEdge[3]);
123 // If coupled modification is set, and this is a
124 // master face, collapse its slaves first.
125 bool localCouple = false, procCouple = false;
127 if (coupledModification_)
128 {
129 const face& fCheck = faces_[fIndex];
131 const label faceEnum = coupleMap::FACE;
132 const label pointEnum = coupleMap::POINT;
134 // Is this a locally coupled edge (either master or slave)?
135 if (locallyCoupledEntity(fIndex, true))
136 {
137 localCouple = true;
138 procCouple = false;
139 }
140 else
141 if (processorCoupledEntity(fIndex))
142 {
143 procCouple = true;
144 localCouple = false;
145 }
147 if (localCouple && !procCouple)
148 {
149 // Determine the slave index.
150 label sIndex = -1, pIndex = -1;
152 forAll(patchCoupling_, patchI)
153 {
154 if (patchCoupling_(patchI))
155 {
156 const coupleMap& cMap = patchCoupling_[patchI].map();
158 if ((sIndex = cMap.findSlave(faceEnum, fIndex)) > -1)
159 {
160 pIndex = patchI;
162 break;
163 }
165 // The following bit happens only during the sliver
166 // exudation process, since slave edges are
167 // usually not added to the coupled edge-stack.
168 if ((sIndex = cMap.findMaster(faceEnum, fIndex)) > -1)
169 {
170 pIndex = patchI;
172 // Notice that we are collapsing a slave edge first.
173 collapsingSlave = true;
175 break;
176 }
177 }
178 }
180 if (sIndex == -1)
181 {
182 FatalErrorIn
183 (
184 "\n"
185 "const changeMap "
186 "dynamicTopoFvMesh::collapseQuadFace\n"
187 "(\n"
188 " const label fIndex,\n"
189 " label overRideCase,\n"
190 " bool checkOnly,\n"
191 " bool forceOp\n"
192 ")\n"
193 )
194 << "Coupled maps were improperly specified." << nl
195 << " Slave index not found for: " << nl
196 << " Face: " << fIndex << ": " << fCheck << nl
197 << abort(FatalError);
198 }
199 else
200 {
201 // If we've found the slave, size up the list
202 meshOps::sizeUpList
203 (
204 changeMap(),
205 slaveMaps
206 );
208 // Save index and patch for posterity
209 slaveMaps[0].index() = sIndex;
210 slaveMaps[0].patchIndex() = pIndex;
211 }
213 if (debug > 1)
214 {
215 Pout<< nl << " >> Collapsing slave face: " << sIndex
216 << " for master face: " << fIndex << endl;
217 }
218 }
219 else
220 if (procCouple && !localCouple)
221 {
222 // If this is a new entity, bail out for now.
223 // This will be handled at the next time-step.
224 if (fIndex >= nOldFaces_)
225 {
226 return map;
227 }
229 // Check slaves
230 forAll(procIndices_, pI)
231 {
232 // Fetch reference to subMeshes
233 const coupledMesh& sendMesh = sendMeshes_[pI];
234 const coupledMesh& recvMesh = recvMeshes_[pI];
236 const coupleMap& scMap = sendMesh.map();
237 const coupleMap& rcMap = recvMesh.map();
239 // If this face was sent to a lower-ranked
240 // processor, skip it.
241 if
242 (
243 priority
244 (
245 procIndices_[pI],
246 lessOp<label>(),
247 Pstream::myProcNo()
248 )
249 )
250 {
251 if (scMap.reverseEntityMap(faceEnum).found(fIndex))
252 {
253 if (debug > 3)
254 {
255 Pout<< "Face: " << fIndex
256 << "::" << fCheck
257 << " was sent to proc: "
258 << procIndices_[pI]
259 << ", so bailing out."
260 << endl;
261 }
263 return map;
264 }
265 }
267 label sIndex = -1;
269 if ((sIndex = rcMap.findSlave(faceEnum, fIndex)) > -1)
270 {
271 // Check if a lower-ranked processor is
272 // handling this face
273 if
274 (
275 priority
276 (
277 procIndices_[pI],
278 lessOp<label>(),
279 Pstream::myProcNo()
280 )
281 )
282 {
283 if (debug > 3)
284 {
285 Pout<< "Face: " << fIndex
286 << "::" << fCheck
287 << " is handled by proc: "
288 << procIndices_[pI]
289 << ", so bailing out."
290 << endl;
291 }
293 return map;
294 }
296 label curIndex = slaveMaps.size();
298 // Size up the list
299 meshOps::sizeUpList
300 (
301 changeMap(),
302 slaveMaps
303 );
305 // Save index and patch for posterity
306 slaveMaps[curIndex].index() = sIndex;
307 slaveMaps[curIndex].patchIndex() = pI;
308 }
309 else
310 if
311 (
312 (
313 rcMap.findSlave(pointEnum, cv0) > -1 &&
314 rcMap.findSlave(pointEnum, cv1) > -1
315 )
316 || (
317 rcMap.findSlave(pointEnum, cv2) > -1 &&
318 rcMap.findSlave(pointEnum, cv3) > -1
319 )
320 )
321 {
322 // An edge-only coupling exists.
324 // Check if a lower-ranked processor is
325 // handling this face
326 if
327 (
328 priority
329 (
330 procIndices_[pI],
331 lessOp<label>(),
332 Pstream::myProcNo()
333 )
334 )
335 {
336 if (debug > 3)
337 {
338 Pout<< "Face edge on: " << fIndex
339 << "::" << fCheck
340 << " is handled by proc: "
341 << procIndices_[pI]
342 << ", so bailing out."
343 << endl;
344 }
346 return map;
347 }
349 label p0 = rcMap.findSlave(pointEnum, cv0);
350 label p1 = rcMap.findSlave(pointEnum, cv1);
352 label p2 = rcMap.findSlave(pointEnum, cv2);
353 label p3 = rcMap.findSlave(pointEnum, cv3);
355 edge cEdge(-1, -1);
356 label edgeCouple = -1;
358 if ((p0 > -1 && p1 > -1) && (p2 == -1 && p3 == -1))
359 {
360 cEdge[0] = p0;
361 cEdge[1] = p1;
363 edgeCouple = 1;
364 sIndex = readLabel(map.lookup("firstEdge"));
365 }
366 else
367 if ((p0 == -1 && p1 == -1) && (p2 > -1 && p3 > -1))
368 {
369 cEdge[0] = p2;
370 cEdge[1] = p3;
372 edgeCouple = 2;
373 sIndex = readLabel(map.lookup("secondEdge"));
374 }
376 label curIndex = slaveMaps.size();
378 // Size up the list
379 meshOps::sizeUpList
380 (
381 changeMap(),
382 slaveMaps
383 );
385 // Unfortunately, since no edge maps are
386 // available in 2D, loop through boundary
387 // faces on subMesh and get the right edge
388 label seIndex = -1;
390 const dynamicTopoFvMesh& sMesh = recvMesh.subMesh();
392 for
393 (
394 label faceI = sMesh.nOldInternalFaces_;
395 faceI < sMesh.faceEdges_.size();
396 faceI++
397 )
398 {
399 const labelList& fEdges = sMesh.faceEdges_[faceI];
401 forAll(fEdges, edgeI)
402 {
403 if (sMesh.edges_[fEdges[edgeI]] == cEdge)
404 {
405 seIndex = fEdges[edgeI];
406 break;
407 }
408 }
410 if (seIndex > -1)
411 {
412 break;
413 }
414 }
416 if (seIndex == -1)
417 {
418 Pout<< "Face edge on: " << fIndex
419 << "::" << fCheck
420 << " sIndex: " << sIndex
421 << " could not be found."
422 << abort(FatalError);
423 }
425 // Save index and patch for posterity
426 // - Negate the index to signify edge coupling
427 slaveMaps[curIndex].index() = -seIndex;
428 slaveMaps[curIndex].patchIndex() = pI;
430 // Save edgeCouple as well, so that
431 // another map comparison is avoided.
432 slaveMaps[curIndex].type() = edgeCouple;
433 }
434 }
435 }
436 else
437 {
438 // Something's wrong with coupling maps
439 FatalErrorIn
440 (
441 "\n"
442 "const changeMap "
443 "dynamicTopoFvMesh::collapseQuadFace\n"
444 "(\n"
445 " const label fIndex,\n"
446 " label overRideCase,\n"
447 " bool checkOnly,\n"
448 " bool forceOp\n"
449 ")\n"
450 )
451 << "Coupled maps were improperly specified." << nl
452 << " localCouple: " << localCouple << nl
453 << " procCouple: " << procCouple << nl
454 << " Face: " << fIndex << ": " << fCheck << nl
455 << abort(FatalError);
456 }
458 // Temporarily turn off coupledModification
459 unsetCoupledModification();
461 // Test the master face for collapse, and decide on a case
462 changeMap masterMap = collapseQuadFace(fIndex, -1, true, forceOp);
464 // Turn it back on.
465 setCoupledModification();
467 // Master couldn't perform collapse
468 if (masterMap.type() <= 0)
469 {
470 return masterMap;
471 }
473 // For edge-only coupling, define the points for checking
474 List<FixedList<point, 2> > slaveMoveNewPoint(slaveMaps.size());
475 List<FixedList<point, 2> > slaveMoveOldPoint(slaveMaps.size());
477 // Now check each of the slaves for collapse feasibility
478 forAll(slaveMaps, slaveI)
479 {
480 // Alias for convenience...
481 changeMap& slaveMap = slaveMaps[slaveI];
483 label slaveOverRide = -1;
484 label sIndex = slaveMap.index();
485 label pI = slaveMap.patchIndex();
486 const coupleMap* cMapPtr = NULL;
488 if (localCouple)
489 {
490 cMapPtr = &(patchCoupling_[pI].map());
491 }
492 else
493 if (procCouple)
494 {
495 const dynamicTopoFvMesh& sMesh =
496 (
497 recvMeshes_[pI].subMesh()
498 );
500 cMapPtr = &(recvMeshes_[pI].map());
502 if (debug > 3)
503 {
504 if (sIndex < 0)
505 {
506 Pout<< "Checking slave edge: " << mag(sIndex)
507 << "::" << sMesh.edges_[mag(sIndex)]
508 << " on proc: " << procIndices_[pI]
509 << " for master face: " << fIndex
510 << " using collapseCase: " << masterMap.type()
511 << endl;
512 }
513 else
514 {
515 Pout<< "Checking slave face: " << sIndex
516 << "::" << sMesh.faces_[sIndex]
517 << " on proc: " << procIndices_[pI]
518 << " for master face: " << fIndex
519 << " using collapseCase: " << masterMap.type()
520 << endl;
521 }
522 }
523 }
525 // Define an overRide case for the slave
526 FixedList<edge, 2> mEdge(edge(-1, -1)), sEdge(edge(-1, -1));
528 if (collapsingSlave)
529 {
530 const Map<label>& rPointMap =
531 (
532 cMapPtr->reverseEntityMap(pointEnum)
533 );
535 if (sIndex < 0)
536 {
537 if (slaveMap.type() == 1)
538 {
539 mEdge[0][0] = rPointMap[cv0];
540 mEdge[0][1] = rPointMap[cv1];
541 }
542 else
543 if (slaveMap.type() == 2)
544 {
545 mEdge[0][0] = rPointMap[cv2];
546 mEdge[0][1] = rPointMap[cv3];
547 }
548 }
549 else
550 {
551 mEdge[0][0] = rPointMap[cv0];
552 mEdge[0][1] = rPointMap[cv1];
554 mEdge[1][0] = rPointMap[cv2];
555 mEdge[1][1] = rPointMap[cv3];
556 }
557 }
558 else
559 {
560 const Map<label>& pointMap =
561 (
562 cMapPtr->entityMap(pointEnum)
563 );
565 if (sIndex < 0)
566 {
567 if (slaveMap.type() == 1)
568 {
569 mEdge[0][0] = pointMap[cv0];
570 mEdge[0][1] = pointMap[cv1];
571 }
572 else
573 if (slaveMap.type() == 2)
574 {
575 mEdge[0][0] = pointMap[cv2];
576 mEdge[0][1] = pointMap[cv3];
577 }
578 }
579 else
580 {
581 mEdge[0][0] = pointMap[cv0];
582 mEdge[0][1] = pointMap[cv1];
584 mEdge[1][0] = pointMap[cv2];
585 mEdge[1][1] = pointMap[cv3];
586 }
587 }
589 // Determine checkEdges for the slave
590 FixedList<edge,4> slaveCheckEdge(edge(-1, -1));
591 FixedList<label,4> slaveCheckEdgeIndex(-1);
593 if (localCouple)
594 {
595 getCheckEdges
596 (
597 sIndex,
598 (*this),
599 slaveMap,
600 slaveCheckEdge,
601 slaveCheckEdgeIndex
602 );
604 sEdge[0] = edges_[readLabel(slaveMap.lookup("firstEdge"))];
605 sEdge[1] = edges_[readLabel(slaveMap.lookup("secondEdge"))];
606 }
607 else
608 if (procCouple)
609 {
610 const dynamicTopoFvMesh& sMesh =
611 (
612 recvMeshes_[pI].subMesh()
613 );
615 if (sIndex < 0)
616 {
617 sEdge[0] = sMesh.edges_[mag(sIndex)];
618 }
619 else
620 {
621 getCheckEdges
622 (
623 sIndex,
624 sMesh,
625 slaveMap,
626 slaveCheckEdge,
627 slaveCheckEdgeIndex
628 );
630 sEdge[0] =
631 (
632 sMesh.edges_[readLabel(slaveMap.lookup("firstEdge"))]
633 );
635 sEdge[1] =
636 (
637 sMesh.edges_[readLabel(slaveMap.lookup("secondEdge"))]
638 );
639 }
640 }
642 // Compare edge orientations for edge-only coupling
643 label compVal = -2;
645 // Perform a topological comparison.
646 switch (masterMap.type())
647 {
648 case 1:
649 {
650 if (sIndex < 0)
651 {
652 slaveMoveNewPoint[slaveI][0] = points_[cv0];
653 slaveMoveNewPoint[slaveI][1] = points_[cv1];
655 slaveMoveOldPoint[slaveI][0] = oldPoints_[cv0];
656 slaveMoveOldPoint[slaveI][1] = oldPoints_[cv1];
657 }
658 else
659 if (mEdge[0] == sEdge[0])
660 {
661 slaveOverRide = 1;
662 }
663 else
664 if (mEdge[1] == sEdge[0])
665 {
666 slaveOverRide = 2;
667 }
668 else
669 {
670 // Write out for for post-processing
671 writeVTK("mFace_" + Foam::name(fIndex), fIndex, 2);
673 FatalErrorIn
674 (
675 "\n"
676 "const changeMap "
677 "dynamicTopoFvMesh::collapseQuadFace\n"
678 "(\n"
679 " const label fIndex,\n"
680 " label overRideCase,\n"
681 " bool checkOnly,\n"
682 " bool forceOp\n"
683 ")\n"
684 )
685 << "Coupled collapse failed." << nl
686 << "Masters: " << nl
687 << checkEdgeIndex[1] << ": "
688 << checkEdge[1] << nl
689 << checkEdgeIndex[2] << ": "
690 << checkEdge[2] << nl
691 << "Slaves: " << nl
692 << readLabel(slaveMap.lookup("firstEdge")) << ": "
693 << sEdge[0] << nl
694 << readLabel(slaveMap.lookup("secondEdge")) << ": "
695 << sEdge[1] << nl
696 << abort(FatalError);
697 }
699 break;
700 }
702 case 2:
703 {
704 if (sIndex < 0)
705 {
706 slaveMoveNewPoint[slaveI][0] = points_[cv2];
707 slaveMoveNewPoint[slaveI][1] = points_[cv3];
709 slaveMoveOldPoint[slaveI][0] = oldPoints_[cv2];
710 slaveMoveOldPoint[slaveI][1] = oldPoints_[cv3];
711 }
712 else
713 if (mEdge[1] == sEdge[1])
714 {
715 slaveOverRide = 2;
716 }
717 else
718 if (mEdge[0] == sEdge[1])
719 {
720 slaveOverRide = 1;
721 }
722 else
723 {
724 // Write out for for post-processing
725 writeVTK("mFace_" + Foam::name(fIndex), fIndex, 2);
727 FatalErrorIn
728 (
729 "\n"
730 "const changeMap "
731 "dynamicTopoFvMesh::collapseQuadFace\n"
732 "(\n"
733 " const label fIndex,\n"
734 " label overRideCase,\n"
735 " bool checkOnly,\n"
736 " bool forceOp\n"
737 ")\n"
738 )
739 << "Coupled collapse failed." << nl
740 << "Masters: " << nl
741 << checkEdgeIndex[1] << ": "
742 << checkEdge[1] << nl
743 << checkEdgeIndex[2] << ": "
744 << checkEdge[2] << nl
745 << "Slaves: " << nl
746 << readLabel(slaveMap.lookup("firstEdge")) << ": "
747 << sEdge[0] << nl
748 << readLabel(slaveMap.lookup("secondEdge")) << ": "
749 << sEdge[1] << nl
750 << abort(FatalError);
751 }
753 break;
754 }
756 case 3:
757 {
758 if (sIndex < 0)
759 {
760 slaveMoveNewPoint[slaveI][0] =
761 (
762 0.5 * (points_[cv0] + points_[cv2])
763 );
765 slaveMoveNewPoint[slaveI][1] =
766 (
767 0.5 * (points_[cv1] + points_[cv3])
768 );
770 slaveMoveOldPoint[slaveI][0] =
771 (
772 0.5 * (oldPoints_[cv0] + oldPoints_[cv2])
773 );
775 slaveMoveOldPoint[slaveI][1] =
776 (
777 0.5 * (oldPoints_[cv1] + oldPoints_[cv3])
778 );
779 }
780 else
781 {
782 overRideCase = 3;
783 }
785 break;
786 }
787 }
789 if (sIndex < 0)
790 {
791 // Check edge orientation
792 compVal = edge::compare(mEdge[0], sEdge[0]);
794 // Swap components if necessary
795 if (compVal == -1)
796 {
797 Foam::Swap
798 (
799 slaveMoveNewPoint[slaveI][0],
800 slaveMoveNewPoint[slaveI][1]
801 );
803 Foam::Swap
804 (
805 slaveMoveOldPoint[slaveI][0],
806 slaveMoveOldPoint[slaveI][1]
807 );
808 }
809 else
810 if (compVal != 1)
811 {
812 FatalErrorIn
813 (
814 "\n"
815 "const changeMap "
816 "dynamicTopoFvMesh::collapseQuadFace\n"
817 "(\n"
818 " const label fIndex,\n"
819 " label overRideCase,\n"
820 " bool checkOnly,\n"
821 " bool forceOp\n"
822 ")\n"
823 )
824 << "Coupled topo-change for slave failed." << nl
825 << " Dissimilar edges: " << nl
826 << " mEdge: " << mEdge << nl
827 << " sEdge: " << sEdge << nl
828 << " masterMap type: " << masterMap.type() << nl
829 << abort(FatalError);
830 }
831 }
833 // Temporarily turn off coupledModification
834 unsetCoupledModification();
836 // Test the slave face
837 if (localCouple)
838 {
839 slaveMap =
840 (
841 collapseQuadFace(sIndex, slaveOverRide, true, forceOp)
842 );
843 }
844 else
845 if (procCouple)
846 {
847 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
849 if (sIndex < 0)
850 {
851 // Edge-based coupling
853 // Build a hull of cells and tri-faces
854 // that are connected to the slave edge
855 labelList slaveHullCells;
856 labelList slaveHullTriFaces;
858 meshOps::constructPrismHull
859 (
860 mag(sIndex),
861 sMesh.faces_,
862 sMesh.cells_,
863 sMesh.owner_,
864 sMesh.neighbour_,
865 sMesh.edgeFaces_,
866 slaveHullTriFaces,
867 slaveHullCells
868 );
870 bool infeasible = false;
872 // Keep track of resulting cell quality,
873 // if collapse is indeed feasible
874 scalar slaveCollapseQuality(GREAT);
876 // Check whether the collapse is possible.
877 if
878 (
879 checkCollapse
880 (
881 sMesh,
882 slaveHullTriFaces,
883 FixedList<label, 2>(-1),
884 FixedList<label, 2>(-1),
885 sEdge[0],
886 slaveMoveNewPoint[slaveI],
887 slaveMoveOldPoint[slaveI],
888 slaveCollapseQuality,
889 false,
890 forceOp
891 )
892 )
893 {
894 infeasible = true;
895 }
897 if (infeasible)
898 {
899 slaveMap.type() = 0;
900 }
901 else
902 {
903 slaveMap.type() = 1;
904 }
905 }
906 else
907 {
908 // Edge-based coupling
909 slaveMap =
910 (
911 sMesh.collapseQuadFace
912 (
913 sIndex,
914 slaveOverRide,
915 true,
916 forceOp
917 )
918 );
919 }
920 }
922 // Turn it back on.
923 setCoupledModification();
925 if (slaveMap.type() <= 0)
926 {
927 // Slave couldn't perform collapse.
928 map.type() = -2;
930 return map;
931 }
933 // Save index and patch for posterity
934 slaveMap.index() = sIndex;
935 slaveMap.patchIndex() = pI;
936 }
938 // Next collapse each slave face (for real this time...)
939 forAll(slaveMaps, slaveI)
940 {
941 // Alias for convenience...
942 changeMap& slaveMap = slaveMaps[slaveI];
944 label sIndex = slaveMap.index();
945 label pI = slaveMap.patchIndex();
946 label slaveOverRide = slaveMap.type();
948 // Temporarily turn off coupledModification
949 unsetCoupledModification();
951 // Collapse the slave
952 if (localCouple)
953 {
954 slaveMap =
955 (
956 collapseQuadFace(sIndex, slaveOverRide, false, forceOp)
957 );
958 }
959 else
960 {
961 const coupleMap& cMap = recvMeshes_[pI].map();
962 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
964 if (sIndex < 0)
965 {
966 // Edge-based coupling
968 // Fetch the slave edge
969 edge sEdge = sMesh.edges_[mag(sIndex)];
971 // Move points to new location,
972 // and update operation into coupleMap
973 sMesh.points_[sEdge[0]] = slaveMoveNewPoint[slaveI][0];
974 sMesh.points_[sEdge[1]] = slaveMoveNewPoint[slaveI][1];
976 sMesh.oldPoints_[sEdge[0]] = slaveMoveOldPoint[slaveI][0];
977 sMesh.oldPoints_[sEdge[1]] = slaveMoveOldPoint[slaveI][1];
979 cMap.pushOperation
980 (
981 sEdge[0],
982 coupleMap::MOVE_POINT,
983 slaveMoveNewPoint[slaveI][0],
984 slaveMoveOldPoint[slaveI][0]
985 );
987 cMap.pushOperation
988 (
989 sEdge[1],
990 coupleMap::MOVE_POINT,
991 slaveMoveNewPoint[slaveI][1],
992 slaveMoveOldPoint[slaveI][1]
993 );
995 // Force operation to succeed
996 slaveMap.type() = 1;
997 }
998 else
999 {
1000 // Face-based coupling
1001 slaveMap =
1002 (
1003 sMesh.collapseQuadFace
1004 (
1005 sIndex,
1006 slaveOverRide,
1007 false,
1008 forceOp
1009 )
1010 );
1012 // Push operation into coupleMap
1013 switch (slaveMap.type())
1014 {
1015 case 1:
1016 {
1017 cMap.pushOperation
1018 (
1019 sIndex,
1020 coupleMap::COLLAPSE_FIRST
1021 );
1023 break;
1024 }
1026 case 2:
1027 {
1028 cMap.pushOperation
1029 (
1030 sIndex,
1031 coupleMap::COLLAPSE_SECOND
1032 );
1034 break;
1035 }
1037 case 3:
1038 {
1039 cMap.pushOperation
1040 (
1041 sIndex,
1042 coupleMap::COLLAPSE_MIDPOINT
1043 );
1045 break;
1046 }
1047 }
1048 }
1049 }
1051 // Turn it back on.
1052 setCoupledModification();
1054 // The final operation has to succeed.
1055 if (slaveMap.type() <= 0)
1056 {
1057 FatalErrorIn
1058 (
1059 "\n"
1060 "const changeMap "
1061 "dynamicTopoFvMesh::collapseQuadFace\n"
1062 "(\n"
1063 " const label fIndex,\n"
1064 " label overRideCase,\n"
1065 " bool checkOnly,\n"
1066 " bool forceOp\n"
1067 ")\n"
1068 )
1069 << "Coupled topo-change for slave failed." << nl
1070 << " Face: " << fIndex << ": " << fCheck << nl
1071 << " Slave index: " << sIndex << nl
1072 << " Patch index: " << pI << nl
1073 << " Type: " << slaveMap.type() << nl
1074 << abort(FatalError);
1075 }
1077 // Save index and patch for posterity
1078 slaveMap.index() = sIndex;
1079 slaveMap.patchIndex() = pI;
1080 }
1081 }
1083 // Build a hull of cells and tri-faces that are connected to each edge
1084 FixedList<labelList, 2> hullCells;
1085 FixedList<labelList, 2> hullTriFaces;
1087 meshOps::constructPrismHull
1088 (
1089 checkEdgeIndex[1],
1090 faces_,
1091 cells_,
1092 owner_,
1093 neighbour_,
1094 edgeFaces_,
1095 hullTriFaces[0],
1096 hullCells[0]
1097 );
1099 meshOps::constructPrismHull
1100 (
1101 checkEdgeIndex[2],
1102 faces_,
1103 cells_,
1104 owner_,
1105 neighbour_,
1106 edgeFaces_,
1107 hullTriFaces[1],
1108 hullCells[1]
1109 );
1111 // Determine the neighbouring cells
1112 label c0 = owner_[fIndex], c1 = neighbour_[fIndex];
1114 // Define variables for the prism-face calculation
1115 FixedList<label,2> c0BdyIndex(-1), c0IntIndex(-1);
1116 FixedList<label,2> c1BdyIndex(-1), c1IntIndex(-1);
1117 FixedList<face,2> c0BdyFace(face(3)), c0IntFace(face(4));
1118 FixedList<face,2> c1BdyFace(face(3)), c1IntFace(face(4));
1120 // Find the prism-faces
1121 meshOps::findPrismFaces
1122 (
1123 fIndex,
1124 c0,
1125 faces_,
1126 cells_,
1127 neighbour_,
1128 c0BdyFace,
1129 c0BdyIndex,
1130 c0IntFace,
1131 c0IntIndex
1132 );
1134 if (c1 != -1)
1135 {
1136 meshOps::findPrismFaces
1137 (
1138 fIndex,
1139 c1,
1140 faces_,
1141 cells_,
1142 neighbour_,
1143 c1BdyFace,
1144 c1BdyIndex,
1145 c1IntFace,
1146 c1IntIndex
1147 );
1148 }
1150 // Determine if either edge belongs to a boundary
1151 FixedList<bool, 2> nBoundCurves(false), edgeBoundary(false);
1153 edgeBoundary[0] = (whichEdgePatch(checkEdgeIndex[1]) > -1);
1154 edgeBoundary[1] = (whichEdgePatch(checkEdgeIndex[2]) > -1);
1156 // Decide on collapseCase
1157 label collapseCase = -1;
1159 if (edgeBoundary[0] && !edgeBoundary[1])
1160 {
1161 collapseCase = 1;
1162 }
1163 else
1164 if (!edgeBoundary[0] && edgeBoundary[1])
1165 {
1166 collapseCase = 2;
1167 }
1168 else
1169 if (edgeBoundary[0] && edgeBoundary[1])
1170 {
1171 if (c1 != -1)
1172 {
1173 if (debug > 2)
1174 {
1175 WarningIn
1176 (
1177 "\n"
1178 "const changeMap "
1179 "dynamicTopoFvMesh::collapseQuadFace\n"
1180 "(\n"
1181 " const label fIndex,\n"
1182 " label overRideCase,\n"
1183 " bool checkOnly,\n"
1184 " bool forceOp\n"
1185 ")\n"
1186 ) << "Collapsing an internal face that "
1187 << "lies on two boundary patches. "
1188 << "Face: " << fIndex << ": " << faces_[fIndex]
1189 << endl;
1190 }
1192 // Bail out for now. If proximity based refinement is
1193 // switched on, mesh may be sliced at this point.
1194 return map;
1195 }
1197 // Check if either edge lies on a bounding curve.
1198 if (checkBoundingCurve(checkEdgeIndex[1]))
1199 {
1200 nBoundCurves[0] = true;
1201 }
1203 if (checkBoundingCurve(checkEdgeIndex[2]))
1204 {
1205 nBoundCurves[1] = true;
1206 }
1208 // Collapse towards a bounding curve
1209 if (nBoundCurves[0] && !nBoundCurves[1])
1210 {
1211 collapseCase = 1;
1212 }
1213 else
1214 if (!nBoundCurves[0] && nBoundCurves[1])
1215 {
1216 collapseCase = 2;
1217 }
1218 else
1219 if (!nBoundCurves[0] && !nBoundCurves[1])
1220 {
1221 // No bounding curves. Collapse to mid-point.
1222 collapseCase = 3;
1223 }
1224 else
1225 {
1226 // Two bounding curves? This might cause pinching.
1227 // Bail out for now.
1228 return map;
1229 }
1230 }
1231 else
1232 {
1233 // Looks like this is an interior face.
1234 // Collapse case [3] by default
1235 collapseCase = 3;
1236 }
1238 // Perform an override if requested.
1239 if (overRideCase != -1)
1240 {
1241 collapseCase = overRideCase;
1242 }
1244 // Configure the new point-positions
1245 FixedList<bool, 2> check(false);
1246 FixedList<FixedList<label, 2>, 2> checkPoints(FixedList<label, 2>(-1));
1247 FixedList<point, 2> newPoint(vector::zero);
1248 FixedList<point, 2> oldPoint(vector::zero);
1250 // Replacement check points
1251 FixedList<label,2> original(-1), replacement(-1);
1253 switch (collapseCase)
1254 {
1255 case 1:
1256 {
1257 // Collapse to the first node
1258 original[0] = cv2; original[1] = cv3;
1259 replacement[0] = cv0; replacement[1] = cv1;
1261 newPoint[0] = points_[replacement[0]];
1262 newPoint[1] = points_[replacement[1]];
1264 oldPoint[0] = oldPoints_[replacement[0]];
1265 oldPoint[1] = oldPoints_[replacement[1]];
1267 // Define check-points
1268 check[1] = true;
1269 checkPoints[1][0] = original[0];
1270 checkPoints[1][1] = original[1];
1272 break;
1273 }
1275 case 2:
1276 {
1277 // Collapse to the second node
1278 original[0] = cv0; original[1] = cv1;
1279 replacement[0] = cv2; replacement[1] = cv3;
1281 newPoint[0] = points_[replacement[0]];
1282 newPoint[1] = points_[replacement[1]];
1284 oldPoint[0] = oldPoints_[replacement[0]];
1285 oldPoint[1] = oldPoints_[replacement[1]];
1287 // Define check-points
1288 check[0] = true;
1289 checkPoints[0][0] = original[0];
1290 checkPoints[0][1] = original[1];
1292 break;
1293 }
1295 case 3:
1296 {
1297 // Collapse to the mid-point
1298 original[0] = cv0; original[1] = cv1;
1299 replacement[0] = cv2; replacement[1] = cv3;
1301 // Define new point-positions
1302 newPoint[0] =
1303 (
1304 0.5 *
1305 (
1306 points_[original[0]]
1307 + points_[replacement[0]]
1308 )
1309 );
1311 newPoint[1] =
1312 (
1313 0.5 *
1314 (
1315 points_[original[1]]
1316 + points_[replacement[1]]
1317 )
1318 );
1320 // Define old point-positions
1321 oldPoint[0] =
1322 (
1323 0.5 *
1324 (
1325 oldPoints_[original[0]]
1326 + oldPoints_[replacement[0]]
1327 )
1328 );
1330 oldPoint[1] =
1331 (
1332 0.5 *
1333 (
1334 oldPoints_[original[1]]
1335 + oldPoints_[replacement[1]]
1336 )
1337 );
1339 // Define check-points
1340 check[0] = true;
1341 checkPoints[0][0] = original[0];
1342 checkPoints[0][1] = original[1];
1344 check[1] = true;
1345 checkPoints[1][0] = replacement[0];
1346 checkPoints[1][1] = replacement[1];
1348 break;
1349 }
1351 default:
1352 {
1353 // Don't think this will ever happen.
1354 FatalErrorIn
1355 (
1356 "\n"
1357 "const changeMap "
1358 "dynamicTopoFvMesh::collapseQuadFace\n"
1359 "(\n"
1360 " const label fIndex,\n"
1361 " label overRideCase,\n"
1362 " bool checkOnly,\n"
1363 " bool forceOp\n"
1364 ")\n"
1365 )
1366 << "Edge: " << fIndex << ": " << faces_[fIndex]
1367 << ". Couldn't decide on collapseCase."
1368 << abort(FatalError);
1370 break;
1371 }
1372 }
1374 // Keep track of resulting cell quality,
1375 // if collapse is indeed feasible
1376 scalar collapseQuality(GREAT);
1378 // Check collapsibility of cells around edges
1379 // with the re-configured point
1380 forAll(check, indexI)
1381 {
1382 if (!check[indexI])
1383 {
1384 continue;
1385 }
1387 // Check whether the collapse is possible.
1388 if
1389 (
1390 checkCollapse
1391 (
1392 (*this),
1393 hullTriFaces[indexI],
1394 c0BdyIndex,
1395 c1BdyIndex,
1396 checkPoints[indexI],
1397 newPoint,
1398 oldPoint,
1399 collapseQuality,
1400 (c1 != -1),
1401 forceOp
1402 )
1403 )
1404 {
1405 map.type() = 0;
1406 return map;
1407 }
1408 }
1410 // Add a map entry of the replacements as 'addedPoints'
1411 // - Used in coupled mapping
1412 map.addPoint(replacement[0]);
1413 map.addPoint(replacement[1]);
1415 // Are we only performing checks?
1416 if (checkOnly)
1417 {
1418 map.type() = collapseCase;
1420 if (debug > 2)
1421 {
1422 Pout<< "Face: " << fIndex
1423 << ":: " << faces_[fIndex] << nl
1424 << " collapseCase determined to be: " << collapseCase << nl
1425 << " Resulting quality: " << collapseQuality << nl
1426 << " edgeBoundary: " << edgeBoundary << nl
1427 << " nBoundCurves: " << nBoundCurves << nl
1428 << " collapsePoints: " << original << nl
1429 << " replacePoints: " << replacement << nl
1430 << endl;
1431 }
1433 return map;
1434 }
1436 if (debug > 1)
1437 {
1438 const labelList& fE = faceEdges_[fIndex];
1440 Pout<< nl << nl
1441 << "Face: " << fIndex << ": " << faces_[fIndex] << nl
1442 << "faceEdges: " << fE << " is to be collapsed. "
1443 << nl;
1445 Pout<< " On SubMesh: " << isSubMesh_ << nl;
1446 Pout<< " coupledModification: " << coupledModification_ << nl;
1448 label epIndex = whichPatch(fIndex);
1450 const polyBoundaryMesh& boundary = boundaryMesh();
1452 Pout<< " Patch: ";
1454 if (epIndex == -1)
1455 {
1456 Pout<< "Internal" << nl;
1457 }
1458 else
1459 if (epIndex < boundary.size())
1460 {
1461 Pout<< boundary[epIndex].name() << nl;
1462 }
1463 else
1464 {
1465 Pout<< " New patch: " << epIndex << endl;
1466 }
1468 if (debug > 2)
1469 {
1470 Pout<< nl
1471 << "~~~~~~~~~~~~~~~~~~~~~~~~~" << nl
1472 << "Hulls before modification" << nl
1473 << "~~~~~~~~~~~~~~~~~~~~~~~~~" << nl;
1475 if (debug > 3)
1476 {
1477 Pout<< "Cell [0] Boundary faces: " << nl
1478 << " Face: "<< c0BdyIndex[0]
1479 << "::" << c0BdyFace[0] << nl
1480 << " Face: "<< c0BdyIndex[1]
1481 << "::" << c0BdyFace[1] << nl;
1483 Pout<< "Cell [0] Interior faces: " << nl
1484 << " Face: "<< c0IntIndex[0]
1485 << "::" << c0IntFace[0] << nl
1486 << " Face: "<< c0IntIndex[1]
1487 << "::" << c0IntFace[1] << nl;
1489 if (c1 != -1)
1490 {
1491 Pout<< "Cell [1] Boundary faces: " << nl
1492 << " Face: "<< c1BdyIndex[0]
1493 << "::" << c1BdyFace[0] << nl
1494 << " Face: "<< c1BdyIndex[1]
1495 << "::" << c1BdyFace[1] << nl;
1497 Pout<< "Cell [1] Interior faces: " << nl
1498 << " Face: "<< c1IntIndex[0]
1499 << "::" << c1IntFace[0] << nl
1500 << " Face: "<< c1IntIndex[1]
1501 << "::" << c1IntFace[1] << nl;
1502 }
1503 }
1505 Pout<< nl << "Cells belonging to first Edge Hull: "
1506 << hullCells[0] << nl;
1508 forAll(hullCells[0],cellI)
1509 {
1510 const cell& firstCurCell = cells_[hullCells[0][cellI]];
1512 Pout<< "Cell: " << hullCells[0][cellI]
1513 << ": " << firstCurCell << nl;
1515 forAll(firstCurCell,faceI)
1516 {
1517 Pout<< " Face: " << firstCurCell[faceI]
1518 << " : " << faces_[firstCurCell[faceI]]
1519 << " fE: " << faceEdges_[firstCurCell[faceI]]
1520 << nl;
1522 if (debug > 3)
1523 {
1524 const labelList& fE = faceEdges_[firstCurCell[faceI]];
1526 forAll(fE, edgeI)
1527 {
1528 Pout<< " Edge: " << fE[edgeI]
1529 << " : " << edges_[fE[edgeI]]
1530 << nl;
1531 }
1532 }
1533 }
1534 }
1536 const labelList& firstEdgeFaces = edgeFaces_[checkEdgeIndex[1]];
1538 Pout<< nl << "First Edge Face Hull: "
1539 << firstEdgeFaces << nl;
1541 forAll(firstEdgeFaces,indexI)
1542 {
1543 Pout<< " Face: " << firstEdgeFaces[indexI]
1544 << " : " << faces_[firstEdgeFaces[indexI]]
1545 << " fE: " << faceEdges_[firstEdgeFaces[indexI]]
1546 << nl;
1548 if (debug > 3)
1549 {
1550 const labelList& fE = faceEdges_[firstEdgeFaces[indexI]];
1552 forAll(fE, edgeI)
1553 {
1554 Pout<< " Edge: " << fE[edgeI]
1555 << " : " << edges_[fE[edgeI]]
1556 << nl;
1557 }
1558 }
1559 }
1561 Pout<< nl << "Cells belonging to second Edge Hull: "
1562 << hullCells[1] << nl;
1564 forAll(hullCells[1], cellI)
1565 {
1566 const cell& secondCurCell = cells_[hullCells[1][cellI]];
1568 Pout<< "Cell: " << hullCells[1][cellI]
1569 << ": " << secondCurCell << nl;
1571 forAll(secondCurCell, faceI)
1572 {
1573 Pout<< " Face: " << secondCurCell[faceI]
1574 << " : " << faces_[secondCurCell[faceI]]
1575 << " fE: " << faceEdges_[secondCurCell[faceI]]
1576 << nl;
1578 if (debug > 3)
1579 {
1580 const labelList& fE = faceEdges_[secondCurCell[faceI]];
1582 forAll(fE, edgeI)
1583 {
1584 Pout<< " Edge: " << fE[edgeI]
1585 << " : " << edges_[fE[edgeI]]
1586 << nl;
1587 }
1588 }
1589 }
1590 }
1592 const labelList& secondEdgeFaces = edgeFaces_[checkEdgeIndex[2]];
1594 Pout<< nl << "Second Edge Face Hull: "
1595 << secondEdgeFaces << nl;
1597 forAll(secondEdgeFaces, indexI)
1598 {
1599 Pout<< " Face: " << secondEdgeFaces[indexI]
1600 << " : " << faces_[secondEdgeFaces[indexI]]
1601 << " fE: " << faceEdges_[secondEdgeFaces[indexI]]
1602 << nl;
1604 if (debug > 3)
1605 {
1606 const labelList& fE = faceEdges_[secondEdgeFaces[indexI]];
1608 forAll(fE, edgeI)
1609 {
1610 Pout<< " Edge: " << fE[edgeI]
1611 << " : " << edges_[fE[edgeI]]
1612 << nl;
1613 }
1614 }
1615 }
1617 // Write out VTK files prior to change
1618 if (debug > 3)
1619 {
1620 labelHashSet vtkCells;
1622 forAll(hullCells[0], cellI)
1623 {
1624 if (!vtkCells.found(hullCells[0][cellI]))
1625 {
1626 vtkCells.insert(hullCells[0][cellI]);
1627 }
1628 }
1630 forAll(hullCells[1], cellI)
1631 {
1632 if (!vtkCells.found(hullCells[1][cellI]))
1633 {
1634 vtkCells.insert(hullCells[1][cellI]);
1635 }
1636 }
1638 writeVTK
1639 (
1640 Foam::name(fIndex)
1641 + "_Collapse_0",
1642 vtkCells.toc(),
1643 3, false, true
1644 );
1645 }
1646 }
1647 }
1649 // Edges / Quad-faces to throw or keep during collapse
1650 FixedList<label,2> ends(-1);
1651 FixedList<label,2> faceToKeep(-1), faceToThrow(-1);
1652 FixedList<label,4> edgeToKeep(-1), edgeToThrow(-1);
1654 // Maintain a list of modified faces for mapping
1655 dynamicLabelList modifiedFaces(10);
1657 // Case 2 & 3 use identical connectivity,
1658 // but different point locations
1659 if (collapseCase == 2 || collapseCase == 3)
1660 {
1661 const labelList& firstEdgeFaces = edgeFaces_[checkEdgeIndex[1]];
1663 // Collapse to the second node...
1664 forAll(firstEdgeFaces,faceI)
1665 {
1666 // Replace point indices on faces.
1667 meshOps::replaceLabel
1668 (
1669 cv0,
1670 cv2,
1671 faces_[firstEdgeFaces[faceI]]
1672 );
1674 meshOps::replaceLabel
1675 (
1676 cv1,
1677 cv3,
1678 faces_[firstEdgeFaces[faceI]]
1679 );
1681 // Add an entry for mapping
1682 if (findIndex(modifiedFaces, firstEdgeFaces[faceI]) == -1)
1683 {
1684 modifiedFaces.append(firstEdgeFaces[faceI]);
1685 }
1687 // Determine the quad-face in cell[0] & cell[1]
1688 // that belongs to firstEdgeFaces
1689 if (firstEdgeFaces[faceI] == c0IntIndex[0])
1690 {
1691 faceToKeep[0] = c0IntIndex[1];
1692 faceToThrow[0] = c0IntIndex[0];
1693 }
1695 if (firstEdgeFaces[faceI] == c0IntIndex[1])
1696 {
1697 faceToKeep[0] = c0IntIndex[0];
1698 faceToThrow[0] = c0IntIndex[1];
1699 }
1701 if (c1 != -1)
1702 {
1703 if (firstEdgeFaces[faceI] == c1IntIndex[0])
1704 {
1705 faceToKeep[1] = c1IntIndex[1];
1706 faceToThrow[1] = c1IntIndex[0];
1707 }
1709 if (firstEdgeFaces[faceI] == c1IntIndex[1])
1710 {
1711 faceToKeep[1] = c1IntIndex[0];
1712 faceToThrow[1] = c1IntIndex[1];
1713 }
1714 }
1715 }
1717 // Find common edges between quad / tri faces...
1718 meshOps::findCommonEdge
1719 (
1720 c0BdyIndex[0],
1721 faceToKeep[0],
1722 faceEdges_,
1723 edgeToKeep[0]
1724 );
1726 meshOps::findCommonEdge
1727 (
1728 c0BdyIndex[1],
1729 faceToKeep[0],
1730 faceEdges_,
1731 edgeToKeep[1]
1732 );
1734 meshOps::findCommonEdge
1735 (
1736 c0BdyIndex[0],
1737 faceToThrow[0],
1738 faceEdges_,
1739 edgeToThrow[0]
1740 );
1742 meshOps::findCommonEdge
1743 (
1744 c0BdyIndex[1],
1745 faceToThrow[0],
1746 faceEdges_,
1747 edgeToThrow[1]
1748 );
1750 // Size down edgeFaces for the ends.
1751 meshOps::findCommonEdge
1752 (
1753 faceToThrow[0],
1754 faceToKeep[0],
1755 faceEdges_,
1756 ends[0]
1757 );
1759 meshOps::sizeDownList
1760 (
1761 faceToThrow[0],
1762 edgeFaces_[ends[0]]
1763 );
1765 if (c1 != -1)
1766 {
1767 meshOps::findCommonEdge
1768 (
1769 c1BdyIndex[0],
1770 faceToKeep[1],
1771 faceEdges_,
1772 edgeToKeep[2]
1773 );
1775 meshOps::findCommonEdge
1776 (
1777 c1BdyIndex[1],
1778 faceToKeep[1],
1779 faceEdges_,
1780 edgeToKeep[3]
1781 );
1783 meshOps::findCommonEdge
1784 (
1785 c1BdyIndex[0],
1786 faceToThrow[1],
1787 faceEdges_,
1788 edgeToThrow[2]
1789 );
1791 meshOps::findCommonEdge
1792 (
1793 c1BdyIndex[1],
1794 faceToThrow[1],
1795 faceEdges_,
1796 edgeToThrow[3]
1797 );
1799 // Size down edgeFaces for the ends.
1800 meshOps::findCommonEdge
1801 (
1802 faceToThrow[1],
1803 faceToKeep[1],
1804 faceEdges_,
1805 ends[1]
1806 );
1808 meshOps::sizeDownList
1809 (
1810 faceToThrow[1],
1811 edgeFaces_[ends[1]]
1812 );
1813 }
1815 // Correct edgeFaces for triangular faces...
1816 forAll(edgeToThrow, indexI)
1817 {
1818 if (edgeToThrow[indexI] == -1)
1819 {
1820 continue;
1821 }
1823 const labelList& eF = edgeFaces_[edgeToThrow[indexI]];
1825 label origTriFace = -1, retTriFace = -1;
1827 // Find the original triangular face index
1828 forAll(eF, faceI)
1829 {
1830 if (eF[faceI] == c0BdyIndex[0])
1831 {
1832 origTriFace = c0BdyIndex[0];
1833 break;
1834 }
1836 if (eF[faceI] == c0BdyIndex[1])
1837 {
1838 origTriFace = c0BdyIndex[1];
1839 break;
1840 }
1842 if (eF[faceI] == c1BdyIndex[0])
1843 {
1844 origTriFace = c1BdyIndex[0];
1845 break;
1846 }
1848 if (eF[faceI] == c1BdyIndex[1])
1849 {
1850 origTriFace = c1BdyIndex[1];
1851 break;
1852 }
1853 }
1855 // Find the retained triangular face index
1856 forAll(eF, faceI)
1857 {
1858 if
1859 (
1860 (eF[faceI] != origTriFace) &&
1861 (eF[faceI] != faceToThrow[0]) &&
1862 (eF[faceI] != faceToThrow[1])
1863 )
1864 {
1865 retTriFace = eF[faceI];
1866 break;
1867 }
1868 }
1870 // Finally replace the face index
1871 if (retTriFace == -1)
1872 {
1873 // Couldn't find a retained face.
1874 // This must be a boundary edge, so size-down instead.
1875 meshOps::sizeDownList
1876 (
1877 origTriFace,
1878 edgeFaces_[edgeToKeep[indexI]]
1879 );
1880 }
1881 else
1882 {
1883 meshOps::replaceLabel
1884 (
1885 origTriFace,
1886 retTriFace,
1887 edgeFaces_[edgeToKeep[indexI]]
1888 );
1890 meshOps::replaceLabel
1891 (
1892 edgeToThrow[indexI],
1893 edgeToKeep[indexI],
1894 faceEdges_[retTriFace]
1895 );
1896 }
1897 }
1899 // Correct faceEdges / edgeFaces for quad-faces...
1900 forAll(firstEdgeFaces,faceI)
1901 {
1902 meshOps::replaceLabel
1903 (
1904 checkEdgeIndex[1],
1905 checkEdgeIndex[2],
1906 faceEdges_[firstEdgeFaces[faceI]]
1907 );
1909 // Renumber the edges on this face
1910 const labelList& fE = faceEdges_[firstEdgeFaces[faceI]];
1912 forAll(fE, edgeI)
1913 {
1914 if (edges_[fE[edgeI]].start() == cv0)
1915 {
1916 edges_[fE[edgeI]].start() = cv2;
1917 }
1919 if (edges_[fE[edgeI]].end() == cv0)
1920 {
1921 edges_[fE[edgeI]].end() = cv2;
1922 }
1924 if (edges_[fE[edgeI]].start() == cv1)
1925 {
1926 edges_[fE[edgeI]].start() = cv3;
1927 }
1929 if (edges_[fE[edgeI]].end() == cv1)
1930 {
1931 edges_[fE[edgeI]].end() = cv3;
1932 }
1933 }
1935 // Size-up edgeFaces for the replacement
1936 if
1937 (
1938 (firstEdgeFaces[faceI] != faceToThrow[0]) &&
1939 (firstEdgeFaces[faceI] != faceToThrow[1]) &&
1940 (firstEdgeFaces[faceI] != fIndex)
1941 )
1942 {
1943 meshOps::sizeUpList
1944 (
1945 firstEdgeFaces[faceI],
1946 edgeFaces_[checkEdgeIndex[2]]
1947 );
1948 }
1949 }
1951 // Remove the current face from the replacement edge
1952 meshOps::sizeDownList
1953 (
1954 fIndex,
1955 edgeFaces_[checkEdgeIndex[2]]
1956 );
1958 // Replace point labels on all triangular boundary faces.
1959 forAll(hullCells[0],cellI)
1960 {
1961 const cell& cellToCheck = cells_[hullCells[0][cellI]];
1963 forAll(cellToCheck,faceI)
1964 {
1965 face& faceToCheck = faces_[cellToCheck[faceI]];
1967 if (faceToCheck.size() == 3)
1968 {
1969 forAll(faceToCheck,pointI)
1970 {
1971 if (faceToCheck[pointI] == cv0)
1972 {
1973 faceToCheck[pointI] = cv2;
1974 }
1976 if (faceToCheck[pointI] == cv1)
1977 {
1978 faceToCheck[pointI] = cv3;
1979 }
1980 }
1981 }
1982 }
1983 }
1985 // Now that we're done with all edges, remove them.
1986 removeEdge(checkEdgeIndex[0]);
1987 removeEdge(checkEdgeIndex[1]);
1988 removeEdge(checkEdgeIndex[3]);
1990 // Update map
1991 map.removeEdge(checkEdgeIndex[0]);
1992 map.removeEdge(checkEdgeIndex[1]);
1993 map.removeEdge(checkEdgeIndex[2]);
1995 forAll(edgeToThrow, indexI)
1996 {
1997 if (edgeToThrow[indexI] == -1)
1998 {
1999 continue;
2000 }
2002 removeEdge(edgeToThrow[indexI]);
2004 // Update map
2005 map.removeEdge(edgeToThrow[indexI]);
2006 }
2008 // Delete the two points...
2009 removePoint(cv0);
2010 removePoint(cv1);
2012 // Update map
2013 map.removePoint(cv0);
2014 map.removePoint(cv1);
2015 }
2016 else
2017 {
2018 const labelList& secondEdgeFaces = edgeFaces_[checkEdgeIndex[2]];
2020 // Collapse to the first node
2021 forAll(secondEdgeFaces,faceI)
2022 {
2023 // Replace point indices on faces.
2024 meshOps::replaceLabel
2025 (
2026 cv2,
2027 cv0,
2028 faces_[secondEdgeFaces[faceI]]
2029 );
2031 meshOps::replaceLabel
2032 (
2033 cv3,
2034 cv1,
2035 faces_[secondEdgeFaces[faceI]]
2036 );
2038 // Add an entry for mapping
2039 if (findIndex(modifiedFaces, secondEdgeFaces[faceI]) == -1)
2040 {
2041 modifiedFaces.append(secondEdgeFaces[faceI]);
2042 }
2044 // Determine the quad-face(s) in cell[0] & cell[1]
2045 // that belongs to secondEdgeFaces
2046 if (secondEdgeFaces[faceI] == c0IntIndex[0])
2047 {
2048 faceToKeep[0] = c0IntIndex[1];
2049 faceToThrow[0] = c0IntIndex[0];
2050 }
2052 if (secondEdgeFaces[faceI] == c0IntIndex[1])
2053 {
2054 faceToKeep[0] = c0IntIndex[0];
2055 faceToThrow[0] = c0IntIndex[1];
2056 }
2058 if (c1 != -1)
2059 {
2060 if (secondEdgeFaces[faceI] == c1IntIndex[0])
2061 {
2062 faceToKeep[1] = c1IntIndex[1];
2063 faceToThrow[1] = c1IntIndex[0];
2064 }
2066 if (secondEdgeFaces[faceI] == c1IntIndex[1])
2067 {
2068 faceToKeep[1] = c1IntIndex[0];
2069 faceToThrow[1] = c1IntIndex[1];
2070 }
2071 }
2072 }
2074 // Find common edges between quad / tri faces...
2075 meshOps::findCommonEdge
2076 (
2077 c0BdyIndex[0],
2078 faceToKeep[0],
2079 faceEdges_,
2080 edgeToKeep[0]
2081 );
2083 meshOps::findCommonEdge
2084 (
2085 c0BdyIndex[1],
2086 faceToKeep[0],
2087 faceEdges_,
2088 edgeToKeep[1]
2089 );
2091 meshOps::findCommonEdge
2092 (
2093 c0BdyIndex[0],
2094 faceToThrow[0],
2095 faceEdges_,
2096 edgeToThrow[0]
2097 );
2099 meshOps::findCommonEdge
2100 (
2101 c0BdyIndex[1],
2102 faceToThrow[0],
2103 faceEdges_,
2104 edgeToThrow[1]
2105 );
2107 // Size down edgeFaces for the ends.
2108 meshOps::findCommonEdge
2109 (
2110 faceToThrow[0],
2111 faceToKeep[0],
2112 faceEdges_,
2113 ends[0]
2114 );
2116 meshOps::sizeDownList
2117 (
2118 faceToThrow[0],
2119 edgeFaces_[ends[0]]
2120 );
2122 if (c1 != -1)
2123 {
2124 meshOps::findCommonEdge
2125 (
2126 c1BdyIndex[0],
2127 faceToKeep[1],
2128 faceEdges_,
2129 edgeToKeep[2]
2130 );
2132 meshOps::findCommonEdge
2133 (
2134 c1BdyIndex[1],
2135 faceToKeep[1],
2136 faceEdges_,
2137 edgeToKeep[3]
2138 );
2140 meshOps::findCommonEdge
2141 (
2142 c1BdyIndex[0],
2143 faceToThrow[1],
2144 faceEdges_,
2145 edgeToThrow[2]
2146 );
2148 meshOps::findCommonEdge
2149 (
2150 c1BdyIndex[1],
2151 faceToThrow[1],
2152 faceEdges_,
2153 edgeToThrow[3]
2154 );
2156 // Size down edgeFaces for the ends.
2157 meshOps::findCommonEdge
2158 (
2159 faceToThrow[1],
2160 faceToKeep[1],
2161 faceEdges_,
2162 ends[1]
2163 );
2165 meshOps::sizeDownList
2166 (
2167 faceToThrow[1],
2168 edgeFaces_[ends[1]]
2169 );
2170 }
2172 // Correct edgeFaces for triangular faces...
2173 forAll(edgeToThrow, indexI)
2174 {
2175 if (edgeToThrow[indexI] == -1)
2176 {
2177 continue;
2178 }
2180 const labelList& eF = edgeFaces_[edgeToThrow[indexI]];
2182 label origTriFace = -1, retTriFace = -1;
2184 // Find the original triangular face index
2185 forAll(eF, faceI)
2186 {
2187 if (eF[faceI] == c0BdyIndex[0])
2188 {
2189 origTriFace = c0BdyIndex[0];
2190 break;
2191 }
2193 if (eF[faceI] == c0BdyIndex[1])
2194 {
2195 origTriFace = c0BdyIndex[1];
2196 break;
2197 }
2199 if (eF[faceI] == c1BdyIndex[0])
2200 {
2201 origTriFace = c1BdyIndex[0];
2202 break;
2203 }
2205 if (eF[faceI] == c1BdyIndex[1])
2206 {
2207 origTriFace = c1BdyIndex[1];
2208 break;
2209 }
2210 }
2212 // Find the retained triangular face index
2213 forAll(eF, faceI)
2214 {
2215 if
2216 (
2217 (eF[faceI] != origTriFace) &&
2218 (eF[faceI] != faceToThrow[0]) &&
2219 (eF[faceI] != faceToThrow[1])
2220 )
2221 {
2222 retTriFace = eF[faceI];
2223 break;
2224 }
2225 }
2227 // Finally replace the face/edge indices
2228 if (retTriFace == -1)
2229 {
2230 // Couldn't find a retained face.
2231 // This must be a boundary edge, so size-down instead.
2232 meshOps::sizeDownList
2233 (
2234 origTriFace,
2235 edgeFaces_[edgeToKeep[indexI]]
2236 );
2237 }
2238 else
2239 {
2240 meshOps::replaceLabel
2241 (
2242 origTriFace,
2243 retTriFace,
2244 edgeFaces_[edgeToKeep[indexI]]
2245 );
2247 meshOps::replaceLabel
2248 (
2249 edgeToThrow[indexI],
2250 edgeToKeep[indexI],
2251 faceEdges_[retTriFace]
2252 );
2253 }
2254 }
2256 // Correct faceEdges / edgeFaces for quad-faces...
2257 forAll(secondEdgeFaces,faceI)
2258 {
2259 meshOps::replaceLabel
2260 (
2261 checkEdgeIndex[2],
2262 checkEdgeIndex[1],
2263 faceEdges_[secondEdgeFaces[faceI]]
2264 );
2266 // Renumber the edges on this face
2267 const labelList& fE = faceEdges_[secondEdgeFaces[faceI]];
2269 forAll(fE, edgeI)
2270 {
2271 if (edges_[fE[edgeI]].start() == cv2)
2272 {
2273 edges_[fE[edgeI]].start() = cv0;
2274 }
2276 if (edges_[fE[edgeI]].end() == cv2)
2277 {
2278 edges_[fE[edgeI]].end() = cv0;
2279 }
2281 if (edges_[fE[edgeI]].start() == cv3)
2282 {
2283 edges_[fE[edgeI]].start() = cv1;
2284 }
2286 if (edges_[fE[edgeI]].end() == cv3)
2287 {
2288 edges_[fE[edgeI]].end() = cv1;
2289 }
2290 }
2292 // Size-up edgeFaces for the replacement
2293 if
2294 (
2295 (secondEdgeFaces[faceI] != faceToThrow[0]) &&
2296 (secondEdgeFaces[faceI] != faceToThrow[1]) &&
2297 (secondEdgeFaces[faceI] != fIndex)
2298 )
2299 {
2300 meshOps::sizeUpList
2301 (
2302 secondEdgeFaces[faceI],
2303 edgeFaces_[checkEdgeIndex[1]]
2304 );
2305 }
2306 }
2308 // Remove the current face from the replacement edge
2309 meshOps::sizeDownList
2310 (
2311 fIndex,
2312 edgeFaces_[checkEdgeIndex[1]]
2313 );
2315 // Replace point labels on all triangular boundary faces.
2316 forAll(hullCells[1], cellI)
2317 {
2318 const cell& cellToCheck = cells_[hullCells[1][cellI]];
2320 forAll(cellToCheck, faceI)
2321 {
2322 face& faceToCheck = faces_[cellToCheck[faceI]];
2324 if (faceToCheck.size() == 3)
2325 {
2326 forAll(faceToCheck, pointI)
2327 {
2328 if (faceToCheck[pointI] == cv2)
2329 {
2330 faceToCheck[pointI] = cv0;
2331 }
2333 if (faceToCheck[pointI] == cv3)
2334 {
2335 faceToCheck[pointI] = cv1;
2336 }
2337 }
2338 }
2339 }
2340 }
2342 // Now that we're done with all edges, remove them.
2343 removeEdge(checkEdgeIndex[0]);
2344 removeEdge(checkEdgeIndex[2]);
2345 removeEdge(checkEdgeIndex[3]);
2347 // Update map
2348 map.removeEdge(checkEdgeIndex[0]);
2349 map.removeEdge(checkEdgeIndex[2]);
2350 map.removeEdge(checkEdgeIndex[3]);
2352 forAll(edgeToThrow, indexI)
2353 {
2354 if (edgeToThrow[indexI] == -1)
2355 {
2356 continue;
2357 }
2359 removeEdge(edgeToThrow[indexI]);
2361 // Update map
2362 map.removeEdge(edgeToThrow[indexI]);
2363 }
2365 // Delete the two points...
2366 removePoint(cv2);
2367 removePoint(cv3);
2369 // Update map
2370 map.removePoint(cv2);
2371 map.removePoint(cv3);
2372 }
2374 if (debug > 2)
2375 {
2376 Pout<< nl
2377 << "~~~~~~~~~~~~~~~~~~~~~~~~" << nl
2378 << "Hulls after modification" << nl
2379 << "~~~~~~~~~~~~~~~~~~~~~~~~" << nl;
2381 Pout<< nl << "Cells belonging to first Edge Hull: "
2382 << hullCells[0] << nl;
2384 forAll(hullCells[0], cellI)
2385 {
2386 const cell& firstCurCell = cells_[hullCells[0][cellI]];
2388 Pout<< "Cell: " << hullCells[0][cellI]
2389 << ": " << firstCurCell << nl;
2391 forAll(firstCurCell, faceI)
2392 {
2393 Pout<< " Face: " << firstCurCell[faceI]
2394 << " : " << faces_[firstCurCell[faceI]]
2395 << " fE: " << faceEdges_[firstCurCell[faceI]]
2396 << nl;
2397 }
2398 }
2400 const labelList& firstEdgeFaces = edgeFaces_[checkEdgeIndex[1]];
2402 Pout<< nl << "First Edge Face Hull: " << firstEdgeFaces << nl;
2404 forAll(firstEdgeFaces, indexI)
2405 {
2406 Pout<< " Face: " << firstEdgeFaces[indexI]
2407 << " : " << faces_[firstEdgeFaces[indexI]]
2408 << " fE: " << faceEdges_[firstEdgeFaces[indexI]]
2409 << nl;
2410 }
2412 Pout<< nl << "Cells belonging to second Edge Hull: "
2413 << hullCells[1] << nl;
2415 forAll(hullCells[1], cellI)
2416 {
2417 const cell& secondCurCell = cells_[hullCells[1][cellI]];
2419 Pout<< "Cell: " << hullCells[1][cellI]
2420 << ": " << secondCurCell << nl;
2422 forAll(secondCurCell, faceI)
2423 {
2424 Pout<< " Face: " << secondCurCell[faceI]
2425 << " : " << faces_[secondCurCell[faceI]]
2426 << " fE: " << faceEdges_[secondCurCell[faceI]]
2427 << nl;
2428 }
2429 }
2431 const labelList& secondEdgeFaces = edgeFaces_[checkEdgeIndex[2]];
2433 Pout<< nl << "Second Edge Face Hull: " << secondEdgeFaces << nl;
2435 forAll(secondEdgeFaces, indexI)
2436 {
2437 Pout<< " Face : " << secondEdgeFaces[indexI]
2438 << " : " << faces_[secondEdgeFaces[indexI]]
2439 << " fE: " << faceEdges_[secondEdgeFaces[indexI]]
2440 << nl;
2441 }
2443 Pout<< "Retained face: "
2444 << faceToKeep[0] << ": "
2445 << " owner: " << owner_[faceToKeep[0]]
2446 << " neighbour: " << neighbour_[faceToKeep[0]]
2447 << nl;
2449 Pout<< "Discarded face: "
2450 << faceToThrow[0] << ": "
2451 << " owner: " << owner_[faceToThrow[0]]
2452 << " neighbour: " << neighbour_[faceToThrow[0]]
2453 << nl;
2455 if (c1 != -1)
2456 {
2457 Pout<< "Retained face: "
2458 << faceToKeep[1] << ": "
2459 << " owner: " << owner_[faceToKeep[1]]
2460 << " neighbour: " << neighbour_[faceToKeep[1]]
2461 << nl;
2463 Pout<< "Discarded face: "
2464 << faceToThrow[1] << ": "
2465 << " owner: " << owner_[faceToThrow[1]]
2466 << " neighbour: " << neighbour_[faceToThrow[1]]
2467 << nl;
2468 }
2470 Pout<< endl;
2471 }
2473 // Ensure proper orientation for the two retained faces
2474 FixedList<label,2> cellCheck(0);
2476 if (owner_[faceToThrow[0]] == c0)
2477 {
2478 cellCheck[0] = neighbour_[faceToThrow[0]];
2480 if (owner_[faceToKeep[0]] == c0)
2481 {
2482 if
2483 (
2484 (neighbour_[faceToThrow[0]] > neighbour_[faceToKeep[0]])
2485 && (neighbour_[faceToKeep[0]] != -1)
2486 )
2487 {
2488 // This face is to be flipped
2489 faces_[faceToKeep[0]] =
2490 (
2491 faces_[faceToKeep[0]].reverseFace()
2492 );
2494 owner_[faceToKeep[0]] = neighbour_[faceToKeep[0]];
2495 neighbour_[faceToKeep[0]] = neighbour_[faceToThrow[0]];
2497 setFlip(faceToKeep[0]);
2498 }
2499 else
2500 {
2501 if (neighbour_[faceToThrow[0]] != -1)
2502 {
2503 // Keep orientation intact, and update the owner
2504 owner_[faceToKeep[0]] = neighbour_[faceToThrow[0]];
2505 }
2506 else
2507 {
2508 // This face will need to be flipped and converted
2509 // to a boundary face. Flip it now, so that conversion
2510 // happens later.
2511 faces_[faceToKeep[0]] =
2512 (
2513 faces_[faceToKeep[0]].reverseFace()
2514 );
2516 owner_[faceToKeep[0]] = neighbour_[faceToKeep[0]];
2517 neighbour_[faceToKeep[0]] = -1;
2519 setFlip(faceToKeep[0]);
2520 }
2521 }
2522 }
2523 else
2524 {
2525 // Keep orientation intact, and update the neighbour
2526 neighbour_[faceToKeep[0]] = neighbour_[faceToThrow[0]];
2527 }
2528 }
2529 else
2530 {
2531 cellCheck[0] = owner_[faceToThrow[0]];
2533 if (neighbour_[faceToKeep[0]] == c0)
2534 {
2535 if (owner_[faceToThrow[0]] < owner_[faceToKeep[0]])
2536 {
2537 // This face is to be flipped
2538 faces_[faceToKeep[0]] =
2539 (
2540 faces_[faceToKeep[0]].reverseFace()
2541 );
2543 neighbour_[faceToKeep[0]] = owner_[faceToKeep[0]];
2544 owner_[faceToKeep[0]] = owner_[faceToThrow[0]];
2546 setFlip(faceToKeep[0]);
2547 }
2548 else
2549 {
2550 // Keep orientation intact, and update the neighbour
2551 neighbour_[faceToKeep[0]] = owner_[faceToThrow[0]];
2552 }
2553 }
2554 else
2555 {
2556 // Keep orientation intact, and update the owner
2557 owner_[faceToKeep[0]] = owner_[faceToThrow[0]];
2558 }
2559 }
2561 if (c1 != -1)
2562 {
2563 if (owner_[faceToThrow[1]] == c1)
2564 {
2565 cellCheck[1] = neighbour_[faceToThrow[1]];
2567 if (owner_[faceToKeep[1]] == c1)
2568 {
2569 if
2570 (
2571 (neighbour_[faceToThrow[1]] > neighbour_[faceToKeep[1]])
2572 && (neighbour_[faceToKeep[1]] != -1)
2573 )
2574 {
2575 // This face is to be flipped
2576 faces_[faceToKeep[1]] =
2577 (
2578 faces_[faceToKeep[1]].reverseFace()
2579 );
2581 owner_[faceToKeep[1]] = neighbour_[faceToKeep[1]];
2582 neighbour_[faceToKeep[1]] = neighbour_[faceToThrow[1]];
2584 setFlip(faceToKeep[1]);
2585 }
2586 else
2587 {
2588 if (neighbour_[faceToThrow[1]] != -1)
2589 {
2590 // Keep orientation intact, and update the owner
2591 owner_[faceToKeep[1]] = neighbour_[faceToThrow[1]];
2592 }
2593 else
2594 {
2595 // This face will need to be flipped and converted
2596 // to a boundary face. Flip it now, so that conversion
2597 // happens later.
2598 faces_[faceToKeep[1]] =
2599 (
2600 faces_[faceToKeep[1]].reverseFace()
2601 );
2603 owner_[faceToKeep[1]] = neighbour_[faceToKeep[1]];
2604 neighbour_[faceToKeep[1]] = -1;
2606 setFlip(faceToKeep[1]);
2607 }
2608 }
2609 }
2610 else
2611 {
2612 // Keep orientation intact, and update the neighbour
2613 neighbour_[faceToKeep[1]] = neighbour_[faceToThrow[1]];
2614 }
2615 }
2616 else
2617 {
2618 cellCheck[1] = owner_[faceToThrow[1]];
2620 if (neighbour_[faceToKeep[1]] == c1)
2621 {
2622 if (owner_[faceToThrow[1]] < owner_[faceToKeep[1]])
2623 {
2624 // This face is to be flipped
2625 faces_[faceToKeep[1]] =
2626 (
2627 faces_[faceToKeep[1]].reverseFace()
2628 );
2630 neighbour_[faceToKeep[1]] = owner_[faceToKeep[1]];
2631 owner_[faceToKeep[1]] = owner_[faceToThrow[1]];
2633 setFlip(faceToKeep[1]);
2634 }
2635 else
2636 {
2637 // Keep orientation intact, and update the neighbour
2638 neighbour_[faceToKeep[1]] = owner_[faceToThrow[1]];
2639 }
2640 }
2641 else
2642 {
2643 // Keep orientation intact, and update the owner
2644 owner_[faceToKeep[1]] = owner_[faceToThrow[1]];
2645 }
2646 }
2647 }
2649 // Remove orphaned faces
2650 if (owner_[faceToKeep[0]] == -1)
2651 {
2652 const face& keepFace = faces_[faceToKeep[0]];
2653 const labelList& rmFE = faceEdges_[faceToKeep[0]];
2655 labelList keepPoints(keepFace.size(), 0);
2657 forAll(rmFE, edgeI)
2658 {
2659 label eIndex = rmFE[edgeI];
2660 labelList& eFaces = edgeFaces_[eIndex];
2661 const edge& checkEdge = edges_[eIndex];
2663 if
2664 (
2665 (eFaces.size() == 1) &&
2666 (eFaces[0] == faceToKeep[0])
2667 )
2668 {
2669 // This edge has to be removed entirely.
2670 removeEdge(eIndex);
2672 // Update map
2673 map.removeEdge(eIndex);
2674 }
2675 else
2676 {
2677 // Mark points
2678 keepPoints[keepFace.which(checkEdge[0])] = 1;
2679 keepPoints[keepFace.which(checkEdge[1])] = 1;
2681 // Size-down edgeFaces
2682 meshOps::sizeDownList
2683 (
2684 faceToKeep[0],
2685 eFaces
2686 );
2687 }
2688 }
2690 // Check if processor patches are being
2691 // converted into physical ones
2692 label neiProc = -1;
2693 label kfPatch = whichPatch(faceToKeep[0]);
2694 label rfPatch = whichPatch(faceToThrow[0]);
2696 if
2697 (
2698 (getNeighbourProcessor(kfPatch) == -1) &&
2699 ((neiProc = getNeighbourProcessor(rfPatch)) > -1)
2700 )
2701 {
2702 if (isSubMesh_)
2703 {
2704 // If this is a subMesh, and faceToKeep is on
2705 // a physical boundary, make a 'special' entry
2706 // for coupled mapping purposes.
2707 map.addFace
2708 (
2709 faceToThrow[0],
2710 labelList(1, (-2 - kfPatch))
2711 );
2712 }
2713 else
2714 {
2715 // This is a wierd overhanging cell on the master
2716 // processor which is being removed entirely.
2717 // Since the processor patch face is being converted
2718 // to a physical one, the slave processor needs to
2719 // be notified of the change.
2720 forAll(procIndices_, pI)
2721 {
2722 if (procIndices_[pI] == neiProc)
2723 {
2724 // Find slave index from the face map
2725 const coupleMap& cMap = recvMeshes_[pI].map();
2727 label replaceFace =
2728 (
2729 cMap.findSlave
2730 (
2731 coupleMap::FACE,
2732 faceToThrow[0]
2733 )
2734 );
2736 if (replaceFace == -1)
2737 {
2738 // Something is wrong here.
2739 Pout<< " Could not find slave face for: " << nl
2740 << " Master face: " << faceToThrow[0]
2741 << " :: " << faces_[faceToThrow[0]] << nl
2742 << " on proc: " << procIndices_[pI] << nl
2743 << " on converting to patch: " << kfPatch
2744 << abort(FatalError);
2745 }
2747 cMap.pushOperation
2748 (
2749 replaceFace,
2750 coupleMap::CONVERT_PHYSICAL,
2751 kfPatch
2752 );
2754 break;
2755 }
2756 }
2757 }
2758 }
2760 // Remove unused points
2761 forAll(keepPoints, pointI)
2762 {
2763 if (keepPoints[pointI] == 0)
2764 {
2765 removePoint(keepFace[pointI]);
2766 map.removePoint(keepFace[pointI]);
2767 }
2768 }
2770 // Remove the face
2771 removeFace(faceToKeep[0]);
2773 // Update map
2774 map.removeFace(faceToKeep[0]);
2775 }
2776 else
2777 if
2778 (
2779 (neighbour_[faceToKeep[0]] == -1)
2780 && (whichPatch(faceToKeep[0]) < 0)
2781 )
2782 {
2783 // Obtain a copy before adding the new face,
2784 // since the reference might become invalid during list resizing.
2785 face newFace = faces_[faceToKeep[0]];
2786 label newOwn = owner_[faceToKeep[0]];
2787 labelList newFaceEdges = faceEdges_[faceToKeep[0]];
2789 // This face is being converted from interior to boundary. Remove
2790 // from the interior list and add as a boundary face to the end.
2791 // Edges don't have to change, since they're all on the boundary anyway.
2792 label newFaceIndex =
2793 (
2794 insertFace
2795 (
2796 whichPatch(faceToThrow[0]),
2797 newFace,
2798 newOwn,
2799 -1,
2800 newFaceEdges
2801 )
2802 );
2804 // Add an entry for mapping
2805 modifiedFaces.append(newFaceIndex);
2807 // Update map
2808 map.addFace(newFaceIndex, labelList(1, faceToThrow[0]));
2810 meshOps::replaceLabel
2811 (
2812 faceToKeep[0],
2813 newFaceIndex,
2814 cells_[newOwn]
2815 );
2817 // Correct edgeFaces with the new face label.
2818 forAll(newFaceEdges, edgeI)
2819 {
2820 meshOps::replaceLabel
2821 (
2822 faceToKeep[0],
2823 newFaceIndex,
2824 edgeFaces_[newFaceEdges[edgeI]]
2825 );
2826 }
2828 // Renumber the neighbour so that this face is removed correctly.
2829 neighbour_[faceToKeep[0]] = 0;
2830 removeFace(faceToKeep[0]);
2832 // Update map
2833 map.removeFace(faceToKeep[0]);
2834 }
2836 // Remove the unwanted faces in the cell(s) adjacent to this face,
2837 // and correct the cells that contain discarded faces
2838 const cell& cell_0 = cells_[c0];
2840 forAll(cell_0,faceI)
2841 {
2842 if (cell_0[faceI] != fIndex && cell_0[faceI] != faceToKeep[0])
2843 {
2844 removeFace(cell_0[faceI]);
2846 // Update map
2847 map.removeFace(cell_0[faceI]);
2848 }
2849 }
2851 if (cellCheck[0] != -1)
2852 {
2853 meshOps::replaceLabel
2854 (
2855 faceToThrow[0],
2856 faceToKeep[0],
2857 cells_[cellCheck[0]]
2858 );
2859 }
2861 // Remove the cell
2862 removeCell(c0);
2864 // Update map
2865 map.removeCell(c0);
2867 if (c1 == -1)
2868 {
2869 // Increment the surface-collapse counter
2870 status(SURFACE_COLLAPSES)++;
2871 }
2872 else
2873 {
2874 // Remove orphaned faces
2875 if (owner_[faceToKeep[1]] == -1)
2876 {
2877 const face& keepFace = faces_[faceToKeep[1]];
2878 const labelList& rmFE = faceEdges_[faceToKeep[1]];
2880 labelList keepPoints(keepFace.size(), 0);
2882 forAll(rmFE, edgeI)
2883 {
2884 label eIndex = rmFE[edgeI];
2885 labelList& eFaces = edgeFaces_[eIndex];
2886 const edge& checkEdge = edges_[eIndex];
2888 if
2889 (
2890 (eFaces.size() == 1) &&
2891 (eFaces[0] == faceToKeep[1])
2892 )
2893 {
2894 // This edge has to be removed entirely.
2895 removeEdge(eIndex);
2897 // Update map
2898 map.removeEdge(eIndex);
2899 }
2900 else
2901 {
2902 // Mark points
2903 keepPoints[keepFace.which(checkEdge[0])] = 1;
2904 keepPoints[keepFace.which(checkEdge[1])] = 1;
2906 // Size-down edgeFaces
2907 meshOps::sizeDownList
2908 (
2909 faceToKeep[1],
2910 eFaces
2911 );
2912 }
2913 }
2915 // Check if processor patches are being
2916 // converted into physical ones
2917 label neiProc = -1;
2918 label kfPatch = whichPatch(faceToKeep[1]);
2919 label rfPatch = whichPatch(faceToThrow[1]);
2921 if
2922 (
2923 (getNeighbourProcessor(kfPatch) == -1) &&
2924 ((neiProc = getNeighbourProcessor(rfPatch)) > -1)
2925 )
2926 {
2927 if (isSubMesh_)
2928 {
2929 // If this is a subMesh, and faceToKeep is on
2930 // a physical boundary, make a 'special' entry
2931 // for coupled mapping purposes.
2932 map.addFace
2933 (
2934 faceToThrow[1],
2935 labelList(1, (-2 - kfPatch))
2936 );
2937 }
2938 else
2939 {
2940 // This is a wierd overhanging cell on the master
2941 // processor which is being removed entirely.
2942 // Since the processor patch face is being converted
2943 // to a physical one, the slave processor needs to
2944 // be notified of the change.
2945 forAll(procIndices_, pI)
2946 {
2947 if (procIndices_[pI] == neiProc)
2948 {
2949 // Find slave index from the face map
2950 const coupleMap& cMap = recvMeshes_[pI].map();
2952 label replaceFace =
2953 (
2954 cMap.findSlave
2955 (
2956 coupleMap::FACE,
2957 faceToThrow[1]
2958 )
2959 );
2961 if (replaceFace == -1)
2962 {
2963 // Something is wrong here.
2964 Pout<< " Could not find slave face for: " << nl
2965 << " Master face: " << faceToThrow[1]
2966 << " :: " << faces_[faceToThrow[1]] << nl
2967 << " on proc: " << procIndices_[pI] << nl
2968 << " on converting to patch: " << kfPatch
2969 << abort(FatalError);
2970 }
2972 cMap.pushOperation
2973 (
2974 replaceFace,
2975 coupleMap::CONVERT_PHYSICAL,
2976 kfPatch
2977 );
2979 break;
2980 }
2981 }
2982 }
2983 }
2985 // Remove unused points
2986 forAll(keepPoints, pointI)
2987 {
2988 if (keepPoints[pointI] == 0)
2989 {
2990 removePoint(keepFace[pointI]);
2991 map.removePoint(keepFace[pointI]);
2992 }
2993 }
2995 // Remove the face
2996 removeFace(faceToKeep[1]);
2998 // Update map
2999 map.removeFace(faceToKeep[1]);
3000 }
3001 else
3002 if
3003 (
3004 (neighbour_[faceToKeep[1]] == -1)
3005 && (whichPatch(faceToKeep[1]) < 0)
3006 )
3007 {
3008 // Obtain a copy before adding the new face,
3009 // since the reference might become invalid during list resizing.
3010 face newFace = faces_[faceToKeep[1]];
3011 label newOwn = owner_[faceToKeep[1]];
3012 labelList newFaceEdges = faceEdges_[faceToKeep[1]];
3014 // This face is being converted from interior to boundary. Remove
3015 // from the interior list and add as a boundary face to the end.
3016 // Edges don't have to change, since they're on the boundary anyway.
3017 label newFaceIndex =
3018 (
3019 insertFace
3020 (
3021 whichPatch(faceToThrow[1]),
3022 newFace,
3023 newOwn,
3024 -1,
3025 newFaceEdges
3026 )
3027 );
3029 // Add an entry for mapping
3030 modifiedFaces.append(newFaceIndex);
3032 // Update map
3033 map.addFace(newFaceIndex, labelList(1, faceToThrow[1]));
3035 meshOps::replaceLabel
3036 (
3037 faceToKeep[1],
3038 newFaceIndex,
3039 cells_[newOwn]
3040 );
3042 // Correct edgeFaces with the new face label.
3043 forAll(newFaceEdges, edgeI)
3044 {
3045 meshOps::replaceLabel
3046 (
3047 faceToKeep[1],
3048 newFaceIndex,
3049 edgeFaces_[newFaceEdges[edgeI]]
3050 );
3051 }
3053 // Renumber the neighbour so that this face is removed correctly.
3054 neighbour_[faceToKeep[1]] = 0;
3055 removeFace(faceToKeep[1]);
3057 // Update map
3058 map.removeFace(faceToKeep[1]);
3059 }
3061 const cell& cell_1 = cells_[c1];
3063 forAll(cell_1, faceI)
3064 {
3065 if (cell_1[faceI] != fIndex && cell_1[faceI] != faceToKeep[1])
3066 {
3067 removeFace(cell_1[faceI]);
3069 // Update map
3070 map.removeFace(cell_1[faceI]);
3071 }
3072 }
3074 if (cellCheck[1] != -1)
3075 {
3076 meshOps::replaceLabel
3077 (
3078 faceToThrow[1],
3079 faceToKeep[1],
3080 cells_[cellCheck[1]]
3081 );
3082 }
3084 // Remove the cell
3085 removeCell(c1);
3087 // Update map
3088 map.removeCell(c1);
3089 }
3091 // Move old / new points
3092 oldPoints_[replacement[0]] = oldPoint[0];
3093 oldPoints_[replacement[1]] = oldPoint[1];
3095 points_[replacement[0]] = newPoint[0];
3096 points_[replacement[1]] = newPoint[1];
3098 // Finally remove the face
3099 removeFace(fIndex);
3101 // Update map
3102 map.removeFace(fIndex);
3104 // Write out VTK files after change
3105 if (debug > 3)
3106 {
3107 labelHashSet vtkCells;
3109 forAll(hullCells[0], cellI)
3110 {
3111 if (hullCells[0][cellI] == c0 || hullCells[0][cellI] == c1)
3112 {
3113 continue;
3114 }
3116 if (!vtkCells.found(hullCells[0][cellI]))
3117 {
3118 vtkCells.insert(hullCells[0][cellI]);
3119 }
3120 }
3122 forAll(hullCells[1], cellI)
3123 {
3124 if (hullCells[1][cellI] == c0 || hullCells[1][cellI] == c1)
3125 {
3126 continue;
3127 }
3129 if (!vtkCells.found(hullCells[1][cellI]))
3130 {
3131 vtkCells.insert(hullCells[1][cellI]);
3132 }
3133 }
3135 writeVTK
3136 (
3137 Foam::name(fIndex)
3138 + "_Collapse_1",
3139 vtkCells.toc(),
3140 3, false, true
3141 );
3142 }
3144 // Fill-in candidate mapping information
3145 labelList mC(2, -1);
3146 mC[0] = c0, mC[1] = c1;
3148 // Now that all old / new cells possess correct connectivity,
3149 // compute mapping information.
3150 forAll(hullCells, indexI)
3151 {
3152 forAll(hullCells[indexI], cellI)
3153 {
3154 label mcIndex = hullCells[indexI][cellI];
3156 // Skip collapsed cells
3157 if (mcIndex == c0 || mcIndex == c1)
3158 {
3159 continue;
3160 }
3162 // Set the mapping for this cell
3163 setCellMapping(mcIndex, mC);
3164 }
3165 }
3167 // Set face mapping information for modified faces
3168 forAll(modifiedFaces, faceI)
3169 {
3170 const label mfIndex = modifiedFaces[faceI];
3172 // Exclude deleted faces
3173 if (faces_[mfIndex].empty())
3174 {
3175 continue;
3176 }
3178 // Decide between default / weighted mapping
3179 // based on boundary information
3180 label fPatch = whichPatch(mfIndex);
3182 if (fPatch == -1)
3183 {
3184 setFaceMapping(mfIndex);
3185 }
3186 else
3187 {
3188 // Fill-in candidate mapping information
3189 labelList faceCandidates;
3191 const labelList& fEdges = faceEdges_[mfIndex];
3193 forAll(fEdges, edgeI)
3194 {
3195 if (whichEdgePatch(fEdges[edgeI]) == fPatch)
3196 {
3197 // Loop through associated edgeFaces
3198 const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
3200 forAll(eFaces, faceI)
3201 {
3202 if
3203 (
3204 (eFaces[faceI] != mfIndex) &&
3205 (whichPatch(eFaces[faceI]) == fPatch)
3206 )
3207 {
3208 faceCandidates.setSize
3209 (
3210 faceCandidates.size() + 1,
3211 eFaces[faceI]
3212 );
3213 }
3214 }
3215 }
3216 }
3218 if (faceCandidates.empty())
3219 {
3220 // Add the face itself
3221 faceCandidates.setSize(1, mfIndex);
3222 }
3224 // Set the mapping for this face
3225 setFaceMapping(mfIndex, faceCandidates);
3226 }
3227 }
3229 // If modification is coupled, update mapping info.
3230 if (coupledModification_)
3231 {
3232 // Check if the collapse points are present
3233 // on a processor not involved in the current
3234 // operation, and update if necessary
3235 if (procCouple && !localCouple)
3236 {
3237 forAll(procIndices_, pI)
3238 {
3239 bool involved = false;
3241 forAll(slaveMaps, slaveI)
3242 {
3243 // Alias for convenience...
3244 const changeMap& slaveMap = slaveMaps[slaveI];
3246 if (slaveMap.patchIndex() == pI && slaveMap.index() >= 0)
3247 {
3248 // Involved in this operation. Break out.
3249 involved = true;
3250 break;
3251 }
3252 }
3254 if (involved)
3255 {
3256 continue;
3257 }
3259 // Check coupleMaps for point coupling
3260 const label pointEnum = coupleMap::POINT;
3262 const coupledMesh& recvMesh = recvMeshes_[pI];
3263 const coupleMap& cMap = recvMesh.map();
3265 // Obtain non-const references
3266 Map<label>& pointMap = cMap.entityMap(pointEnum);
3267 Map<label>& rPointMap = cMap.reverseEntityMap(pointEnum);
3269 label sI0 = -1, sI1 = -1;
3271 if (collapsingSlave)
3272 {
3273 if ((sI0 = cMap.findMaster(pointEnum, original[0])) > -1)
3274 {
3275 if (rPointMap.found(replacement[0]))
3276 {
3277 rPointMap[replacement[0]] = sI0;
3278 }
3279 else
3280 {
3281 rPointMap.insert(replacement[0], sI0);
3282 }
3284 pointMap[sI0] = replacement[0];
3285 }
3287 if ((sI1 = cMap.findMaster(pointEnum, original[1])) > -1)
3288 {
3289 if (rPointMap.found(replacement[1]))
3290 {
3291 rPointMap[replacement[1]] = sI1;
3292 }
3293 else
3294 {
3295 rPointMap.insert(replacement[1], sI1);
3296 }
3298 pointMap[sI1] = replacement[1];
3299 }
3300 }
3301 else
3302 {
3303 if ((sI0 = cMap.findSlave(pointEnum, original[0])) > -1)
3304 {
3305 if (pointMap.found(replacement[0]))
3306 {
3307 pointMap[replacement[0]] = sI0;
3308 }
3309 else
3310 {
3311 pointMap.insert(replacement[0], sI0);
3312 }
3314 rPointMap[sI0] = replacement[0];
3315 }
3317 if ((sI1 = cMap.findSlave(pointEnum, original[1])) > -1)
3318 {
3319 if (pointMap.found(replacement[1]))
3320 {
3321 pointMap[replacement[1]] = sI1;
3322 }
3323 else
3324 {
3325 pointMap.insert(replacement[1], sI1);
3326 }
3328 rPointMap[sI1] = replacement[1];
3329 }
3330 }
3332 if (sI0 > -1 && sI1 > -1 && debug > 2)
3333 {
3334 Pout<< " Found " << original[0] << " and " << original[1]
3335 << " on proc: " << procIndices_[pI]
3336 << endl;
3337 }
3338 }
3339 }
3341 forAll(slaveMaps, slaveI)
3342 {
3343 // Alias for convenience...
3344 const changeMap& slaveMap = slaveMaps[slaveI];
3346 // Skip updates for edge-based coupling
3347 if (slaveMap.index() < 0)
3348 {
3349 continue;
3350 }
3352 label pI = slaveMap.patchIndex();
3354 // Fetch the appropriate coupleMap
3355 const coupleMap* cMapPtr = NULL;
3357 if (localCouple && !procCouple)
3358 {
3359 cMapPtr = &(patchCoupling_[pI].map());
3360 }
3361 else
3362 if (procCouple && !localCouple)
3363 {
3364 cMapPtr = &(recvMeshes_[pI].map());
3365 }
3367 // Configure the slave replacement points.
3368 // - collapseQuadFace stores this as 'addedPoints'
3369 label scP0 = slaveMap.removedPointList()[0];
3370 label scP1 = slaveMap.removedPointList()[1];
3372 label srP0 = slaveMap.addedPointList()[0].index();
3373 label srP1 = slaveMap.addedPointList()[1].index();
3375 // Alias for convenience
3376 const coupleMap& cMap = *cMapPtr;
3378 // Remove the point entries.
3379 const label pointEnum = coupleMap::POINT;
3381 // Obtain references
3382 Map<label>& pointMap = cMap.entityMap(pointEnum);
3383 Map<label>& rPointMap = cMap.reverseEntityMap(pointEnum);
3385 if (collapsingSlave)
3386 {
3387 if (rPointMap[replacement[0]] == scP0)
3388 {
3389 pointMap[srP0] = replacement[0];
3390 rPointMap[replacement[0]] = srP0;
3391 }
3392 else
3393 if (rPointMap[replacement[0]] == scP1)
3394 {
3395 pointMap[srP1] = replacement[0];
3396 rPointMap[replacement[0]] = srP1;
3397 }
3399 if (rPointMap[original[0]] == scP0)
3400 {
3401 pointMap.erase(scP0);
3402 }
3403 else
3404 if (rPointMap[original[0]] == scP1)
3405 {
3406 pointMap.erase(scP1);
3407 }
3409 rPointMap.erase(original[0]);
3410 }
3411 else
3412 {
3413 if (pointMap[replacement[0]] == scP0)
3414 {
3415 rPointMap[srP0] = replacement[0];
3416 pointMap[replacement[0]] = srP0;
3417 }
3418 else
3419 if (pointMap[replacement[0]] == scP1)
3420 {
3421 rPointMap[srP1] = replacement[0];
3422 pointMap[replacement[0]] = srP1;
3423 }
3425 if (pointMap[original[0]] == scP0)
3426 {
3427 rPointMap.erase(scP0);
3428 }
3429 else
3430 if (pointMap[original[0]] == scP1)
3431 {
3432 rPointMap.erase(scP1);
3433 }
3435 pointMap.erase(original[0]);
3436 }
3438 if (collapsingSlave)
3439 {
3440 if (rPointMap[replacement[1]] == scP0)
3441 {
3442 pointMap[srP0] = replacement[1];
3443 rPointMap[replacement[1]] = srP0;
3444 }
3445 else
3446 if (rPointMap[replacement[1]] == scP1)
3447 {
3448 pointMap[srP1] = replacement[1];
3449 rPointMap[replacement[1]] = srP1;
3450 }
3452 if (rPointMap[original[1]] == scP0)
3453 {
3454 pointMap.erase(scP0);
3455 }
3456 else
3457 if (rPointMap[original[1]] == scP1)
3458 {
3459 pointMap.erase(scP1);
3460 }
3462 rPointMap.erase(original[1]);
3463 }
3464 else
3465 {
3466 if (pointMap[replacement[1]] == scP0)
3467 {
3468 rPointMap[srP0] = replacement[1];
3469 pointMap[replacement[1]] = srP0;
3470 }
3471 else
3472 if (pointMap[replacement[1]] == scP1)
3473 {
3474 rPointMap[srP1] = replacement[1];
3475 pointMap[replacement[1]] = srP1;
3476 }
3478 if (pointMap[original[1]] == scP0)
3479 {
3480 rPointMap.erase(scP0);
3481 }
3482 else
3483 if (pointMap[original[1]] == scP1)
3484 {
3485 rPointMap.erase(scP1);
3486 }
3488 pointMap.erase(original[1]);
3489 }
3491 // Remove the face entries
3492 const label faceEnum = coupleMap::FACE;
3494 // Obtain references
3495 Map<label>& faceMap = cMap.entityMap(faceEnum);
3496 Map<label>& rFaceMap = cMap.reverseEntityMap(faceEnum);
3498 if (collapsingSlave)
3499 {
3500 faceMap.erase(faceMap[fIndex]);
3501 rFaceMap.erase(fIndex);
3502 }
3503 else
3504 {
3505 rFaceMap.erase(faceMap[fIndex]);
3506 faceMap.erase(fIndex);
3507 }
3509 // Configure a comparison face
3510 face cFace(4);
3512 // If any interior faces in the master map were
3513 // converted to boundaries, account for it
3514 const List<objectMap>& madF = map.addedFaceList();
3516 forAll(madF, faceI)
3517 {
3518 label mfIndex = madF[faceI].index();
3519 const face& mFace = faces_[mfIndex];
3521 // Select appropriate mesh
3522 const dynamicTopoFvMesh* meshPtr = NULL;
3524 // Fetch the appropriate coupleMap
3525 const coupleMap* crMapPtr = NULL;
3527 // Fetch patch info
3528 label ofPatch = whichPatch(fIndex);
3529 label mfPatch = whichPatch(mfIndex);
3531 if (localCouple && !procCouple)
3532 {
3533 // Local coupling. Use this mesh itself
3534 meshPtr = this;
3535 crMapPtr = &(patchCoupling_[pI].map());
3536 }
3537 else
3538 if (procCouple && !localCouple)
3539 {
3540 // Occasionally, this face might talk to
3541 // a processor other than the slave
3542 if (ofPatch == mfPatch)
3543 {
3544 meshPtr = &(recvMeshes_[pI].subMesh());
3545 crMapPtr = &(recvMeshes_[pI].map());
3546 }
3547 else
3548 {
3549 label neiProcNo = getNeighbourProcessor(mfPatch);
3551 if (neiProcNo == -1)
3552 {
3553 // Not a processor patch. No mapping required.
3554 continue;
3555 }
3556 else
3557 {
3558 // Find an appropriate subMesh
3559 label prI = findIndex(procIndices_, neiProcNo);
3561 meshPtr = &(recvMeshes_[prI].subMesh());
3562 crMapPtr = &(recvMeshes_[prI].map());
3563 }
3564 }
3565 }
3567 bool matchedFace = false;
3569 // Alias for convenience
3570 const coupleMap& crMap = *crMapPtr;
3571 const dynamicTopoFvMesh& sMesh = *meshPtr;
3573 // Obtain references
3574 Map<label>& pMap = crMap.entityMap(pointEnum);
3575 Map<label>& fMap = crMap.entityMap(faceEnum);
3576 Map<label>& rfMap = crMap.reverseEntityMap(faceEnum);
3578 // Configure the face
3579 forAll(cFace, pointI)
3580 {
3581 cFace[pointI] =
3582 (
3583 pMap.found(mFace[pointI]) ?
3584 pMap[mFace[pointI]] : -1
3585 );
3586 }
3588 // Loop through all boundary faces on the subMesh
3589 for
3590 (
3591 label faceJ = sMesh.nOldInternalFaces_;
3592 faceJ < sMesh.faces_.size();
3593 faceJ++
3594 )
3595 {
3596 const face& sFace = sMesh.faces_[faceJ];
3598 if (face::compare(sFace, cFace))
3599 {
3600 if (debug > 2)
3601 {
3602 Pout<< " Found face: " << faceJ
3603 << "::" << sFace
3604 << " with mfIndex: " << mfIndex
3605 << "::" << mFace
3606 << endl;
3607 }
3609 // Update faceMaps
3610 if (collapsingSlave)
3611 {
3612 if (fMap.found(faceJ))
3613 {
3614 fMap[faceJ] = mfIndex;
3615 }
3616 else
3617 {
3618 fMap.insert(faceJ, mfIndex);
3619 }
3621 if (rfMap.found(mfIndex))
3622 {
3623 rfMap[mfIndex] = faceJ;
3624 }
3625 else
3626 {
3627 rfMap.insert(mfIndex, faceJ);
3628 }
3629 }
3630 else
3631 {
3632 if (rfMap.found(faceJ))
3633 {
3634 rfMap[faceJ] = mfIndex;
3635 }
3636 else
3637 {
3638 rfMap.insert(faceJ, mfIndex);
3639 }
3641 if (fMap.found(mfIndex))
3642 {
3643 fMap[mfIndex] = faceJ;
3644 }
3645 else
3646 {
3647 fMap.insert(mfIndex, faceJ);
3648 }
3649 }
3651 matchedFace = true;
3653 break;
3654 }
3655 }
3657 if (!matchedFace)
3658 {
3659 // Write out for post-processing
3660 writeVTK("masterFace_" + Foam::name(mfIndex), mfIndex, 2);
3662 FatalErrorIn
3663 (
3664 "\n"
3665 "const changeMap "
3666 "dynamicTopoFvMesh::collapseQuadFace\n"
3667 "(\n"
3668 " const label fIndex,\n"
3669 " label overRideCase,\n"
3670 " bool checkOnly,\n"
3671 " bool forceOp\n"
3672 ")\n"
3673 )
3674 << " Master face: " << mfIndex
3675 << ": " << mFace << " could not be matched." << nl
3676 << " cFace: " << cFace << nl
3677 << abort(FatalError);
3678 }
3679 }
3681 // If any interior faces in the slave map were
3682 // converted to boundaries, account for it
3683 const List<objectMap>& sadF = slaveMap.addedFaceList();
3685 forAll(sadF, faceI)
3686 {
3687 label sIndex = slaveMap.index();
3688 label sfIndex = sadF[faceI].index();
3690 // Select appropriate mesh
3691 const dynamicTopoFvMesh* meshPtr = NULL;
3693 if (localCouple && !procCouple)
3694 {
3695 // Local coupling. Use this mesh itself
3696 meshPtr = this;
3697 }
3698 else
3699 if (procCouple && !localCouple)
3700 {
3701 meshPtr = &(recvMeshes_[pI].subMesh());
3702 }
3704 // Alias for convenience
3705 const dynamicTopoFvMesh& sMesh = *meshPtr;
3707 label ofPatch = sMesh.whichPatch(sIndex);
3708 label sfPatch = sMesh.whichPatch(sfIndex);
3710 // Skip dissimilar patches
3711 if (ofPatch != sfPatch && localCouple)
3712 {
3713 continue;
3714 }
3716 const face& sFace = sMesh.faces_[sfIndex];
3718 if (sFace.empty())
3719 {
3720 // Slave face was removed. Update map.
3721 Map<label>::iterator sit = rFaceMap.find(sfIndex);
3723 label mfIndex = -1;
3725 if (sit != rFaceMap.end())
3726 {
3727 mfIndex = sit();
3728 faceMap.erase(sit());
3729 rFaceMap.erase(sit);
3730 }
3732 // Check if this is a special entry
3733 label mo =
3734 (
3735 sadF[faceI].masterObjects().size() ?
3736 sadF[faceI].masterObjects()[0] : 0
3737 );
3739 // Check if a patch conversion is necessary
3740 label newPatch = -1;
3742 // Back out the physical patch ID
3743 if (mo < 0 && mfIndex > -1)
3744 {
3745 newPatch = Foam::mag(mo + 2);
3746 }
3747 else
3748 {
3749 continue;
3750 }
3752 // Obtain a copy before adding the new face,
3753 // since the reference might become
3754 // invalid during list resizing.
3755 face newFace = faces_[mfIndex];
3756 label newOwn = owner_[mfIndex];
3757 labelList newFaceEdges = faceEdges_[mfIndex];
3759 label newFaceIndex =
3760 (
3761 insertFace
3762 (
3763 newPatch,
3764 newFace,
3765 newOwn,
3766 -1,
3767 newFaceEdges
3768 )
3769 );
3771 meshOps::replaceLabel
3772 (
3773 mfIndex,
3774 newFaceIndex,
3775 cells_[newOwn]
3776 );
3778 // Correct edgeFaces with the new face label.
3779 forAll(newFaceEdges, edgeI)
3780 {
3781 meshOps::replaceLabel
3782 (
3783 mfIndex,
3784 newFaceIndex,
3785 edgeFaces_[newFaceEdges[edgeI]]
3786 );
3787 }
3789 // Finally remove the face
3790 removeFace(mfIndex);
3792 // Update map
3793 map.removeFace(mfIndex);
3794 map.addFace(newFaceIndex, labelList(1, mfIndex));
3796 continue;
3797 }
3798 else
3799 {
3800 forAll(cFace, pointI)
3801 {
3802 cFace[pointI] =
3803 (
3804 rPointMap.found(sFace[pointI]) ?
3805 rPointMap[sFace[pointI]] : -1
3806 );
3807 }
3808 }
3810 bool matchedFace = false;
3812 // Loop through all boundary faces on this mesh
3813 for
3814 (
3815 label faceJ = nOldInternalFaces_;
3816 faceJ < faces_.size();
3817 faceJ++
3818 )
3819 {
3820 const face& mFace = faces_[faceJ];
3822 if (face::compare(mFace, cFace))
3823 {
3824 if (debug > 2)
3825 {
3826 Pout<< " Found face: " << faceJ
3827 << "::" << mFace
3828 << " with sfIndex: " << sfIndex
3829 << "::" << sFace
3830 << endl;
3831 }
3833 // Update faceMaps
3834 if (collapsingSlave)
3835 {
3836 if (rFaceMap.found(faceJ))
3837 {
3838 rFaceMap[faceJ] = sfIndex;
3839 }
3840 else
3841 {
3842 rFaceMap.insert(faceJ, sfIndex);
3843 }
3845 if (faceMap.found(sfIndex))
3846 {
3847 faceMap[sfIndex] = faceJ;
3848 }
3849 else
3850 {
3851 faceMap.insert(sfIndex, faceJ);
3852 }
3853 }
3854 else
3855 {
3856 if (faceMap.found(faceJ))
3857 {
3858 faceMap[faceJ] = sfIndex;
3859 }
3860 else
3861 {
3862 faceMap.insert(faceJ, sfIndex);
3863 }
3865 if (rFaceMap.found(sfIndex))
3866 {
3867 rFaceMap[sfIndex] = faceJ;
3868 }
3869 else
3870 {
3871 rFaceMap.insert(sfIndex, faceJ);
3872 }
3873 }
3875 matchedFace = true;
3877 break;
3878 }
3879 }
3881 if (!matchedFace)
3882 {
3883 FatalErrorIn
3884 (
3885 "\n"
3886 "const changeMap "
3887 "dynamicTopoFvMesh::collapseQuadFace\n"
3888 "(\n"
3889 " const label fIndex,\n"
3890 " label overRideCase,\n"
3891 " bool checkOnly,\n"
3892 " bool forceOp\n"
3893 ")\n"
3894 )
3895 << " Slave face: " << sfIndex
3896 << ": " << sFace << " could not be matched." << nl
3897 << " cFace: " << cFace << nl
3898 << abort(FatalError);
3899 }
3900 }
3901 }
3902 }
3904 // Set the flag
3905 topoChangeFlag_ = true;
3907 // Increment the counter
3908 status(TOTAL_COLLAPSES)++;
3910 // Increment the number of modifications
3911 status(TOTAL_MODIFICATIONS)++;
3913 // Return a succesful collapse
3914 map.type() = collapseCase;
3916 return map;
3917 }
3920 // Method for the collapse of an edge in 3D
3921 // - Returns a changeMap with a type specifying:
3922 // -3: Collapse failed since edge was on a noRefinement patch.
3923 // -1: Collapse failed since max number of topo-changes was reached.
3924 // 0: Collapse could not be performed.
3925 // 1: Collapsed to first node.
3926 // 2: Collapsed to second node.
3927 // 3: Collapsed to mid-point (default).
3928 // - overRideCase is used to force a certain collapse configuration.
3929 // -1: Use this value to let collapseEdge decide a case.
3930 // 1: Force collapse to first node.
3931 // 2: Force collapse to second node.
3932 // 3: Force collapse to mid-point.
3933 // - checkOnly performs a feasibility check and returns without modifications.
3934 // - forceOp to force the collapse.
3935 const changeMap dynamicTopoFvMesh::collapseEdge
3936 (
3937 const label eIndex,
3938 label overRideCase,
3939 bool checkOnly,
3940 bool forceOp
3941 )
3942 {
3943 // Edge collapse performs the following operations:
3944 // [1] Checks if either vertex of the edge is on a boundary
3945 // [2] Checks whether cells attached to deleted vertices will be valid
3946 // after the edge-collapse operation
3947 // [3] Deletes all cells surrounding this edge
3948 // [4] Deletes all faces surrounding this edge
3949 // [5] Deletes all faces surrounding the deleted vertex attached
3950 // to the cells in [3]
3951 // [6] Checks the orientation of faces connected to the retained
3952 // vertices
3953 // [7] Remove one of the vertices of the edge
3954 // Update faceEdges and edgeFaces information
3956 // For 2D meshes, perform face-collapse
3957 if (is2D())
3958 {
3959 return collapseQuadFace(eIndex, overRideCase, checkOnly);
3960 }
3962 // Figure out which thread this is...
3963 label tIndex = self();
3965 // Prepare the changeMaps
3966 changeMap map;
3967 List<changeMap> slaveMaps;
3968 bool collapsingSlave = false;
3970 if
3971 (
3972 (status(TOTAL_MODIFICATIONS) > maxModifications_)
3973 && (maxModifications_ > -1)
3974 )
3975 {
3976 // Reached the max allowable topo-changes.
3977 stack(tIndex).clear();
3979 return map;
3980 }
3982 // Check if edgeRefinements are to be avoided on patch.
3983 const labelList& eF = edgeFaces_[eIndex];
3985 forAll(eF, fI)
3986 {
3987 label fPatch = whichPatch(eF[fI]);
3989 if (baseMesh_.lengthEstimator().checkRefinementPatch(fPatch))
3990 {
3991 map.type() = -3;
3993 return map;
3994 }
3995 }
3997 // Sanity check: Is the index legitimate?
3998 if (eIndex < 0)
3999 {
4000 FatalErrorIn
4001 (
4002 "\n"
4003 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
4004 "(\n"
4005 " const label eIndex,\n"
4006 " label overRideCase,\n"
4007 " bool checkOnly,\n"
4008 " bool forceOp\n"
4009 ")\n"
4010 )
4011 << " Invalid index: " << eIndex
4012 << abort(FatalError);
4013 }
4015 // If coupled modification is set, and this is a
4016 // master edge, collapse its slaves as well.
4017 bool localCouple = false, procCouple = false;
4019 if (coupledModification_)
4020 {
4021 const edge& eCheck = edges_[eIndex];
4023 const label edgeEnum = coupleMap::EDGE;
4024 const label pointEnum = coupleMap::POINT;
4026 // Is this a locally coupled edge (either master or slave)?
4027 if (locallyCoupledEntity(eIndex, true))
4028 {
4029 localCouple = true;
4030 procCouple = false;
4031 }
4032 else
4033 if (processorCoupledEntity(eIndex))
4034 {
4035 procCouple = true;
4036 localCouple = false;
4037 }
4039 if (localCouple && !procCouple)
4040 {
4041 // Determine the slave index.
4042 label sIndex = -1, pIndex = -1;
4044 forAll(patchCoupling_, patchI)
4045 {
4046 if (patchCoupling_(patchI))
4047 {
4048 const coupleMap& cMap = patchCoupling_[patchI].map();
4050 if ((sIndex = cMap.findSlave(edgeEnum, eIndex)) > -1)
4051 {
4052 pIndex = patchI;
4054 break;
4055 }
4057 // The following bit happens only during the sliver
4058 // exudation process, since slave edges are
4059 // usually not added to the coupled edge-stack.
4060 if ((sIndex = cMap.findMaster(edgeEnum, eIndex)) > -1)
4061 {
4062 pIndex = patchI;
4064 // Notice that we are collapsing a slave edge first.
4065 collapsingSlave = true;
4067 break;
4068 }
4069 }
4070 }
4072 if (sIndex == -1)
4073 {
4074 FatalErrorIn
4075 (
4076 "\n"
4077 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
4078 "(\n"
4079 " const label eIndex,\n"
4080 " label overRideCase,\n"
4081 " bool checkOnly,\n"
4082 " bool forceOp\n"
4083 ")\n"
4084 )
4085 << "Coupled maps were improperly specified." << nl
4086 << " Slave index not found for: " << nl
4087 << " Edge: " << eIndex << ": " << eCheck << nl
4088 << abort(FatalError);
4089 }
4090 else
4091 {
4092 // If we've found the slave, size up the list
4093 meshOps::sizeUpList
4094 (
4095 changeMap(),
4096 slaveMaps
4097 );
4099 // Save index and patch for posterity
4100 slaveMaps[0].index() = sIndex;
4101 slaveMaps[0].patchIndex() = pIndex;
4102 }
4104 if (debug > 1)
4105 {
4106 Pout<< nl << " >> Collapsing slave edge: " << sIndex
4107 << " for master edge: " << eIndex << endl;
4108 }
4109 }
4110 else
4111 if (procCouple && !localCouple)
4112 {
4113 // If this is a new entity, bail out for now.
4114 // This will be handled at the next time-step.
4115 if (eIndex >= nOldEdges_)
4116 {
4117 return map;
4118 }
4120 // Check slaves
4121 forAll(procIndices_, pI)
4122 {
4123 // Fetch reference to subMeshes
4124 const coupledMesh& sendMesh = sendMeshes_[pI];
4125 const coupledMesh& recvMesh = recvMeshes_[pI];
4127 const coupleMap& scMap = sendMesh.map();
4128 const coupleMap& rcMap = recvMesh.map();
4130 // If this edge was sent to a lower-ranked
4131 // processor, skip it.
4132 if
4133 (
4134 priority
4135 (
4136 procIndices_[pI],
4137 lessOp<label>(),
4138 Pstream::myProcNo()
4139 )
4140 )
4141 {
4142 if (scMap.reverseEntityMap(edgeEnum).found(eIndex))
4143 {
4144 if (debug > 3)
4145 {
4146 Pout<< "Edge: " << eIndex
4147 << "::" << eCheck
4148 << " was sent to proc: "
4149 << procIndices_[pI]
4150 << ", so bailing out."
4151 << endl;
4152 }
4154 return map;
4155 }
4156 }
4158 label sIndex = -1;
4160 if ((sIndex = rcMap.findSlave(edgeEnum, eIndex)) > -1)
4161 {
4162 // Check if a lower-ranked processor is
4163 // handling this edge
4164 if
4165 (
4166 priority
4167 (
4168 procIndices_[pI],
4169 lessOp<label>(),
4170 Pstream::myProcNo()
4171 )
4172 )
4173 {
4174 if (debug > 3)
4175 {
4176 Pout<< "Edge: " << eIndex
4177 << "::" << eCheck
4178 << " is handled by proc: "
4179 << procIndices_[pI]
4180 << ", so bailing out."
4181 << endl;
4182 }
4184 return map;
4185 }
4187 label curIndex = slaveMaps.size();
4189 // Size up the list
4190 meshOps::sizeUpList
4191 (
4192 changeMap(),
4193 slaveMaps
4194 );
4196 // Save index and patch for posterity
4197 slaveMaps[curIndex].index() = sIndex;
4198 slaveMaps[curIndex].patchIndex() = pI;
4199 }
4200 else
4201 if
4202 (
4203 (rcMap.findSlave(pointEnum, eCheck[0]) > -1) ||
4204 (rcMap.findSlave(pointEnum, eCheck[1]) > -1)
4205 )
4206 {
4207 // A point-only coupling exists.
4209 // Check if a lower-ranked processor is
4210 // handling this edge
4211 if
4212 (
4213 priority
4214 (
4215 procIndices_[pI],
4216 lessOp<label>(),
4217 Pstream::myProcNo()
4218 )
4219 )
4220 {
4221 if (debug > 3)
4222 {
4223 Pout<< "Edge point on: " << eIndex
4224 << "::" << eCheck
4225 << " is handled by proc: "
4226 << procIndices_[pI]
4227 << ", so bailing out."
4228 << endl;
4229 }
4231 return map;
4232 }
4234 label p0Index = rcMap.findSlave(pointEnum, eCheck[0]);
4235 label p1Index = rcMap.findSlave(pointEnum, eCheck[1]);
4237 if (p0Index > -1 && p1Index == -1)
4238 {
4239 sIndex = p0Index;
4240 }
4241 else
4242 if (p0Index == -1 && p1Index > -1)
4243 {
4244 sIndex = p1Index;
4245 }
4247 label curIndex = slaveMaps.size();
4249 // Size up the list
4250 meshOps::sizeUpList
4251 (
4252 changeMap(),
4253 slaveMaps
4254 );
4256 // Save index and patch for posterity
4257 // - Negate the index to signify point coupling
4258 slaveMaps[curIndex].index() = -sIndex;
4259 slaveMaps[curIndex].patchIndex() = pI;
4260 }
4261 }
4262 }
4263 else
4264 {
4265 // Something's wrong with coupling maps
4266 FatalErrorIn
4267 (
4268 "\n"
4269 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
4270 "(\n"
4271 " const label eIndex,\n"
4272 " label overRideCase,\n"
4273 " bool checkOnly,\n"
4274 " bool forceOp\n"
4275 ")\n"
4276 )
4277 << "Coupled maps were improperly specified." << nl
4278 << " localCouple: " << localCouple << nl
4279 << " procCouple: " << procCouple << nl
4280 << " Edge: " << eIndex << ": " << eCheck << nl
4281 << abort(FatalError);
4282 }
4284 // Temporarily turn off coupledModification
4285 unsetCoupledModification();
4287 // Test the master edge for collapse, and decide on a case
4288 changeMap masterMap = collapseEdge(eIndex, -1, true, forceOp);
4290 // Turn it back on.
4291 setCoupledModification();
4293 // Master couldn't perform collapse
4294 if (masterMap.type() <= 0)
4295 {
4296 return masterMap;
4297 }
4299 // For point-only coupling, define the points for checking
4300 pointField slaveMoveNewPoint(slaveMaps.size(), vector::zero);
4301 pointField slaveMoveOldPoint(slaveMaps.size(), vector::zero);
4303 // Now check each of the slaves for collapse feasibility
4304 forAll(slaveMaps, slaveI)
4305 {
4306 // Alias for convenience...
4307 changeMap& slaveMap = slaveMaps[slaveI];
4309 label slaveOverRide = -1;
4310 label sIndex = slaveMap.index();
4311 label pI = slaveMap.patchIndex();
4313 // If the edge is mapped onto itself, skip check
4314 // (occurs for cyclic edges)
4315 if ((sIndex == eIndex) && localCouple)
4316 {
4317 continue;
4318 }
4320 const coupleMap* cMapPtr = NULL;
4322 edge mEdge(eCheck), sEdge(-1, -1);
4324 if (localCouple)
4325 {
4326 sEdge = edges_[sIndex];
4328 cMapPtr = &(patchCoupling_[pI].map());
4329 }
4330 else
4331 if (procCouple)
4332 {
4333 const dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
4335 cMapPtr = &(recvMeshes_[pI].map());
4337 if (sIndex < 0)
4338 {
4339 if (debug > 3)
4340 {
4341 Pout<< "Checking slave point: " << mag(sIndex)
4342 << "::" << sMesh.points_[mag(sIndex)]
4343 << " on proc: " << procIndices_[pI]
4344 << " for master edge: " << eIndex
4345 << " using collapseCase: " << masterMap.type()
4346 << endl;
4347 }
4348 }
4349 else
4350 {
4351 sEdge = sMesh.edges_[sIndex];
4353 if (debug > 3)
4354 {
4355 Pout<< "Checking slave edge: " << sIndex
4356 << "::" << sMesh.edges_[sIndex]
4357 << " on proc: " << procIndices_[pI]
4358 << " for master edge: " << eIndex
4359 << " using collapseCase: " << masterMap.type()
4360 << endl;
4361 }
4362 }
4363 }
4365 const Map<label>& pointMap = cMapPtr->entityMap(pointEnum);
4367 // Perform a topological comparison.
4368 switch (masterMap.type())
4369 {
4370 case 1:
4371 {
4372 if (sIndex < 0)
4373 {
4374 slaveMoveNewPoint[slaveI] = points_[mEdge[0]];
4375 slaveMoveOldPoint[slaveI] = oldPoints_[mEdge[0]];
4376 }
4377 else
4378 if (pointMap[mEdge[0]] == sEdge[0])
4379 {
4380 slaveOverRide = 1;
4381 }
4382 else
4383 if (pointMap[mEdge[1]] == sEdge[0])
4384 {
4385 slaveOverRide = 2;
4386 }
4387 else
4388 {
4389 // Write out for for post-processing
4390 writeVTK("mEdge_" + Foam::name(eIndex), eIndex, 1);
4392 FatalErrorIn
4393 (
4394 "\n"
4395 "const changeMap dynamicTopoFvMesh"
4396 "::collapseEdge\n"
4397 "(\n"
4398 " const label eIndex,\n"
4399 " label overRideCase,\n"
4400 " bool checkOnly,\n"
4401 " bool forceOp\n"
4402 ")\n"
4403 )
4404 << "Coupled collapse failed." << nl
4405 << "Master: " << eIndex << " : " << mEdge << nl
4406 << "Slave: " << sIndex << " : " << sEdge << nl
4407 << abort(FatalError);
4408 }
4410 break;
4411 }
4413 case 2:
4414 {
4415 if (sIndex < 0)
4416 {
4417 slaveMoveNewPoint[slaveI] = points_[mEdge[1]];
4418 slaveMoveOldPoint[slaveI] = oldPoints_[mEdge[1]];
4419 }
4420 else
4421 if (pointMap[mEdge[1]] == sEdge[1])
4422 {
4423 slaveOverRide = 2;
4424 }
4425 else
4426 if (pointMap[mEdge[0]] == sEdge[1])
4427 {
4428 slaveOverRide = 1;
4429 }
4430 else
4431 {
4432 // Write out for for post-processing
4433 writeVTK("mEdge_" + Foam::name(eIndex), eIndex, 1);
4435 FatalErrorIn
4436 (
4437 "\n"
4438 "const changeMap dynamicTopoFvMesh"
4439 "::collapseEdge\n"
4440 "(\n"
4441 " const label eIndex,\n"
4442 " label overRideCase,\n"
4443 " bool checkOnly,\n"
4444 " bool forceOp\n"
4445 ")\n"
4446 )
4447 << "Coupled collapse failed." << nl
4448 << "Master: " << eIndex << " : " << mEdge << nl
4449 << "Slave: " << sIndex << " : " << sEdge << nl
4450 << abort(FatalError);
4451 }
4453 break;
4454 }
4456 case 3:
4457 {
4458 if (sIndex < 0)
4459 {
4460 slaveMoveNewPoint[slaveI] =
4461 (
4462 0.5 * (points_[mEdge[0]] + points_[mEdge[1]])
4463 );
4465 slaveMoveOldPoint[slaveI] =
4466 (
4467 0.5 * (oldPoints_[mEdge[0]] + oldPoints_[mEdge[1]])
4468 );
4469 }
4470 else
4471 {
4472 slaveOverRide = 3;
4473 }
4475 break;
4476 }
4477 }
4479 // Temporarily turn off coupledModification
4480 unsetCoupledModification();
4482 // Test the slave edge
4483 if (localCouple)
4484 {
4485 slaveMap = collapseEdge(sIndex, slaveOverRide, true, forceOp);
4486 }
4487 else
4488 if (procCouple)
4489 {
4490 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
4492 if (sIndex < 0)
4493 {
4494 // Point-based coupling
4495 scalar slaveCollapseQuality(GREAT);
4496 dynamicLabelList cellsChecked(10);
4498 // Check cells connected to coupled point
4499 const labelList& pEdges = sMesh.pointEdges_[mag(sIndex)];
4501 bool infeasible = false;
4503 forAll(pEdges, edgeI)
4504 {
4505 const labelList& eFaces =
4506 (
4507 sMesh.edgeFaces_[pEdges[edgeI]]
4508 );
4510 // Build a list of cells to check
4511 forAll(eFaces, faceI)
4512 {
4513 label own = sMesh.owner_[eFaces[faceI]];
4514 label nei = sMesh.neighbour_[eFaces[faceI]];
4516 // Check owner cell
4517 if (findIndex(cellsChecked, own) == -1)
4518 {
4519 // Check if point movement is feasible
4520 if
4521 (
4522 sMesh.checkCollapse
4523 (
4524 slaveMoveNewPoint[slaveI],
4525 slaveMoveOldPoint[slaveI],
4526 mag(sIndex),
4527 own,
4528 cellsChecked,
4529 slaveCollapseQuality,
4530 forceOp
4531 )
4532 )
4533 {
4534 infeasible = true;
4535 break;
4536 }
4537 }
4539 if (nei == -1)
4540 {
4541 continue;
4542 }
4544 // Check neighbour cell
4545 if (findIndex(cellsChecked, nei) == -1)
4546 {
4547 // Check if point movement is feasible
4548 if
4549 (
4550 sMesh.checkCollapse
4551 (
4552 slaveMoveNewPoint[slaveI],
4553 slaveMoveOldPoint[slaveI],
4554 mag(sIndex),
4555 nei,
4556 cellsChecked,
4557 slaveCollapseQuality,
4558 forceOp
4559 )
4560 )
4561 {
4562 infeasible = true;
4563 break;
4564 }
4565 }
4566 }
4568 if (infeasible)
4569 {
4570 break;
4571 }
4572 }
4574 if (infeasible)
4575 {
4576 slaveMap.type() = 0;
4577 }
4578 else
4579 {
4580 slaveMap.type() = 1;
4581 }
4582 }
4583 else
4584 {
4585 // Edge-based coupling
4586 slaveMap =
4587 (
4588 sMesh.collapseEdge
4589 (
4590 sIndex,
4591 slaveOverRide,
4592 true,
4593 forceOp
4594 )
4595 );
4596 }
4597 }
4599 // Turn it back on.
4600 setCoupledModification();
4602 if (slaveMap.type() <= 0)
4603 {
4604 // Slave couldn't perform collapse.
4605 map.type() = -2;
4607 return map;
4608 }
4610 // Save index and patch for posterity
4611 slaveMap.index() = sIndex;
4612 slaveMap.patchIndex() = pI;
4613 }
4615 // Next collapse each slave edge (for real this time...)
4616 forAll(slaveMaps, slaveI)
4617 {
4618 // Alias for convenience...
4619 changeMap& slaveMap = slaveMaps[slaveI];
4621 label sIndex = slaveMap.index();
4622 label pI = slaveMap.patchIndex();
4624 // If the edge is mapped onto itself, skip modification
4625 // (occurs for cyclic edges)
4626 if ((sIndex == eIndex) && localCouple)
4627 {
4628 continue;
4629 }
4631 label slaveOverRide = slaveMap.type();
4633 // Temporarily turn off coupledModification
4634 unsetCoupledModification();
4636 // Collapse the slave
4637 if (localCouple)
4638 {
4639 slaveMap = collapseEdge(sIndex, slaveOverRide, false, forceOp);
4640 }
4641 else
4642 {
4643 const coupleMap& cMap = recvMeshes_[pI].map();
4644 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
4646 if (sIndex < 0)
4647 {
4648 // Point based coupling
4650 // Move points to new location,
4651 // and update operation into coupleMap
4652 sMesh.points_[mag(sIndex)] = slaveMoveNewPoint[slaveI];
4653 sMesh.oldPoints_[mag(sIndex)] = slaveMoveOldPoint[slaveI];
4655 cMap.pushOperation
4656 (
4657 mag(sIndex),
4658 coupleMap::MOVE_POINT,
4659 slaveMoveNewPoint[slaveI],
4660 slaveMoveOldPoint[slaveI]
4661 );
4663 // Force operation to succeed
4664 slaveMap.type() = 1;
4665 }
4666 else
4667 {
4668 // Edge-based coupling
4669 slaveMap =
4670 (
4671 sMesh.collapseEdge
4672 (
4673 sIndex,
4674 slaveOverRide,
4675 false,
4676 forceOp
4677 )
4678 );
4680 // Push operation into coupleMap
4681 switch (slaveMap.type())
4682 {
4683 case 1:
4684 {
4685 cMap.pushOperation
4686 (
4687 sIndex,
4688 coupleMap::COLLAPSE_FIRST
4689 );
4691 break;
4692 }
4694 case 2:
4695 {
4696 cMap.pushOperation
4697 (
4698 sIndex,
4699 coupleMap::COLLAPSE_SECOND
4700 );
4702 break;
4703 }
4705 case 3:
4706 {
4707 cMap.pushOperation
4708 (
4709 sIndex,
4710 coupleMap::COLLAPSE_MIDPOINT
4711 );
4713 break;
4714 }
4715 }
4716 }
4717 }
4719 // Turn it back on.
4720 setCoupledModification();
4722 // The final operation has to succeed.
4723 if (slaveMap.type() <= 0)
4724 {
4725 FatalErrorIn
4726 (
4727 "\n"
4728 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
4729 "(\n"
4730 " const label eIndex,\n"
4731 " label overRideCase,\n"
4732 " bool checkOnly,\n"
4733 " bool forceOp\n"
4734 ")\n"
4735 )
4736 << "Coupled topo-change for slave failed." << nl
4737 << " Edge: " << eIndex << ": " << eCheck << nl
4738 << " Slave index: " << sIndex << nl
4739 << " Patch index: " << pI << nl
4740 << " Type: " << slaveMap.type() << nl
4741 << abort(FatalError);
4742 }
4744 // Save index and patch for posterity
4745 slaveMap.index() = sIndex;
4746 slaveMap.patchIndex() = pI;
4747 }
4748 }
4750 // Build hullVertices for this edge
4751 labelList vertexHull;
4752 buildVertexHull(eIndex, vertexHull);
4754 // Hull variables
4755 bool found = false;
4756 label replaceIndex = -1, m = vertexHull.size();
4758 // Size up the hull lists
4759 labelList cellHull(m, -1);
4760 labelList faceHull(m, -1);
4761 labelList edgeHull(m, -1);
4762 labelListList ringEntities(4, labelList(m, -1));
4764 // Construct a hull around this edge
4765 meshOps::constructHull
4766 (
4767 eIndex,
4768 faces_,
4769 edges_,
4770 cells_,
4771 owner_,
4772 neighbour_,
4773 faceEdges_,
4774 edgeFaces_,
4775 vertexHull,
4776 edgeHull,
4777 faceHull,
4778 cellHull,
4779 ringEntities
4780 );
4782 // Check whether points of the edge lies on a boundary
4783 const FixedList<bool,2> edgeBoundary = checkEdgeBoundary(eIndex);
4784 FixedList<label, 2> nBoundCurves(0), nProcCurves(0), checkPoints(-1);
4786 // Decide on collapseCase
4787 label collapseCase = -1;
4789 if (edgeBoundary[0] && !edgeBoundary[1])
4790 {
4791 collapseCase = 1;
4792 }
4793 else
4794 if (!edgeBoundary[0] && edgeBoundary[1])
4795 {
4796 collapseCase = 2;
4797 }
4798 else
4799 if (edgeBoundary[0] && edgeBoundary[1])
4800 {
4801 // If this is an interior edge with two boundary points.
4802 // Bail out for now. If proximity based refinement is
4803 // switched on, mesh may be sliced at this point.
4804 label edgePatch = whichEdgePatch(eIndex);
4806 if (edgePatch == -1)
4807 {
4808 return map;
4809 }
4811 // Check if either point lies on a bounding curve
4812 // Used to ensure that collapses happen towards bounding curves.
4813 // If the edge itself is on a bounding curve, collapse is valid.
4814 const edge& edgeCheck = edges_[eIndex];
4816 forAll(edgeCheck, pointI)
4817 {
4818 const labelList& pEdges = pointEdges_[edgeCheck[pointI]];
4820 forAll(pEdges, edgeI)
4821 {
4822 if
4823 (
4824 checkBoundingCurve
4825 (
4826 pEdges[edgeI],
4827 false,
4828 &(nProcCurves[pointI])
4829 )
4830 )
4831 {
4832 nBoundCurves[pointI]++;
4833 }
4834 }
4835 }
4837 // Check for procCurves first
4838 if (!coupledModification_ && !isSubMesh_)
4839 {
4840 // Bail out for now
4841 if (nProcCurves[0] && nProcCurves[1])
4842 {
4843 return map;
4844 }
4846 if (nProcCurves[0] && !nProcCurves[1])
4847 {
4848 if (nBoundCurves[1])
4849 {
4850 // Bail out for now
4851 return map;
4852 }
4853 else
4854 {
4855 collapseCase = 1;
4856 }
4857 }
4859 if (!nProcCurves[0] && nProcCurves[1])
4860 {
4861 if (nBoundCurves[0])
4862 {
4863 // Bail out for now
4864 return map;
4865 }
4866 else
4867 {
4868 collapseCase = 2;
4869 }
4870 }
4871 }
4873 // If still undecided, choose based on bounding curves
4874 if (collapseCase == -1)
4875 {
4876 // Pick the point which is connected to more bounding curves
4877 if (nBoundCurves[0] > nBoundCurves[1])
4878 {
4879 collapseCase = 1;
4880 }
4881 else
4882 if (nBoundCurves[1] > nBoundCurves[0])
4883 {
4884 collapseCase = 2;
4885 }
4886 else
4887 {
4888 // Bounding edge: collapseEdge can collapse this edge
4889 collapseCase = 3;
4890 }
4891 }
4892 }
4893 else
4894 {
4895 // Looks like this is an interior edge.
4896 // Collapse case [3] by default
4897 collapseCase = 3;
4898 }
4900 // Perform an override if requested.
4901 if (overRideCase != -1)
4902 {
4903 collapseCase = overRideCase;
4904 }
4906 // Configure the new point-position
4907 point newPoint = vector::zero;
4908 point oldPoint = vector::zero;
4910 label collapsePoint = -1, replacePoint = -1;
4912 switch (collapseCase)
4913 {
4914 case 1:
4915 {
4916 // Collapse to the first node
4917 replacePoint = edges_[eIndex][0];
4918 collapsePoint = edges_[eIndex][1];
4920 newPoint = points_[replacePoint];
4921 oldPoint = oldPoints_[replacePoint];
4923 checkPoints[0] = collapsePoint;
4925 break;
4926 }
4928 case 2:
4929 {
4930 // Collapse to the second node
4931 replacePoint = edges_[eIndex][1];
4932 collapsePoint = edges_[eIndex][0];
4934 newPoint = points_[replacePoint];
4935 oldPoint = oldPoints_[replacePoint];
4937 checkPoints[0] = collapsePoint;
4939 break;
4940 }
4942 case 3:
4943 {
4944 // Collapse to the mid-point
4945 replacePoint = edges_[eIndex][1];
4946 collapsePoint = edges_[eIndex][0];
4948 newPoint =
4949 (
4950 0.5 *
4951 (
4952 points_[replacePoint]
4953 + points_[collapsePoint]
4954 )
4955 );
4957 oldPoint =
4958 (
4959 0.5 *
4960 (
4961 oldPoints_[replacePoint]
4962 + oldPoints_[collapsePoint]
4963 )
4964 );
4966 checkPoints[0] = replacePoint;
4967 checkPoints[1] = collapsePoint;
4969 break;
4970 }
4972 default:
4973 {
4974 // Don't think this will ever happen.
4975 FatalErrorIn
4976 (
4977 "\n"
4978 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
4979 "(\n"
4980 " const label eIndex,\n"
4981 " label overRideCase,\n"
4982 " bool checkOnly,\n"
4983 " bool forceOp\n"
4984 ")\n"
4985 )
4986 << "Edge: " << eIndex << ": " << edges_[eIndex]
4987 << ". Couldn't decide on collapseCase."
4988 << abort(FatalError);
4990 break;
4991 }
4992 }
4994 // Loop through edges and check for feasibility of collapse
4995 // Also, keep track of resulting cell quality,
4996 // if collapse is indeed feasible
4997 scalar collapseQuality(GREAT);
4998 dynamicLabelList cellsChecked(10);
5000 // Add all hull cells as 'checked',
5001 // and therefore, feasible
5002 forAll(cellHull, cellI)
5003 {
5004 if (cellHull[cellI] == -1)
5005 {
5006 continue;
5007 }
5009 cellsChecked.append(cellHull[cellI]);
5010 }
5012 // Check collapsibility of cells around edges
5013 // with the re-configured point
5014 forAll(checkPoints, pointI)
5015 {
5016 if (checkPoints[pointI] == -1)
5017 {
5018 continue;
5019 }
5021 const labelList& checkPointEdges = pointEdges_[checkPoints[pointI]];
5023 forAll(checkPointEdges, edgeI)
5024 {
5025 const labelList& eFaces = edgeFaces_[checkPointEdges[edgeI]];
5027 // Build a list of cells to check
5028 forAll(eFaces, faceI)
5029 {
5030 label own = owner_[eFaces[faceI]];
5031 label nei = neighbour_[eFaces[faceI]];
5033 // Check owner cell
5034 if (findIndex(cellsChecked, own) == -1)
5035 {
5036 // Check if a collapse is feasible
5037 if
5038 (
5039 checkCollapse
5040 (
5041 newPoint,
5042 oldPoint,
5043 checkPoints[pointI],
5044 own,
5045 cellsChecked,
5046 collapseQuality,
5047 forceOp
5048 )
5049 )
5050 {
5051 map.type() = 0;
5052 return map;
5053 }
5054 }
5056 if (nei == -1)
5057 {
5058 continue;
5059 }
5061 // Check neighbour cell
5062 if (findIndex(cellsChecked, nei) == -1)
5063 {
5064 // Check if a collapse is feasible
5065 if
5066 (
5067 checkCollapse
5068 (
5069 newPoint,
5070 oldPoint,
5071 checkPoints[pointI],
5072 nei,
5073 cellsChecked,
5074 collapseQuality,
5075 forceOp
5076 )
5077 )
5078 {
5079 map.type() = 0;
5080 return map;
5081 }
5082 }
5083 }
5084 }
5085 }
5087 // Add a map entry of the replacePoint as an 'addedPoint'
5088 // - Used in coupled mapping
5089 map.addPoint(replacePoint);
5091 // Are we only performing checks?
5092 if (checkOnly)
5093 {
5094 // Return a succesful collapse possibility
5095 map.type() = collapseCase;
5097 // Make note of the removed point
5098 map.removePoint(collapsePoint);
5100 if (debug > 2)
5101 {
5102 Pout<< nl << "Edge: " << eIndex
5103 << ":: " << edges_[eIndex] << nl
5104 << " collapseCase determined to be: " << collapseCase << nl
5105 << " Resulting quality: " << collapseQuality << nl
5106 << " collapsePoint: " << collapsePoint << nl
5107 << " nBoundCurves: " << nBoundCurves << nl
5108 << " nProcCurves: " << nProcCurves << nl
5109 << endl;
5110 }
5112 return map;
5113 }
5115 // Define indices on the hull for removal / replacement
5116 label removeEdgeIndex = -1, replaceEdgeIndex = -1;
5117 label removeFaceIndex = -1, replaceFaceIndex = -1;
5119 if (replacePoint == edges_[eIndex][0])
5120 {
5121 replaceEdgeIndex = 0;
5122 replaceFaceIndex = 1;
5123 removeEdgeIndex = 2;
5124 removeFaceIndex = 3;
5125 }
5126 else
5127 if (replacePoint == edges_[eIndex][1])
5128 {
5129 removeEdgeIndex = 0;
5130 removeFaceIndex = 1;
5131 replaceEdgeIndex = 2;
5132 replaceFaceIndex = 3;
5133 }
5134 else
5135 {
5136 // Don't think this will ever happen.
5137 FatalErrorIn
5138 (
5139 "\n"
5140 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
5141 "(\n"
5142 " const label eIndex,\n"
5143 " label overRideCase,\n"
5144 " bool checkOnly,\n"
5145 " bool forceOp\n"
5146 ")\n"
5147 )
5148 << "Edge: " << eIndex << ": " << edges_[eIndex]
5149 << ". Couldn't decide on removal / replacement indices."
5150 << abort(FatalError);
5151 }
5153 if (debug > 1)
5154 {
5155 Pout<< nl << nl
5156 << "Edge: " << eIndex << ": " << edges_[eIndex]
5157 << " is to be collapsed. " << nl;
5159 Pout<< " On SubMesh: " << isSubMesh_ << nl;
5160 Pout<< " coupledModification: " << coupledModification_ << nl;
5162 label epIndex = whichEdgePatch(eIndex);
5164 const polyBoundaryMesh& boundary = boundaryMesh();
5166 Pout<< " Patch: ";
5168 if (epIndex == -1)
5169 {
5170 Pout<< "Internal" << nl;
5171 }
5172 else
5173 if (epIndex < boundary.size())
5174 {
5175 Pout<< boundary[epIndex].name() << nl;
5176 }
5177 else
5178 {
5179 Pout<< " New patch: " << epIndex << endl;
5180 }
5182 Pout<< " nBoundCurves: " << nBoundCurves << nl
5183 << " nProcCurves: " << nProcCurves << nl
5184 << " collapseCase: " << collapseCase << nl
5185 << " Resulting quality: " << collapseQuality << endl;
5187 if (debug > 2)
5188 {
5189 Pout<< " Edges: " << edgeHull << nl
5190 << " Faces: " << faceHull << nl
5191 << " Cells: " << cellHull << nl
5192 << " replacePoint: " << replacePoint << nl
5193 << " collapsePoint: " << collapsePoint << nl
5194 << " checkPoints: " << checkPoints << nl;
5196 Pout<< " ringEntities (removed faces): " << nl;
5198 forAll(ringEntities[removeFaceIndex], faceI)
5199 {
5200 label fIndex = ringEntities[removeFaceIndex][faceI];
5202 if (fIndex == -1)
5203 {
5204 continue;
5205 }
5207 Pout<< fIndex << ": " << faces_[fIndex] << nl;
5208 }
5210 Pout<< " ringEntities (removed edges): " << nl;
5212 forAll(ringEntities[removeEdgeIndex], edgeI)
5213 {
5214 label ieIndex = ringEntities[removeEdgeIndex][edgeI];
5216 if (ieIndex == -1)
5217 {
5218 continue;
5219 }
5221 Pout<< ieIndex << ": " << edges_[ieIndex] << nl;
5222 }
5224 Pout<< " ringEntities (replacement faces): " << nl;
5226 forAll(ringEntities[replaceFaceIndex], faceI)
5227 {
5228 label fIndex = ringEntities[replaceFaceIndex][faceI];
5230 if (fIndex == -1)
5231 {
5232 continue;
5233 }
5235 Pout<< fIndex << ": " << faces_[fIndex] << nl;
5236 }
5238 Pout<< " ringEntities (replacement edges): " << nl;
5240 forAll(ringEntities[replaceEdgeIndex], edgeI)
5241 {
5242 label ieIndex = ringEntities[replaceEdgeIndex][edgeI];
5244 if (ieIndex == -1)
5245 {
5246 continue;
5247 }
5249 Pout<< ieIndex << ": " << edges_[ieIndex] << nl;
5250 }
5252 labelList& collapsePointEdges = pointEdges_[collapsePoint];
5254 Pout<< " pointEdges (collapsePoint): ";
5256 forAll(collapsePointEdges, edgeI)
5257 {
5258 Pout<< collapsePointEdges[edgeI] << " ";
5259 }
5261 Pout<< nl;
5263 // Write out VTK files prior to change
5264 // - Using old-points for convenience in post-processing
5265 if (debug > 3)
5266 {
5267 writeVTK
5268 (
5269 Foam::name(eIndex)
5270 + '(' + Foam::name(collapsePoint)
5271 + ',' + Foam::name(replacePoint) + ')'
5272 + "_Collapse_0",
5273 cellsChecked,
5274 3, false, true
5275 );
5276 }
5277 }
5278 }
5280 if (whichEdgePatch(eIndex) > -1)
5281 {
5282 // Update number of surface collapses, if necessary.
5283 status(SURFACE_COLLAPSES)++;
5284 }
5286 // Maintain a list of modified faces for mapping
5287 dynamicLabelList modifiedFaces(10);
5289 // Renumber all hull faces and edges
5290 forAll(faceHull, indexI)
5291 {
5292 // Loop through all faces of the edge to be removed
5293 // and reassign them to the replacement edge
5294 label edgeToRemove = ringEntities[removeEdgeIndex][indexI];
5295 label faceToRemove = ringEntities[removeFaceIndex][indexI];
5296 label cellToRemove = cellHull[indexI];
5297 label replaceEdge = ringEntities[replaceEdgeIndex][indexI];
5298 label replaceFace = ringEntities[replaceFaceIndex][indexI];
5300 label replacePatch = whichEdgePatch(replaceEdge);
5301 label removePatch = whichEdgePatch(edgeToRemove);
5303 // Check if a patch conversion is necessary
5304 if (replacePatch == -1 && removePatch > -1)
5305 {
5306 if (debug > 2)
5307 {
5308 Pout<< " Edge: " << replaceEdge
5309 << " :: " << edges_[replaceEdge]
5310 << " is interior, but edge: " << edgeToRemove
5311 << " :: " << edges_[edgeToRemove]
5312 << " is on a boundary patch."
5313 << endl;
5314 }
5316 // Convert patch for edge
5317 edge newEdge = edges_[replaceEdge];
5318 labelList newEdgeFaces = edgeFaces_[replaceEdge];
5320 // Insert the new edge
5321 label newEdgeIndex =
5322 (
5323 insertEdge
5324 (
5325 removePatch,
5326 newEdge,
5327 newEdgeFaces
5328 )
5329 );
5331 // Replace faceEdges for all
5332 // connected faces.
5333 forAll(newEdgeFaces, faceI)
5334 {
5335 meshOps::replaceLabel
5336 (
5337 replaceEdge,
5338 newEdgeIndex,
5339 faceEdges_[newEdgeFaces[faceI]]
5340 );
5341 }
5343 // Remove the edge
5344 removeEdge(replaceEdge);
5346 // Update map
5347 map.removeEdge(replaceEdge);
5348 map.addEdge(newEdgeIndex, labelList(1, replaceEdge));
5350 // Replace index and patch
5351 replaceEdge = newEdgeIndex;
5353 // Modify ringEntities
5354 ringEntities[replaceEdgeIndex][indexI] = newEdgeIndex;
5355 }
5357 const labelList& rmvEdgeFaces = edgeFaces_[edgeToRemove];
5359 forAll(rmvEdgeFaces, faceI)
5360 {
5361 // Replace edge labels for faces
5362 meshOps::replaceLabel
5363 (
5364 edgeToRemove,
5365 replaceEdge,
5366 faceEdges_[rmvEdgeFaces[faceI]]
5367 );
5369 // Loop through faces associated with this edge,
5370 // and renumber them as well.
5371 const face& faceToCheck = faces_[rmvEdgeFaces[faceI]];
5373 if ((replaceIndex = faceToCheck.which(collapsePoint)) > -1)
5374 {
5375 // Renumber the face...
5376 faces_[rmvEdgeFaces[faceI]][replaceIndex] = replacePoint;
5378 // Add an entry for mapping
5379 if (findIndex(modifiedFaces, rmvEdgeFaces[faceI]) == -1)
5380 {
5381 modifiedFaces.append(rmvEdgeFaces[faceI]);
5382 }
5383 }
5385 // Hull faces should be removed for the replacement edge
5386 if (rmvEdgeFaces[faceI] == faceHull[indexI])
5387 {
5388 meshOps::sizeDownList
5389 (
5390 faceHull[indexI],
5391 edgeFaces_[replaceEdge]
5392 );
5394 continue;
5395 }
5397 found = false;
5399 // Need to avoid ring faces as well.
5400 forAll(ringEntities[removeFaceIndex], faceII)
5401 {
5402 if
5403 (
5404 rmvEdgeFaces[faceI]
5405 == ringEntities[removeFaceIndex][faceII]
5406 )
5407 {
5408 found = true;
5409 break;
5410 }
5411 }
5413 // Size-up the replacement edge list if the face hasn't been found.
5414 // These faces are connected to the edge slated for
5415 // removal, but do not belong to the hull.
5416 if (!found)
5417 {
5418 meshOps::sizeUpList
5419 (
5420 rmvEdgeFaces[faceI],
5421 edgeFaces_[replaceEdge]
5422 );
5423 }
5424 }
5426 if (cellToRemove == -1)
5427 {
5428 continue;
5429 }
5431 // Size down edgeFaces for the ring edges
5432 meshOps::sizeDownList
5433 (
5434 faceToRemove,
5435 edgeFaces_[edgeHull[indexI]]
5436 );
5438 // Ensure proper orientation of retained faces
5439 if (owner_[faceToRemove] == cellToRemove)
5440 {
5441 if (owner_[replaceFace] == cellToRemove)
5442 {
5443 if
5444 (
5445 (neighbour_[faceToRemove] > neighbour_[replaceFace])
5446 && (neighbour_[replaceFace] != -1)
5447 )
5448 {
5449 // This face is to be flipped
5450 faces_[replaceFace] = faces_[replaceFace].reverseFace();
5451 owner_[replaceFace] = neighbour_[replaceFace];
5452 neighbour_[replaceFace] = neighbour_[faceToRemove];
5454 setFlip(replaceFace);
5455 }
5456 else
5457 if
5458 (
5459 (neighbour_[faceToRemove] == -1) &&
5460 (neighbour_[replaceFace] != -1)
5461 )
5462 {
5463 // This interior face would need to be converted
5464 // to a boundary one, and flipped as well.
5465 face newFace = faces_[replaceFace].reverseFace();
5466 label newOwner = neighbour_[replaceFace];
5467 label newNeighbour = neighbour_[faceToRemove];
5468 labelList newFE = faceEdges_[replaceFace];
5470 label newFaceIndex =
5471 (
5472 insertFace
5473 (
5474 whichPatch(faceToRemove),
5475 newFace,
5476 newOwner,
5477 newNeighbour,
5478 newFE
5479 )
5480 );
5482 // Set this face aside for mapping
5483 modifiedFaces.append(newFaceIndex);
5485 // Update map.
5486 map.addFace(newFaceIndex, labelList(1, faceToRemove));
5488 // Ensure that all edges of this face are
5489 // on the boundary.
5490 forAll(newFE, edgeI)
5491 {
5492 if (whichEdgePatch(newFE[edgeI]) == -1)
5493 {
5494 edge newEdge = edges_[newFE[edgeI]];
5495 labelList newEF = edgeFaces_[newFE[edgeI]];
5497 // Need patch information for the new edge.
5498 // Find the corresponding edge in ringEntities.
5499 // Note that hullEdges doesn't need to be checked,
5500 // since they are common to both faces.
5501 label i =
5502 (
5503 findIndex
5504 (
5505 ringEntities[replaceEdgeIndex],
5506 newFE[edgeI]
5507 )
5508 );
5510 label repIndex =
5511 (
5512 whichEdgePatch
5513 (
5514 ringEntities[removeEdgeIndex][i]
5515 )
5516 );
5518 // Insert the new edge
5519 label newEdgeIndex =
5520 (
5521 insertEdge
5522 (
5523 repIndex,
5524 newEdge,
5525 newEF
5526 )
5527 );
5529 // Replace faceEdges for all
5530 // connected faces.
5531 forAll(newEF, faceI)
5532 {
5533 meshOps::replaceLabel
5534 (
5535 newFE[edgeI],
5536 newEdgeIndex,
5537 faceEdges_[newEF[faceI]]
5538 );
5539 }
5541 // Remove the edge
5542 removeEdge(newFE[edgeI]);
5544 // Update map
5545 map.removeEdge(newFE[edgeI]);
5547 map.addEdge
5548 (
5549 newEdgeIndex,
5550 labelList(1, newFE[edgeI])
5551 );
5553 // Replace faceEdges with new edge index
5554 newFE[edgeI] = newEdgeIndex;
5556 // Modify ringEntities
5557 ringEntities[replaceEdgeIndex][i] = newEdgeIndex;
5558 }
5559 }
5561 // Add the new faceEdges
5562 faceEdges_[newFaceIndex] = newFE;
5564 // Replace edgeFaces with the new face index
5565 const labelList& newFEdges = faceEdges_[newFaceIndex];
5567 forAll(newFEdges, edgeI)
5568 {
5569 meshOps::replaceLabel
5570 (
5571 replaceFace,
5572 newFaceIndex,
5573 edgeFaces_[newFEdges[edgeI]]
5574 );
5575 }
5577 // Remove the face
5578 removeFace(replaceFace);
5580 // Update map
5581 map.removeFace(replaceFace);
5583 // Replace label for the new owner
5584 meshOps::replaceLabel
5585 (
5586 replaceFace,
5587 newFaceIndex,
5588 cells_[newOwner]
5589 );
5591 // Modify ringEntities and replaceFace
5592 replaceFace = newFaceIndex;
5593 ringEntities[replaceFaceIndex][indexI] = newFaceIndex;
5594 }
5595 else
5596 if
5597 (
5598 (neighbour_[faceToRemove] == -1) &&
5599 (neighbour_[replaceFace] == -1)
5600 )
5601 {
5602 // Wierd overhanging cell. Since replaceFace
5603 // would be an orphan if this continued, remove
5604 // the face entirely.
5605 labelList rmFE = faceEdges_[replaceFace];
5607 forAll(rmFE, edgeI)
5608 {
5609 if
5610 (
5611 (edgeFaces_[rmFE[edgeI]].size() == 1) &&
5612 (edgeFaces_[rmFE[edgeI]][0] == replaceFace)
5613 )
5614 {
5615 // This edge has to be removed entirely.
5616 removeEdge(rmFE[edgeI]);
5618 // Update map
5619 map.removeEdge(rmFE[edgeI]);
5621 label i =
5622 (
5623 findIndex
5624 (
5625 ringEntities[replaceEdgeIndex],
5626 rmFE[edgeI]
5627 )
5628 );
5630 if (i > -1)
5631 {
5632 // Modify ringEntities
5633 ringEntities[replaceEdgeIndex][i] = -1;
5634 }
5635 }
5636 else
5637 {
5638 // Size-down edgeFaces
5639 meshOps::sizeDownList
5640 (
5641 replaceFace,
5642 edgeFaces_[rmFE[edgeI]]
5643 );
5644 }
5645 }
5647 // If this is a subMesh, and replaceFace is on
5648 // a physical boundary, make a 'special' entry
5649 // for coupled mapping purposes.
5650 if (isSubMesh_)
5651 {
5652 label kfPatch = whichPatch(replaceFace);
5653 label rfPatch = whichPatch(faceToRemove);
5655 if
5656 (
5657 (getNeighbourProcessor(rfPatch) > -1) &&
5658 (getNeighbourProcessor(kfPatch) == -1)
5659 )
5660 {
5661 map.addFace
5662 (
5663 faceToRemove,
5664 labelList(1, (-2 - kfPatch))
5665 );
5666 }
5667 }
5669 // Remove the face
5670 removeFace(replaceFace);
5672 // Update map
5673 map.removeFace(replaceFace);
5675 // Modify ringEntities and replaceFace
5676 replaceFace = -1;
5677 ringEntities[replaceFaceIndex][indexI] = -1;
5678 }
5679 else
5680 {
5681 // Keep orientation intact, and update the owner
5682 owner_[replaceFace] = neighbour_[faceToRemove];
5683 }
5684 }
5685 else
5686 if (neighbour_[faceToRemove] == -1)
5687 {
5688 // This interior face would need to be converted
5689 // to a boundary one, but with orientation intact.
5690 face newFace = faces_[replaceFace];
5691 label newOwner = owner_[replaceFace];
5692 label newNeighbour = neighbour_[faceToRemove];
5693 labelList newFE = faceEdges_[replaceFace];
5695 label newFaceIndex =
5696 (
5697 insertFace
5698 (
5699 whichPatch(faceToRemove),
5700 newFace,
5701 newOwner,
5702 newNeighbour,
5703 newFE
5704 )
5705 );
5707 // Set this face aside for mapping
5708 modifiedFaces.append(newFaceIndex);
5710 // Update map
5711 map.addFace(newFaceIndex, labelList(1, faceToRemove));
5713 // Ensure that all edges of this face are
5714 // on the boundary.
5715 forAll(newFE, edgeI)
5716 {
5717 if (whichEdgePatch(newFE[edgeI]) == -1)
5718 {
5719 edge newEdge = edges_[newFE[edgeI]];
5720 labelList newEF = edgeFaces_[newFE[edgeI]];
5722 // Need patch information for the new edge.
5723 // Find the corresponding edge in ringEntities.
5724 // Note that hullEdges doesn't need to be checked,
5725 // since they are common to both faces.
5726 label i =
5727 (
5728 findIndex
5729 (
5730 ringEntities[replaceEdgeIndex],
5731 newFE[edgeI]
5732 )
5733 );
5735 label repIndex =
5736 (
5737 whichEdgePatch
5738 (
5739 ringEntities[removeEdgeIndex][i]
5740 )
5741 );
5743 // Insert the new edge
5744 label newEdgeIndex =
5745 (
5746 insertEdge
5747 (
5748 repIndex,
5749 newEdge,
5750 newEF
5751 )
5752 );
5754 // Replace faceEdges for all
5755 // connected faces.
5756 forAll(newEF, faceI)
5757 {
5758 meshOps::replaceLabel
5759 (
5760 newFE[edgeI],
5761 newEdgeIndex,
5762 faceEdges_[newEF[faceI]]
5763 );
5764 }
5766 // Remove the edge
5767 removeEdge(newFE[edgeI]);
5769 // Update map
5770 map.removeEdge(newFE[edgeI]);
5772 map.addEdge
5773 (
5774 newEdgeIndex,
5775 labelList(1, newFE[edgeI])
5776 );
5778 // Replace faceEdges with new edge index
5779 newFE[edgeI] = newEdgeIndex;
5781 // Modify ringEntities
5782 ringEntities[replaceEdgeIndex][i] = newEdgeIndex;
5783 }
5784 }
5786 // Add the new faceEdges
5787 faceEdges_[newFaceIndex] = newFE;
5789 // Replace edgeFaces with the new face index
5790 const labelList& newFEdges = faceEdges_[newFaceIndex];
5792 forAll(newFEdges, edgeI)
5793 {
5794 meshOps::replaceLabel
5795 (
5796 replaceFace,
5797 newFaceIndex,
5798 edgeFaces_[newFEdges[edgeI]]
5799 );
5800 }
5802 // Remove the face
5803 removeFace(replaceFace);
5805 // Update map
5806 map.removeFace(replaceFace);
5808 // Replace label for the new owner
5809 meshOps::replaceLabel
5810 (
5811 replaceFace,
5812 newFaceIndex,
5813 cells_[newOwner]
5814 );
5816 // Modify ringEntities and replaceFace
5817 replaceFace = newFaceIndex;
5818 ringEntities[replaceFaceIndex][indexI] = newFaceIndex;
5819 }
5820 else
5821 {
5822 // Keep orientation intact, and update the neighbour
5823 neighbour_[replaceFace] = neighbour_[faceToRemove];
5824 }
5826 // Update the cell
5827 if (neighbour_[faceToRemove] != -1)
5828 {
5829 meshOps::replaceLabel
5830 (
5831 faceToRemove,
5832 replaceFace,
5833 cells_[neighbour_[faceToRemove]]
5834 );
5835 }
5836 }
5837 else
5838 {
5839 if (neighbour_[replaceFace] == cellToRemove)
5840 {
5841 if (owner_[faceToRemove] < owner_[replaceFace])
5842 {
5843 // This face is to be flipped
5844 faces_[replaceFace] = faces_[replaceFace].reverseFace();
5845 neighbour_[replaceFace] = owner_[replaceFace];
5846 owner_[replaceFace] = owner_[faceToRemove];
5848 setFlip(replaceFace);
5849 }
5850 else
5851 {
5852 // Keep orientation intact, and update the neighbour
5853 neighbour_[replaceFace] = owner_[faceToRemove];
5854 }
5855 }
5856 else
5857 {
5858 // Keep orientation intact, and update the owner
5859 owner_[replaceFace] = owner_[faceToRemove];
5860 }
5862 // Update the cell
5863 meshOps::replaceLabel
5864 (
5865 faceToRemove,
5866 replaceFace,
5867 cells_[owner_[faceToRemove]]
5868 );
5869 }
5870 }
5872 // Remove all hull entities
5873 forAll(faceHull, indexI)
5874 {
5875 label edgeToRemove = ringEntities[removeEdgeIndex][indexI];
5876 label faceToRemove = ringEntities[removeFaceIndex][indexI];
5877 label cellToRemove = cellHull[indexI];
5879 if (cellToRemove != -1)
5880 {
5881 // Remove faceToRemove and associated faceEdges
5882 removeFace(faceToRemove);
5884 // Update map
5885 map.removeFace(faceToRemove);
5887 // Remove the hull cell
5888 removeCell(cellToRemove);
5890 // Update map
5891 map.removeCell(cellToRemove);
5892 }
5894 // Remove the hull edge and associated edgeFaces
5895 removeEdge(edgeToRemove);
5897 // Update map
5898 map.removeEdge(edgeToRemove);
5900 // Remove the hull face
5901 removeFace(faceHull[indexI]);
5903 // Update map
5904 map.removeFace(faceHull[indexI]);
5905 }
5907 // Loop through pointEdges for the collapsePoint,
5908 // and replace all occurrences with replacePoint.
5909 // Size-up pointEdges for the replacePoint as well.
5910 const labelList& pEdges = pointEdges_[collapsePoint];
5912 forAll(pEdges, edgeI)
5913 {
5914 // Renumber edges
5915 label edgeIndex = pEdges[edgeI];
5917 if (edgeIndex != eIndex)
5918 {
5919 if (edges_[edgeIndex][0] == collapsePoint)
5920 {
5921 edges_[edgeIndex][0] = replacePoint;
5923 meshOps::sizeUpList
5924 (
5925 edgeIndex,
5926 pointEdges_[replacePoint]
5927 );
5928 }
5929 else
5930 if (edges_[edgeIndex][1] == collapsePoint)
5931 {
5932 edges_[edgeIndex][1] = replacePoint;
5934 meshOps::sizeUpList
5935 (
5936 edgeIndex,
5937 pointEdges_[replacePoint]
5938 );
5939 }
5940 else
5941 {
5942 // Looks like pointEdges is inconsistent
5943 FatalErrorIn
5944 (
5945 "\n"
5946 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
5947 "(\n"
5948 " const label eIndex,\n"
5949 " label overRideCase,\n"
5950 " bool checkOnly,\n"
5951 " bool forceOp\n"
5952 ")\n"
5953 )
5954 << "pointEdges is inconsistent." << nl
5955 << "Point: " << collapsePoint << nl
5956 << "pointEdges: " << pEdges << nl
5957 << abort(FatalError);
5958 }
5960 // Loop through faces associated with this edge,
5961 // and renumber them as well.
5962 const labelList& eFaces = edgeFaces_[edgeIndex];
5964 forAll(eFaces, faceI)
5965 {
5966 const face& faceToCheck = faces_[eFaces[faceI]];
5968 if ((replaceIndex = faceToCheck.which(collapsePoint)) > -1)
5969 {
5970 // Renumber the face...
5971 faces_[eFaces[faceI]][replaceIndex] = replacePoint;
5973 // Set this face aside for mapping
5974 if (findIndex(modifiedFaces, eFaces[faceI]) == -1)
5975 {
5976 modifiedFaces.append(eFaces[faceI]);
5977 }
5978 }
5979 }
5980 }
5981 }
5983 // Set old / new points
5984 oldPoints_[replacePoint] = oldPoint;
5985 points_[replacePoint] = newPoint;
5987 // Remove the collapse point
5988 removePoint(collapsePoint);
5990 // Update map
5991 map.removePoint(collapsePoint);
5993 // Remove the edge
5994 removeEdge(eIndex);
5996 // Update map
5997 map.removeEdge(eIndex);
5999 // Check for duplicate edges connected to the replacePoint
6000 const labelList& rpEdges = pointEdges_[replacePoint];
6002 dynamicLabelList mergeFaces(10);
6004 forAll(rpEdges, edgeI)
6005 {
6006 const edge& eCheckI = edges_[rpEdges[edgeI]];
6007 const label vCheckI = eCheckI.otherVertex(replacePoint);
6009 for (label edgeJ = edgeI + 1; edgeJ < rpEdges.size(); edgeJ++)
6010 {
6011 const edge& eCheckJ = edges_[rpEdges[edgeJ]];
6012 const label vCheckJ = eCheckJ.otherVertex(replacePoint);
6014 if (vCheckJ == vCheckI)
6015 {
6016 labelPair efCheck(rpEdges[edgeI], rpEdges[edgeJ]);
6018 forAll(efCheck, edgeI)
6019 {
6020 const labelList& eF = edgeFaces_[efCheck[edgeI]];
6022 forAll(eF, faceI)
6023 {
6024 label patch = whichPatch(eF[faceI]);
6026 if (patch == -1)
6027 {
6028 continue;
6029 }
6031 if (findIndex(mergeFaces, eF[faceI]) == -1)
6032 {
6033 mergeFaces.append(eF[faceI]);
6034 }
6035 }
6036 }
6038 if (debug > 2)
6039 {
6040 Pout<< " Found duplicate: " << eCheckI << nl
6041 << " Merge faces: " << mergeFaces << nl
6042 << endl;
6043 }
6044 }
6045 }
6046 }
6048 // Merge faces if necessary
6049 if (mergeFaces.size())
6050 {
6051 mergeBoundaryFaces(mergeFaces);
6052 }
6054 // Check and remove edges with an empty edgeFaces list
6055 const labelList& replaceEdges = ringEntities[replaceEdgeIndex];
6057 forAll(replaceEdges, edgeI)
6058 {
6059 label replaceEdge = replaceEdges[edgeI];
6061 // If the ring edge was removed, don't bother.
6062 if (replaceEdge > -1)
6063 {
6064 // Account for merged edges as well
6065 if (edgeFaces_[replaceEdge].empty())
6066 {
6067 // Is this edge truly removed? If not, remove it.
6068 if (edges_[replaceEdge] != edge(-1, -1))
6069 {
6070 if (debug > 2)
6071 {
6072 Pout<< " Edge: " << replaceEdge
6073 << " :: " << edges_[replaceEdge]
6074 << " has empty edgeFaces."
6075 << endl;
6076 }
6078 // Remove the edge
6079 removeEdge(replaceEdge);
6081 // Update map
6082 map.removeEdge(replaceEdge);
6083 }
6084 }
6085 }
6086 }
6088 // Create a list of candidates for mapping
6089 // using the list of checked cells
6090 labelList mC(cellsChecked.size());
6092 // For cell-mapping, exclude all hull-cells
6093 forAll(cellsChecked, indexI)
6094 {
6095 mC[indexI] = cellsChecked[indexI];
6097 if (cells_[cellsChecked[indexI]].empty())
6098 {
6099 cellsChecked[indexI] = -1;
6100 }
6101 else
6102 if (findIndex(cellHull, cellsChecked[indexI]) > 0)
6103 {
6104 cellsChecked[indexI] = -1;
6105 }
6106 }
6108 // Write out VTK files after change
6109 // - Using old-points for convenience in post-processing
6110 if (debug > 3)
6111 {
6112 writeVTK
6113 (
6114 Foam::name(eIndex)
6115 + '(' + Foam::name(collapsePoint)
6116 + ',' + Foam::name(replacePoint) + ')'
6117 + "_Collapse_1",
6118 cellsChecked,
6119 3, false, true
6120 );
6121 }
6123 // Now that all old / new cells possess correct connectivity,
6124 // compute mapping information.
6125 forAll(cellsChecked, cellI)
6126 {
6127 if (cellsChecked[cellI] < 0)
6128 {
6129 continue;
6130 }
6132 // Set the mapping for this cell
6133 setCellMapping(cellsChecked[cellI], mC);
6134 }
6136 // Set face mapping information for modified faces
6137 forAll(modifiedFaces, faceI)
6138 {
6139 const label mfIndex = modifiedFaces[faceI];
6141 // Exclude deleted faces
6142 if (faces_[mfIndex].empty())
6143 {
6144 continue;
6145 }
6147 // Decide between default / weighted mapping
6148 // based on boundary information
6149 label fPatch = whichPatch(mfIndex);
6151 if (fPatch == -1)
6152 {
6153 setFaceMapping(mfIndex);
6154 }
6155 else
6156 {
6157 // Fill-in candidate mapping information
6158 labelList faceCandidates;
6160 const labelList& fEdges = faceEdges_[mfIndex];
6162 forAll(fEdges, edgeI)
6163 {
6164 if (whichEdgePatch(fEdges[edgeI]) == fPatch)
6165 {
6166 // Loop through associated edgeFaces
6167 const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
6169 forAll(eFaces, faceI)
6170 {
6171 if
6172 (
6173 (eFaces[faceI] != mfIndex) &&
6174 (whichPatch(eFaces[faceI]) == fPatch)
6175 )
6176 {
6177 faceCandidates.setSize
6178 (
6179 faceCandidates.size() + 1,
6180 eFaces[faceI]
6181 );
6182 }
6183 }
6184 }
6185 }
6187 if (faceCandidates.empty())
6188 {
6189 // Add the face itself
6190 faceCandidates.setSize(1, mfIndex);
6191 }
6193 // Set the mapping for this face
6194 setFaceMapping(mfIndex, faceCandidates);
6195 }
6196 }
6198 // If modification is coupled, update mapping info.
6199 if (coupledModification_)
6200 {
6201 // Build a list of boundary edges / faces for mapping
6202 dynamicLabelList checkEdges(10), checkFaces(10);
6204 const labelList& pEdges = pointEdges_[replacePoint];
6206 forAll(pEdges, edgeI)
6207 {
6208 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
6210 forAll(eFaces, faceI)
6211 {
6212 label fPatch = whichPatch(eFaces[faceI]);
6214 if (localCouple && !procCouple)
6215 {
6216 if (!locallyCoupledEntity(eFaces[faceI], true, false, true))
6217 {
6218 continue;
6219 }
6221 // Check for cyclics
6222 if (boundaryMesh()[fPatch].type() == "cyclic")
6223 {
6224 // Check if this is a master face
6225 const coupleMap& cM = patchCoupling_[fPatch].map();
6226 const Map<label>& fM = cM.entityMap(coupleMap::FACE);
6228 // Only add master faces
6229 // (belonging to the first half)
6230 if (!fM.found(eFaces[faceI]))
6231 {
6232 continue;
6233 }
6234 }
6235 }
6236 else
6237 if (procCouple && !localCouple)
6238 {
6239 if (getNeighbourProcessor(fPatch) == -1)
6240 {
6241 continue;
6242 }
6243 }
6245 // Add face and its edges for checking
6246 if (findIndex(checkFaces, eFaces[faceI]) == -1)
6247 {
6248 // Add this face
6249 checkFaces.append(eFaces[faceI]);
6251 const labelList& fEdges = faceEdges_[eFaces[faceI]];
6253 forAll(fEdges, edgeJ)
6254 {
6255 if (findIndex(checkEdges, fEdges[edgeJ]) == -1)
6256 {
6257 checkEdges.append(fEdges[edgeJ]);
6258 }
6259 }
6260 }
6261 }
6262 }
6264 // Prepare a checklist
6265 boolList matchedFaces(checkFaces.size(), false);
6266 boolList matchedEdges(checkEdges.size(), false);
6268 // Output check entities
6269 if (debug > 4)
6270 {
6271 writeVTK
6272 (
6273 "checkEdges_" + Foam::name(eIndex),
6274 checkEdges, 1, false, true
6275 );
6277 writeVTK
6278 (
6279 "checkFaces_" + Foam::name(eIndex),
6280 checkFaces, 2, false, true
6281 );
6282 }
6284 if (localCouple && !procCouple)
6285 {
6286 // Alias for convenience...
6287 const changeMap& slaveMap = slaveMaps[0];
6289 const label pI = slaveMap.patchIndex();
6290 const coupleMap& cMap = patchCoupling_[pI].map();
6292 // Obtain references
6293 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
6294 Map<label>& rPointMap = cMap.reverseEntityMap(coupleMap::POINT);
6296 const labelList& rpList = slaveMap.removedPointList();
6297 const List<objectMap>& apList = slaveMap.addedPointList();
6299 // Configure the slave replacement point.
6300 // - collapseEdge stores this as an 'addedPoint'
6301 label scPoint = -1, srPoint = -1;
6303 if ((slaveMap.index() == eIndex) && localCouple)
6304 {
6305 // Cyclic edge at axis
6306 scPoint = collapsePoint;
6307 srPoint = replacePoint;
6308 }
6309 else
6310 {
6311 scPoint = rpList[0];
6312 srPoint = apList[0].index();
6313 }
6315 if (collapsingSlave)
6316 {
6317 if (rPointMap[replacePoint] == scPoint)
6318 {
6319 pointMap[srPoint] = replacePoint;
6320 rPointMap[replacePoint] = srPoint;
6321 }
6322 }
6323 else
6324 {
6325 if (pointMap[replacePoint] == scPoint)
6326 {
6327 rPointMap[srPoint] = replacePoint;
6328 pointMap[replacePoint] = srPoint;
6329 }
6330 }
6332 // Update face maps
6333 const label faceEnum = coupleMap::FACE;
6335 // Obtain references
6336 Map<label>& faceMap = cMap.entityMap(faceEnum);
6337 Map<label>& rFaceMap = cMap.reverseEntityMap(faceEnum);
6339 forAll(checkFaces, faceI)
6340 {
6341 label mfIndex = checkFaces[faceI];
6342 label mfPatch = whichPatch(mfIndex);
6344 const face& mF = faces_[mfIndex];
6346 triFace cF
6347 (
6348 pointMap.found(mF[0]) ? pointMap[mF[0]] : -1,
6349 pointMap.found(mF[1]) ? pointMap[mF[1]] : -1,
6350 pointMap.found(mF[2]) ? pointMap[mF[2]] : -1
6351 );
6353 if (cF[0] == -1 || cF[1] == -1 || cF[2] == -1)
6354 {
6355 writeVTK
6356 (
6357 "failedFace_"
6358 + Foam::name(mfIndex),
6359 mfIndex,
6360 2, false, true
6361 );
6363 Pout<< " Failed to configure face for: "
6364 << mfIndex << " :: " << faces_[mfIndex]
6365 << " using comparison face: " << cF
6366 << abort(FatalError);
6367 }
6369 // Fetch edges connected to first slave point
6370 const labelList& spEdges = pointEdges_[cF[0]];
6372 forAll(spEdges, edgeJ)
6373 {
6374 label seIndex = spEdges[edgeJ];
6376 if (whichEdgePatch(seIndex) == -1)
6377 {
6378 continue;
6379 }
6381 const labelList& seFaces = edgeFaces_[seIndex];
6383 forAll(seFaces, faceJ)
6384 {
6385 label sfIndex = seFaces[faceJ];
6387 if (whichPatch(sfIndex) == -1)
6388 {
6389 continue;
6390 }
6392 const face& sF = faces_[sfIndex];
6394 if
6395 (
6396 triFace::compare
6397 (
6398 triFace(sF[0], sF[1], sF[2]), cF
6399 )
6400 )
6401 {
6402 if (debug > 2)
6403 {
6404 word pN(boundaryMesh()[mfPatch].name());
6406 Pout<< " Found face: " << sfIndex
6407 << " :: " << sF
6408 << " with mfIndex: " << mfIndex
6409 << " :: " << faces_[mfIndex]
6410 << " Patch: " << pN
6411 << endl;
6412 }
6414 if (rFaceMap.found(sfIndex))
6415 {
6416 rFaceMap[sfIndex] = mfIndex;
6417 }
6418 else
6419 {
6420 rFaceMap.insert(sfIndex, mfIndex);
6421 }
6423 if (faceMap.found(mfIndex))
6424 {
6425 faceMap[mfIndex] = sfIndex;
6426 }
6427 else
6428 {
6429 faceMap.insert(mfIndex, sfIndex);
6430 }
6432 matchedFaces[faceI] = true;
6434 break;
6435 }
6436 }
6438 if (matchedFaces[faceI])
6439 {
6440 break;
6441 }
6442 }
6444 if (!matchedFaces[faceI])
6445 {
6446 writeVTK
6447 (
6448 "failedFacePoints_"
6449 + Foam::name(mfIndex),
6450 labelList(cF), 0, false, true
6451 );
6453 Pout<< " Failed to match face: "
6454 << mfIndex << " :: " << faces_[mfIndex]
6455 << " using comparison face: " << cF
6456 << abort(FatalError);
6457 }
6458 }
6460 // Update edge maps
6461 const label edgeEnum = coupleMap::EDGE;
6463 // Obtain references
6464 Map<label>& edgeMap = cMap.entityMap(edgeEnum);
6465 Map<label>& rEdgeMap = cMap.reverseEntityMap(edgeEnum);
6467 forAll(checkEdges, edgeI)
6468 {
6469 label meIndex = checkEdges[edgeI];
6471 const edge& mE = edges_[meIndex];
6473 edge cE
6474 (
6475 pointMap.found(mE[0]) ? pointMap[mE[0]] : -1,
6476 pointMap.found(mE[1]) ? pointMap[mE[1]] : -1
6477 );
6479 if (cE[0] == -1 || cE[1] == -1)
6480 {
6481 writeVTK
6482 (
6483 "failedEdge_"
6484 + Foam::name(meIndex),
6485 meIndex,
6486 1, false, true
6487 );
6489 Pout<< " Failed to configure edge for: "
6490 << meIndex << " :: " << edges_[meIndex]
6491 << " using comparison edge: " << cE
6492 << abort(FatalError);
6493 }
6495 // Fetch edges connected to first slave point
6496 const labelList& spEdges = pointEdges_[cE[0]];
6498 forAll(spEdges, edgeJ)
6499 {
6500 label seIndex = spEdges[edgeJ];
6502 const edge& sE = edges_[seIndex];
6504 if (sE == cE)
6505 {
6506 if (debug > 2)
6507 {
6508 Pout<< " Found edge: " << seIndex
6509 << " :: " << sE
6510 << " with meIndex: " << meIndex
6511 << " :: " << edges_[meIndex]
6512 << endl;
6513 }
6515 // Update reverse map
6516 if (rEdgeMap.found(seIndex))
6517 {
6518 rEdgeMap[seIndex] = meIndex;
6519 }
6520 else
6521 {
6522 rEdgeMap.insert(seIndex, meIndex);
6523 }
6525 // Update map
6526 if (edgeMap.found(meIndex))
6527 {
6528 edgeMap[meIndex] = seIndex;
6529 }
6530 else
6531 {
6532 edgeMap.insert(meIndex, seIndex);
6533 }
6535 matchedEdges[edgeI] = true;
6537 break;
6538 }
6539 }
6541 if (!matchedEdges[edgeI])
6542 {
6543 writeVTK
6544 (
6545 "failedEdgePoints_"
6546 + Foam::name(meIndex),
6547 labelList(cE), 0, false, true
6548 );
6550 Pout<< " Failed to match edge: "
6551 << meIndex << " :: " << edges_[meIndex]
6552 << " using comparison edge: " << cE
6553 << abort(FatalError);
6554 }
6555 }
6556 }
6557 else
6558 if (procCouple && !localCouple)
6559 {
6560 // Update point mapping
6561 forAll(procIndices_, pI)
6562 {
6563 const coupleMap& cMap = recvMeshes_[pI].map();
6565 // Obtain references
6566 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
6567 Map<label>& rPointMap = cMap.reverseEntityMap(coupleMap::POINT);
6569 const changeMap* slaveMapPtr = NULL;
6570 const label pointEnum = coupleMap::POINT;
6572 forAll(slaveMaps, slaveI)
6573 {
6574 const changeMap& slaveMap = slaveMaps[slaveI];
6576 if (slaveMap.patchIndex() == pI)
6577 {
6578 if (slaveMap.index() < 0)
6579 {
6580 // Point-coupling
6581 label sI = -1;
6583 if (collapsingSlave)
6584 {
6585 sI = cMap.findMaster(pointEnum, collapsePoint);
6587 if (sI > -1)
6588 {
6589 if (rPointMap.found(replacePoint))
6590 {
6591 rPointMap[replacePoint] = sI;
6592 }
6593 else
6594 {
6595 rPointMap.insert(replacePoint, sI);
6596 }
6598 pointMap[sI] = replacePoint;
6599 }
6600 }
6601 else
6602 {
6603 sI = cMap.findSlave(pointEnum, collapsePoint);
6605 if (sI > -1)
6606 {
6607 if (pointMap.found(replacePoint))
6608 {
6609 pointMap[replacePoint] = sI;
6610 }
6611 else
6612 {
6613 pointMap.insert(replacePoint, sI);
6614 }
6616 rPointMap[sI] = replacePoint;
6617 }
6618 }
6620 if (sI > -1 && debug > 2)
6621 {
6622 Pout<< " Found point: " << collapsePoint
6623 << " on proc: " << procIndices_[pI]
6624 << endl;
6625 }
6626 }
6627 else
6628 {
6629 // Edge-coupling. Fetch address for later.
6630 slaveMapPtr = &slaveMap;
6631 break;
6632 }
6633 }
6634 }
6636 if (slaveMapPtr)
6637 {
6638 // Alias for convenience...
6639 const changeMap& slaveMap = *slaveMapPtr;
6641 const labelList& rpList = slaveMap.removedPointList();
6642 const List<objectMap>& apList = slaveMap.addedPointList();
6644 // Configure the slave replacement point.
6645 // - collapseEdge stores this as an 'addedPoint'
6646 label scPoint = rpList[0];
6647 label srPoint = apList[0].index();
6649 if (collapsingSlave)
6650 {
6651 if (rPointMap[replacePoint] == scPoint)
6652 {
6653 pointMap[srPoint] = replacePoint;
6654 rPointMap[replacePoint] = srPoint;
6655 }
6657 pointMap.erase(rPointMap[collapsePoint]);
6658 rPointMap.erase(collapsePoint);
6659 }
6660 else
6661 {
6662 if (pointMap[replacePoint] == scPoint)
6663 {
6664 rPointMap[srPoint] = replacePoint;
6665 pointMap[replacePoint] = srPoint;
6666 }
6668 rPointMap.erase(pointMap[collapsePoint]);
6669 pointMap.erase(collapsePoint);
6670 }
6672 // If any other points were removed, update map
6673 for (label pointI = 1; pointI < rpList.size(); pointI++)
6674 {
6675 if (collapsingSlave)
6676 {
6677 if (pointMap.found(rpList[pointI]))
6678 {
6679 rPointMap.erase(pointMap[rpList[pointI]]);
6680 pointMap.erase(rpList[pointI]);
6681 }
6682 }
6683 else
6684 {
6685 if (rPointMap.found(rpList[pointI]))
6686 {
6687 if (debug > 2)
6688 {
6689 Pout<< " Found removed point: "
6690 << rpList[pointI]
6691 << " on proc: " << procIndices_[pI]
6692 << " for point on this proc: "
6693 << rPointMap[rpList[pointI]]
6694 << endl;
6695 }
6697 pointMap.erase(rPointMap[rpList[pointI]]);
6698 rPointMap.erase(rpList[pointI]);
6699 }
6700 }
6701 }
6702 }
6703 }
6705 // Update face mapping
6706 forAll(procIndices_, pI)
6707 {
6708 const coupleMap& cMap = recvMeshes_[pI].map();
6709 const dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
6711 // Obtain point maps
6712 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
6714 // Update face maps
6715 const label faceEnum = coupleMap::FACE;
6717 // Obtain references
6718 Map<label>& faceMap = cMap.entityMap(faceEnum);
6719 Map<label>& rFaceMap = cMap.reverseEntityMap(faceEnum);
6721 const changeMap* slaveMapPtr = NULL;
6723 forAll(slaveMaps, slaveI)
6724 {
6725 const changeMap& slaveMap = slaveMaps[slaveI];
6727 if (slaveMap.patchIndex() == pI)
6728 {
6729 if (slaveMap.index() < 0)
6730 {
6731 // Point-coupling
6732 continue;
6733 }
6735 // Edge-coupling. Fetch address for later.
6736 slaveMapPtr = &slaveMap;
6737 break;
6738 }
6739 }
6741 forAll(checkFaces, faceI)
6742 {
6743 label mfIndex = checkFaces[faceI];
6745 const face& mF = faces_[mfIndex];
6747 // Skip checks if a conversion occurred,
6748 // and this face was removed as a result
6749 if (mF.empty())
6750 {
6751 continue;
6752 }
6754 label mfPatch = whichPatch(mfIndex);
6755 label neiProc = getNeighbourProcessor(mfPatch);
6757 triFace cF
6758 (
6759 pointMap.found(mF[0]) ? pointMap[mF[0]] : -1,
6760 pointMap.found(mF[1]) ? pointMap[mF[1]] : -1,
6761 pointMap.found(mF[2]) ? pointMap[mF[2]] : -1
6762 );
6764 // Check if a patch conversion is necessary
6765 label newPatch = -1;
6766 bool requireConversion = false, physConvert = false;
6768 // Skip mapping if all points were not found
6769 if (cF[0] == -1 || cF[1] == -1 || cF[2] == -1)
6770 {
6771 // Is this actually a non-processor patch?
6772 if (neiProc == -1)
6773 {
6774 continue;
6775 }
6776 else
6777 {
6778 physConvert = true;
6779 }
6780 }
6782 // Has this face been converted to a physical boundary?
6783 if (faceMap.found(mfIndex) && slaveMapPtr)
6784 {
6785 // Alias for convenience...
6786 const changeMap& sMap = *slaveMapPtr;
6787 const labelList& rfList = sMap.removedFaceList();
6788 const List<objectMap>& afList = sMap.addedFaceList();
6790 // Was the face removed by the slave?
6791 label sfIndex = faceMap[mfIndex];
6793 if (findIndex(rfList, sfIndex) > -1)
6794 {
6795 if (debug > 2)
6796 {
6797 Pout<< " Found removed face: " << sfIndex
6798 << " with mfIndex: " << mfIndex
6799 << " :: " << faces_[mfIndex]
6800 << " on proc: " << procIndices_[pI]
6801 << endl;
6802 }
6804 // Remove map entry
6805 rFaceMap.erase(sfIndex);
6806 faceMap.erase(mfIndex);
6807 }
6809 // Check added faces for special entries
6810 forAll(afList, indexI)
6811 {
6812 const objectMap& af = afList[indexI];
6814 label mo =
6815 (
6816 af.masterObjects().size() ?
6817 af.masterObjects()[0] : 0
6818 );
6820 // Back out the physical patch ID
6821 if ((af.index() == sfIndex) && (mo < 0))
6822 {
6823 newPatch = Foam::mag(mo + 2);
6824 requireConversion = true;
6825 break;
6826 }
6827 }
6828 }
6830 // Was an appropriate physical patch found?
6831 if (physConvert && !requireConversion)
6832 {
6833 continue;
6834 }
6836 // Are we talking to a different processor?
6837 if (neiProc != procIndices_[pI])
6838 {
6839 // This face needs to be converted
6840 const polyBoundaryMesh& boundary = boundaryMesh();
6842 forAll(boundary, pi)
6843 {
6844 if (getNeighbourProcessor(pi) == procIndices_[pI])
6845 {
6846 newPatch = pi;
6847 requireConversion = true;
6848 break;
6849 }
6850 }
6851 }
6853 if (requireConversion)
6854 {
6855 // Obtain a copy before adding the new face,
6856 // since the reference might become
6857 // invalid during list resizing.
6858 face newFace = faces_[mfIndex];
6859 label newOwn = owner_[mfIndex];
6860 labelList newFaceEdges = faceEdges_[mfIndex];
6862 label newFaceIndex =
6863 (
6864 insertFace
6865 (
6866 newPatch,
6867 newFace,
6868 newOwn,
6869 -1,
6870 newFaceEdges
6871 )
6872 );
6874 meshOps::replaceLabel
6875 (
6876 mfIndex,
6877 newFaceIndex,
6878 cells_[newOwn]
6879 );
6881 // Correct edgeFaces with the new face label.
6882 forAll(newFaceEdges, edgeI)
6883 {
6884 meshOps::replaceLabel
6885 (
6886 mfIndex,
6887 newFaceIndex,
6888 edgeFaces_[newFaceEdges[edgeI]]
6889 );
6890 }
6892 // Finally remove the face
6893 removeFace(mfIndex);
6895 // Update map
6896 map.removeFace(mfIndex);
6897 map.addFace(newFaceIndex, labelList(1, mfIndex));
6899 // Replace index and patch
6900 mfIndex = newFaceIndex;
6901 mfPatch = newPatch;
6903 // If conversion was to a physical patch,
6904 // skip the remaining face mapping steps
6905 if (getNeighbourProcessor(newPatch) == -1)
6906 {
6907 continue;
6908 }
6910 // Fail for now
6911 Pout<< " Conversion required... "
6912 << " Face: " << newFace << " : "
6913 << newFace.centre(points_)
6914 << abort(FatalError);
6916 // Push a patch-conversion operation
6917 cMap.pushOperation
6918 (
6919 newFaceIndex,
6920 coupleMap::CONVERT_PATCH,
6921 newFace.centre(points_),
6922 newFace.centre(oldPoints_)
6923 );
6924 }
6926 // Fetch edges connected to first slave point
6927 const labelList& spEdges = sMesh.pointEdges_[cF[0]];
6929 forAll(spEdges, edgeJ)
6930 {
6931 label seIndex = spEdges[edgeJ];
6933 if (sMesh.whichEdgePatch(seIndex) == -1)
6934 {
6935 continue;
6936 }
6938 const labelList& seFaces = sMesh.edgeFaces_[seIndex];
6940 forAll(seFaces, faceJ)
6941 {
6942 label sfIndex = seFaces[faceJ];
6944 if (sMesh.whichPatch(sfIndex) == -1)
6945 {
6946 continue;
6947 }
6949 const face& sF = sMesh.faces_[sfIndex];
6951 if
6952 (
6953 triFace::compare
6954 (
6955 triFace(sF[0], sF[1], sF[2]), cF
6956 )
6957 )
6958 {
6959 if (debug > 2)
6960 {
6961 word pN(boundaryMesh()[mfPatch].name());
6963 Pout<< " Found face: " << sfIndex
6964 << " :: " << sF
6965 << " with mfIndex: " << mfIndex
6966 << " :: " << faces_[mfIndex]
6967 << " on proc: " << procIndices_[pI]
6968 << " Patch: " << pN
6969 << endl;
6970 }
6972 if (rFaceMap.found(sfIndex))
6973 {
6974 rFaceMap[sfIndex] = mfIndex;
6975 }
6976 else
6977 {
6978 rFaceMap.insert(sfIndex, mfIndex);
6979 }
6981 if (faceMap.found(mfIndex))
6982 {
6983 faceMap[mfIndex] = sfIndex;
6984 }
6985 else
6986 {
6987 faceMap.insert(mfIndex, sfIndex);
6988 }
6990 matchedFaces[faceI] = true;
6992 break;
6993 }
6994 }
6996 if (matchedFaces[faceI])
6997 {
6998 break;
6999 }
7000 }
7002 if ((debug > 4) && !matchedFaces[faceI])
7003 {
7004 sMesh.writeVTK
7005 (
7006 "failedFacePoints_"
7007 + Foam::name(mfIndex),
7008 labelList(cF), 0, false, true
7009 );
7011 Pout<< " Failed to match face: "
7012 << mfIndex << " :: " << faces_[mfIndex]
7013 << " using comparison face: " << cF
7014 << " on proc: " << procIndices_[pI]
7015 << endl;
7016 }
7017 }
7018 }
7020 // Update edge mapping
7021 forAll(procIndices_, pI)
7022 {
7023 const coupleMap& cMap = recvMeshes_[pI].map();
7024 const dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
7026 // Obtain point maps
7027 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
7029 // Update edge maps
7030 const label edgeEnum = coupleMap::EDGE;
7032 // Obtain references
7033 Map<label>& edgeMap = cMap.entityMap(edgeEnum);
7034 Map<label>& rEdgeMap = cMap.reverseEntityMap(edgeEnum);
7036 const changeMap* slaveMapPtr = NULL;
7038 forAll(slaveMaps, slaveI)
7039 {
7040 const changeMap& slaveMap = slaveMaps[slaveI];
7042 if (slaveMap.patchIndex() == pI)
7043 {
7044 if (slaveMap.index() < 0)
7045 {
7046 // Point-coupling
7047 continue;
7048 }
7050 // Edge-coupling. Fetch address for later.
7051 slaveMapPtr = &slaveMap;
7052 break;
7053 }
7054 }
7056 forAll(checkEdges, edgeI)
7057 {
7058 label meIndex = checkEdges[edgeI];
7060 const edge& mE = edges_[meIndex];
7062 // Skip checks if a conversion occurred,
7063 // and this edge was removed as a result
7064 if (edgeFaces_[meIndex].empty())
7065 {
7066 continue;
7067 }
7069 label mePatch = whichEdgePatch(meIndex);
7070 label neiProc = getNeighbourProcessor(mePatch);
7072 edge cE
7073 (
7074 pointMap.found(mE[0]) ? pointMap[mE[0]] : -1,
7075 pointMap.found(mE[1]) ? pointMap[mE[1]] : -1
7076 );
7078 // Check if a patch conversion is necessary
7079 label newPatch = -1;
7080 bool requireConversion = true, physConvert = false;
7082 // Skip mapping if all points were not found
7083 if (cE[0] == -1 || cE[1] == -1)
7084 {
7085 // Is this actually a non-processor patch?
7086 if (neiProc == -1)
7087 {
7088 continue;
7089 }
7090 else
7091 {
7092 physConvert = true;
7093 }
7094 }
7096 // Has this edge been converted to a physical boundary?
7097 const labelList& meFaces = edgeFaces_[meIndex];
7099 forAll(meFaces, faceI)
7100 {
7101 label mfPatch = whichPatch(meFaces[faceI]);
7103 if (mfPatch == -1)
7104 {
7105 continue;
7106 }
7108 if (getNeighbourProcessor(mfPatch) == -1)
7109 {
7110 // Store physical patch for conversion
7111 newPatch = mfPatch;
7112 }
7113 else
7114 {
7115 requireConversion = false;
7116 }
7117 }
7119 // Was an appropriate physical patch found?
7120 if (physConvert && !requireConversion)
7121 {
7122 continue;
7123 }
7125 if (requireConversion)
7126 {
7127 edge newEdge = edges_[meIndex];
7128 labelList newEdgeFaces = edgeFaces_[meIndex];
7130 // Insert the new edge
7131 label newEdgeIndex =
7132 (
7133 insertEdge
7134 (
7135 newPatch,
7136 newEdge,
7137 newEdgeFaces
7138 )
7139 );
7141 // Replace faceEdges for all
7142 // connected faces.
7143 forAll(newEdgeFaces, faceI)
7144 {
7145 meshOps::replaceLabel
7146 (
7147 meIndex,
7148 newEdgeIndex,
7149 faceEdges_[newEdgeFaces[faceI]]
7150 );
7151 }
7153 // Remove the edge
7154 removeEdge(meIndex);
7156 // Update map
7157 map.removeEdge(meIndex);
7158 map.addEdge(newEdgeIndex, labelList(1, meIndex));
7160 // Replace index and patch
7161 meIndex = newEdgeIndex;
7162 mePatch = newPatch;
7164 // If conversion was to a physical patch,
7165 // skip the remaining face mapping steps
7166 if (getNeighbourProcessor(newPatch) == -1)
7167 {
7168 continue;
7169 }
7170 }
7172 // Fetch edges connected to first slave point
7173 const labelList& spEdges = sMesh.pointEdges_[cE[0]];
7175 forAll(spEdges, edgeJ)
7176 {
7177 label seIndex = spEdges[edgeJ];
7179 const edge& sE = sMesh.edges_[seIndex];
7181 if (sE == cE)
7182 {
7183 if (debug > 2)
7184 {
7185 Pout<< " Found edge: " << seIndex
7186 << " :: " << sE
7187 << " with meIndex: " << meIndex
7188 << " :: " << edges_[meIndex]
7189 << " on proc: " << procIndices_[pI]
7190 << endl;
7191 }
7193 // Update reverse map
7194 if (rEdgeMap.found(seIndex))
7195 {
7196 rEdgeMap[seIndex] = meIndex;
7197 }
7198 else
7199 {
7200 rEdgeMap.insert(seIndex, meIndex);
7201 }
7203 // Update map
7204 if (edgeMap.found(meIndex))
7205 {
7206 edgeMap[meIndex] = seIndex;
7207 }
7208 else
7209 {
7210 edgeMap.insert(meIndex, seIndex);
7211 }
7213 matchedEdges[edgeI] = true;
7215 break;
7216 }
7217 }
7219 if (!matchedEdges[edgeI])
7220 {
7221 // Check if a coupling existed before
7222 if (edgeMap.found(meIndex) && slaveMapPtr)
7223 {
7224 // Alias for convenience...
7225 const changeMap& sMap = *slaveMapPtr;
7226 const labelList& reList = sMap.removedEdgeList();
7228 // Was the edge removed by the slave?
7229 if (findIndex(reList, edgeMap[meIndex]) > -1)
7230 {
7231 if (debug > 2)
7232 {
7233 Pout<< " Found removed edge: "
7234 << edgeMap[meIndex]
7235 << " with meIndex: " << meIndex
7236 << " :: " << edges_[meIndex]
7237 << " on proc: " << procIndices_[pI]
7238 << endl;
7239 }
7241 // Remove map entry
7242 rEdgeMap.erase(edgeMap[meIndex]);
7243 edgeMap.erase(meIndex);
7245 matchedEdges[edgeI] = true;
7246 }
7247 }
7248 }
7250 if ((debug > 4) && !matchedEdges[edgeI])
7251 {
7252 sMesh.writeVTK
7253 (
7254 "failedEdgePoints_"
7255 + Foam::name(meIndex),
7256 labelList(cE), 0, false, true
7257 );
7259 Pout<< " Failed to match edge: "
7260 << meIndex << " :: " << edges_[meIndex]
7261 << " using comparison edge: " << cE
7262 << " on proc: " << procIndices_[pI]
7263 << endl;
7264 }
7265 }
7266 }
7267 }
7269 // Ensure that all entities were matched
7270 label nFailFace = 0, nFailEdge = 0;
7272 forAll(matchedFaces, faceI)
7273 {
7274 if (!matchedFaces[faceI])
7275 {
7276 ++nFailFace;
7277 }
7278 }
7280 forAll(matchedEdges, edgeI)
7281 {
7282 if (!matchedEdges[edgeI])
7283 {
7284 ++nFailEdge;
7285 }
7286 }
7288 if (nFailFace || nFailEdge)
7289 {
7290 Pout<< " Failed to match all entities. " << nl
7291 << " Faces: " << nFailFace << nl
7292 << " Edges: " << nFailEdge << nl
7293 << abort(FatalError);
7294 }
7295 }
7297 // Set the flag
7298 topoChangeFlag_ = true;
7300 // Increment the counter
7301 status(TOTAL_COLLAPSES)++;
7303 // Increment the number of modifications
7304 status(TOTAL_MODIFICATIONS)++;
7306 // Return a succesful collapse
7307 map.type() = collapseCase;
7309 return map;
7310 }
7313 // Remove cell layer above specified patch
7314 const changeMap dynamicTopoFvMesh::removeCellLayer
7315 (
7316 const label patchID
7317 )
7318 {
7319 changeMap map;
7321 labelHashSet edgesToRemove, facesToRemove;
7322 Map<labelPair> pointsToRemove, edgesToKeep;
7324 dynamicLabelList patchFaces(patchSizes_[patchID]);
7325 DynamicList<labelPair> cellsToRemove(patchSizes_[patchID]);
7326 DynamicList<labelPair> hFacesToRemove(patchSizes_[patchID]);
7328 for (label faceI = nOldInternalFaces_; faceI < faces_.size(); faceI++)
7329 {
7330 label pIndex = whichPatch(faceI);
7332 if (pIndex != patchID)
7333 {
7334 continue;
7335 }
7337 // Fetch owner cell
7338 label cIndex = owner_[faceI];
7340 // Add face to the list
7341 patchFaces.append(faceI);
7343 // Fetch appropriate face / cell
7344 const face& bFace = faces_[faceI];
7345 const cell& ownCell = cells_[cIndex];
7347 // Get the opposing face from the cell
7348 oppositeFace oFace = ownCell.opposingFace(faceI, faces_);
7350 if (!oFace.found())
7351 {
7352 // Something's wrong here.
7353 FatalErrorIn
7354 (
7355 "const changeMap dynamicTopoFvMesh::removeCellLayer"
7356 "(const label patchID)"
7357 )
7358 << " Face: " << faceI << " :: " << bFace << nl
7359 << " has no opposing face in cell: "
7360 << cIndex << " :: " << ownCell << nl
7361 << abort(FatalError);
7362 }
7364 // Fetch cell on the other-side of the opposite face
7365 label otherCellIndex =
7366 (
7367 (owner_[oFace.oppositeIndex()] == cIndex) ?
7368 neighbour_[oFace.oppositeIndex()] :
7369 owner_[oFace.oppositeIndex()]
7370 );
7372 if (otherCellIndex == -1)
7373 {
7374 // Opposite face is on a boundary, and layer
7375 // removal would be invalid if we continued.
7376 map.type() = -2;
7378 return map;
7379 }
7381 // Fetch reference to other cell
7382 const cell& otherCell = cells_[otherCellIndex];
7384 // Find opposite face on the other cell
7385 oppositeFace otFace =
7386 (
7387 otherCell.opposingFace
7388 (
7389 oFace.oppositeIndex(),
7390 faces_
7391 )
7392 );
7394 if (!otFace.found())
7395 {
7396 // Something's wrong here.
7397 FatalErrorIn
7398 (
7399 "const changeMap dynamicTopoFvMesh::removeCellLayer"
7400 "(const label patchID)"
7401 )
7402 << " Face: " << oFace.oppositeIndex()
7403 << " :: " << oFace << nl
7404 << " has no opposing face in cell: "
7405 << otherCellIndex << " :: " << otherCell << nl
7406 << abort(FatalError);
7407 }
7409 // All edges on the boundary face are to be retained
7410 const labelList& fEdges = faceEdges_[faceI];
7411 const labelList& ofEdges = faceEdges_[oFace.oppositeIndex()];
7412 const labelList& otfEdges = faceEdges_[otFace.oppositeIndex()];
7414 forAll(fEdges, edgeI)
7415 {
7416 label eIndex = fEdges[edgeI];
7418 if (!edgesToKeep.found(eIndex))
7419 {
7420 // Find equivalent edge on opposite face
7421 label oeIndex = -1, oteIndex = -1;
7422 const edge& bEdge = edges_[eIndex];
7424 // Build edges for comparison
7425 label startLoc = bFace.which(bEdge[0]);
7426 label endLoc = bFace.which(bEdge[1]);
7428 edge cEdge(oFace[startLoc], oFace[endLoc]);
7429 edge ctEdge(otFace[startLoc], otFace[endLoc]);
7431 forAll(ofEdges, edgeJ)
7432 {
7433 const edge& ofEdge = edges_[ofEdges[edgeJ]];
7435 if (cEdge == ofEdge)
7436 {
7437 oeIndex = ofEdges[edgeJ];
7438 break;
7439 }
7440 }
7442 forAll(otfEdges, edgeJ)
7443 {
7444 const edge& otfEdge = edges_[otfEdges[edgeJ]];
7446 if (ctEdge == otfEdge)
7447 {
7448 oteIndex = otfEdges[edgeJ];
7449 break;
7450 }
7451 }
7453 if (oeIndex < 0 || oteIndex < 0)
7454 {
7455 FatalErrorIn
7456 (
7457 "const changeMap dynamicTopoFvMesh::removeCellLayer"
7458 "(const label patchID)"
7459 )
7460 << " Could not find comparison edge: " << nl
7461 << " cEdge: " << cEdge
7462 << " oeIndex: " << oeIndex << nl
7463 << " ctEdge: " << ctEdge
7464 << " oteIndex: " << oteIndex << nl
7465 << " for edge: " << bEdge
7466 << abort(FatalError);
7467 }
7469 // Make entry
7470 edgesToKeep.insert(eIndex, labelPair(oeIndex, oteIndex));
7471 }
7472 }
7474 // Add information to removal lists
7475 cellsToRemove.append
7476 (
7477 labelPair
7478 (
7479 cIndex,
7480 otherCellIndex
7481 )
7482 );
7484 hFacesToRemove.append
7485 (
7486 labelPair
7487 (
7488 oFace.oppositeIndex(),
7489 otFace.oppositeIndex()
7490 )
7491 );
7493 // Mark points for removal
7494 forAll(oFace, pointI)
7495 {
7496 label pIndex = oFace[pointI];
7498 if (!pointsToRemove.found(pIndex))
7499 {
7500 // Make entry
7501 pointsToRemove.insert
7502 (
7503 pIndex,
7504 labelPair(bFace[pointI], otFace[pointI])
7505 );
7506 }
7507 }
7509 // Loop through cell faces and mark
7510 // faces / edges for removal
7511 forAll(ownCell, faceJ)
7512 {
7513 label fIndex = ownCell[faceJ];
7515 if (fIndex == faceI || fIndex == oFace.oppositeIndex())
7516 {
7517 continue;
7518 }
7520 if (!facesToRemove.found(fIndex))
7521 {
7522 facesToRemove.insert(fIndex);
7523 }
7525 const labelList& checkEdges = faceEdges_[fIndex];
7527 forAll(checkEdges, edgeI)
7528 {
7529 label eIndex = checkEdges[edgeI];
7531 if (edgesToKeep.found(eIndex))
7532 {
7533 continue;
7534 }
7536 if (!edgesToRemove.found(eIndex))
7537 {
7538 edgesToRemove.insert(eIndex);
7539 }
7540 }
7541 }
7542 }
7544 // Correct edgeFaces / faceEdges for retained edges
7545 forAllConstIter(Map<labelPair>, edgesToKeep, eIter)
7546 {
7547 const labelList& rmeFaces = edgeFaces_[eIter().first()];
7549 forAll(rmeFaces, faceI)
7550 {
7551 labelList& fE = faceEdges_[rmeFaces[faceI]];
7553 bool foundRp = (findIndex(fE, eIter.key()) > -1);
7554 bool foundRn = (findIndex(fE, eIter().second()) > -1);
7556 if (foundRp)
7557 {
7558 // Size-down edgeFaces for replacement
7559 meshOps::sizeDownList
7560 (
7561 rmeFaces[faceI],
7562 edgeFaces_[eIter.key()]
7563 );
7564 }
7566 if (foundRn)
7567 {
7568 // Size-up edgeFaces for replacement
7569 meshOps::sizeUpList
7570 (
7571 rmeFaces[faceI],
7572 edgeFaces_[eIter.key()]
7573 );
7575 // Replace edge index
7576 meshOps::replaceLabel
7577 (
7578 eIter().first(),
7579 eIter.key(),
7580 fE
7581 );
7582 }
7583 }
7584 }
7586 // Remove unwanted faces
7587 forAllConstIter(labelHashSet, facesToRemove, fIter)
7588 {
7589 // Remove the face
7590 removeFace(fIter.key());
7592 // Update map
7593 map.removeFace(fIter.key());
7594 }
7596 // Remove unwanted edges
7597 forAllConstIter(labelHashSet, edgesToRemove, eIter)
7598 {
7599 // Remove the edge
7600 removeEdge(eIter.key());
7602 // Update map
7603 map.removeEdge(eIter.key());
7604 }
7606 // Remove unwanted points
7607 forAllConstIter(Map<labelPair>, pointsToRemove, pIter)
7608 {
7609 // Update pointEdges information first
7610 if (is3D())
7611 {
7612 const labelList& pEdges = pointEdges_[pIter.key()];
7614 // Configure edge for comparison
7615 edge cEdge
7616 (
7617 pIter.key(),
7618 pIter().second()
7619 );
7621 label replaceEdge = -1;
7623 forAll(pEdges, edgeI)
7624 {
7625 const edge& checkEdge = edges_[pEdges[edgeI]];
7627 if (checkEdge == cEdge)
7628 {
7629 replaceEdge = pEdges[edgeI];
7630 break;
7631 }
7632 }
7634 if (replaceEdge == -1)
7635 {
7636 FatalErrorIn
7637 (
7638 "const changeMap dynamicTopoFvMesh::removeCellLayer"
7639 "(const label patchID)"
7640 )
7641 << " Could not find comparison edge: " << nl
7642 << " cEdge: " << cEdge
7643 << " for point: " << pIter.key() << nl
7644 << " pointEdges: " << pEdges
7645 << abort(FatalError);
7646 }
7648 // Size-up pointEdges
7649 meshOps::sizeUpList
7650 (
7651 replaceEdge,
7652 pointEdges_[pIter().first()]
7653 );
7654 }
7656 // Remove the point
7657 removePoint(pIter.key());
7659 // Update map
7660 map.removePoint(pIter.key());
7661 }
7663 // Track modified faces for mapping
7664 labelHashSet modifiedFaces;
7666 // Remove all cells
7667 forAll(cellsToRemove, indexI)
7668 {
7669 // Replace face label on the other cell
7670 meshOps::replaceLabel
7671 (
7672 hFacesToRemove[indexI].first(),
7673 patchFaces[indexI],
7674 cells_[cellsToRemove[indexI].second()]
7675 );
7677 // Set owner information
7678 owner_[patchFaces[indexI]] = cellsToRemove[indexI].second();
7680 // Replace points on faces / edges
7681 const cell& otherCell = cells_[cellsToRemove[indexI].second()];
7683 forAll(otherCell, faceI)
7684 {
7685 face& faceToCheck = faces_[otherCell[faceI]];
7687 bool modified = false;
7689 forAll(faceToCheck, pointI)
7690 {
7691 if (pointsToRemove.found(faceToCheck[pointI]))
7692 {
7693 faceToCheck[pointI] =
7694 (
7695 pointsToRemove[faceToCheck[pointI]].first()
7696 );
7697 }
7699 modified = true;
7700 }
7702 // Add face to list
7703 if (modified && !modifiedFaces.found(otherCell[faceI]))
7704 {
7705 modifiedFaces.insert(otherCell[faceI]);
7706 }
7708 const labelList& fEdges = faceEdges_[otherCell[faceI]];
7710 forAll(fEdges, edgeI)
7711 {
7712 edge& edgeToCheck = edges_[fEdges[edgeI]];
7714 if (pointsToRemove.found(edgeToCheck[0]))
7715 {
7716 edgeToCheck[0] =
7717 (
7718 pointsToRemove[edgeToCheck[0]].first()
7719 );
7720 }
7722 if (pointsToRemove.found(edgeToCheck[1]))
7723 {
7724 edgeToCheck[1] =
7725 (
7726 pointsToRemove[edgeToCheck[1]].first()
7727 );
7728 }
7729 }
7730 }
7732 // Remove horizontal interior face
7733 removeFace(hFacesToRemove[indexI].first());
7735 // Update map
7736 map.removeFace(hFacesToRemove[indexI].first());
7738 // Remove the cell
7739 removeCell(cellsToRemove[indexI].first());
7741 // Update map
7742 map.removeCell(cellsToRemove[indexI].first());
7743 }
7745 // Now that all old / new cells possess correct connectivity,
7746 // compute mapping information.
7747 forAll(cellsToRemove, indexI)
7748 {
7749 // Set mapping for the modified cell
7750 setCellMapping
7751 (
7752 cellsToRemove[indexI].second(),
7753 labelList(cellsToRemove[indexI])
7754 );
7755 }
7757 // Set face mapping information for modified faces
7758 forAllConstIter(labelHashSet, modifiedFaces, fIter)
7759 {
7760 // Decide between default / weighted mapping
7761 // based on boundary information
7762 label fPatch = whichPatch(fIter.key());
7764 if (fPatch == -1)
7765 {
7766 // Default mapping for interior faces
7767 setFaceMapping(fIter.key());
7768 }
7769 else
7770 {
7771 // Map boundary face from itself
7772 setFaceMapping(fIter.key(), labelList(1, fIter.key()));
7773 }
7774 }
7776 // Set the flag
7777 topoChangeFlag_ = true;
7779 // Specify that the operation was successful
7780 map.type() = 1;
7782 // Return the changeMap
7783 return map;
7784 }
7787 // Merge a set of boundary faces into internal
7788 const changeMap dynamicTopoFvMesh::mergeBoundaryFaces
7789 (
7790 const labelList& mergeFaces
7791 )
7792 {
7793 changeMap map;
7795 if (debug > 2)
7796 {
7797 Pout << "Merging faces: " << mergeFaces << endl;
7798 }
7800 List<labelPair> mergePairs(0);
7802 forAll(mergeFaces, faceI)
7803 {
7804 const face& fI = faces_[mergeFaces[faceI]];
7806 for (label faceJ = faceI + 1; faceJ < mergeFaces.size(); faceJ++)
7807 {
7808 const face& fJ = faces_[mergeFaces[faceJ]];
7810 bool merge = false;
7812 if (fI.size() == fJ.size() && fI.size() == 3)
7813 {
7814 if
7815 (
7816 triFace::compare
7817 (
7818 triFace(fI[0], fI[1], fI[2]),
7819 triFace(fJ[0], fJ[1], fJ[2])
7820 )
7821 )
7822 {
7823 merge = true;
7824 }
7825 }
7826 else
7827 if (face::compare(fI, fJ))
7828 {
7829 merge = true;
7830 }
7832 if (merge)
7833 {
7834 meshOps::sizeUpList
7835 (
7836 labelPair(mergeFaces[faceI], mergeFaces[faceJ]),
7837 mergePairs
7838 );
7840 break;
7841 }
7842 }
7843 }
7845 if (debug > 2)
7846 {
7847 Pout<< " mergePairs: " << mergePairs << endl;
7848 }
7850 dynamicLabelList checkEdges(10);
7852 forAll(mergePairs, pairI)
7853 {
7854 label firstFace = mergePairs[pairI].first();
7855 label secondFace = mergePairs[pairI].second();
7857 // Obtain owners for both faces, and compare their labels
7858 label newOwner = -1, newNeighbour = -1;
7859 label removedFace = -1, retainedFace = -1;
7861 if (owner_[firstFace] < owner_[secondFace])
7862 {
7863 // Retain the first face
7864 removedFace = secondFace;
7865 retainedFace = firstFace;
7867 newOwner = owner_[firstFace];
7868 newNeighbour = owner_[secondFace];
7869 }
7870 else
7871 {
7872 // Retain the second face
7873 removedFace = firstFace;
7874 retainedFace = secondFace;
7876 newOwner = owner_[secondFace];
7877 newNeighbour = owner_[firstFace];
7878 }
7880 // Insert a new interior face
7881 label newFaceIndex =
7882 (
7883 insertFace
7884 (
7885 -1,
7886 face(faces_[retainedFace]),
7887 newOwner,
7888 newNeighbour,
7889 labelList(faceEdges_[retainedFace])
7890 )
7891 );
7893 // Update map
7894 map.addFace(newFaceIndex);
7896 // Replace cell with the new face label
7897 meshOps::replaceLabel
7898 (
7899 removedFace,
7900 newFaceIndex,
7901 cells_[owner_[removedFace]]
7902 );
7904 meshOps::replaceLabel
7905 (
7906 retainedFace,
7907 newFaceIndex,
7908 cells_[owner_[retainedFace]]
7909 );
7911 const labelList& fEdges = faceEdges_[newFaceIndex];
7912 const labelList& rfEdges = faceEdges_[removedFace];
7914 // Check for common edges on the removed face
7915 forAll(rfEdges, edgeI)
7916 {
7917 label reIndex = rfEdges[edgeI];
7919 if (findIndex(fEdges, reIndex) == -1)
7920 {
7921 // Find the equivalent edge
7922 const edge& rEdge = edges_[reIndex];
7923 const labelList& reFaces = edgeFaces_[reIndex];
7925 label keIndex = -1;
7927 forAll(fEdges, edgeJ)
7928 {
7929 if (edges_[fEdges[edgeJ]] == rEdge)
7930 {
7931 keIndex = fEdges[edgeJ];
7932 break;
7933 }
7934 }
7936 if (keIndex == -1)
7937 {
7938 Pout<< " Could not find edge for: "
7939 << reIndex << " :: " << rEdge
7940 << abort(FatalError);
7941 }
7943 // Add edgeFaces from this edge
7944 forAll(reFaces, faceI)
7945 {
7946 if (reFaces[faceI] == removedFace)
7947 {
7948 continue;
7949 }
7951 meshOps::sizeUpList
7952 (
7953 reFaces[faceI],
7954 edgeFaces_[keIndex]
7955 );
7957 meshOps::replaceLabel
7958 (
7959 reIndex,
7960 keIndex,
7961 faceEdges_[reFaces[faceI]]
7962 );
7963 }
7965 // Remove the old edge
7966 removeEdge(reIndex);
7968 // Update map
7969 map.removeEdge(reIndex);
7970 }
7971 }
7973 // Replace edgeFaces with the new face label
7974 forAll(fEdges, edgeI)
7975 {
7976 label eIndex = fEdges[edgeI];
7978 if (findIndex(edgeFaces_[eIndex], removedFace) > -1)
7979 {
7980 meshOps::sizeDownList
7981 (
7982 removedFace,
7983 edgeFaces_[eIndex]
7984 );
7985 }
7987 if (findIndex(edgeFaces_[eIndex], retainedFace) > -1)
7988 {
7989 meshOps::sizeDownList
7990 (
7991 retainedFace,
7992 edgeFaces_[eIndex]
7993 );
7994 }
7996 // Size-up list with the new face index
7997 meshOps::sizeUpList
7998 (
7999 newFaceIndex,
8000 edgeFaces_[eIndex]
8001 );
8003 if (findIndex(checkEdges, eIndex) == -1)
8004 {
8005 checkEdges.append(eIndex);
8006 }
8007 }
8009 // Remove the boundary faces
8010 removeFace(removedFace);
8011 removeFace(retainedFace);
8013 // Update map
8014 map.removeFace(removedFace);
8015 map.removeFace(retainedFace);
8016 }
8018 forAll(checkEdges, edgeI)
8019 {
8020 bool allInterior = true;
8021 label eIndex = checkEdges[edgeI];
8023 const labelList& eFaces = edgeFaces_[eIndex];
8025 forAll(eFaces, faceI)
8026 {
8027 if (whichPatch(eFaces[faceI]) > -1)
8028 {
8029 allInterior = false;
8030 break;
8031 }
8032 }
8034 if (allInterior)
8035 {
8036 // This edge needs to be converted to an interior one
8037 label newEdgeIndex =
8038 (
8039 insertEdge
8040 (
8041 -1,
8042 edge(edges_[eIndex]),
8043 labelList(edgeFaces_[eIndex])
8044 )
8045 );
8047 // Update map
8048 map.addEdge(newEdgeIndex, labelList(1, eIndex));
8050 // Update faceEdges information for all connected faces
8051 const labelList& neFaces = edgeFaces_[newEdgeIndex];
8053 forAll(neFaces, faceI)
8054 {
8055 meshOps::replaceLabel
8056 (
8057 eIndex,
8058 newEdgeIndex,
8059 faceEdges_[neFaces[faceI]]
8060 );
8061 }
8063 // Remove the old boundary edge
8064 removeEdge(eIndex);
8066 // Update map
8067 map.removeEdge(eIndex);
8069 // Replace index
8070 checkEdges[edgeI] = newEdgeIndex;
8071 }
8072 }
8074 if (debug > 2)
8075 {
8076 Pout<< "Merge complete." << nl << endl;
8077 }
8079 // Return a succesful merge
8080 map.type() = 1;
8082 return map;
8083 }
8086 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
8088 } // End namespace Foam
8090 // ************************************************************************* //