| blob | fd5e67e189602e98b3bd94a2532f19b40b39e87c |
1 /*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | Copyright held by original author
6 \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19 for more details.
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "Stack.H"
28 #include "objectRegistry.H"
29 #include "triFace.H"
30 #include "objectMap.H"
31 #include "changeMap.H"
32 #include "multiThreader.H"
33 #include "coupledInfo.H"
34 #include "dynamicTopoFvMesh.H"
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 namespace Foam
39 {
41 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
43 // Method to collapse a quad-face in 2D
44 // - Returns a changeMap with a type specifying:
45 // -3: Collapse failed since face was on a noRefinement patch.
46 // -1: Collapse failed since max number of topo-changes was reached.
47 // 0: Collapse could not be performed.
48 // 1: Collapsed to first node.
49 // 2: Collapsed to second node.
50 // 3: Collapse to mid-point.
51 // - overRideCase is used to force a certain collapse configuration.
52 // -1: Use this value to let collapseQuadFace decide a case.
53 // 1: Force collapse to first node.
54 // 2: Force collapse to second node.
55 // 3: Force collapse to mid-point.
56 // - checkOnly performs a feasibility check and returns without modifications.
57 const changeMap dynamicTopoFvMesh::collapseQuadFace
58 (
59 const label fIndex,
60 label overRideCase,
61 bool checkOnly,
62 bool forceOp
63 )
64 {
65 // Figure out which thread this is...
66 label tIndex = self();
68 // Prepare the changeMaps
69 changeMap map;
70 List<changeMap> slaveMaps;
71 bool collapsingSlave = false;
73 if
74 (
75 (statistics_[0] > maxModifications_)
76 && (maxModifications_ > -1)
77 )
78 {
79 // Reached the max allowable topo-changes.
80 stack(tIndex).clear();
82 return map;
83 }
85 // Check if edgeRefinements are to be avoided on patch.
86 if (baseMesh_.lengthEstimator().checkRefinementPatch(whichPatch(fIndex)))
87 {
88 map.type() = -3;
90 return map;
91 }
93 // Sanity check: Is the index legitimate?
94 if (fIndex < 0)
95 {
96 FatalErrorIn
97 (
98 "\n"
99 "const changeMap "
100 "dynamicTopoFvMesh::collapseQuadFace\n"
101 "(\n"
102 " const label fIndex,\n"
103 " label overRideCase,\n"
104 " bool checkOnly,\n"
105 " bool forceOp\n"
106 ")\n"
107 )
108 << " Invalid index: " << fIndex
109 << abort(FatalError);
110 }
112 // Define the edges on the face to be collapsed
113 FixedList<edge,4> checkEdge(edge(-1, -1));
114 FixedList<label,4> checkEdgeIndex(-1);
116 // Define checkEdges
117 getCheckEdges(fIndex, (*this), map, checkEdge, checkEdgeIndex);
119 // Determine the common vertices for the first and second edges
120 label cv0 = checkEdge[1].commonVertex(checkEdge[0]);
121 label cv1 = checkEdge[1].commonVertex(checkEdge[3]);
122 label cv2 = checkEdge[2].commonVertex(checkEdge[0]);
123 label cv3 = checkEdge[2].commonVertex(checkEdge[3]);
125 // If coupled modification is set, and this is a
126 // master face, collapse its slaves first.
127 bool localCouple = false, procCouple = false;
129 if (coupledModification_)
130 {
131 const face& fCheck = faces_[fIndex];
133 const label faceEnum = coupleMap::FACE;
134 const label pointEnum = coupleMap::POINT;
136 // Is this a locally coupled edge (either master or slave)?
137 if (locallyCoupledEntity(fIndex, true))
138 {
139 localCouple = true;
140 procCouple = false;
141 }
142 else
143 if (processorCoupledEntity(fIndex))
144 {
145 procCouple = true;
146 localCouple = false;
147 }
149 if (localCouple && !procCouple)
150 {
151 // Determine the slave index.
152 label sIndex = -1, pIndex = -1;
154 forAll(patchCoupling_, patchI)
155 {
156 if (patchCoupling_(patchI))
157 {
158 const coupleMap& cMap = patchCoupling_[patchI].map();
160 if ((sIndex = cMap.findSlave(faceEnum, fIndex)) > -1)
161 {
162 pIndex = patchI;
164 break;
165 }
167 // The following bit happens only during the sliver
168 // exudation process, since slave edges are
169 // usually not added to the coupled edge-stack.
170 if ((sIndex = cMap.findMaster(faceEnum, fIndex)) > -1)
171 {
172 pIndex = patchI;
174 // Notice that we are collapsing a slave edge first.
175 collapsingSlave = true;
177 break;
178 }
179 }
180 }
182 if (sIndex == -1)
183 {
184 FatalErrorIn
185 (
186 "\n"
187 "const changeMap "
188 "dynamicTopoFvMesh::collapseQuadFace\n"
189 "(\n"
190 " const label fIndex,\n"
191 " label overRideCase,\n"
192 " bool checkOnly,\n"
193 " bool forceOp\n"
194 ")\n"
195 )
196 << "Coupled maps were improperly specified." << nl
197 << " Slave index not found for: " << nl
198 << " Face: " << fIndex << ": " << fCheck << nl
199 << abort(FatalError);
200 }
201 else
202 {
203 // If we've found the slave, size up the list
204 meshOps::sizeUpList
205 (
206 changeMap(),
207 slaveMaps
208 );
210 // Save index and patch for posterity
211 slaveMaps[0].index() = sIndex;
212 slaveMaps[0].patchIndex() = pIndex;
213 }
215 if (debug > 1)
216 {
217 Pout<< nl << " >> Collapsing slave face: " << sIndex
218 << " for master face: " << fIndex << endl;
219 }
220 }
221 else
222 if (procCouple && !localCouple)
223 {
224 // If this is a new entity, bail out for now.
225 // This will be handled at the next time-step.
226 if (fIndex >= nOldFaces_)
227 {
228 return map;
229 }
231 // Check slaves
232 forAll(procIndices_, pI)
233 {
234 // Fetch reference to subMeshes
235 const coupledInfo& sendMesh = sendMeshes_[pI];
236 const coupledInfo& recvMesh = recvMeshes_[pI];
238 const coupleMap& scMap = sendMesh.map();
239 const coupleMap& rcMap = recvMesh.map();
241 // If this face was sent to a lower-ranked
242 // processor, skip it.
243 if (procIndices_[pI] < Pstream::myProcNo())
244 {
245 if (scMap.reverseEntityMap(faceEnum).found(fIndex))
246 {
247 if (debug > 3)
248 {
249 Pout<< "Face: " << fIndex
250 << "::" << fCheck
251 << " was sent to proc: "
252 << procIndices_[pI]
253 << ", so bailing out."
254 << endl;
255 }
257 return map;
258 }
259 }
261 label sIndex = -1;
263 if ((sIndex = rcMap.findSlave(faceEnum, fIndex)) > -1)
264 {
265 // Check if a lower-ranked processor is
266 // handling this face
267 if (procIndices_[pI] < Pstream::myProcNo())
268 {
269 if (debug > 3)
270 {
271 Pout<< "Face: " << fIndex
272 << "::" << fCheck
273 << " is handled by proc: "
274 << procIndices_[pI]
275 << ", so bailing out."
276 << endl;
277 }
279 return map;
280 }
282 label curIndex = slaveMaps.size();
284 // Size up the list
285 meshOps::sizeUpList
286 (
287 changeMap(),
288 slaveMaps
289 );
291 // Save index and patch for posterity
292 slaveMaps[curIndex].index() = sIndex;
293 slaveMaps[curIndex].patchIndex() = pI;
294 }
295 else
296 if
297 (
298 (
299 rcMap.findSlave(pointEnum, cv0) > -1 &&
300 rcMap.findSlave(pointEnum, cv1) > -1
301 )
302 || (
303 rcMap.findSlave(pointEnum, cv2) > -1 &&
304 rcMap.findSlave(pointEnum, cv3) > -1
305 )
306 )
307 {
308 // An edge-only coupling exists.
310 // Check if a lower-ranked processor is
311 // handling this face
312 if (procIndices_[pI] < Pstream::myProcNo())
313 {
314 if (debug > 3)
315 {
316 Pout<< "Face edge on: " << fIndex
317 << "::" << fCheck
318 << " is handled by proc: "
319 << procIndices_[pI]
320 << ", so bailing out."
321 << endl;
322 }
324 return map;
325 }
327 label p0 = rcMap.findSlave(pointEnum, cv0);
328 label p1 = rcMap.findSlave(pointEnum, cv1);
330 label p2 = rcMap.findSlave(pointEnum, cv2);
331 label p3 = rcMap.findSlave(pointEnum, cv3);
333 edge cEdge(-1, -1);
334 label edgeCouple = -1;
336 if ((p0 > -1 && p1 > -1) && (p2 == -1 && p3 == -1))
337 {
338 cEdge[0] = p0;
339 cEdge[1] = p1;
341 edgeCouple = 1;
342 sIndex = readLabel(map.lookup("firstEdge"));
343 }
344 else
345 if ((p0 == -1 && p1 == -1) && (p2 > -1 && p3 > -1))
346 {
347 cEdge[0] = p2;
348 cEdge[1] = p3;
350 edgeCouple = 2;
351 sIndex = readLabel(map.lookup("secondEdge"));
352 }
354 label curIndex = slaveMaps.size();
356 // Size up the list
357 meshOps::sizeUpList
358 (
359 changeMap(),
360 slaveMaps
361 );
363 // Unfortunately, since no edge maps are
364 // available in 2D, loop through boundary
365 // faces on subMesh and get the right edge
366 label seIndex = -1;
368 const dynamicTopoFvMesh& sMesh = recvMesh.subMesh();
370 for
371 (
372 label faceI = sMesh.nOldInternalFaces_;
373 faceI < sMesh.faceEdges_.size();
374 faceI++
375 )
376 {
377 const labelList& fEdges = sMesh.faceEdges_[faceI];
379 forAll(fEdges, edgeI)
380 {
381 if (sMesh.edges_[fEdges[edgeI]] == cEdge)
382 {
383 seIndex = fEdges[edgeI];
384 break;
385 }
386 }
388 if (seIndex > -1)
389 {
390 break;
391 }
392 }
394 if (seIndex == -1)
395 {
396 Pout<< "Face edge on: " << fIndex
397 << "::" << fCheck
398 << " sIndex: " << sIndex
399 << " could not be found."
400 << abort(FatalError);
401 }
403 // Save index and patch for posterity
404 // - Negate the index to signify edge coupling
405 slaveMaps[curIndex].index() = -seIndex;
406 slaveMaps[curIndex].patchIndex() = pI;
408 // Save edgeCouple as well, so that
409 // another map comparison is avoided.
410 slaveMaps[curIndex].type() = edgeCouple;
411 }
412 }
413 }
414 else
415 {
416 // Something's wrong with coupling maps
417 FatalErrorIn
418 (
419 "\n"
420 "const changeMap "
421 "dynamicTopoFvMesh::collapseQuadFace\n"
422 "(\n"
423 " const label fIndex,\n"
424 " label overRideCase,\n"
425 " bool checkOnly,\n"
426 " bool forceOp\n"
427 ")\n"
428 )
429 << "Coupled maps were improperly specified." << nl
430 << " localCouple: " << localCouple << nl
431 << " procCouple: " << procCouple << nl
432 << " Face: " << fIndex << ": " << fCheck << nl
433 << abort(FatalError);
434 }
436 // Temporarily turn off coupledModification
437 unsetCoupledModification();
439 // Test the master face for collapse, and decide on a case
440 changeMap masterMap = collapseQuadFace(fIndex, -1, true, forceOp);
442 // Turn it back on.
443 setCoupledModification();
445 // Master couldn't perform collapse
446 if (masterMap.type() <= 0)
447 {
448 return masterMap;
449 }
451 // For edge-only coupling, define the points for checking
452 List<FixedList<point, 2> > slaveMoveNewPoint(slaveMaps.size());
453 List<FixedList<point, 2> > slaveMoveOldPoint(slaveMaps.size());
455 // Now check each of the slaves for collapse feasibility
456 forAll(slaveMaps, slaveI)
457 {
458 // Alias for convenience...
459 changeMap& slaveMap = slaveMaps[slaveI];
461 label slaveOverRide = -1;
462 label sIndex = slaveMap.index();
463 label pI = slaveMap.patchIndex();
464 const coupleMap* cMapPtr = NULL;
466 if (localCouple)
467 {
468 cMapPtr = &(patchCoupling_[pI].map());
469 }
470 else
471 if (procCouple)
472 {
473 const dynamicTopoFvMesh& sMesh =
474 (
475 recvMeshes_[pI].subMesh()
476 );
478 cMapPtr = &(recvMeshes_[pI].map());
480 if (debug > 3)
481 {
482 if (sIndex < 0)
483 {
484 Pout<< "Checking slave edge: " << mag(sIndex)
485 << "::" << sMesh.edges_[mag(sIndex)]
486 << " on proc: " << procIndices_[pI]
487 << " for master face: " << fIndex
488 << " using collapseCase: " << masterMap.type()
489 << endl;
490 }
491 else
492 {
493 Pout<< "Checking slave face: " << sIndex
494 << "::" << sMesh.faces_[sIndex]
495 << " on proc: " << procIndices_[pI]
496 << " for master face: " << fIndex
497 << " using collapseCase: " << masterMap.type()
498 << endl;
499 }
500 }
501 }
503 // Define an overRide case for the slave
504 FixedList<edge, 2> mEdge(edge(-1, -1)), sEdge(edge(-1, -1));
506 if (collapsingSlave)
507 {
508 const Map<label>& rPointMap =
509 (
510 cMapPtr->reverseEntityMap(pointEnum)
511 );
513 if (sIndex < 0)
514 {
515 if (slaveMap.type() == 1)
516 {
517 mEdge[0][0] = rPointMap[cv0];
518 mEdge[0][1] = rPointMap[cv1];
519 }
520 else
521 if (slaveMap.type() == 2)
522 {
523 mEdge[0][0] = rPointMap[cv2];
524 mEdge[0][1] = rPointMap[cv3];
525 }
526 }
527 else
528 {
529 mEdge[0][0] = rPointMap[cv0];
530 mEdge[0][1] = rPointMap[cv1];
532 mEdge[1][0] = rPointMap[cv2];
533 mEdge[1][1] = rPointMap[cv3];
534 }
535 }
536 else
537 {
538 const Map<label>& pointMap =
539 (
540 cMapPtr->entityMap(pointEnum)
541 );
543 if (sIndex < 0)
544 {
545 if (slaveMap.type() == 1)
546 {
547 mEdge[0][0] = pointMap[cv0];
548 mEdge[0][1] = pointMap[cv1];
549 }
550 else
551 if (slaveMap.type() == 2)
552 {
553 mEdge[0][0] = pointMap[cv2];
554 mEdge[0][1] = pointMap[cv3];
555 }
556 }
557 else
558 {
559 mEdge[0][0] = pointMap[cv0];
560 mEdge[0][1] = pointMap[cv1];
562 mEdge[1][0] = pointMap[cv2];
563 mEdge[1][1] = pointMap[cv3];
564 }
565 }
567 // Determine checkEdges for the slave
568 FixedList<edge,4> slaveCheckEdge(edge(-1, -1));
569 FixedList<label,4> slaveCheckEdgeIndex(-1);
571 if (localCouple)
572 {
573 getCheckEdges
574 (
575 sIndex,
576 (*this),
577 slaveMap,
578 slaveCheckEdge,
579 slaveCheckEdgeIndex
580 );
582 sEdge[0] = edges_[readLabel(slaveMap.lookup("firstEdge"))];
583 sEdge[1] = edges_[readLabel(slaveMap.lookup("secondEdge"))];
584 }
585 else
586 if (procCouple)
587 {
588 const dynamicTopoFvMesh& sMesh =
589 (
590 recvMeshes_[pI].subMesh()
591 );
593 if (sIndex < 0)
594 {
595 sEdge[0] = sMesh.edges_[mag(sIndex)];
596 }
597 else
598 {
599 getCheckEdges
600 (
601 sIndex,
602 sMesh,
603 slaveMap,
604 slaveCheckEdge,
605 slaveCheckEdgeIndex
606 );
608 sEdge[0] =
609 (
610 sMesh.edges_[readLabel(slaveMap.lookup("firstEdge"))]
611 );
613 sEdge[1] =
614 (
615 sMesh.edges_[readLabel(slaveMap.lookup("secondEdge"))]
616 );
617 }
618 }
620 // Compare edge orientations for edge-only coupling
621 label compVal = -2;
623 // Perform a topological comparison.
624 switch (masterMap.type())
625 {
626 case 1:
627 {
628 if (sIndex < 0)
629 {
630 slaveMoveNewPoint[slaveI][0] = points_[cv0];
631 slaveMoveNewPoint[slaveI][1] = points_[cv1];
633 slaveMoveOldPoint[slaveI][0] = oldPoints_[cv0];
634 slaveMoveOldPoint[slaveI][1] = oldPoints_[cv1];
635 }
636 else
637 if (mEdge[0] == sEdge[0])
638 {
639 slaveOverRide = 1;
640 }
641 else
642 if (mEdge[1] == sEdge[0])
643 {
644 slaveOverRide = 2;
645 }
646 else
647 {
648 // Write out for for post-processing
649 writeVTK("mFace_" + Foam::name(fIndex), fIndex, 2);
651 FatalErrorIn
652 (
653 "\n"
654 "const changeMap "
655 "dynamicTopoFvMesh::collapseQuadFace\n"
656 "(\n"
657 " const label fIndex,\n"
658 " label overRideCase,\n"
659 " bool checkOnly,\n"
660 " bool forceOp\n"
661 ")\n"
662 )
663 << "Coupled collapse failed." << nl
664 << "Masters: " << nl
665 << checkEdgeIndex[1] << ": "
666 << checkEdge[1] << nl
667 << checkEdgeIndex[2] << ": "
668 << checkEdge[2] << nl
669 << "Slaves: " << nl
670 << readLabel(slaveMap.lookup("firstEdge")) << ": "
671 << sEdge[0] << nl
672 << readLabel(slaveMap.lookup("secondEdge")) << ": "
673 << sEdge[1] << nl
674 << abort(FatalError);
675 }
677 break;
678 }
680 case 2:
681 {
682 if (sIndex < 0)
683 {
684 slaveMoveNewPoint[slaveI][0] = points_[cv2];
685 slaveMoveNewPoint[slaveI][1] = points_[cv3];
687 slaveMoveOldPoint[slaveI][0] = oldPoints_[cv2];
688 slaveMoveOldPoint[slaveI][1] = oldPoints_[cv3];
689 }
690 else
691 if (mEdge[1] == sEdge[1])
692 {
693 slaveOverRide = 2;
694 }
695 else
696 if (mEdge[0] == sEdge[1])
697 {
698 slaveOverRide = 1;
699 }
700 else
701 {
702 // Write out for for post-processing
703 writeVTK("mFace_" + Foam::name(fIndex), fIndex, 2);
705 FatalErrorIn
706 (
707 "\n"
708 "const changeMap "
709 "dynamicTopoFvMesh::collapseQuadFace\n"
710 "(\n"
711 " const label fIndex,\n"
712 " label overRideCase,\n"
713 " bool checkOnly,\n"
714 " bool forceOp\n"
715 ")\n"
716 )
717 << "Coupled collapse failed." << nl
718 << "Masters: " << nl
719 << checkEdgeIndex[1] << ": "
720 << checkEdge[1] << nl
721 << checkEdgeIndex[2] << ": "
722 << checkEdge[2] << nl
723 << "Slaves: " << nl
724 << readLabel(slaveMap.lookup("firstEdge")) << ": "
725 << sEdge[0] << nl
726 << readLabel(slaveMap.lookup("secondEdge")) << ": "
727 << sEdge[1] << nl
728 << abort(FatalError);
729 }
731 break;
732 }
734 case 3:
735 {
736 if (sIndex < 0)
737 {
738 slaveMoveNewPoint[slaveI][0] =
739 (
740 0.5 * (points_[cv0] + points_[cv2])
741 );
743 slaveMoveNewPoint[slaveI][1] =
744 (
745 0.5 * (points_[cv1] + points_[cv3])
746 );
748 slaveMoveOldPoint[slaveI][0] =
749 (
750 0.5 * (oldPoints_[cv0] + oldPoints_[cv2])
751 );
753 slaveMoveOldPoint[slaveI][1] =
754 (
755 0.5 * (oldPoints_[cv1] + oldPoints_[cv3])
756 );
757 }
758 else
759 {
760 overRideCase = 3;
761 }
763 break;
764 }
765 }
767 if (sIndex < 0)
768 {
769 // Check edge orientation
770 compVal = edge::compare(mEdge[0], sEdge[0]);
772 // Swap components if necessary
773 if (compVal == -1)
774 {
775 Foam::Swap
776 (
777 slaveMoveNewPoint[slaveI][0],
778 slaveMoveNewPoint[slaveI][1]
779 );
781 Foam::Swap
782 (
783 slaveMoveOldPoint[slaveI][0],
784 slaveMoveOldPoint[slaveI][1]
785 );
786 }
787 else
788 if (compVal != 1)
789 {
790 FatalErrorIn
791 (
792 "\n"
793 "const changeMap "
794 "dynamicTopoFvMesh::collapseQuadFace\n"
795 "(\n"
796 " const label fIndex,\n"
797 " label overRideCase,\n"
798 " bool checkOnly,\n"
799 " bool forceOp\n"
800 ")\n"
801 )
802 << "Coupled topo-change for slave failed." << nl
803 << " Dissimilar edges: " << nl
804 << " mEdge: " << mEdge << nl
805 << " sEdge: " << sEdge << nl
806 << " masterMap type: " << masterMap.type() << nl
807 << abort(FatalError);
808 }
809 }
811 // Temporarily turn off coupledModification
812 unsetCoupledModification();
814 // Test the slave face
815 if (localCouple)
816 {
817 slaveMap =
818 (
819 collapseQuadFace(sIndex, slaveOverRide, true, forceOp)
820 );
821 }
822 else
823 if (procCouple)
824 {
825 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
827 if (sIndex < 0)
828 {
829 // Edge-based coupling
831 // Build a hull of cells and tri-faces
832 // that are connected to the slave edge
833 labelList slaveHullCells;
834 labelList slaveHullTriFaces;
836 meshOps::constructPrismHull
837 (
838 mag(sIndex),
839 sMesh.faces_,
840 sMesh.cells_,
841 sMesh.owner_,
842 sMesh.neighbour_,
843 sMesh.edgeFaces_,
844 slaveHullTriFaces,
845 slaveHullCells
846 );
848 bool infeasible = false;
850 // Keep track of resulting cell quality,
851 // if collapse is indeed feasible
852 scalar slaveCollapseQuality(GREAT);
854 // Check whether the collapse is possible.
855 if
856 (
857 checkCollapse
858 (
859 sMesh,
860 slaveHullTriFaces,
861 FixedList<label, 2>(-1),
862 FixedList<label, 2>(-1),
863 sEdge[0],
864 slaveMoveNewPoint[slaveI],
865 slaveMoveOldPoint[slaveI],
866 slaveCollapseQuality,
867 false,
868 forceOp
869 )
870 )
871 {
872 infeasible = true;
873 }
875 if (infeasible)
876 {
877 slaveMap.type() = 0;
878 }
879 else
880 {
881 slaveMap.type() = 1;
882 }
883 }
884 else
885 {
886 // Edge-based coupling
887 slaveMap =
888 (
889 sMesh.collapseQuadFace
890 (
891 sIndex,
892 slaveOverRide,
893 true,
894 forceOp
895 )
896 );
897 }
898 }
900 // Turn it back on.
901 setCoupledModification();
903 if (slaveMap.type() <= 0)
904 {
905 // Slave couldn't perform collapse.
906 map.type() = -2;
908 return map;
909 }
911 // Save index and patch for posterity
912 slaveMap.index() = sIndex;
913 slaveMap.patchIndex() = pI;
914 }
916 // Next collapse each slave face (for real this time...)
917 forAll(slaveMaps, slaveI)
918 {
919 // Alias for convenience...
920 changeMap& slaveMap = slaveMaps[slaveI];
922 label sIndex = slaveMap.index();
923 label pI = slaveMap.patchIndex();
924 label slaveOverRide = slaveMap.type();
926 // Temporarily turn off coupledModification
927 unsetCoupledModification();
929 // Collapse the slave
930 if (localCouple)
931 {
932 slaveMap =
933 (
934 collapseQuadFace(sIndex, slaveOverRide, false, forceOp)
935 );
936 }
937 else
938 {
939 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
941 if (sIndex < 0)
942 {
943 // Edge-based coupling
944 const coupleMap& cMap = recvMeshes_[pI].map();
946 // Fetch the slave edge
947 edge sEdge = sMesh.edges_[mag(sIndex)];
949 // Move points to new location,
950 // and update operation into coupleMap
951 sMesh.points_[sEdge[0]] = slaveMoveNewPoint[slaveI][0];
952 sMesh.points_[sEdge[1]] = slaveMoveNewPoint[slaveI][1];
954 sMesh.oldPoints_[sEdge[0]] = slaveMoveOldPoint[slaveI][0];
955 sMesh.oldPoints_[sEdge[1]] = slaveMoveOldPoint[slaveI][1];
957 cMap.pushOperation
958 (
959 sEdge[0],
960 coupleMap::MOVE_POINT,
961 slaveMoveNewPoint[slaveI][0],
962 slaveMoveOldPoint[slaveI][0]
963 );
965 cMap.pushOperation
966 (
967 sEdge[1],
968 coupleMap::MOVE_POINT,
969 slaveMoveNewPoint[slaveI][1],
970 slaveMoveOldPoint[slaveI][1]
971 );
973 // Force operation to succeed
974 slaveMap.type() = 1;
975 }
976 else
977 {
978 // Face-based coupling
979 slaveMap =
980 (
981 sMesh.collapseQuadFace
982 (
983 sIndex,
984 slaveOverRide,
985 false,
986 forceOp
987 )
988 );
989 }
990 }
992 // Turn it back on.
993 setCoupledModification();
995 // The final operation has to succeed.
996 if (slaveMap.type() <= 0)
997 {
998 FatalErrorIn
999 (
1000 "\n"
1001 "const changeMap "
1002 "dynamicTopoFvMesh::collapseQuadFace\n"
1003 "(\n"
1004 " const label fIndex,\n"
1005 " label overRideCase,\n"
1006 " bool checkOnly,\n"
1007 " bool forceOp\n"
1008 ")\n"
1009 )
1010 << "Coupled topo-change for slave failed." << nl
1011 << " Face: " << fIndex << ": " << fCheck << nl
1012 << " Slave index: " << sIndex << nl
1013 << " Patch index: " << pI << nl
1014 << " Type: " << slaveMap.type() << nl
1015 << abort(FatalError);
1016 }
1018 // Save index and patch for posterity
1019 slaveMap.index() = sIndex;
1020 slaveMap.patchIndex() = pI;
1021 }
1022 }
1024 // Build a hull of cells and tri-faces that are connected to each edge
1025 FixedList<labelList, 2> hullCells;
1026 FixedList<labelList, 2> hullTriFaces;
1028 meshOps::constructPrismHull
1029 (
1030 checkEdgeIndex[1],
1031 faces_,
1032 cells_,
1033 owner_,
1034 neighbour_,
1035 edgeFaces_,
1036 hullTriFaces[0],
1037 hullCells[0]
1038 );
1040 meshOps::constructPrismHull
1041 (
1042 checkEdgeIndex[2],
1043 faces_,
1044 cells_,
1045 owner_,
1046 neighbour_,
1047 edgeFaces_,
1048 hullTriFaces[1],
1049 hullCells[1]
1050 );
1052 // Determine the neighbouring cells
1053 label c0 = owner_[fIndex], c1 = neighbour_[fIndex];
1055 // Define variables for the prism-face calculation
1056 FixedList<label,2> c0BdyIndex(-1), c0IntIndex(-1);
1057 FixedList<label,2> c1BdyIndex(-1), c1IntIndex(-1);
1058 FixedList<face,2> c0BdyFace(face(3)), c0IntFace(face(4));
1059 FixedList<face,2> c1BdyFace(face(3)), c1IntFace(face(4));
1061 // Find the prism-faces
1062 meshOps::findPrismFaces
1063 (
1064 fIndex,
1065 c0,
1066 faces_,
1067 cells_,
1068 neighbour_,
1069 c0BdyFace,
1070 c0BdyIndex,
1071 c0IntFace,
1072 c0IntIndex
1073 );
1075 if (c1 != -1)
1076 {
1077 meshOps::findPrismFaces
1078 (
1079 fIndex,
1080 c1,
1081 faces_,
1082 cells_,
1083 neighbour_,
1084 c1BdyFace,
1085 c1BdyIndex,
1086 c1IntFace,
1087 c1IntIndex
1088 );
1089 }
1091 // Determine if either edge belongs to a boundary
1092 FixedList<bool, 2> nBoundCurves(false), edgeBoundary(false);
1094 edgeBoundary[0] = (whichEdgePatch(checkEdgeIndex[1]) > -1);
1095 edgeBoundary[1] = (whichEdgePatch(checkEdgeIndex[2]) > -1);
1097 // Decide on collapseCase
1098 label collapseCase = -1;
1100 if (edgeBoundary[0] && !edgeBoundary[1])
1101 {
1102 collapseCase = 1;
1103 }
1104 else
1105 if (!edgeBoundary[0] && edgeBoundary[1])
1106 {
1107 collapseCase = 2;
1108 }
1109 else
1110 if (edgeBoundary[0] && edgeBoundary[1])
1111 {
1112 if (c1 != -1)
1113 {
1114 if (debug > 2)
1115 {
1116 WarningIn
1117 (
1118 "\n"
1119 "const changeMap "
1120 "dynamicTopoFvMesh::collapseQuadFace\n"
1121 "(\n"
1122 " const label fIndex,\n"
1123 " label overRideCase,\n"
1124 " bool checkOnly,\n"
1125 " bool forceOp\n"
1126 ")\n"
1127 ) << "Collapsing an internal face that "
1128 << "lies on two boundary patches. "
1129 << "Face: " << fIndex << ": " << faces_[fIndex]
1130 << endl;
1131 }
1133 // Bail out for now. If proximity based refinement is
1134 // switched on, mesh may be sliced at this point.
1135 return map;
1136 }
1138 // Check if either edge lies on a bounding curve.
1139 if (checkBoundingCurve(checkEdgeIndex[1]))
1140 {
1141 nBoundCurves[0] = true;
1142 }
1144 if (checkBoundingCurve(checkEdgeIndex[2]))
1145 {
1146 nBoundCurves[1] = true;
1147 }
1149 // Collapse towards a bounding curve
1150 if (nBoundCurves[0] && !nBoundCurves[1])
1151 {
1152 collapseCase = 1;
1153 }
1154 else
1155 if (!nBoundCurves[0] && nBoundCurves[1])
1156 {
1157 collapseCase = 2;
1158 }
1159 else
1160 if (!nBoundCurves[0] && !nBoundCurves[1])
1161 {
1162 // No bounding curves. Collapse to mid-point.
1163 collapseCase = 3;
1164 }
1165 else
1166 {
1167 // Two bounding curves? This might cause pinching.
1168 // Bail out for now.
1169 return map;
1170 }
1171 }
1172 else
1173 {
1174 // Looks like this is an interior face.
1175 // Collapse case [3] by default
1176 collapseCase = 3;
1177 }
1179 // Perform an override if requested.
1180 if (overRideCase != -1)
1181 {
1182 collapseCase = overRideCase;
1183 }
1185 // Configure the new point-positions
1186 FixedList<bool, 2> check(false);
1187 FixedList<FixedList<label, 2>, 2> checkPoints(FixedList<label, 2>(-1));
1188 FixedList<point, 2> newPoint(vector::zero);
1189 FixedList<point, 2> oldPoint(vector::zero);
1191 // Replacement check points
1192 FixedList<label,2> original(-1), replacement(-1);
1194 switch (collapseCase)
1195 {
1196 case 1:
1197 {
1198 // Collapse to the first node
1199 original[0] = cv2; original[1] = cv3;
1200 replacement[0] = cv0; replacement[1] = cv1;
1202 newPoint[0] = points_[replacement[0]];
1203 newPoint[1] = points_[replacement[1]];
1205 oldPoint[0] = oldPoints_[replacement[0]];
1206 oldPoint[1] = oldPoints_[replacement[1]];
1208 // Define check-points
1209 check[1] = true;
1210 checkPoints[1][0] = original[0];
1211 checkPoints[1][1] = original[1];
1213 break;
1214 }
1216 case 2:
1217 {
1218 // Collapse to the second node
1219 original[0] = cv0; original[1] = cv1;
1220 replacement[0] = cv2; replacement[1] = cv3;
1222 newPoint[0] = points_[replacement[0]];
1223 newPoint[1] = points_[replacement[1]];
1225 oldPoint[0] = oldPoints_[replacement[0]];
1226 oldPoint[1] = oldPoints_[replacement[1]];
1228 // Define check-points
1229 check[0] = true;
1230 checkPoints[0][0] = original[0];
1231 checkPoints[0][1] = original[1];
1233 break;
1234 }
1236 case 3:
1237 {
1238 // Collapse to the mid-point
1239 original[0] = cv0; original[1] = cv1;
1240 replacement[0] = cv2; replacement[1] = cv3;
1242 // Define new point-positions
1243 newPoint[0] =
1244 (
1245 0.5 *
1246 (
1247 points_[original[0]]
1248 + points_[replacement[0]]
1249 )
1250 );
1252 newPoint[1] =
1253 (
1254 0.5 *
1255 (
1256 points_[original[1]]
1257 + points_[replacement[1]]
1258 )
1259 );
1261 // Define old point-positions
1262 oldPoint[0] =
1263 (
1264 0.5 *
1265 (
1266 oldPoints_[original[0]]
1267 + oldPoints_[replacement[0]]
1268 )
1269 );
1271 oldPoint[1] =
1272 (
1273 0.5 *
1274 (
1275 oldPoints_[original[1]]
1276 + oldPoints_[replacement[1]]
1277 )
1278 );
1280 // Define check-points
1281 check[0] = true;
1282 checkPoints[0][0] = original[0];
1283 checkPoints[0][1] = original[1];
1285 check[1] = true;
1286 checkPoints[1][0] = replacement[0];
1287 checkPoints[1][1] = replacement[1];
1289 break;
1290 }
1292 default:
1293 {
1294 // Don't think this will ever happen.
1295 FatalErrorIn
1296 (
1297 "\n"
1298 "const changeMap "
1299 "dynamicTopoFvMesh::collapseQuadFace\n"
1300 "(\n"
1301 " const label fIndex,\n"
1302 " label overRideCase,\n"
1303 " bool checkOnly,\n"
1304 " bool forceOp\n"
1305 ")\n"
1306 )
1307 << "Edge: " << fIndex << ": " << faces_[fIndex]
1308 << ". Couldn't decide on collapseCase."
1309 << abort(FatalError);
1311 break;
1312 }
1313 }
1315 // Keep track of resulting cell quality,
1316 // if collapse is indeed feasible
1317 scalar collapseQuality(GREAT);
1319 // Check collapsibility of cells around edges
1320 // with the re-configured point
1321 forAll(check, indexI)
1322 {
1323 if (!check[indexI])
1324 {
1325 continue;
1326 }
1328 // Check whether the collapse is possible.
1329 if
1330 (
1331 checkCollapse
1332 (
1333 (*this),
1334 hullTriFaces[indexI],
1335 c0BdyIndex,
1336 c1BdyIndex,
1337 checkPoints[indexI],
1338 newPoint,
1339 oldPoint,
1340 collapseQuality,
1341 (c1 != -1),
1342 forceOp
1343 )
1344 )
1345 {
1346 map.type() = 0;
1347 return map;
1348 }
1349 }
1351 // Add a map entry of the replacements as 'addedPoints'
1352 // - Used in coupled mapping
1353 map.addPoint(replacement[0]);
1354 map.addPoint(replacement[1]);
1356 // Are we only performing checks?
1357 if (checkOnly)
1358 {
1359 map.type() = collapseCase;
1361 if (debug > 2)
1362 {
1363 Pout<< "Face: " << fIndex
1364 << ":: " << faces_[fIndex] << nl
1365 << " collapseCase determined to be: " << collapseCase << nl
1366 << " Resulting quality: " << collapseQuality << nl
1367 << " edgeBoundary: " << edgeBoundary << nl
1368 << " nBoundCurves: " << nBoundCurves << nl
1369 << " collapsePoints: " << original << nl
1370 << " replacePoints: " << replacement << nl
1371 << endl;
1372 }
1374 return map;
1375 }
1377 if (debug > 1)
1378 {
1379 const labelList& fE = faceEdges_[fIndex];
1381 Pout<< nl << nl
1382 << "Face: " << fIndex << ": " << faces_[fIndex] << nl
1383 << "faceEdges: " << fE << " is to be collapsed. "
1384 << nl;
1386 Pout<< " On SubMesh: " << Switch::asText(isSubMesh_) << nl;
1387 Pout<< " coupledModification: " << coupledModification_ << nl;
1389 label epIndex = whichPatch(fIndex);
1391 const polyBoundaryMesh& boundary = boundaryMesh();
1393 Pout<< " Patch: ";
1395 if (epIndex == -1)
1396 {
1397 Pout<< "Internal" << nl;
1398 }
1399 else
1400 if (epIndex < boundary.size())
1401 {
1402 Pout<< boundary[epIndex].name() << nl;
1403 }
1404 else
1405 {
1406 Pout<< " New patch: " << epIndex << endl;
1407 }
1409 if (debug > 2)
1410 {
1411 Pout<< nl
1412 << "~~~~~~~~~~~~~~~~~~~~~~~~~" << nl
1413 << "Hulls before modification" << nl
1414 << "~~~~~~~~~~~~~~~~~~~~~~~~~" << nl;
1416 if (debug > 3)
1417 {
1418 Pout<< "Cell [0] Boundary faces: " << nl
1419 << " Face: "<< c0BdyIndex[0]
1420 << "::" << c0BdyFace[0] << nl
1421 << " Face: "<< c0BdyIndex[1]
1422 << "::" << c0BdyFace[1] << nl;
1424 Pout<< "Cell [0] Interior faces: " << nl
1425 << " Face: "<< c0IntIndex[0]
1426 << "::" << c0IntFace[0] << nl
1427 << " Face: "<< c0IntIndex[1]
1428 << "::" << c0IntFace[1] << nl;
1430 if (c1 != -1)
1431 {
1432 Pout<< "Cell [1] Boundary faces: " << nl
1433 << " Face: "<< c1BdyIndex[0]
1434 << "::" << c1BdyFace[0] << nl
1435 << " Face: "<< c1BdyIndex[1]
1436 << "::" << c1BdyFace[1] << nl;
1438 Pout<< "Cell [1] Interior faces: " << nl
1439 << " Face: "<< c1IntIndex[0]
1440 << "::" << c1IntFace[0] << nl
1441 << " Face: "<< c1IntIndex[1]
1442 << "::" << c1IntFace[1] << nl;
1443 }
1444 }
1446 Pout<< nl << "Cells belonging to first Edge Hull: "
1447 << hullCells[0] << nl;
1449 forAll(hullCells[0],cellI)
1450 {
1451 const cell& firstCurCell = cells_[hullCells[0][cellI]];
1453 Pout<< "Cell: " << hullCells[0][cellI]
1454 << ": " << firstCurCell << nl;
1456 forAll(firstCurCell,faceI)
1457 {
1458 Pout<< " Face: " << firstCurCell[faceI]
1459 << " : " << faces_[firstCurCell[faceI]]
1460 << " fE: " << faceEdges_[firstCurCell[faceI]]
1461 << nl;
1463 if (debug > 3)
1464 {
1465 const labelList& fE = faceEdges_[firstCurCell[faceI]];
1467 forAll(fE, edgeI)
1468 {
1469 Pout<< " Edge: " << fE[edgeI]
1470 << " : " << edges_[fE[edgeI]]
1471 << nl;
1472 }
1473 }
1474 }
1475 }
1477 const labelList& firstEdgeFaces = edgeFaces_[checkEdgeIndex[1]];
1479 Pout<< nl << "First Edge Face Hull: "
1480 << firstEdgeFaces << nl;
1482 forAll(firstEdgeFaces,indexI)
1483 {
1484 Pout<< " Face: " << firstEdgeFaces[indexI]
1485 << " : " << faces_[firstEdgeFaces[indexI]]
1486 << " fE: " << faceEdges_[firstEdgeFaces[indexI]]
1487 << nl;
1489 if (debug > 3)
1490 {
1491 const labelList& fE = faceEdges_[firstEdgeFaces[indexI]];
1493 forAll(fE, edgeI)
1494 {
1495 Pout<< " Edge: " << fE[edgeI]
1496 << " : " << edges_[fE[edgeI]]
1497 << nl;
1498 }
1499 }
1500 }
1502 Pout<< nl << "Cells belonging to second Edge Hull: "
1503 << hullCells[1] << nl;
1505 forAll(hullCells[1], cellI)
1506 {
1507 const cell& secondCurCell = cells_[hullCells[1][cellI]];
1509 Pout<< "Cell: " << hullCells[1][cellI]
1510 << ": " << secondCurCell << nl;
1512 forAll(secondCurCell, faceI)
1513 {
1514 Pout<< " Face: " << secondCurCell[faceI]
1515 << " : " << faces_[secondCurCell[faceI]]
1516 << " fE: " << faceEdges_[secondCurCell[faceI]]
1517 << nl;
1519 if (debug > 3)
1520 {
1521 const labelList& fE = faceEdges_[secondCurCell[faceI]];
1523 forAll(fE, edgeI)
1524 {
1525 Pout<< " Edge: " << fE[edgeI]
1526 << " : " << edges_[fE[edgeI]]
1527 << nl;
1528 }
1529 }
1530 }
1531 }
1533 const labelList& secondEdgeFaces = edgeFaces_[checkEdgeIndex[2]];
1535 Pout<< nl << "Second Edge Face Hull: "
1536 << secondEdgeFaces << nl;
1538 forAll(secondEdgeFaces, indexI)
1539 {
1540 Pout<< " Face: " << secondEdgeFaces[indexI]
1541 << " : " << faces_[secondEdgeFaces[indexI]]
1542 << " fE: " << faceEdges_[secondEdgeFaces[indexI]]
1543 << nl;
1545 if (debug > 3)
1546 {
1547 const labelList& fE = faceEdges_[secondEdgeFaces[indexI]];
1549 forAll(fE, edgeI)
1550 {
1551 Pout<< " Edge: " << fE[edgeI]
1552 << " : " << edges_[fE[edgeI]]
1553 << nl;
1554 }
1555 }
1556 }
1558 // Write out VTK files prior to change
1559 if (debug > 3)
1560 {
1561 labelHashSet vtkCells;
1563 forAll(hullCells[0], cellI)
1564 {
1565 if (!vtkCells.found(hullCells[0][cellI]))
1566 {
1567 vtkCells.insert(hullCells[0][cellI]);
1568 }
1569 }
1571 forAll(hullCells[1], cellI)
1572 {
1573 if (!vtkCells.found(hullCells[1][cellI]))
1574 {
1575 vtkCells.insert(hullCells[1][cellI]);
1576 }
1577 }
1579 writeVTK
1580 (
1581 Foam::name(fIndex)
1582 + "_Collapse_0",
1583 vtkCells.toc()
1584 );
1585 }
1586 }
1587 }
1589 // Edges / Quad-faces to throw or keep during collapse
1590 FixedList<label,2> ends(-1);
1591 FixedList<label,2> faceToKeep(-1), faceToThrow(-1);
1592 FixedList<label,4> edgeToKeep(-1), edgeToThrow(-1);
1594 // Maintain a list of modified faces for mapping
1595 DynamicList<label> modifiedFaces(10);
1597 // Case 2 & 3 use identical connectivity,
1598 // but different point locations
1599 if (collapseCase == 2 || collapseCase == 3)
1600 {
1601 const labelList& firstEdgeFaces = edgeFaces_[checkEdgeIndex[1]];
1603 // Collapse to the second node...
1604 forAll(firstEdgeFaces,faceI)
1605 {
1606 // Replace point indices on faces.
1607 meshOps::replaceLabel
1608 (
1609 cv0,
1610 cv2,
1611 faces_[firstEdgeFaces[faceI]]
1612 );
1614 meshOps::replaceLabel
1615 (
1616 cv1,
1617 cv3,
1618 faces_[firstEdgeFaces[faceI]]
1619 );
1621 // Add an entry for mapping
1622 if (findIndex(modifiedFaces, firstEdgeFaces[faceI]) == -1)
1623 {
1624 modifiedFaces.append(firstEdgeFaces[faceI]);
1625 }
1627 // Determine the quad-face in cell[0] & cell[1]
1628 // that belongs to firstEdgeFaces
1629 if (firstEdgeFaces[faceI] == c0IntIndex[0])
1630 {
1631 faceToKeep[0] = c0IntIndex[1];
1632 faceToThrow[0] = c0IntIndex[0];
1633 }
1635 if (firstEdgeFaces[faceI] == c0IntIndex[1])
1636 {
1637 faceToKeep[0] = c0IntIndex[0];
1638 faceToThrow[0] = c0IntIndex[1];
1639 }
1641 if (c1 != -1)
1642 {
1643 if (firstEdgeFaces[faceI] == c1IntIndex[0])
1644 {
1645 faceToKeep[1] = c1IntIndex[1];
1646 faceToThrow[1] = c1IntIndex[0];
1647 }
1649 if (firstEdgeFaces[faceI] == c1IntIndex[1])
1650 {
1651 faceToKeep[1] = c1IntIndex[0];
1652 faceToThrow[1] = c1IntIndex[1];
1653 }
1654 }
1655 }
1657 // Find common edges between quad / tri faces...
1658 meshOps::findCommonEdge
1659 (
1660 c0BdyIndex[0],
1661 faceToKeep[0],
1662 faceEdges_,
1663 edgeToKeep[0]
1664 );
1666 meshOps::findCommonEdge
1667 (
1668 c0BdyIndex[1],
1669 faceToKeep[0],
1670 faceEdges_,
1671 edgeToKeep[1]
1672 );
1674 meshOps::findCommonEdge
1675 (
1676 c0BdyIndex[0],
1677 faceToThrow[0],
1678 faceEdges_,
1679 edgeToThrow[0]
1680 );
1682 meshOps::findCommonEdge
1683 (
1684 c0BdyIndex[1],
1685 faceToThrow[0],
1686 faceEdges_,
1687 edgeToThrow[1]
1688 );
1690 // Size down edgeFaces for the ends.
1691 meshOps::findCommonEdge
1692 (
1693 faceToThrow[0],
1694 faceToKeep[0],
1695 faceEdges_,
1696 ends[0]
1697 );
1699 meshOps::sizeDownList
1700 (
1701 faceToThrow[0],
1702 edgeFaces_[ends[0]]
1703 );
1705 if (c1 != -1)
1706 {
1707 meshOps::findCommonEdge
1708 (
1709 c1BdyIndex[0],
1710 faceToKeep[1],
1711 faceEdges_,
1712 edgeToKeep[2]
1713 );
1715 meshOps::findCommonEdge
1716 (
1717 c1BdyIndex[1],
1718 faceToKeep[1],
1719 faceEdges_,
1720 edgeToKeep[3]
1721 );
1723 meshOps::findCommonEdge
1724 (
1725 c1BdyIndex[0],
1726 faceToThrow[1],
1727 faceEdges_,
1728 edgeToThrow[2]
1729 );
1731 meshOps::findCommonEdge
1732 (
1733 c1BdyIndex[1],
1734 faceToThrow[1],
1735 faceEdges_,
1736 edgeToThrow[3]
1737 );
1739 // Size down edgeFaces for the ends.
1740 meshOps::findCommonEdge
1741 (
1742 faceToThrow[1],
1743 faceToKeep[1],
1744 faceEdges_,
1745 ends[1]
1746 );
1748 meshOps::sizeDownList
1749 (
1750 faceToThrow[1],
1751 edgeFaces_[ends[1]]
1752 );
1753 }
1755 // Correct edgeFaces for triangular faces...
1756 forAll(edgeToThrow, indexI)
1757 {
1758 if (edgeToThrow[indexI] == -1)
1759 {
1760 continue;
1761 }
1763 const labelList& eF = edgeFaces_[edgeToThrow[indexI]];
1765 label origTriFace = -1, retTriFace = -1;
1767 // Find the original triangular face index
1768 forAll(eF, faceI)
1769 {
1770 if (eF[faceI] == c0BdyIndex[0])
1771 {
1772 origTriFace = c0BdyIndex[0];
1773 break;
1774 }
1776 if (eF[faceI] == c0BdyIndex[1])
1777 {
1778 origTriFace = c0BdyIndex[1];
1779 break;
1780 }
1782 if (eF[faceI] == c1BdyIndex[0])
1783 {
1784 origTriFace = c1BdyIndex[0];
1785 break;
1786 }
1788 if (eF[faceI] == c1BdyIndex[1])
1789 {
1790 origTriFace = c1BdyIndex[1];
1791 break;
1792 }
1793 }
1795 // Find the retained triangular face index
1796 forAll(eF, faceI)
1797 {
1798 if
1799 (
1800 (eF[faceI] != origTriFace) &&
1801 (eF[faceI] != faceToThrow[0]) &&
1802 (eF[faceI] != faceToThrow[1])
1803 )
1804 {
1805 retTriFace = eF[faceI];
1806 break;
1807 }
1808 }
1810 // Finally replace the face index
1811 if (retTriFace == -1)
1812 {
1813 // Couldn't find a retained face.
1814 // This must be a boundary edge, so size-down instead.
1815 meshOps::sizeDownList
1816 (
1817 origTriFace,
1818 edgeFaces_[edgeToKeep[indexI]]
1819 );
1820 }
1821 else
1822 {
1823 meshOps::replaceLabel
1824 (
1825 origTriFace,
1826 retTriFace,
1827 edgeFaces_[edgeToKeep[indexI]]
1828 );
1830 meshOps::replaceLabel
1831 (
1832 edgeToThrow[indexI],
1833 edgeToKeep[indexI],
1834 faceEdges_[retTriFace]
1835 );
1836 }
1837 }
1839 // Correct faceEdges / edgeFaces for quad-faces...
1840 forAll(firstEdgeFaces,faceI)
1841 {
1842 meshOps::replaceLabel
1843 (
1844 checkEdgeIndex[1],
1845 checkEdgeIndex[2],
1846 faceEdges_[firstEdgeFaces[faceI]]
1847 );
1849 // Renumber the edges on this face
1850 const labelList& fE = faceEdges_[firstEdgeFaces[faceI]];
1852 forAll(fE, edgeI)
1853 {
1854 if (edges_[fE[edgeI]].start() == cv0)
1855 {
1856 edges_[fE[edgeI]].start() = cv2;
1857 }
1859 if (edges_[fE[edgeI]].end() == cv0)
1860 {
1861 edges_[fE[edgeI]].end() = cv2;
1862 }
1864 if (edges_[fE[edgeI]].start() == cv1)
1865 {
1866 edges_[fE[edgeI]].start() = cv3;
1867 }
1869 if (edges_[fE[edgeI]].end() == cv1)
1870 {
1871 edges_[fE[edgeI]].end() = cv3;
1872 }
1873 }
1875 // Size-up edgeFaces for the replacement
1876 if
1877 (
1878 (firstEdgeFaces[faceI] != faceToThrow[0]) &&
1879 (firstEdgeFaces[faceI] != faceToThrow[1]) &&
1880 (firstEdgeFaces[faceI] != fIndex)
1881 )
1882 {
1883 meshOps::sizeUpList
1884 (
1885 firstEdgeFaces[faceI],
1886 edgeFaces_[checkEdgeIndex[2]]
1887 );
1888 }
1889 }
1891 // Remove the current face from the replacement edge
1892 meshOps::sizeDownList
1893 (
1894 fIndex,
1895 edgeFaces_[checkEdgeIndex[2]]
1896 );
1898 // Replace point labels on all triangular boundary faces.
1899 forAll(hullCells[0],cellI)
1900 {
1901 const cell& cellToCheck = cells_[hullCells[0][cellI]];
1903 forAll(cellToCheck,faceI)
1904 {
1905 face& faceToCheck = faces_[cellToCheck[faceI]];
1907 if (faceToCheck.size() == 3)
1908 {
1909 forAll(faceToCheck,pointI)
1910 {
1911 if (faceToCheck[pointI] == cv0)
1912 {
1913 faceToCheck[pointI] = cv2;
1914 }
1916 if (faceToCheck[pointI] == cv1)
1917 {
1918 faceToCheck[pointI] = cv3;
1919 }
1920 }
1921 }
1922 }
1923 }
1925 // Now that we're done with all edges, remove them.
1926 removeEdge(checkEdgeIndex[0]);
1927 removeEdge(checkEdgeIndex[1]);
1928 removeEdge(checkEdgeIndex[3]);
1930 // Update map
1931 map.removeEdge(checkEdgeIndex[0]);
1932 map.removeEdge(checkEdgeIndex[1]);
1933 map.removeEdge(checkEdgeIndex[2]);
1935 forAll(edgeToThrow, indexI)
1936 {
1937 if (edgeToThrow[indexI] == -1)
1938 {
1939 continue;
1940 }
1942 removeEdge(edgeToThrow[indexI]);
1944 // Update map
1945 map.removeEdge(edgeToThrow[indexI]);
1946 }
1948 // Delete the two points...
1949 removePoint(cv0);
1950 removePoint(cv1);
1952 // Update map
1953 map.removePoint(cv0);
1954 map.removePoint(cv1);
1955 }
1956 else
1957 {
1958 const labelList& secondEdgeFaces = edgeFaces_[checkEdgeIndex[2]];
1960 // Collapse to the first node
1961 forAll(secondEdgeFaces,faceI)
1962 {
1963 // Replace point indices on faces.
1964 meshOps::replaceLabel
1965 (
1966 cv2,
1967 cv0,
1968 faces_[secondEdgeFaces[faceI]]
1969 );
1971 meshOps::replaceLabel
1972 (
1973 cv3,
1974 cv1,
1975 faces_[secondEdgeFaces[faceI]]
1976 );
1978 // Add an entry for mapping
1979 if (findIndex(modifiedFaces, secondEdgeFaces[faceI]) == -1)
1980 {
1981 modifiedFaces.append(secondEdgeFaces[faceI]);
1982 }
1984 // Determine the quad-face(s) in cell[0] & cell[1]
1985 // that belongs to secondEdgeFaces
1986 if (secondEdgeFaces[faceI] == c0IntIndex[0])
1987 {
1988 faceToKeep[0] = c0IntIndex[1];
1989 faceToThrow[0] = c0IntIndex[0];
1990 }
1992 if (secondEdgeFaces[faceI] == c0IntIndex[1])
1993 {
1994 faceToKeep[0] = c0IntIndex[0];
1995 faceToThrow[0] = c0IntIndex[1];
1996 }
1998 if (c1 != -1)
1999 {
2000 if (secondEdgeFaces[faceI] == c1IntIndex[0])
2001 {
2002 faceToKeep[1] = c1IntIndex[1];
2003 faceToThrow[1] = c1IntIndex[0];
2004 }
2006 if (secondEdgeFaces[faceI] == c1IntIndex[1])
2007 {
2008 faceToKeep[1] = c1IntIndex[0];
2009 faceToThrow[1] = c1IntIndex[1];
2010 }
2011 }
2012 }
2014 // Find common edges between quad / tri faces...
2015 meshOps::findCommonEdge
2016 (
2017 c0BdyIndex[0],
2018 faceToKeep[0],
2019 faceEdges_,
2020 edgeToKeep[0]
2021 );
2023 meshOps::findCommonEdge
2024 (
2025 c0BdyIndex[1],
2026 faceToKeep[0],
2027 faceEdges_,
2028 edgeToKeep[1]
2029 );
2031 meshOps::findCommonEdge
2032 (
2033 c0BdyIndex[0],
2034 faceToThrow[0],
2035 faceEdges_,
2036 edgeToThrow[0]
2037 );
2039 meshOps::findCommonEdge
2040 (
2041 c0BdyIndex[1],
2042 faceToThrow[0],
2043 faceEdges_,
2044 edgeToThrow[1]
2045 );
2047 // Size down edgeFaces for the ends.
2048 meshOps::findCommonEdge
2049 (
2050 faceToThrow[0],
2051 faceToKeep[0],
2052 faceEdges_,
2053 ends[0]
2054 );
2056 meshOps::sizeDownList
2057 (
2058 faceToThrow[0],
2059 edgeFaces_[ends[0]]
2060 );
2062 if (c1 != -1)
2063 {
2064 meshOps::findCommonEdge
2065 (
2066 c1BdyIndex[0],
2067 faceToKeep[1],
2068 faceEdges_,
2069 edgeToKeep[2]
2070 );
2072 meshOps::findCommonEdge
2073 (
2074 c1BdyIndex[1],
2075 faceToKeep[1],
2076 faceEdges_,
2077 edgeToKeep[3]
2078 );
2080 meshOps::findCommonEdge
2081 (
2082 c1BdyIndex[0],
2083 faceToThrow[1],
2084 faceEdges_,
2085 edgeToThrow[2]
2086 );
2088 meshOps::findCommonEdge
2089 (
2090 c1BdyIndex[1],
2091 faceToThrow[1],
2092 faceEdges_,
2093 edgeToThrow[3]
2094 );
2096 // Size down edgeFaces for the ends.
2097 meshOps::findCommonEdge
2098 (
2099 faceToThrow[1],
2100 faceToKeep[1],
2101 faceEdges_,
2102 ends[1]
2103 );
2105 meshOps::sizeDownList
2106 (
2107 faceToThrow[1],
2108 edgeFaces_[ends[1]]
2109 );
2110 }
2112 // Correct edgeFaces for triangular faces...
2113 forAll(edgeToThrow, indexI)
2114 {
2115 if (edgeToThrow[indexI] == -1)
2116 {
2117 continue;
2118 }
2120 const labelList& eF = edgeFaces_[edgeToThrow[indexI]];
2122 label origTriFace = -1, retTriFace = -1;
2124 // Find the original triangular face index
2125 forAll(eF, faceI)
2126 {
2127 if (eF[faceI] == c0BdyIndex[0])
2128 {
2129 origTriFace = c0BdyIndex[0];
2130 break;
2131 }
2133 if (eF[faceI] == c0BdyIndex[1])
2134 {
2135 origTriFace = c0BdyIndex[1];
2136 break;
2137 }
2139 if (eF[faceI] == c1BdyIndex[0])
2140 {
2141 origTriFace = c1BdyIndex[0];
2142 break;
2143 }
2145 if (eF[faceI] == c1BdyIndex[1])
2146 {
2147 origTriFace = c1BdyIndex[1];
2148 break;
2149 }
2150 }
2152 // Find the retained triangular face index
2153 forAll(eF, faceI)
2154 {
2155 if
2156 (
2157 (eF[faceI] != origTriFace) &&
2158 (eF[faceI] != faceToThrow[0]) &&
2159 (eF[faceI] != faceToThrow[1])
2160 )
2161 {
2162 retTriFace = eF[faceI];
2163 break;
2164 }
2165 }
2167 // Finally replace the face/edge indices
2168 if (retTriFace == -1)
2169 {
2170 // Couldn't find a retained face.
2171 // This must be a boundary edge, so size-down instead.
2172 meshOps::sizeDownList
2173 (
2174 origTriFace,
2175 edgeFaces_[edgeToKeep[indexI]]
2176 );
2177 }
2178 else
2179 {
2180 meshOps::replaceLabel
2181 (
2182 origTriFace,
2183 retTriFace,
2184 edgeFaces_[edgeToKeep[indexI]]
2185 );
2187 meshOps::replaceLabel
2188 (
2189 edgeToThrow[indexI],
2190 edgeToKeep[indexI],
2191 faceEdges_[retTriFace]
2192 );
2193 }
2194 }
2196 // Correct faceEdges / edgeFaces for quad-faces...
2197 forAll(secondEdgeFaces,faceI)
2198 {
2199 meshOps::replaceLabel
2200 (
2201 checkEdgeIndex[2],
2202 checkEdgeIndex[1],
2203 faceEdges_[secondEdgeFaces[faceI]]
2204 );
2206 // Renumber the edges on this face
2207 const labelList& fE = faceEdges_[secondEdgeFaces[faceI]];
2209 forAll(fE, edgeI)
2210 {
2211 if (edges_[fE[edgeI]].start() == cv2)
2212 {
2213 edges_[fE[edgeI]].start() = cv0;
2214 }
2216 if (edges_[fE[edgeI]].end() == cv2)
2217 {
2218 edges_[fE[edgeI]].end() = cv0;
2219 }
2221 if (edges_[fE[edgeI]].start() == cv3)
2222 {
2223 edges_[fE[edgeI]].start() = cv1;
2224 }
2226 if (edges_[fE[edgeI]].end() == cv3)
2227 {
2228 edges_[fE[edgeI]].end() = cv1;
2229 }
2230 }
2232 // Size-up edgeFaces for the replacement
2233 if
2234 (
2235 (secondEdgeFaces[faceI] != faceToThrow[0]) &&
2236 (secondEdgeFaces[faceI] != faceToThrow[1]) &&
2237 (secondEdgeFaces[faceI] != fIndex)
2238 )
2239 {
2240 meshOps::sizeUpList
2241 (
2242 secondEdgeFaces[faceI],
2243 edgeFaces_[checkEdgeIndex[1]]
2244 );
2245 }
2246 }
2248 // Remove the current face from the replacement edge
2249 meshOps::sizeDownList
2250 (
2251 fIndex,
2252 edgeFaces_[checkEdgeIndex[1]]
2253 );
2255 // Replace point labels on all triangular boundary faces.
2256 forAll(hullCells[1], cellI)
2257 {
2258 const cell& cellToCheck = cells_[hullCells[1][cellI]];
2260 forAll(cellToCheck, faceI)
2261 {
2262 face& faceToCheck = faces_[cellToCheck[faceI]];
2264 if (faceToCheck.size() == 3)
2265 {
2266 forAll(faceToCheck, pointI)
2267 {
2268 if (faceToCheck[pointI] == cv2)
2269 {
2270 faceToCheck[pointI] = cv0;
2271 }
2273 if (faceToCheck[pointI] == cv3)
2274 {
2275 faceToCheck[pointI] = cv1;
2276 }
2277 }
2278 }
2279 }
2280 }
2282 // Now that we're done with all edges, remove them.
2283 removeEdge(checkEdgeIndex[0]);
2284 removeEdge(checkEdgeIndex[2]);
2285 removeEdge(checkEdgeIndex[3]);
2287 // Update map
2288 map.removeEdge(checkEdgeIndex[0]);
2289 map.removeEdge(checkEdgeIndex[2]);
2290 map.removeEdge(checkEdgeIndex[3]);
2292 forAll(edgeToThrow, indexI)
2293 {
2294 if (edgeToThrow[indexI] == -1)
2295 {
2296 continue;
2297 }
2299 removeEdge(edgeToThrow[indexI]);
2301 // Update map
2302 map.removeEdge(edgeToThrow[indexI]);
2303 }
2305 // Delete the two points...
2306 removePoint(cv2);
2307 removePoint(cv3);
2309 // Update map
2310 map.removePoint(cv2);
2311 map.removePoint(cv3);
2312 }
2314 if (debug > 2)
2315 {
2316 Pout<< nl
2317 << "~~~~~~~~~~~~~~~~~~~~~~~~" << nl
2318 << "Hulls after modification" << nl
2319 << "~~~~~~~~~~~~~~~~~~~~~~~~" << nl;
2321 Pout<< nl << "Cells belonging to first Edge Hull: "
2322 << hullCells[0] << nl;
2324 forAll(hullCells[0], cellI)
2325 {
2326 const cell& firstCurCell = cells_[hullCells[0][cellI]];
2328 Pout<< "Cell: " << hullCells[0][cellI]
2329 << ": " << firstCurCell << nl;
2331 forAll(firstCurCell, faceI)
2332 {
2333 Pout<< " Face: " << firstCurCell[faceI]
2334 << " : " << faces_[firstCurCell[faceI]]
2335 << " fE: " << faceEdges_[firstCurCell[faceI]]
2336 << nl;
2337 }
2338 }
2340 const labelList& firstEdgeFaces = edgeFaces_[checkEdgeIndex[1]];
2342 Pout<< nl << "First Edge Face Hull: " << firstEdgeFaces << nl;
2344 forAll(firstEdgeFaces, indexI)
2345 {
2346 Pout<< " Face: " << firstEdgeFaces[indexI]
2347 << " : " << faces_[firstEdgeFaces[indexI]]
2348 << " fE: " << faceEdges_[firstEdgeFaces[indexI]]
2349 << nl;
2350 }
2352 Pout<< nl << "Cells belonging to second Edge Hull: "
2353 << hullCells[1] << nl;
2355 forAll(hullCells[1], cellI)
2356 {
2357 const cell& secondCurCell = cells_[hullCells[1][cellI]];
2359 Pout<< "Cell: " << hullCells[1][cellI]
2360 << ": " << secondCurCell << nl;
2362 forAll(secondCurCell, faceI)
2363 {
2364 Pout<< " Face: " << secondCurCell[faceI]
2365 << " : " << faces_[secondCurCell[faceI]]
2366 << " fE: " << faceEdges_[secondCurCell[faceI]]
2367 << nl;
2368 }
2369 }
2371 const labelList& secondEdgeFaces = edgeFaces_[checkEdgeIndex[2]];
2373 Pout<< nl << "Second Edge Face Hull: " << secondEdgeFaces << nl;
2375 forAll(secondEdgeFaces, indexI)
2376 {
2377 Pout<< " Face : " << secondEdgeFaces[indexI]
2378 << " : " << faces_[secondEdgeFaces[indexI]]
2379 << " fE: " << faceEdges_[secondEdgeFaces[indexI]]
2380 << nl;
2381 }
2383 Pout<< "Retained face: "
2384 << faceToKeep[0] << ": "
2385 << " owner: " << owner_[faceToKeep[0]]
2386 << " neighbour: " << neighbour_[faceToKeep[0]]
2387 << nl;
2389 Pout<< "Discarded face: "
2390 << faceToThrow[0] << ": "
2391 << " owner: " << owner_[faceToThrow[0]]
2392 << " neighbour: " << neighbour_[faceToThrow[0]]
2393 << nl;
2395 if (c1 != -1)
2396 {
2397 Pout<< "Retained face: "
2398 << faceToKeep[1] << ": "
2399 << " owner: " << owner_[faceToKeep[1]]
2400 << " neighbour: " << neighbour_[faceToKeep[1]]
2401 << nl;
2403 Pout<< "Discarded face: "
2404 << faceToThrow[1] << ": "
2405 << " owner: " << owner_[faceToThrow[1]]
2406 << " neighbour: " << neighbour_[faceToThrow[1]]
2407 << nl;
2408 }
2410 Pout<< endl;
2411 }
2413 // Ensure proper orientation for the two retained faces
2414 FixedList<label,2> cellCheck(0);
2416 if (owner_[faceToThrow[0]] == c0)
2417 {
2418 cellCheck[0] = neighbour_[faceToThrow[0]];
2420 if (owner_[faceToKeep[0]] == c0)
2421 {
2422 if
2423 (
2424 (neighbour_[faceToThrow[0]] > neighbour_[faceToKeep[0]])
2425 && (neighbour_[faceToKeep[0]] != -1)
2426 )
2427 {
2428 // This face is to be flipped
2429 faces_[faceToKeep[0]] =
2430 (
2431 faces_[faceToKeep[0]].reverseFace()
2432 );
2434 owner_[faceToKeep[0]] = neighbour_[faceToKeep[0]];
2435 neighbour_[faceToKeep[0]] = neighbour_[faceToThrow[0]];
2437 setFlip(faceToKeep[0]);
2438 }
2439 else
2440 {
2441 if (neighbour_[faceToThrow[0]] != -1)
2442 {
2443 // Keep orientation intact, and update the owner
2444 owner_[faceToKeep[0]] = neighbour_[faceToThrow[0]];
2445 }
2446 else
2447 {
2448 // This face will need to be flipped and converted
2449 // to a boundary face. Flip it now, so that conversion
2450 // happens later.
2451 faces_[faceToKeep[0]] =
2452 (
2453 faces_[faceToKeep[0]].reverseFace()
2454 );
2456 owner_[faceToKeep[0]] = neighbour_[faceToKeep[0]];
2457 neighbour_[faceToKeep[0]] = -1;
2459 setFlip(faceToKeep[0]);
2460 }
2461 }
2462 }
2463 else
2464 {
2465 // Keep orientation intact, and update the neighbour
2466 neighbour_[faceToKeep[0]] = neighbour_[faceToThrow[0]];
2467 }
2468 }
2469 else
2470 {
2471 cellCheck[0] = owner_[faceToThrow[0]];
2473 if (neighbour_[faceToKeep[0]] == c0)
2474 {
2475 if (owner_[faceToThrow[0]] < owner_[faceToKeep[0]])
2476 {
2477 // This face is to be flipped
2478 faces_[faceToKeep[0]] =
2479 (
2480 faces_[faceToKeep[0]].reverseFace()
2481 );
2483 neighbour_[faceToKeep[0]] = owner_[faceToKeep[0]];
2484 owner_[faceToKeep[0]] = owner_[faceToThrow[0]];
2486 setFlip(faceToKeep[0]);
2487 }
2488 else
2489 {
2490 // Keep orientation intact, and update the neighbour
2491 neighbour_[faceToKeep[0]] = owner_[faceToThrow[0]];
2492 }
2493 }
2494 else
2495 {
2496 // Keep orientation intact, and update the owner
2497 owner_[faceToKeep[0]] = owner_[faceToThrow[0]];
2498 }
2499 }
2501 if (c1 != -1)
2502 {
2503 if (owner_[faceToThrow[1]] == c1)
2504 {
2505 cellCheck[1] = neighbour_[faceToThrow[1]];
2507 if (owner_[faceToKeep[1]] == c1)
2508 {
2509 if
2510 (
2511 (neighbour_[faceToThrow[1]] > neighbour_[faceToKeep[1]])
2512 && (neighbour_[faceToKeep[1]] != -1)
2513 )
2514 {
2515 // This face is to be flipped
2516 faces_[faceToKeep[1]] =
2517 (
2518 faces_[faceToKeep[1]].reverseFace()
2519 );
2521 owner_[faceToKeep[1]] = neighbour_[faceToKeep[1]];
2522 neighbour_[faceToKeep[1]] = neighbour_[faceToThrow[1]];
2524 setFlip(faceToKeep[1]);
2525 }
2526 else
2527 {
2528 if (neighbour_[faceToThrow[1]] != -1)
2529 {
2530 // Keep orientation intact, and update the owner
2531 owner_[faceToKeep[1]] = neighbour_[faceToThrow[1]];
2532 }
2533 else
2534 {
2535 // This face will need to be flipped and converted
2536 // to a boundary face. Flip it now, so that conversion
2537 // happens later.
2538 faces_[faceToKeep[1]] =
2539 (
2540 faces_[faceToKeep[1]].reverseFace()
2541 );
2543 owner_[faceToKeep[1]] = neighbour_[faceToKeep[1]];
2544 neighbour_[faceToKeep[1]] = -1;
2546 setFlip(faceToKeep[1]);
2547 }
2548 }
2549 }
2550 else
2551 {
2552 // Keep orientation intact, and update the neighbour
2553 neighbour_[faceToKeep[1]] = neighbour_[faceToThrow[1]];
2554 }
2555 }
2556 else
2557 {
2558 cellCheck[1] = owner_[faceToThrow[1]];
2560 if (neighbour_[faceToKeep[1]] == c1)
2561 {
2562 if (owner_[faceToThrow[1]] < owner_[faceToKeep[1]])
2563 {
2564 // This face is to be flipped
2565 faces_[faceToKeep[1]] =
2566 (
2567 faces_[faceToKeep[1]].reverseFace()
2568 );
2570 neighbour_[faceToKeep[1]] = owner_[faceToKeep[1]];
2571 owner_[faceToKeep[1]] = owner_[faceToThrow[1]];
2573 setFlip(faceToKeep[1]);
2574 }
2575 else
2576 {
2577 // Keep orientation intact, and update the neighbour
2578 neighbour_[faceToKeep[1]] = owner_[faceToThrow[1]];
2579 }
2580 }
2581 else
2582 {
2583 // Keep orientation intact, and update the owner
2584 owner_[faceToKeep[1]] = owner_[faceToThrow[1]];
2585 }
2586 }
2587 }
2589 // Remove orphaned faces
2590 if (owner_[faceToKeep[0]] == -1)
2591 {
2592 removeFace(faceToKeep[0]);
2594 // Update map
2595 map.removeFace(faceToKeep[0]);
2596 }
2597 else
2598 if
2599 (
2600 (neighbour_[faceToKeep[0]] == -1)
2601 && (whichPatch(faceToKeep[0]) < 0)
2602 )
2603 {
2604 // Obtain a copy before adding the new face,
2605 // since the reference might become invalid during list resizing.
2606 face newFace = faces_[faceToKeep[0]];
2607 label newOwn = owner_[faceToKeep[0]];
2608 labelList newFaceEdges = faceEdges_[faceToKeep[0]];
2610 // This face is being converted from interior to boundary. Remove
2611 // from the interior list and add as a boundary face to the end.
2612 label newFaceIndex =
2613 (
2614 insertFace
2615 (
2616 whichPatch(faceToThrow[0]),
2617 newFace,
2618 newOwn,
2619 -1
2620 )
2621 );
2623 // Add an entry for mapping
2624 modifiedFaces.append(newFaceIndex);
2626 // Update map
2627 map.addFace(newFaceIndex, labelList(1, faceToThrow[0]));
2629 // Add a faceEdges entry as well.
2630 // Edges don't have to change, since they're
2631 // all on the boundary anyway.
2632 faceEdges_.append(newFaceEdges);
2634 meshOps::replaceLabel
2635 (
2636 faceToKeep[0],
2637 newFaceIndex,
2638 cells_[newOwn]
2639 );
2641 // Correct edgeFaces with the new face label.
2642 forAll(newFaceEdges, edgeI)
2643 {
2644 meshOps::replaceLabel
2645 (
2646 faceToKeep[0],
2647 newFaceIndex,
2648 edgeFaces_[newFaceEdges[edgeI]]
2649 );
2650 }
2652 // Renumber the neighbour so that this face is removed correctly.
2653 neighbour_[faceToKeep[0]] = 0;
2654 removeFace(faceToKeep[0]);
2656 // Update map
2657 map.removeFace(faceToKeep[0]);
2658 }
2660 // Remove the unwanted faces in the cell(s) adjacent to this face,
2661 // and correct the cells that contain discarded faces
2662 const cell& cell_0 = cells_[c0];
2664 forAll(cell_0,faceI)
2665 {
2666 if (cell_0[faceI] != fIndex && cell_0[faceI] != faceToKeep[0])
2667 {
2668 removeFace(cell_0[faceI]);
2670 // Update map
2671 map.removeFace(cell_0[faceI]);
2672 }
2673 }
2675 if (cellCheck[0] != -1)
2676 {
2677 meshOps::replaceLabel
2678 (
2679 faceToThrow[0],
2680 faceToKeep[0],
2681 cells_[cellCheck[0]]
2682 );
2683 }
2685 // Remove the cell
2686 removeCell(c0);
2688 // Update map
2689 map.removeCell(c0);
2691 if (c1 == -1)
2692 {
2693 // Increment the surface-collapse counter
2694 statistics_[6]++;
2695 }
2696 else
2697 {
2698 // Remove orphaned faces
2699 if (owner_[faceToKeep[1]] == -1)
2700 {
2701 removeFace(faceToKeep[1]);
2703 // Update map
2704 map.removeFace(faceToKeep[1]);
2705 }
2706 else
2707 if
2708 (
2709 (neighbour_[faceToKeep[1]] == -1)
2710 && (whichPatch(faceToKeep[1]) < 0)
2711 )
2712 {
2713 // Obtain a copy before adding the new face,
2714 // since the reference might become invalid during list resizing.
2715 face newFace = faces_[faceToKeep[1]];
2716 label newOwn = owner_[faceToKeep[1]];
2717 labelList newFaceEdges = faceEdges_[faceToKeep[1]];
2719 // This face is being converted from interior to boundary. Remove
2720 // from the interior list and add as a boundary face to the end.
2721 label newFaceIndex =
2722 (
2723 insertFace
2724 (
2725 whichPatch(faceToThrow[1]),
2726 newFace,
2727 newOwn,
2728 -1
2729 )
2730 );
2732 // Add an entry for mapping
2733 modifiedFaces.append(newFaceIndex);
2735 // Update map
2736 map.addFace(newFaceIndex, labelList(1, faceToThrow[1]));
2738 // Add a faceEdges entry as well.
2739 // Edges don't have to change, since they're
2740 // all on the boundary anyway.
2741 faceEdges_.append(newFaceEdges);
2743 meshOps::replaceLabel
2744 (
2745 faceToKeep[1],
2746 newFaceIndex,
2747 cells_[newOwn]
2748 );
2750 // Correct edgeFaces with the new face label.
2751 forAll(newFaceEdges, edgeI)
2752 {
2753 meshOps::replaceLabel
2754 (
2755 faceToKeep[1],
2756 newFaceIndex,
2757 edgeFaces_[newFaceEdges[edgeI]]
2758 );
2759 }
2761 // Renumber the neighbour so that this face is removed correctly.
2762 neighbour_[faceToKeep[1]] = 0;
2763 removeFace(faceToKeep[1]);
2765 // Update map
2766 map.removeFace(faceToKeep[1]);
2767 }
2769 const cell& cell_1 = cells_[c1];
2771 forAll(cell_1, faceI)
2772 {
2773 if (cell_1[faceI] != fIndex && cell_1[faceI] != faceToKeep[1])
2774 {
2775 removeFace(cell_1[faceI]);
2777 // Update map
2778 map.removeFace(cell_1[faceI]);
2779 }
2780 }
2782 if (cellCheck[1] != -1)
2783 {
2784 meshOps::replaceLabel
2785 (
2786 faceToThrow[1],
2787 faceToKeep[1],
2788 cells_[cellCheck[1]]
2789 );
2790 }
2792 // Remove the cell
2793 removeCell(c1);
2795 // Update map
2796 map.removeCell(c1);
2797 }
2799 // Move old / new points
2800 oldPoints_[replacement[0]] = oldPoint[0];
2801 oldPoints_[replacement[1]] = oldPoint[1];
2803 points_[replacement[0]] = newPoint[0];
2804 points_[replacement[1]] = newPoint[1];
2806 // Finally remove the face
2807 removeFace(fIndex);
2809 // Update map
2810 map.removeFace(fIndex);
2812 // Write out VTK files after change
2813 if (debug > 3)
2814 {
2815 labelHashSet vtkCells;
2817 forAll(hullCells[0], cellI)
2818 {
2819 if (hullCells[0][cellI] == c0 || hullCells[0][cellI] == c1)
2820 {
2821 continue;
2822 }
2824 if (!vtkCells.found(hullCells[0][cellI]))
2825 {
2826 vtkCells.insert(hullCells[0][cellI]);
2827 }
2828 }
2830 forAll(hullCells[1], cellI)
2831 {
2832 if (hullCells[1][cellI] == c0 || hullCells[1][cellI] == c1)
2833 {
2834 continue;
2835 }
2837 if (!vtkCells.found(hullCells[1][cellI]))
2838 {
2839 vtkCells.insert(hullCells[1][cellI]);
2840 }
2841 }
2843 writeVTK
2844 (
2845 Foam::name(fIndex)
2846 + "_Collapse_1",
2847 vtkCells.toc()
2848 );
2849 }
2851 // Fill-in candidate mapping information
2852 labelList mC(2, -1);
2853 mC[0] = c0, mC[1] = c1;
2855 // Now that all old / new cells possess correct connectivity,
2856 // compute mapping information.
2857 forAll(hullCells, indexI)
2858 {
2859 forAll(hullCells[indexI], cellI)
2860 {
2861 label mcIndex = hullCells[indexI][cellI];
2863 // Skip collapsed cells
2864 if (mcIndex == c0 || mcIndex == c1)
2865 {
2866 continue;
2867 }
2869 // Set the mapping for this cell
2870 setCellMapping(mcIndex, mC);
2871 }
2872 }
2874 // Set face mapping information for modified faces
2875 forAll(modifiedFaces, faceI)
2876 {
2877 const label mfIndex = modifiedFaces[faceI];
2879 // Exclude deleted faces
2880 if (faces_[mfIndex].empty())
2881 {
2882 continue;
2883 }
2885 // Decide between default / weighted mapping
2886 // based on boundary information
2887 label fPatch = whichPatch(mfIndex);
2889 if (fPatch == -1)
2890 {
2891 setFaceMapping(mfIndex);
2892 }
2893 else
2894 {
2895 // Fill-in candidate mapping information
2896 labelList faceCandidates;
2898 const labelList& fEdges = faceEdges_[mfIndex];
2900 forAll(fEdges, edgeI)
2901 {
2902 if (whichEdgePatch(fEdges[edgeI]) == fPatch)
2903 {
2904 // Loop through associated edgeFaces
2905 const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
2907 forAll(eFaces, faceI)
2908 {
2909 if
2910 (
2911 (eFaces[faceI] != mfIndex) &&
2912 (whichPatch(eFaces[faceI]) == fPatch)
2913 )
2914 {
2915 faceCandidates.setSize
2916 (
2917 faceCandidates.size() + 1,
2918 eFaces[faceI]
2919 );
2920 }
2921 }
2922 }
2923 }
2925 // Set the mapping for this face
2926 setFaceMapping(mfIndex, faceCandidates);
2927 }
2928 }
2930 // If modification is coupled, update mapping info.
2931 if (coupledModification_)
2932 {
2933 // Check if the collapse points are present
2934 // on a processor not involved in the current
2935 // operation, and update if necessary
2936 if (procCouple && !localCouple)
2937 {
2938 forAll(procIndices_, pI)
2939 {
2940 bool involved = false;
2942 forAll(slaveMaps, slaveI)
2943 {
2944 // Alias for convenience...
2945 const changeMap& slaveMap = slaveMaps[slaveI];
2947 if (slaveMap.patchIndex() == pI && slaveMap.index() >= 0)
2948 {
2949 // Involved in this operation. Break out.
2950 involved = true;
2951 break;
2952 }
2953 }
2955 if (involved)
2956 {
2957 continue;
2958 }
2960 // Check coupleMaps for point coupling
2961 const label pointEnum = coupleMap::POINT;
2963 const coupledInfo& recvMesh = recvMeshes_[pI];
2964 const coupleMap& cMap = recvMesh.map();
2966 // Obtain non-const references
2967 Map<label>& pointMap = cMap.entityMap(pointEnum);
2968 Map<label>& rPointMap = cMap.reverseEntityMap(pointEnum);
2970 label sI0 = -1, sI1 = -1;
2972 if (collapsingSlave)
2973 {
2974 if ((sI0 = cMap.findMaster(pointEnum, original[0])) > -1)
2975 {
2976 if (rPointMap.found(replacement[0]))
2977 {
2978 rPointMap[replacement[0]] = sI0;
2979 }
2980 else
2981 {
2982 rPointMap.insert(replacement[0], sI0);
2983 }
2985 pointMap[sI0] = replacement[0];
2986 }
2988 if ((sI1 = cMap.findMaster(pointEnum, original[1])) > -1)
2989 {
2990 if (rPointMap.found(replacement[1]))
2991 {
2992 rPointMap[replacement[1]] = sI1;
2993 }
2994 else
2995 {
2996 rPointMap.insert(replacement[1], sI1);
2997 }
2999 pointMap[sI1] = replacement[1];
3000 }
3001 }
3002 else
3003 {
3004 if ((sI0 = cMap.findSlave(pointEnum, original[0])) > -1)
3005 {
3006 if (pointMap.found(replacement[0]))
3007 {
3008 pointMap[replacement[0]] = sI0;
3009 }
3010 else
3011 {
3012 pointMap.insert(replacement[0], sI0);
3013 }
3015 rPointMap[sI0] = replacement[0];
3016 }
3018 if ((sI1 = cMap.findSlave(pointEnum, original[1])) > -1)
3019 {
3020 if (pointMap.found(replacement[1]))
3021 {
3022 pointMap[replacement[1]] = sI1;
3023 }
3024 else
3025 {
3026 pointMap.insert(replacement[1], sI1);
3027 }
3029 rPointMap[sI1] = replacement[1];
3030 }
3031 }
3033 if (sI0 > -1 && sI1 > -1 && debug > 2)
3034 {
3035 Pout<< " Found " << original[0] << " and " << original[1]
3036 << " on proc: " << procIndices_[pI]
3037 << endl;
3038 }
3039 }
3040 }
3042 forAll(slaveMaps, slaveI)
3043 {
3044 // Alias for convenience...
3045 const changeMap& slaveMap = slaveMaps[slaveI];
3047 // Skip updates for edge-based coupling
3048 if (slaveMap.index() < 0)
3049 {
3050 continue;
3051 }
3053 label pI = slaveMap.patchIndex();
3055 // Fetch the appropriate coupleMap
3056 const coupleMap* cMapPtr = NULL;
3058 if (localCouple && !procCouple)
3059 {
3060 cMapPtr = &(patchCoupling_[pI].map());
3061 }
3062 else
3063 if (procCouple && !localCouple)
3064 {
3065 cMapPtr = &(recvMeshes_[pI].map());
3066 }
3068 // Configure the slave replacement points.
3069 // - collapseQuadFace stores this as 'addedPoints'
3070 label scP0 = slaveMap.removedPointList()[0];
3071 label scP1 = slaveMap.removedPointList()[1];
3073 label srP0 = slaveMap.addedPointList()[0].index();
3074 label srP1 = slaveMap.addedPointList()[1].index();
3076 // Alias for convenience
3077 const coupleMap& cMap = *cMapPtr;
3079 // Remove the point entries.
3080 const label pointEnum = coupleMap::POINT;
3082 // Obtain references
3083 Map<label>& pointMap = cMap.entityMap(pointEnum);
3084 Map<label>& rPointMap = cMap.reverseEntityMap(pointEnum);
3086 if (collapsingSlave)
3087 {
3088 if (rPointMap[replacement[0]] == scP0)
3089 {
3090 pointMap[srP0] = replacement[0];
3091 rPointMap[replacement[0]] = srP0;
3092 }
3093 else
3094 if (rPointMap[replacement[0]] == scP1)
3095 {
3096 pointMap[srP1] = replacement[0];
3097 rPointMap[replacement[0]] = srP1;
3098 }
3100 pointMap.erase(rPointMap[original[0]]);
3101 rPointMap.erase(original[0]);
3102 }
3103 else
3104 {
3105 if (pointMap[replacement[0]] == scP0)
3106 {
3107 rPointMap[srP0] = replacement[0];
3108 pointMap[replacement[0]] = srP0;
3109 }
3110 else
3111 if (pointMap[replacement[0]] == scP1)
3112 {
3113 rPointMap[srP1] = replacement[0];
3114 pointMap[replacement[0]] = srP1;
3115 }
3117 rPointMap.erase(pointMap[original[0]]);
3118 pointMap.erase(original[0]);
3119 }
3121 if (collapsingSlave)
3122 {
3123 if (rPointMap[replacement[1]] == scP0)
3124 {
3125 pointMap[srP0] = replacement[1];
3126 rPointMap[replacement[1]] = srP0;
3127 }
3128 else
3129 if (rPointMap[replacement[1]] == scP1)
3130 {
3131 pointMap[srP1] = replacement[1];
3132 rPointMap[replacement[1]] = srP1;
3133 }
3135 pointMap.erase(rPointMap[original[1]]);
3136 rPointMap.erase(original[1]);
3137 }
3138 else
3139 {
3140 if (pointMap[replacement[1]] == scP0)
3141 {
3142 rPointMap[srP0] = replacement[1];
3143 pointMap[replacement[1]] = srP0;
3144 }
3145 else
3146 if (pointMap[replacement[1]] == scP1)
3147 {
3148 rPointMap[srP1] = replacement[1];
3149 pointMap[replacement[1]] = srP1;
3150 }
3152 rPointMap.erase(pointMap[original[1]]);
3153 pointMap.erase(original[1]);
3154 }
3156 // Remove the face entries
3157 const label faceEnum = coupleMap::FACE;
3159 // Obtain references
3160 Map<label>& faceMap = cMap.entityMap(faceEnum);
3161 Map<label>& rFaceMap = cMap.reverseEntityMap(faceEnum);
3163 if (collapsingSlave)
3164 {
3165 faceMap.erase(faceMap[fIndex]);
3166 rFaceMap.erase(fIndex);
3167 }
3168 else
3169 {
3170 rFaceMap.erase(faceMap[fIndex]);
3171 faceMap.erase(fIndex);
3172 }
3174 // Configure a comparison face
3175 face cFace(4);
3177 // If any interior faces in the master map were
3178 // converted to boundaries, account for it
3179 const List<objectMap>& madF = map.addedFaceList();
3181 forAll(madF, faceI)
3182 {
3183 label mfIndex = madF[faceI].index();
3184 const face& mFace = faces_[mfIndex];
3186 // Select appropriate mesh
3187 const dynamicTopoFvMesh* meshPtr = NULL;
3189 // Fetch the appropriate coupleMap
3190 const coupleMap* crMapPtr = NULL;
3192 // Fetch patch info
3193 label ofPatch = whichPatch(fIndex);
3194 label mfPatch = whichPatch(mfIndex);
3196 if (localCouple && !procCouple)
3197 {
3198 // Local coupling. Use this mesh itself
3199 meshPtr = this;
3200 crMapPtr = &(patchCoupling_[pI].map());
3201 }
3202 else
3203 if (procCouple && !localCouple)
3204 {
3205 // Occasionally, this face might talk to
3206 // a processor other than the slave
3207 if (ofPatch == mfPatch)
3208 {
3209 meshPtr = &(recvMeshes_[pI].subMesh());
3210 crMapPtr = &(recvMeshes_[pI].map());
3211 }
3212 else
3213 {
3214 label neiProcNo = getNeighbourProcessor(mfPatch);
3216 if (neiProcNo == -1)
3217 {
3218 // Not a processor patch. No mapping required.
3219 continue;
3220 }
3221 else
3222 {
3223 // Find an appropriate subMesh
3224 label prI = findIndex(procIndices_, neiProcNo);
3226 meshPtr = &(recvMeshes_[prI].subMesh());
3227 crMapPtr = &(recvMeshes_[prI].map());
3228 }
3229 }
3230 }
3232 bool matchedFace = false;
3234 // Alias for convenience
3235 const coupleMap& crMap = *crMapPtr;
3236 const dynamicTopoFvMesh& sMesh = *meshPtr;
3238 // Obtain references
3239 Map<label>& pMap = crMap.entityMap(pointEnum);
3240 Map<label>& fMap = crMap.entityMap(faceEnum);
3241 Map<label>& rfMap = crMap.reverseEntityMap(faceEnum);
3243 // Configure the face
3244 forAll(cFace, pointI)
3245 {
3246 cFace[pointI] =
3247 (
3248 pMap.found(mFace[pointI]) ?
3249 pMap[mFace[pointI]] : -1
3250 );
3251 }
3253 // Loop through all boundary faces on the subMesh
3254 for
3255 (
3256 label faceJ = sMesh.nOldInternalFaces_;
3257 faceJ < sMesh.faces_.size();
3258 faceJ++
3259 )
3260 {
3261 const face& sFace = sMesh.faces_[faceJ];
3263 if (face::compare(sFace, cFace))
3264 {
3265 if (debug > 2)
3266 {
3267 Pout<< " Found face: " << faceJ
3268 << "::" << sFace
3269 << " with mfIndex: " << mfIndex
3270 << "::" << mFace
3271 << endl;
3272 }
3274 // Update faceMaps
3275 if (collapsingSlave)
3276 {
3277 if (fMap.found(faceJ))
3278 {
3279 fMap[faceJ] = mfIndex;
3280 }
3281 else
3282 {
3283 fMap.insert(faceJ, mfIndex);
3284 }
3286 if (rfMap.found(mfIndex))
3287 {
3288 rfMap[mfIndex] = faceJ;
3289 }
3290 else
3291 {
3292 rfMap.insert(mfIndex, faceJ);
3293 }
3294 }
3295 else
3296 {
3297 if (rfMap.found(faceJ))
3298 {
3299 rfMap[faceJ] = mfIndex;
3300 }
3301 else
3302 {
3303 rfMap.insert(faceJ, mfIndex);
3304 }
3306 if (fMap.found(mfIndex))
3307 {
3308 fMap[mfIndex] = faceJ;
3309 }
3310 else
3311 {
3312 fMap.insert(mfIndex, faceJ);
3313 }
3314 }
3316 matchedFace = true;
3318 break;
3319 }
3320 }
3322 if (!matchedFace)
3323 {
3324 // Write out for post-processing
3325 writeVTK("masterFace_" + Foam::name(mfIndex), mfIndex, 2);
3327 FatalErrorIn
3328 (
3329 "\n"
3330 "const changeMap "
3331 "dynamicTopoFvMesh::collapseQuadFace\n"
3332 "(\n"
3333 " const label fIndex,\n"
3334 " label overRideCase,\n"
3335 " bool checkOnly,\n"
3336 " bool forceOp\n"
3337 ")\n"
3338 )
3339 << " Master face: " << mfIndex
3340 << ": " << mFace << " could not be matched." << nl
3341 << " cFace: " << cFace << nl
3342 << abort(FatalError);
3343 }
3344 }
3346 // If any interior faces in the slave map were
3347 // converted to boundaries, account for it
3348 const List<objectMap>& sadF = slaveMap.addedFaceList();
3350 forAll(sadF, faceI)
3351 {
3352 label sIndex = slaveMap.index();
3353 label sfIndex = sadF[faceI].index();
3355 // Select appropriate mesh
3356 const dynamicTopoFvMesh* meshPtr = NULL;
3358 if (localCouple && !procCouple)
3359 {
3360 // Local coupling. Use this mesh itself
3361 meshPtr = this;
3362 }
3363 else
3364 if (procCouple && !localCouple)
3365 {
3366 meshPtr = &(recvMeshes_[pI].subMesh());
3367 }
3369 // Alias for convenience
3370 const dynamicTopoFvMesh& sMesh = *meshPtr;
3372 label ofPatch = sMesh.whichPatch(sIndex);
3373 label sfPatch = sMesh.whichPatch(sfIndex);
3375 // Skip dissimilar patches
3376 if (ofPatch != sfPatch)
3377 {
3378 continue;
3379 }
3381 const face& sFace = sMesh.faces_[sfIndex];
3383 forAll(cFace, pointI)
3384 {
3385 cFace[pointI] =
3386 (
3387 rPointMap.found(sFace[pointI]) ?
3388 rPointMap[sFace[pointI]] : -1
3389 );
3390 }
3392 bool matchedFace = false;
3394 // Loop through all boundary faces on this mesh
3395 for
3396 (
3397 label faceJ = nOldInternalFaces_;
3398 faceJ < faces_.size();
3399 faceJ++
3400 )
3401 {
3402 const face& mFace = faces_[faceJ];
3404 if (face::compare(mFace, cFace))
3405 {
3406 if (debug > 2)
3407 {
3408 Pout<< " Found face: " << faceJ
3409 << "::" << mFace
3410 << " with sfIndex: " << sfIndex
3411 << "::" << sFace
3412 << endl;
3413 }
3415 // Update faceMaps
3416 if (collapsingSlave)
3417 {
3418 if (rFaceMap.found(faceJ))
3419 {
3420 rFaceMap[faceJ] = sfIndex;
3421 }
3422 else
3423 {
3424 rFaceMap.insert(faceJ, sfIndex);
3425 }
3427 if (faceMap.found(sfIndex))
3428 {
3429 faceMap[sfIndex] = faceJ;
3430 }
3431 else
3432 {
3433 faceMap.insert(sfIndex, faceJ);
3434 }
3435 }
3436 else
3437 {
3438 if (faceMap.found(faceJ))
3439 {
3440 faceMap[faceJ] = sfIndex;
3441 }
3442 else
3443 {
3444 faceMap.insert(faceJ, sfIndex);
3445 }
3447 if (rFaceMap.found(sfIndex))
3448 {
3449 rFaceMap[sfIndex] = faceJ;
3450 }
3451 else
3452 {
3453 rFaceMap.insert(sfIndex, faceJ);
3454 }
3455 }
3457 matchedFace = true;
3459 break;
3460 }
3461 }
3463 if (!matchedFace)
3464 {
3465 FatalErrorIn
3466 (
3467 "\n"
3468 "const changeMap "
3469 "dynamicTopoFvMesh::collapseQuadFace\n"
3470 "(\n"
3471 " const label fIndex,\n"
3472 " label overRideCase,\n"
3473 " bool checkOnly,\n"
3474 " bool forceOp\n"
3475 ")\n"
3476 )
3477 << " Slave face: " << sfIndex
3478 << ": " << sFace << " could not be matched." << nl
3479 << " cFace: " << cFace << nl
3480 << abort(FatalError);
3481 }
3482 }
3484 // Push operation into coupleMap
3485 switch (slaveMap.type())
3486 {
3487 case 1:
3488 {
3489 cMap.pushOperation
3490 (
3491 slaveMap.index(),
3492 coupleMap::COLLAPSE_FIRST
3493 );
3495 break;
3496 }
3498 case 2:
3499 {
3500 cMap.pushOperation
3501 (
3502 slaveMap.index(),
3503 coupleMap::COLLAPSE_SECOND
3504 );
3506 break;
3507 }
3509 case 3:
3510 {
3511 cMap.pushOperation
3512 (
3513 slaveMap.index(),
3514 coupleMap::COLLAPSE_MIDPOINT
3515 );
3517 break;
3518 }
3519 }
3520 }
3521 }
3523 // Set the flag
3524 topoChangeFlag_ = true;
3526 // Increment the counter
3527 statistics_[4]++;
3529 // Increment the number of modifications
3530 statistics_[0]++;
3532 // Return a succesful collapse
3533 map.type() = collapseCase;
3535 return map;
3536 }
3539 // Method for the collapse of an edge in 3D
3540 // - Returns a changeMap with a type specifying:
3541 // -3: Collapse failed since edge was on a noRefinement patch.
3542 // -1: Collapse failed since max number of topo-changes was reached.
3543 // 0: Collapse could not be performed.
3544 // 1: Collapsed to first node.
3545 // 2: Collapsed to second node.
3546 // 3: Collapsed to mid-point (default).
3547 // - overRideCase is used to force a certain collapse configuration.
3548 // -1: Use this value to let collapseEdge decide a case.
3549 // 1: Force collapse to first node.
3550 // 2: Force collapse to second node.
3551 // 3: Force collapse to mid-point.
3552 // - checkOnly performs a feasibility check and returns without modifications.
3553 // - forceOp to force the collapse.
3554 const changeMap dynamicTopoFvMesh::collapseEdge
3555 (
3556 const label eIndex,
3557 label overRideCase,
3558 bool checkOnly,
3559 bool forceOp
3560 )
3561 {
3562 // Edge collapse performs the following operations:
3563 // [1] Checks if either vertex of the edge is on a boundary
3564 // [2] Checks whether cells attached to deleted vertices will be valid
3565 // after the edge-collapse operation
3566 // [3] Deletes all cells surrounding this edge
3567 // [4] Deletes all faces surrounding this edge
3568 // [5] Deletes all faces surrounding the deleted vertex attached
3569 // to the cells in [3]
3570 // [6] Checks the orientation of faces connected to the retained
3571 // vertices
3572 // [7] Remove one of the vertices of the edge
3573 // Update faceEdges and edgeFaces information
3575 // For 2D meshes, perform face-collapse
3576 if (twoDMesh_)
3577 {
3578 return collapseQuadFace(eIndex, overRideCase, checkOnly);
3579 }
3581 // Figure out which thread this is...
3582 label tIndex = self();
3584 // Prepare the changeMaps
3585 changeMap map;
3586 List<changeMap> slaveMaps;
3587 bool collapsingSlave = false;
3589 if
3590 (
3591 (statistics_[0] > maxModifications_)
3592 && (maxModifications_ > -1)
3593 )
3594 {
3595 // Reached the max allowable topo-changes.
3596 stack(tIndex).clear();
3598 return map;
3599 }
3601 // Check if edgeRefinements are to be avoided on patch.
3602 const labelList& eF = edgeFaces_[eIndex];
3604 forAll(eF, fI)
3605 {
3606 label fPatch = whichPatch(eF[fI]);
3608 if (baseMesh_.lengthEstimator().checkRefinementPatch(fPatch))
3609 {
3610 map.type() = -3;
3612 return map;
3613 }
3614 }
3616 // Sanity check: Is the index legitimate?
3617 if (eIndex < 0)
3618 {
3619 FatalErrorIn
3620 (
3621 "\n"
3622 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
3623 "(\n"
3624 " const label eIndex,\n"
3625 " label overRideCase,\n"
3626 " bool checkOnly,\n"
3627 " bool forceOp\n"
3628 ")\n"
3629 )
3630 << " Invalid index: " << eIndex
3631 << abort(FatalError);
3632 }
3634 // If coupled modification is set, and this is a
3635 // master edge, collapse its slaves as well.
3636 bool localCouple = false, procCouple = false;
3638 if (coupledModification_)
3639 {
3640 const edge& eCheck = edges_[eIndex];
3642 const label edgeEnum = coupleMap::EDGE;
3643 const label pointEnum = coupleMap::POINT;
3645 // Is this a locally coupled edge (either master or slave)?
3646 if (locallyCoupledEntity(eIndex, true))
3647 {
3648 localCouple = true;
3649 procCouple = false;
3650 }
3651 else
3652 if (processorCoupledEntity(eIndex))
3653 {
3654 procCouple = true;
3655 localCouple = false;
3656 }
3658 if (localCouple && !procCouple)
3659 {
3660 // Determine the slave index.
3661 label sIndex = -1, pIndex = -1;
3663 forAll(patchCoupling_, patchI)
3664 {
3665 if (patchCoupling_(patchI))
3666 {
3667 const coupleMap& cMap = patchCoupling_[patchI].map();
3669 if ((sIndex = cMap.findSlave(edgeEnum, eIndex)) > -1)
3670 {
3671 pIndex = patchI;
3673 break;
3674 }
3676 // The following bit happens only during the sliver
3677 // exudation process, since slave edges are
3678 // usually not added to the coupled edge-stack.
3679 if ((sIndex = cMap.findMaster(edgeEnum, eIndex)) > -1)
3680 {
3681 pIndex = patchI;
3683 // Notice that we are collapsing a slave edge first.
3684 collapsingSlave = true;
3686 break;
3687 }
3688 }
3689 }
3691 if (sIndex == -1)
3692 {
3693 FatalErrorIn
3694 (
3695 "\n"
3696 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
3697 "(\n"
3698 " const label eIndex,\n"
3699 " label overRideCase,\n"
3700 " bool checkOnly,\n"
3701 " bool forceOp\n"
3702 ")\n"
3703 )
3704 << "Coupled maps were improperly specified." << nl
3705 << " Slave index not found for: " << nl
3706 << " Edge: " << eIndex << ": " << eCheck << nl
3707 << abort(FatalError);
3708 }
3709 else
3710 {
3711 // If we've found the slave, size up the list
3712 meshOps::sizeUpList
3713 (
3714 changeMap(),
3715 slaveMaps
3716 );
3718 // Save index and patch for posterity
3719 slaveMaps[0].index() = sIndex;
3720 slaveMaps[0].patchIndex() = pIndex;
3721 }
3723 if (debug > 1)
3724 {
3725 Pout<< nl << " >> Collapsing slave edge: " << sIndex
3726 << " for master edge: " << eIndex << endl;
3727 }
3728 }
3729 else
3730 if (procCouple && !localCouple)
3731 {
3732 // If this is a new entity, bail out for now.
3733 // This will be handled at the next time-step.
3734 if (eIndex >= nOldEdges_)
3735 {
3736 return map;
3737 }
3739 // Check slaves
3740 forAll(procIndices_, pI)
3741 {
3742 // Fetch reference to subMeshes
3743 const coupledInfo& sendMesh = sendMeshes_[pI];
3744 const coupledInfo& recvMesh = recvMeshes_[pI];
3746 const coupleMap& scMap = sendMesh.map();
3747 const coupleMap& rcMap = recvMesh.map();
3749 // If this edge was sent to a lower-ranked
3750 // processor, skip it.
3751 if (procIndices_[pI] < Pstream::myProcNo())
3752 {
3753 if (scMap.reverseEntityMap(edgeEnum).found(eIndex))
3754 {
3755 if (debug > 3)
3756 {
3757 Pout<< "Edge: " << eIndex
3758 << "::" << eCheck
3759 << " was sent to proc: "
3760 << procIndices_[pI]
3761 << ", so bailing out."
3762 << endl;
3763 }
3765 return map;
3766 }
3767 }
3769 label sIndex = -1;
3771 if ((sIndex = rcMap.findSlave(edgeEnum, eIndex)) > -1)
3772 {
3773 // Check if a lower-ranked processor is
3774 // handling this edge
3775 if (procIndices_[pI] < Pstream::myProcNo())
3776 {
3777 if (debug > 3)
3778 {
3779 Pout<< "Edge: " << eIndex
3780 << "::" << eCheck
3781 << " is handled by proc: "
3782 << procIndices_[pI]
3783 << ", so bailing out."
3784 << endl;
3785 }
3787 return map;
3788 }
3790 label curIndex = slaveMaps.size();
3792 // Size up the list
3793 meshOps::sizeUpList
3794 (
3795 changeMap(),
3796 slaveMaps
3797 );
3799 // Save index and patch for posterity
3800 slaveMaps[curIndex].index() = sIndex;
3801 slaveMaps[curIndex].patchIndex() = pI;
3802 }
3803 else
3804 if
3805 (
3806 (rcMap.findSlave(pointEnum, eCheck[0]) > -1) ||
3807 (rcMap.findSlave(pointEnum, eCheck[1]) > -1)
3808 )
3809 {
3810 // A point-only coupling exists.
3812 // Check if a lower-ranked processor is
3813 // handling this edge
3814 if (procIndices_[pI] < Pstream::myProcNo())
3815 {
3816 if (debug > 3)
3817 {
3818 Pout<< "Edge point on: " << eIndex
3819 << "::" << eCheck
3820 << " is handled by proc: "
3821 << procIndices_[pI]
3822 << ", so bailing out."
3823 << endl;
3824 }
3826 return map;
3827 }
3829 label p0Index = rcMap.findSlave(pointEnum, eCheck[0]);
3830 label p1Index = rcMap.findSlave(pointEnum, eCheck[1]);
3832 if (p0Index > -1 && p1Index == -1)
3833 {
3834 sIndex = p0Index;
3835 }
3836 else
3837 if (p0Index == -1 && p1Index > -1)
3838 {
3839 sIndex = p1Index;
3840 }
3842 label curIndex = slaveMaps.size();
3844 // Size up the list
3845 meshOps::sizeUpList
3846 (
3847 changeMap(),
3848 slaveMaps
3849 );
3851 // Save index and patch for posterity
3852 // - Negate the index to signify point coupling
3853 slaveMaps[curIndex].index() = -sIndex;
3854 slaveMaps[curIndex].patchIndex() = pI;
3855 }
3856 }
3857 }
3858 else
3859 {
3860 // Something's wrong with coupling maps
3861 FatalErrorIn
3862 (
3863 "\n"
3864 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
3865 "(\n"
3866 " const label eIndex,\n"
3867 " label overRideCase,\n"
3868 " bool checkOnly,\n"
3869 " bool forceOp\n"
3870 ")\n"
3871 )
3872 << "Coupled maps were improperly specified." << nl
3873 << " localCouple: " << localCouple << nl
3874 << " procCouple: " << procCouple << nl
3875 << " Edge: " << eIndex << ": " << eCheck << nl
3876 << abort(FatalError);
3877 }
3879 // Temporarily turn off coupledModification
3880 unsetCoupledModification();
3882 // Test the master edge for collapse, and decide on a case
3883 changeMap masterMap = collapseEdge(eIndex, -1, true, forceOp);
3885 // Turn it back on.
3886 setCoupledModification();
3888 // Master couldn't perform collapse
3889 if (masterMap.type() <= 0)
3890 {
3891 return masterMap;
3892 }
3894 // For point-only coupling, define the points for checking
3895 pointField slaveMoveNewPoint(slaveMaps.size(), vector::zero);
3896 pointField slaveMoveOldPoint(slaveMaps.size(), vector::zero);
3898 // Now check each of the slaves for collapse feasibility
3899 forAll(slaveMaps, slaveI)
3900 {
3901 // Alias for convenience...
3902 changeMap& slaveMap = slaveMaps[slaveI];
3904 label slaveOverRide = -1;
3905 label sIndex = slaveMap.index();
3906 label pI = slaveMap.patchIndex();
3908 // If the edge is mapped onto itself, skip check
3909 // (occurs for cyclic edges)
3910 if ((sIndex == eIndex) && localCouple)
3911 {
3912 continue;
3913 }
3915 const coupleMap* cMapPtr = NULL;
3917 edge mEdge(eCheck), sEdge(-1, -1);
3919 if (localCouple)
3920 {
3921 sEdge = edges_[sIndex];
3923 cMapPtr = &(patchCoupling_[pI].map());
3924 }
3925 else
3926 if (procCouple)
3927 {
3928 const dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
3930 cMapPtr = &(recvMeshes_[pI].map());
3932 if (sIndex < 0)
3933 {
3934 if (debug > 3)
3935 {
3936 Pout<< "Checking slave point: " << mag(sIndex)
3937 << "::" << sMesh.points_[mag(sIndex)]
3938 << " on proc: " << procIndices_[pI]
3939 << " for master edge: " << eIndex
3940 << " using collapseCase: " << masterMap.type()
3941 << endl;
3942 }
3943 }
3944 else
3945 {
3946 sEdge = sMesh.edges_[sIndex];
3948 if (debug > 3)
3949 {
3950 Pout<< "Checking slave edge: " << sIndex
3951 << "::" << sMesh.edges_[sIndex]
3952 << " on proc: " << procIndices_[pI]
3953 << " for master edge: " << eIndex
3954 << " using collapseCase: " << masterMap.type()
3955 << endl;
3956 }
3957 }
3958 }
3960 const Map<label>& pointMap = cMapPtr->entityMap(pointEnum);
3962 // Perform a topological comparison.
3963 switch (masterMap.type())
3964 {
3965 case 1:
3966 {
3967 if (sIndex < 0)
3968 {
3969 slaveMoveNewPoint[slaveI] = points_[mEdge[0]];
3970 slaveMoveOldPoint[slaveI] = oldPoints_[mEdge[0]];
3971 }
3972 else
3973 if (pointMap[mEdge[0]] == sEdge[0])
3974 {
3975 slaveOverRide = 1;
3976 }
3977 else
3978 if (pointMap[mEdge[1]] == sEdge[0])
3979 {
3980 slaveOverRide = 2;
3981 }
3982 else
3983 {
3984 // Write out for for post-processing
3985 writeVTK("mEdge_" + Foam::name(eIndex), eIndex, 1);
3987 FatalErrorIn
3988 (
3989 "\n"
3990 "const changeMap dynamicTopoFvMesh"
3991 "::collapseEdge\n"
3992 "(\n"
3993 " const label eIndex,\n"
3994 " label overRideCase,\n"
3995 " bool checkOnly,\n"
3996 " bool forceOp\n"
3997 ")\n"
3998 )
3999 << "Coupled collapse failed." << nl
4000 << "Master: " << eIndex << " : " << mEdge << nl
4001 << "Slave: " << sIndex << " : " << sEdge << nl
4002 << abort(FatalError);
4003 }
4005 break;
4006 }
4008 case 2:
4009 {
4010 if (sIndex < 0)
4011 {
4012 slaveMoveNewPoint[slaveI] = points_[mEdge[1]];
4013 slaveMoveOldPoint[slaveI] = oldPoints_[mEdge[1]];
4014 }
4015 else
4016 if (pointMap[mEdge[1]] == sEdge[1])
4017 {
4018 slaveOverRide = 2;
4019 }
4020 else
4021 if (pointMap[mEdge[0]] == sEdge[1])
4022 {
4023 slaveOverRide = 1;
4024 }
4025 else
4026 {
4027 // Write out for for post-processing
4028 writeVTK("mEdge_" + Foam::name(eIndex), eIndex, 1);
4030 FatalErrorIn
4031 (
4032 "\n"
4033 "const changeMap dynamicTopoFvMesh"
4034 "::collapseEdge\n"
4035 "(\n"
4036 " const label eIndex,\n"
4037 " label overRideCase,\n"
4038 " bool checkOnly,\n"
4039 " bool forceOp\n"
4040 ")\n"
4041 )
4042 << "Coupled collapse failed." << nl
4043 << "Master: " << eIndex << " : " << mEdge << nl
4044 << "Slave: " << sIndex << " : " << sEdge << nl
4045 << abort(FatalError);
4046 }
4048 break;
4049 }
4051 case 3:
4052 {
4053 if (sIndex < 0)
4054 {
4055 slaveMoveNewPoint[slaveI] =
4056 (
4057 0.5 * (points_[mEdge[0]] + points_[mEdge[1]])
4058 );
4060 slaveMoveOldPoint[slaveI] =
4061 (
4062 0.5 * (oldPoints_[mEdge[0]] + oldPoints_[mEdge[1]])
4063 );
4064 }
4065 else
4066 {
4067 slaveOverRide = 3;
4068 }
4070 break;
4071 }
4072 }
4074 // Temporarily turn off coupledModification
4075 unsetCoupledModification();
4077 // Test the slave edge
4078 if (localCouple)
4079 {
4080 slaveMap = collapseEdge(sIndex, slaveOverRide, true, forceOp);
4081 }
4082 else
4083 if (procCouple)
4084 {
4085 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
4087 if (sIndex < 0)
4088 {
4089 // Point-based coupling
4090 scalar slaveCollapseQuality(GREAT);
4091 DynamicList<label> cellsChecked(10);
4093 // Check cells connected to coupled point
4094 const labelList& pEdges = sMesh.pointEdges_[mag(sIndex)];
4096 bool infeasible = false;
4098 forAll(pEdges, edgeI)
4099 {
4100 const labelList& eFaces =
4101 (
4102 sMesh.edgeFaces_[pEdges[edgeI]]
4103 );
4105 // Build a list of cells to check
4106 forAll(eFaces, faceI)
4107 {
4108 label own = sMesh.owner_[eFaces[faceI]];
4109 label nei = sMesh.neighbour_[eFaces[faceI]];
4111 // Check owner cell
4112 if (findIndex(cellsChecked, own) == -1)
4113 {
4114 // Check if point movement is feasible
4115 if
4116 (
4117 sMesh.checkCollapse
4118 (
4119 slaveMoveNewPoint[slaveI],
4120 slaveMoveOldPoint[slaveI],
4121 mag(sIndex),
4122 own,
4123 cellsChecked,
4124 slaveCollapseQuality,
4125 forceOp
4126 )
4127 )
4128 {
4129 infeasible = true;
4130 break;
4131 }
4132 }
4134 if (nei == -1)
4135 {
4136 continue;
4137 }
4139 // Check neighbour cell
4140 if (findIndex(cellsChecked, nei) == -1)
4141 {
4142 // Check if point movement is feasible
4143 if
4144 (
4145 sMesh.checkCollapse
4146 (
4147 slaveMoveNewPoint[slaveI],
4148 slaveMoveOldPoint[slaveI],
4149 mag(sIndex),
4150 nei,
4151 cellsChecked,
4152 slaveCollapseQuality,
4153 forceOp
4154 )
4155 )
4156 {
4157 infeasible = true;
4158 break;
4159 }
4160 }
4161 }
4163 if (infeasible)
4164 {
4165 break;
4166 }
4167 }
4169 if (infeasible)
4170 {
4171 slaveMap.type() = 0;
4172 }
4173 else
4174 {
4175 slaveMap.type() = 1;
4176 }
4177 }
4178 else
4179 {
4180 // Edge-based coupling
4181 slaveMap =
4182 (
4183 sMesh.collapseEdge
4184 (
4185 sIndex,
4186 slaveOverRide,
4187 true,
4188 forceOp
4189 )
4190 );
4191 }
4192 }
4194 // Turn it back on.
4195 setCoupledModification();
4197 if (slaveMap.type() <= 0)
4198 {
4199 // Slave couldn't perform collapse.
4200 map.type() = -2;
4202 return map;
4203 }
4205 // Save index and patch for posterity
4206 slaveMap.index() = sIndex;
4207 slaveMap.patchIndex() = pI;
4208 }
4210 // Next collapse each slave edge (for real this time...)
4211 forAll(slaveMaps, slaveI)
4212 {
4213 // Alias for convenience...
4214 changeMap& slaveMap = slaveMaps[slaveI];
4216 label sIndex = slaveMap.index();
4217 label pI = slaveMap.patchIndex();
4219 // If the edge is mapped onto itself, skip modification
4220 // (occurs for cyclic edges)
4221 if ((sIndex == eIndex) && localCouple)
4222 {
4223 continue;
4224 }
4226 label slaveOverRide = slaveMap.type();
4228 // Temporarily turn off coupledModification
4229 unsetCoupledModification();
4231 // Collapse the slave
4232 if (localCouple)
4233 {
4234 slaveMap = collapseEdge(sIndex, slaveOverRide, false, forceOp);
4235 }
4236 else
4237 {
4238 const coupleMap& cMap = recvMeshes_[pI].map();
4239 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
4241 if (sIndex < 0)
4242 {
4243 // Point based coupling
4245 // Move points to new location,
4246 // and update operation into coupleMap
4247 sMesh.points_[mag(sIndex)] = slaveMoveNewPoint[slaveI];
4248 sMesh.oldPoints_[mag(sIndex)] = slaveMoveOldPoint[slaveI];
4250 cMap.pushOperation
4251 (
4252 mag(sIndex),
4253 coupleMap::MOVE_POINT,
4254 slaveMoveNewPoint[slaveI],
4255 slaveMoveOldPoint[slaveI]
4256 );
4258 // Force operation to succeed
4259 slaveMap.type() = 1;
4260 }
4261 else
4262 {
4263 // Edge-based coupling
4264 slaveMap =
4265 (
4266 sMesh.collapseEdge
4267 (
4268 sIndex,
4269 slaveOverRide,
4270 false,
4271 forceOp
4272 )
4273 );
4275 // Push operation into coupleMap
4276 switch (slaveMap.type())
4277 {
4278 case 1:
4279 {
4280 cMap.pushOperation
4281 (
4282 sIndex,
4283 coupleMap::COLLAPSE_FIRST
4284 );
4286 break;
4287 }
4289 case 2:
4290 {
4291 cMap.pushOperation
4292 (
4293 sIndex,
4294 coupleMap::COLLAPSE_SECOND
4295 );
4297 break;
4298 }
4300 case 3:
4301 {
4302 cMap.pushOperation
4303 (
4304 sIndex,
4305 coupleMap::COLLAPSE_MIDPOINT
4306 );
4308 break;
4309 }
4310 }
4311 }
4312 }
4314 // Turn it back on.
4315 setCoupledModification();
4317 // The final operation has to succeed.
4318 if (slaveMap.type() <= 0)
4319 {
4320 FatalErrorIn
4321 (
4322 "\n"
4323 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
4324 "(\n"
4325 " const label eIndex,\n"
4326 " label overRideCase,\n"
4327 " bool checkOnly,\n"
4328 " bool forceOp\n"
4329 ")\n"
4330 )
4331 << "Coupled topo-change for slave failed." << nl
4332 << " Edge: " << eIndex << ": " << eCheck << nl
4333 << " Slave index: " << sIndex << nl
4334 << " Patch index: " << pI << nl
4335 << " Type: " << slaveMap.type() << nl
4336 << abort(FatalError);
4337 }
4339 // Save index and patch for posterity
4340 slaveMap.index() = sIndex;
4341 slaveMap.patchIndex() = pI;
4342 }
4343 }
4345 // Build hullVertices for this edge
4346 labelList vertexHull;
4347 buildVertexHull(eIndex, vertexHull);
4349 // Hull variables
4350 bool found = false;
4351 label replaceIndex = -1, m = vertexHull.size();
4353 // Size up the hull lists
4354 labelList cellHull(m, -1);
4355 labelList faceHull(m, -1);
4356 labelList edgeHull(m, -1);
4357 labelListList ringEntities(4, labelList(m, -1));
4359 // Construct a hull around this edge
4360 meshOps::constructHull
4361 (
4362 eIndex,
4363 faces_,
4364 edges_,
4365 cells_,
4366 owner_,
4367 neighbour_,
4368 faceEdges_,
4369 edgeFaces_,
4370 vertexHull,
4371 edgeHull,
4372 faceHull,
4373 cellHull,
4374 ringEntities
4375 );
4377 // Check whether points of the edge lies on a boundary
4378 const FixedList<bool,2> edgeBoundary = checkEdgeBoundary(eIndex);
4379 FixedList<label, 2> nBoundCurves(0), nProcCurves(0), checkPoints(-1);
4381 // Decide on collapseCase
4382 label collapseCase = -1;
4384 if (edgeBoundary[0] && !edgeBoundary[1])
4385 {
4386 collapseCase = 1;
4387 }
4388 else
4389 if (!edgeBoundary[0] && edgeBoundary[1])
4390 {
4391 collapseCase = 2;
4392 }
4393 else
4394 if (edgeBoundary[0] && edgeBoundary[1])
4395 {
4396 // If this is an interior edge with two boundary points.
4397 // Bail out for now. If proximity based refinement is
4398 // switched on, mesh may be sliced at this point.
4399 label edgePatch = whichEdgePatch(eIndex);
4401 if (edgePatch == -1)
4402 {
4403 return map;
4404 }
4406 // Check if either point lies on a bounding curve
4407 // Used to ensure that collapses happen towards bounding curves.
4408 // If the edge itself is on a bounding curve, collapse is valid.
4409 const edge& edgeCheck = edges_[eIndex];
4411 forAll(edgeCheck, pointI)
4412 {
4413 const labelList& pEdges = pointEdges_[edgeCheck[pointI]];
4415 forAll(pEdges, edgeI)
4416 {
4417 if
4418 (
4419 checkBoundingCurve
4420 (
4421 pEdges[edgeI],
4422 false,
4423 &(nProcCurves[pointI])
4424 )
4425 )
4426 {
4427 nBoundCurves[pointI]++;
4428 }
4429 }
4430 }
4432 // Check for procCurves first
4433 if (!coupledModification_ && !isSubMesh_)
4434 {
4435 // Bail out for now
4436 if (nProcCurves[0] && nProcCurves[1])
4437 {
4438 return map;
4439 }
4441 if (nProcCurves[0] && !nProcCurves[1])
4442 {
4443 if (nBoundCurves[1])
4444 {
4445 // Bail out for now
4446 return map;
4447 }
4448 else
4449 {
4450 collapseCase = 1;
4451 }
4452 }
4454 if (!nProcCurves[0] && nProcCurves[1])
4455 {
4456 if (nBoundCurves[0])
4457 {
4458 // Bail out for now
4459 return map;
4460 }
4461 else
4462 {
4463 collapseCase = 2;
4464 }
4465 }
4466 }
4468 // If still undecided, choose based on bounding curves
4469 if (collapseCase == -1)
4470 {
4471 // Pick the point which is connected to more bounding curves
4472 if (nBoundCurves[0] > nBoundCurves[1])
4473 {
4474 collapseCase = 1;
4475 }
4476 else
4477 if (nBoundCurves[1] > nBoundCurves[0])
4478 {
4479 collapseCase = 2;
4480 }
4481 else
4482 {
4483 // Bounding edge: collapseEdge can collapse this edge
4484 collapseCase = 3;
4485 }
4486 }
4487 }
4488 else
4489 {
4490 // Looks like this is an interior edge.
4491 // Collapse case [3] by default
4492 collapseCase = 3;
4493 }
4495 // Perform an override if requested.
4496 if (overRideCase != -1)
4497 {
4498 collapseCase = overRideCase;
4499 }
4501 // Configure the new point-position
4502 point newPoint = vector::zero;
4503 point oldPoint = vector::zero;
4505 label collapsePoint = -1, replacePoint = -1;
4507 switch (collapseCase)
4508 {
4509 case 1:
4510 {
4511 // Collapse to the first node
4512 replacePoint = edges_[eIndex][0];
4513 collapsePoint = edges_[eIndex][1];
4515 newPoint = points_[replacePoint];
4516 oldPoint = oldPoints_[replacePoint];
4518 checkPoints[0] = collapsePoint;
4520 break;
4521 }
4523 case 2:
4524 {
4525 // Collapse to the second node
4526 replacePoint = edges_[eIndex][1];
4527 collapsePoint = edges_[eIndex][0];
4529 newPoint = points_[replacePoint];
4530 oldPoint = oldPoints_[replacePoint];
4532 checkPoints[0] = collapsePoint;
4534 break;
4535 }
4537 case 3:
4538 {
4539 // Collapse to the mid-point
4540 replacePoint = edges_[eIndex][1];
4541 collapsePoint = edges_[eIndex][0];
4543 newPoint =
4544 (
4545 0.5 *
4546 (
4547 points_[replacePoint]
4548 + points_[collapsePoint]
4549 )
4550 );
4552 oldPoint =
4553 (
4554 0.5 *
4555 (
4556 oldPoints_[replacePoint]
4557 + oldPoints_[collapsePoint]
4558 )
4559 );
4561 checkPoints[0] = replacePoint;
4562 checkPoints[1] = collapsePoint;
4564 break;
4565 }
4567 default:
4568 {
4569 // Don't think this will ever happen.
4570 FatalErrorIn
4571 (
4572 "\n"
4573 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
4574 "(\n"
4575 " const label eIndex,\n"
4576 " label overRideCase,\n"
4577 " bool checkOnly,\n"
4578 " bool forceOp\n"
4579 ")\n"
4580 )
4581 << "Edge: " << eIndex << ": " << edges_[eIndex]
4582 << ". Couldn't decide on collapseCase."
4583 << abort(FatalError);
4585 break;
4586 }
4587 }
4589 // Loop through edges and check for feasibility of collapse
4590 // Also, keep track of resulting cell quality,
4591 // if collapse is indeed feasible
4592 scalar collapseQuality(GREAT);
4593 DynamicList<label> cellsChecked(10);
4595 // Add all hull cells as 'checked',
4596 // and therefore, feasible
4597 forAll(cellHull, cellI)
4598 {
4599 if (cellHull[cellI] == -1)
4600 {
4601 continue;
4602 }
4604 cellsChecked.append(cellHull[cellI]);
4605 }
4607 // Check collapsibility of cells around edges
4608 // with the re-configured point
4609 forAll(checkPoints, pointI)
4610 {
4611 if (checkPoints[pointI] == -1)
4612 {
4613 continue;
4614 }
4616 const labelList& checkPointEdges = pointEdges_[checkPoints[pointI]];
4618 forAll(checkPointEdges, edgeI)
4619 {
4620 const labelList& eFaces = edgeFaces_[checkPointEdges[edgeI]];
4622 // Build a list of cells to check
4623 forAll(eFaces, faceI)
4624 {
4625 label own = owner_[eFaces[faceI]];
4626 label nei = neighbour_[eFaces[faceI]];
4628 // Check owner cell
4629 if (findIndex(cellsChecked, own) == -1)
4630 {
4631 // Check if a collapse is feasible
4632 if
4633 (
4634 checkCollapse
4635 (
4636 newPoint,
4637 oldPoint,
4638 checkPoints[pointI],
4639 own,
4640 cellsChecked,
4641 collapseQuality,
4642 forceOp
4643 )
4644 )
4645 {
4646 map.type() = 0;
4647 return map;
4648 }
4649 }
4651 if (nei == -1)
4652 {
4653 continue;
4654 }
4656 // Check neighbour cell
4657 if (findIndex(cellsChecked, nei) == -1)
4658 {
4659 // Check if a collapse is feasible
4660 if
4661 (
4662 checkCollapse
4663 (
4664 newPoint,
4665 oldPoint,
4666 checkPoints[pointI],
4667 nei,
4668 cellsChecked,
4669 collapseQuality,
4670 forceOp
4671 )
4672 )
4673 {
4674 map.type() = 0;
4675 return map;
4676 }
4677 }
4678 }
4679 }
4680 }
4682 // Add a map entry of the replacePoint as an 'addedPoint'
4683 // - Used in coupled mapping
4684 map.addPoint(replacePoint);
4686 // Are we only performing checks?
4687 if (checkOnly)
4688 {
4689 // Return a succesful collapse possibility
4690 map.type() = collapseCase;
4692 // Make note of the removed point
4693 map.removePoint(collapsePoint);
4695 if (debug > 2)
4696 {
4697 Pout<< nl << "Edge: " << eIndex
4698 << ":: " << edges_[eIndex] << nl
4699 << " collapseCase determined to be: " << collapseCase << nl
4700 << " Resulting quality: " << collapseQuality << nl
4701 << " collapsePoint: " << collapsePoint << nl
4702 << " nBoundCurves: " << nBoundCurves << nl
4703 << " nProcCurves: " << nProcCurves << nl
4704 << endl;
4705 }
4707 return map;
4708 }
4710 // Define indices on the hull for removal / replacement
4711 label removeEdgeIndex = -1, replaceEdgeIndex = -1;
4712 label removeFaceIndex = -1, replaceFaceIndex = -1;
4714 if (replacePoint == edges_[eIndex][0])
4715 {
4716 replaceEdgeIndex = 0;
4717 replaceFaceIndex = 1;
4718 removeEdgeIndex = 2;
4719 removeFaceIndex = 3;
4720 }
4721 else
4722 if (replacePoint == edges_[eIndex][1])
4723 {
4724 removeEdgeIndex = 0;
4725 removeFaceIndex = 1;
4726 replaceEdgeIndex = 2;
4727 replaceFaceIndex = 3;
4728 }
4729 else
4730 {
4731 // Don't think this will ever happen.
4732 FatalErrorIn
4733 (
4734 "\n"
4735 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
4736 "(\n"
4737 " const label eIndex,\n"
4738 " label overRideCase,\n"
4739 " bool checkOnly,\n"
4740 " bool forceOp\n"
4741 ")\n"
4742 )
4743 << "Edge: " << eIndex << ": " << edges_[eIndex]
4744 << ". Couldn't decide on removal / replacement indices."
4745 << abort(FatalError);
4746 }
4748 if (debug > 1)
4749 {
4750 Pout<< nl << nl
4751 << "Edge: " << eIndex << ": " << edges_[eIndex]
4752 << " is to be collapsed. " << nl;
4754 Pout<< " On SubMesh: " << Switch::asText(isSubMesh_) << nl;
4755 Pout<< " coupledModification: " << coupledModification_ << nl;
4757 label epIndex = whichEdgePatch(eIndex);
4759 const polyBoundaryMesh& boundary = boundaryMesh();
4761 Pout<< " Patch: ";
4763 if (epIndex == -1)
4764 {
4765 Pout<< "Internal" << nl;
4766 }
4767 else
4768 if (epIndex < boundary.size())
4769 {
4770 Pout<< boundary[epIndex].name() << nl;
4771 }
4772 else
4773 {
4774 Pout<< " New patch: " << epIndex << endl;
4775 }
4777 Pout<< " nBoundCurves: " << nBoundCurves << nl
4778 << " nProcCurves: " << nProcCurves << nl
4779 << " collapseCase: " << collapseCase << nl
4780 << " Resulting quality: " << collapseQuality << endl;
4782 if (debug > 2)
4783 {
4784 Pout<< " Edges: " << edgeHull << nl
4785 << " Faces: " << faceHull << nl
4786 << " Cells: " << cellHull << nl
4787 << " replacePoint: " << replacePoint << nl
4788 << " collapsePoint: " << collapsePoint << nl
4789 << " checkPoints: " << checkPoints << nl;
4791 Pout<< " ringEntities (removed faces): " << nl;
4793 forAll(ringEntities[removeFaceIndex], faceI)
4794 {
4795 label fIndex = ringEntities[removeFaceIndex][faceI];
4797 if (fIndex == -1)
4798 {
4799 continue;
4800 }
4802 Pout<< fIndex << ": " << faces_[fIndex] << nl;
4803 }
4805 Pout<< " ringEntities (removed edges): " << nl;
4807 forAll(ringEntities[removeEdgeIndex], edgeI)
4808 {
4809 label ieIndex = ringEntities[removeEdgeIndex][edgeI];
4811 if (ieIndex == -1)
4812 {
4813 continue;
4814 }
4816 Pout<< ieIndex << ": " << edges_[ieIndex] << nl;
4817 }
4819 Pout<< " ringEntities (replacement faces): " << nl;
4821 forAll(ringEntities[replaceFaceIndex], faceI)
4822 {
4823 label fIndex = ringEntities[replaceFaceIndex][faceI];
4825 if (fIndex == -1)
4826 {
4827 continue;
4828 }
4830 Pout<< fIndex << ": " << faces_[fIndex] << nl;
4831 }
4833 Pout<< " ringEntities (replacement edges): " << nl;
4835 forAll(ringEntities[replaceEdgeIndex], edgeI)
4836 {
4837 label ieIndex = ringEntities[replaceEdgeIndex][edgeI];
4839 if (ieIndex == -1)
4840 {
4841 continue;
4842 }
4844 Pout<< ieIndex << ": " << edges_[ieIndex] << nl;
4845 }
4847 labelList& collapsePointEdges = pointEdges_[collapsePoint];
4849 Pout<< " pointEdges (collapsePoint): ";
4851 forAll(collapsePointEdges, edgeI)
4852 {
4853 Pout<< collapsePointEdges[edgeI] << " ";
4854 }
4856 Pout<< nl;
4858 // Write out VTK files prior to change
4859 // - Using old-points for convenience in post-processing
4860 if (debug > 3)
4861 {
4862 writeVTK
4863 (
4864 Foam::name(eIndex)
4865 + '(' + Foam::name(collapsePoint)
4866 + ',' + Foam::name(replacePoint) + ')'
4867 + "_Collapse_0",
4868 cellsChecked,
4869 3, false, true
4870 );
4871 }
4872 }
4873 }
4875 if (whichEdgePatch(eIndex) > -1)
4876 {
4877 // Update number of surface collapses, if necessary.
4878 statistics_[6]++;
4879 }
4881 // Maintain a list of modified faces for mapping
4882 DynamicList<label> modifiedFaces(10);
4884 // Renumber all hull faces and edges
4885 forAll(faceHull, indexI)
4886 {
4887 // Loop through all faces of the edge to be removed
4888 // and reassign them to the replacement edge
4889 label edgeToRemove = ringEntities[removeEdgeIndex][indexI];
4890 label faceToRemove = ringEntities[removeFaceIndex][indexI];
4891 label cellToRemove = cellHull[indexI];
4892 label replaceEdge = ringEntities[replaceEdgeIndex][indexI];
4893 label replaceFace = ringEntities[replaceFaceIndex][indexI];
4895 label replacePatch = whichEdgePatch(replaceEdge);
4896 label removePatch = whichEdgePatch(edgeToRemove);
4898 // Check if a patch conversion is necessary
4899 if (replacePatch == -1 && removePatch > -1)
4900 {
4901 if (debug > 2)
4902 {
4903 Pout<< " Edge: " << replaceEdge
4904 << " :: " << edges_[replaceEdge]
4905 << " is interior, but edge: " << edgeToRemove
4906 << " :: " << edges_[edgeToRemove]
4907 << " is on a boundary patch."
4908 << endl;
4909 }
4911 // Convert patch for edge
4912 edge newEdge = edges_[replaceEdge];
4913 labelList newEdgeFaces = edgeFaces_[replaceEdge];
4915 // Insert the new edge
4916 label newEdgeIndex =
4917 (
4918 insertEdge
4919 (
4920 removePatch,
4921 newEdge,
4922 newEdgeFaces
4923 )
4924 );
4926 // Replace faceEdges for all
4927 // connected faces.
4928 forAll(newEdgeFaces, faceI)
4929 {
4930 meshOps::replaceLabel
4931 (
4932 replaceEdge,
4933 newEdgeIndex,
4934 faceEdges_[newEdgeFaces[faceI]]
4935 );
4936 }
4938 // Remove the edge
4939 removeEdge(replaceEdge);
4941 // Update map
4942 map.removeEdge(replaceEdge);
4943 map.addEdge(newEdgeIndex, labelList(1, replaceEdge));
4945 // Replace index and patch
4946 replaceEdge = newEdgeIndex;
4948 // Modify ringEntities
4949 ringEntities[replaceEdgeIndex][indexI] = newEdgeIndex;
4950 }
4952 const labelList& rmvEdgeFaces = edgeFaces_[edgeToRemove];
4954 forAll(rmvEdgeFaces, faceI)
4955 {
4956 // Replace edge labels for faces
4957 meshOps::replaceLabel
4958 (
4959 edgeToRemove,
4960 replaceEdge,
4961 faceEdges_[rmvEdgeFaces[faceI]]
4962 );
4964 // Loop through faces associated with this edge,
4965 // and renumber them as well.
4966 const face& faceToCheck = faces_[rmvEdgeFaces[faceI]];
4968 if ((replaceIndex = faceToCheck.which(collapsePoint)) > -1)
4969 {
4970 // Renumber the face...
4971 faces_[rmvEdgeFaces[faceI]][replaceIndex] = replacePoint;
4973 // Add an entry for mapping
4974 if (findIndex(modifiedFaces, rmvEdgeFaces[faceI]) == -1)
4975 {
4976 modifiedFaces.append(rmvEdgeFaces[faceI]);
4977 }
4978 }
4980 // Hull faces should be removed for the replacement edge
4981 if (rmvEdgeFaces[faceI] == faceHull[indexI])
4982 {
4983 meshOps::sizeDownList
4984 (
4985 faceHull[indexI],
4986 edgeFaces_[replaceEdge]
4987 );
4989 continue;
4990 }
4992 found = false;
4994 // Need to avoid ring faces as well.
4995 forAll(ringEntities[removeFaceIndex], faceII)
4996 {
4997 if
4998 (
4999 rmvEdgeFaces[faceI]
5000 == ringEntities[removeFaceIndex][faceII]
5001 )
5002 {
5003 found = true;
5004 break;
5005 }
5006 }
5008 // Size-up the replacement edge list if the face hasn't been found.
5009 // These faces are connected to the edge slated for
5010 // removal, but do not belong to the hull.
5011 if (!found)
5012 {
5013 meshOps::sizeUpList
5014 (
5015 rmvEdgeFaces[faceI],
5016 edgeFaces_[replaceEdge]
5017 );
5018 }
5019 }
5021 if (cellToRemove == -1)
5022 {
5023 continue;
5024 }
5026 // Size down edgeFaces for the ring edges
5027 meshOps::sizeDownList
5028 (
5029 faceToRemove,
5030 edgeFaces_[edgeHull[indexI]]
5031 );
5033 // Ensure proper orientation of retained faces
5034 if (owner_[faceToRemove] == cellToRemove)
5035 {
5036 if (owner_[replaceFace] == cellToRemove)
5037 {
5038 if
5039 (
5040 (neighbour_[faceToRemove] > neighbour_[replaceFace])
5041 && (neighbour_[replaceFace] != -1)
5042 )
5043 {
5044 // This face is to be flipped
5045 faces_[replaceFace] = faces_[replaceFace].reverseFace();
5046 owner_[replaceFace] = neighbour_[replaceFace];
5047 neighbour_[replaceFace] = neighbour_[faceToRemove];
5049 setFlip(replaceFace);
5050 }
5051 else
5052 if
5053 (
5054 (neighbour_[faceToRemove] == -1) &&
5055 (neighbour_[replaceFace] != -1)
5056 )
5057 {
5058 // This interior face would need to be converted
5059 // to a boundary one, and flipped as well.
5060 face newFace = faces_[replaceFace].reverseFace();
5061 label newOwner = neighbour_[replaceFace];
5062 label newNeighbour = neighbour_[faceToRemove];
5063 labelList newFE = faceEdges_[replaceFace];
5065 label newFaceIndex =
5066 (
5067 insertFace
5068 (
5069 whichPatch(faceToRemove),
5070 newFace,
5071 newOwner,
5072 newNeighbour
5073 )
5074 );
5076 // Set this face aside for mapping
5077 modifiedFaces.append(newFaceIndex);
5079 // Update map.
5080 map.addFace(newFaceIndex, labelList(1, faceToRemove));
5082 // Ensure that all edges of this face are
5083 // on the boundary.
5084 forAll(newFE, edgeI)
5085 {
5086 if (whichEdgePatch(newFE[edgeI]) == -1)
5087 {
5088 edge newEdge = edges_[newFE[edgeI]];
5089 labelList newEF = edgeFaces_[newFE[edgeI]];
5091 // Need patch information for the new edge.
5092 // Find the corresponding edge in ringEntities.
5093 // Note that hullEdges doesn't need to be checked,
5094 // since they are common to both faces.
5095 label i =
5096 (
5097 findIndex
5098 (
5099 ringEntities[replaceEdgeIndex],
5100 newFE[edgeI]
5101 )
5102 );
5104 label repIndex =
5105 (
5106 whichEdgePatch
5107 (
5108 ringEntities[removeEdgeIndex][i]
5109 )
5110 );
5112 // Insert the new edge
5113 label newEdgeIndex =
5114 (
5115 insertEdge
5116 (
5117 repIndex,
5118 newEdge,
5119 newEF
5120 )
5121 );
5123 // Replace faceEdges for all
5124 // connected faces.
5125 forAll(newEF, faceI)
5126 {
5127 meshOps::replaceLabel
5128 (
5129 newFE[edgeI],
5130 newEdgeIndex,
5131 faceEdges_[newEF[faceI]]
5132 );
5133 }
5135 // Remove the edge
5136 removeEdge(newFE[edgeI]);
5138 // Update map
5139 map.removeEdge(newFE[edgeI]);
5141 map.addEdge
5142 (
5143 newEdgeIndex,
5144 labelList(1, newFE[edgeI])
5145 );
5147 // Replace faceEdges with new edge index
5148 newFE[edgeI] = newEdgeIndex;
5150 // Modify ringEntities
5151 ringEntities[replaceEdgeIndex][i] = newEdgeIndex;
5152 }
5153 }
5155 // Add the new faceEdges
5156 faceEdges_.append(newFE);
5158 // Replace edgeFaces with the new face index
5159 const labelList& newFEdges = faceEdges_[newFaceIndex];
5161 forAll(newFEdges, edgeI)
5162 {
5163 meshOps::replaceLabel
5164 (
5165 replaceFace,
5166 newFaceIndex,
5167 edgeFaces_[newFEdges[edgeI]]
5168 );
5169 }
5171 // Remove the face
5172 removeFace(replaceFace);
5174 // Update map
5175 map.removeFace(replaceFace);
5177 // Replace label for the new owner
5178 meshOps::replaceLabel
5179 (
5180 replaceFace,
5181 newFaceIndex,
5182 cells_[newOwner]
5183 );
5185 // Modify ringEntities and replaceFace
5186 replaceFace = newFaceIndex;
5187 ringEntities[replaceFaceIndex][indexI] = newFaceIndex;
5188 }
5189 else
5190 if
5191 (
5192 (neighbour_[faceToRemove] == -1) &&
5193 (neighbour_[replaceFace] == -1)
5194 )
5195 {
5196 // Wierd overhanging cell. Since replaceFace
5197 // would be an orphan if this continued, remove
5198 // the face entirely.
5199 labelList rmFE = faceEdges_[replaceFace];
5201 forAll(rmFE, edgeI)
5202 {
5203 if
5204 (
5205 (edgeFaces_[rmFE[edgeI]].size() == 1) &&
5206 (edgeFaces_[rmFE[edgeI]][0] == replaceFace)
5207 )
5208 {
5209 // This edge has to be removed entirely.
5210 removeEdge(rmFE[edgeI]);
5212 // Update map
5213 map.removeEdge(rmFE[edgeI]);
5215 label i =
5216 (
5217 findIndex
5218 (
5219 ringEntities[replaceEdgeIndex],
5220 rmFE[edgeI]
5221 )
5222 );
5224 if (i > -1)
5225 {
5226 // Modify ringEntities
5227 ringEntities[replaceEdgeIndex][i] = -1;
5228 }
5229 }
5230 else
5231 {
5232 // Size-down edgeFaces
5233 meshOps::sizeDownList
5234 (
5235 replaceFace,
5236 edgeFaces_[rmFE[edgeI]]
5237 );
5238 }
5239 }
5241 // If this is a subMesh, and replaceFace is on
5242 // a physical boundary, make a 'special' entry
5243 // for coupled mapping purposes.
5244 if (isSubMesh_)
5245 {
5246 label kfPatch = whichPatch(replaceFace);
5247 label rfPatch = whichPatch(faceToRemove);
5249 if
5250 (
5251 (getNeighbourProcessor(rfPatch) > -1) &&
5252 (getNeighbourProcessor(kfPatch) == -1)
5253 )
5254 {
5255 map.addFace
5256 (
5257 faceToRemove,
5258 labelList(1, (-2 - kfPatch))
5259 );
5260 }
5261 }
5263 // Remove the face
5264 removeFace(replaceFace);
5266 // Update map
5267 map.removeFace(replaceFace);
5269 // Modify ringEntities and replaceFace
5270 replaceFace = -1;
5271 ringEntities[replaceFaceIndex][indexI] = -1;
5272 }
5273 else
5274 {
5275 // Keep orientation intact, and update the owner
5276 owner_[replaceFace] = neighbour_[faceToRemove];
5277 }
5278 }
5279 else
5280 if (neighbour_[faceToRemove] == -1)
5281 {
5282 // This interior face would need to be converted
5283 // to a boundary one, but with orientation intact.
5284 face newFace = faces_[replaceFace];
5285 label newOwner = owner_[replaceFace];
5286 label newNeighbour = neighbour_[faceToRemove];
5287 labelList newFE = faceEdges_[replaceFace];
5289 label newFaceIndex =
5290 (
5291 insertFace
5292 (
5293 whichPatch(faceToRemove),
5294 newFace,
5295 newOwner,
5296 newNeighbour
5297 )
5298 );
5300 // Set this face aside for mapping
5301 modifiedFaces.append(newFaceIndex);
5303 // Update map
5304 map.addFace(newFaceIndex, labelList(1, faceToRemove));
5306 // Ensure that all edges of this face are
5307 // on the boundary.
5308 forAll(newFE, edgeI)
5309 {
5310 if (whichEdgePatch(newFE[edgeI]) == -1)
5311 {
5312 edge newEdge = edges_[newFE[edgeI]];
5313 labelList newEF = edgeFaces_[newFE[edgeI]];
5315 // Need patch information for the new edge.
5316 // Find the corresponding edge in ringEntities.
5317 // Note that hullEdges doesn't need to be checked,
5318 // since they are common to both faces.
5319 label i =
5320 (
5321 findIndex
5322 (
5323 ringEntities[replaceEdgeIndex],
5324 newFE[edgeI]
5325 )
5326 );
5328 label repIndex =
5329 (
5330 whichEdgePatch
5331 (
5332 ringEntities[removeEdgeIndex][i]
5333 )
5334 );
5336 // Insert the new edge
5337 label newEdgeIndex =
5338 (
5339 insertEdge
5340 (
5341 repIndex,
5342 newEdge,
5343 newEF
5344 )
5345 );
5347 // Replace faceEdges for all
5348 // connected faces.
5349 forAll(newEF, faceI)
5350 {
5351 meshOps::replaceLabel
5352 (
5353 newFE[edgeI],
5354 newEdgeIndex,
5355 faceEdges_[newEF[faceI]]
5356 );
5357 }
5359 // Remove the edge
5360 removeEdge(newFE[edgeI]);
5362 // Update map
5363 map.removeEdge(newFE[edgeI]);
5365 map.addEdge
5366 (
5367 newEdgeIndex,
5368 labelList(1, newFE[edgeI])
5369 );
5371 // Replace faceEdges with new edge index
5372 newFE[edgeI] = newEdgeIndex;
5374 // Modify ringEntities
5375 ringEntities[replaceEdgeIndex][i] = newEdgeIndex;
5376 }
5377 }
5379 // Add the new faceEdges
5380 faceEdges_.append(newFE);
5382 // Replace edgeFaces with the new face index
5383 const labelList& newFEdges = faceEdges_[newFaceIndex];
5385 forAll(newFEdges, edgeI)
5386 {
5387 meshOps::replaceLabel
5388 (
5389 replaceFace,
5390 newFaceIndex,
5391 edgeFaces_[newFEdges[edgeI]]
5392 );
5393 }
5395 // Remove the face
5396 removeFace(replaceFace);
5398 // Update map
5399 map.removeFace(replaceFace);
5401 // Replace label for the new owner
5402 meshOps::replaceLabel
5403 (
5404 replaceFace,
5405 newFaceIndex,
5406 cells_[newOwner]
5407 );
5409 // Modify ringEntities and replaceFace
5410 replaceFace = newFaceIndex;
5411 ringEntities[replaceFaceIndex][indexI] = newFaceIndex;
5412 }
5413 else
5414 {
5415 // Keep orientation intact, and update the neighbour
5416 neighbour_[replaceFace] = neighbour_[faceToRemove];
5417 }
5419 // Update the cell
5420 if (neighbour_[faceToRemove] != -1)
5421 {
5422 meshOps::replaceLabel
5423 (
5424 faceToRemove,
5425 replaceFace,
5426 cells_[neighbour_[faceToRemove]]
5427 );
5428 }
5429 }
5430 else
5431 {
5432 if (neighbour_[replaceFace] == cellToRemove)
5433 {
5434 if (owner_[faceToRemove] < owner_[replaceFace])
5435 {
5436 // This face is to be flipped
5437 faces_[replaceFace] = faces_[replaceFace].reverseFace();
5438 neighbour_[replaceFace] = owner_[replaceFace];
5439 owner_[replaceFace] = owner_[faceToRemove];
5441 setFlip(replaceFace);
5442 }
5443 else
5444 {
5445 // Keep orientation intact, and update the neighbour
5446 neighbour_[replaceFace] = owner_[faceToRemove];
5447 }
5448 }
5449 else
5450 {
5451 // Keep orientation intact, and update the owner
5452 owner_[replaceFace] = owner_[faceToRemove];
5453 }
5455 // Update the cell
5456 meshOps::replaceLabel
5457 (
5458 faceToRemove,
5459 replaceFace,
5460 cells_[owner_[faceToRemove]]
5461 );
5462 }
5463 }
5465 // Remove all hull entities
5466 forAll(faceHull, indexI)
5467 {
5468 label edgeToRemove = ringEntities[removeEdgeIndex][indexI];
5469 label faceToRemove = ringEntities[removeFaceIndex][indexI];
5470 label cellToRemove = cellHull[indexI];
5472 if (cellToRemove != -1)
5473 {
5474 // Remove faceToRemove and associated faceEdges
5475 removeFace(faceToRemove);
5477 // Update map
5478 map.removeFace(faceToRemove);
5480 // Remove the hull cell
5481 removeCell(cellToRemove);
5483 // Update map
5484 map.removeCell(cellToRemove);
5485 }
5487 // Remove the hull edge and associated edgeFaces
5488 removeEdge(edgeToRemove);
5490 // Update map
5491 map.removeEdge(edgeToRemove);
5493 // Remove the hull face
5494 removeFace(faceHull[indexI]);
5496 // Update map
5497 map.removeFace(faceHull[indexI]);
5498 }
5500 // Loop through pointEdges for the collapsePoint,
5501 // and replace all occurrences with replacePoint.
5502 // Size-up pointEdges for the replacePoint as well.
5503 const labelList& pEdges = pointEdges_[collapsePoint];
5505 forAll(pEdges, edgeI)
5506 {
5507 // Renumber edges
5508 label edgeIndex = pEdges[edgeI];
5510 if (edgeIndex != eIndex)
5511 {
5512 if (edges_[edgeIndex][0] == collapsePoint)
5513 {
5514 edges_[edgeIndex][0] = replacePoint;
5516 meshOps::sizeUpList
5517 (
5518 edgeIndex,
5519 pointEdges_[replacePoint]
5520 );
5521 }
5522 else
5523 if (edges_[edgeIndex][1] == collapsePoint)
5524 {
5525 edges_[edgeIndex][1] = replacePoint;
5527 meshOps::sizeUpList
5528 (
5529 edgeIndex,
5530 pointEdges_[replacePoint]
5531 );
5532 }
5533 else
5534 {
5535 // Looks like pointEdges is inconsistent
5536 FatalErrorIn
5537 (
5538 "\n"
5539 "const changeMap dynamicTopoFvMesh::collapseEdge\n"
5540 "(\n"
5541 " const label eIndex,\n"
5542 " label overRideCase,\n"
5543 " bool checkOnly,\n"
5544 " bool forceOp\n"
5545 ")\n"
5546 )
5547 << "pointEdges is inconsistent." << nl
5548 << "Point: " << collapsePoint << nl
5549 << "pointEdges: " << pEdges << nl
5550 << abort(FatalError);
5551 }
5553 // Loop through faces associated with this edge,
5554 // and renumber them as well.
5555 const labelList& eFaces = edgeFaces_[edgeIndex];
5557 forAll(eFaces, faceI)
5558 {
5559 const face& faceToCheck = faces_[eFaces[faceI]];
5561 if ((replaceIndex = faceToCheck.which(collapsePoint)) > -1)
5562 {
5563 // Renumber the face...
5564 faces_[eFaces[faceI]][replaceIndex] = replacePoint;
5566 // Set this face aside for mapping
5567 if (findIndex(modifiedFaces, eFaces[faceI]) == -1)
5568 {
5569 modifiedFaces.append(eFaces[faceI]);
5570 }
5571 }
5572 }
5573 }
5574 }
5576 // Set old / new points
5577 oldPoints_[replacePoint] = oldPoint;
5578 points_[replacePoint] = newPoint;
5580 // Remove the collapse point
5581 removePoint(collapsePoint);
5583 // Update map
5584 map.removePoint(collapsePoint);
5586 // Remove the edge
5587 removeEdge(eIndex);
5589 // Update map
5590 map.removeEdge(eIndex);
5592 // Check for duplicate edges connected to the replacePoint
5593 const labelList& rpEdges = pointEdges_[replacePoint];
5595 DynamicList<label> mergeFaces(10);
5597 forAll(rpEdges, edgeI)
5598 {
5599 const edge& eCheckI = edges_[rpEdges[edgeI]];
5600 const label vCheckI = eCheckI.otherVertex(replacePoint);
5602 for (label edgeJ = edgeI + 1; edgeJ < rpEdges.size(); edgeJ++)
5603 {
5604 const edge& eCheckJ = edges_[rpEdges[edgeJ]];
5605 const label vCheckJ = eCheckJ.otherVertex(replacePoint);
5607 if (vCheckJ == vCheckI)
5608 {
5609 labelPair efCheck(rpEdges[edgeI], rpEdges[edgeJ]);
5611 forAll(efCheck, edgeI)
5612 {
5613 const labelList& eF = edgeFaces_[efCheck[edgeI]];
5615 forAll(eF, faceI)
5616 {
5617 label patch = whichPatch(eF[faceI]);
5619 if (patch == -1)
5620 {
5621 continue;
5622 }
5624 if (findIndex(mergeFaces, eF[faceI]) == -1)
5625 {
5626 mergeFaces.append(eF[faceI]);
5627 }
5628 }
5629 }
5631 if (debug > 2)
5632 {
5633 Pout<< " Found duplicate: " << eCheckI << nl
5634 << " Merge faces: " << mergeFaces << nl
5635 << endl;
5636 }
5637 }
5638 }
5639 }
5641 // Merge faces if necessary
5642 if (mergeFaces.size())
5643 {
5644 mergeBoundaryFaces(mergeFaces);
5645 }
5647 // Check and remove edges with an empty edgeFaces list
5648 const labelList& replaceEdges = ringEntities[replaceEdgeIndex];
5650 forAll(replaceEdges, edgeI)
5651 {
5652 label replaceEdge = replaceEdges[edgeI];
5654 // If the ring edge was removed, don't bother.
5655 if (replaceEdge > -1)
5656 {
5657 // Account for merged edges as well
5658 if (edgeFaces_[replaceEdge].empty())
5659 {
5660 // Is this edge truly removed? If not, remove it.
5661 if (edges_[replaceEdge] != edge(-1, -1))
5662 {
5663 if (debug > 2)
5664 {
5665 Pout<< " Edge: " << replaceEdge
5666 << " :: " << edges_[replaceEdge]
5667 << " has empty edgeFaces."
5668 << endl;
5669 }
5671 // Remove the edge
5672 removeEdge(replaceEdge);
5674 // Update map
5675 map.removeEdge(replaceEdge);
5676 }
5677 }
5678 }
5679 }
5681 // For cell-mapping, exclude all hull-cells
5682 forAll(cellsChecked, indexI)
5683 {
5684 if (cells_[cellsChecked[indexI]].empty())
5685 {
5686 cellsChecked[indexI] = -1;
5687 }
5688 else
5689 if (findIndex(cellHull, cellsChecked[indexI]) > 0)
5690 {
5691 cellsChecked[indexI] = -1;
5692 }
5693 }
5695 // Write out VTK files after change
5696 // - Using old-points for convenience in post-processing
5697 if (debug > 3)
5698 {
5699 writeVTK
5700 (
5701 Foam::name(eIndex)
5702 + '(' + Foam::name(collapsePoint)
5703 + ',' + Foam::name(replacePoint) + ')'
5704 + "_Collapse_1",
5705 cellsChecked,
5706 3, false, true
5707 );
5708 }
5710 // Now that all old / new cells possess correct connectivity,
5711 // compute mapping information.
5712 forAll(cellsChecked, cellI)
5713 {
5714 if (cellsChecked[cellI] < 0)
5715 {
5716 continue;
5717 }
5719 // Fill-in candidate mapping information
5720 labelList mC(1, cellsChecked[cellI]);
5722 // Set the mapping for this cell
5723 setCellMapping(cellsChecked[cellI], mC);
5724 }
5726 // Set face mapping information for modified faces
5727 forAll(modifiedFaces, faceI)
5728 {
5729 const label mfIndex = modifiedFaces[faceI];
5731 // Exclude deleted faces
5732 if (faces_[mfIndex].empty())
5733 {
5734 continue;
5735 }
5737 // Decide between default / weighted mapping
5738 // based on boundary information
5739 label fPatch = whichPatch(mfIndex);
5741 if (fPatch == -1)
5742 {
5743 setFaceMapping(mfIndex);
5744 }
5745 else
5746 {
5747 // Fill-in candidate mapping information
5748 labelList faceCandidates;
5750 const labelList& fEdges = faceEdges_[mfIndex];
5752 forAll(fEdges, edgeI)
5753 {
5754 if (whichEdgePatch(fEdges[edgeI]) == fPatch)
5755 {
5756 // Loop through associated edgeFaces
5757 const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
5759 forAll(eFaces, faceI)
5760 {
5761 if
5762 (
5763 (eFaces[faceI] != mfIndex) &&
5764 (whichPatch(eFaces[faceI]) == fPatch)
5765 )
5766 {
5767 faceCandidates.setSize
5768 (
5769 faceCandidates.size() + 1,
5770 eFaces[faceI]
5771 );
5772 }
5773 }
5774 }
5775 }
5777 // Set the mapping for this face
5778 setFaceMapping(mfIndex, faceCandidates);
5779 }
5780 }
5782 // If modification is coupled, update mapping info.
5783 if (coupledModification_)
5784 {
5785 // Build a list of boundary edges / faces for mapping
5786 DynamicList<label> checkEdges(10), checkFaces(10);
5788 const labelList& pEdges = pointEdges_[replacePoint];
5790 forAll(pEdges, edgeI)
5791 {
5792 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
5794 forAll(eFaces, faceI)
5795 {
5796 label fPatch = whichPatch(eFaces[faceI]);
5798 if (localCouple && !procCouple)
5799 {
5800 if (!locallyCoupledEntity(eFaces[faceI], true, false, true))
5801 {
5802 continue;
5803 }
5805 // Check for cyclics
5806 if (boundaryMesh()[fPatch].type() == "cyclic")
5807 {
5808 // Check if this is a master face
5809 const coupleMap& cM = patchCoupling_[fPatch].map();
5810 const Map<label>& fM = cM.entityMap(coupleMap::FACE);
5812 // Only add master faces
5813 // (belonging to the first half)
5814 if (!fM.found(eFaces[faceI]))
5815 {
5816 continue;
5817 }
5818 }
5819 }
5820 else
5821 if (procCouple && !localCouple)
5822 {
5823 if (getNeighbourProcessor(fPatch) == -1)
5824 {
5825 continue;
5826 }
5827 }
5829 // Add face and its edges for checking
5830 if (findIndex(checkFaces, eFaces[faceI]) == -1)
5831 {
5832 // Add this face
5833 checkFaces.append(eFaces[faceI]);
5835 const labelList& fEdges = faceEdges_[eFaces[faceI]];
5837 forAll(fEdges, edgeJ)
5838 {
5839 if (findIndex(checkEdges, fEdges[edgeJ]) == -1)
5840 {
5841 checkEdges.append(fEdges[edgeJ]);
5842 }
5843 }
5844 }
5845 }
5846 }
5848 // Output check entities
5849 if (debug > 4)
5850 {
5851 writeVTK
5852 (
5853 "checkEdges_" + Foam::name(eIndex),
5854 checkEdges, 1, false, true
5855 );
5857 writeVTK
5858 (
5859 "checkFaces_" + Foam::name(eIndex),
5860 checkFaces, 2, false, true
5861 );
5862 }
5864 if (localCouple && !procCouple)
5865 {
5866 // Alias for convenience...
5867 const changeMap& slaveMap = slaveMaps[0];
5869 const label pI = slaveMap.patchIndex();
5870 const coupleMap& cMap = patchCoupling_[pI].map();
5872 // Obtain references
5873 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
5874 Map<label>& rPointMap = cMap.reverseEntityMap(coupleMap::POINT);
5876 const labelList& rpList = slaveMap.removedPointList();
5877 const List<objectMap>& apList = slaveMap.addedPointList();
5879 // Configure the slave replacement point.
5880 // - collapseEdge stores this as an 'addedPoint'
5881 label scPoint = -1, srPoint = -1;
5883 if ((slaveMap.index() == eIndex) && localCouple)
5884 {
5885 // Cyclic edge at axis
5886 scPoint = collapsePoint;
5887 srPoint = replacePoint;
5888 }
5889 else
5890 {
5891 scPoint = rpList[0];
5892 srPoint = apList[0].index();
5893 }
5895 if (collapsingSlave)
5896 {
5897 if (rPointMap[replacePoint] == scPoint)
5898 {
5899 pointMap[srPoint] = replacePoint;
5900 rPointMap[replacePoint] = srPoint;
5901 }
5902 }
5903 else
5904 {
5905 if (pointMap[replacePoint] == scPoint)
5906 {
5907 rPointMap[srPoint] = replacePoint;
5908 pointMap[replacePoint] = srPoint;
5909 }
5910 }
5912 // Update face maps
5913 const label faceEnum = coupleMap::FACE;
5915 // Obtain references
5916 Map<label>& faceMap = cMap.entityMap(faceEnum);
5917 Map<label>& rFaceMap = cMap.reverseEntityMap(faceEnum);
5919 forAll(checkFaces, faceI)
5920 {
5921 label mfIndex = checkFaces[faceI];
5922 label mfPatch = whichPatch(mfIndex);
5924 const face& mF = faces_[mfIndex];
5926 triFace cF
5927 (
5928 pointMap.found(mF[0]) ? pointMap[mF[0]] : -1,
5929 pointMap.found(mF[1]) ? pointMap[mF[1]] : -1,
5930 pointMap.found(mF[2]) ? pointMap[mF[2]] : -1
5931 );
5933 if (cF[0] == -1 || cF[1] == -1 || cF[2] == -1)
5934 {
5935 writeVTK
5936 (
5937 "failedFace_"
5938 + Foam::name(mfIndex),
5939 mfIndex,
5940 2, false, true
5941 );
5943 Pout<< " Failed to configure face for: "
5944 << mfIndex << " :: " << faces_[mfIndex]
5945 << " using comparison face: " << cF
5946 << abort(FatalError);
5947 }
5949 bool matchedFace = false;
5951 // Fetch edges connected to first slave point
5952 const labelList& spEdges = pointEdges_[cF[0]];
5954 forAll(spEdges, edgeJ)
5955 {
5956 label seIndex = spEdges[edgeJ];
5958 if (whichEdgePatch(seIndex) == -1)
5959 {
5960 continue;
5961 }
5963 const labelList& seFaces = edgeFaces_[seIndex];
5965 forAll(seFaces, faceI)
5966 {
5967 label sfIndex = seFaces[faceI];
5969 if (whichPatch(sfIndex) == -1)
5970 {
5971 continue;
5972 }
5974 const face& sF = faces_[sfIndex];
5976 if
5977 (
5978 triFace::compare
5979 (
5980 triFace(sF[0], sF[1], sF[2]), cF
5981 )
5982 )
5983 {
5984 if (debug > 2)
5985 {
5986 word pN(boundaryMesh()[mfPatch].name());
5988 Pout<< " Found face: " << sfIndex
5989 << " :: " << sF
5990 << " with mfIndex: " << mfIndex
5991 << " :: " << faces_[mfIndex]
5992 << " Patch: " << pN
5993 << endl;
5994 }
5996 if (rFaceMap.found(sfIndex))
5997 {
5998 rFaceMap[sfIndex] = mfIndex;
5999 }
6000 else
6001 {
6002 rFaceMap.insert(sfIndex, mfIndex);
6003 }
6005 if (faceMap.found(mfIndex))
6006 {
6007 faceMap[mfIndex] = sfIndex;
6008 }
6009 else
6010 {
6011 faceMap.insert(mfIndex, sfIndex);
6012 }
6014 matchedFace = true;
6016 break;
6017 }
6018 }
6020 if (matchedFace)
6021 {
6022 break;
6023 }
6024 }
6026 if (!matchedFace)
6027 {
6028 writeVTK
6029 (
6030 "failedFacePoints_"
6031 + Foam::name(mfIndex),
6032 cF, 0, false, true
6033 );
6035 Pout<< " Failed to match face: "
6036 << mfIndex << " :: " << faces_[mfIndex]
6037 << " using comparison face: " << cF
6038 << abort(FatalError);
6039 }
6040 }
6042 // Update edge maps
6043 const label edgeEnum = coupleMap::EDGE;
6045 // Obtain references
6046 Map<label>& edgeMap = cMap.entityMap(edgeEnum);
6047 Map<label>& rEdgeMap = cMap.reverseEntityMap(edgeEnum);
6049 forAll(checkEdges, edgeI)
6050 {
6051 label meIndex = checkEdges[edgeI];
6053 const edge& mE = edges_[meIndex];
6055 edge cE
6056 (
6057 pointMap.found(mE[0]) ? pointMap[mE[0]] : -1,
6058 pointMap.found(mE[1]) ? pointMap[mE[1]] : -1
6059 );
6061 if (cE[0] == -1 || cE[1] == -1)
6062 {
6063 writeVTK
6064 (
6065 "failedEdge_"
6066 + Foam::name(meIndex),
6067 meIndex,
6068 1, false, true
6069 );
6071 Pout<< " Failed to configure edge for: "
6072 << meIndex << " :: " << edges_[meIndex]
6073 << " using comparison edge: " << cE
6074 << abort(FatalError);
6075 }
6077 bool matchedEdge = false;
6079 // Fetch edges connected to first slave point
6080 const labelList& spEdges = pointEdges_[cE[0]];
6082 forAll(spEdges, edgeJ)
6083 {
6084 label seIndex = spEdges[edgeJ];
6086 const edge& sE = edges_[seIndex];
6088 if (sE == cE)
6089 {
6090 if (debug > 2)
6091 {
6092 Pout<< " Found edge: " << seIndex
6093 << " :: " << sE
6094 << " with meIndex: " << meIndex
6095 << " :: " << edges_[meIndex]
6096 << endl;
6097 }
6099 // Update reverse map
6100 if (rEdgeMap.found(seIndex))
6101 {
6102 rEdgeMap[seIndex] = meIndex;
6103 }
6104 else
6105 {
6106 rEdgeMap.insert(seIndex, meIndex);
6107 }
6109 // Update map
6110 if (edgeMap.found(meIndex))
6111 {
6112 edgeMap[meIndex] = seIndex;
6113 }
6114 else
6115 {
6116 edgeMap.insert(meIndex, seIndex);
6117 }
6119 matchedEdge = true;
6121 break;
6122 }
6123 }
6125 if (!matchedEdge)
6126 {
6127 writeVTK
6128 (
6129 "failedEdgePoints_"
6130 + Foam::name(meIndex),
6131 cE, 0, false, true
6132 );
6134 Pout<< " Failed to match edge: "
6135 << meIndex << " :: " << edges_[meIndex]
6136 << " using comparison edge: " << cE
6137 << abort(FatalError);
6138 }
6139 }
6140 }
6141 else
6142 if (procCouple && !localCouple)
6143 {
6144 // Update point mapping
6145 forAll(procIndices_, pI)
6146 {
6147 const coupleMap& cMap = recvMeshes_[pI].map();
6149 // Obtain references
6150 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
6151 Map<label>& rPointMap = cMap.reverseEntityMap(coupleMap::POINT);
6153 const changeMap* slaveMapPtr = NULL;
6154 const label pointEnum = coupleMap::POINT;
6156 forAll(slaveMaps, slaveI)
6157 {
6158 const changeMap& slaveMap = slaveMaps[slaveI];
6160 if (slaveMap.patchIndex() == pI)
6161 {
6162 if (slaveMap.index() < 0)
6163 {
6164 // Point-coupling
6165 label sI = -1;
6167 if (collapsingSlave)
6168 {
6169 sI = cMap.findMaster(pointEnum, collapsePoint);
6171 if (sI > -1)
6172 {
6173 if (rPointMap.found(replacePoint))
6174 {
6175 rPointMap[replacePoint] = sI;
6176 }
6177 else
6178 {
6179 rPointMap.insert(replacePoint, sI);
6180 }
6182 pointMap[sI] = replacePoint;
6183 }
6184 }
6185 else
6186 {
6187 sI = cMap.findSlave(pointEnum, collapsePoint);
6189 if (sI > -1)
6190 {
6191 if (pointMap.found(replacePoint))
6192 {
6193 pointMap[replacePoint] = sI;
6194 }
6195 else
6196 {
6197 pointMap.insert(replacePoint, sI);
6198 }
6200 rPointMap[sI] = replacePoint;
6201 }
6202 }
6204 if (sI > -1 && debug > 2)
6205 {
6206 Pout<< " Found point: " << collapsePoint
6207 << " on proc: " << procIndices_[pI]
6208 << endl;
6209 }
6210 }
6211 else
6212 {
6213 // Edge-coupling. Fetch address for later.
6214 slaveMapPtr = &slaveMap;
6215 break;
6216 }
6217 }
6218 }
6220 if (slaveMapPtr)
6221 {
6222 // Alias for convenience...
6223 const changeMap& slaveMap = *slaveMapPtr;
6225 const labelList& rpList = slaveMap.removedPointList();
6226 const List<objectMap>& apList = slaveMap.addedPointList();
6228 // Configure the slave replacement point.
6229 // - collapseEdge stores this as an 'addedPoint'
6230 label scPoint = rpList[0];
6231 label srPoint = apList[0].index();
6233 if (collapsingSlave)
6234 {
6235 if (rPointMap[replacePoint] == scPoint)
6236 {
6237 pointMap[srPoint] = replacePoint;
6238 rPointMap[replacePoint] = srPoint;
6239 }
6241 pointMap.erase(rPointMap[collapsePoint]);
6242 rPointMap.erase(collapsePoint);
6243 }
6244 else
6245 {
6246 if (pointMap[replacePoint] == scPoint)
6247 {
6248 rPointMap[srPoint] = replacePoint;
6249 pointMap[replacePoint] = srPoint;
6250 }
6252 rPointMap.erase(pointMap[collapsePoint]);
6253 pointMap.erase(collapsePoint);
6254 }
6256 // If any other points were removed, update map
6257 for (label pointI = 1; pointI < rpList.size(); pointI++)
6258 {
6259 if (collapsingSlave)
6260 {
6261 if (pointMap.found(rpList[pointI]))
6262 {
6263 rPointMap.erase(pointMap[rpList[pointI]]);
6264 pointMap.erase(rpList[pointI]);
6265 }
6266 }
6267 else
6268 {
6269 if (rPointMap.found(rpList[pointI]))
6270 {
6271 if (debug > 2)
6272 {
6273 Pout<< " Found removed point: "
6274 << rpList[pointI]
6275 << " on proc: " << procIndices_[pI]
6276 << " for point on this proc: "
6277 << rPointMap[rpList[pointI]]
6278 << endl;
6279 }
6281 pointMap.erase(rPointMap[rpList[pointI]]);
6282 rPointMap.erase(rpList[pointI]);
6283 }
6284 }
6285 }
6286 }
6287 }
6289 // Update face mapping
6290 forAll(procIndices_, pI)
6291 {
6292 const coupleMap& cMap = recvMeshes_[pI].map();
6293 const dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
6295 // Obtain point maps
6296 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
6298 // Update face maps
6299 const label faceEnum = coupleMap::FACE;
6301 // Obtain references
6302 Map<label>& faceMap = cMap.entityMap(faceEnum);
6303 Map<label>& rFaceMap = cMap.reverseEntityMap(faceEnum);
6305 const changeMap* slaveMapPtr = NULL;
6307 forAll(slaveMaps, slaveI)
6308 {
6309 const changeMap& slaveMap = slaveMaps[slaveI];
6311 if (slaveMap.patchIndex() == pI)
6312 {
6313 if (slaveMap.index() < 0)
6314 {
6315 // Point-coupling
6316 continue;
6317 }
6319 // Edge-coupling. Fetch address for later.
6320 slaveMapPtr = &slaveMap;
6321 break;
6322 }
6323 }
6325 forAll(checkFaces, faceI)
6326 {
6327 label mfIndex = checkFaces[faceI];
6329 const face& mF = faces_[mfIndex];
6331 // Skip checks if a conversion occurred,
6332 // and this face was removed as a result
6333 if (mF.empty())
6334 {
6335 continue;
6336 }
6338 label mfPatch = whichPatch(mfIndex);
6339 label neiProc = getNeighbourProcessor(mfPatch);
6341 triFace cF
6342 (
6343 pointMap.found(mF[0]) ? pointMap[mF[0]] : -1,
6344 pointMap.found(mF[1]) ? pointMap[mF[1]] : -1,
6345 pointMap.found(mF[2]) ? pointMap[mF[2]] : -1
6346 );
6348 // Check if a patch conversion is necessary
6349 label newPatch = -1;
6350 bool requireConversion = false, physConvert = false;
6352 // Skip mapping if all points were not found
6353 if (cF[0] == -1 || cF[1] == -1 || cF[2] == -1)
6354 {
6355 // Is this actually a non-processor patch?
6356 if (neiProc == -1)
6357 {
6358 continue;
6359 }
6360 else
6361 {
6362 physConvert = true;
6363 }
6364 }
6366 // Has this face been converted to a physical boundary?
6367 if (faceMap.found(mfIndex) && slaveMapPtr)
6368 {
6369 // Alias for convenience...
6370 const changeMap& sMap = *slaveMapPtr;
6371 const labelList& rfList = sMap.removedFaceList();
6372 const List<objectMap>& afList = sMap.addedFaceList();
6374 // Was the face removed by the slave?
6375 label sfIndex = faceMap[mfIndex];
6377 if (findIndex(rfList, sfIndex) > -1)
6378 {
6379 if (debug > 2)
6380 {
6381 Pout<< " Found removed face: " << sfIndex
6382 << " with mfIndex: " << mfIndex
6383 << " :: " << faces_[mfIndex]
6384 << " on proc: " << procIndices_[pI]
6385 << endl;
6386 }
6388 // Remove map entry
6389 rFaceMap.erase(sfIndex);
6390 faceMap.erase(mfIndex);
6391 }
6393 // Check added faces for special entries
6394 forAll(afList, indexI)
6395 {
6396 const objectMap& af = afList[indexI];
6398 label mo =
6399 (
6400 af.masterObjects().size() ?
6401 af.masterObjects()[0] : 0
6402 );
6404 // Back out the physical patch ID
6405 if ((af.index() == sfIndex) && (mo < 0))
6406 {
6407 newPatch = Foam::mag(mo + 2);
6408 requireConversion = true;
6409 break;
6410 }
6411 }
6412 }
6414 // Was an appropriate physical patch found?
6415 if (physConvert && !requireConversion)
6416 {
6417 continue;
6418 }
6420 // Are we talking to a different processor?
6421 if (neiProc != procIndices_[pI])
6422 {
6423 // This face needs to be converted
6424 const polyBoundaryMesh& boundary = boundaryMesh();
6426 forAll(boundary, pi)
6427 {
6428 if (getNeighbourProcessor(pi) == procIndices_[pI])
6429 {
6430 newPatch = pi;
6431 requireConversion = true;
6432 break;
6433 }
6434 }
6435 }
6437 if (requireConversion)
6438 {
6439 // Obtain a copy before adding the new face,
6440 // since the reference might become
6441 // invalid during list resizing.
6442 face newFace = faces_[mfIndex];
6443 label newOwn = owner_[mfIndex];
6444 labelList newFaceEdges = faceEdges_[mfIndex];
6446 label newFaceIndex =
6447 (
6448 insertFace
6449 (
6450 newPatch,
6451 newFace,
6452 newOwn,
6453 -1
6454 )
6455 );
6457 faceEdges_.append(newFaceEdges);
6459 meshOps::replaceLabel
6460 (
6461 mfIndex,
6462 newFaceIndex,
6463 cells_[newOwn]
6464 );
6466 // Correct edgeFaces with the new face label.
6467 forAll(newFaceEdges, edgeI)
6468 {
6469 meshOps::replaceLabel
6470 (
6471 mfIndex,
6472 newFaceIndex,
6473 edgeFaces_[newFaceEdges[edgeI]]
6474 );
6475 }
6477 // Finally remove the face
6478 removeFace(mfIndex);
6480 // Update map
6481 map.removeFace(mfIndex);
6482 map.addFace(newFaceIndex, labelList(1, mfIndex));
6484 // Replace index and patch
6485 mfIndex = newFaceIndex;
6486 mfPatch = newPatch;
6488 // If conversion was to a physical patch,
6489 // skip the remaining face mapping steps
6490 if (getNeighbourProcessor(newPatch) == -1)
6491 {
6492 continue;
6493 }
6495 // Fail for now
6496 Pout<< " Conversion required... "
6497 << " Face: " << newFace << " : "
6498 << newFace.centre(points_)
6499 << abort(FatalError);
6501 // Push a patch-conversion operation
6502 cMap.pushOperation
6503 (
6504 newFaceIndex,
6505 coupleMap::CONVERT_PATCH,
6506 newFace.centre(points_),
6507 newFace.centre(oldPoints_)
6508 );
6509 }
6511 bool matchedFace = false;
6513 // Fetch edges connected to first slave point
6514 const labelList& spEdges = sMesh.pointEdges_[cF[0]];
6516 forAll(spEdges, edgeJ)
6517 {
6518 label seIndex = spEdges[edgeJ];
6520 if (sMesh.whichEdgePatch(seIndex) == -1)
6521 {
6522 continue;
6523 }
6525 const labelList& seFaces = sMesh.edgeFaces_[seIndex];
6527 forAll(seFaces, faceI)
6528 {
6529 label sfIndex = seFaces[faceI];
6531 if (sMesh.whichPatch(sfIndex) == -1)
6532 {
6533 continue;
6534 }
6536 const face& sF = sMesh.faces_[sfIndex];
6538 if
6539 (
6540 triFace::compare
6541 (
6542 triFace(sF[0], sF[1], sF[2]), cF
6543 )
6544 )
6545 {
6546 if (debug > 2)
6547 {
6548 word pN(boundaryMesh()[mfPatch].name());
6550 Pout<< " Found face: " << sfIndex
6551 << " :: " << sF
6552 << " with mfIndex: " << mfIndex
6553 << " :: " << faces_[mfIndex]
6554 << " on proc: " << procIndices_[pI]
6555 << " Patch: " << pN
6556 << endl;
6557 }
6559 if (rFaceMap.found(sfIndex))
6560 {
6561 rFaceMap[sfIndex] = mfIndex;
6562 }
6563 else
6564 {
6565 rFaceMap.insert(sfIndex, mfIndex);
6566 }
6568 if (faceMap.found(mfIndex))
6569 {
6570 faceMap[mfIndex] = sfIndex;
6571 }
6572 else
6573 {
6574 faceMap.insert(mfIndex, sfIndex);
6575 }
6577 matchedFace = true;
6579 break;
6580 }
6581 }
6583 if (matchedFace)
6584 {
6585 break;
6586 }
6587 }
6589 if (!matchedFace)
6590 {
6591 sMesh.writeVTK
6592 (
6593 "failedFacePoints_"
6594 + Foam::name(mfIndex),
6595 cF, 0, false, true
6596 );
6598 Pout<< " Failed to match face: "
6599 << mfIndex << " :: " << faces_[mfIndex]
6600 << " using comparison face: " << cF
6601 << " on proc: " << procIndices_[pI]
6602 << abort(FatalError);
6603 }
6604 }
6605 }
6607 // Update edge mapping
6608 forAll(procIndices_, pI)
6609 {
6610 const coupleMap& cMap = recvMeshes_[pI].map();
6611 const dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
6613 // Obtain point maps
6614 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
6616 // Update edge maps
6617 const label edgeEnum = coupleMap::EDGE;
6619 // Obtain references
6620 Map<label>& edgeMap = cMap.entityMap(edgeEnum);
6621 Map<label>& rEdgeMap = cMap.reverseEntityMap(edgeEnum);
6623 const changeMap* slaveMapPtr = NULL;
6625 forAll(slaveMaps, slaveI)
6626 {
6627 const changeMap& slaveMap = slaveMaps[slaveI];
6629 if (slaveMap.patchIndex() == pI)
6630 {
6631 if (slaveMap.index() < 0)
6632 {
6633 // Point-coupling
6634 continue;
6635 }
6637 // Edge-coupling. Fetch address for later.
6638 slaveMapPtr = &slaveMap;
6639 break;
6640 }
6641 }
6643 forAll(checkEdges, edgeI)
6644 {
6645 label meIndex = checkEdges[edgeI];
6647 const edge& mE = edges_[meIndex];
6649 // Skip checks if a conversion occurred,
6650 // and this edge was removed as a result
6651 if (edgeFaces_[meIndex].empty())
6652 {
6653 continue;
6654 }
6656 label mePatch = whichEdgePatch(meIndex);
6657 label neiProc = getNeighbourProcessor(mePatch);
6659 edge cE
6660 (
6661 pointMap.found(mE[0]) ? pointMap[mE[0]] : -1,
6662 pointMap.found(mE[1]) ? pointMap[mE[1]] : -1
6663 );
6665 // Check if a patch conversion is necessary
6666 label newPatch = -1;
6667 bool requireConversion = true, physConvert = false;
6669 // Skip mapping if all points were not found
6670 if (cE[0] == -1 || cE[1] == -1)
6671 {
6672 // Is this actually a non-processor patch?
6673 if (neiProc == -1)
6674 {
6675 continue;
6676 }
6677 else
6678 {
6679 physConvert = true;
6680 }
6681 }
6683 // Has this edge been converted to a physical boundary?
6684 const labelList& meFaces = edgeFaces_[meIndex];
6686 forAll(meFaces, faceI)
6687 {
6688 label mfPatch = whichPatch(meFaces[faceI]);
6690 if (mfPatch == -1)
6691 {
6692 continue;
6693 }
6695 if (getNeighbourProcessor(mfPatch) == -1)
6696 {
6697 // Store physical patch for conversion
6698 newPatch = mfPatch;
6699 }
6700 else
6701 {
6702 requireConversion = false;
6703 }
6704 }
6706 // Was an appropriate physical patch found?
6707 if (physConvert && !requireConversion)
6708 {
6709 continue;
6710 }
6712 if (requireConversion)
6713 {
6714 edge newEdge = edges_[meIndex];
6715 labelList newEdgeFaces = edgeFaces_[meIndex];
6717 // Insert the new edge
6718 label newEdgeIndex =
6719 (
6720 insertEdge
6721 (
6722 newPatch,
6723 newEdge,
6724 newEdgeFaces
6725 )
6726 );
6728 // Replace faceEdges for all
6729 // connected faces.
6730 forAll(newEdgeFaces, faceI)
6731 {
6732 meshOps::replaceLabel
6733 (
6734 meIndex,
6735 newEdgeIndex,
6736 faceEdges_[newEdgeFaces[faceI]]
6737 );
6738 }
6740 // Remove the edge
6741 removeEdge(meIndex);
6743 // Update map
6744 map.removeEdge(meIndex);
6745 map.addEdge(newEdgeIndex, labelList(1, meIndex));
6747 // Replace index and patch
6748 meIndex = newEdgeIndex;
6749 mePatch = newPatch;
6751 // If conversion was to a physical patch,
6752 // skip the remaining face mapping steps
6753 if (getNeighbourProcessor(newPatch) == -1)
6754 {
6755 continue;
6756 }
6757 }
6759 bool matchedEdge = false;
6761 // Fetch edges connected to first slave point
6762 const labelList& spEdges = sMesh.pointEdges_[cE[0]];
6764 forAll(spEdges, edgeJ)
6765 {
6766 label seIndex = spEdges[edgeJ];
6768 const edge& sE = sMesh.edges_[seIndex];
6770 if (sE == cE)
6771 {
6772 if (debug > 2)
6773 {
6774 Pout<< " Found edge: " << seIndex
6775 << " :: " << sE
6776 << " with meIndex: " << meIndex
6777 << " :: " << edges_[meIndex]
6778 << " on proc: " << procIndices_[pI]
6779 << endl;
6780 }
6782 // Update reverse map
6783 if (rEdgeMap.found(seIndex))
6784 {
6785 rEdgeMap[seIndex] = meIndex;
6786 }
6787 else
6788 {
6789 rEdgeMap.insert(seIndex, meIndex);
6790 }
6792 // Update map
6793 if (edgeMap.found(meIndex))
6794 {
6795 edgeMap[meIndex] = seIndex;
6796 }
6797 else
6798 {
6799 edgeMap.insert(meIndex, seIndex);
6800 }
6802 matchedEdge = true;
6804 break;
6805 }
6806 }
6808 if (!matchedEdge)
6809 {
6810 // Check if a coupling existed before
6811 if (edgeMap.found(meIndex) && slaveMapPtr)
6812 {
6813 // Alias for convenience...
6814 const changeMap& sMap = *slaveMapPtr;
6815 const labelList& reList = sMap.removedEdgeList();
6817 // Was the edge removed by the slave?
6818 if (findIndex(reList, edgeMap[meIndex]) > -1)
6819 {
6820 if (debug > 2)
6821 {
6822 Pout<< " Found removed edge: "
6823 << edgeMap[meIndex]
6824 << " with meIndex: " << meIndex
6825 << " :: " << edges_[meIndex]
6826 << " on proc: " << procIndices_[pI]
6827 << endl;
6828 }
6830 // Remove map entry
6831 rEdgeMap.erase(edgeMap[meIndex]);
6832 edgeMap.erase(meIndex);
6834 matchedEdge = true;
6835 }
6836 }
6837 }
6839 if (!matchedEdge)
6840 {
6841 sMesh.writeVTK
6842 (
6843 "failedEdgePoints_"
6844 + Foam::name(meIndex),
6845 cE, 0, false, true
6846 );
6848 Pout<< " Failed to match edge: "
6849 << meIndex << " :: " << edges_[meIndex]
6850 << " using comparison edge: " << cE
6851 << " on proc: " << procIndices_[pI]
6852 << abort(FatalError);
6853 }
6854 }
6855 }
6856 }
6857 }
6859 // Set the flag
6860 topoChangeFlag_ = true;
6862 // Increment the counter
6863 statistics_[4]++;
6865 // Increment the number of modifications
6866 statistics_[0]++;
6868 // Return a succesful collapse
6869 map.type() = collapseCase;
6871 return map;
6872 }
6875 // Remove cell layer above specified patch
6876 const changeMap dynamicTopoFvMesh::removeCellLayer
6877 (
6878 const label patchID
6879 )
6880 {
6881 changeMap map;
6883 labelHashSet edgesToRemove, facesToRemove;
6884 Map<labelPair> pointsToRemove, edgesToKeep;
6886 DynamicList<label> patchFaces(patchSizes_[patchID]);
6887 DynamicList<labelPair> cellsToRemove(patchSizes_[patchID]);
6888 DynamicList<labelPair> hFacesToRemove(patchSizes_[patchID]);
6890 for (label faceI = nOldInternalFaces_; faceI < faces_.size(); faceI++)
6891 {
6892 label pIndex = whichPatch(faceI);
6894 if (pIndex != patchID)
6895 {
6896 continue;
6897 }
6899 // Fetch owner cell
6900 label cIndex = owner_[faceI];
6902 // Add face to the list
6903 patchFaces.append(faceI);
6905 // Fetch appropriate face / cell
6906 const face& bFace = faces_[faceI];
6907 const cell& ownCell = cells_[cIndex];
6909 // Get the opposing face from the cell
6910 oppositeFace oFace = ownCell.opposingFace(faceI, faces_);
6912 if (!oFace.found())
6913 {
6914 // Something's wrong here.
6915 FatalErrorIn
6916 (
6917 "const changeMap dynamicTopoFvMesh::removeCellLayer"
6918 "(const label patchID)"
6919 )
6920 << " Face: " << faceI << " :: " << bFace << nl
6921 << " has no opposing face in cell: "
6922 << cIndex << " :: " << ownCell << nl
6923 << abort(FatalError);
6924 }
6926 // Fetch cell on the other-side of the opposite face
6927 label otherCellIndex =
6928 (
6929 (owner_[oFace.oppositeIndex()] == cIndex) ?
6930 neighbour_[oFace.oppositeIndex()] :
6931 owner_[oFace.oppositeIndex()]
6932 );
6934 if (otherCellIndex == -1)
6935 {
6936 // Opposite face is on a boundary, and layer
6937 // removal would be invalid if we continued.
6938 map.type() = -2;
6940 return map;
6941 }
6943 // Fetch reference to other cell
6944 const cell& otherCell = cells_[otherCellIndex];
6946 // Find opposite face on the other cell
6947 oppositeFace otFace =
6948 (
6949 otherCell.opposingFace
6950 (
6951 oFace.oppositeIndex(),
6952 faces_
6953 )
6954 );
6956 if (!otFace.found())
6957 {
6958 // Something's wrong here.
6959 FatalErrorIn
6960 (
6961 "const changeMap dynamicTopoFvMesh::removeCellLayer"
6962 "(const label patchID)"
6963 )
6964 << " Face: " << oFace.oppositeIndex()
6965 << " :: " << oFace << nl
6966 << " has no opposing face in cell: "
6967 << otherCellIndex << " :: " << otherCell << nl
6968 << abort(FatalError);
6969 }
6971 // All edges on the boundary face are to be retained
6972 const labelList& fEdges = faceEdges_[faceI];
6973 const labelList& ofEdges = faceEdges_[oFace.oppositeIndex()];
6974 const labelList& otfEdges = faceEdges_[otFace.oppositeIndex()];
6976 forAll(fEdges, edgeI)
6977 {
6978 label eIndex = fEdges[edgeI];
6980 if (!edgesToKeep.found(eIndex))
6981 {
6982 // Find equivalent edge on opposite face
6983 label oeIndex = -1, oteIndex = -1;
6984 const edge& bEdge = edges_[eIndex];
6986 // Build edges for comparison
6987 label startLoc = bFace.which(bEdge[0]);
6988 label endLoc = bFace.which(bEdge[1]);
6990 edge cEdge(oFace[startLoc], oFace[endLoc]);
6991 edge ctEdge(otFace[startLoc], otFace[endLoc]);
6993 forAll(ofEdges, edgeJ)
6994 {
6995 const edge& ofEdge = edges_[ofEdges[edgeJ]];
6997 if (cEdge == ofEdge)
6998 {
6999 oeIndex = ofEdges[edgeJ];
7000 break;
7001 }
7002 }
7004 forAll(otfEdges, edgeJ)
7005 {
7006 const edge& otfEdge = edges_[otfEdges[edgeJ]];
7008 if (ctEdge == otfEdge)
7009 {
7010 oteIndex = otfEdges[edgeJ];
7011 break;
7012 }
7013 }
7015 if (oeIndex < 0 || oteIndex < 0)
7016 {
7017 FatalErrorIn
7018 (
7019 "const changeMap dynamicTopoFvMesh::removeCellLayer"
7020 "(const label patchID)"
7021 )
7022 << " Could not find comparison edge: " << nl
7023 << " cEdge: " << cEdge
7024 << " oeIndex: " << oeIndex << nl
7025 << " ctEdge: " << ctEdge
7026 << " oteIndex: " << oteIndex << nl
7027 << " for edge: " << bEdge
7028 << abort(FatalError);
7029 }
7031 // Make entry
7032 edgesToKeep.insert(eIndex, labelPair(oeIndex, oteIndex));
7033 }
7034 }
7036 // Add information to removal lists
7037 cellsToRemove.append
7038 (
7039 labelPair
7040 (
7041 cIndex,
7042 otherCellIndex
7043 )
7044 );
7046 hFacesToRemove.append
7047 (
7048 labelPair
7049 (
7050 oFace.oppositeIndex(),
7051 otFace.oppositeIndex()
7052 )
7053 );
7055 // Mark points for removal
7056 forAll(oFace, pointI)
7057 {
7058 label pIndex = oFace[pointI];
7060 if (!pointsToRemove.found(pIndex))
7061 {
7062 // Make entry
7063 pointsToRemove.insert
7064 (
7065 pIndex,
7066 labelPair(bFace[pointI], otFace[pointI])
7067 );
7068 }
7069 }
7071 // Loop through cell faces and mark
7072 // faces / edges for removal
7073 forAll(ownCell, faceJ)
7074 {
7075 label fIndex = ownCell[faceJ];
7077 if (fIndex == faceI || fIndex == oFace.oppositeIndex())
7078 {
7079 continue;
7080 }
7082 if (!facesToRemove.found(fIndex))
7083 {
7084 facesToRemove.insert(fIndex);
7085 }
7087 const labelList& checkEdges = faceEdges_[fIndex];
7089 forAll(checkEdges, edgeI)
7090 {
7091 label eIndex = checkEdges[edgeI];
7093 if (edgesToKeep.found(eIndex))
7094 {
7095 continue;
7096 }
7098 if (!edgesToRemove.found(eIndex))
7099 {
7100 edgesToRemove.insert(eIndex);
7101 }
7102 }
7103 }
7104 }
7106 // Correct edgeFaces / faceEdges for retained edges
7107 forAllConstIter(Map<labelPair>, edgesToKeep, eIter)
7108 {
7109 const labelList& rmeFaces = edgeFaces_[eIter().first()];
7111 forAll(rmeFaces, faceI)
7112 {
7113 labelList& fE = faceEdges_[rmeFaces[faceI]];
7115 bool foundRp = (findIndex(fE, eIter.key()) > -1);
7116 bool foundRn = (findIndex(fE, eIter().second()) > -1);
7118 if (foundRp)
7119 {
7120 // Size-down edgeFaces for replacement
7121 meshOps::sizeDownList
7122 (
7123 rmeFaces[faceI],
7124 edgeFaces_[eIter.key()]
7125 );
7126 }
7128 if (foundRn)
7129 {
7130 // Size-up edgeFaces for replacement
7131 meshOps::sizeUpList
7132 (
7133 rmeFaces[faceI],
7134 edgeFaces_[eIter.key()]
7135 );
7137 // Replace edge index
7138 meshOps::replaceLabel
7139 (
7140 eIter().first(),
7141 eIter.key(),
7142 fE
7143 );
7144 }
7145 }
7146 }
7148 // Remove unwanted faces
7149 forAllConstIter(labelHashSet, facesToRemove, fIter)
7150 {
7151 // Remove the face
7152 removeFace(fIter.key());
7154 // Update map
7155 map.removeFace(fIter.key());
7156 }
7158 // Remove unwanted edges
7159 forAllConstIter(labelHashSet, edgesToRemove, eIter)
7160 {
7161 // Remove the edge
7162 removeEdge(eIter.key());
7164 // Update map
7165 map.removeEdge(eIter.key());
7166 }
7168 // Remove unwanted points
7169 forAllConstIter(Map<labelPair>, pointsToRemove, pIter)
7170 {
7171 // Update pointEdges information first
7172 if (!twoDMesh_)
7173 {
7174 const labelList& pEdges = pointEdges_[pIter.key()];
7176 // Configure edge for comparison
7177 edge cEdge
7178 (
7179 pIter.key(),
7180 pIter().second()
7181 );
7183 label replaceEdge = -1;
7185 forAll(pEdges, edgeI)
7186 {
7187 const edge& checkEdge = edges_[pEdges[edgeI]];
7189 if (checkEdge == cEdge)
7190 {
7191 replaceEdge = pEdges[edgeI];
7192 break;
7193 }
7194 }
7196 if (replaceEdge == -1)
7197 {
7198 FatalErrorIn
7199 (
7200 "const changeMap dynamicTopoFvMesh::removeCellLayer"
7201 "(const label patchID)"
7202 )
7203 << " Could not find comparison edge: " << nl
7204 << " cEdge: " << cEdge
7205 << " for point: " << pIter.key() << nl
7206 << " pointEdges: " << pEdges
7207 << abort(FatalError);
7208 }
7210 // Size-up pointEdges
7211 meshOps::sizeUpList
7212 (
7213 replaceEdge,
7214 pointEdges_[pIter().first()]
7215 );
7216 }
7218 // Remove the point
7219 removePoint(pIter.key());
7221 // Update map
7222 map.removePoint(pIter.key());
7223 }
7225 // Track modified faces for mapping
7226 labelHashSet modifiedFaces;
7228 // Remove all cells
7229 forAll(cellsToRemove, indexI)
7230 {
7231 // Replace face label on the other cell
7232 meshOps::replaceLabel
7233 (
7234 hFacesToRemove[indexI].first(),
7235 patchFaces[indexI],
7236 cells_[cellsToRemove[indexI].second()]
7237 );
7239 // Set owner information
7240 owner_[patchFaces[indexI]] = cellsToRemove[indexI].second();
7242 // Replace points on faces / edges
7243 const cell& otherCell = cells_[cellsToRemove[indexI].second()];
7245 forAll(otherCell, faceI)
7246 {
7247 face& faceToCheck = faces_[otherCell[faceI]];
7249 bool modified = false;
7251 forAll(faceToCheck, pointI)
7252 {
7253 if (pointsToRemove.found(faceToCheck[pointI]))
7254 {
7255 faceToCheck[pointI] =
7256 (
7257 pointsToRemove[faceToCheck[pointI]].first()
7258 );
7259 }
7261 modified = true;
7262 }
7264 // Add face to list
7265 if (modified && !modifiedFaces.found(otherCell[faceI]))
7266 {
7267 modifiedFaces.insert(otherCell[faceI]);
7268 }
7270 const labelList& fEdges = faceEdges_[otherCell[faceI]];
7272 forAll(fEdges, edgeI)
7273 {
7274 edge& edgeToCheck = edges_[fEdges[edgeI]];
7276 if (pointsToRemove.found(edgeToCheck[0]))
7277 {
7278 edgeToCheck[0] =
7279 (
7280 pointsToRemove[edgeToCheck[0]].first()
7281 );
7282 }
7284 if (pointsToRemove.found(edgeToCheck[1]))
7285 {
7286 edgeToCheck[1] =
7287 (
7288 pointsToRemove[edgeToCheck[1]].first()
7289 );
7290 }
7291 }
7292 }
7294 // Remove horizontal interior face
7295 removeFace(hFacesToRemove[indexI].first());
7297 // Update map
7298 map.removeFace(hFacesToRemove[indexI].first());
7300 // Remove the cell
7301 removeCell(cellsToRemove[indexI].first());
7303 // Update map
7304 map.removeCell(cellsToRemove[indexI].first());
7305 }
7307 // Now that all old / new cells possess correct connectivity,
7308 // compute mapping information.
7309 forAll(cellsToRemove, indexI)
7310 {
7311 // Set mapping for the modified cell
7312 setCellMapping
7313 (
7314 cellsToRemove[indexI].second(),
7315 cellsToRemove[indexI]
7316 );
7317 }
7319 // Set face mapping information for modified faces
7320 forAllConstIter(labelHashSet, modifiedFaces, fIter)
7321 {
7322 // Decide between default / weighted mapping
7323 // based on boundary information
7324 label fPatch = whichPatch(fIter.key());
7326 if (fPatch == -1)
7327 {
7328 // Default mapping for interior faces
7329 setFaceMapping(fIter.key());
7330 }
7331 else
7332 {
7333 // Map boundary face from itself
7334 setFaceMapping(fIter.key(), labelList(1, fIter.key()));
7335 }
7336 }
7338 // Set the flag
7339 topoChangeFlag_ = true;
7341 // Specify that the operation was successful
7342 map.type() = 1;
7344 // Return the changeMap
7345 return map;
7346 }
7349 // Merge a set of boundary faces into internal
7350 const changeMap dynamicTopoFvMesh::mergeBoundaryFaces
7351 (
7352 const labelList& mergeFaces
7353 )
7354 {
7355 changeMap map;
7357 if (debug > 2)
7358 {
7359 Pout << "Merging faces: " << mergeFaces << endl;
7360 }
7362 List<labelPair> mergePairs(0);
7364 forAll(mergeFaces, faceI)
7365 {
7366 const face& fI = faces_[mergeFaces[faceI]];
7368 for (label faceJ = faceI + 1; faceJ < mergeFaces.size(); faceJ++)
7369 {
7370 const face& fJ = faces_[mergeFaces[faceJ]];
7372 bool merge = false;
7374 if (fI.size() == fJ.size() && fI.size() == 3)
7375 {
7376 if
7377 (
7378 triFace::compare
7379 (
7380 triFace(fI[0], fI[1], fI[2]),
7381 triFace(fJ[0], fJ[1], fJ[2])
7382 )
7383 )
7384 {
7385 merge = true;
7386 }
7387 }
7388 else
7389 if (face::compare(fI, fJ))
7390 {
7391 merge = true;
7392 }
7394 if (merge)
7395 {
7396 meshOps::sizeUpList
7397 (
7398 labelPair(mergeFaces[faceI], mergeFaces[faceJ]),
7399 mergePairs
7400 );
7402 break;
7403 }
7404 }
7405 }
7407 if (debug > 2)
7408 {
7409 Pout<< " mergePairs: " << mergePairs << endl;
7410 }
7412 DynamicList<label> checkEdges(10);
7414 forAll(mergePairs, pairI)
7415 {
7416 label firstFace = mergePairs[pairI].first();
7417 label secondFace = mergePairs[pairI].second();
7419 // Obtain owners for both faces, and compare their labels
7420 label newOwner = -1, newNeighbour = -1;
7421 label removedFace = -1, retainedFace = -1;
7423 if (owner_[firstFace] < owner_[secondFace])
7424 {
7425 // Retain the first face
7426 removedFace = secondFace;
7427 retainedFace = firstFace;
7429 newOwner = owner_[firstFace];
7430 newNeighbour = owner_[secondFace];
7431 }
7432 else
7433 {
7434 // Retain the second face
7435 removedFace = firstFace;
7436 retainedFace = secondFace;
7438 newOwner = owner_[secondFace];
7439 newNeighbour = owner_[firstFace];
7440 }
7442 // Insert a new interior face
7443 label newFaceIndex =
7444 (
7445 insertFace
7446 (
7447 -1,
7448 face(faces_[retainedFace]),
7449 newOwner,
7450 newNeighbour
7451 )
7452 );
7454 // Add the faceEdges entry
7455 faceEdges_.append(labelList(faceEdges_[retainedFace]));
7457 // Update map
7458 map.addFace(newFaceIndex);
7460 // Replace cell with the new face label
7461 meshOps::replaceLabel
7462 (
7463 removedFace,
7464 newFaceIndex,
7465 cells_[owner_[removedFace]]
7466 );
7468 meshOps::replaceLabel
7469 (
7470 retainedFace,
7471 newFaceIndex,
7472 cells_[owner_[retainedFace]]
7473 );
7475 const labelList& fEdges = faceEdges_[newFaceIndex];
7476 const labelList& rfEdges = faceEdges_[removedFace];
7478 // Check for common edges on the removed face
7479 forAll(rfEdges, edgeI)
7480 {
7481 label reIndex = rfEdges[edgeI];
7483 if (findIndex(fEdges, reIndex) == -1)
7484 {
7485 // Find the equivalent edge
7486 const edge& rEdge = edges_[reIndex];
7487 const labelList& reFaces = edgeFaces_[reIndex];
7489 label keIndex = -1;
7491 forAll(fEdges, edgeJ)
7492 {
7493 if (edges_[fEdges[edgeJ]] == rEdge)
7494 {
7495 keIndex = fEdges[edgeJ];
7496 break;
7497 }
7498 }
7500 if (keIndex == -1)
7501 {
7502 Pout<< " Could not find edge for: "
7503 << reIndex << " :: " << rEdge
7504 << abort(FatalError);
7505 }
7507 // Add edgeFaces from this edge
7508 forAll(reFaces, faceI)
7509 {
7510 if (reFaces[faceI] == removedFace)
7511 {
7512 continue;
7513 }
7515 meshOps::sizeUpList
7516 (
7517 reFaces[faceI],
7518 edgeFaces_[keIndex]
7519 );
7521 meshOps::replaceLabel
7522 (
7523 reIndex,
7524 keIndex,
7525 faceEdges_[reFaces[faceI]]
7526 );
7527 }
7529 // Remove the old edge
7530 removeEdge(reIndex);
7532 // Update map
7533 map.removeEdge(reIndex);
7534 }
7535 }
7537 // Replace edgeFaces with the new face label
7538 forAll(fEdges, edgeI)
7539 {
7540 label eIndex = fEdges[edgeI];
7542 if (findIndex(edgeFaces_[eIndex], removedFace) > -1)
7543 {
7544 meshOps::sizeDownList
7545 (
7546 removedFace,
7547 edgeFaces_[eIndex]
7548 );
7549 }
7551 if (findIndex(edgeFaces_[eIndex], retainedFace) > -1)
7552 {
7553 meshOps::sizeDownList
7554 (
7555 retainedFace,
7556 edgeFaces_[eIndex]
7557 );
7558 }
7560 // Size-up list with the new face index
7561 meshOps::sizeUpList
7562 (
7563 newFaceIndex,
7564 edgeFaces_[eIndex]
7565 );
7567 if (findIndex(checkEdges, eIndex) == -1)
7568 {
7569 checkEdges.append(eIndex);
7570 }
7571 }
7573 // Remove the boundary faces
7574 removeFace(removedFace);
7575 removeFace(retainedFace);
7577 // Update map
7578 map.removeFace(removedFace);
7579 map.removeFace(retainedFace);
7580 }
7582 forAll(checkEdges, edgeI)
7583 {
7584 bool allInterior = true;
7585 label eIndex = checkEdges[edgeI];
7587 const labelList& eFaces = edgeFaces_[eIndex];
7589 forAll(eFaces, faceI)
7590 {
7591 if (whichPatch(eFaces[faceI]) > -1)
7592 {
7593 allInterior = false;
7594 break;
7595 }
7596 }
7598 if (allInterior)
7599 {
7600 // This edge needs to be converted to an interior one
7601 label newEdgeIndex =
7602 (
7603 insertEdge
7604 (
7605 -1,
7606 edge(edges_[eIndex]),
7607 labelList(edgeFaces_[eIndex])
7608 )
7609 );
7611 // Update map
7612 map.addEdge(newEdgeIndex, labelList(1, eIndex));
7614 // Update faceEdges information for all connected faces
7615 const labelList& neFaces = edgeFaces_[newEdgeIndex];
7617 forAll(neFaces, faceI)
7618 {
7619 meshOps::replaceLabel
7620 (
7621 eIndex,
7622 newEdgeIndex,
7623 faceEdges_[neFaces[faceI]]
7624 );
7625 }
7627 // Remove the old boundary edge
7628 removeEdge(eIndex);
7630 // Update map
7631 map.removeEdge(eIndex);
7633 // Replace index
7634 checkEdges[edgeI] = newEdgeIndex;
7635 }
7636 }
7638 if (debug > 2)
7639 {
7640 Pout<< "Merge complete." << nl << endl;
7641 }
7643 // Return a succesful merge
7644 map.type() = 1;
7646 return map;
7647 }
7650 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
7652 } // End namespace Foam
7654 // ************************************************************************* //