| blob | 1c90538f42543f73c9399d2f4d23a0bccb466d56 |
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 for the bisection of a quad-face in 2D
44 // - Returns a changeMap with a type specifying:
45 // 1: Bisection was successful
46 // -1: Bisection failed since max number of topo-changes was reached.
47 // -2: Bisection failed since resulting quality would be unacceptable.
48 // -3: Bisection failed since edge was on a noRefinement patch.
49 const changeMap dynamicTopoFvMesh::bisectQuadFace
50 (
51 const label fIndex,
52 const changeMap& masterMap,
53 bool checkOnly,
54 bool forceOp
55 )
56 {
57 // Quad-face bisection performs the following operations:
58 // [1] Add two points at the middle of the face
59 // [2] Create a new internal face for each bisected cell
60 // [3] Modify existing face and create a new half-face
61 // [4] Modify triangular boundary faces, and create new ones as well
62 // [5] Create edges for new faces
63 // Update faceEdges and edgeFaces information
65 // Figure out which thread this is...
66 label tIndex = self();
68 // Prepare the changeMaps
69 changeMap map;
70 List<changeMap> slaveMaps;
71 bool bisectingSlave = 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 dynamicTopoFvMesh::bisectQuadFace\n"
100 "(\n"
101 " const label fIndex,\n"
102 " const changeMap& masterMap,\n"
103 " bool checkOnly,\n"
104 " bool forceOp\n"
105 ")\n"
106 )
107 << " Invalid index: " << fIndex << nl
108 << " nFaces: " << nFaces_
109 << abort(FatalError);
110 }
112 bool found;
113 label replaceFace = -1, retainFace = -1;
114 face tmpQuadFace(4), tmpTriFace(3);
115 labelList tmpQFEdges(4, -1), tmpTFEdges(3, -1);
116 FixedList<label,7> newFaceIndex(-1), newEdgeIndex(-1);
117 FixedList<edge,4> commonEdges(edge(-1, -1));
118 FixedList<label,4> cornerEdgeIndex(-1), commonEdgeIndex(-1);
119 FixedList<label,4> commonFaceIndex(-1);
120 FixedList<label,2> newPointIndex(-1), newCellIndex(-1);
121 FixedList<label,4> otherEdgeIndex(-1), otherEdgePoint(-1);
122 FixedList<label,4> otherPointIndex(-1), nextToOtherPoint(-1);
123 FixedList<label,2> c0BdyIndex(-1), c0IntIndex(-1);
124 FixedList<label,2> c1BdyIndex(-1), c1IntIndex(-1);
125 FixedList<face,2> c0BdyFace(face(3)), c0IntFace(face(4));
126 FixedList<face,2> c1BdyFace(face(3)), c1IntFace(face(4));
128 // Get the two cells on either side...
129 label c0 = owner_[fIndex], c1 = neighbour_[fIndex];
131 // Keep track of old / new cells
132 FixedList<cell, 2> oldCells(cell(5));
133 FixedList<cell, 2> newCells(cell(5));
135 // Find the prism faces for cell[0].
136 oldCells[0] = cells_[c0];
138 meshOps::findPrismFaces
139 (
140 fIndex,
141 c0,
142 faces_,
143 cells_,
144 neighbour_,
145 c0BdyFace,
146 c0BdyIndex,
147 c0IntFace,
148 c0IntIndex
149 );
151 // Check for resulting quality
152 if (checkBisection(fIndex, c0BdyIndex[0], forceOp))
153 {
154 map.type() = -2;
155 return map;
156 }
158 if (c1 != -1)
159 {
160 // Find the prism faces for cell[1].
161 meshOps::findPrismFaces
162 (
163 fIndex,
164 c1,
165 faces_,
166 cells_,
167 neighbour_,
168 c1BdyFace,
169 c1BdyIndex,
170 c1IntFace,
171 c1IntIndex
172 );
174 // Check for resulting quality
175 if (checkBisection(fIndex, c1BdyIndex[0], forceOp))
176 {
177 map.type() = -2;
178 return map;
179 }
180 }
182 // Find the common-edge between the triangular boundary faces
183 // and the face under consideration.
184 meshOps::findCommonEdge
185 (
186 c0BdyIndex[0],
187 fIndex,
188 faceEdges_,
189 commonEdgeIndex[0]
190 );
192 meshOps::findCommonEdge
193 (
194 c0BdyIndex[1],
195 fIndex,
196 faceEdges_,
197 commonEdgeIndex[1]
198 );
200 commonEdges[0] = edges_[commonEdgeIndex[0]];
201 commonEdges[1] = edges_[commonEdgeIndex[1]];
203 // If coupled modification is set, and this is a
204 // master edge, bisect its slaves as well.
205 bool localCouple = false, procCouple = false;
207 if (coupledModification_)
208 {
209 const label faceEnum = coupleMap::FACE;
210 const label pointEnum = coupleMap::POINT;
212 // Is this a locally coupled face (either master or slave)?
213 if (locallyCoupledEntity(fIndex, true))
214 {
215 localCouple = true;
216 procCouple = false;
217 }
218 else
219 if (processorCoupledEntity(fIndex))
220 {
221 procCouple = true;
222 localCouple = false;
223 }
225 if (localCouple && !procCouple)
226 {
227 // Determine the slave index.
228 label sIndex = -1, pIndex = -1;
230 forAll(patchCoupling_, patchI)
231 {
232 if (patchCoupling_(patchI))
233 {
234 const coupleMap& cMap = patchCoupling_[patchI].map();
236 if ((sIndex = cMap.findSlave(faceEnum, fIndex)) > -1)
237 {
238 pIndex = patchI;
240 break;
241 }
243 // The following bit happens only during the sliver
244 // exudation process, since slave faces are
245 // usually not added to the coupled face-stack.
246 if ((sIndex = cMap.findMaster(faceEnum, fIndex)) > -1)
247 {
248 pIndex = patchI;
250 // Notice that we are collapsing a slave face first.
251 bisectingSlave = true;
253 break;
254 }
255 }
256 }
258 if (sIndex == -1)
259 {
260 FatalErrorIn
261 (
262 "\n"
263 "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
264 "(\n"
265 " const label fIndex,\n"
266 " const changeMap& masterMap,\n"
267 " bool checkOnly,\n"
268 " bool forceOp"
269 ")\n"
270 )
271 << "Coupled maps were improperly specified." << nl
272 << " Slave index not found for: " << nl
273 << " Face: " << fIndex << ": " << faces_[fIndex] << nl
274 << abort(FatalError);
275 }
276 else
277 {
278 // If we've found the slave, size up the list
279 meshOps::sizeUpList
280 (
281 changeMap(),
282 slaveMaps
283 );
285 // Save index and patch for posterity
286 slaveMaps[0].index() = sIndex;
287 slaveMaps[0].patchIndex() = pIndex;
288 }
290 if (debug > 1)
291 {
292 Pout<< nl << " >> Bisecting slave face: " << sIndex
293 << " for master face: " << fIndex << endl;
294 }
295 }
296 else
297 if (procCouple && !localCouple)
298 {
299 // If this is a new entity, bail out for now.
300 // This will be handled at the next time-step.
301 if (fIndex >= nOldFaces_)
302 {
303 return map;
304 }
306 // Check slaves
307 forAll(procIndices_, pI)
308 {
309 // Fetch reference to subMesh
310 const coupledInfo& recvMesh = recvMeshes_[pI];
311 const coupleMap& cMap = recvMesh.map();
313 label sIndex = -1;
315 if ((sIndex = cMap.findSlave(faceEnum, fIndex)) > -1)
316 {
317 if (debug > 3)
318 {
319 Pout<< "Checking slave face: " << sIndex
320 << " on proc: " << procIndices_[pI]
321 << " for master face: " << fIndex
322 << endl;
323 }
325 // Check if a lower-ranked processor is
326 // handling this edge
327 if (procIndices_[pI] < Pstream::myProcNo())
328 {
329 if (debug > 3)
330 {
331 Pout<< "Face: " << fIndex
332 << " is handled by proc: "
333 << procIndices_[pI]
334 << ", so bailing out."
335 << endl;
336 }
338 return map;
339 }
341 label curIndex = slaveMaps.size();
343 // Size up the list
344 meshOps::sizeUpList
345 (
346 changeMap(),
347 slaveMaps
348 );
350 // Save index and patch for posterity
351 slaveMaps[curIndex].index() = sIndex;
352 slaveMaps[curIndex].patchIndex() = pI;
354 // Only one slave coupling is possible, so bail out
355 break;
356 }
357 }
358 }
359 else
360 {
361 // Something's wrong with coupling maps
362 FatalErrorIn
363 (
364 "\n"
365 "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
366 "(\n"
367 " const label fIndex,\n"
368 " const changeMap& masterMap,\n"
369 " bool checkOnly,\n"
370 " bool forceOp"
371 ")\n"
372 )
373 << "Coupled maps were improperly specified." << nl
374 << " localCouple: " << localCouple << nl
375 << " procCouple: " << procCouple << nl
376 << " Face: " << fIndex << ": " << faces_[fIndex] << nl
377 << abort(FatalError);
378 }
380 // Alias for convenience...
381 changeMap& slaveMap = slaveMaps[0];
383 label sIndex = slaveMap.index();
384 label pI = slaveMap.patchIndex();
385 const coupleMap* cMapPtr = NULL;
387 // Temporarily turn off coupledModification
388 unsetCoupledModification();
390 if (localCouple)
391 {
392 cMapPtr = &(patchCoupling_[pI].map());
394 // First check the slave for bisection feasibility.
395 slaveMap = bisectQuadFace(sIndex, changeMap(), true, forceOp);
396 }
397 else
398 if (procCouple)
399 {
400 cMapPtr = &(recvMeshes_[pI].map());
402 coupledInfo& recvMesh = recvMeshes_[pI];
404 // First check the slave for bisection feasibility.
405 slaveMap =
406 (
407 recvMesh.subMesh().bisectQuadFace
408 (
409 sIndex,
410 changeMap(),
411 true,
412 forceOp
413 )
414 );
415 }
417 // Turn it back on.
418 setCoupledModification();
420 if (slaveMap.type() <= 0)
421 {
422 // Slave couldn't perform bisection.
423 map.type() = -2;
425 return map;
426 }
428 // Save index and patch for posterity
429 slaveMap.index() = sIndex;
430 slaveMap.patchIndex() = pI;
432 // Alias for convenience..
433 const coupleMap& cMap = *cMapPtr;
435 // Temporarily turn off coupledModification
436 unsetCoupledModification();
438 // First check the master for bisection feasibility.
439 changeMap masterMap = bisectQuadFace(fIndex, changeMap(), true);
441 // Turn it back on.
442 setCoupledModification();
444 // Master couldn't perform bisection
445 if (masterMap.type() <= 0)
446 {
447 return masterMap;
448 }
450 // Fill the masterMap with points that
451 // we seek maps for...
452 FixedList<labelList, 2> slaveEdges(labelList(2, -1));
454 forAll(slaveEdges, edgeI)
455 {
456 labelList& slaveEdge = slaveEdges[edgeI];
458 // Renumber to slave indices
459 forAll(slaveEdge, pointI)
460 {
461 slaveEdge[pointI] =
462 (
463 cMap.findSlave
464 (
465 pointEnum,
466 commonEdges[edgeI][pointI]
467 )
468 );
469 }
471 masterMap.addPoint(-1, slaveEdge);
472 }
474 // Temporarily turn off coupledModification
475 unsetCoupledModification();
477 if (localCouple)
478 {
479 // Bisect the local slave.
480 slaveMap = bisectQuadFace(sIndex, masterMap, false, forceOp);
481 }
482 else
483 {
484 coupledInfo& recvMesh = recvMeshes_[pI];
486 // Bisect the slave face
487 slaveMap =
488 (
489 recvMesh.subMesh().bisectQuadFace
490 (
491 sIndex,
492 masterMap,
493 false,
494 forceOp
495 )
496 );
497 }
499 // Turn coupledModification back on.
500 setCoupledModification();
502 // The final operation has to succeed.
503 if (slaveMap.type() <= 0)
504 {
505 FatalErrorIn
506 (
507 "\n"
508 "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
509 "(\n"
510 " const label fIndex,\n"
511 " const changeMap& masterMap,\n"
512 " bool checkOnly,\n"
513 " bool forceOp"
514 ")\n"
515 )
516 << "Coupled topo-change for slave failed."
517 << " Master face: " << fIndex << nl
518 << " Slave face: " << sIndex << nl
519 << " Patch index: " << pI << nl
520 << " Type: " << slaveMap.type() << nl
521 << abort(FatalError);
522 }
524 // Save index and patch for posterity
525 slaveMap.index() = sIndex;
526 slaveMap.patchIndex() = pI;
527 }
529 // Are we performing only checks?
530 if (checkOnly)
531 {
532 map.type() = 1;
533 return map;
534 }
536 if (debug > 1)
537 {
538 Pout<< nl << nl
539 << "Face: " << fIndex
540 << ": " << faces_[fIndex]
541 << " is to be bisected. " << endl;
543 if (debug > 2)
544 {
545 Pout<< " On SubMesh: " << Switch::asText(isSubMesh_) << nl;
546 Pout<< " coupledModification: " << coupledModification_ << nl;
548 const polyBoundaryMesh& boundary = boundaryMesh();
550 label epIndex = whichPatch(fIndex);
552 Pout<< " Patch: ";
554 if (epIndex == -1)
555 {
556 Pout<< "Internal" << endl;
557 }
558 else
559 if (epIndex < boundary.size())
560 {
561 Pout<< boundary[epIndex].name() << endl;
562 }
563 else
564 {
565 Pout<< " New patch: " << epIndex << endl;
566 }
568 Pout<< "Cell[0]: " << c0 << ": " << oldCells[0] << endl;
570 forAll(oldCells[0], faceI)
571 {
572 const labelList& fE = faceEdges_[oldCells[0][faceI]];
574 Pout<< oldCells[0][faceI] << ": "
575 << faces_[oldCells[0][faceI]]
576 << " fE: " << fE
577 << endl;
579 forAll(fE, edgeI)
580 {
581 const labelList& eF = edgeFaces_[fE[edgeI]];
583 Pout<< '\t' << fE[edgeI]
584 << ": " << edges_[fE[edgeI]]
585 << " eF: " << eF
586 << endl;
587 }
588 }
589 }
591 // Write out VTK files prior to change
592 if (debug > 3)
593 {
594 labelList cellHull(2, -1);
596 cellHull[0] = owner_[fIndex];
597 cellHull[1] = neighbour_[fIndex];
599 writeVTK
600 (
601 Foam::name(fIndex)
602 + "_Bisect_0",
603 cellHull
604 );
605 }
606 }
608 // Find the isolated point on both boundary faces of cell[0]
609 meshOps::findIsolatedPoint
610 (
611 c0BdyFace[0],
612 commonEdges[0],
613 otherPointIndex[0],
614 nextToOtherPoint[0]
615 );
617 meshOps::findIsolatedPoint
618 (
619 c0BdyFace[1],
620 commonEdges[1],
621 otherPointIndex[1],
622 nextToOtherPoint[1]
623 );
625 // For convenience...
626 otherEdgePoint[0] = commonEdges[0].otherVertex(nextToOtherPoint[0]);
627 otherEdgePoint[1] = commonEdges[1].otherVertex(nextToOtherPoint[1]);
629 labelList mP(2, -1);
631 // Set mapping for this point
632 mP[0] = commonEdges[0][0];
633 mP[1] = commonEdges[0][1];
635 // Add two new points to the end of the list
636 newPointIndex[0] =
637 (
638 insertPoint
639 (
640 0.5 * (points_[mP[0]] + points_[mP[1]]),
641 0.5 * (oldPoints_[mP[0]] + oldPoints_[mP[1]]),
642 mP
643 )
644 );
646 // Set mapping for this point
647 mP[0] = commonEdges[1][0];
648 mP[1] = commonEdges[1][1];
650 newPointIndex[1] =
651 (
652 insertPoint
653 (
654 0.5 * (points_[mP[0]] + points_[mP[1]]),
655 0.5 * (oldPoints_[mP[0]] + oldPoints_[mP[1]]),
656 mP
657 )
658 );
660 // Add the points to the map. Since this might require master mapping,
661 // first check to see if a slave is being bisected.
662 if (masterMap.addedPointList().size())
663 {
664 const List<objectMap>& pMap =
665 (
666 masterMap.addedPointList()
667 );
669 edge firstEdge
670 (
671 pMap[0].masterObjects()[0],
672 pMap[0].masterObjects()[1]
673 );
675 edge secondEdge
676 (
677 pMap[1].masterObjects()[0],
678 pMap[1].masterObjects()[1]
679 );
681 if (firstEdge == commonEdges[0] && secondEdge == commonEdges[1])
682 {
683 // Update in conventional order
684 map.addPoint(newPointIndex[0]);
685 map.addPoint(newPointIndex[1]);
686 }
687 else
688 if (firstEdge == commonEdges[1] && secondEdge == commonEdges[0])
689 {
690 // Update in reverse order
691 map.addPoint(newPointIndex[1]);
692 map.addPoint(newPointIndex[0]);
693 }
694 else
695 {
696 // We have a problem
697 FatalErrorIn
698 (
699 "\n"
700 "const changeMap "
701 "dynamicTopoFvMesh::bisectQuadFace\n"
702 "(\n"
703 " const label fIndex,\n"
704 " const changeMap& masterMap,\n"
705 " bool checkOnly,\n"
706 " bool forceOp\n"
707 ")\n"
708 )
709 << "Coupled topo-change for slave failed."
710 << " firstEdge: " << firstEdge << nl
711 << " secondEdge: " << secondEdge << nl
712 << " commonEdges[0]: " << commonEdges[0] << nl
713 << " commonEdges[1]: " << commonEdges[1] << nl
714 << abort(FatalError);
715 }
716 }
717 else
718 {
719 map.addPoint(newPointIndex[0]);
720 map.addPoint(newPointIndex[1]);
721 }
723 // Add a new prism cell to the end of the list.
724 // Currently invalid, but will be updated later.
725 newCellIndex[0] = insertCell(newCells[0], lengthScale_[c0]);
727 // Add this cell to the map.
728 map.addCell(newCellIndex[0]);
730 // Modify the two existing triangle boundary faces
732 // Zeroth boundary face - Owner = c[0] & Neighbour [-1] (unchanged)
733 meshOps::replaceLabel
734 (
735 otherEdgePoint[0],
736 newPointIndex[0],
737 c0BdyFace[0]
738 );
740 // First boundary face - Owner = newCell[0], Neighbour = -1
741 meshOps::replaceLabel
742 (
743 otherEdgePoint[1],
744 newPointIndex[1],
745 c0BdyFace[1]
746 );
748 // Update faces.
749 faces_[c0BdyIndex[0]] = c0BdyFace[0];
750 faces_[c0BdyIndex[1]] = c0BdyFace[1];
752 owner_[c0BdyIndex[1]] = newCellIndex[0];
754 meshOps::replaceLabel
755 (
756 c0BdyIndex[1],
757 -1,
758 oldCells[0]
759 );
761 // Detect edges other than commonEdges
762 const labelList& fEdges = faceEdges_[fIndex];
764 forAll(fEdges, edgeI)
765 {
766 if
767 (
768 fEdges[edgeI] != commonEdgeIndex[0] &&
769 fEdges[edgeI] != commonEdgeIndex[1]
770 )
771 {
772 // Obtain a reference to this edge
773 const edge& eThis = edges_[fEdges[edgeI]];
775 if
776 (
777 eThis[0] == nextToOtherPoint[0]
778 || eThis[1] == nextToOtherPoint[0]
779 )
780 {
781 otherEdgeIndex[0] = fEdges[edgeI];
782 }
783 else
784 {
785 otherEdgeIndex[1] = fEdges[edgeI];
786 }
787 }
788 }
790 // Modify point-labels on the quad face under consideration
791 meshOps::replaceLabel
792 (
793 otherEdgePoint[0],
794 newPointIndex[0],
795 faces_[fIndex]
796 );
798 meshOps::replaceLabel
799 (
800 nextToOtherPoint[1],
801 newPointIndex[1],
802 faces_[fIndex]
803 );
805 // Add this face to the map.
806 // Although this face isn't technically 'added', it's
807 // required for coupled patch mapping.
808 map.addFace(fIndex);
810 if (debug > 1)
811 {
812 Pout<< "Modified face: " << fIndex
813 << ": " << faces_[fIndex] << endl;
815 if (debug > 2)
816 {
817 Pout<< "Common edges: " << nl
818 << commonEdgeIndex[0] << ": " << commonEdges[0] << nl
819 << commonEdgeIndex[1] << ": " << commonEdges[1]
820 << endl;
821 }
822 }
824 // Find the quad face that contains otherEdgeIndex[1]
825 found = false;
827 const labelList& e1 = faceEdges_[c0IntIndex[0]];
829 forAll(e1, edgeI)
830 {
831 if (e1[edgeI] == otherEdgeIndex[1])
832 {
833 meshOps::replaceLabel
834 (
835 c0IntIndex[0],
836 -1,
837 oldCells[0]
838 );
840 replaceFace = c0IntIndex[0];
841 retainFace = c0IntIndex[1];
842 found = true;
843 break;
844 }
845 }
847 if (!found)
848 {
849 // The edge was not found before
850 meshOps::replaceLabel
851 (
852 c0IntIndex[1],
853 -1,
854 oldCells[0]
855 );
857 replaceFace = c0IntIndex[1];
858 retainFace = c0IntIndex[0];
859 }
861 // Check if face reversal is necessary for the replacement
862 if (owner_[replaceFace] == c0)
863 {
864 if (neighbour_[replaceFace] == -1)
865 {
866 // Change the owner
867 owner_[replaceFace] = newCellIndex[0];
868 }
869 else
870 {
871 // This face has to be reversed
872 faces_[replaceFace] = faces_[replaceFace].reverseFace();
873 owner_[replaceFace] = neighbour_[replaceFace];
874 neighbour_[replaceFace] = newCellIndex[0];
876 setFlip(replaceFace);
877 }
878 }
879 else
880 {
881 // Keep owner, but change neighbour
882 neighbour_[replaceFace] = newCellIndex[0];
883 }
885 // Define the faces for the new cell
886 newCells[0][0] = c0BdyIndex[1];
887 newCells[0][1] = replaceFace;
889 // Define the set of new faces and insert...
891 // New interior face; Owner = cell[0] & Neighbour = newCell[0]
892 tmpQuadFace[0] = otherPointIndex[0];
893 tmpQuadFace[1] = newPointIndex[0];
894 tmpQuadFace[2] = newPointIndex[1];
895 tmpQuadFace[3] = otherPointIndex[1];
897 newFaceIndex[0] =
898 (
899 insertFace
900 (
901 -1,
902 tmpQuadFace,
903 c0,
904 newCellIndex[0]
905 )
906 );
908 // Add a faceEdges entry as well
909 faceEdges_.append(tmpQFEdges);
911 // Add this face to the map.
912 map.addFace(newFaceIndex[0]);
914 // Find the common edge between quad/quad faces...
915 meshOps::findCommonEdge
916 (
917 c0IntIndex[0],
918 c0IntIndex[1],
919 faceEdges_,
920 otherEdgeIndex[2]
921 );
923 // ... and size-up edgeFaces for the edge
924 meshOps::sizeUpList
925 (
926 newFaceIndex[0],
927 edgeFaces_[otherEdgeIndex[2]]
928 );
930 meshOps::replaceLabel
931 (
932 -1,
933 newFaceIndex[0],
934 newCells[0]
935 );
937 meshOps::replaceLabel
938 (
939 -1,
940 newFaceIndex[0],
941 oldCells[0]
942 );
944 // Second boundary face; Owner = newCell[0] & Neighbour = [-1]
945 tmpTriFace[0] = otherPointIndex[0];
946 tmpTriFace[1] = newPointIndex[0];
947 tmpTriFace[2] = otherEdgePoint[0];
949 newFaceIndex[1] =
950 (
951 insertFace
952 (
953 whichPatch(c0BdyIndex[0]),
954 tmpTriFace,
955 newCellIndex[0],
956 -1
957 )
958 );
960 // Add a faceEdges entry as well
961 faceEdges_.append(tmpTFEdges);
963 // Add this face to the map.
964 map.addFace(newFaceIndex[1]);
966 meshOps::replaceLabel
967 (
968 -1,
969 newFaceIndex[1],
970 newCells[0]
971 );
973 // Third boundary face; Owner = c[0] & Neighbour = [-1]
974 tmpTriFace[0] = otherPointIndex[1];
975 tmpTriFace[1] = newPointIndex[1];
976 tmpTriFace[2] = otherEdgePoint[1];
978 newFaceIndex[2] =
979 (
980 insertFace
981 (
982 whichPatch(c0BdyIndex[1]),
983 tmpTriFace,
984 c0,
985 -1
986 )
987 );
989 // Add a faceEdges entry as well
990 faceEdges_.append(tmpTFEdges);
992 // Add this face to the map.
993 map.addFace(newFaceIndex[2]);
995 meshOps::replaceLabel
996 (
997 -1,
998 newFaceIndex[2],
999 oldCells[0]
1000 );
1002 // Create / modify edges...
1003 labelList tmpTriEdgeFaces(3, -1);
1005 // The edge bisecting the zeroth boundary triangular face
1006 tmpTriEdgeFaces[0] = c0BdyIndex[0];
1007 tmpTriEdgeFaces[1] = newFaceIndex[0];
1008 tmpTriEdgeFaces[2] = newFaceIndex[1];
1010 newEdgeIndex[1] =
1011 (
1012 insertEdge
1013 (
1014 whichEdgePatch(commonEdgeIndex[0]),
1015 edge(newPointIndex[0], otherPointIndex[0]),
1016 tmpTriEdgeFaces
1017 )
1018 );
1020 // Add this edge to the map.
1021 map.addEdge(newEdgeIndex[1]);
1023 // Find the common edge between the quad/tri faces...
1024 meshOps::findCommonEdge
1025 (
1026 c0BdyIndex[0],
1027 replaceFace,
1028 faceEdges_,
1029 cornerEdgeIndex[0]
1030 );
1032 // ...and correct faceEdges / edgeFaces
1033 meshOps::replaceLabel
1034 (
1035 cornerEdgeIndex[0],
1036 newEdgeIndex[1],
1037 faceEdges_[c0BdyIndex[0]]
1038 );
1040 meshOps::replaceLabel
1041 (
1042 c0BdyIndex[0],
1043 newFaceIndex[1],
1044 edgeFaces_[cornerEdgeIndex[0]]
1045 );
1047 // The edge bisecting the first boundary triangular face
1048 tmpTriEdgeFaces[0] = c0BdyIndex[1];
1049 tmpTriEdgeFaces[1] = newFaceIndex[0];
1050 tmpTriEdgeFaces[2] = newFaceIndex[2];
1052 newEdgeIndex[2] =
1053 (
1054 insertEdge
1055 (
1056 whichEdgePatch(commonEdgeIndex[1]),
1057 edge(newPointIndex[1], otherPointIndex[1]),
1058 tmpTriEdgeFaces
1059 )
1060 );
1062 // Add this edge to the map.
1063 map.addEdge(newEdgeIndex[2]);
1065 // Find the common edge between the quad/tri faces...
1066 meshOps::findCommonEdge
1067 (
1068 c0BdyIndex[1],
1069 retainFace,
1070 faceEdges_,
1071 cornerEdgeIndex[1]
1072 );
1074 // ...and correct faceEdges / edgeFaces
1075 meshOps::replaceLabel
1076 (
1077 cornerEdgeIndex[1],
1078 newEdgeIndex[2],
1079 faceEdges_[c0BdyIndex[1]]
1080 );
1082 meshOps::replaceLabel
1083 (
1084 c0BdyIndex[1],
1085 newFaceIndex[2],
1086 edgeFaces_[cornerEdgeIndex[1]]
1087 );
1089 if (c1 == -1)
1090 {
1091 // The quad boundary face resulting from bisection;
1092 // Owner = newCell[0] & Neighbour = [-1]
1093 tmpQuadFace[0] = newPointIndex[1];
1094 tmpQuadFace[1] = nextToOtherPoint[1];
1095 tmpQuadFace[2] = otherEdgePoint[0];
1096 tmpQuadFace[3] = newPointIndex[0];
1098 newFaceIndex[3] =
1099 (
1100 insertFace
1101 (
1102 whichPatch(fIndex),
1103 tmpQuadFace,
1104 newCellIndex[0],
1105 -1
1106 )
1107 );
1109 // Add a faceEdges entry as well
1110 faceEdges_.append(tmpQFEdges);
1112 // Add this face to the map.
1113 map.addFace(newFaceIndex[3]);
1115 // Correct edgeFaces for otherEdgeIndex[1]
1116 meshOps::replaceLabel
1117 (
1118 fIndex,
1119 newFaceIndex[3],
1120 edgeFaces_[otherEdgeIndex[1]]
1121 );
1123 meshOps::replaceLabel
1124 (
1125 -1,
1126 newFaceIndex[3],
1127 newCells[0]
1128 );
1130 labelList tmpBiEdgeFaces(2, -1);
1132 // The edge bisecting the face
1133 tmpTriEdgeFaces[0] = newFaceIndex[3];
1134 tmpTriEdgeFaces[1] = newFaceIndex[0];
1135 tmpTriEdgeFaces[2] = fIndex;
1137 newEdgeIndex[0] =
1138 (
1139 insertEdge
1140 (
1141 whichEdgePatch(otherEdgeIndex[0]),
1142 edge(newPointIndex[0], newPointIndex[1]),
1143 tmpTriEdgeFaces
1144 )
1145 );
1147 // Add this edge to the map.
1148 map.addEdge(newEdgeIndex[0]);
1150 // Replace an edge on the bisected face
1151 meshOps::replaceLabel
1152 (
1153 otherEdgeIndex[1],
1154 newEdgeIndex[0],
1155 faceEdges_[fIndex]
1156 );
1158 // Create / replace side edges created from face bisection
1159 tmpBiEdgeFaces[0] = newFaceIndex[1];
1160 tmpBiEdgeFaces[1] = newFaceIndex[3];
1162 newEdgeIndex[3] =
1163 (
1164 insertEdge
1165 (
1166 whichEdgePatch(commonEdgeIndex[0]),
1167 edge(newPointIndex[0], otherEdgePoint[0]),
1168 tmpBiEdgeFaces
1169 )
1170 );
1172 // Add this edge to the map.
1173 map.addEdge(newEdgeIndex[3]);
1175 tmpBiEdgeFaces[0] = c0BdyIndex[1];
1176 tmpBiEdgeFaces[1] = newFaceIndex[3];
1178 newEdgeIndex[4] =
1179 (
1180 insertEdge
1181 (
1182 whichEdgePatch(commonEdgeIndex[1]),
1183 edge(newPointIndex[1], nextToOtherPoint[1]),
1184 tmpBiEdgeFaces
1185 )
1186 );
1188 // Add this edge to the map.
1189 map.addEdge(newEdgeIndex[4]);
1191 // Now that edges are defined, configure faceEdges
1192 // for all new faces
1194 // The quad interior face; Owner = cell[0] & Neighbour = newCell[0]
1195 tmpQFEdges[0] = newEdgeIndex[0];
1196 tmpQFEdges[1] = newEdgeIndex[1];
1197 tmpQFEdges[2] = otherEdgeIndex[2];
1198 tmpQFEdges[3] = newEdgeIndex[2];
1199 faceEdges_[newFaceIndex[0]] = tmpQFEdges;
1201 // Second boundary face; Owner = newCell[0] & Neighbour = [-1]
1202 tmpTFEdges[0] = newEdgeIndex[3];
1203 tmpTFEdges[1] = cornerEdgeIndex[0];
1204 tmpTFEdges[2] = newEdgeIndex[1];
1205 faceEdges_[newFaceIndex[1]] = tmpTFEdges;
1207 // Third boundary face; Owner = c[0] & Neighbour = [-1]
1208 tmpTFEdges[0] = newEdgeIndex[2];
1209 tmpTFEdges[1] = cornerEdgeIndex[1];
1210 tmpTFEdges[2] = commonEdgeIndex[1];
1211 faceEdges_[newFaceIndex[2]] = tmpTFEdges;
1213 // The quad face from bisection:
1214 tmpQFEdges[0] = otherEdgeIndex[1];
1215 tmpQFEdges[1] = newEdgeIndex[3];
1216 tmpQFEdges[2] = newEdgeIndex[0];
1217 tmpQFEdges[3] = newEdgeIndex[4];
1218 faceEdges_[newFaceIndex[3]] = tmpQFEdges;
1220 meshOps::replaceLabel
1221 (
1222 commonEdgeIndex[1],
1223 newEdgeIndex[4],
1224 faceEdges_[c0BdyIndex[1]]
1225 );
1227 meshOps::replaceLabel
1228 (
1229 c0BdyIndex[1],
1230 newFaceIndex[2],
1231 edgeFaces_[commonEdgeIndex[1]]
1232 );
1234 if (debug > 2)
1235 {
1236 Pout<< "Modified Cell[0]: "
1237 << c0 << ": " << oldCells[0] << endl;
1239 forAll(oldCells[0], faceI)
1240 {
1241 const labelList& fE = faceEdges_[oldCells[0][faceI]];
1243 Pout<< oldCells[0][faceI]
1244 << ": " << faces_[oldCells[0][faceI]]
1245 << " fE: " << fE
1246 << endl;
1248 forAll(fE, edgeI)
1249 {
1250 const labelList& eF = edgeFaces_[fE[edgeI]];
1252 Pout<< '\t' << fE[edgeI]
1253 << ": " << edges_[fE[edgeI]]
1254 << " eF: " << eF
1255 << endl;
1256 }
1257 }
1259 Pout<< "New Cell[0]: " << newCellIndex[0]
1260 << ": " << newCells[0] << endl;
1262 forAll(newCells[0], faceI)
1263 {
1264 const labelList& fE = faceEdges_[newCells[0][faceI]];
1266 Pout<< newCells[0][faceI]
1267 << ": " << faces_[newCells[0][faceI]]
1268 << " fE: " << fE
1269 << endl;
1271 forAll(fE, edgeI)
1272 {
1273 const labelList& eF = edgeFaces_[fE[edgeI]];
1275 Pout<< '\t' << fE[edgeI]
1276 << ": " << edges_[fE[edgeI]]
1277 << " eF: " << eF
1278 << endl;
1279 }
1280 }
1281 }
1282 }
1283 else
1284 {
1285 oldCells[1] = cells_[c1];
1287 // Add a new prism cell to the end of the list.
1288 // Currently invalid, but will be updated later.
1289 newCellIndex[1] = insertCell(newCells[1], lengthScale_[c1]);
1291 // Add this cell to the map.
1292 map.addCell(newCellIndex[1]);
1294 if (debug > 2)
1295 {
1296 Pout<< "Cell[1]: " << c1 << ": " << oldCells[1] << endl;
1298 forAll(oldCells[1], faceI)
1299 {
1300 const labelList& fE = faceEdges_[oldCells[1][faceI]];
1302 Pout<< oldCells[1][faceI] << ": "
1303 << faces_[oldCells[1][faceI]]
1304 << " fE: " << fE
1305 << endl;
1307 forAll(fE, edgeI)
1308 {
1309 const labelList& eF = edgeFaces_[fE[edgeI]];
1311 Pout<< '\t' << fE[edgeI]
1312 << ": " << edges_[fE[edgeI]]
1313 << " eF: " << eF
1314 << endl;
1315 }
1316 }
1317 }
1319 // Find the interior face that contains otherEdgeIndex[1]
1320 found = false;
1322 const labelList& e2 = faceEdges_[c1IntIndex[0]];
1324 forAll(e2, edgeI)
1325 {
1326 if (e2[edgeI] == otherEdgeIndex[1])
1327 {
1328 meshOps::replaceLabel
1329 (
1330 c1IntIndex[0],
1331 -1,
1332 oldCells[1]
1333 );
1335 replaceFace = c1IntIndex[0];
1336 retainFace = c1IntIndex[1];
1337 found = true;
1338 break;
1339 }
1340 }
1342 if (!found)
1343 {
1344 // The edge was not found before
1345 meshOps::replaceLabel
1346 (
1347 c1IntIndex[1],
1348 -1,
1349 oldCells[1]
1350 );
1352 replaceFace = c1IntIndex[1];
1353 retainFace = c1IntIndex[0];
1354 }
1356 // Check if face reversal is necessary for the replacement
1357 if (owner_[replaceFace] == c1)
1358 {
1359 if (neighbour_[replaceFace] == -1)
1360 {
1361 // Change the owner
1362 owner_[replaceFace] = newCellIndex[1];
1363 }
1364 else
1365 {
1366 // This face has to be reversed
1367 faces_[replaceFace] = faces_[replaceFace].reverseFace();
1368 owner_[replaceFace] = neighbour_[replaceFace];
1369 neighbour_[replaceFace] = newCellIndex[1];
1371 setFlip(replaceFace);
1372 }
1373 }
1374 else
1375 {
1376 // Keep owner, but change neighbour
1377 neighbour_[replaceFace] = newCellIndex[1];
1378 }
1380 // Define attributes for the new prism cell
1381 newCells[1][0] = replaceFace;
1383 // The quad interior face resulting from bisection;
1384 // Owner = newCell[0] & Neighbour = newCell[1]
1385 tmpQuadFace[0] = newPointIndex[1];
1386 tmpQuadFace[1] = nextToOtherPoint[1];
1387 tmpQuadFace[2] = otherEdgePoint[0];
1388 tmpQuadFace[3] = newPointIndex[0];
1390 newFaceIndex[3] =
1391 (
1392 insertFace
1393 (
1394 -1,
1395 tmpQuadFace,
1396 newCellIndex[0],
1397 newCellIndex[1]
1398 )
1399 );
1401 // Add a faceEdges entry as well
1402 faceEdges_.append(tmpQFEdges);
1404 // Add this face to the map.
1405 map.addFace(newFaceIndex[3]);
1407 // Correct edgeFaces for otherEdgeIndex[1]
1408 meshOps::replaceLabel
1409 (
1410 fIndex,
1411 newFaceIndex[3],
1412 edgeFaces_[otherEdgeIndex[1]]
1413 );
1415 meshOps::replaceLabel
1416 (
1417 -1,
1418 newFaceIndex[3],
1419 newCells[0]
1420 );
1422 meshOps::replaceLabel
1423 (
1424 -1,
1425 newFaceIndex[3],
1426 newCells[1]
1427 );
1429 newCells[1][1] = newFaceIndex[3];
1431 // Check for common edges among the two boundary faces
1432 // Find the isolated point on both boundary faces of cell[1]
1433 if
1434 (
1435 meshOps::findCommonEdge
1436 (
1437 c1BdyIndex[0],
1438 c0BdyIndex[0],
1439 faceEdges_,
1440 commonEdgeIndex[2]
1441 )
1442 )
1443 {
1444 meshOps::findCommonEdge
1445 (
1446 c1BdyIndex[1],
1447 c0BdyIndex[1],
1448 faceEdges_,
1449 commonEdgeIndex[3]
1450 );
1452 commonFaceIndex[2] = c1BdyIndex[0];
1453 commonFaceIndex[3] = c1BdyIndex[1];
1454 }
1455 else
1456 {
1457 meshOps::findCommonEdge
1458 (
1459 c1BdyIndex[0],
1460 c0BdyIndex[1],
1461 faceEdges_,
1462 commonEdgeIndex[3]
1463 );
1465 meshOps::findCommonEdge
1466 (
1467 c1BdyIndex[1],
1468 c0BdyIndex[0],
1469 faceEdges_,
1470 commonEdgeIndex[2]
1471 );
1473 commonFaceIndex[2] = c1BdyIndex[1];
1474 commonFaceIndex[3] = c1BdyIndex[0];
1475 }
1477 commonEdges[2] = edges_[commonEdgeIndex[2]];
1478 commonEdges[3] = edges_[commonEdgeIndex[3]];
1480 if (debug > 2)
1481 {
1482 Pout<< "Common edges: " << nl
1483 << commonEdgeIndex[2] << ": " << commonEdges[2] << nl
1484 << commonEdgeIndex[3] << ": " << commonEdges[3]
1485 << endl;
1486 }
1488 meshOps::findIsolatedPoint
1489 (
1490 faces_[commonFaceIndex[2]],
1491 commonEdges[2],
1492 otherPointIndex[2],
1493 nextToOtherPoint[2]
1494 );
1496 meshOps::findIsolatedPoint
1497 (
1498 faces_[commonFaceIndex[3]],
1499 commonEdges[3],
1500 otherPointIndex[3],
1501 nextToOtherPoint[3]
1502 );
1504 // For convenience...
1505 otherEdgePoint[2] = commonEdges[2].otherVertex(nextToOtherPoint[2]);
1506 otherEdgePoint[3] = commonEdges[3].otherVertex(nextToOtherPoint[3]);
1508 // Modify the two existing triangle boundary faces
1510 // Zeroth boundary face - Owner = newCell[1], Neighbour = -1
1511 meshOps::replaceLabel
1512 (
1513 otherEdgePoint[2],
1514 newPointIndex[0],
1515 faces_[commonFaceIndex[2]]
1516 );
1518 owner_[commonFaceIndex[2]] = newCellIndex[1];
1520 meshOps::replaceLabel
1521 (
1522 commonFaceIndex[2],
1523 -1,
1524 oldCells[1]
1525 );
1527 newCells[1][2] = commonFaceIndex[2];
1529 // First boundary face - Owner = c[1] & Neighbour [-1] (unchanged)
1530 meshOps::replaceLabel
1531 (
1532 otherEdgePoint[3],
1533 newPointIndex[1],
1534 faces_[commonFaceIndex[3]]
1535 );
1537 // New interior face; Owner = cell[1] & Neighbour = newCell[1]
1538 tmpQuadFace[0] = newPointIndex[0];
1539 tmpQuadFace[1] = otherPointIndex[2];
1540 tmpQuadFace[2] = otherPointIndex[3];
1541 tmpQuadFace[3] = newPointIndex[1];
1543 newFaceIndex[4] =
1544 (
1545 insertFace
1546 (
1547 -1,
1548 tmpQuadFace,
1549 c1,
1550 newCellIndex[1]
1551 )
1552 );
1554 // Add a faceEdges entry as well
1555 faceEdges_.append(tmpQFEdges);
1557 // Add this face to the map.
1558 map.addFace(newFaceIndex[4]);
1560 // Find the common edge between quad/quad faces...
1561 meshOps::findCommonEdge
1562 (
1563 c1IntIndex[0],
1564 c1IntIndex[1],
1565 faceEdges_,
1566 otherEdgeIndex[3]
1567 );
1569 // ... and size-up edgeFaces for the edge
1570 meshOps::sizeUpList
1571 (
1572 newFaceIndex[4],
1573 edgeFaces_[otherEdgeIndex[3]]
1574 );
1576 meshOps::replaceLabel
1577 (
1578 -1,
1579 newFaceIndex[4],
1580 newCells[1]
1581 );
1583 meshOps::replaceLabel
1584 (
1585 -1,
1586 newFaceIndex[4],
1587 oldCells[1]
1588 );
1590 // Second boundary face; Owner = cell[1] & Neighbour [-1]
1591 tmpTriFace[0] = otherPointIndex[2];
1592 tmpTriFace[1] = newPointIndex[0];
1593 tmpTriFace[2] = otherEdgePoint[2];
1595 newFaceIndex[5] =
1596 (
1597 insertFace
1598 (
1599 whichPatch(commonFaceIndex[2]),
1600 tmpTriFace,
1601 c1,
1602 -1
1603 )
1604 );
1606 // Add a faceEdges entry as well
1607 faceEdges_.append(tmpTFEdges);
1609 // Add this face to the map.
1610 map.addFace(newFaceIndex[5]);
1612 meshOps::replaceLabel
1613 (
1614 -1,
1615 newFaceIndex[5],
1616 oldCells[1]
1617 );
1619 // Third boundary face; Owner = newCell[1] & Neighbour [-1]
1620 tmpTriFace[0] = otherPointIndex[3];
1621 tmpTriFace[1] = newPointIndex[1];
1622 tmpTriFace[2] = otherEdgePoint[3];
1624 newFaceIndex[6] =
1625 (
1626 insertFace
1627 (
1628 whichPatch(commonFaceIndex[3]),
1629 tmpTriFace,
1630 newCellIndex[1],
1631 -1
1632 )
1633 );
1635 // Add a faceEdges entry as well
1636 faceEdges_.append(tmpTFEdges);
1638 // Add this face to the map.
1639 map.addFace(newFaceIndex[6]);
1641 meshOps::replaceLabel
1642 (
1643 -1,
1644 newFaceIndex[6],
1645 newCells[1]
1646 );
1648 // Create / modify edges...
1649 labelList tmpQuadEdgeFaces(4, -1);
1651 // The internal edge bisecting the face
1652 tmpQuadEdgeFaces[0] = fIndex;
1653 tmpQuadEdgeFaces[1] = newFaceIndex[0];
1654 tmpQuadEdgeFaces[2] = newFaceIndex[3];
1655 tmpQuadEdgeFaces[3] = newFaceIndex[4];
1657 newEdgeIndex[0] =
1658 (
1659 insertEdge
1660 (
1661 -1,
1662 edge(newPointIndex[0], newPointIndex[1]),
1663 tmpQuadEdgeFaces
1664 )
1665 );
1667 // Add this edge to the map.
1668 map.addEdge(newEdgeIndex[0]);
1670 // Replace an edge on the bisected face
1671 meshOps::replaceLabel
1672 (
1673 otherEdgeIndex[1],
1674 newEdgeIndex[0],
1675 faceEdges_[fIndex]
1676 );
1678 // Create / replace side edges created from face bisection
1679 tmpTriEdgeFaces[0] = commonFaceIndex[2];
1680 tmpTriEdgeFaces[1] = newFaceIndex[3];
1681 tmpTriEdgeFaces[2] = newFaceIndex[1];
1683 newEdgeIndex[3] =
1684 (
1685 insertEdge
1686 (
1687 whichEdgePatch(commonEdgeIndex[2]),
1688 edge(newPointIndex[0], nextToOtherPoint[2]),
1689 tmpTriEdgeFaces
1690 )
1691 );
1693 // Add this edge to the map.
1694 map.addEdge(newEdgeIndex[3]);
1696 tmpTriEdgeFaces[0] = c0BdyIndex[1];
1697 tmpTriEdgeFaces[1] = newFaceIndex[3];
1698 tmpTriEdgeFaces[2] = newFaceIndex[6];
1700 newEdgeIndex[4] =
1701 (
1702 insertEdge
1703 (
1704 whichEdgePatch(commonEdgeIndex[3]),
1705 edge(newPointIndex[1], otherEdgePoint[3]),
1706 tmpTriEdgeFaces
1707 )
1708 );
1710 // Add this edge to the map.
1711 map.addEdge(newEdgeIndex[4]);
1713 // The edge bisecting the second boundary triangular face
1714 tmpTriEdgeFaces[0] = commonFaceIndex[2];
1715 tmpTriEdgeFaces[1] = newFaceIndex[4];
1716 tmpTriEdgeFaces[2] = newFaceIndex[5];
1718 newEdgeIndex[5] =
1719 (
1720 insertEdge
1721 (
1722 whichEdgePatch(commonEdgeIndex[2]),
1723 edge(newPointIndex[0], otherPointIndex[2]),
1724 tmpTriEdgeFaces
1725 )
1726 );
1728 // Add this edge to the map.
1729 map.addEdge(newEdgeIndex[5]);
1731 // Find the common edge between the quad/tri faces...
1732 meshOps::findCommonEdge
1733 (
1734 commonFaceIndex[2],
1735 retainFace,
1736 faceEdges_,
1737 cornerEdgeIndex[2]
1738 );
1740 // ...and correct faceEdges / edgeFaces
1741 meshOps::replaceLabel
1742 (
1743 cornerEdgeIndex[2],
1744 newEdgeIndex[5],
1745 faceEdges_[commonFaceIndex[2]]
1746 );
1748 meshOps::replaceLabel
1749 (
1750 commonFaceIndex[2],
1751 newFaceIndex[5],
1752 edgeFaces_[cornerEdgeIndex[2]]
1753 );
1755 // The edge bisecting the third boundary triangular face
1756 tmpTriEdgeFaces[0] = commonFaceIndex[3];
1757 tmpTriEdgeFaces[1] = newFaceIndex[4];
1758 tmpTriEdgeFaces[2] = newFaceIndex[6];
1760 newEdgeIndex[6] =
1761 (
1762 insertEdge
1763 (
1764 whichEdgePatch(commonEdgeIndex[3]),
1765 edge(newPointIndex[1], otherPointIndex[3]),
1766 tmpTriEdgeFaces
1767 )
1768 );
1770 // Add this edge to the map.
1771 map.addEdge(newEdgeIndex[6]);
1773 // Find the common edge between the quad/tri faces...
1774 meshOps::findCommonEdge
1775 (
1776 commonFaceIndex[3],
1777 replaceFace,
1778 faceEdges_,
1779 cornerEdgeIndex[3]
1780 );
1782 // ...and correct faceEdges / edgeFaces
1783 meshOps::replaceLabel
1784 (
1785 cornerEdgeIndex[3],
1786 newEdgeIndex[6],
1787 faceEdges_[commonFaceIndex[3]]
1788 );
1790 meshOps::replaceLabel
1791 (
1792 commonFaceIndex[3],
1793 newFaceIndex[6],
1794 edgeFaces_[cornerEdgeIndex[3]]
1795 );
1797 // Now that edges are defined, configure faceEdges
1798 // for all new faces
1800 // The quad interior face; Owner = c[0] & Neighbour = newCell[0]
1801 tmpQFEdges[0] = newEdgeIndex[0];
1802 tmpQFEdges[1] = newEdgeIndex[1];
1803 tmpQFEdges[2] = otherEdgeIndex[2];
1804 tmpQFEdges[3] = newEdgeIndex[2];
1805 faceEdges_[newFaceIndex[0]] = tmpQFEdges;
1807 // Second boundary face; Owner = newCell[0] & Neighbour = [-1]
1808 tmpTFEdges[0] = newEdgeIndex[3];
1809 tmpTFEdges[1] = cornerEdgeIndex[0];
1810 tmpTFEdges[2] = newEdgeIndex[1];
1811 faceEdges_[newFaceIndex[1]] = tmpTFEdges;
1813 // Third boundary face; Owner = c[0] & Neighbour = [-1]
1814 tmpTFEdges[0] = newEdgeIndex[2];
1815 tmpTFEdges[1] = cornerEdgeIndex[1];
1816 tmpTFEdges[2] = commonEdgeIndex[3];
1817 faceEdges_[newFaceIndex[2]] = tmpTFEdges;
1819 // The quad face from bisection:
1820 tmpQFEdges[0] = otherEdgeIndex[1];
1821 tmpQFEdges[1] = newEdgeIndex[3];
1822 tmpQFEdges[2] = newEdgeIndex[0];
1823 tmpQFEdges[3] = newEdgeIndex[4];
1824 faceEdges_[newFaceIndex[3]] = tmpQFEdges;
1826 // The quad interior face; Owner = c[1] & Neighbour = newCell[1]
1827 tmpQFEdges[0] = newEdgeIndex[0];
1828 tmpQFEdges[1] = newEdgeIndex[5];
1829 tmpQFEdges[2] = otherEdgeIndex[3];
1830 tmpQFEdges[3] = newEdgeIndex[6];
1831 faceEdges_[newFaceIndex[4]] = tmpQFEdges;
1833 // Second boundary face; Owner = c[1] & Neighbour = [-1]
1834 tmpTFEdges[0] = commonEdgeIndex[2];
1835 tmpTFEdges[1] = cornerEdgeIndex[2];
1836 tmpTFEdges[2] = newEdgeIndex[5];
1837 faceEdges_[newFaceIndex[5]] = tmpTFEdges;
1839 // Third boundary face; Owner = newCell[1] & Neighbour = [-1]
1840 tmpTFEdges[0] = newEdgeIndex[4];
1841 tmpTFEdges[1] = cornerEdgeIndex[3];
1842 tmpTFEdges[2] = newEdgeIndex[6];
1843 faceEdges_[newFaceIndex[6]] = tmpTFEdges;
1845 meshOps::replaceLabel
1846 (
1847 commonEdgeIndex[1],
1848 newEdgeIndex[4],
1849 faceEdges_[c0BdyIndex[1]]
1850 );
1852 meshOps::replaceLabel
1853 (
1854 c0BdyIndex[1],
1855 newFaceIndex[2],
1856 edgeFaces_[commonEdgeIndex[1]]
1857 );
1859 meshOps::replaceLabel
1860 (
1861 commonEdgeIndex[2],
1862 newEdgeIndex[3],
1863 faceEdges_[commonFaceIndex[2]]
1864 );
1866 meshOps::replaceLabel
1867 (
1868 commonFaceIndex[2],
1869 newFaceIndex[5],
1870 edgeFaces_[commonEdgeIndex[2]]
1871 );
1873 if (debug > 2)
1874 {
1875 Pout<< nl << "Modified Cell[0]: "
1876 << c0 << ": " << oldCells[0] << endl;
1878 forAll(oldCells[0], faceI)
1879 {
1880 const labelList& fE = faceEdges_[oldCells[0][faceI]];
1882 Pout<< oldCells[0][faceI]
1883 << ": " << faces_[oldCells[0][faceI]]
1884 << " fE: " << fE
1885 << endl;
1887 forAll(fE, edgeI)
1888 {
1889 const labelList& eF = edgeFaces_[fE[edgeI]];
1891 Pout<< '\t' << fE[edgeI]
1892 << ": " << edges_[fE[edgeI]]
1893 << " eF: " << eF
1894 << endl;
1895 }
1896 }
1898 Pout<< "New Cell[0]: "
1899 << newCellIndex[0] << ": " << newCells[0] << endl;
1901 forAll(newCells[0], faceI)
1902 {
1903 const labelList& fE = faceEdges_[newCells[0][faceI]];
1905 Pout<< newCells[0][faceI] << ": "
1906 << faces_[newCells[0][faceI]]
1907 << " fE: " << fE
1908 << endl;
1910 forAll(fE, edgeI)
1911 {
1912 const labelList& eF = edgeFaces_[fE[edgeI]];
1914 Pout<< '\t' << fE[edgeI]
1915 << ": " << edges_[fE[edgeI]]
1916 << " eF: " << eF
1917 << endl;
1918 }
1919 }
1921 Pout<< nl << "Modified Cell[1]: "
1922 << c1 << ": " << oldCells[1] << endl;
1924 forAll(oldCells[1], faceI)
1925 {
1926 const labelList& fE = faceEdges_[oldCells[1][faceI]];
1928 Pout<< oldCells[1][faceI] << ": "
1929 << faces_[oldCells[1][faceI]]
1930 << " fE: " << fE
1931 << endl;
1933 forAll(fE, edgeI)
1934 {
1935 const labelList& eF = edgeFaces_[fE[edgeI]];
1937 Pout<< '\t' << fE[edgeI]
1938 << ": " << edges_[fE[edgeI]]
1939 << " eF: " << eF
1940 << endl;
1941 }
1942 }
1944 Pout<< "New Cell[1]: "
1945 << newCellIndex[1] << ": " << newCells[1] << endl;
1947 forAll(newCells[1], faceI)
1948 {
1949 const labelList& fE = faceEdges_[newCells[1][faceI]];
1951 Pout<< newCells[1][faceI] << ": "
1952 << faces_[newCells[1][faceI]]
1953 << " fE: " << fE
1954 << endl;
1956 forAll(fE, edgeI)
1957 {
1958 const labelList& eF = edgeFaces_[fE[edgeI]];
1960 Pout<< '\t' << fE[edgeI]
1961 << ": " << edges_[fE[edgeI]]
1962 << " eF: " << eF
1963 << endl;
1964 }
1965 }
1966 }
1968 // Update the cell list.
1969 cells_[c1] = oldCells[1];
1970 cells_[newCellIndex[1]] = newCells[1];
1971 }
1973 // Update the cell list.
1974 cells_[c0] = oldCells[0];
1975 cells_[newCellIndex[0]] = newCells[0];
1977 // Modify point labels for common edges
1978 if (edges_[commonEdgeIndex[0]].start() == otherEdgePoint[0])
1979 {
1980 edges_[commonEdgeIndex[0]].start() = newPointIndex[0];
1981 }
1982 else
1983 {
1984 edges_[commonEdgeIndex[0]].end() = newPointIndex[0];
1985 }
1987 if (edges_[commonEdgeIndex[1]].start() == nextToOtherPoint[1])
1988 {
1989 edges_[commonEdgeIndex[1]].start() = newPointIndex[1];
1990 }
1991 else
1992 {
1993 edges_[commonEdgeIndex[1]].end() = newPointIndex[1];
1994 }
1996 if (coupledModification_)
1997 {
1998 // Alias for convenience...
1999 const changeMap& slaveMap = slaveMaps[0];
2001 const label pI = slaveMap.patchIndex();
2003 // Fetch the appropriate coupleMap
2004 const coupleMap* cMapPtr = NULL;
2006 if (localCouple && !procCouple)
2007 {
2008 cMapPtr = &(patchCoupling_[pI].map());
2009 }
2010 else
2011 if (procCouple && !localCouple)
2012 {
2013 cMapPtr = &(recvMeshes_[pI].map());
2014 }
2016 // Alias for convenience
2017 const coupleMap& cMap = *cMapPtr;
2019 // Add the new points to the coupling map
2020 const List<objectMap>& apList = slaveMap.addedPointList();
2022 if (bisectingSlave)
2023 {
2024 // Update reverse pointMap
2025 cMap.mapMaster
2026 (
2027 coupleMap::POINT,
2028 newPointIndex[0],
2029 apList[0].index()
2030 );
2032 cMap.mapMaster
2033 (
2034 coupleMap::POINT,
2035 newPointIndex[1],
2036 apList[1].index()
2037 );
2039 // Update pointMap
2040 cMap.mapSlave
2041 (
2042 coupleMap::POINT,
2043 apList[0].index(),
2044 newPointIndex[0]
2045 );
2047 cMap.mapSlave
2048 (
2049 coupleMap::POINT,
2050 apList[1].index(),
2051 newPointIndex[1]
2052 );
2053 }
2054 else
2055 {
2056 // Update pointMap
2057 cMap.mapSlave
2058 (
2059 coupleMap::POINT,
2060 newPointIndex[0],
2061 apList[0].index()
2062 );
2064 cMap.mapSlave
2065 (
2066 coupleMap::POINT,
2067 newPointIndex[1],
2068 apList[1].index()
2069 );
2071 // Update reverse pointMap
2072 cMap.mapMaster
2073 (
2074 coupleMap::POINT,
2075 apList[0].index(),
2076 newPointIndex[0]
2077 );
2079 cMap.mapMaster
2080 (
2081 coupleMap::POINT,
2082 apList[1].index(),
2083 newPointIndex[1]
2084 );
2085 }
2087 // Create a master/slave entry for the new face on the patch.
2088 FixedList<bool, 2> foundMatch(false);
2089 FixedList<label, 2> checkFaces(-1);
2091 // Fill in the faces to check for...
2092 checkFaces[0] = fIndex;
2093 checkFaces[1] = newFaceIndex[3];
2095 const List<objectMap>& afList = slaveMap.addedFaceList();
2097 forAll (checkFaces, indexI)
2098 {
2099 const face& mFace = faces_[checkFaces[indexI]];
2101 label sFaceIndex = -1;
2103 forAll(afList, sfI)
2104 {
2105 const face* facePtr = NULL;
2107 if (localCouple && !procCouple)
2108 {
2109 facePtr = &(faces_[afList[sfI].index()]);
2110 }
2111 else
2112 if (procCouple && !localCouple)
2113 {
2114 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2116 facePtr = &(mesh.faces_[afList[sfI].index()]);
2117 }
2119 const face& tFace = *facePtr;
2120 FixedList<bool, 4> cP(false);
2122 forAll(mFace, pointI)
2123 {
2124 const Map<label>& pointMap =
2125 (
2126 cMap.entityMap(coupleMap::POINT)
2127 );
2129 if (tFace.which(pointMap[mFace[pointI]]) > -1)
2130 {
2131 cP[pointI] = true;
2132 }
2133 }
2135 if (cP[0] && cP[1] && cP[2] && cP[3])
2136 {
2137 sFaceIndex = afList[sfI].index();
2138 foundMatch[indexI] = true;
2139 break;
2140 }
2141 }
2143 if (foundMatch[indexI])
2144 {
2145 if (bisectingSlave)
2146 {
2147 cMap.mapMaster
2148 (
2149 coupleMap::FACE,
2150 checkFaces[indexI],
2151 sFaceIndex
2152 );
2154 cMap.mapSlave
2155 (
2156 coupleMap::FACE,
2157 sFaceIndex,
2158 checkFaces[indexI]
2159 );
2160 }
2161 else
2162 {
2163 cMap.mapSlave
2164 (
2165 coupleMap::FACE,
2166 checkFaces[indexI],
2167 sFaceIndex
2168 );
2170 cMap.mapMaster
2171 (
2172 coupleMap::FACE,
2173 sFaceIndex,
2174 checkFaces[indexI]
2175 );
2176 }
2177 }
2178 else
2179 {
2180 FatalErrorIn
2181 (
2182 "\n"
2183 "const changeMap "
2184 "dynamicTopoFvMesh::bisectQuadFace\n"
2185 "(\n"
2186 " const label fIndex,\n"
2187 " const changeMap& masterMap,\n"
2188 " bool checkOnly,\n"
2189 " bool forceOp\n"
2190 ")\n"
2191 )
2192 << "Failed to build coupled maps." << nl
2193 << " masterFace: " << mFace << nl
2194 << " Added slave faces: " << afList << nl
2195 << abort(FatalError);
2196 }
2197 }
2199 // Push operation into coupleMap
2200 cMap.pushOperation
2201 (
2202 slaveMap.index(),
2203 coupleMap::BISECTION
2204 );
2205 }
2207 // Write out VTK files after change
2208 if (debug > 3)
2209 {
2210 labelList cellHull(4, -1);
2212 cellHull[0] = owner_[fIndex];
2213 cellHull[1] = neighbour_[fIndex];
2214 cellHull[2] = owner_[newFaceIndex[3]];
2215 cellHull[3] = neighbour_[newFaceIndex[3]];
2217 writeVTK
2218 (
2219 Foam::name(fIndex)
2220 + "_Bisect_1",
2221 cellHull
2222 );
2223 }
2225 // Fill-in mapping information
2226 FixedList<label, 4> mapCells(-1);
2228 mapCells[0] = c0;
2229 mapCells[1] = newCellIndex[0];
2231 if (c1 != -1)
2232 {
2233 mapCells[2] = c1;
2234 mapCells[3] = newCellIndex[1];
2235 }
2237 labelList mC(1, c0);
2239 forAll(mapCells, cellI)
2240 {
2241 if (mapCells[cellI] == -1)
2242 {
2243 continue;
2244 }
2246 // Set the mapping for this cell
2247 setCellMapping(mapCells[cellI], mC);
2248 }
2250 // Set fill-in mapping information for the modified face.
2251 if (c1 == -1)
2252 {
2253 // Set the mapping for this face
2254 setFaceMapping(fIndex, labelList(1, fIndex));
2255 }
2256 else
2257 {
2258 // Internal face. Default mapping.
2259 setFaceMapping(fIndex);
2260 }
2262 forAll(newFaceIndex, faceI)
2263 {
2264 if (newFaceIndex[faceI] == -1)
2265 {
2266 continue;
2267 }
2269 // Check for boundary faces
2270 if (neighbour_[newFaceIndex[faceI]] == -1)
2271 {
2272 // Boundary face. Compute mapping.
2273 labelList mC;
2275 if (faces_[newFaceIndex[faceI]].size() == 4)
2276 {
2277 // Quad-face on boundary
2278 mC.setSize(1, fIndex);
2279 }
2280 else
2281 if (faces_[newFaceIndex[faceI]].size() == 3)
2282 {
2283 label triFacePatch = whichPatch(newFaceIndex[faceI]);
2285 // Fetch face-normals
2286 vector tfNorm, f0Norm, f1Norm;
2288 tfNorm = faces_[newFaceIndex[faceI]].normal(oldPoints_);
2289 f0Norm = faces_[c0BdyIndex[0]].normal(oldPoints_);
2290 f1Norm = faces_[c0BdyIndex[1]].normal(oldPoints_);
2292 // Tri-face on boundary. Perform normal checks
2293 // also, because of empty patches.
2294 if
2295 (
2296 (whichPatch(c0BdyIndex[0]) == triFacePatch) &&
2297 ((tfNorm & f0Norm) > 0.0)
2298 )
2299 {
2300 mC.setSize(1, c0BdyIndex[0]);
2301 }
2302 else
2303 if
2304 (
2305 (whichPatch(c0BdyIndex[1]) == triFacePatch) &&
2306 ((tfNorm & f1Norm) > 0.0)
2307 )
2308 {
2309 mC.setSize(1, c0BdyIndex[1]);
2310 }
2311 else
2312 {
2313 FatalErrorIn
2314 (
2315 "\n"
2316 "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
2317 "(\n"
2318 " const label fIndex,\n"
2319 " const changeMap& masterMap,\n"
2320 " bool checkOnly\n"
2321 ")\n"
2322 )
2323 << " Unable to find patch for face: "
2324 << newFaceIndex[faceI] << ":: "
2325 << faces_[newFaceIndex[faceI]] << nl
2326 << " Patch: " << triFacePatch << nl
2327 << abort(FatalError);
2328 }
2329 }
2331 // Set the mapping for this face
2332 setFaceMapping(newFaceIndex[faceI], mC);
2333 }
2334 else
2335 {
2336 // Internal quad-faces get default mapping.
2337 setFaceMapping(newFaceIndex[faceI]);
2338 }
2339 }
2341 // Set the flag
2342 topoChangeFlag_ = true;
2344 // Increment the counter
2345 statistics_[3]++;
2347 // Increment surface-counter
2348 if (c1 == -1)
2349 {
2350 // Do not update stats for processor patches
2351 if (!processorCoupledEntity(fIndex))
2352 {
2353 statistics_[5]++;
2354 }
2355 }
2357 // Increment the number of modifications
2358 statistics_[0]++;
2360 // Specify that the operation was successful
2361 map.type() = 1;
2363 // Return the changeMap
2364 return map;
2365 }
2368 // Method for the bisection of an edge in 3D
2369 // - Returns a changeMap with a type specifying:
2370 // 1: Bisection was successful
2371 // -1: Bisection failed since max number of topo-changes was reached.
2372 // -2: Bisection failed since resulting quality would be unacceptable.
2373 // -3: Bisection failed since edge was on a noRefinement patch.
2374 // - AddedPoints contain the index of the newly added point.
2375 const changeMap dynamicTopoFvMesh::bisectEdge
2376 (
2377 const label eIndex,
2378 bool checkOnly,
2379 bool forceOp
2380 )
2381 {
2382 // Edge bisection performs the following operations:
2383 // [1] Add a point at middle of the edge
2384 // [2] Bisect all faces surrounding this edge
2385 // [3] Bisect all cells surrounding this edge
2386 // [4] Create internal/external edges for each bisected face
2387 // [5] Create internal faces for each bisected cell
2388 // Update faceEdges and edgeFaces information
2390 // For 2D meshes, perform face-bisection
2391 if (twoDMesh_)
2392 {
2393 return bisectQuadFace(eIndex, changeMap(), checkOnly);
2394 }
2396 // Figure out which thread this is...
2397 label tIndex = self();
2399 // Prepare the changeMaps
2400 changeMap map;
2401 List<changeMap> slaveMaps;
2402 bool bisectingSlave = false;
2404 if
2405 (
2406 (statistics_[0] > maxModifications_) &&
2407 (maxModifications_ > -1)
2408 )
2409 {
2410 // Reached the max allowable topo-changes.
2411 stack(tIndex).clear();
2413 return map;
2414 }
2416 // Check if edgeRefinements are to be avoided on patch.
2417 const labelList& eF = edgeFaces_[eIndex];
2419 forAll(eF, fI)
2420 {
2421 label fPatch = whichPatch(eF[fI]);
2423 if (baseMesh_.lengthEstimator().checkRefinementPatch(fPatch))
2424 {
2425 map.type() = -3;
2427 return map;
2428 }
2429 }
2431 // Sanity check: Is the index legitimate?
2432 if (eIndex < 0)
2433 {
2434 FatalErrorIn
2435 (
2436 "\n"
2437 "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2438 "(\n"
2439 " const label eIndex,\n"
2440 " bool checkOnly,\n"
2441 " bool forceOp\n"
2442 ")\n"
2443 )
2444 << " Invalid index: " << eIndex
2445 << abort(FatalError);
2446 }
2448 // If coupled modification is set, and this is a
2449 // master edge, bisect its slaves as well.
2450 bool localCouple = false, procCouple = false;
2452 if (coupledModification_)
2453 {
2454 const edge& eCheck = edges_[eIndex];
2456 const label edgeEnum = coupleMap::EDGE;
2458 // Is this a locally coupled edge (either master or slave)?
2459 if (locallyCoupledEntity(eIndex, true))
2460 {
2461 localCouple = true;
2462 procCouple = false;
2463 }
2464 else
2465 if (processorCoupledEntity(eIndex))
2466 {
2467 procCouple = true;
2468 localCouple = false;
2469 }
2471 if (localCouple && !procCouple)
2472 {
2473 // Determine the slave index.
2474 label sIndex = -1, pIndex = -1;
2476 forAll(patchCoupling_, patchI)
2477 {
2478 if (patchCoupling_(patchI))
2479 {
2480 const coupleMap& cMap = patchCoupling_[patchI].map();
2482 if ((sIndex = cMap.findSlave(edgeEnum, eIndex)) > -1)
2483 {
2484 pIndex = patchI;
2486 break;
2487 }
2489 // The following bit happens only during the sliver
2490 // exudation process, since slave edges are
2491 // usually not added to the coupled edge-stack.
2492 if ((sIndex = cMap.findMaster(edgeEnum, eIndex)) > -1)
2493 {
2494 pIndex = patchI;
2496 // Notice that we are bisecting a slave edge first.
2497 bisectingSlave = true;
2499 break;
2500 }
2501 }
2502 }
2504 if (sIndex == -1)
2505 {
2506 FatalErrorIn
2507 (
2508 "\n"
2509 "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2510 "(\n"
2511 " const label eIndex,\n"
2512 " bool checkOnly,\n"
2513 " bool forceOp\n"
2514 ")\n"
2515 )
2516 << "Coupled maps were improperly specified." << nl
2517 << " Slave index not found for: " << nl
2518 << " Edge: " << eIndex << ": " << edges_[eIndex] << nl
2519 << abort(FatalError);
2520 }
2521 else
2522 {
2523 // If we've found the slave, size up the list
2524 meshOps::sizeUpList
2525 (
2526 changeMap(),
2527 slaveMaps
2528 );
2530 // Save index and patch for posterity
2531 slaveMaps[0].index() = sIndex;
2532 slaveMaps[0].patchIndex() = pIndex;
2533 }
2535 if (debug > 1)
2536 {
2537 Pout<< nl << " >> Bisecting slave edge: " << sIndex
2538 << " for master edge: " << eIndex << endl;
2539 }
2540 }
2541 else
2542 if (procCouple && !localCouple)
2543 {
2544 // If this is a new entity, bail out for now.
2545 // This will be handled at the next time-step.
2546 if (eIndex >= nOldEdges_)
2547 {
2548 return map;
2549 }
2551 // Check slaves
2552 forAll(procIndices_, pI)
2553 {
2554 // Fetch reference to subMeshes
2555 const coupledInfo& sendMesh = sendMeshes_[pI];
2556 const coupledInfo& recvMesh = recvMeshes_[pI];
2558 const coupleMap& scMap = sendMesh.map();
2559 const coupleMap& rcMap = recvMesh.map();
2561 // If this edge was sent to a lower-ranked
2562 // processor, skip it.
2563 if (procIndices_[pI] < Pstream::myProcNo())
2564 {
2565 if (scMap.reverseEntityMap(edgeEnum).found(eIndex))
2566 {
2567 if (debug > 3)
2568 {
2569 Pout<< "Edge: " << eIndex
2570 << "::" << eCheck
2571 << " was sent to proc: "
2572 << procIndices_[pI]
2573 << ", so bailing out."
2574 << endl;
2575 }
2577 return map;
2578 }
2579 }
2581 label sIndex = -1;
2583 if ((sIndex = rcMap.findSlave(edgeEnum, eIndex)) > -1)
2584 {
2585 if (debug > 3)
2586 {
2587 Pout<< "Checking slave edge: " << sIndex
2588 << " on proc: " << procIndices_[pI]
2589 << " for master edge: " << eIndex
2590 << endl;
2591 }
2593 // Check if a lower-ranked processor is
2594 // handling this edge
2595 if (procIndices_[pI] < Pstream::myProcNo())
2596 {
2597 if (debug > 3)
2598 {
2599 Pout<< "Edge: " << eIndex
2600 << " is handled by proc: "
2601 << procIndices_[pI]
2602 << ", so bailing out."
2603 << endl;
2604 }
2606 return map;
2607 }
2609 label curIndex = slaveMaps.size();
2611 // Size up the list
2612 meshOps::sizeUpList
2613 (
2614 changeMap(),
2615 slaveMaps
2616 );
2618 // Save index and patch for posterity
2619 slaveMaps[curIndex].index() = sIndex;
2620 slaveMaps[curIndex].patchIndex() = pI;
2621 }
2622 }
2623 }
2624 else
2625 {
2626 // Something's wrong with coupling maps
2627 FatalErrorIn
2628 (
2629 "\n"
2630 "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2631 "(\n"
2632 " const label eIndex,\n"
2633 " bool checkOnly,\n"
2634 " bool forceOp\n"
2635 ")\n"
2636 )
2637 << "Coupled maps were improperly specified." << nl
2638 << " localCouple: " << localCouple << nl
2639 << " procCouple: " << procCouple << nl
2640 << " Edge: " << eIndex << ": " << edges_[eIndex] << nl
2641 << abort(FatalError);
2642 }
2644 // Temporarily turn off coupledModification
2645 unsetCoupledModification();
2647 // First check the master for bisection feasibility.
2648 changeMap masterMap = bisectEdge(eIndex, true, forceOp);
2650 // Turn it back on.
2651 setCoupledModification();
2653 // Master couldn't perform bisection
2654 if (masterMap.type() <= 0)
2655 {
2656 return masterMap;
2657 }
2659 // Now check each of the slaves for bisection feasibility
2660 forAll(slaveMaps, slaveI)
2661 {
2662 // Alias for convenience...
2663 changeMap& slaveMap = slaveMaps[slaveI];
2665 label sIndex = slaveMap.index();
2666 label pI = slaveMap.patchIndex();
2668 // If the edge is mapped onto itself, skip check
2669 // (occurs for cyclic edges)
2670 if ((sIndex == eIndex) && localCouple)
2671 {
2672 continue;
2673 }
2675 // Temporarily turn off coupledModification
2676 unsetCoupledModification();
2678 // Test the slave edge
2679 if (localCouple)
2680 {
2681 slaveMap = bisectEdge(sIndex, true, forceOp);
2682 }
2683 else
2684 if (procCouple)
2685 {
2686 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
2688 if (debug > 3)
2689 {
2690 Pout<< "Checking slave edge: " << sIndex
2691 << "::" << sMesh.edges_[sIndex]
2692 << " on proc: " << procIndices_[pI]
2693 << " for master edge: " << eIndex
2694 << endl;
2695 }
2697 slaveMap = sMesh.bisectEdge(sIndex, true, forceOp);
2698 }
2700 // Turn it back on.
2701 setCoupledModification();
2703 if (slaveMap.type() <= 0)
2704 {
2705 // Slave couldn't perform bisection.
2706 map.type() = -2;
2708 return map;
2709 }
2711 // Save index and patch for posterity
2712 slaveMap.index() = sIndex;
2713 slaveMap.patchIndex() = pI;
2714 }
2716 // Next bisect each slave edge..
2717 forAll(slaveMaps, slaveI)
2718 {
2719 // Alias for convenience...
2720 changeMap& slaveMap = slaveMaps[slaveI];
2722 label sIndex = slaveMap.index();
2723 label pI = slaveMap.patchIndex();
2725 // If the edge is mapped onto itself, skip modification
2726 // (occurs for cyclic edges)
2727 if ((sIndex == eIndex) && localCouple)
2728 {
2729 continue;
2730 }
2732 // Temporarily turn off coupledModification
2733 unsetCoupledModification();
2735 // Bisect the slave edge
2736 if (localCouple)
2737 {
2738 slaveMap = bisectEdge(sIndex, false, forceOp);
2739 }
2740 else
2741 {
2742 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
2744 slaveMap = sMesh.bisectEdge(sIndex, false, forceOp);
2745 }
2747 // Turn it back on.
2748 setCoupledModification();
2750 // The final operation has to succeed.
2751 if (slaveMap.type() <= 0)
2752 {
2753 FatalErrorIn
2754 (
2755 "\n"
2756 "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2757 "(\n"
2758 " const label eIndex,\n"
2759 " bool checkOnly,\n"
2760 " bool forceOp\n"
2761 ")\n"
2762 )
2763 << "Coupled topo-change for slave failed." << nl
2764 << " Edge: " << eIndex << ": " << edges_[eIndex] << nl
2765 << " Slave index: " << sIndex << nl
2766 << " Patch index: " << pI << nl
2767 << " Type: " << slaveMap.type() << nl
2768 << abort(FatalError);
2769 }
2771 // Save index and patch for posterity
2772 slaveMap.index() = sIndex;
2773 slaveMap.patchIndex() = pI;
2774 }
2775 }
2777 // Before we bisect this edge, check whether the operation will
2778 // yield an acceptable cell-quality.
2779 scalar minQ = 0.0;
2781 if ((minQ = computeBisectionQuality(eIndex)) < sliverThreshold_)
2782 {
2783 // Check if the quality is actually valid before forcing it.
2784 if (forceOp && (minQ < 0.0))
2785 {
2786 FatalErrorIn
2787 (
2788 "\n"
2789 "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2790 "(\n"
2791 " const label eIndex,\n"
2792 " bool checkOnly,\n"
2793 " bool forceOp\n"
2794 ")\n"
2795 )
2796 << " Forcing bisection on edge: " << eIndex
2797 << " will yield an invalid cell."
2798 << abort(FatalError);
2799 }
2800 else
2801 if (!forceOp)
2802 {
2803 map.type() = -2;
2804 return map;
2805 }
2806 }
2808 // Are we performing only checks?
2809 if (checkOnly)
2810 {
2811 if (debug > 3 && isSubMesh_)
2812 {
2813 Pout<< " Slave edge: " << eIndex
2814 << " can be bisected."
2815 << endl;
2816 }
2818 map.type() = 1;
2819 return map;
2820 }
2822 // Update number of surface bisections, if necessary.
2823 if (whichEdgePatch(eIndex) > -1)
2824 {
2825 // Do not update stats for processor patches
2826 if (!processorCoupledEntity(eIndex))
2827 {
2828 statistics_[5]++;
2829 }
2830 }
2832 // Hull variables
2833 face tmpTriFace(3);
2834 edge origEdge(edges_[eIndex]);
2835 labelList tmpEdgeFaces(3,-1);
2836 labelList tmpIntEdgeFaces(4,-1);
2837 labelList tmpFaceEdges(3,-1);
2839 // Build vertexHull for this edge
2840 labelList vertexHull;
2841 buildVertexHull(eIndex, vertexHull);
2843 // Size up the hull lists
2844 label m = vertexHull.size();
2845 labelList cellHull(m, -1);
2846 labelList faceHull(m, -1);
2847 labelList edgeHull(m, -1);
2848 labelListList ringEntities(4, labelList(m, -1));
2850 // Construct a hull around this edge
2851 meshOps::constructHull
2852 (
2853 eIndex,
2854 faces_,
2855 edges_,
2856 cells_,
2857 owner_,
2858 neighbour_,
2859 faceEdges_,
2860 edgeFaces_,
2861 vertexHull,
2862 edgeHull,
2863 faceHull,
2864 cellHull,
2865 ringEntities
2866 );
2868 if (debug > 1)
2869 {
2870 Pout<< nl << nl
2871 << "Edge: " << eIndex
2872 << ": " << origEdge
2873 << " is to be bisected. " << endl;
2875 Pout<< " On SubMesh: " << Switch::asText(isSubMesh_) << nl;
2876 Pout<< " coupledModification: " << coupledModification_ << nl;
2878 label epIndex = whichEdgePatch(eIndex);
2880 const polyBoundaryMesh& boundary = boundaryMesh();
2882 Pout<< " Patch: ";
2884 if (epIndex == -1)
2885 {
2886 Pout<< "Internal" << endl;
2887 }
2888 else
2889 if (epIndex < boundary.size())
2890 {
2891 Pout<< boundary[epIndex].name() << endl;
2892 }
2893 else
2894 {
2895 Pout<< " New patch: " << epIndex << endl;
2896 }
2898 // Write out VTK files prior to change
2899 // - Using old-points for convenience in post-processing
2900 if (debug > 3)
2901 {
2902 writeVTK
2903 (
2904 Foam::name(eIndex)
2905 + '(' + Foam::name(origEdge[0])
2906 + ',' + Foam::name(origEdge[1]) + ')'
2907 + "_Bisect_0",
2908 cellHull,
2909 3, false, true
2910 );
2911 }
2912 }
2914 labelList mP(2, -1);
2916 // Set mapping for this point
2917 mP[0] = edges_[eIndex][0];
2918 mP[1] = edges_[eIndex][1];
2920 // Add a new point to the end of the list
2921 label newPointIndex =
2922 (
2923 insertPoint
2924 (
2925 0.5 * (points_[mP[0]] + points_[mP[1]]),
2926 0.5 * (oldPoints_[mP[0]] + oldPoints_[mP[1]]),
2927 mP
2928 )
2929 );
2931 // Add this point to the map.
2932 map.addPoint(newPointIndex);
2934 // New edges can lie on a bounding curve between
2935 // coupled and non-coupled faces. Preferentially
2936 // add edges to coupled-patches, if they exist.
2937 // This makes it convenient for coupled patch matching.
2938 label nePatch = -1;
2940 if (locallyCoupledEntity(eIndex, true))
2941 {
2942 nePatch = locallyCoupledEdgePatch(eIndex);
2943 }
2944 else
2945 {
2946 nePatch = whichEdgePatch(eIndex);
2947 }
2949 // Add a new edge to the end of the list
2950 label newEdgeIndex =
2951 (
2952 insertEdge
2953 (
2954 nePatch,
2955 edge(newPointIndex,edges_[eIndex][1]),
2956 labelList(faceHull.size(),-1)
2957 )
2958 );
2960 // Add this edge to the map.
2961 map.addEdge(newEdgeIndex);
2963 // Remove the existing edge from the pointEdges list
2964 // of the modified point, and add it to the new point
2965 meshOps::sizeDownList(eIndex, pointEdges_[edges_[eIndex][1]]);
2966 meshOps::sizeUpList(eIndex, pointEdges_[newPointIndex]);
2968 // Modify the existing edge
2969 edges_[eIndex][1] = newPointIndex;
2971 // Add this edge to the map.
2972 // Although this edge isn't technically 'added', it's
2973 // required for coupled patch mapping.
2974 map.addEdge(eIndex);
2976 // Keep track of added entities
2977 labelList addedCellIndices(cellHull.size(),-1);
2978 labelList addedFaceIndices(faceHull.size(),-1);
2979 labelList addedEdgeIndices(faceHull.size(),-1);
2980 labelList addedIntFaceIndices(faceHull.size(),-1);
2982 // Now loop through the hull and bisect individual entities
2983 forAll(vertexHull, indexI)
2984 {
2985 // Modify the existing face
2986 meshOps::replaceLabel
2987 (
2988 edges_[newEdgeIndex][1],
2989 newPointIndex,
2990 faces_[faceHull[indexI]]
2991 );
2993 // Obtain circular indices
2994 label nextI = vertexHull.fcIndex(indexI);
2995 label prevI = vertexHull.rcIndex(indexI);
2997 // Check if this is an interior/boundary face
2998 if (cellHull[indexI] != -1)
2999 {
3000 // Create a new cell. Add it for now, but update later.
3001 cell newCell(4);
3003 addedCellIndices[indexI] =
3004 (
3005 insertCell
3006 (
3007 newCell,
3008 lengthScale_[cellHull[indexI]]
3009 )
3010 );
3012 // Add this cell to the map.
3013 map.addCell(addedCellIndices[indexI]);
3015 // Configure the interior face
3016 tmpTriFace[0] = vertexHull[nextI];
3017 tmpTriFace[1] = vertexHull[indexI];
3018 tmpTriFace[2] = newPointIndex;
3020 // Insert the face
3021 addedIntFaceIndices[indexI] =
3022 (
3023 insertFace
3024 (
3025 -1,
3026 tmpTriFace,
3027 cellHull[indexI],
3028 addedCellIndices[indexI]
3029 )
3030 );
3032 // Add a faceEdges entry as well
3033 faceEdges_.append(tmpFaceEdges);
3035 // Add this face to the map.
3036 map.addFace(addedIntFaceIndices[indexI]);
3038 // Add to the new cell
3039 newCell[0] = addedIntFaceIndices[indexI];
3041 // Modify the existing ring face connected to newEdge[1]
3042 label replaceFace = ringEntities[3][indexI];
3044 // Check if face reversal is necessary
3045 if (owner_[replaceFace] == cellHull[indexI])
3046 {
3047 if (neighbour_[replaceFace] == -1)
3048 {
3049 // Change the owner
3050 owner_[replaceFace] = addedCellIndices[indexI];
3051 }
3052 else
3053 {
3054 // This face has to be reversed
3055 faces_[replaceFace] = faces_[replaceFace].reverseFace();
3056 owner_[replaceFace] = neighbour_[replaceFace];
3057 neighbour_[replaceFace] = addedCellIndices[indexI];
3059 setFlip(replaceFace);
3060 }
3061 }
3062 else
3063 {
3064 // Keep owner, but change neighbour
3065 neighbour_[replaceFace] = addedCellIndices[indexI];
3066 }
3068 // Modify the edge on the ring.
3069 // Add the new interior face to edgeFaces.
3070 meshOps::sizeUpList
3071 (
3072 addedIntFaceIndices[indexI],
3073 edgeFaces_[edgeHull[indexI]]
3074 );
3076 // Add this edge to faceEdges for the new interior face
3077 faceEdges_[addedIntFaceIndices[indexI]][0] = edgeHull[indexI];
3079 // Replace face labels
3080 meshOps::replaceLabel
3081 (
3082 replaceFace,
3083 addedIntFaceIndices[indexI],
3084 cells_[cellHull[indexI]]
3085 );
3087 // Add to the new cell
3088 newCell[1] = replaceFace;
3090 // Check if this is a boundary face
3091 if (cellHull[prevI] == -1)
3092 {
3093 // Configure the boundary face
3094 tmpTriFace[0] = newPointIndex;
3095 tmpTriFace[1] = edges_[newEdgeIndex][1];
3096 tmpTriFace[2] = vertexHull[indexI];
3098 // Insert the face
3099 addedFaceIndices[indexI] =
3100 (
3101 insertFace
3102 (
3103 whichPatch(faceHull[indexI]),
3104 tmpTriFace,
3105 addedCellIndices[indexI],
3106 -1
3107 )
3108 );
3110 // Add this face to the map.
3111 map.addFace(addedFaceIndices[indexI]);
3113 // Configure edgeFaces
3114 tmpEdgeFaces[0] = faceHull[indexI];
3115 tmpEdgeFaces[1] = addedIntFaceIndices[indexI];
3116 tmpEdgeFaces[2] = addedFaceIndices[indexI];
3118 // Add an edge
3119 addedEdgeIndices[indexI] =
3120 (
3121 insertEdge
3122 (
3123 whichPatch(faceHull[indexI]),
3124 edge(newPointIndex,vertexHull[indexI]),
3125 tmpEdgeFaces
3126 )
3127 );
3129 // Add this edge to the map.
3130 map.addEdge(addedEdgeIndices[indexI]);
3132 // Add this edge to the interior-face faceEdges entry
3133 faceEdges_[addedIntFaceIndices[indexI]][1] =
3134 (
3135 addedEdgeIndices[indexI]
3136 );
3138 // Configure faceEdges for this boundary face
3139 tmpFaceEdges[0] = addedEdgeIndices[indexI];
3140 tmpFaceEdges[1] = newEdgeIndex;
3141 tmpFaceEdges[2] = ringEntities[2][indexI];
3143 // Modify faceEdges for the hull face
3144 meshOps::replaceLabel
3145 (
3146 ringEntities[2][indexI],
3147 addedEdgeIndices[indexI],
3148 faceEdges_[faceHull[indexI]]
3149 );
3151 // Modify edgeFaces for the edge connected to newEdge[1]
3152 meshOps::replaceLabel
3153 (
3154 faceHull[indexI],
3155 addedFaceIndices[indexI],
3156 edgeFaces_[ringEntities[2][indexI]]
3157 );
3159 // Add the faceEdges entry
3160 faceEdges_.append(tmpFaceEdges);
3162 // Add an entry to edgeFaces
3163 edgeFaces_[newEdgeIndex][indexI] = addedFaceIndices[indexI];
3165 // Add an entry for this cell
3166 newCell[2] = addedFaceIndices[indexI];
3167 }
3168 else
3169 // Check if a cell was added before this
3170 if (addedCellIndices[prevI] != -1)
3171 {
3172 // Configure the interior face
3173 tmpTriFace[0] = vertexHull[indexI];
3174 tmpTriFace[1] = edges_[newEdgeIndex][1];
3175 tmpTriFace[2] = newPointIndex;
3177 // Insert the face
3178 addedFaceIndices[indexI] =
3179 (
3180 insertFace
3181 (
3182 -1,
3183 tmpTriFace,
3184 addedCellIndices[prevI],
3185 addedCellIndices[indexI]
3186 )
3187 );
3189 // Add this face to the map.
3190 map.addFace(addedFaceIndices[indexI]);
3192 // Configure edgeFaces
3193 tmpIntEdgeFaces[0] = faceHull[indexI];
3194 tmpIntEdgeFaces[1] = addedIntFaceIndices[indexI];
3195 tmpIntEdgeFaces[2] = addedFaceIndices[indexI];
3196 tmpIntEdgeFaces[3] = addedIntFaceIndices[prevI];
3198 // Add an internal edge
3199 addedEdgeIndices[indexI] =
3200 (
3201 insertEdge
3202 (
3203 -1,
3204 edge(newPointIndex,vertexHull[indexI]),
3205 tmpIntEdgeFaces
3206 )
3207 );
3209 // Add this edge to the map.
3210 map.addEdge(addedEdgeIndices[indexI]);
3212 // Add this edge to the interior-face faceEdges entry..
3213 faceEdges_[addedIntFaceIndices[indexI]][1] =
3214 (
3215 addedEdgeIndices[indexI]
3216 );
3218 // ... and to the previous interior face as well
3219 faceEdges_[addedIntFaceIndices[prevI]][2] =
3220 (
3221 addedEdgeIndices[indexI]
3222 );
3224 // Configure faceEdges for this split interior face
3225 tmpFaceEdges[0] = addedEdgeIndices[indexI];
3226 tmpFaceEdges[1] = newEdgeIndex;
3227 tmpFaceEdges[2] = ringEntities[2][indexI];
3229 // Modify faceEdges for the hull face
3230 meshOps::replaceLabel
3231 (
3232 ringEntities[2][indexI],
3233 addedEdgeIndices[indexI],
3234 faceEdges_[faceHull[indexI]]
3235 );
3237 // Modify edgeFaces for the edge connected to newEdge[1]
3238 meshOps::replaceLabel
3239 (
3240 faceHull[indexI],
3241 addedFaceIndices[indexI],
3242 edgeFaces_[ringEntities[2][indexI]]
3243 );
3245 // Add the faceEdges entry
3246 faceEdges_.append(tmpFaceEdges);
3248 // Add an entry to edgeFaces
3249 edgeFaces_[newEdgeIndex][indexI] = addedFaceIndices[indexI];
3251 // Add an entry for this cell
3252 newCell[2] = addedFaceIndices[indexI];
3254 // Make the final entry for the previous cell
3255 cells_[addedCellIndices[prevI]][3] = addedFaceIndices[indexI];
3256 }
3258 // Do the first interior face at the end
3259 if (indexI == vertexHull.size() - 1)
3260 {
3261 // Configure the interior face
3262 tmpTriFace[0] = newPointIndex;
3263 tmpTriFace[1] = edges_[newEdgeIndex][1];
3264 tmpTriFace[2] = vertexHull[0];
3266 // Insert the face
3267 addedFaceIndices[0] =
3268 (
3269 insertFace
3270 (
3271 -1,
3272 tmpTriFace,
3273 addedCellIndices[0],
3274 addedCellIndices[indexI]
3275 )
3276 );
3278 // Add this face to the map.
3279 map.addFace(addedFaceIndices[0]);
3281 // Configure edgeFaces
3282 tmpIntEdgeFaces[0] = faceHull[0];
3283 tmpIntEdgeFaces[1] = addedIntFaceIndices[0];
3284 tmpIntEdgeFaces[2] = addedFaceIndices[0];
3285 tmpIntEdgeFaces[3] = addedIntFaceIndices[indexI];
3287 // Add an internal edge
3288 addedEdgeIndices[0] =
3289 (
3290 insertEdge
3291 (
3292 -1,
3293 edge(newPointIndex,vertexHull[0]),
3294 tmpIntEdgeFaces
3295 )
3296 );
3298 // Add this edge to the map.
3299 map.addEdge(addedEdgeIndices[0]);
3301 // Add this edge to the interior-face faceEdges entry..
3302 faceEdges_[addedIntFaceIndices[0]][1] =
3303 (
3304 addedEdgeIndices[0]
3305 );
3307 // ... and to the previous interior face as well
3308 faceEdges_[addedIntFaceIndices[indexI]][2] =
3309 (
3310 addedEdgeIndices[0]
3311 );
3313 // Configure faceEdges for the first split face
3314 tmpFaceEdges[0] = addedEdgeIndices[0];
3315 tmpFaceEdges[1] = newEdgeIndex;
3316 tmpFaceEdges[2] = ringEntities[2][0];
3318 // Modify faceEdges for the hull face
3319 meshOps::replaceLabel
3320 (
3321 ringEntities[2][0],
3322 addedEdgeIndices[0],
3323 faceEdges_[faceHull[0]]
3324 );
3326 // Modify edgeFaces for the edge connected to newEdge[1]
3327 meshOps::replaceLabel
3328 (
3329 faceHull[0],
3330 addedFaceIndices[0],
3331 edgeFaces_[ringEntities[2][0]]
3332 );
3334 // Add the faceEdges entry
3335 faceEdges_.append(tmpFaceEdges);
3337 // Add an entry to edgeFaces
3338 edgeFaces_[newEdgeIndex][0] = addedFaceIndices[0];
3340 // Add an entry for this cell
3341 newCell[3] = addedFaceIndices[0];
3343 // Make the final entry for the first cell
3344 cells_[addedCellIndices[0]][2] = addedFaceIndices[0];
3345 }
3347 // Update the cell list with the new cell.
3348 cells_[addedCellIndices[indexI]] = newCell;
3349 }
3350 else
3351 {
3352 // Configure the final boundary face
3353 tmpTriFace[0] = vertexHull[indexI];
3354 tmpTriFace[1] = edges_[newEdgeIndex][1];
3355 tmpTriFace[2] = newPointIndex;
3357 // Insert the face
3358 addedFaceIndices[indexI] =
3359 (
3360 insertFace
3361 (
3362 whichPatch(faceHull[indexI]),
3363 tmpTriFace,
3364 addedCellIndices[prevI],
3365 -1
3366 )
3367 );
3369 // Add this face to the map.
3370 map.addFace(addedFaceIndices[indexI]);
3372 // Configure edgeFaces
3373 tmpEdgeFaces[0] = addedFaceIndices[indexI];
3374 tmpEdgeFaces[1] = addedIntFaceIndices[prevI];
3375 tmpEdgeFaces[2] = faceHull[indexI];
3377 // Add an edge
3378 addedEdgeIndices[indexI] =
3379 (
3380 insertEdge
3381 (
3382 whichPatch(faceHull[indexI]),
3383 edge(newPointIndex,vertexHull[indexI]),
3384 tmpEdgeFaces
3385 )
3386 );
3388 // Add this edge to the map.
3389 map.addEdge(addedEdgeIndices[indexI]);
3391 // Add a faceEdges entry to the previous interior face
3392 faceEdges_[addedIntFaceIndices[prevI]][2] =
3393 (
3394 addedEdgeIndices[indexI]
3395 );
3397 // Configure faceEdges for the final boundary face
3398 tmpFaceEdges[0] = addedEdgeIndices[indexI];
3399 tmpFaceEdges[1] = newEdgeIndex;
3400 tmpFaceEdges[2] = ringEntities[2][indexI];
3402 // Modify faceEdges for the hull face
3403 meshOps::replaceLabel
3404 (
3405 ringEntities[2][indexI],
3406 addedEdgeIndices[indexI],
3407 faceEdges_[faceHull[indexI]]
3408 );
3410 // Modify edgeFaces for the edge connected to newEdge[1]
3411 meshOps::replaceLabel
3412 (
3413 faceHull[indexI],
3414 addedFaceIndices[indexI],
3415 edgeFaces_[ringEntities[2][indexI]]
3416 );
3418 // Add the faceEdges entry
3419 faceEdges_.append(tmpFaceEdges);
3421 // Add an entry to edgeFaces
3422 edgeFaces_[newEdgeIndex][indexI] = addedFaceIndices[indexI];
3424 // Make the final entry for the previous cell
3425 cells_[addedCellIndices[prevI]][3] = addedFaceIndices[indexI];
3426 }
3427 }
3429 // Now that all old / new cells possess correct connectivity,
3430 // compute mapping information.
3431 forAll(cellHull, indexI)
3432 {
3433 if (cellHull[indexI] == -1)
3434 {
3435 continue;
3436 }
3438 // Set mapping for both new and modified cells.
3439 FixedList<label, 2> cmIndex;
3441 cmIndex[0] = cellHull[indexI];
3442 cmIndex[1] = addedCellIndices[indexI];
3444 // Fill-in candidate mapping information
3445 labelList mC(1, cellHull[indexI]);
3447 forAll(cmIndex, cmI)
3448 {
3449 // Set the mapping for this cell
3450 setCellMapping(cmIndex[cmI], mC);
3451 }
3452 }
3454 // Set mapping information for old / new faces
3455 forAll(faceHull, indexI)
3456 {
3457 // Interior faces get default mapping
3458 if (addedIntFaceIndices[indexI] > -1)
3459 {
3460 setFaceMapping(addedIntFaceIndices[indexI]);
3461 }
3463 // Decide between default / weighted mapping
3464 // based on boundary information
3465 if (whichPatch(faceHull[indexI]) == -1)
3466 {
3467 // Interior faces get default mapping
3468 setFaceMapping(faceHull[indexI]);
3469 setFaceMapping(addedFaceIndices[indexI]);
3470 }
3471 else
3472 {
3473 // Compute mapping weights for boundary faces
3474 FixedList<label, 2> fmIndex;
3476 fmIndex[0] = faceHull[indexI];
3477 fmIndex[1] = addedFaceIndices[indexI];
3479 // Fill-in candidate mapping information
3480 labelList mF(1, faceHull[indexI]);
3482 forAll(fmIndex, fmI)
3483 {
3484 // Set the mapping for this face
3485 setFaceMapping(fmIndex[fmI], mF);
3486 }
3487 }
3488 }
3490 // If modification is coupled, update mapping info.
3491 if (coupledModification_)
3492 {
3493 // Build a list of boundary edges / faces for mapping
3494 DynamicList<label> checkEdges(8), checkFaces(4);
3496 const labelList& oeFaces = edgeFaces_[eIndex];
3497 const labelList& neFaces = edgeFaces_[newEdgeIndex];
3499 forAll(oeFaces, faceI)
3500 {
3501 FixedList<label, 2> check(-1);
3503 check[0] = oeFaces[faceI];
3504 check[1] = neFaces[faceI];
3506 forAll(check, indexI)
3507 {
3508 label cPatch = whichPatch(check[indexI]);
3510 if (localCouple && !procCouple)
3511 {
3512 if (!locallyCoupledEntity(check[indexI], true, false, true))
3513 {
3514 continue;
3515 }
3517 // Check for cyclics
3518 if (boundaryMesh()[cPatch].type() == "cyclic")
3519 {
3520 // Check if this is a master face
3521 const coupleMap& cM = patchCoupling_[cPatch].map();
3522 const Map<label>& fM = cM.entityMap(coupleMap::FACE);
3524 // Only add master faces
3525 // (belonging to the first half)
3526 // - Only check[0] will be present,
3527 // so check for that alone
3528 if (!fM.found(check[0]))
3529 {
3530 continue;
3531 }
3532 }
3533 }
3534 else
3535 if (procCouple && !localCouple)
3536 {
3537 if (getNeighbourProcessor(cPatch) == -1)
3538 {
3539 continue;
3540 }
3541 }
3543 // Add face and its edges for checking
3544 if (findIndex(checkFaces, check[indexI]) == -1)
3545 {
3546 // Add this face
3547 checkFaces.append(check[indexI]);
3549 const labelList& fEdges = faceEdges_[check[indexI]];
3551 forAll(fEdges, edgeI)
3552 {
3553 if (findIndex(checkEdges, fEdges[edgeI]) == -1)
3554 {
3555 checkEdges.append(fEdges[edgeI]);
3556 }
3557 }
3558 }
3559 }
3560 }
3562 // Output check entities
3563 if (debug > 4)
3564 {
3565 writeVTK
3566 (
3567 "checkEdges_" + Foam::name(eIndex),
3568 checkEdges, 1, false, true
3569 );
3571 writeVTK
3572 (
3573 "checkFaces_" + Foam::name(eIndex),
3574 checkFaces, 2, false, true
3575 );
3576 }
3578 forAll(slaveMaps, slaveI)
3579 {
3580 // Alias for convenience...
3581 const changeMap& slaveMap = slaveMaps[slaveI];
3583 const label pI = slaveMap.patchIndex();
3585 // Fetch the appropriate coupleMap / mesh
3586 const coupleMap* cMapPtr = NULL;
3587 const dynamicTopoFvMesh* sMeshPtr = NULL;
3589 if (localCouple && !procCouple)
3590 {
3591 cMapPtr = &(patchCoupling_[pI].map());
3592 sMeshPtr = this;
3593 }
3594 else
3595 if (procCouple && !localCouple)
3596 {
3597 cMapPtr = &(recvMeshes_[pI].map());
3598 sMeshPtr = &(recvMeshes_[pI].subMesh());
3599 }
3601 // Alias for convenience
3602 const coupleMap& cMap = *cMapPtr;
3603 const dynamicTopoFvMesh& sMesh = *sMeshPtr;
3605 // Add the new points to the coupling map
3606 const List<objectMap>& apList = slaveMap.addedPointList();
3608 // Fetch the slave point
3609 label slavePoint = -1;
3611 if ((slaveMap.index() == eIndex) && localCouple)
3612 {
3613 // Cyclic edge at axis
3614 slavePoint = newPointIndex;
3615 }
3616 else
3617 {
3618 slavePoint = apList[0].index();
3619 }
3621 // Map points
3622 if (bisectingSlave)
3623 {
3624 // Update reverse pointMap
3625 cMap.mapMaster
3626 (
3627 coupleMap::POINT,
3628 newPointIndex,
3629 slavePoint
3630 );
3632 // Update pointMap
3633 cMap.mapSlave
3634 (
3635 coupleMap::POINT,
3636 slavePoint,
3637 newPointIndex
3638 );
3639 }
3640 else
3641 {
3642 // Update pointMap
3643 cMap.mapSlave
3644 (
3645 coupleMap::POINT,
3646 newPointIndex,
3647 slavePoint
3648 );
3650 // Update reverse pointMap
3651 cMap.mapMaster
3652 (
3653 coupleMap::POINT,
3654 slavePoint,
3655 newPointIndex
3656 );
3657 }
3659 if (debug > 2)
3660 {
3661 Pout<< " Adding point: " << slavePoint
3662 << " for point: " << newPointIndex
3663 << endl;
3664 }
3666 // Obtain point maps
3667 const Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
3669 // Update face mapping
3670 const label faceEnum = coupleMap::FACE;
3672 // Obtain references
3673 Map<label>& faceMap = cMap.entityMap(faceEnum);
3674 Map<label>& rFaceMap = cMap.reverseEntityMap(faceEnum);
3676 forAll(checkFaces, faceI)
3677 {
3678 label mfIndex = checkFaces[faceI];
3680 const face& mF = faces_[mfIndex];
3682 label mfPatch = whichPatch(mfIndex);
3684 // Check for processor match
3685 label neiProc = getNeighbourProcessor(mfPatch);
3687 if (procCouple && !localCouple)
3688 {
3689 if (neiProc != procIndices_[pI])
3690 {
3691 continue;
3692 }
3693 }
3695 triFace cF
3696 (
3697 pointMap.found(mF[0]) ? pointMap[mF[0]] : -1,
3698 pointMap.found(mF[1]) ? pointMap[mF[1]] : -1,
3699 pointMap.found(mF[2]) ? pointMap[mF[2]] : -1
3700 );
3702 // Skip mapping if all points were not found
3703 if (cF[0] == -1 || cF[1] == -1 || cF[2] == -1)
3704 {
3705 continue;
3706 }
3708 bool matchedFace = false;
3710 // Fetch edges connected to the slave point
3711 const labelList& spEdges = sMesh.pointEdges_[slavePoint];
3713 forAll(spEdges, edgeI)
3714 {
3715 label seIndex = spEdges[edgeI];
3717 if (sMesh.whichEdgePatch(seIndex) == -1)
3718 {
3719 continue;
3720 }
3722 const labelList& seFaces = sMesh.edgeFaces_[seIndex];
3724 forAll(seFaces, faceI)
3725 {
3726 label sfIndex = seFaces[faceI];
3728 if (sMesh.whichPatch(sfIndex) == -1)
3729 {
3730 continue;
3731 }
3733 const face& sF = sMesh.faces_[sfIndex];
3735 if
3736 (
3737 triFace::compare
3738 (
3739 triFace(sF[0], sF[1], sF[2]), cF
3740 )
3741 )
3742 {
3743 if (debug > 2)
3744 {
3745 word pN(boundaryMesh()[mfPatch].name());
3747 Pout<< " Found face: " << sfIndex
3748 << " :: " << sF
3749 << " with mfIndex: " << mfIndex
3750 << " :: " << mF
3751 << " Patch: " << pN
3752 << endl;
3753 }
3755 if (rFaceMap.found(sfIndex))
3756 {
3757 rFaceMap[sfIndex] = mfIndex;
3758 }
3759 else
3760 {
3761 rFaceMap.insert(sfIndex, mfIndex);
3762 }
3764 if (faceMap.found(mfIndex))
3765 {
3766 faceMap[mfIndex] = sfIndex;
3767 }
3768 else
3769 {
3770 faceMap.insert(mfIndex, sfIndex);
3771 }
3773 matchedFace = true;
3775 break;
3776 }
3777 }
3779 if (matchedFace)
3780 {
3781 break;
3782 }
3783 }
3785 if (!matchedFace)
3786 {
3787 sMesh.writeVTK
3788 (
3789 "failedFacePoints_"
3790 + Foam::name(mfIndex),
3791 cF, 0, false, true
3792 );
3794 writeVTK
3795 (
3796 "checkFaces_" + Foam::name(eIndex),
3797 checkFaces, 2, false, true
3798 );
3800 Pout<< " Failed to match face: "
3801 << mfIndex << " :: " << mF
3802 << " masterPatch: " << mfPatch
3803 << " using comparison face: " << cF
3804 << abort(FatalError);
3805 }
3806 }
3808 // Update edge mapping
3809 const label edgeEnum = coupleMap::EDGE;
3811 // Obtain references
3812 Map<label>& edgeMap = cMap.entityMap(edgeEnum);
3813 Map<label>& rEdgeMap = cMap.reverseEntityMap(edgeEnum);
3815 forAll(checkEdges, edgeI)
3816 {
3817 label meIndex = checkEdges[edgeI];
3819 const edge& mE = edges_[meIndex];
3821 label mePatch = whichPatch(meIndex);
3822 label neiProc = getNeighbourProcessor(mePatch);
3824 edge cE
3825 (
3826 pointMap.found(mE[0]) ? pointMap[mE[0]] : -1,
3827 pointMap.found(mE[1]) ? pointMap[mE[1]] : -1
3828 );
3830 // Skip mapping if all points were not found
3831 if (cE[0] == -1 || cE[1] == -1)
3832 {
3833 continue;
3834 }
3836 bool matchedEdge = false;
3838 // Fetch edges connected to the slave point
3839 const labelList& spEdges = sMesh.pointEdges_[cE[0]];
3841 forAll(spEdges, edgeI)
3842 {
3843 label seIndex = spEdges[edgeI];
3845 const edge& sE = sMesh.edges_[seIndex];
3847 if (sE == cE)
3848 {
3849 if (debug > 2)
3850 {
3851 Pout<< " Found edge: " << seIndex
3852 << " :: " << sE
3853 << " with meIndex: " << meIndex
3854 << " :: " << mE
3855 << endl;
3856 }
3858 // Update reverse map
3859 if (rEdgeMap.found(seIndex))
3860 {
3861 rEdgeMap[seIndex] = meIndex;
3862 }
3863 else
3864 {
3865 rEdgeMap.insert(seIndex, meIndex);
3866 }
3868 // Update map
3869 if (edgeMap.found(meIndex))
3870 {
3871 edgeMap[meIndex] = seIndex;
3872 }
3873 else
3874 {
3875 edgeMap.insert(meIndex, seIndex);
3876 }
3878 matchedEdge = true;
3880 break;
3881 }
3882 }
3884 if (!matchedEdge)
3885 {
3886 if (procCouple && !localCouple)
3887 {
3888 // Rare occassion where both points
3889 // of the edge lie on processor, but
3890 // not the edge itself.
3891 if (neiProc != procIndices_[pI])
3892 {
3893 if (debug > 2)
3894 {
3895 Pout<< " Edge: " << meIndex
3896 << " :: " << mE
3897 << " with comparison: " << cE
3898 << " has points on processor: " << neiProc
3899 << " but no edge. Marking as matched."
3900 << endl;
3901 }
3903 matchedEdge = true;
3904 }
3905 }
3906 }
3908 if (!matchedEdge)
3909 {
3910 sMesh.writeVTK
3911 (
3912 "failedEdge_"
3913 + Foam::name(meIndex),
3914 cE, 0, false, true
3915 );
3917 writeVTK
3918 (
3919 "checkEdges_" + Foam::name(eIndex),
3920 checkEdges, 1, false, true
3921 );
3923 Pout<< " Failed to match edge: "
3924 << meIndex << " :: " << mE
3925 << " using comparison edge: " << cE
3926 << abort(FatalError);
3927 }
3928 }
3930 // Push operation into coupleMap
3931 if (procCouple && !localCouple)
3932 {
3933 cMap.pushOperation
3934 (
3935 slaveMap.index(),
3936 coupleMap::BISECTION
3937 );
3938 }
3939 }
3940 }
3942 if (debug > 2)
3943 {
3944 label bPatch = whichEdgePatch(eIndex);
3946 const polyBoundaryMesh& boundary = boundaryMesh();
3948 if (bPatch == -1)
3949 {
3950 Pout<< "Patch: Internal" << nl;
3951 }
3952 else
3953 if (bPatch < boundary.size())
3954 {
3955 Pout<< "Patch: " << boundary[bPatch].name() << nl;
3956 }
3957 else
3958 {
3959 Pout<< " New patch: " << bPatch << endl;
3960 }
3962 Pout<< "vertexHull: " << vertexHull << nl
3963 << "Edges: " << edgeHull << nl
3964 << "Faces: " << faceHull << nl
3965 << "Cells: " << cellHull << nl;
3967 Pout<< "Modified cells: " << nl;
3969 forAll(cellHull, cellI)
3970 {
3971 if (cellHull[cellI] == -1)
3972 {
3973 continue;
3974 }
3976 Pout<< cellHull[cellI] << ":: "
3977 << cells_[cellHull[cellI]]
3978 << nl;
3979 }
3981 Pout<< nl << "Added cells: " << nl;
3983 forAll(addedCellIndices, cellI)
3984 {
3985 if (addedCellIndices[cellI] == -1)
3986 {
3987 continue;
3988 }
3990 Pout<< addedCellIndices[cellI] << ":: "
3991 << cells_[addedCellIndices[cellI]] << nl
3992 << "lengthScale: " << lengthScale_[addedCellIndices[cellI]]
3993 << nl;
3994 }
3996 Pout<< nl << "Modified faces: " << nl;
3998 forAll(faceHull, faceI)
3999 {
4000 Pout<< faceHull[faceI] << ":: "
4001 << faces_[faceHull[faceI]] << ": "
4002 << owner_[faceHull[faceI]] << ": "
4003 << neighbour_[faceHull[faceI]] << " "
4004 << "faceEdges:: " << faceEdges_[faceHull[faceI]]
4005 << nl;
4006 }
4008 Pout<< nl << "Added faces: " << nl;
4010 forAll(addedFaceIndices, faceI)
4011 {
4012 Pout<< addedFaceIndices[faceI] << ":: "
4013 << faces_[addedFaceIndices[faceI]] << ": "
4014 << owner_[addedFaceIndices[faceI]] << ": "
4015 << neighbour_[addedFaceIndices[faceI]] << " "
4016 << "faceEdges:: " << faceEdges_[addedFaceIndices[faceI]]
4017 << nl;
4018 }
4020 forAll(addedIntFaceIndices, faceI)
4021 {
4022 if (addedIntFaceIndices[faceI] == -1)
4023 {
4024 continue;
4025 }
4027 Pout<< addedIntFaceIndices[faceI] << ":: "
4028 << faces_[addedIntFaceIndices[faceI]] << ": "
4029 << owner_[addedIntFaceIndices[faceI]] << ": "
4030 << neighbour_[addedIntFaceIndices[faceI]] << " "
4031 << "faceEdges:: " << faceEdges_[addedIntFaceIndices[faceI]]
4032 << nl;
4033 }
4035 Pout<< nl << "New edge:: " << newEdgeIndex
4036 << ": " << edges_[newEdgeIndex] << nl
4037 << " edgeFaces:: " << edgeFaces_[newEdgeIndex]
4038 << nl;
4040 Pout<< nl << "Added edges: " << nl;
4042 forAll(addedEdgeIndices, edgeI)
4043 {
4044 Pout<< addedEdgeIndices[edgeI]
4045 << ":: " << edges_[addedEdgeIndices[edgeI]] << nl
4046 << " edgeFaces:: " << edgeFaces_[addedEdgeIndices[edgeI]]
4047 << nl;
4048 }
4050 Pout<< "New Point:: " << newPointIndex << nl
4051 << "pointEdges:: " << pointEdges_[newPointIndex] << nl;
4053 // Flush buffer
4054 Pout<< endl;
4056 // Write out VTK files after change
4057 if (debug > 3)
4058 {
4059 labelList newHull(cellHull.size() + addedCellIndices.size(), 0);
4061 // Combine both lists into one.
4062 forAll(cellHull, i)
4063 {
4064 newHull[i] = cellHull[i];
4065 }
4067 label start = cellHull.size();
4069 for(label i = start; i < newHull.size(); i++)
4070 {
4071 newHull[i] = addedCellIndices[i - start];
4072 }
4074 // Write out VTK files after change
4075 // - Using old-points for convenience in post-processing
4076 writeVTK
4077 (
4078 Foam::name(eIndex)
4079 + '(' + Foam::name(origEdge[0])
4080 + ',' + Foam::name(origEdge[1]) + ')'
4081 + "_Bisect_1",
4082 newHull,
4083 3, false, true
4084 );
4085 }
4086 }
4088 // Set the flag
4089 topoChangeFlag_ = true;
4091 // Increment the counter
4092 statistics_[3]++;
4094 // Increment the number of modifications
4095 statistics_[0]++;
4097 // Specify that the operation was successful
4098 map.type() = 1;
4100 // Return the changeMap
4101 return map;
4102 }
4105 // Method for the trisection of a face in 3D
4106 // - Returns a changeMap with a type specifying:
4107 // 1: Trisection was successful
4108 // -1: Trisection failed since max number of topo-changes was reached.
4109 // -2: Trisection failed since resulting quality would be really bad.
4110 // - AddedPoint is the index of the newly added point.
4111 const changeMap dynamicTopoFvMesh::trisectFace
4112 (
4113 const label fIndex,
4114 bool checkOnly,
4115 bool forceOp
4116 )
4117 {
4118 // Face trisection performs the following operations:
4119 // [1] Add a point at middle of the face
4120 // [2] Remove the face and add three new faces in place.
4121 // [3] Add three cells for each trisected cell (remove the originals).
4122 // [4] Create one internal edge for each trisected cell.
4123 // [5] Create three edges for the trisected face.
4124 // [6] Create three internal faces for each trisected cell.
4125 // Update faceEdges and edgeFaces information.
4127 // Figure out which thread this is...
4128 label tIndex = self(), pIndex = -1;
4130 // Prepare the changeMaps
4131 changeMap map, slaveMap;
4132 bool trisectingSlave = false;
4134 if
4135 (
4136 (statistics_[0] > maxModifications_) &&
4137 (maxModifications_ > -1)
4138 )
4139 {
4140 // Reached the max allowable topo-changes.
4141 stack(tIndex).clear();
4143 return map;
4144 }
4146 // Sanity check: Is the index legitimate?
4147 if (fIndex < 0)
4148 {
4149 FatalErrorIn
4150 (
4151 "const changeMap dynamicTopoFvMesh::trisectFace\n"
4152 "(\n"
4153 " const label fIndex,\n"
4154 " bool checkOnly,\n"
4155 " bool forceOp\n"
4156 ")\n"
4157 )
4158 << " Invalid index: " << fIndex
4159 << abort(FatalError);
4160 }
4162 if (coupledModification_)
4163 {
4164 // Is this a locally coupled face (either master or slave)?
4165 if (locallyCoupledEntity(fIndex, true))
4166 {
4167 label sIndex = -1;
4169 // Loop through master/slave maps
4170 // and determine the coupled edge index.
4171 forAll(patchCoupling_, patchI)
4172 {
4173 if (!patchCoupling_(patchI))
4174 {
4175 continue;
4176 }
4178 const label faceEnum = coupleMap::FACE;
4179 const coupleMap& cMap = patchCoupling_[patchI].map();
4181 if ((sIndex = cMap.findSlave(faceEnum, fIndex)) > -1)
4182 {
4183 // Keep this index for master/slave mapping.
4184 pIndex = patchI;
4186 break;
4187 }
4189 // The following bit happens only during the sliver
4190 // exudation process.
4191 if ((sIndex = cMap.findMaster(faceEnum, fIndex)) > -1)
4192 {
4193 // Keep this index for master/slave mapping.
4194 pIndex = patchI;
4196 // Notice that we are trisecting a slave edge.
4197 trisectingSlave = true;
4199 break;
4200 }
4201 }
4203 if (sIndex == -1)
4204 {
4205 FatalErrorIn
4206 (
4207 "const changeMap dynamicTopoFvMesh::trisectFace\n"
4208 "(\n"
4209 " const label fIndex,\n"
4210 " bool checkOnly,\n"
4211 " bool forceOp\n"
4212 ")\n"
4213 )
4214 << "Coupled maps were improperly specified." << nl
4215 << " Slave index not found for: " << nl
4216 << " Face: " << fIndex << nl
4217 << abort(FatalError);
4218 }
4220 if (debug > 1)
4221 {
4222 Pout<< nl << "Trisecting slave face: " << sIndex
4223 << " for master face: " << fIndex << endl;
4224 }
4226 // Temporarily turn off coupledModification.
4227 unsetCoupledModification();
4229 // First check the slave for trisection feasibility.
4230 slaveMap = trisectFace(sIndex, true, forceOp);
4232 if (slaveMap.type() == 1)
4233 {
4234 // Can the master be trisected as well?
4235 changeMap masterMap = trisectFace(fIndex, true, forceOp);
4237 // Master couldn't perform trisection
4238 if (masterMap.type() != 1)
4239 {
4240 setCoupledModification();
4242 return masterMap;
4243 }
4245 // Trisect the slave face
4246 slaveMap = trisectFace(sIndex, false, forceOp);
4248 // The final operation has to succeed.
4249 if (slaveMap.type() <= 0)
4250 {
4251 FatalErrorIn
4252 (
4253 "const changeMap dynamicTopoFvMesh::trisectFace\n"
4254 "(\n"
4255 " const label fIndex,\n"
4256 " bool checkOnly,\n"
4257 " bool forceOp\n"
4258 ")\n"
4259 )
4260 << "Coupled topo-change for slave failed." << nl
4261 << " Type: " << slaveMap.type()
4262 << abort(FatalError);
4263 }
4264 }
4265 else
4266 {
4267 // Slave couldn't perform trisection.
4268 setCoupledModification();
4270 map.type() = -2;
4272 return map;
4273 }
4275 // Turn it back on.
4276 setCoupledModification();
4277 }
4278 else
4279 if (processorCoupledEntity(fIndex))
4280 {
4281 // Trisect face on the patchSubMesh.
4283 }
4284 }
4286 // Before we trisect this face, check whether the operation will
4287 // yield an acceptable cell-quality.
4288 scalar minQ = 0.0;
4290 if ((minQ = computeTrisectionQuality(fIndex)) < sliverThreshold_)
4291 {
4292 // Check if the quality is actually valid before forcing it.
4293 if (forceOp && (minQ < 0.0))
4294 {
4295 FatalErrorIn
4296 (
4297 "const changeMap dynamicTopoFvMesh::trisectFace\n"
4298 "(\n"
4299 " const label fIndex,\n"
4300 " bool checkOnly,\n"
4301 " bool forceOp\n"
4302 ")\n"
4303 )
4304 << " Forcing trisection on face: " << fIndex
4305 << " will yield an invalid cell."
4306 << abort(FatalError);
4307 }
4308 else
4309 if (!forceOp)
4310 {
4311 map.type() = -2;
4312 return map;
4313 }
4314 }
4316 // Are we performing only checks?
4317 if (checkOnly)
4318 {
4319 map.type() = 1;
4320 return map;
4321 }
4323 // Update number of surface bisections, if necessary.
4324 if (whichPatch(fIndex) > -1)
4325 {
4326 statistics_[5]++;
4327 }
4329 // Hull variables
4330 face tmpTriFace(3);
4331 labelList newTriEdgeFaces(3);
4332 labelList newQuadEdgeFaces(4);
4334 FixedList<label,2> apexPoint(-1);
4335 FixedList<face, 3> checkFace(face(3));
4336 FixedList<label,5> newEdgeIndex(-1);
4337 FixedList<label,9> newFaceIndex(-1);
4338 FixedList<label,6> newCellIndex(-1);
4339 FixedList<cell, 6> newTetCell(cell(4));
4340 FixedList<labelList, 9> newFaceEdges(labelList(3));
4342 // Counters for entities
4343 FixedList<label, 9> nE(0);
4344 FixedList<label, 6> nF(0);
4346 // Determine the two cells to be removed
4347 FixedList<label,2> cellsForRemoval;
4348 cellsForRemoval[0] = owner_[fIndex];
4349 cellsForRemoval[1] = neighbour_[fIndex];
4351 if (debug > 1)
4352 {
4353 Pout<< nl << nl
4354 << "Face: " << fIndex
4355 << ": " << faces_[fIndex]
4356 << " is to be trisected. " << endl;
4358 // Write out VTK files prior to change
4359 if (debug > 3)
4360 {
4361 labelList vtkCells;
4363 if (neighbour_[fIndex] == -1)
4364 {
4365 vtkCells.setSize(1);
4366 vtkCells[0] = owner_[fIndex];
4367 }
4368 else
4369 {
4370 vtkCells.setSize(2);
4371 vtkCells[0] = owner_[fIndex];
4372 vtkCells[1] = neighbour_[fIndex];
4373 }
4375 writeVTK
4376 (
4377 Foam::name(fIndex)
4378 + "Trisect_0",
4379 vtkCells
4380 );
4381 }
4382 }
4384 labelList mP(3, -1);
4386 // Fill in mapping information
4387 mP[0] = faces_[fIndex][0];
4388 mP[1] = faces_[fIndex][1];
4389 mP[2] = faces_[fIndex][2];
4391 // Add a new point to the end of the list
4392 scalar oT = (1.0/3.0);
4394 label newPointIndex =
4395 (
4396 insertPoint
4397 (
4398 oT * (points_[mP[0]] + points_[mP[1]] + points_[mP[2]]),
4399 oT * (oldPoints_[mP[0]] + oldPoints_[mP[1]] + oldPoints_[mP[2]]),
4400 mP
4401 )
4402 );
4404 // Add this point to the map.
4405 map.addPoint(newPointIndex);
4407 // Add three new cells to the end of the cell list
4408 for (label i = 0; i < 3; i++)
4409 {
4410 scalar parentScale = -1.0;
4412 if (edgeRefinement_)
4413 {
4414 parentScale = lengthScale_[cellsForRemoval[0]];
4415 }
4417 newCellIndex[i] = insertCell(newTetCell[i], parentScale);
4419 // Add cells to the map
4420 map.addCell(newCellIndex[i]);
4421 }
4423 // Find the apex point for this cell
4424 apexPoint[0] =
4425 (
4426 meshOps::tetApexPoint
4427 (
4428 owner_[fIndex],
4429 fIndex,
4430 faces_,
4431 cells_
4432 )
4433 );
4435 // Insert three new internal faces
4437 // First face: Owner: newCellIndex[0], Neighbour: newCellIndex[1]
4438 tmpTriFace[0] = newPointIndex;
4439 tmpTriFace[1] = faces_[fIndex][0];
4440 tmpTriFace[2] = apexPoint[0];
4442 newFaceIndex[0] =
4443 (
4444 insertFace
4445 (
4446 -1,
4447 tmpTriFace,
4448 newCellIndex[0],
4449 newCellIndex[1]
4450 )
4451 );
4453 // Second face: Owner: newCellIndex[1], Neighbour: newCellIndex[2]
4454 tmpTriFace[0] = newPointIndex;
4455 tmpTriFace[1] = faces_[fIndex][1];
4456 tmpTriFace[2] = apexPoint[0];
4458 newFaceIndex[1] =
4459 (
4460 insertFace
4461 (
4462 -1,
4463 tmpTriFace,
4464 newCellIndex[1],
4465 newCellIndex[2]
4466 )
4467 );
4469 // Third face: Owner: newCellIndex[0], Neighbour: newCellIndex[2]
4470 tmpTriFace[0] = newPointIndex;
4471 tmpTriFace[1] = apexPoint[0];
4472 tmpTriFace[2] = faces_[fIndex][2];
4474 newFaceIndex[2] =
4475 (
4476 insertFace
4477 (
4478 -1,
4479 tmpTriFace,
4480 newCellIndex[0],
4481 newCellIndex[2]
4482 )
4483 );
4485 // Add an entry to edgeFaces
4486 newTriEdgeFaces[0] = newFaceIndex[0];
4487 newTriEdgeFaces[1] = newFaceIndex[1];
4488 newTriEdgeFaces[2] = newFaceIndex[2];
4490 // Add a new internal edge to the mesh
4491 newEdgeIndex[0] =
4492 (
4493 insertEdge
4494 (
4495 -1,
4496 edge
4497 (
4498 newPointIndex,
4499 apexPoint[0]
4500 ),
4501 newTriEdgeFaces
4502 )
4503 );
4505 // Configure faceEdges with the new internal edge
4506 newFaceEdges[0][nE[0]++] = newEdgeIndex[0];
4507 newFaceEdges[1][nE[1]++] = newEdgeIndex[0];
4508 newFaceEdges[2][nE[2]++] = newEdgeIndex[0];
4510 // Add the newly created faces to cells
4511 newTetCell[0][nF[0]++] = newFaceIndex[0];
4512 newTetCell[0][nF[0]++] = newFaceIndex[2];
4513 newTetCell[1][nF[1]++] = newFaceIndex[0];
4514 newTetCell[1][nF[1]++] = newFaceIndex[1];
4515 newTetCell[2][nF[2]++] = newFaceIndex[1];
4516 newTetCell[2][nF[2]++] = newFaceIndex[2];
4518 // Define the three faces to check for orientation:
4519 checkFace[0][0] = faces_[fIndex][2];
4520 checkFace[0][1] = apexPoint[0];
4521 checkFace[0][2] = faces_[fIndex][0];
4523 checkFace[1][0] = faces_[fIndex][0];
4524 checkFace[1][1] = apexPoint[0];
4525 checkFace[1][2] = faces_[fIndex][1];
4527 checkFace[2][0] = faces_[fIndex][1];
4528 checkFace[2][1] = apexPoint[0];
4529 checkFace[2][2] = faces_[fIndex][2];
4531 // Check the orientation of faces on the first cell.
4532 forAll(cells_[owner_[fIndex]], faceI)
4533 {
4534 label faceIndex = cells_[owner_[fIndex]][faceI];
4536 if (faceIndex == fIndex)
4537 {
4538 continue;
4539 }
4541 const face& faceToCheck = faces_[faceIndex];
4542 label cellIndex = cellsForRemoval[0];
4543 label newIndex = -1;
4545 // Check against faces.
4546 if (triFace::compare(triFace(faceToCheck), triFace(checkFace[0])))
4547 {
4548 newIndex = newCellIndex[0];
4549 newTetCell[0][nF[0]++] = faceIndex;
4550 }
4551 else
4552 if (triFace::compare(triFace(faceToCheck), triFace(checkFace[1])))
4553 {
4554 newIndex = newCellIndex[1];
4555 newTetCell[1][nF[1]++] = faceIndex;
4556 }
4557 else
4558 if (triFace::compare(triFace(faceToCheck), triFace(checkFace[2])))
4559 {
4560 newIndex = newCellIndex[2];
4561 newTetCell[2][nF[2]++] = faceIndex;
4562 }
4563 else
4564 {
4565 // Something's terribly wrong.
4566 FatalErrorIn
4567 (
4568 "const changeMap dynamicTopoFvMesh::trisectFace\n"
4569 "(\n"
4570 " const label fIndex,\n"
4571 " bool checkOnly,\n"
4572 " bool forceOp\n"
4573 ")\n"
4574 )
4575 << "Failed to determine a face match."
4576 << abort(FatalError);
4577 }
4579 // Check if a face-flip is necessary
4580 if (owner_[faceIndex] == cellIndex)
4581 {
4582 if (neighbour_[faceIndex] == -1)
4583 {
4584 // Change the owner
4585 owner_[faceIndex] = newIndex;
4586 }
4587 else
4588 {
4589 // Flip this face
4590 faces_[faceIndex] = faceToCheck.reverseFace();
4591 owner_[faceIndex] = neighbour_[faceIndex];
4592 neighbour_[faceIndex] = newIndex;
4594 setFlip(faceIndex);
4595 }
4596 }
4597 else
4598 {
4599 // Flip is unnecessary. Just update neighbour
4600 neighbour_[faceIndex] = newIndex;
4601 }
4602 }
4604 if (cellsForRemoval[1] == -1)
4605 {
4606 // Boundary face. Determine its patch.
4607 label facePatch = whichPatch(fIndex);
4609 // Add three new boundary faces.
4611 // Fourth face: Owner: newCellIndex[0], Neighbour: -1
4612 tmpTriFace[0] = newPointIndex;
4613 tmpTriFace[1] = faces_[fIndex][2];
4614 tmpTriFace[2] = faces_[fIndex][0];
4616 newFaceIndex[3] =
4617 (
4618 insertFace
4619 (
4620 facePatch,
4621 tmpTriFace,
4622 newCellIndex[0],
4623 -1
4624 )
4625 );
4627 // Fifth face: Owner: newCellIndex[1], Neighbour: -1
4628 tmpTriFace[0] = newPointIndex;
4629 tmpTriFace[1] = faces_[fIndex][0];
4630 tmpTriFace[2] = faces_[fIndex][1];
4632 newFaceIndex[4] =
4633 (
4634 insertFace
4635 (
4636 facePatch,
4637 tmpTriFace,
4638 newCellIndex[1],
4639 -1
4640 )
4641 );
4643 // Sixth face: Owner: newCellIndex[2], Neighbour: -1
4644 tmpTriFace[0] = newPointIndex;
4645 tmpTriFace[1] = faces_[fIndex][1];
4646 tmpTriFace[2] = faces_[fIndex][2];
4648 newFaceIndex[5] =
4649 (
4650 insertFace
4651 (
4652 facePatch,
4653 tmpTriFace,
4654 newCellIndex[2],
4655 -1
4656 )
4657 );
4659 // Add the newly created faces to cells
4660 newTetCell[0][nF[0]++] = newFaceIndex[3];
4661 newTetCell[1][nF[1]++] = newFaceIndex[4];
4662 newTetCell[2][nF[2]++] = newFaceIndex[5];
4664 // Configure edgeFaces for three new boundary edges.
4665 newTriEdgeFaces[0] = newFaceIndex[4];
4666 newTriEdgeFaces[1] = newFaceIndex[0];
4667 newTriEdgeFaces[2] = newFaceIndex[3];
4669 newEdgeIndex[1] =
4670 (
4671 insertEdge
4672 (
4673 facePatch,
4674 edge
4675 (
4676 newPointIndex,
4677 faces_[fIndex][0]
4678 ),
4679 newTriEdgeFaces
4680 )
4681 );
4683 newTriEdgeFaces[0] = newFaceIndex[5];
4684 newTriEdgeFaces[1] = newFaceIndex[1];
4685 newTriEdgeFaces[2] = newFaceIndex[4];
4687 newEdgeIndex[2] =
4688 (
4689 insertEdge
4690 (
4691 facePatch,
4692 edge
4693 (
4694 newPointIndex,
4695 faces_[fIndex][1]
4696 ),
4697 newTriEdgeFaces
4698 )
4699 );
4701 newTriEdgeFaces[0] = newFaceIndex[3];
4702 newTriEdgeFaces[1] = newFaceIndex[2];
4703 newTriEdgeFaces[2] = newFaceIndex[5];
4705 newEdgeIndex[3] =
4706 (
4707 insertEdge
4708 (
4709 facePatch,
4710 edge
4711 (
4712 newPointIndex,
4713 faces_[fIndex][2]
4714 ),
4715 newTriEdgeFaces
4716 )
4717 );
4719 // Configure faceEdges with the three new edges.
4720 newFaceEdges[0][nE[0]++] = newEdgeIndex[1];
4721 newFaceEdges[1][nE[1]++] = newEdgeIndex[2];
4722 newFaceEdges[2][nE[2]++] = newEdgeIndex[3];
4724 newFaceEdges[3][nE[3]++] = newEdgeIndex[1];
4725 newFaceEdges[3][nE[3]++] = newEdgeIndex[3];
4726 newFaceEdges[4][nE[4]++] = newEdgeIndex[1];
4727 newFaceEdges[4][nE[4]++] = newEdgeIndex[2];
4728 newFaceEdges[5][nE[5]++] = newEdgeIndex[2];
4729 newFaceEdges[5][nE[5]++] = newEdgeIndex[3];
4731 // Define the six edges to check while building faceEdges:
4732 FixedList<edge,6> check;
4734 check[0][0] = apexPoint[0]; check[0][1] = faces_[fIndex][0];
4735 check[1][0] = apexPoint[0]; check[1][1] = faces_[fIndex][1];
4736 check[2][0] = apexPoint[0]; check[2][1] = faces_[fIndex][2];
4738 check[3][0] = faces_[fIndex][2]; check[3][1] = faces_[fIndex][0];
4739 check[4][0] = faces_[fIndex][0]; check[4][1] = faces_[fIndex][1];
4740 check[5][0] = faces_[fIndex][1]; check[5][1] = faces_[fIndex][2];
4742 // Build a list of cellEdges
4743 DynamicList<label> cellEdges(6);
4745 forAll(cells_[owner_[fIndex]], faceI)
4746 {
4747 const labelList& fEdges =
4748 (
4749 faceEdges_[cells_[owner_[fIndex]][faceI]]
4750 );
4752 forAll(fEdges, edgeI)
4753 {
4754 if (findIndex(cellEdges, fEdges[edgeI]) == -1)
4755 {
4756 cellEdges.append(fEdges[edgeI]);
4757 }
4758 }
4759 }
4761 // Loop through cellEdges, and perform appropriate actions.
4762 forAll(cellEdges, edgeI)
4763 {
4764 const label ceIndex = cellEdges[edgeI];
4765 const edge& edgeToCheck = edges_[ceIndex];
4767 // Check against the specified edges.
4768 if (edgeToCheck == check[0])
4769 {
4770 meshOps::sizeUpList(newFaceIndex[0], edgeFaces_[ceIndex]);
4771 newFaceEdges[0][nE[0]++] = ceIndex;
4772 }
4774 if (edgeToCheck == check[1])
4775 {
4776 meshOps::sizeUpList(newFaceIndex[1], edgeFaces_[ceIndex]);
4777 newFaceEdges[1][nE[1]++] = ceIndex;
4778 }
4780 if (edgeToCheck == check[2])
4781 {
4782 meshOps::sizeUpList(newFaceIndex[2], edgeFaces_[ceIndex]);
4783 newFaceEdges[2][nE[2]++] = ceIndex;
4784 }
4786 if (edgeToCheck == check[3])
4787 {
4788 meshOps::replaceLabel
4789 (
4790 fIndex,
4791 newFaceIndex[3],
4792 edgeFaces_[ceIndex]
4793 );
4795 newFaceEdges[3][nE[3]++] = ceIndex;
4796 }
4798 if (edgeToCheck == check[4])
4799 {
4800 meshOps::replaceLabel
4801 (
4802 fIndex,
4803 newFaceIndex[4],
4804 edgeFaces_[ceIndex]
4805 );
4807 newFaceEdges[4][nE[4]++] = ceIndex;
4808 }
4810 if (edgeToCheck == check[5])
4811 {
4812 meshOps::replaceLabel
4813 (
4814 fIndex,
4815 newFaceIndex[5],
4816 edgeFaces_[ceIndex]
4817 );
4819 newFaceEdges[5][nE[5]++] = ceIndex;
4820 }
4821 }
4823 // Now that faceEdges has been configured, append them to the list.
4824 for (label i = 0; i < 6; i++)
4825 {
4826 faceEdges_.append(newFaceEdges[i]);
4828 // Add faces to the map.
4829 map.addFace(newFaceIndex[i]);
4830 }
4832 // If modification is coupled, generate mapping info.
4833 if (coupledModification_)
4834 {
4835 // Create a master/slave entry for the new edges/faces on the patch.
4836 if (locallyCoupledEntity(fIndex, true))
4837 {
4838 if (patchCoupling_(pIndex))
4839 {
4840 // Add the new point to the coupling map
4841 const coupleMap& cMap = patchCoupling_[pIndex].map();
4843 if (trisectingSlave)
4844 {
4845 cMap.mapMaster
4846 (
4847 coupleMap::POINT,
4848 newPointIndex,
4849 slaveMap.addedPointList()[0].index()
4850 );
4852 cMap.mapSlave
4853 (
4854 coupleMap::POINT,
4855 slaveMap.addedPointList()[0].index(),
4856 newPointIndex
4857 );
4858 }
4859 else
4860 {
4861 cMap.mapSlave
4862 (
4863 coupleMap::POINT,
4864 newPointIndex,
4865 slaveMap.addedPointList()[0].index()
4866 );
4868 cMap.mapMaster
4869 (
4870 coupleMap::POINT,
4871 slaveMap.addedPointList()[0].index(),
4872 newPointIndex
4873 );
4874 }
4875 }
4877 const List<objectMap>& ameList = map.addedEdgeList();
4878 const List<objectMap>& aseList = slaveMap.addedEdgeList();
4880 // Compare with all check entries.
4881 forAll(ameList, meI)
4882 {
4883 label epIndex = whichEdgePatch(ameList[meI].index());
4885 if (epIndex != pIndex && !trisectingSlave)
4886 {
4887 continue;
4888 }
4890 if (patchCoupling_(pIndex))
4891 {
4892 const coupleMap& cM = patchCoupling_[pIndex].map();
4894 if (trisectingSlave && epIndex != cM.slaveIndex())
4895 {
4896 continue;
4897 }
4898 }
4899 else
4900 {
4901 continue;
4902 }
4904 const coupleMap& cMap = patchCoupling_[pIndex].map();
4906 // Configure an edge for comparison.
4907 edge cE(-1, -1);
4909 const edge& mE = edges_[ameList[meI].index()];
4911 if (trisectingSlave)
4912 {
4913 cE[0] = cMap.reverseEntityMap(coupleMap::POINT)[mE[0]];
4914 cE[1] = cMap.reverseEntityMap(coupleMap::POINT)[mE[1]];
4915 }
4916 else
4917 {
4918 cE[0] = cMap.entityMap(coupleMap::POINT)[mE[0]];
4919 cE[1] = cMap.entityMap(coupleMap::POINT)[mE[1]];
4920 }
4922 bool matched = false;
4924 forAll(aseList, seI)
4925 {
4926 const edge& sE = edges_[aseList[seI].index()];
4928 if (cE == sE)
4929 {
4930 if (trisectingSlave)
4931 {
4932 cMap.mapMaster
4933 (
4934 coupleMap::EDGE,
4935 ameList[meI].index(),
4936 aseList[seI].index()
4937 );
4939 cMap.mapSlave
4940 (
4941 coupleMap::EDGE,
4942 aseList[seI].index(),
4943 ameList[meI].index()
4944 );
4945 }
4946 else
4947 {
4948 cMap.mapSlave
4949 (
4950 coupleMap::EDGE,
4951 ameList[meI].index(),
4952 aseList[seI].index()
4953 );
4955 cMap.mapMaster
4956 (
4957 coupleMap::EDGE,
4958 aseList[seI].index(),
4959 ameList[meI].index()
4960 );
4961 }
4963 matched = true;
4965 break;
4966 }
4967 }
4969 if (!matched)
4970 {
4971 Pout<< "masterEdges: " << nl
4972 << ameList << endl;
4974 Pout<< "slaveEdges: " << nl
4975 << aseList << endl;
4977 forAll(ameList, meI)
4978 {
4979 Pout<< ameList[meI].index() << ": "
4980 << edges_[ameList[meI].index()]
4981 << endl;
4982 }
4984 forAll(aseList, seI)
4985 {
4986 Pout<< aseList[seI].index() << ": "
4987 << edges_[aseList[seI].index()]
4988 << endl;
4989 }
4991 FatalErrorIn
4992 (
4993 "const changeMap dynamicTopoFvMesh::trisectFace\n"
4994 "(\n"
4995 " const label fIndex,\n"
4996 " bool checkOnly,\n"
4997 " bool forceOp\n"
4998 ")\n"
4999 )
5000 << "Failed to build coupled edge maps."
5001 << abort(FatalError);
5002 }
5003 }
5005 // Add a mapping entry for three new faces as well.
5006 face cF(3);
5008 const List<objectMap>& amfList = map.addedFaceList();
5009 const List<objectMap>& asfList = slaveMap.addedFaceList();
5011 forAll(amfList, mfI)
5012 {
5013 label fpIndex = whichPatch(amfList[mfI].index());
5015 if (fpIndex != pIndex && !trisectingSlave)
5016 {
5017 continue;
5018 }
5020 if (patchCoupling_(pIndex))
5021 {
5022 const coupleMap& cM = patchCoupling_[pIndex].map();
5024 if (trisectingSlave && fpIndex != cM.slaveIndex())
5025 {
5026 continue;
5027 }
5028 }
5029 else
5030 {
5031 continue;
5032 }
5034 const coupleMap& cMap = patchCoupling_[pIndex].map();
5036 // Configure a face for comparison.
5037 const face& mF = faces_[amfList[mfI].index()];
5039 forAll(mF, pI)
5040 {
5041 if (trisectingSlave)
5042 {
5043 cF[pI] =
5044 (
5045 cMap.reverseEntityMap(coupleMap::POINT)[mF[pI]]
5046 );
5047 }
5048 else
5049 {
5050 cF[pI] =
5051 (
5052 cMap.entityMap(coupleMap::POINT)[mF[pI]]
5053 );
5054 }
5055 }
5057 bool matched = false;
5059 forAll(asfList, sfI)
5060 {
5061 const face& sF = faces_[asfList[sfI].index()];
5063 if (triFace::compare(triFace(cF), triFace(sF)))
5064 {
5065 if (trisectingSlave)
5066 {
5067 cMap.mapMaster
5068 (
5069 coupleMap::FACE,
5070 amfList[mfI].index(),
5071 asfList[sfI].index()
5072 );
5074 cMap.mapSlave
5075 (
5076 coupleMap::FACE,
5077 asfList[sfI].index(),
5078 amfList[mfI].index()
5079 );
5080 }
5081 else
5082 {
5083 cMap.mapSlave
5084 (
5085 coupleMap::FACE,
5086 amfList[mfI].index(),
5087 asfList[sfI].index()
5088 );
5090 cMap.mapMaster
5091 (
5092 coupleMap::FACE,
5093 asfList[sfI].index(),
5094 amfList[mfI].index()
5095 );
5096 }
5098 matched = true;
5100 break;
5101 }
5102 }
5104 if (!matched)
5105 {
5106 Pout<< "masterFaces: " << nl
5107 << amfList << endl;
5109 Pout<< "slaveFaces: " << nl
5110 << asfList << endl;
5112 forAll(amfList, mfI)
5113 {
5114 Pout<< amfList[mfI].index() << ": "
5115 << faces_[amfList[mfI].index()]
5116 << endl;
5117 }
5119 forAll(asfList, sfI)
5120 {
5121 Pout<< asfList[sfI].index() << ": "
5122 << faces_[asfList[sfI].index()]
5123 << endl;
5124 }
5126 FatalErrorIn
5127 (
5128 "const changeMap dynamicTopoFvMesh::trisectFace\n"
5129 "(\n"
5130 " const label fIndex,\n"
5131 " bool checkOnly,\n"
5132 " bool forceOp\n"
5133 ")\n"
5134 )
5135 << "Failed to build coupled face maps."
5136 << abort(FatalError);
5137 }
5138 }
5139 }
5140 else
5141 if (processorCoupledEntity(fIndex))
5142 {
5143 // Look for matching slave edges on the patchSubMesh.
5145 }
5146 }
5147 }
5148 else
5149 {
5150 // Add three new cells to the end of the cell list
5151 for (label i = 3; i < 6; i++)
5152 {
5153 scalar parentScale = -1.0;
5155 if (edgeRefinement_)
5156 {
5157 parentScale = lengthScale_[cellsForRemoval[1]];
5158 }
5160 newCellIndex[i] = insertCell(newTetCell[i], parentScale);
5162 // Add to the map.
5163 map.addCell(newCellIndex[i]);
5164 }
5166 // Find the apex point for this cell
5167 apexPoint[1] =
5168 (
5169 meshOps::tetApexPoint
5170 (
5171 neighbour_[fIndex],
5172 fIndex,
5173 faces_,
5174 cells_
5175 )
5176 );
5178 // Add six new interior faces.
5180 // Fourth face: Owner: newCellIndex[0], Neighbour: newCellIndex[3]
5181 tmpTriFace[0] = newPointIndex;
5182 tmpTriFace[1] = faces_[fIndex][2];
5183 tmpTriFace[2] = faces_[fIndex][0];
5185 newFaceIndex[3] =
5186 (
5187 insertFace
5188 (
5189 -1,
5190 tmpTriFace,
5191 newCellIndex[0],
5192 newCellIndex[3]
5193 )
5194 );
5196 // Fifth face: Owner: newCellIndex[1], Neighbour: newCellIndex[4]
5197 tmpTriFace[0] = newPointIndex;
5198 tmpTriFace[1] = faces_[fIndex][0];
5199 tmpTriFace[2] = faces_[fIndex][1];
5201 newFaceIndex[4] =
5202 (
5203 insertFace
5204 (
5205 -1,
5206 tmpTriFace,
5207 newCellIndex[1],
5208 newCellIndex[4]
5209 )
5210 );
5212 // Sixth face: Owner: newCellIndex[2], Neighbour: newCellIndex[5]
5213 tmpTriFace[0] = newPointIndex;
5214 tmpTriFace[1] = faces_[fIndex][1];
5215 tmpTriFace[2] = faces_[fIndex][2];
5217 newFaceIndex[5] =
5218 (
5219 insertFace
5220 (
5221 -1,
5222 tmpTriFace,
5223 newCellIndex[2],
5224 newCellIndex[5]
5225 )
5226 );
5228 // Seventh face: Owner: newCellIndex[3], Neighbour: newCellIndex[4]
5229 tmpTriFace[0] = newPointIndex;
5230 tmpTriFace[1] = apexPoint[1];
5231 tmpTriFace[2] = faces_[fIndex][0];
5233 newFaceIndex[6] =
5234 (
5235 insertFace
5236 (
5237 -1,
5238 tmpTriFace,
5239 newCellIndex[3],
5240 newCellIndex[4]
5241 )
5242 );
5244 // Eighth face: Owner: newCellIndex[4], Neighbour: newCellIndex[5]
5245 tmpTriFace[0] = newPointIndex;
5246 tmpTriFace[1] = apexPoint[1];
5247 tmpTriFace[2] = faces_[fIndex][1];
5249 newFaceIndex[7] =
5250 (
5251 insertFace
5252 (
5253 -1,
5254 tmpTriFace,
5255 newCellIndex[4],
5256 newCellIndex[5]
5257 )
5258 );
5260 // Ninth face: Owner: newCellIndex[3], Neighbour: newCellIndex[5]
5261 tmpTriFace[0] = newPointIndex;
5262 tmpTriFace[1] = faces_[fIndex][2];
5263 tmpTriFace[2] = apexPoint[1];
5265 newFaceIndex[8] =
5266 (
5267 insertFace
5268 (
5269 -1,
5270 tmpTriFace,
5271 newCellIndex[3],
5272 newCellIndex[5]
5273 )
5274 );
5276 // Add the newly created faces to cells
5277 newTetCell[3][nF[3]++] = newFaceIndex[6];
5278 newTetCell[3][nF[3]++] = newFaceIndex[8];
5279 newTetCell[4][nF[4]++] = newFaceIndex[6];
5280 newTetCell[4][nF[4]++] = newFaceIndex[7];
5281 newTetCell[5][nF[5]++] = newFaceIndex[7];
5282 newTetCell[5][nF[5]++] = newFaceIndex[8];
5284 newTetCell[0][nF[0]++] = newFaceIndex[3];
5285 newTetCell[1][nF[1]++] = newFaceIndex[4];
5286 newTetCell[2][nF[2]++] = newFaceIndex[5];
5288 newTetCell[3][nF[3]++] = newFaceIndex[3];
5289 newTetCell[4][nF[4]++] = newFaceIndex[4];
5290 newTetCell[5][nF[5]++] = newFaceIndex[5];
5292 // Define the three faces to check for orientation:
5293 checkFace[0][0] = faces_[fIndex][2];
5294 checkFace[0][1] = apexPoint[1];
5295 checkFace[0][2] = faces_[fIndex][0];
5297 checkFace[1][0] = faces_[fIndex][0];
5298 checkFace[1][1] = apexPoint[1];
5299 checkFace[1][2] = faces_[fIndex][1];
5301 checkFace[2][0] = faces_[fIndex][1];
5302 checkFace[2][1] = apexPoint[1];
5303 checkFace[2][2] = faces_[fIndex][2];
5305 // Check the orientation of faces on the second cell.
5306 forAll(cells_[neighbour_[fIndex]], faceI)
5307 {
5308 label faceIndex = cells_[neighbour_[fIndex]][faceI];
5310 if (faceIndex == fIndex)
5311 {
5312 continue;
5313 }
5315 const face& faceToCheck = faces_[faceIndex];
5316 label cellIndex = cellsForRemoval[1];
5317 label newIndex = -1;
5319 // Check against faces.
5320 if (triFace::compare(triFace(faceToCheck), triFace(checkFace[0])))
5321 {
5322 newIndex = newCellIndex[3];
5323 newTetCell[3][nF[3]++] = faceIndex;
5324 }
5325 else
5326 if (triFace::compare(triFace(faceToCheck), triFace(checkFace[1])))
5327 {
5328 newIndex = newCellIndex[4];
5329 newTetCell[4][nF[4]++] = faceIndex;
5330 }
5331 else
5332 if (triFace::compare(triFace(faceToCheck), triFace(checkFace[2])))
5333 {
5334 newIndex = newCellIndex[5];
5335 newTetCell[5][nF[5]++] = faceIndex;
5336 }
5337 else
5338 {
5339 // Something's terribly wrong.
5340 FatalErrorIn
5341 (
5342 "const changeMap dynamicTopoFvMesh::trisectFace\n"
5343 "(\n"
5344 " const label fIndex,\n"
5345 " bool checkOnly,\n"
5346 " bool forceOp\n"
5347 ")\n"
5348 )
5349 << "Failed to determine a face match."
5350 << abort(FatalError);
5351 }
5353 // Check if a face-flip is necessary
5354 if (owner_[faceIndex] == cellIndex)
5355 {
5356 if (neighbour_[faceIndex] == -1)
5357 {
5358 // Change the owner
5359 owner_[faceIndex] = newIndex;
5360 }
5361 else
5362 {
5363 // Flip this face
5364 faces_[faceIndex] = faceToCheck.reverseFace();
5365 owner_[faceIndex] = neighbour_[faceIndex];
5366 neighbour_[faceIndex] = newIndex;
5368 setFlip(faceIndex);
5369 }
5370 }
5371 else
5372 {
5373 // Flip is unnecessary. Just update neighbour
5374 neighbour_[faceIndex] = newIndex;
5375 }
5376 }
5378 // Configure edgeFaces for four new interior edges.
5379 newQuadEdgeFaces[0] = newFaceIndex[4];
5380 newQuadEdgeFaces[1] = newFaceIndex[0];
5381 newQuadEdgeFaces[2] = newFaceIndex[3];
5382 newQuadEdgeFaces[3] = newFaceIndex[6];
5384 newEdgeIndex[1] =
5385 (
5386 insertEdge
5387 (
5388 -1,
5389 edge
5390 (
5391 newPointIndex,
5392 faces_[fIndex][0]
5393 ),
5394 newQuadEdgeFaces
5395 )
5396 );
5398 newQuadEdgeFaces[0] = newFaceIndex[5];
5399 newQuadEdgeFaces[1] = newFaceIndex[1];
5400 newQuadEdgeFaces[2] = newFaceIndex[4];
5401 newQuadEdgeFaces[3] = newFaceIndex[7];
5403 newEdgeIndex[2] =
5404 (
5405 insertEdge
5406 (
5407 -1,
5408 edge
5409 (
5410 newPointIndex,
5411 faces_[fIndex][1]
5412 ),
5413 newQuadEdgeFaces
5414 )
5415 );
5417 newQuadEdgeFaces[0] = newFaceIndex[3];
5418 newQuadEdgeFaces[1] = newFaceIndex[2];
5419 newQuadEdgeFaces[2] = newFaceIndex[5];
5420 newQuadEdgeFaces[3] = newFaceIndex[8];
5422 newEdgeIndex[3] =
5423 (
5424 insertEdge
5425 (
5426 -1,
5427 edge
5428 (
5429 newPointIndex,
5430 faces_[fIndex][2]
5431 ),
5432 newQuadEdgeFaces
5433 )
5434 );
5436 newTriEdgeFaces[0] = newFaceIndex[6];
5437 newTriEdgeFaces[1] = newFaceIndex[7];
5438 newTriEdgeFaces[2] = newFaceIndex[8];
5440 newEdgeIndex[4] =
5441 (
5442 insertEdge
5443 (
5444 -1,
5445 edge
5446 (
5447 apexPoint[1],
5448 newPointIndex
5449 ),
5450 newTriEdgeFaces
5451 )
5452 );
5454 // Configure faceEdges with the new internal edges
5455 newFaceEdges[0][nE[0]++] = newEdgeIndex[1];
5456 newFaceEdges[1][nE[1]++] = newEdgeIndex[2];
5457 newFaceEdges[2][nE[2]++] = newEdgeIndex[3];
5459 newFaceEdges[3][nE[3]++] = newEdgeIndex[1];
5460 newFaceEdges[3][nE[3]++] = newEdgeIndex[3];
5461 newFaceEdges[4][nE[4]++] = newEdgeIndex[1];
5462 newFaceEdges[4][nE[4]++] = newEdgeIndex[2];
5463 newFaceEdges[5][nE[5]++] = newEdgeIndex[2];
5464 newFaceEdges[5][nE[5]++] = newEdgeIndex[3];
5466 newFaceEdges[6][nE[6]++] = newEdgeIndex[1];
5467 newFaceEdges[7][nE[7]++] = newEdgeIndex[2];
5468 newFaceEdges[8][nE[8]++] = newEdgeIndex[3];
5470 newFaceEdges[6][nE[6]++] = newEdgeIndex[4];
5471 newFaceEdges[7][nE[7]++] = newEdgeIndex[4];
5472 newFaceEdges[8][nE[8]++] = newEdgeIndex[4];
5474 // Define the nine edges to check while building faceEdges:
5475 FixedList<edge,9> check;
5477 check[0][0] = apexPoint[0]; check[0][1] = faces_[fIndex][0];
5478 check[1][0] = apexPoint[0]; check[1][1] = faces_[fIndex][1];
5479 check[2][0] = apexPoint[0]; check[2][1] = faces_[fIndex][2];
5481 check[3][0] = faces_[fIndex][2]; check[3][1] = faces_[fIndex][0];
5482 check[4][0] = faces_[fIndex][0]; check[4][1] = faces_[fIndex][1];
5483 check[5][0] = faces_[fIndex][1]; check[5][1] = faces_[fIndex][2];
5485 check[6][0] = apexPoint[1]; check[6][1] = faces_[fIndex][0];
5486 check[7][0] = apexPoint[1]; check[7][1] = faces_[fIndex][1];
5487 check[8][0] = apexPoint[1]; check[8][1] = faces_[fIndex][2];
5489 // Build a list of cellEdges
5490 DynamicList<label> cellEdges(12);
5492 forAll(cells_[owner_[fIndex]], faceI)
5493 {
5494 const labelList& fEdges =
5495 (
5496 faceEdges_[cells_[owner_[fIndex]][faceI]]
5497 );
5499 forAll(fEdges, edgeI)
5500 {
5501 if (findIndex(cellEdges, fEdges[edgeI]) == -1)
5502 {
5503 cellEdges.append(fEdges[edgeI]);
5504 }
5505 }
5506 }
5508 forAll(cells_[neighbour_[fIndex]], faceI)
5509 {
5510 const labelList& fEdges =
5511 (
5512 faceEdges_[cells_[neighbour_[fIndex]][faceI]]
5513 );
5515 forAll(fEdges, edgeI)
5516 {
5517 if (findIndex(cellEdges, fEdges[edgeI]) == -1)
5518 {
5519 cellEdges.append(fEdges[edgeI]);
5520 }
5521 }
5522 }
5524 // Loop through cellEdges, and perform appropriate actions.
5525 forAll(cellEdges, edgeI)
5526 {
5527 const label ceIndex = cellEdges[edgeI];
5528 const edge& edgeToCheck = edges_[ceIndex];
5530 // Check against the specified edges.
5531 if (edgeToCheck == check[0])
5532 {
5533 meshOps::sizeUpList(newFaceIndex[0], edgeFaces_[ceIndex]);
5534 newFaceEdges[0][nE[0]++] = ceIndex;
5535 }
5537 if (edgeToCheck == check[1])
5538 {
5539 meshOps::sizeUpList(newFaceIndex[1], edgeFaces_[ceIndex]);
5540 newFaceEdges[1][nE[1]++] = ceIndex;
5541 }
5543 if (edgeToCheck == check[2])
5544 {
5545 meshOps::sizeUpList(newFaceIndex[2], edgeFaces_[ceIndex]);
5546 newFaceEdges[2][nE[2]++] = ceIndex;
5547 }
5549 if (edgeToCheck == check[3])
5550 {
5551 meshOps::replaceLabel
5552 (
5553 fIndex,
5554 newFaceIndex[3],
5555 edgeFaces_[ceIndex]
5556 );
5558 newFaceEdges[3][nE[3]++] = ceIndex;
5559 }
5561 if (edgeToCheck == check[4])
5562 {
5563 meshOps::replaceLabel
5564 (
5565 fIndex,
5566 newFaceIndex[4],
5567 edgeFaces_[ceIndex]
5568 );
5570 newFaceEdges[4][nE[4]++] = ceIndex;
5571 }
5573 if (edgeToCheck == check[5])
5574 {
5575 meshOps::replaceLabel
5576 (
5577 fIndex,
5578 newFaceIndex[5],
5579 edgeFaces_[ceIndex]
5580 );
5582 newFaceEdges[5][nE[5]++] = ceIndex;
5583 }
5585 if (edgeToCheck == check[6])
5586 {
5587 meshOps::sizeUpList(newFaceIndex[6], edgeFaces_[ceIndex]);
5588 newFaceEdges[6][nE[6]++] = ceIndex;
5589 }
5591 if (edgeToCheck == check[7])
5592 {
5593 meshOps::sizeUpList(newFaceIndex[7], edgeFaces_[ceIndex]);
5594 newFaceEdges[7][nE[7]++] = ceIndex;
5595 }
5597 if (edgeToCheck == check[8])
5598 {
5599 meshOps::sizeUpList(newFaceIndex[8], edgeFaces_[ceIndex]);
5600 newFaceEdges[8][nE[8]++] = ceIndex;
5601 }
5602 }
5604 // Now that faceEdges has been configured, append them to the list.
5605 for (label i = 0; i < 9; i++)
5606 {
5607 faceEdges_.append(newFaceEdges[i]);
5609 // Add faces to the map.
5610 map.addFace(newFaceIndex[i]);
5611 }
5612 }
5614 // Added edges are those connected to the new point
5615 const labelList& pointEdges = pointEdges_[newPointIndex];
5617 forAll(pointEdges, edgeI)
5618 {
5619 map.addEdge(pointEdges[edgeI]);
5620 }
5622 // Now generate mapping info and remove entities.
5623 forAll(cellsForRemoval, cellI)
5624 {
5625 label cIndex = cellsForRemoval[cellI];
5627 if (cIndex == -1)
5628 {
5629 continue;
5630 }
5632 // Fill-in mapping information
5633 labelList mC(1, cellsForRemoval[cellI]);
5635 if (cellI == 0)
5636 {
5637 for (label i = 0; i < 3; i++)
5638 {
5639 // Update the cell list with newly configured cells.
5640 cells_[newCellIndex[i]] = newTetCell[i];
5642 setCellMapping(newCellIndex[i], mC);
5643 }
5644 }
5645 else
5646 {
5647 for (label i = 3; i < 6; i++)
5648 {
5649 // Update the cell list with newly configured cells.
5650 cells_[newCellIndex[i]] = newTetCell[i];
5652 setCellMapping(newCellIndex[i], mC);
5653 }
5654 }
5656 removeCell(cIndex);
5657 }
5659 // Set default mapping for interior faces.
5660 for (label i = 0; i < 3; i++)
5661 {
5662 setFaceMapping(newFaceIndex[i]);
5663 }
5665 if (cellsForRemoval[1] == -1)
5666 {
5667 // Set mapping for boundary faces.
5668 for (label i = 3; i < 6; i++)
5669 {
5670 setFaceMapping(newFaceIndex[i], labelList(1, fIndex));
5671 }
5672 }
5673 else
5674 {
5675 // Set default mapping for interior faces.
5676 for (label i = 3; i < 9; i++)
5677 {
5678 setFaceMapping(newFaceIndex[i]);
5679 }
5680 }
5682 // Now finally remove the face...
5683 removeFace(fIndex);
5685 if (debug > 2)
5686 {
5687 Pout<< "New Point:: " << newPointIndex << endl;
5689 const labelList& pEdges = pointEdges_[newPointIndex];
5691 Pout<< "pointEdges:: " << pEdges << endl;
5693 Pout<< "Added edges: " << endl;
5695 forAll(pEdges, edgeI)
5696 {
5697 Pout<< pEdges[edgeI]
5698 << ":: " << edges_[pEdges[edgeI]] << nl
5699 << " edgeFaces:: " << edgeFaces_[pEdges[edgeI]] << nl
5700 << endl;
5701 }
5703 Pout<< "Added faces: " << endl;
5705 forAll(newFaceIndex, faceI)
5706 {
5707 if (newFaceIndex[faceI] == -1)
5708 {
5709 continue;
5710 }
5712 Pout<< newFaceIndex[faceI] << ":: "
5713 << faces_[newFaceIndex[faceI]]
5714 << endl;
5715 }
5717 Pout<< "Added cells: " << endl;
5719 forAll(newCellIndex, cellI)
5720 {
5721 if (newCellIndex[cellI] == -1)
5722 {
5723 continue;
5724 }
5726 Pout<< newCellIndex[cellI] << ":: "
5727 << cells_[newCellIndex[cellI]]
5728 << endl;
5729 }
5731 // Write out VTK files after change
5732 if (debug > 3)
5733 {
5734 labelList vtkCells;
5736 if (cellsForRemoval[1] == -1)
5737 {
5738 vtkCells.setSize(3);
5740 // Fill in cell indices
5741 vtkCells[0] = newCellIndex[0];
5742 vtkCells[1] = newCellIndex[1];
5743 vtkCells[2] = newCellIndex[2];
5744 }
5745 else
5746 {
5747 vtkCells.setSize(6);
5749 // Fill in cell indices
5750 forAll(newCellIndex, indexI)
5751 {
5752 vtkCells[indexI] = newCellIndex[indexI];
5753 }
5754 }
5756 writeVTK
5757 (
5758 Foam::name(fIndex)
5759 + "Trisect_1",
5760 vtkCells
5761 );
5762 }
5763 }
5765 // Set the flag
5766 topoChangeFlag_ = true;
5768 // Increment the counter
5769 statistics_[3]++;
5771 // Increment the number of modifications
5772 statistics_[0]++;
5774 // Specify that the operation was successful
5775 map.type() = 1;
5777 // Return the changeMap
5778 return map;
5779 }
5782 // Utility method to compute the quality of a
5783 // vertex hull around an edge after bisection.
5784 scalar dynamicTopoFvMesh::computeBisectionQuality
5785 (
5786 const label eIndex
5787 ) const
5788 {
5789 scalar minQuality = GREAT, minVolume = GREAT;
5790 scalar cQuality = 0.0, oldVolume = 0.0;
5792 // Obtain a reference to this edge
5793 const edge& edgeToCheck = edges_[eIndex];
5795 // Obtain point references
5796 const point& a = points_[edgeToCheck[0]];
5797 const point& c = points_[edgeToCheck[1]];
5799 const point& aOld = oldPoints_[edgeToCheck[0]];
5800 const point& cOld = oldPoints_[edgeToCheck[1]];
5802 // Compute the mid-point of the edge
5803 point midPoint = 0.5 * (a + c);
5804 point oldPoint = 0.5 * (aOld + cOld);
5806 DynamicList<label> eCells(10);
5808 const labelList& eFaces = edgeFaces_[eIndex];
5810 // Accumulate cells connected to this edge
5811 forAll(eFaces, faceI)
5812 {
5813 label own = owner_[eFaces[faceI]];
5814 label nei = neighbour_[eFaces[faceI]];
5816 if (findIndex(eCells, own) == -1)
5817 {
5818 eCells.append(own);
5819 }
5821 if (nei == -1)
5822 {
5823 continue;
5824 }
5826 if (findIndex(eCells, nei) == -1)
5827 {
5828 eCells.append(nei);
5829 }
5830 }
5832 // Loop through all cells and compute quality
5833 forAll(eCells, cellI)
5834 {
5835 label cellIndex = eCells[cellI];
5837 const cell& checkCell = cells_[cellIndex];
5839 // Find two faces that don't contain the edge
5840 forAll(checkCell, faceI)
5841 {
5842 const face& checkFace = faces_[checkCell[faceI]];
5844 if
5845 (
5846 (findIndex(checkFace, edgeToCheck[0]) == -1) &&
5847 (findIndex(checkFace, edgeToCheck[1]) == -1)
5848 )
5849 {
5850 // Check orientation
5851 if (owner_[checkCell[faceI]] == cellIndex)
5852 {
5853 cQuality =
5854 (
5855 tetMetric_
5856 (
5857 points_[checkFace[2]],
5858 points_[checkFace[1]],
5859 points_[checkFace[0]],
5860 midPoint
5861 )
5862 );
5864 oldVolume =
5865 (
5866 tetPointRef
5867 (
5868 oldPoints_[checkFace[2]],
5869 oldPoints_[checkFace[1]],
5870 oldPoints_[checkFace[0]],
5871 oldPoint
5872 ).mag()
5873 );
5874 }
5875 else
5876 {
5877 cQuality =
5878 (
5879 tetMetric_
5880 (
5881 points_[checkFace[0]],
5882 points_[checkFace[1]],
5883 points_[checkFace[2]],
5884 midPoint
5885 )
5886 );
5888 oldVolume =
5889 (
5890 tetPointRef
5891 (
5892 oldPoints_[checkFace[0]],
5893 oldPoints_[checkFace[1]],
5894 oldPoints_[checkFace[2]],
5895 oldPoint
5896 ).mag()
5897 );
5898 }
5900 // Check if the quality is worse
5901 minQuality = Foam::min(cQuality, minQuality);
5902 minVolume = Foam::min(oldVolume, minVolume);
5903 }
5904 }
5905 }
5907 // Ensure that the mesh is valid
5908 if (minQuality < sliverThreshold_)
5909 {
5910 if (debug > 3 && minQuality < 0.0)
5911 {
5912 writeVTK(Foam::name(eIndex) + "_iCells", eCells);
5913 }
5915 if (debug > 2)
5916 {
5917 InfoIn
5918 (
5919 "scalar dynamicTopoFvMesh::computeBisectionQuality"
5920 "(const label eIndex) const"
5921 )
5922 << " Bisecting edge will fall below the"
5923 << " sliver threshold of: " << sliverThreshold_ << nl
5924 << " Edge: " << eIndex << ": " << edgeToCheck << nl
5925 << " Minimum Quality: " << minQuality << nl
5926 << " Minimum Volume: " << minVolume << nl
5927 << " Mid point: " << midPoint << nl
5928 << " Old point: " << oldPoint << nl
5929 << endl;
5930 }
5931 }
5933 // If a negative old-volume was encountered,
5934 // return an invalid quality.
5935 if (minVolume < 0.0)
5936 {
5937 return minVolume;
5938 }
5940 return minQuality;
5941 }
5944 // Utility method to compute the quality of cells
5945 // around a face after trisection.
5946 scalar dynamicTopoFvMesh::computeTrisectionQuality
5947 (
5948 const label fIndex
5949 ) const
5950 {
5951 scalar minQuality = GREAT;
5952 scalar cQuality = 0.0;
5954 point midPoint;
5956 // Fetch the midPoint
5957 midPoint = faces_[fIndex].centre(points_);
5959 FixedList<label,2> apexPoint(-1);
5961 // Find the apex point
5962 apexPoint[0] =
5963 (
5964 meshOps::tetApexPoint
5965 (
5966 owner_[fIndex],
5967 fIndex,
5968 faces_,
5969 cells_
5970 )
5971 );
5973 const face& faceToCheck = faces_[fIndex];
5975 forAll(faceToCheck, pointI)
5976 {
5977 // Pick vertices off the list
5978 const point& b = points_[faceToCheck[pointI]];
5979 const point& c = points_[apexPoint[0]];
5980 const point& d = points_[faceToCheck[faceToCheck.fcIndex(pointI)]];
5982 // Compute the quality of the upper half.
5983 cQuality = tetMetric_(midPoint, b, c, d);
5985 // Check if the quality is worse
5986 minQuality = Foam::min(cQuality, minQuality);
5987 }
5989 if (whichPatch(fIndex) == -1)
5990 {
5991 apexPoint[1] =
5992 (
5993 meshOps::tetApexPoint
5994 (
5995 neighbour_[fIndex],
5996 fIndex,
5997 faces_,
5998 cells_
5999 )
6000 );
6002 forAll(faceToCheck, pointI)
6003 {
6004 // Pick vertices off the list
6005 const point& b = points_[faceToCheck[pointI]];
6006 const point& c = points_[apexPoint[1]];
6007 const point& d = points_[faceToCheck[faceToCheck.rcIndex(pointI)]];
6009 // Compute the quality of the upper half.
6010 cQuality = tetMetric_(midPoint, b, c, d);
6012 // Check if the quality is worse
6013 minQuality = Foam::min(cQuality, minQuality);
6014 }
6015 }
6017 return minQuality;
6018 }
6021 // Slice the mesh at a particular location
6022 void dynamicTopoFvMesh::sliceMesh
6023 (
6024 const labelPair& pointPair
6025 )
6026 {
6027 if (debug > 1)
6028 {
6029 Pout<< nl << nl
6030 << "Pair: " << pointPair
6031 << " is to be used for mesh slicing. " << endl;
6032 }
6034 label patchIndex = -1;
6035 scalar dx = 0.0;
6036 vector gCentre = vector::zero;
6037 FixedList<vector, 2> fC(vector::zero);
6039 if (twoDMesh_)
6040 {
6041 patchIndex = whichPatch(pointPair.first());
6043 // Fetch face centres
6044 fC[0] = faces_[pointPair.first()].centre(points_);
6045 fC[1] = faces_[pointPair.second()].centre(points_);
6046 }
6047 else
6048 {
6049 // Find the patch that the edge-vertex is connected to.
6050 const labelList& pEdges = pointEdges_[pointPair.first()];
6052 forAll(pEdges, edgeI)
6053 {
6054 if ((patchIndex = whichEdgePatch(pEdges[edgeI])) > -1)
6055 {
6056 break;
6057 }
6058 }
6060 fC[0] = points_[pointPair.first()];
6061 fC[1] = points_[pointPair.second()];
6062 }
6064 linePointRef lpr(fC[0], fC[1]);
6066 // Specify the centre.
6067 gCentre = lpr.centre();
6069 // Specify a search distance
6070 dx = lpr.mag();
6072 // Is this edge in the vicinity of a previous slice-point?
6073 if (lengthEstimator().checkOldSlices(gCentre))
6074 {
6075 if (debug > 1)
6076 {
6077 Pout<< nl << nl
6078 << "Pair: " << pointPair
6079 << " is too close to another slice point. "
6080 << endl;
6081 }
6083 // Too close to another slice-point. Bail out.
6084 return;
6085 }
6087 // Choose a box around the centre and scan all
6088 // surface entities that fall into this region.
6089 boundBox bBox
6090 (
6091 gCentre - vector(dx, dx, dx),
6092 gCentre + vector(dx, dx, dx)
6093 );
6095 vector p = vector::zero, N = vector::zero;
6096 Map<vector> checkPoints, surfFaces;
6097 Map<edge> checkEdges;
6099 if (twoDMesh_)
6100 {
6101 // Assign plane point / normal
6102 p = gCentre;
6104 vector gNorm = faces_[pointPair.first()].normal(points_);
6106 gNorm /= (mag(gNorm) + VSMALL);
6108 // Since this is 2D, assume XY-plane here.
6109 N = (gNorm ^ vector(0.0, 0.0, 1.0));
6110 }
6111 else
6112 {
6113 // Prepare surface points / edges for Dijkstra's algorithm
6114 for (label edgeI = nOldInternalEdges_; edgeI < nEdges_; edgeI++)
6115 {
6116 if (edgeFaces_[edgeI].empty())
6117 {
6118 continue;
6119 }
6121 if (whichEdgePatch(edgeI) == patchIndex)
6122 {
6123 const edge& surfaceEdge = edges_[edgeI];
6125 if
6126 (
6127 (bBox.contains(points_[surfaceEdge[0]])) &&
6128 (bBox.contains(points_[surfaceEdge[1]]))
6129 )
6130 {
6131 checkEdges.insert(edgeI, surfaceEdge);
6133 if (!checkPoints.found(surfaceEdge[0]))
6134 {
6135 checkPoints.insert
6136 (
6137 surfaceEdge[0],
6138 points_[surfaceEdge[0]]
6139 );
6140 }
6142 if (!checkPoints.found(surfaceEdge[1]))
6143 {
6144 checkPoints.insert
6145 (
6146 surfaceEdge[1],
6147 points_[surfaceEdge[1]]
6148 );
6149 }
6151 // Add surface faces as well.
6152 const labelList& eFaces = edgeFaces_[edgeI];
6154 forAll(eFaces, faceI)
6155 {
6156 if
6157 (
6158 (neighbour_[eFaces[faceI]] == -1) &&
6159 (!surfFaces.found(eFaces[faceI]))
6160 )
6161 {
6162 vector surfNorm =
6163 (
6164 faces_[eFaces[faceI]].normal(points_)
6165 );
6167 surfFaces.insert(eFaces[faceI], surfNorm);
6168 }
6169 }
6170 }
6171 }
6172 }
6174 if (debug > 1)
6175 {
6176 Pout<< nl << nl
6177 << " Point [0]: " << points_[pointPair.first()] << nl
6178 << " Point [1]: " << points_[pointPair.second()] << endl;
6180 if (debug > 3)
6181 {
6182 writeVTK("slicePoints", checkPoints.toc(), 0);
6183 writeVTK("sliceEdges", checkEdges.toc(), 1);
6184 }
6185 }
6187 // Find the shortest path using Dijkstra's algorithm.
6188 Map<label> shortestPath;
6190 bool foundPath =
6191 (
6192 meshOps::Dijkstra
6193 (
6194 checkPoints,
6195 checkEdges,
6196 pointPair.first(),
6197 pointPair.second(),
6198 shortestPath
6199 )
6200 );
6202 // First normalize all face-normals
6203 forAllIter(Map<vector>, surfFaces, sIter)
6204 {
6205 sIter() /= (mag(sIter()) + VSMALL);
6206 }
6208 if (foundPath)
6209 {
6210 // Next, take cross-products with every other
6211 // vector in the list, and accumulate.
6212 forAllIter(Map<vector>, surfFaces, sIterI)
6213 {
6214 forAllIter(Map<vector>, surfFaces, sIterJ)
6215 {
6216 if (sIterI.key() != sIterJ.key())
6217 {
6218 vector n = (sIterI() ^ sIterJ());
6220 n /= (mag(n) + VSMALL);
6222 // Reverse the vector if necessary
6223 if ((N & n) < 0.0)
6224 {
6225 n = -n;
6226 }
6228 N += n;
6229 }
6230 }
6231 }
6233 N /= (mag(N) + VSMALL);
6235 // Obtain point position
6236 p = gCentre;
6237 }
6238 else
6239 {
6240 // Probably a membrane-type configuration.
6241 labelHashSet checkCells;
6243 // Prepare a bounding cylinder with radius dx.
6244 forAllIter(Map<vector>, surfFaces, sIter)
6245 {
6246 const face& thisFace = faces_[sIter.key()];
6248 if (thisFace.which(pointPair.first()) > -1)
6249 {
6250 N += sIter();
6251 }
6252 }
6254 // Normalize and reverse.
6255 N /= -(mag(N) + VSMALL);
6257 vector a0 = points_[pointPair.first()];
6258 vector a1 = points_[pointPair.second()];
6259 scalar dist = mag(a1 - a0);
6261 forAll(cells_, cellI)
6262 {
6263 if (cells_[cellI].empty())
6264 {
6265 continue;
6266 }
6268 scalar v = 0.0;
6269 vector x = vector::zero;
6271 meshOps::cellCentreAndVolume
6272 (
6273 cellI,
6274 points_,
6275 faces_,
6276 cells_,
6277 owner_,
6278 x,
6279 v
6280 );
6282 vector rx = (x - a0);
6283 vector ra = (rx & N)*N;
6285 // Check if point falls off cylinder ends.
6286 if (mag(ra) > dist || mag(ra) < 0.0)
6287 {
6288 continue;
6289 }
6291 vector r = (rx - ra);
6293 // Check if the magnitude of 'r' is within radius.
6294 if (mag(r) < dx)
6295 {
6296 checkCells.insert(cellI);
6297 }
6298 }
6300 labelList cList = checkCells.toc();
6302 if (debug > 1)
6303 {
6304 Pout<< "Dijkstra's algorithm could not find a path." << endl;
6306 if (debug > 3)
6307 {
6308 writeVTK("checkCells", cList, 3);
6309 }
6310 }
6312 changeMap sliceMap =
6313 (
6314 removeCells
6315 (
6316 cList,
6317 patchIndex,
6318 "Slice_"
6319 + Foam::name(pointPair.first()) + '_'
6320 + Foam::name(pointPair.second())
6321 )
6322 );
6324 if (debug)
6325 {
6326 checkConnectivity(10);
6327 }
6329 // Accumulate a list of projection points
6330 labelHashSet pPoints;
6332 const List<objectMap>& afList = sliceMap.addedFaceList();
6334 forAll(afList, faceI)
6335 {
6336 const face& thisFace = faces_[afList[faceI].index()];
6338 forAll(thisFace, pointI)
6339 {
6340 if (!pPoints.found(thisFace[pointI]))
6341 {
6342 pPoints.insert(thisFace[pointI]);
6343 }
6344 }
6345 }
6347 // Now project points in a series of intermediate steps.
6349 // Add an entry to sliceBoxes.
6350 lengthEstimator().appendBox(bBox);
6352 return;
6353 }
6354 }
6356 if (debug > 1)
6357 {
6358 Pout<< nl << nl
6359 << " Plane point: " << p << nl
6360 << " Plane normal: " << N << endl;
6361 }
6363 // Mark cells and interior faces that fall
6364 // within the bounding box.
6365 labelHashSet checkCells, checkFaces, splitFaces;
6366 Map<bool> cellColors;
6368 forAll(faces_, faceI)
6369 {
6370 if (faces_[faceI].empty())
6371 {
6372 continue;
6373 }
6375 if (twoDMesh_ && faces_[faceI].size() == 3)
6376 {
6377 continue;
6378 }
6380 vector fCentre = faces_[faceI].centre(points_);
6382 FixedList<label, 2> cellsToCheck(-1);
6383 cellsToCheck[0] = owner_[faceI];
6384 cellsToCheck[1] = neighbour_[faceI];
6386 if (bBox.contains(fCentre) && cellsToCheck[1] != -1)
6387 {
6388 // Add this internal face to the list.
6389 checkFaces.insert(faceI);
6391 scalar volume = 0.0;
6392 vector centre = vector::zero;
6394 forAll(cellsToCheck, cellI)
6395 {
6396 if (!checkCells.found(cellsToCheck[cellI]))
6397 {
6398 meshOps::cellCentreAndVolume
6399 (
6400 cellsToCheck[cellI],
6401 points_,
6402 faces_,
6403 cells_,
6404 owner_,
6405 centre,
6406 volume
6407 );
6409 checkCells.insert(cellsToCheck[cellI]);
6411 if (((centre - p) & N) > 0.0)
6412 {
6413 cellColors.insert(cellsToCheck[cellI], true);
6414 }
6415 else
6416 {
6417 cellColors.insert(cellsToCheck[cellI], false);
6418 }
6419 }
6420 }
6421 }
6422 }
6424 // Prepare a list of internal faces for mesh splitting.
6425 forAllIter(labelHashSet, checkFaces, fIter)
6426 {
6427 if
6428 (
6429 cellColors[owner_[fIter.key()]]
6430 != cellColors[neighbour_[fIter.key()]]
6431 )
6432 {
6433 splitFaces.insert(fIter.key());
6434 }
6436 // Loop through all points (and associated pointEdges)
6437 // for this face, and check if connected cells are also
6438 // present in the checkCells/cellColors list
6439 if (twoDMesh_)
6440 {
6441 const labelList& fEdges = faceEdges_[fIter.key()];
6443 forAll(fEdges, edgeI)
6444 {
6445 const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
6447 forAll(eFaces, faceI)
6448 {
6449 label own = owner_[eFaces[faceI]];
6450 label nei = neighbour_[eFaces[faceI]];
6452 if (!checkCells.found(own))
6453 {
6454 scalar volume = 0.0;
6455 vector centre = vector::zero;
6457 meshOps::cellCentreAndVolume
6458 (
6459 own,
6460 points_,
6461 faces_,
6462 cells_,
6463 owner_,
6464 centre,
6465 volume
6466 );
6468 checkCells.insert(own);
6470 if (((centre - p) & N) > 0.0)
6471 {
6472 cellColors.insert(own, true);
6473 }
6474 else
6475 {
6476 cellColors.insert(own, false);
6477 }
6478 }
6480 if (!checkCells.found(nei) && nei != -1)
6481 {
6482 scalar volume = 0.0;
6483 vector centre = vector::zero;
6485 meshOps::cellCentreAndVolume
6486 (
6487 nei,
6488 points_,
6489 faces_,
6490 cells_,
6491 owner_,
6492 centre,
6493 volume
6494 );
6496 checkCells.insert(nei);
6498 if (((centre - p) & N) > 0.0)
6499 {
6500 cellColors.insert(nei, true);
6501 }
6502 else
6503 {
6504 cellColors.insert(nei, false);
6505 }
6506 }
6507 }
6508 }
6509 }
6510 else
6511 {
6512 const face& faceToCheck = faces_[fIter.key()];
6514 forAll(faceToCheck, pointI)
6515 {
6516 const labelList& pEdges = pointEdges_[faceToCheck[pointI]];
6518 forAll(pEdges, edgeI)
6519 {
6520 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
6522 forAll(eFaces, faceI)
6523 {
6524 label own = owner_[eFaces[faceI]];
6525 label nei = neighbour_[eFaces[faceI]];
6527 if (!checkCells.found(own))
6528 {
6529 scalar volume = 0.0;
6530 vector centre = vector::zero;
6532 meshOps::cellCentreAndVolume
6533 (
6534 own,
6535 points_,
6536 faces_,
6537 cells_,
6538 owner_,
6539 centre,
6540 volume
6541 );
6543 checkCells.insert(own);
6545 if (((centre - p) & N) > 0.0)
6546 {
6547 cellColors.insert(own, true);
6548 }
6549 else
6550 {
6551 cellColors.insert(own, false);
6552 }
6553 }
6555 if (!checkCells.found(nei) && nei != -1)
6556 {
6557 scalar volume = 0.0;
6558 vector centre = vector::zero;
6560 meshOps::cellCentreAndVolume
6561 (
6562 nei,
6563 points_,
6564 faces_,
6565 cells_,
6566 owner_,
6567 centre,
6568 volume
6569 );
6571 checkCells.insert(nei);
6573 if (((centre - p) & N) > 0.0)
6574 {
6575 cellColors.insert(nei, true);
6576 }
6577 else
6578 {
6579 cellColors.insert(nei, false);
6580 }
6581 }
6582 }
6583 }
6584 }
6585 }
6586 }
6588 if (debug > 3)
6589 {
6590 writeVTK("splitFaces", splitFaces.toc(), 2);
6591 writeVTK("checkCells", checkCells.toc(), 3);
6592 }
6594 // Pass this info into the splitInternalFaces routine.
6595 splitInternalFaces
6596 (
6597 patchIndex,
6598 splitFaces.toc(),
6599 cellColors
6600 );
6602 // Add an entry to sliceBoxes.
6603 lengthEstimator().appendBox(bBox);
6604 }
6607 // Add cell layer above specified patch
6608 const changeMap dynamicTopoFvMesh::addCellLayer
6609 (
6610 const label patchID
6611 )
6612 {
6613 changeMap map;
6615 // Maps for added entities
6616 Map<label> addedPoints;
6617 Map<label> addedHEdges, addedVEdges, currentVEdges;
6618 Map<label> addedHFaces, addedVFaces, currentVFaces;
6619 Map<labelPair> addedCells;
6621 // Allocate a list of patch faces
6622 DynamicList<label> patchFaces(patchSizes_[patchID]);
6624 // Loop through all patch faces and create a cell for each
6625 for (label faceI = nOldInternalFaces_; faceI < faces_.size(); faceI++)
6626 {
6627 label pIndex = whichPatch(faceI);
6629 if (pIndex != patchID)
6630 {
6631 continue;
6632 }
6634 // Add face to the list
6635 patchFaces.append(faceI);
6637 // Add a new cell
6638 label cIndex = owner_[faceI];
6639 scalar newLengthScale = -1.0;
6640 const cell& ownCell = cells_[cIndex];
6642 if (edgeRefinement_)
6643 {
6644 newLengthScale = lengthScale_[cIndex];
6645 }
6647 label newCellIndex =
6648 (
6649 insertCell
6650 (
6651 cell(ownCell.size()),
6652 newLengthScale
6653 )
6654 );
6656 // Update maps
6657 map.addCell(newCellIndex, labelList(1, cIndex));
6658 addedCells.insert(cIndex, labelPair(newCellIndex, 0));
6659 }
6661 FixedList<label, 2> mP(-1);
6663 forAll(patchFaces, indexI)
6664 {
6665 label faceI = patchFaces[indexI];
6666 label cIndex = owner_[faceI];
6668 // Fetch appropriate face / cell
6669 // - Make copies, since holding references
6670 // to data within this loop is unsafe.
6671 const face bFace = faces_[faceI];
6672 const cell ownCell = cells_[cIndex];
6674 // Configure a new face for insertion
6675 face newHFace(bFace);
6676 labelList newHFaceEdges(bFace.size(), -1);
6678 // Get the opposing face from the cell
6679 oppositeFace oFace = ownCell.opposingFace(faceI, faces_);
6681 if (!oFace.found())
6682 {
6683 // Something's wrong here.
6684 FatalErrorIn
6685 (
6686 "const changeMap dynamicTopoFvMesh::addCellLayer"
6687 "(const label patchID)"
6688 )
6689 << " Face: " << faceI << " :: " << bFace << nl
6690 << " has no opposing face in cell: "
6691 << cIndex << " :: " << ownCell << nl
6692 << abort(FatalError);
6693 }
6695 // Create points
6696 forAll(bFace, pointI)
6697 {
6698 label pIndex = bFace[pointI];
6700 // Skip if we've added this already
6701 if (addedPoints.found(pIndex))
6702 {
6703 continue;
6704 }
6706 // Set master points
6707 mP[0] = pIndex;
6708 mP[1] = oFace[pointI];
6710 label newPointIndex =
6711 (
6712 insertPoint
6713 (
6714 0.5 * (points_[mP[0]] + points_[mP[1]]),
6715 oldPoints_[mP[0]],
6716 mP
6717 )
6718 );
6720 // Update maps
6721 map.addPoint(newPointIndex, mP);
6722 addedPoints.insert(pIndex, newPointIndex);
6723 }
6725 // Fetch faceEdges from opposite faces.
6726 // - Make copies, since holding references is unsafe
6727 const labelList bfEdges = faceEdges_[faceI];
6728 const labelList ofEdges = faceEdges_[oFace.oppositeIndex()];
6730 // Create edges for each edge of the new horizontal face
6731 forAll(bfEdges, edgeI)
6732 {
6733 label beIndex = bfEdges[edgeI];
6735 // Skip if we've added this already
6736 if (addedHEdges.found(beIndex))
6737 {
6738 // Update face edges for the new horizontal face
6739 newHFaceEdges[edgeI] = addedHEdges[beIndex];
6741 continue;
6742 }
6744 // Configure the new edge
6745 label oeIndex = -1;
6746 const edge bEdge = edges_[beIndex];
6748 // Build an edge for comparison
6749 edge cEdge
6750 (
6751 oFace[bFace.which(bEdge[0])],
6752 oFace[bFace.which(bEdge[1])]
6753 );
6755 forAll(ofEdges, edgeJ)
6756 {
6757 const edge& ofEdge = edges_[ofEdges[edgeJ]];
6759 if (cEdge == ofEdge)
6760 {
6761 oeIndex = ofEdges[edgeJ];
6762 break;
6763 }
6764 }
6766 if (oeIndex < 0)
6767 {
6768 FatalErrorIn
6769 (
6770 "const changeMap dynamicTopoFvMesh::addCellLayer"
6771 "(const label patchID)"
6772 )
6773 << " Could not find comparison edge: " << cEdge
6774 << " for edge: " << bEdge
6775 << abort(FatalError);
6776 }
6778 // Fetch patch information
6779 label hEdgePatch = whichEdgePatch(oeIndex);
6781 // Set indices
6782 edge newHEdge
6783 (
6784 addedPoints[bEdge[0]],
6785 addedPoints[bEdge[1]]
6786 );
6788 // Insert a new edge with empty edgeFaces
6789 label newHEdgeIndex =
6790 (
6791 insertEdge
6792 (
6793 hEdgePatch,
6794 newHEdge,
6795 labelList(0)
6796 )
6797 );
6799 // Update maps
6800 map.addEdge(newHEdgeIndex);
6801 addedHEdges.insert(beIndex, newHEdgeIndex);
6803 // Update face edges for the new horizontal face
6804 newHFaceEdges[edgeI] = newHEdgeIndex;
6806 // Add a new vertical face for this edge
6807 label vFaceIndex = -1;
6809 // Find a vertical face that contains both edges
6810 const labelList& beFaces = edgeFaces_[beIndex];
6812 forAll(beFaces, faceJ)
6813 {
6814 const labelList& testEdges = faceEdges_[beFaces[faceJ]];
6816 if
6817 (
6818 (findIndex(testEdges, beIndex) > -1) &&
6819 (findIndex(testEdges, oeIndex) > -1)
6820 )
6821 {
6822 vFaceIndex = beFaces[faceJ];
6823 break;
6824 }
6825 }
6827 if (vFaceIndex < 0)
6828 {
6829 FatalErrorIn
6830 (
6831 "const changeMap dynamicTopoFvMesh::addCellLayer"
6832 "(const label patchID)"
6833 )
6834 << " Could not find an appropriate vertical face"
6835 << " containing edges: " << oeIndex
6836 << " and " << beIndex
6837 << abort(FatalError);
6838 }
6840 // Find two vertical edges on this face
6841 const labelList& vfEdges = faceEdges_[vFaceIndex];
6843 forAll(vfEdges, edgeJ)
6844 {
6845 const edge& vfEdge = edges_[vfEdges[edgeJ]];
6847 forAll(bEdge, i)
6848 {
6849 if (vfEdge == edge(bEdge[i], cEdge[i]))
6850 {
6851 // Skip if we've added this already
6852 if (addedVEdges.found(bEdge[i]))
6853 {
6854 continue;
6855 }
6857 label veIndex = vfEdges[edgeJ];
6859 // Fetch edge patch information
6860 label vEdgePatch = whichEdgePatch(veIndex);
6862 // Set indices
6863 edge newVEdge
6864 (
6865 bEdge[i],
6866 addedPoints[bEdge[i]]
6867 );
6869 // Insert a new edge with empty edgeFaces
6870 label newVEdgeIndex =
6871 (
6872 insertEdge
6873 (
6874 vEdgePatch,
6875 newVEdge,
6876 labelList(0)
6877 )
6878 );
6880 // Update maps
6881 map.addEdge(newVEdgeIndex);
6882 addedVEdges.insert(bEdge[i], newVEdgeIndex);
6884 // Note edge indices for later renumbering
6885 currentVEdges.insert(bEdge[i], veIndex);
6886 }
6887 }
6888 }
6890 // Configure the new vertical face
6891 face newVFace(faces_[vFaceIndex]);
6892 label newOwner = -1, newNeighbour = -1;
6894 label oldOwner = owner_[vFaceIndex];
6895 label oldNeighbour = neighbour_[vFaceIndex];
6897 // Fetch owner / neighbour
6898 newOwner = addedCells[oldOwner].first();
6900 if (oldNeighbour > -1)
6901 {
6902 newNeighbour = addedCells[oldNeighbour].first();
6903 }
6905 // Replace point indices on the new face
6906 forAll(bEdge, i)
6907 {
6908 meshOps::replaceLabel
6909 (
6910 cEdge[i],
6911 addedPoints[bEdge[i]],
6912 newVFace
6913 );
6914 }
6916 // Note face indices for later renumbering
6917 currentVFaces.insert(beIndex, vFaceIndex);
6919 // Check if reversal is necessary
6920 if ((newNeighbour < newOwner) && (newNeighbour > -1))
6921 {
6922 // Flip face
6923 newVFace = newVFace.reverseFace();
6925 // Swap addressing
6926 Foam::Swap(newOwner, newNeighbour);
6927 Foam::Swap(oldOwner, oldNeighbour);
6928 }
6930 // Configure faceEdges for the new vertical face
6931 labelList newVFaceEdges(4, -1);
6933 newVFaceEdges[0] = beIndex;
6934 newVFaceEdges[1] = newHEdgeIndex;
6935 newVFaceEdges[2] = addedVEdges[bEdge[0]];
6936 newVFaceEdges[3] = addedVEdges[bEdge[1]];
6938 // Add the new vertical face
6939 label newVFaceIndex =
6940 (
6941 insertFace
6942 (
6943 whichPatch(vFaceIndex),
6944 newVFace,
6945 newOwner,
6946 newNeighbour
6947 )
6948 );
6950 // Add a faceEdges entry
6951 faceEdges_.append(newVFaceEdges);
6953 // Update maps
6954 map.addFace(newVFaceIndex, labelList(1, vFaceIndex));
6955 addedVFaces.insert(beIndex, newVFaceIndex);
6957 // Update face count on the new cells
6958 cells_[newOwner][addedCells[oldOwner].second()++] =
6959 (
6960 newVFaceIndex
6961 );
6963 if (newNeighbour > -1)
6964 {
6965 cells_[newNeighbour][addedCells[oldNeighbour].second()++] =
6966 (
6967 newVFaceIndex
6968 );
6969 }
6971 // Size up edgeFaces for each edge
6972 forAll(newVFaceEdges, edgeJ)
6973 {
6974 label vfeIndex = newVFaceEdges[edgeJ];
6976 meshOps::sizeUpList
6977 (
6978 newVFaceIndex,
6979 edgeFaces_[vfeIndex]
6980 );
6981 }
6982 }
6984 // Add a new interior face, with identical orientation
6985 forAll(newHFace, pointI)
6986 {
6987 newHFace[pointI] = addedPoints[newHFace[pointI]];
6988 }
6990 // Add the new horizontal face
6991 label newHFaceIndex =
6992 (
6993 insertFace
6994 (
6995 -1,
6996 newHFace,
6997 cIndex,
6998 addedCells[cIndex].first()
6999 )
7000 );
7002 // Add a faceEdges entry
7003 faceEdges_.append(newHFaceEdges);
7005 // Update maps
7006 map.addFace(newHFaceIndex, labelList(1, faceI));
7007 addedHFaces.insert(faceI, newHFaceIndex);
7009 // Replace index on the old cell
7010 meshOps::replaceLabel
7011 (
7012 faceI,
7013 newHFaceIndex,
7014 cells_[cIndex]
7015 );
7017 // Update face count on the new cell
7018 label newCellIndex = addedCells[cIndex].first();
7020 // Modify owner for the existing boundary face
7021 owner_[faceI] = newCellIndex;
7023 cells_[newCellIndex][addedCells[cIndex].second()++] = faceI;
7024 cells_[newCellIndex][addedCells[cIndex].second()++] = newHFaceIndex;
7026 // Size up edgeFaces for each edge
7027 forAll(newHFaceEdges, edgeI)
7028 {
7029 label hfeIndex = newHFaceEdges[edgeI];
7031 meshOps::sizeUpList
7032 (
7033 newHFaceIndex,
7034 edgeFaces_[hfeIndex]
7035 );
7036 }
7037 }
7039 // Renumber vertical edges
7040 forAllConstIter(Map<label>, currentVEdges, eIter)
7041 {
7042 // Fetch reference to edge
7043 edge& curEdge = edges_[eIter()];
7045 if (curEdge[0] == eIter.key())
7046 {
7047 curEdge[0] = addedPoints[eIter.key()];
7048 }
7050 if (curEdge[1] == eIter.key())
7051 {
7052 curEdge[1] = addedPoints[eIter.key()];
7053 }
7055 // Size down pointEdges
7056 if (!twoDMesh_)
7057 {
7058 meshOps::sizeDownList
7059 (
7060 eIter(),
7061 pointEdges_[eIter.key()]
7062 );
7064 meshOps::sizeUpList
7065 (
7066 eIter(),
7067 pointEdges_[addedPoints[eIter.key()]]
7068 );
7069 }
7070 }
7072 // Renumber vertical faces
7073 forAllConstIter(Map<label>, currentVFaces, fIter)
7074 {
7075 // Fetch reference to existing edge
7076 const edge& bEdge = edges_[fIter.key()];
7078 // Replace point indices on vertical face
7079 forAll(bEdge, i)
7080 {
7081 meshOps::replaceLabel
7082 (
7083 bEdge[i],
7084 addedPoints[bEdge[i]],
7085 faces_[fIter()]
7086 );
7087 }
7089 // Replace edge on the existing vertical face
7090 meshOps::replaceLabel
7091 (
7092 fIter.key(),
7093 addedHEdges[fIter.key()],
7094 faceEdges_[fIter()]
7095 );
7097 // Remove old face on existing boundary edge
7098 meshOps::sizeDownList
7099 (
7100 fIter(),
7101 edgeFaces_[fIter.key()]
7102 );
7104 // Add old face to new horizontal edge
7105 meshOps::sizeUpList
7106 (
7107 fIter(),
7108 edgeFaces_[addedHEdges[fIter.key()]]
7109 );
7110 }
7112 // Now that all old / new cells possess correct connectivity,
7113 // compute mapping information.
7114 const List<objectMap>& afList = map.addedFaceList();
7116 forAll(afList, indexI)
7117 {
7118 label parent = afList[indexI].masterObjects()[0];
7120 if (whichPatch(afList[indexI].index()) == -1)
7121 {
7122 // Interior faces get default mapping
7123 if (whichPatch(parent) == -1)
7124 {
7125 setFaceMapping(parent);
7126 }
7128 setFaceMapping(afList[indexI].index());
7129 }
7130 }
7132 // Set the flag
7133 topoChangeFlag_ = true;
7135 // Specify that the operation was successful
7136 map.type() = 1;
7138 // Return the changeMap
7139 return map;
7140 }
7143 // Split a set of internal faces into boundary faces
7144 // - Add boundary faces and edges to the patch specified by 'patchIndex'
7145 // - Cell color should specify a binary value dictating either side
7146 // of the split face.
7147 void dynamicTopoFvMesh::splitInternalFaces
7148 (
7149 const label patchIndex,
7150 const labelList& internalFaces,
7151 const Map<bool>& cellColors
7152 )
7153 {
7154 Map<label> mirrorPointLabels;
7155 FixedList<Map<label>, 2> mirrorEdgeLabels, mirrorFaceLabels;
7157 // First loop through the list and accumulate a list of
7158 // points and edges that need to be duplicated.
7159 forAll(internalFaces, faceI)
7160 {
7161 const face& faceToCheck = faces_[internalFaces[faceI]];
7163 forAll(faceToCheck, pointI)
7164 {
7165 if (!mirrorPointLabels.found(faceToCheck[pointI]))
7166 {
7167 mirrorPointLabels.insert(faceToCheck[pointI], -1);
7168 }
7169 }
7171 const labelList& fEdges = faceEdges_[internalFaces[faceI]];
7173 forAll(fEdges, edgeI)
7174 {
7175 if (!mirrorEdgeLabels[0].found(fEdges[edgeI]))
7176 {
7177 mirrorEdgeLabels[0].insert(fEdges[edgeI], -1);
7178 }
7179 }
7180 }
7182 // Now for every point in the list, add a new one.
7183 // Add a mapping entry as well.
7184 forAllIter(Map<label>, mirrorPointLabels, pIter)
7185 {
7186 // Obtain a copy of the point before adding it,
7187 // since the reference might become invalid during list resizing.
7188 point newPoint = points_[pIter.key()];
7189 point oldPoint = oldPoints_[pIter.key()];
7191 pIter() = insertPoint(newPoint, oldPoint, labelList(1, pIter.key()));
7193 if (!twoDMesh_)
7194 {
7195 const labelList& pEdges = pointEdges_[pIter.key()];
7197 labelHashSet edgesToRemove;
7199 forAll(pEdges, edgeI)
7200 {
7201 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
7203 bool allTrue = true;
7205 forAll(eFaces, faceI)
7206 {
7207 label own = owner_[eFaces[faceI]];
7208 label nei = neighbour_[eFaces[faceI]];
7210 // Check if an owner/neighbour cell is false
7211 if (!cellColors[own])
7212 {
7213 allTrue = false;
7214 break;
7215 }
7217 if (nei != -1)
7218 {
7219 if (!cellColors[nei])
7220 {
7221 allTrue = false;
7222 break;
7223 }
7224 }
7225 }
7227 if (allTrue)
7228 {
7229 // Mark this edge label to be discarded later
7230 edgesToRemove.insert(pEdges[edgeI]);
7231 }
7232 }
7234 // It is dangerous to use the pointEdges references,
7235 // so call it using array-lookup instead.
7236 forAllIter(labelHashSet, edgesToRemove, hsIter)
7237 {
7238 // Add the edge to the mirror point list
7239 meshOps::sizeUpList
7240 (
7241 hsIter.key(),
7242 pointEdges_[pIter()]
7243 );
7245 // Remove the edge from the original point list
7246 meshOps::sizeDownList
7247 (
7248 hsIter.key(),
7249 pointEdges_[pIter.key()]
7250 );
7251 }
7252 }
7253 }
7255 if (debug > 3)
7256 {
7257 label i = 0;
7258 labelList mPoints(mirrorPointLabels.size());
7260 if (!twoDMesh_)
7261 {
7262 forAllIter(Map<label>, mirrorPointLabels, pIter)
7263 {
7264 writeVTK
7265 (
7266 "pEdges_o_" + Foam::name(pIter.key()) + '_',
7267 pointEdges_[pIter.key()],
7268 1
7269 );
7271 writeVTK
7272 (
7273 "pEdges_m_" + Foam::name(pIter()) + '_',
7274 pointEdges_[pIter()],
7275 1
7276 );
7278 mPoints[i++] = pIter();
7279 }
7281 writeVTK
7282 (
7283 "points_o_",
7284 mirrorPointLabels.toc(),
7285 0
7286 );
7288 writeVTK
7289 (
7290 "points_m_",
7291 mPoints,
7292 0
7293 );
7294 }
7295 }
7297 // For every internal face, add a new one.
7298 // - Stick to the rule:
7299 // [1] Every cell marked false keeps the existing entities.
7300 // [2] Every cell marked true gets new points/edges/faces.
7301 // - If faces are improperly oriented, reverse them.
7302 forAll(internalFaces, faceI)
7303 {
7304 FixedList<face, 2> newFace;
7305 FixedList<label, 2> newFaceIndex(-1);
7306 FixedList<label, 2> newOwner(-1);
7308 label oldOwn = owner_[internalFaces[faceI]];
7309 label oldNei = neighbour_[internalFaces[faceI]];
7311 if (cellColors[oldOwn] && !cellColors[oldNei])
7312 {
7313 // The owner gets a new boundary face.
7314 // Note that orientation is already correct.
7315 newFace[0] = faces_[internalFaces[faceI]];
7317 // The neighbour needs to have its face reversed
7318 // and moved to the boundary patch, thereby getting
7319 // deleted in the process.
7320 newFace[1] = newFace[0].reverseFace();
7322 newOwner[0] = oldOwn;
7323 newOwner[1] = oldNei;
7324 }
7325 else
7326 if (!cellColors[oldOwn] && cellColors[oldNei])
7327 {
7328 // The neighbour gets a new boundary face.
7329 // The face is oriented in the opposite sense, however.
7330 newFace[0] = faces_[internalFaces[faceI]].reverseFace();
7332 // The owner keeps the existing face and orientation.
7333 // But it also needs to be moved to the boundary.
7334 newFace[1] = faces_[internalFaces[faceI]];
7336 newOwner[0] = oldNei;
7337 newOwner[1] = oldOwn;
7338 }
7339 else
7340 {
7341 // Something's wrong here.
7342 FatalErrorIn
7343 (
7344 "dynamicTopoFvMesh::splitInternalFaces\n"
7345 "(\n"
7346 " const label patchIndex,\n"
7347 " const labelList& internalFaces,\n"
7348 " const Map<bool>& cellColors\n"
7349 ")\n"
7350 )
7351 << nl << " Face: "
7352 << internalFaces[faceI]
7353 << " has cells which are improperly marked: " << nl
7354 << oldOwn << ":: " << cellColors[oldOwn] << nl
7355 << oldNei << ":: " << cellColors[oldNei]
7356 << abort(FatalError);
7357 }
7359 // Renumber point labels for the first new face.
7360 forAll(newFace[0], pointI)
7361 {
7362 newFace[0][pointI] = mirrorPointLabels[newFace[0][pointI]];
7363 }
7365 // Insert the new boundary faces.
7366 forAll(newFace, indexI)
7367 {
7368 newFaceIndex[indexI] =
7369 (
7370 insertFace
7371 (
7372 patchIndex,
7373 newFace[indexI],
7374 newOwner[indexI],
7375 -1
7376 )
7377 );
7379 // Make an identical faceEdges entry.
7380 // This will be renumbered once new edges are added.
7381 labelList newFaceEdges(faceEdges_[internalFaces[faceI]]);
7383 faceEdges_.append(newFaceEdges);
7385 // Replace face labels on cells
7386 meshOps::replaceLabel
7387 (
7388 internalFaces[faceI],
7389 newFaceIndex[indexI],
7390 cells_[newOwner[indexI]]
7391 );
7393 // Make face mapping entries for posterity.
7394 mirrorFaceLabels[indexI].insert
7395 (
7396 internalFaces[faceI],
7397 newFaceIndex[indexI]
7398 );
7399 }
7400 }
7402 // For every edge in the list, add a new one.
7403 // We'll deal with correcting edgeFaces later.
7404 forAllIter(Map<label>, mirrorEdgeLabels[0], eIter)
7405 {
7406 // Obtain copies for the append method
7407 edge origEdge = edges_[eIter.key()];
7408 labelList newEdgeFaces(edgeFaces_[eIter.key()]);
7410 eIter() =
7411 (
7412 insertEdge
7413 (
7414 patchIndex,
7415 edge
7416 (
7417 mirrorPointLabels[origEdge[0]],
7418 mirrorPointLabels[origEdge[1]]
7419 ),
7420 newEdgeFaces
7421 )
7422 );
7424 // Is the original edge an internal one?
7425 // If it is, we need to move it to the boundary.
7426 if (whichEdgePatch(eIter.key()) == -1)
7427 {
7428 label rplEdgeIndex =
7429 (
7430 insertEdge
7431 (
7432 patchIndex,
7433 origEdge,
7434 newEdgeFaces
7435 )
7436 );
7438 // Map the new entry.
7439 mirrorEdgeLabels[1].insert(eIter.key(), rplEdgeIndex);
7440 }
7441 else
7442 {
7443 // This is already a boundary edge.
7444 // Make an identical map.
7445 mirrorEdgeLabels[1].insert(eIter.key(), eIter.key());
7446 }
7447 }
7449 // Correct edgeFaces for all new edges
7450 forAll(mirrorEdgeLabels, indexI)
7451 {
7452 forAllIter(Map<label>, mirrorEdgeLabels[indexI], eIter)
7453 {
7454 labelList& eFaces = edgeFaces_[eIter()];
7456 labelHashSet facesToRemove;
7458 forAll(eFaces, faceI)
7459 {
7460 bool remove = false;
7462 label own = owner_[eFaces[faceI]];
7463 label nei = neighbour_[eFaces[faceI]];
7465 if
7466 (
7467 (!cellColors[own] && !indexI) ||
7468 ( cellColors[own] && indexI)
7469 )
7470 {
7471 remove = true;
7472 }
7474 if (nei != -1)
7475 {
7476 if
7477 (
7478 (!cellColors[nei] && !indexI) ||
7479 ( cellColors[nei] && indexI)
7480 )
7481 {
7482 remove = true;
7483 }
7484 }
7486 if (mirrorFaceLabels[indexI].found(eFaces[faceI]))
7487 {
7488 // Perform a replacement instead of a removal.
7489 eFaces[faceI] = mirrorFaceLabels[indexI][eFaces[faceI]];
7491 remove = false;
7492 }
7494 if (remove)
7495 {
7496 facesToRemove.insert(eFaces[faceI]);
7497 }
7498 }
7500 // Now successively size down edgeFaces.
7501 // It is dangerous to use the eFaces reference,
7502 // so call it using array-lookup.
7503 forAllIter(labelHashSet, facesToRemove, hsIter)
7504 {
7505 meshOps::sizeDownList(hsIter.key(), edgeFaces_[eIter()]);
7506 }
7507 }
7508 }
7510 // Renumber faceEdges for all faces connected to new edges
7511 forAll(mirrorEdgeLabels, indexI)
7512 {
7513 forAllIter(Map<label>, mirrorEdgeLabels[indexI], eIter)
7514 {
7515 const labelList& eFaces = edgeFaces_[eIter()];
7517 forAll(eFaces, faceI)
7518 {
7519 labelList& fEdges = faceEdges_[eFaces[faceI]];
7521 forAll(fEdges, edgeI)
7522 {
7523 if (mirrorEdgeLabels[indexI].found(fEdges[edgeI]))
7524 {
7525 fEdges[edgeI] =
7526 (
7527 mirrorEdgeLabels[indexI][fEdges[edgeI]]
7528 );
7529 }
7530 }
7531 }
7532 }
7533 }
7535 if (twoDMesh_)
7536 {
7537 // Renumber edges and faces
7538 forAllIter(Map<label>, mirrorEdgeLabels[0], eIter)
7539 {
7540 const labelList& eFaces = edgeFaces_[eIter()];
7542 // Two levels of indirection to ensure
7543 // that all entities we renumbered.
7544 // A flip-side for the lack of a pointEdges list in 2D.
7545 forAll(eFaces, faceI)
7546 {
7547 const labelList& fEdges = faceEdges_[eFaces[faceI]];
7549 forAll(fEdges, edgeI)
7550 {
7551 // Renumber this edge.
7552 edge& edgeToCheck = edges_[fEdges[edgeI]];
7554 forAll(edgeToCheck, pointI)
7555 {
7556 if (mirrorPointLabels.found(edgeToCheck[pointI]))
7557 {
7558 edgeToCheck[pointI] =
7559 (
7560 mirrorPointLabels[edgeToCheck[pointI]]
7561 );
7562 }
7563 }
7565 // Also renumber faces connected to this edge.
7566 const labelList& efFaces = edgeFaces_[fEdges[edgeI]];
7568 forAll(efFaces, faceJ)
7569 {
7570 face& faceToCheck = faces_[efFaces[faceJ]];
7572 forAll(faceToCheck, pointI)
7573 {
7574 if
7575 (
7576 mirrorPointLabels.found(faceToCheck[pointI])
7577 )
7578 {
7579 faceToCheck[pointI] =
7580 (
7581 mirrorPointLabels[faceToCheck[pointI]]
7582 );
7583 }
7584 }
7585 }
7586 }
7587 }
7588 }
7589 }
7590 else
7591 {
7592 // Point renumbering of entities connected to mirror points
7593 forAllIter(Map<label>, mirrorPointLabels, pIter)
7594 {
7595 const labelList& pEdges = pointEdges_[pIter()];
7597 forAll(pEdges, edgeI)
7598 {
7599 // Renumber this edge.
7600 edge& edgeToCheck = edges_[pEdges[edgeI]];
7602 forAll(edgeToCheck, pointI)
7603 {
7604 if (edgeToCheck[pointI] == pIter.key())
7605 {
7606 edgeToCheck[pointI] = pIter();
7607 }
7608 }
7610 // Also renumber faces connected to this edge.
7611 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
7613 forAll(eFaces, faceI)
7614 {
7615 face& faceToCheck = faces_[eFaces[faceI]];
7617 forAll(faceToCheck, pointI)
7618 {
7619 if (faceToCheck[pointI] == pIter.key())
7620 {
7621 faceToCheck[pointI] = pIter();
7622 }
7623 }
7624 }
7625 }
7626 }
7627 }
7629 // Now that we're done with the internal faces, remove them.
7630 forAll(internalFaces, faceI)
7631 {
7632 removeFace(internalFaces[faceI]);
7633 }
7635 // Remove old internal edges as well.
7636 forAllIter(Map<label>, mirrorEdgeLabels[1], eIter)
7637 {
7638 if (eIter.key() != eIter())
7639 {
7640 removeEdge(eIter.key());
7641 }
7642 }
7644 // Set the flag
7645 topoChangeFlag_ = true;
7646 }
7648 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
7650 } // End namespace Foam
7652 // ************************************************************************* //