[foam-extend-3.2.git] / src / dynamicMesh / dynamicFvMesh / dynamicTopoFvMesh / dynamicTopoFvMeshCoupled.C
| blob | 4f0281f0ab63a435cc4007a3784c86df171ae44e |
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 Class
26 dynamicTopoFvMesh
28 Description
29 Functions specific to coupled connectivity
31 Author
32 Sandeep Menon
33 University of Massachusetts Amherst
34 All rights reserved
36 \*---------------------------------------------------------------------------*/
38 #include "dynamicTopoFvMesh.H"
40 #include "Time.H"
41 #include "eMesh.H"
42 #include "triFace.H"
43 #include "volFields.H"
44 #include "changeMap.H"
45 #include "topoMapper.H"
46 #include "coupledInfo.H"
47 #include "matchPoints.H"
48 #include "SortableList.H"
49 #include "motionSolver.H"
50 #include "surfaceFields.H"
51 #include "globalMeshData.H"
52 #include "cyclicPolyPatch.H"
53 #include "fvMeshDistribute.H"
54 #include "faceSetAlgorithm.H"
55 #include "cellSetAlgorithm.H"
56 #include "coordinateSystem.H"
57 #include "processorPolyPatch.H"
58 #include "decompositionMethod.H"
59 #include "mapDistributePolyMesh.H"
61 namespace Foam
62 {
64 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
66 defineTypeNameAndDebug(coupleMap, 0);
68 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
70 // Set coupled modification
71 void dynamicTopoFvMesh::setCoupledModification() const
72 {
73 coupledModification_ = true;
74 }
77 // Unset coupled modification
78 void dynamicTopoFvMesh::unsetCoupledModification() const
79 {
80 coupledModification_ = false;
81 }
84 // Initialize the coupled stack
85 void dynamicTopoFvMesh::initCoupledStack
86 (
87 const labelHashSet& entities,
88 bool useEntities
89 )
90 {
91 // Clear existing lists/stacks.
92 stack(0).clear();
94 if (useEntities)
95 {
96 bool emptyEntity;
98 // Initialize the stack with entries
99 // in the labelHashSet and return
100 forAllConstIter(labelHashSet, entities, eIter)
101 {
102 // Add only valid entities
103 emptyEntity =
104 (
105 twoDMesh_ ?
106 faces_[eIter.key()].empty() :
107 edgeFaces_[eIter.key()].empty()
108 );
110 if (emptyEntity)
111 {
112 continue;
113 }
115 stack(0).insert(eIter.key());
116 }
118 if (debug > 3 && Pstream::parRun())
119 {
120 Pout<< nl << "Entity stack size: " << stack(0).size() << endl;
122 if (debug > 4)
123 {
124 // Write out stack entities
125 labelList stackElements(stack(0).size(), -1);
127 forAll(stackElements, elemI)
128 {
129 stackElements[elemI] = stack(0)[elemI];
130 }
132 label elemType = twoDMesh_ ? 2 : 1;
134 writeVTK
135 (
136 "entityStack_"
137 + Foam::name(Pstream::myProcNo()),
138 stackElements,
139 elemType
140 );
141 }
142 }
144 return;
145 }
147 // Loop though boundary faces and check whether
148 // they belong to master/slave coupled patches.
149 for (label faceI = nOldInternalFaces_; faceI < faces_.size(); faceI++)
150 {
151 // Add only valid faces
152 if (faces_[faceI].empty())
153 {
154 continue;
155 }
157 label pIndex = whichPatch(faceI);
159 if (pIndex == -1)
160 {
161 continue;
162 }
164 // Check if this is a locally coupled master face.
165 if (patchCoupling_.size())
166 {
167 if (patchCoupling_(pIndex))
168 {
169 // Is patch cyclic?
170 bool addIndex = true;
172 if (isA<cyclicPolyPatch>(boundaryMesh()[pIndex]))
173 {
174 // Check if this is a cyclic master face
175 const coupleMap& cMap = patchCoupling_[pIndex].map();
176 const Map<label>& fMap = cMap.entityMap(coupleMap::FACE);
178 if (!fMap.found(faceI))
179 {
180 addIndex = false;
181 }
182 }
184 // Add this to the coupled modification stack.
185 if (addIndex)
186 {
187 if (twoDMesh_)
188 {
189 stack(0).push(faceI);
190 }
191 else
192 {
193 const labelList& mfEdges = faceEdges_[faceI];
195 forAll(mfEdges, edgeI)
196 {
197 // Add this to the coupled modification stack.
198 stack(0).push(mfEdges[edgeI]);
199 }
200 }
201 }
202 }
203 }
205 // Check if this is a processor patch.
206 label neiProcID = getNeighbourProcessor(pIndex);
208 if (neiProcID > -1)
209 {
210 // Check if this is a master processor patch.
211 if (neiProcID > Pstream::myProcNo())
212 {
213 // Add this to the coupled modification stack.
214 if (twoDMesh_)
215 {
216 stack(0).push(faceI);
217 }
218 else
219 {
220 const label edgeEnum = coupleMap::EDGE;
221 const labelList& mfEdges = faceEdges_[faceI];
223 forAll(mfEdges, edgeI)
224 {
225 label eIndex = mfEdges[edgeI];
227 // Need to avoid this edge if it is
228 // talking to lower-ranked processors.
229 bool permissible = true;
231 forAll(procIndices_, pI)
232 {
233 if (procIndices_[pI] > Pstream::myProcNo())
234 {
235 continue;
236 }
238 // Fetch reference to subMeshes
239 const coupledInfo& recvMesh = recvMeshes_[pI];
240 const coupleMap& rcMap = recvMesh.map();
242 if (rcMap.findSlave(edgeEnum, eIndex) > -1)
243 {
244 permissible = false;
245 break;
246 }
247 }
249 if (permissible)
250 {
251 stack(0).push(eIndex);
252 }
253 }
254 }
255 }
256 }
257 }
259 if (debug > 3 && Pstream::parRun())
260 {
261 Pout<< nl << "Coupled stack size: " << stack(0).size() << endl;
263 if (debug > 4)
264 {
265 // Write out stack entities
266 labelList stackElements(stack(0).size(), -1);
268 forAll(stackElements, elemI)
269 {
270 stackElements[elemI] = stack(0)[elemI];
271 }
273 label elemType = twoDMesh_ ? 2 : 1;
275 writeVTK
276 (
277 "coupledStack_"
278 + Foam::name(Pstream::myProcNo()),
279 stackElements,
280 elemType
281 );
282 }
283 }
284 }
287 // Execute load balancing operations on the mesh
288 void dynamicTopoFvMesh::executeLoadBalancing
289 (
290 const dictionary& balanceDict
291 )
292 {
293 if (!Pstream::parRun())
294 {
295 return;
296 }
298 bool enabled = readBool(balanceDict.lookup("enabled"));
300 if (!enabled)
301 {
302 return;
303 }
305 // Read interval
306 label interval = readLabel(balanceDict.lookup("interval"));
308 // Are we at the interval?
309 if ((interval < 0) || (time().timeIndex() % interval != 0))
310 {
311 return;
312 }
314 if (debug)
315 {
316 Info<< " void dynamicTopoFvMesh::executeLoadBalancing() :"
317 << " Executing load-balancing operations on the mesh."
318 << endl;
320 // Write out patch boundaries before re-distribution
321 for (label pI = 0; pI < nPatches_; pI++)
322 {
323 label neiProcNo = getNeighbourProcessor(pI);
325 if (neiProcNo == -1)
326 {
327 continue;
328 }
330 writeVTK
331 (
332 "procPatchFacesNew_"
333 + Foam::name(Pstream::myProcNo())
334 + '_'
335 + Foam::name(neiProcNo),
336 identity(patchSizes_[pI]) + patchStarts_[pI],
337 2
338 );
340 writeVTK
341 (
342 "procPatchFacesOld_"
343 + Foam::name(Pstream::myProcNo())
344 + '_'
345 + Foam::name(neiProcNo),
346 identity(patchSizes_[pI]) + patchStarts_[pI],
347 2, false, true
348 );
349 }
350 }
352 // Read decomposition options and initialize
353 autoPtr<decompositionMethod> decomposerPtr
354 (
355 decompositionMethod::New
356 (
357 balanceDict,
358 (*this)
359 )
360 );
362 // Alias for convenience...
363 decompositionMethod& decomposer = decomposerPtr();
365 if (!decomposer.parallelAware())
366 {
367 FatalErrorIn("void dynamicTopoFvMesh::executeLoadBalancing()")
368 << "You have selected decomposition method "
369 << decomposer.typeName
370 << " which is not parallel aware." << endl
371 << "Please select one that is (hierarchical, parMetis)"
372 << exit(FatalError);
373 }
375 // Calculate a merge-distance
376 const boundBox& box = polyMesh::bounds();
377 scalar mergeTol = readScalar(balanceDict.lookup("mergeTol"));
379 // Compute mergeDist
380 scalar mergeDist = mergeTol * box.mag();
382 // Compute write tolerance
383 scalar writeTol =
384 (
385 std::pow
386 (
387 scalar(10.0),
388 -scalar(IOstream::defaultPrecision())
389 )
390 );
392 Info<< nl
393 << " Overall mesh bounding box : " << box << nl
394 << " Relative tolerance : " << mergeTol << nl
395 << " Absolute matching distance : " << mergeDist << nl
396 << endl;
398 // Check for compatibility
399 if (time().writeFormat() == IOstream::ASCII && mergeTol < writeTol)
400 {
401 FatalErrorIn("void dynamicTopoFvMesh::executeLoadBalancing()")
402 << " Your current settings specify ASCII writing with "
403 << IOstream::defaultPrecision() << " digits precision." << nl
404 << " Your mergeTol (" << mergeTol << ") is finer than this." << nl
405 << " Remedies: " << nl
406 << " - Change writeFormat to binary" << nl
407 << " - Increase writePrecision" << nl
408 << " - Adjust the merge tolerance (mergeTol)" << endl
409 << exit(FatalError);
410 }
412 // Decompose the mesh and obtain distribution
413 labelList cellDistribution
414 (
415 decomposer.decompose
416 (
417 primitiveMesh::cellCentres(),
418 scalarField(nCells_, 1)
419 )
420 );
422 // Set the coupledModification switch
423 setCoupledModification();
425 // Initialize the mesh distribution engine
426 fvMeshDistribute distributor(*this, mergeDist);
428 // Re-distribute mesh according to new decomposition
429 distributor.distribute(cellDistribution);
431 // Reset the coupledModification switch
432 unsetCoupledModification();
434 // Set old / new sizes
435 nPoints_ = nOldPoints_ = primitiveMesh::nPoints();
436 nEdges_ = nOldEdges_ = 0;
437 nFaces_ = nOldFaces_ = primitiveMesh::nFaces();
438 nCells_ = nOldCells_ = primitiveMesh::nCells();
439 nInternalFaces_ = nOldInternalFaces_ = primitiveMesh::nInternalFaces();
440 nInternalEdges_ = nOldInternalEdges_ = 0;
442 // Now re-initialize all connectivity structures
443 oldPoints_ = polyMesh::points();
444 points_ = polyMesh::points();
445 faces_ = polyMesh::faces();
446 owner_ = polyMesh::faceOwner();
447 neighbour_ = polyMesh::faceNeighbour();
448 cells_ = primitiveMesh::cells();
450 // Check the size of owner / neighbour
451 if (owner_.size() != neighbour_.size())
452 {
453 // Size up to number of faces
454 neighbour_.setSize(nFaces_, -1);
455 }
457 const polyBoundaryMesh& boundary = polyMesh::boundaryMesh();
459 // Initialize patch-size information
460 nPatches_ = boundary.size();
462 oldPatchSizes_.setSize(nPatches_, 0);
463 oldPatchStarts_.setSize(nPatches_, -1);
464 oldEdgePatchSizes_.setSize(nPatches_, 0);
465 oldEdgePatchStarts_.setSize(nPatches_, -1);
466 oldPatchNMeshPoints_.setSize(nPatches_, -1);
468 patchSizes_.setSize(nPatches_, 0);
469 patchStarts_.setSize(nPatches_, -1);
470 edgePatchSizes_.setSize(nPatches_, 0);
471 edgePatchStarts_.setSize(nPatches_, -1);
472 patchNMeshPoints_.setSize(nPatches_, -1);
474 for (label i = 0; i < nPatches_; i++)
475 {
476 patchNMeshPoints_[i] = boundary[i].meshPoints().size();
477 oldPatchSizes_[i] = patchSizes_[i] = boundary[i].size();
478 oldPatchStarts_[i] = patchStarts_[i] = boundary[i].start();
479 oldPatchNMeshPoints_[i] = patchNMeshPoints_[i];
480 }
482 // Clear pointers and re-initialize
483 mapper_.clear();
484 eMeshPtr_.clear();
485 motionSolver_.clear();
486 lengthEstimator_.clear();
488 // Initialize edge-related connectivity structures
489 initEdges();
491 // Load the mesh-motion solver
492 loadMotionSolver();
494 // Load the field-mapper
495 loadFieldMapper();
497 // Set sizes for the reverse maps
498 reversePointMap_.setSize(nPoints_, -7);
499 reverseEdgeMap_.setSize(nEdges_, -7);
500 reverseFaceMap_.setSize(nFaces_, -7);
501 reverseCellMap_.setSize(nCells_, -7);
503 // Load the length-scale estimator,
504 // and read refinement options
505 loadLengthScaleEstimator();
507 // Clear parallel structures
508 procIndices_.clear();
509 sendMeshes_.clear();
510 recvMeshes_.clear();
511 }
514 // Identify coupled patches.
515 // - Also builds global shared point information.
516 // - Returns true if no coupled patches were found.
517 bool dynamicTopoFvMesh::identifyCoupledPatches()
518 {
519 bool coupledPatchesAbsent = true;
521 // Check if patches are explicitly coupled
522 if (patchCoupling_.size())
523 {
524 coupledPatchesAbsent = false;
525 }
527 // Maintain a separate list of processor IDs in procIndices.
528 // This is done because this sub-domain may talk to processors
529 // that share only edges/points.
530 if (!Pstream::parRun() || procIndices_.size())
531 {
532 return coupledPatchesAbsent;
533 }
535 // Prepare a list of points for sub-mesh creation.
536 // - Obtain global shared-points information, if necessary.
537 const polyBoundaryMesh& boundary = boundaryMesh();
539 // Set flag as default for parallel runs
540 coupledPatchesAbsent = false;
542 // Fetch the list of global points from polyMesh.
543 // - This should be the first evaluation of globalData,
544 // so this involves global communication
545 const globalMeshData& gData = polyMesh::globalData();
547 const labelList& spA = gData.sharedPointAddr();
548 const labelList& spL = gData.sharedPointLabels();
550 labelListList spB(Pstream::nProcs(), labelList(0));
552 if (gData.nGlobalPoints())
553 {
554 if (debug)
555 {
556 Info<< " dynamicTopoFvMesh::identifyCoupledPatches :"
557 << " Found " << gData.nGlobalPoints()
558 << " global points." << endl;
559 }
561 // Send others my addressing.
562 for (label proc = 0; proc < Pstream::nProcs(); proc++)
563 {
564 if (proc != Pstream::myProcNo())
565 {
566 // Send number of entities first.
567 meshOps::pWrite(proc, spA.size());
569 // Send the buffer.
570 if (spA.size())
571 {
572 meshOps::pWrite(proc, spA);
573 }
574 }
575 }
577 // Receive addressing from others
578 for (label proc = 0; proc < Pstream::nProcs(); proc++)
579 {
580 if (proc != Pstream::myProcNo())
581 {
582 label procInfoSize = -1;
584 // How many entities am I going to be receiving?
585 meshOps::pRead(proc, procInfoSize);
587 if (procInfoSize)
588 {
589 // Size the receive buffer.
590 spB[proc].setSize(procInfoSize, -1);
592 // Schedule for receipt.
593 meshOps::pRead(proc, spB[proc]);
594 }
595 }
596 }
597 }
598 else
599 if (debug)
600 {
601 Info<< " dynamicTopoFvMesh::identifyCoupledPatches :"
602 << " Did not find any global points." << endl;
603 }
605 Map<label> immNeighbours;
606 labelListList procPoints(Pstream::nProcs());
608 // Track globally shared points
609 List<DynamicList<labelPair> > globalProcPoints
610 (
611 Pstream::nProcs(),
612 DynamicList<labelPair>(5)
613 );
615 // Insert my immediate neighbours into the list.
616 forAll(boundary, pI)
617 {
618 if (isA<processorPolyPatch>(boundary[pI]))
619 {
620 const processorPolyPatch& pp =
621 (
622 refCast<const processorPolyPatch>(boundary[pI])
623 );
625 label neiProcNo = pp.neighbProcNo();
627 // Insert all boundary points.
628 procPoints[neiProcNo] = pp.meshPoints();
630 // Keep track of immediate neighbours.
631 immNeighbours.insert(neiProcNo, pI);
632 }
633 }
635 if (gData.nGlobalPoints())
636 {
637 // Wait for all transfers to complete.
638 meshOps::waitForBuffers();
640 // Now loop through all processor addressing, and check if
641 // any labels coincide with my global shared points.
642 // If this is true, we need to be talking to that neighbour
643 // as well (if not already).
644 for (label proc = 0; proc < Pstream::nProcs(); proc++)
645 {
646 if (proc == Pstream::myProcNo())
647 {
648 continue;
649 }
651 bool foundGlobalMatch = false;
653 // Fetch reference to buffer
654 const labelList& procBuffer = spB[proc];
656 forAll(procBuffer, pointI)
657 {
658 forAll(spA, pointJ)
659 {
660 if (spA[pointJ] == procBuffer[pointI])
661 {
662 // Make an entry, if one wasn't made already
663 if (findIndex(procPoints[proc], spL[pointJ]) == -1)
664 {
665 globalProcPoints[proc].append
666 (
667 labelPair(spL[pointJ], spA[pointJ])
668 );
669 }
671 foundGlobalMatch = true;
673 break;
674 }
675 }
676 }
678 if (!immNeighbours.found(proc))
679 {
680 if (foundGlobalMatch && debug)
681 {
682 Pout<< " dynamicTopoFvMesh::identifyCoupledPatches :"
683 << " Additionally talking to processor: "
684 << proc << endl;
685 }
686 }
687 }
688 }
690 // Estimate an initial size
691 procIndices_.setSize(Pstream::nProcs());
693 // Patch sub meshes need to be prepared in ascending
694 // order of neighbouring processors.
695 label nTotalProcs = 0;
697 forAll(procPoints, procI)
698 {
699 if (procPoints[procI].size())
700 {
701 procIndices_[nTotalProcs++] = procI;
702 }
704 // Check for point / edge coupling
705 if (globalProcPoints[procI].size() && !procPoints[procI].size())
706 {
707 procIndices_[nTotalProcs++] = procI;
708 }
709 }
711 // Shrink to actual size
712 procIndices_.setSize(nTotalProcs);
714 // Size the PtrLists.
715 sendMeshes_.setSize(nTotalProcs);
716 recvMeshes_.setSize(nTotalProcs);
718 // Create send/recv patch meshes, and copy
719 // the list of points for each processor.
720 forAll(procIndices_, pI)
721 {
722 label proc = procIndices_[pI], master = -1, slave = -1;
724 // For processors that have only point / edge coupling
725 // specify an invalid patch ID for now
726 label patchID = immNeighbours.found(proc) ? immNeighbours[proc] : -1;
728 if (proc < Pstream::myProcNo())
729 {
730 master = proc;
731 slave = Pstream::myProcNo();
732 }
733 else
734 {
735 master = Pstream::myProcNo();
736 slave = proc;
737 }
739 sendMeshes_.set
740 (
741 pI,
742 new coupledInfo
743 (
744 *this, // Reference to this mesh
745 twoDMesh_, // 2D or 3D
746 false, // Not local
747 true, // Sent to neighbour
748 patchID, // Patch index
749 master, // Master index
750 slave // Slave index
751 )
752 );
754 sendMeshes_[pI].map().subMeshPoints() =
755 (
756 procPoints[procIndices_[pI]]
757 );
759 sendMeshes_[pI].map().globalProcPoints() =
760 (
761 globalProcPoints[procIndices_[pI]]
762 );
764 recvMeshes_.set
765 (
766 pI,
767 new coupledInfo
768 (
769 *this, // Reference to this mesh
770 twoDMesh_, // 2D or 3D
771 false, // Not local
772 false, // Not sent to neighbour
773 patchID, // Patch index
774 master, // Master index
775 slave // Slave index
776 )
777 );
779 recvMeshes_[pI].map().subMeshPoints() =
780 (
781 procPoints[procIndices_[pI]]
782 );
784 recvMeshes_[pI].map().globalProcPoints() =
785 (
786 globalProcPoints[procIndices_[pI]]
787 );
788 }
790 if (debug > 3)
791 {
792 Pout<< " identifyCoupledPatches :"
793 << " Talking to processors: "
794 << procIndices_ << endl;
796 forAll(procIndices_, pI)
797 {
798 label proc = procIndices_[pI];
800 const List<labelPair>* gppPtr = NULL;
802 // Write out subMeshPoints as a VTK
803 if (proc < Pstream::myProcNo())
804 {
805 const coupleMap& cMap = sendMeshes_[pI].map();
807 writeVTK
808 (
809 "subMeshPoints_" +
810 Foam::name(Pstream::myProcNo()) + "to" +
811 Foam::name(proc),
812 cMap.subMeshPoints(),
813 0
814 );
816 // Fetch reference
817 gppPtr = &(cMap.globalProcPoints());
818 }
819 else
820 {
821 const coupleMap& cMap = sendMeshes_[pI].map();
823 writeVTK
824 (
825 "subMeshPoints_" +
826 Foam::name(Pstream::myProcNo()) + "to" +
827 Foam::name(proc),
828 cMap.subMeshPoints(),
829 0
830 );
832 // Fetch reference
833 gppPtr = &(cMap.globalProcPoints());
834 }
836 // Write out globalProcPoints as a VTK
837 const List<labelPair>& gpp = *gppPtr;
839 if (gpp.size())
840 {
841 labelList gpPoints(gpp.size(), -1);
843 forAll(gpp, pointI)
844 {
845 gpPoints[pointI] = gpp[pointI].first();
846 }
848 writeVTK
849 (
850 "globalProcPoints_" +
851 Foam::name(Pstream::myProcNo()) + "to" +
852 Foam::name(proc),
853 gpPoints,
854 0
855 );
856 }
857 }
858 }
860 return coupledPatchesAbsent;
861 }
864 // Read coupled patch information from dictionary.
865 void dynamicTopoFvMesh::readCoupledPatches()
866 {
867 const polyBoundaryMesh& boundary = boundaryMesh();
869 // Clear list
870 patchCoupling_.clear();
872 if (dict_.found("coupledPatches") || mandatory_)
873 {
874 // Size the initial list
875 patchCoupling_.setSize(boundary.size());
877 const dictionary& coupledPatches = dict_.subDict("coupledPatches");
879 // Determine master and slave patches
880 forAllConstIter(dictionary, coupledPatches, dIter)
881 {
882 const dictionary& dictI = dIter().dict();
884 // Lookup the master / slave patches
885 word masterPatch = dictI.lookup("master");
886 word slavePatch = dictI.lookup("slave");
888 // Determine patch indices
889 label mPatch = boundary.findPatchID(masterPatch);
890 label sPatch = boundary.findPatchID(slavePatch);
892 if (debug)
893 {
894 Info<< " Adding master: " << masterPatch << " : " << mPatch
895 << " with slave: " << slavePatch << " : " << sPatch
896 << endl;
897 }
899 // Add to the list if entries are legitimate
900 if
901 (
902 mPatch != sPatch &&
903 boundary[mPatch].size() == boundary[sPatch].size()
904 )
905 {
906 // Check whether patches are associated with zones.
907 Switch specifyZones
908 (
909 dictI.lookup("specifyZones")
910 );
912 label mZone = -1, sZone = -1;
914 if (specifyZones)
915 {
916 const faceZoneMesh& faceZones = polyMesh::faceZones();
918 mZone = faceZones.findZoneID
919 (
920 dictI.lookup("masterZone")
921 );
923 sZone = faceZones.findZoneID
924 (
925 dictI.lookup("slaveZone")
926 );
927 }
929 // Configure with regIOobject for check-in
930 coupleMap cMap
931 (
932 IOobject
933 (
934 "coupleMap_"
935 + Foam::name(mPatch)
936 + "_To_"
937 + Foam::name(sPatch)
938 + "_Local",
939 time().timeName(),
940 *this,
941 IOobject::READ_IF_PRESENT,
942 IOobject::AUTO_WRITE,
943 true
944 ),
945 twoDMesh_,
946 true,
947 false,
948 mPatch,
949 mPatch,
950 sPatch
951 );
953 // Set the pointer
954 patchCoupling_.set
955 (
956 mPatch,
957 new coupledInfo
958 (
959 *this,
960 cMap,
961 mZone,
962 sZone
963 )
964 );
965 }
966 else
967 {
968 FatalErrorIn("void dynamicTopoFvMesh::readCoupledPatches()")
969 << " Coupled patches are either wrongly specified,"
970 << " or the sizes don't match." << nl
971 << " Master: " << mPatch << ":" << masterPatch
972 << " Size: " << boundary[mPatch].size() << nl
973 << " Slave: " << sPatch << ":" << slavePatch
974 << " Size: " << boundary[sPatch].size() << nl
975 << abort(FatalError);
976 }
977 }
978 }
980 // Loop through boundaries and add any cyclic patches
981 forAll(boundary, patchI)
982 {
983 if (!isA<cyclicPolyPatch>(boundary[patchI]))
984 {
985 continue;
986 }
988 // Size up patchCoupling, if necessary
989 if (patchCoupling_.empty())
990 {
991 patchCoupling_.setSize(boundary.size());
992 }
994 if (debug)
995 {
996 Info<< " Adding cyclic: " << patchI
997 << " : " << boundary[patchI].name()
998 << endl;
999 }
1001 // Configure with regIOobject for check-in
1002 coupleMap cMap
1003 (
1004 IOobject
1005 (
1006 "coupleMap_"
1007 + Foam::name(patchI)
1008 + "_To_"
1009 + Foam::name(patchI)
1010 + "_Local",
1011 time().timeName(),
1012 *this,
1013 IOobject::READ_IF_PRESENT,
1014 IOobject::AUTO_WRITE,
1015 true
1016 ),
1017 twoDMesh_,
1018 true,
1019 false,
1020 patchI,
1021 patchI,
1022 patchI
1023 );
1025 // Set the pointer
1026 patchCoupling_.set
1027 (
1028 patchI,
1029 new coupledInfo(*this, cMap, -1, -1)
1030 );
1031 }
1032 }
1035 // Initialize coupled patch connectivity for topology modifications.
1036 // - Send and receive sub meshes for processor patches.
1037 // - Made static because this function may be run in a separate thread.
1038 void dynamicTopoFvMesh::initCoupledConnectivity
1039 (
1040 void *argument
1041 )
1042 {
1043 // Recast the argument
1044 dynamicTopoFvMesh *mesh = static_cast<dynamicTopoFvMesh*>(argument);
1046 // Identify coupled patches.
1047 if (mesh->identifyCoupledPatches())
1048 {
1049 return;
1050 }
1052 // Build and send patch sub-meshes.
1053 mesh->buildProcessorPatchMeshes();
1055 // Build inital maps for locally coupled patches.
1056 mesh->buildLocalCoupledMaps();
1058 // Build maps for coupled processor patches.
1059 mesh->buildProcessorCoupledMaps();
1060 }
1063 // Move coupled subMesh points
1064 void dynamicTopoFvMesh::moveCoupledSubMeshes()
1065 {
1066 if (!Pstream::parRun())
1067 {
1068 return;
1069 }
1071 if (debug)
1072 {
1073 Info<< " void dynamicTopoFvMesh::moveCoupledSubMeshes() :"
1074 << " Moving points for coupled subMeshes."
1075 << endl;
1076 }
1078 forAll(procIndices_, pI)
1079 {
1080 label proc = procIndices_[pI];
1082 const coupledInfo& sPM = sendMeshes_[pI];
1083 const coupledInfo& rPM = recvMeshes_[pI];
1085 // Fetch the coupleMap
1086 const coupleMap& scMap = sPM.map();
1087 const coupleMap& rcMap = rPM.map();
1089 Map<label>& pointMap = scMap.entityMap(coupleMap::POINT);
1091 // Fill point buffers
1092 pointField& pBuffer = scMap.pointBuffer();
1093 pointField& opBuffer = scMap.oldPointBuffer();
1095 forAllConstIter(Map<label>, pointMap, pIter)
1096 {
1097 pBuffer[pIter.key()] = points_[pIter()];
1098 opBuffer[pIter.key()] = oldPoints_[pIter()];
1099 }
1101 // Buffers have already been allocated
1102 // to the right size, so just transfer points
1104 // Send point buffers to neighbour
1105 meshOps::pWrite(proc, scMap.pointBuffer());
1106 meshOps::pWrite(proc, scMap.oldPointBuffer());
1108 // Receive point buffers from neighbour
1109 meshOps::pRead(proc, rcMap.pointBuffer());
1110 meshOps::pRead(proc, rcMap.oldPointBuffer());
1111 }
1113 // Wait for transfers to complete before moving on
1114 meshOps::waitForBuffers();
1116 // Set points in mesh
1117 forAll(procIndices_, pI)
1118 {
1119 // Fetch non-const reference to patchSubMesh
1120 coupledInfo& rPM = recvMeshes_[pI];
1122 // Fetch the coupleMap
1123 const coupleMap& rcMap = rPM.map();
1125 dynamicTopoFvMesh& mesh = rPM.subMesh();
1126 const polyBoundaryMesh& boundary = mesh.boundaryMesh();
1128 // Set points / oldPoints
1129 mesh.points_ = rcMap.pointBuffer();
1130 mesh.oldPoints_ = rcMap.oldPointBuffer();
1132 labelList patchSizes(boundary.size(), -1);
1133 labelList patchStarts(boundary.size(), -1);
1135 forAll(boundary, patchI)
1136 {
1137 patchSizes[patchI] = boundary[patchI].size();
1138 patchStarts[patchI] = boundary[patchI].start();
1139 }
1141 // Reset underlying mesh.
1142 // - Use null lists for addressing to avoid over-writes
1143 // - Specify non-valid boundary to avoid globalData creation
1144 mesh.resetPrimitives
1145 (
1146 xferCopy(rcMap.oldPointBuffer()),
1147 Xfer<faceList>::null(),
1148 Xfer<labelList>::null(),
1149 Xfer<labelList>::null(),
1150 patchSizes,
1151 patchStarts,
1152 false
1153 );
1154 }
1155 }
1158 // Insert the cells around the coupled master entity to the mesh
1159 // - Returns a changeMap with a type specifying:
1160 // 1: Insertion was successful
1161 // -1: Insertion failed
1162 // -2: Failed because entity was being handled elsewhere
1163 // - The changeMap index specifies the converted mIndex.
1164 const changeMap dynamicTopoFvMesh::insertCells(const label mIndex)
1165 {
1166 // Prepare the changeMaps
1167 changeMap map;
1169 // Maintain face counts for each inserted cell
1170 Map<label> nCellFaces;
1172 // Track inserted entities
1173 List<Map<label> > pointsToInsert(procIndices_.size());
1174 List<Map<label> > edgesToInsert(procIndices_.size());
1175 List<Map<label> > facesToInsert(procIndices_.size());
1176 List<Map<label> > cellsToInsert(procIndices_.size());
1178 // Track converted entities
1179 List<Map<label> > edgesToConvert(procIndices_.size());
1180 List<Map<label> > facesToConvert(procIndices_.size());
1182 // Track edge / face patch information
1183 List<Map<label> > masterConvertEdgePatch(procIndices_.size());
1184 List<Map<label> > masterConvertFacePatch(procIndices_.size());
1186 Map<label> createPatch;
1187 labelList slaveConvertPatch(procIndices_.size(), -1);
1188 List<Map<Pair<point> > > convertPatchPoints(procIndices_.size());
1190 // First check to ensure that this case can be handled
1191 forAll(procIndices_, pI)
1192 {
1193 const label cplEnum = twoDMesh_ ? coupleMap::FACE : coupleMap::EDGE;
1195 // Fetch reference to maps
1196 const coupleMap& cMap = recvMeshes_[pI].map();
1198 // Does a coupling exist?
1199 if (cMap.findSlave(cplEnum, mIndex) == -1)
1200 {
1201 continue;
1202 }
1204 // If this entity is being handled elsewhere, bail out
1205 if (procIndices_[pI] < Pstream::myProcNo())
1206 {
1207 map.type() = -2;
1209 return map;
1210 }
1211 }
1213 const polyBoundaryMesh& boundary = boundaryMesh();
1215 // Agglomerate cells from surrounding subMeshes
1216 forAll(procIndices_, pI)
1217 {
1218 const label cplEnum = twoDMesh_ ? coupleMap::FACE : coupleMap::EDGE;
1220 // Fetch reference to maps
1221 const coupleMap& cMap = recvMeshes_[pI].map();
1222 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
1224 label sIndex = -1;
1226 if ((sIndex = cMap.findSlave(cplEnum, mIndex)) == -1)
1227 {
1228 continue;
1229 }
1231 // Fetch references
1232 Map<label>& procCellMap = cellsToInsert[pI];
1234 if (twoDMesh_)
1235 {
1236 // Insert the owner cell
1237 procCellMap.insert(mesh.owner_[sIndex], -1);
1238 }
1239 else
1240 {
1241 // Insert all cells connected to this edge
1242 const labelList& eFaces = mesh.edgeFaces_[sIndex];
1244 forAll(eFaces, faceI)
1245 {
1246 const label own = mesh.owner_[eFaces[faceI]];
1247 const label nei = mesh.neighbour_[eFaces[faceI]];
1249 if (!procCellMap.found(own))
1250 {
1251 procCellMap.insert(own, -1);
1252 }
1254 if (!procCellMap.found(nei) && (nei != -1))
1255 {
1256 procCellMap.insert(nei, -1);
1257 }
1258 }
1259 }
1261 // Find a patch that talks to this processor
1262 label neiProcPatch = -1;
1264 forAll(boundary, patchI)
1265 {
1266 if (getNeighbourProcessor(patchI) == procIndices_[pI])
1267 {
1268 neiProcPatch = patchI;
1269 break;
1270 }
1271 }
1273 // Prepare information about inserted / converted faces and edges
1274 const polyBoundaryMesh& slaveBoundary = mesh.boundaryMesh();
1276 // Set the slave conversion patch
1277 forAll(slaveBoundary, patchI)
1278 {
1279 if (mesh.getNeighbourProcessor(patchI) == Pstream::myProcNo())
1280 {
1281 slaveConvertPatch[pI] = patchI;
1282 break;
1283 }
1284 }
1286 forAllIter(Map<label>, procCellMap, cIter)
1287 {
1288 const cell& checkCell = mesh.cells_[cIter.key()];
1290 forAll(checkCell, fI)
1291 {
1292 label mfIndex = -1, sfIndex = checkCell[fI];
1293 label sfPatch = mesh.whichPatch(sfIndex);
1295 // Check whether a face mapping exists for this face
1296 if ((mfIndex = cMap.findMaster(coupleMap::FACE, sfIndex)) > -1)
1297 {
1298 // Mark as a candidate for conversion
1299 facesToInsert[pI].insert(sfIndex, mfIndex);
1300 facesToConvert[pI].insert(sfIndex, mfIndex);
1301 }
1302 else
1303 {
1304 // Mark for insertion. Some faces may be duplicated
1305 // across processors, but this will be dealt
1306 // with at a later stage.
1307 facesToInsert[pI].insert(sfIndex, -1);
1308 }
1310 // Check face points for coupling
1311 const face& sFace = mesh.faces_[sfIndex];
1313 // Configure points which need to be inserted
1314 forAll(sFace, pointI)
1315 {
1316 // Skip if added already
1317 if (pointsToInsert[pI].found(sFace[pointI]))
1318 {
1319 continue;
1320 }
1322 // Check for coupled points
1323 pointsToInsert[pI].insert
1324 (
1325 sFace[pointI],
1326 cMap.findMaster
1327 (
1328 coupleMap::POINT,
1329 sFace[pointI]
1330 )
1331 );
1332 }
1334 // Loop through edges and check whether edge-mapping exists
1335 const labelList& sfEdges = mesh.faceEdges_[sfIndex];
1337 forAll(sfEdges, edgeI)
1338 {
1339 const label seIndex = sfEdges[edgeI];
1340 const edge& sEdge = mesh.edges_[seIndex];
1342 // Meshes in 2D don't have edge-mapping, so check
1343 // point maps instead. If either point doesn't exist,
1344 // this is an edge that needs to be inserted.
1345 label cMs = cMap.findMaster(coupleMap::POINT, sEdge[0]);
1346 label cMe = cMap.findMaster(coupleMap::POINT, sEdge[1]);
1348 if (!edgesToInsert[pI].found(seIndex))
1349 {
1350 edgesToInsert[pI].insert(seIndex, -1);
1351 }
1353 // If both points have maps, this is a conversion edge
1354 if (cMs > -1 && cMe > -1)
1355 {
1356 if (!edgesToConvert[pI].found(seIndex))
1357 {
1358 edgesToConvert[pI].insert(seIndex, -1);
1359 }
1360 }
1362 // Add a patch entry as well
1363 if (masterConvertEdgePatch[pI].found(seIndex))
1364 {
1365 continue;
1366 }
1368 label edgePatch = -1;
1369 label sePatch = mesh.whichEdgePatch(seIndex);
1370 label neiProcNo = mesh.getNeighbourProcessor(sePatch);
1372 // Determine patch
1373 if (sePatch == -1)
1374 {
1375 // Slave edge was an interior one
1376 edgePatch = neiProcPatch;
1377 }
1378 else
1379 if (sePatch == (slaveBoundary.size() - 1))
1380 {
1381 // The 'defaultPatch'
1382 edgePatch = neiProcPatch;
1383 }
1384 else
1385 if (neiProcNo > -1)
1386 {
1387 if (neiProcNo == Pstream::myProcNo())
1388 {
1389 // Set the master edge patch
1390 edgePatch = neiProcPatch;
1391 }
1392 else
1393 {
1394 // If this other processor is
1395 // lesser-ranked, bail out.
1396 if (neiProcNo < Pstream::myProcNo())
1397 {
1398 map.type() = -2;
1400 return map;
1401 }
1402 else
1403 if (debug > 3)
1404 {
1405 Pout<< nl << nl
1406 << " Edge: " << seIndex
1407 << " :: " << mesh.edges_[seIndex]
1408 << " is talking to processor: "
1409 << neiProcNo << endl;
1410 }
1412 // Find an appropriate boundary patch
1413 forAll(boundary, patchI)
1414 {
1415 label eP = getNeighbourProcessor(patchI);
1417 if (eP == neiProcNo)
1418 {
1419 edgePatch = patchI;
1420 break;
1421 }
1422 }
1424 if (edgePatch == -1)
1425 {
1426 if (debug > 3)
1427 {
1428 Pout<< nl << nl
1429 << " No direct contact with"
1430 << " processor: " << neiProcNo
1431 << " so adding to: " << neiProcPatch
1432 << endl;
1433 }
1435 // Add this to neiProcPatch instead.
1436 edgePatch = neiProcPatch;
1437 }
1438 }
1439 }
1440 else
1441 {
1442 // Physical type
1443 edgePatch = sePatch;
1444 }
1446 if (edgePatch == -1)
1447 {
1448 // Write out the edge
1449 mesh.writeVTK("seEdge", seIndex, 1);
1451 Pout<< " Could not find correct patch info: " << nl
1452 << " seIndex: " << seIndex << nl
1453 << " sePatch: " << sePatch << nl
1454 << " neiProcPatch: " << neiProcPatch << nl
1455 << " Patch Name: " <<
1456 (
1457 sePatch > -1 ?
1458 slaveBoundary[sePatch].name() :
1459 "Internal"
1460 )
1461 << abort(FatalError);
1462 }
1464 // Add the patch entry
1465 masterConvertEdgePatch[pI].insert
1466 (
1467 seIndex,
1468 edgePatch
1469 );
1470 }
1472 // Determine patch
1473 if (masterConvertFacePatch[pI].found(sfIndex))
1474 {
1475 continue;
1476 }
1478 label facePatch = -1;
1479 label neiProcNo = mesh.getNeighbourProcessor(sfPatch);
1481 if (sfPatch == -1)
1482 {
1483 // Slave face was an interior one
1484 facePatch = neiProcPatch;
1485 }
1486 else
1487 if (sfPatch == (slaveBoundary.size() - 1))
1488 {
1489 // The 'defaultPatch'
1490 facePatch = neiProcPatch;
1491 }
1492 else
1493 if (neiProcNo > -1)
1494 {
1495 if (neiProcNo == Pstream::myProcNo())
1496 {
1497 // Check to see if this face was converted before
1498 if
1499 (
1500 findIndex
1501 (
1502 cMap.entityOperations(),
1503 coupleMap::CONVERT_PATCH
1504 ) > -1
1505 )
1506 {
1507 // Loop through points and check for match
1508 const pointField& pts = cMap.moveNewPoints();
1509 const face& fCheck = mesh.faces_[sfIndex];
1510 const point fC = fCheck.centre(mesh.points_);
1512 scalar tol = mag(fC - mesh.points_[fCheck[0]]);
1514 forAll(pts, ptI)
1515 {
1516 scalar dist = mag(fC - pts[ptI]);
1518 if (dist < (1e-4 * tol))
1519 {
1520 // Face was converted before
1521 if (debug > 3)
1522 {
1523 Pout<< nl << nl
1524 << " Face: " << sfIndex
1525 << " :: " << mesh.faces_[sfIndex]
1526 << " was converted before."
1527 << endl;
1528 }
1530 map.type() = -2;
1532 return map;
1533 }
1534 }
1535 }
1537 // Set the master face patch
1538 facePatch = neiProcPatch;
1539 }
1540 else
1541 {
1542 // If this other processor is
1543 // lesser-ranked, bail out.
1544 if (neiProcNo < Pstream::myProcNo())
1545 {
1546 map.type() = -2;
1548 return map;
1549 }
1550 else
1551 if (debug > 3)
1552 {
1553 Pout<< nl << nl
1554 << " Face: " << sfIndex
1555 << " :: " << mesh.faces_[sfIndex]
1556 << " is talking to processor: "
1557 << neiProcNo << endl;
1558 }
1560 // Find an appropriate boundary patch
1561 forAll(boundary, patchI)
1562 {
1563 label fP = getNeighbourProcessor(patchI);
1565 if (fP == neiProcNo)
1566 {
1567 facePatch = patchI;
1568 break;
1569 }
1570 }
1572 // Specify a convert-patch operation
1573 // for later, if this operation succeeds
1574 const face& sfCheck = mesh.faces_[sfIndex];
1576 point nfC = sfCheck.centre(mesh.points_);
1577 point ofC = sfCheck.centre(mesh.oldPoints_);
1579 // Are we talking to this processor already?
1580 label prI = findIndex(procIndices_, neiProcNo);
1582 // Check for face-contact
1583 if (facePatch == -1)
1584 {
1585 if (debug > 3)
1586 {
1587 Pout<< nl << nl
1588 << " No face contact with"
1589 << " processor: " << neiProcNo
1590 << endl;
1591 }
1593 if (prI == -1)
1594 {
1595 Pout<< " No contact with: " << neiProcNo
1596 << abort(FatalError);
1597 }
1599 label pIdx =
1600 (
1601 sendMeshes_[prI].map().patchIndex()
1602 );
1604 // Add a patch creation order, if necessary
1605 if (pIdx == -1 && !createPatch.found(neiProcNo))
1606 {
1607 createPatch.insert(neiProcNo, -1);
1608 }
1610 // Specify a value for facePatch:
1611 // - Specify a value that we can use
1612 // to back out the patch after creation
1613 // - Also needs to bypass failure check
1614 facePatch = (-2 - neiProcNo);
1615 }
1617 // Add to the list
1618 convertPatchPoints[prI].insert
1619 (
1620 sfIndex,
1621 Pair<point>(nfC, ofC)
1622 );
1623 }
1624 }
1625 else
1626 {
1627 // Physical type
1628 facePatch = sfPatch;
1629 }
1631 if (facePatch == -1)
1632 {
1633 // Write out the face
1634 mesh.writeVTK("sfFace", sfIndex, 2);
1636 Pout<< " Could not find correct patch info: " << nl
1637 << " sfIndex: " << sfIndex << nl
1638 << " sfPatch: " << sfPatch << nl
1639 << " neiProcPatch: " << neiProcPatch << nl
1640 << " slavePatch: " <<
1641 (
1642 sfPatch > -1 ?
1643 slaveBoundary[sfPatch].name() :
1644 "Internal"
1645 )
1646 << abort(FatalError);
1647 }
1649 // Add the patch entry
1650 masterConvertFacePatch[pI].insert
1651 (
1652 sfIndex,
1653 facePatch
1654 );
1655 }
1656 }
1657 }
1659 // Perform a few debug calls before modifications
1660 if (debug > 1)
1661 {
1662 Pout<< nl << nl
1663 << "Inserting cell(s) around coupled "
1664 << (twoDMesh_ ? "face: " : "edge: ") << mIndex
1665 << endl;
1666 }
1668 // Check to see if any new processor
1669 // patches need to be created
1670 forAllIter(Map<label>, createPatch, procIter)
1671 {
1672 label neiProcNo = procIter.key();
1674 // Create a new processor patch
1675 procIter() = createProcessorPatch(neiProcNo);
1677 forAll(masterConvertFacePatch, pI)
1678 {
1679 Map<label>& convertMap = masterConvertFacePatch[pI];
1681 forAllIter(Map<label>, convertMap, mIter)
1682 {
1683 if (mIter() < 0)
1684 {
1685 // Back out the neighbouring processor ID
1686 label proc = Foam::mag(mIter() + 2);
1688 if (proc == neiProcNo)
1689 {
1690 // Set the index
1691 mIter() = procIter();
1692 }
1693 }
1694 }
1695 }
1697 // Find index in processor list
1698 label pI = findIndex(procIndices_, neiProcNo);
1700 if (pI == -1)
1701 {
1702 Pout<< " Could not find index for processor: " << neiProcNo
1703 << " in indices: " << procIndices_
1704 << abort(FatalError);
1705 }
1707 // Create patch on subMesh
1708 dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
1710 mesh.createProcessorPatch(Pstream::myProcNo());
1711 }
1713 // Build a list of mapping entities on this processor
1714 DynamicList<label> mapCells(10);
1716 if (twoDMesh_)
1717 {
1718 label own = owner_[mIndex];
1720 mapCells.append(own);
1721 }
1722 else
1723 {
1724 const labelList& eFaces = edgeFaces_[mIndex];
1726 forAll(eFaces, faceI)
1727 {
1728 label own = owner_[eFaces[faceI]];
1729 label nei = neighbour_[eFaces[faceI]];
1731 if (findIndex(mapCells, own) == -1)
1732 {
1733 mapCells.append(own);
1734 }
1736 if (nei > -1)
1737 {
1738 if (findIndex(mapCells, nei) == -1)
1739 {
1740 mapCells.append(nei);
1741 }
1742 }
1743 }
1744 }
1746 // Loop through insertion cells and
1747 // create an equivalent on this mesh
1748 forAll(procIndices_, pI)
1749 {
1750 const label cplEnum = twoDMesh_ ? coupleMap::FACE : coupleMap::EDGE;
1752 // Fetch reference to maps
1753 const coupleMap& cMap = recvMeshes_[pI].map();
1754 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
1756 label sIndex = -1;
1758 if ((sIndex = cMap.findSlave(cplEnum, mIndex)) == -1)
1759 {
1760 continue;
1761 }
1763 // Fetch references
1764 Map<label>& procCellMap = cellsToInsert[pI];
1766 forAllIter(Map<label>, procCellMap, cIter)
1767 {
1768 const cell& checkCell = mesh.cells_[cIter.key()];
1770 scalar sLengthScale = -1.0;
1772 if (edgeRefinement_)
1773 {
1774 sLengthScale = mesh.lengthScale_[cIter.key()];
1775 }
1777 // Add an empty cell for now, and update
1778 // with face information at a later stage.
1779 cIter() = insertCell(cell(checkCell.size()), sLengthScale);
1781 // Initialize a face counter
1782 nCellFaces.insert(cIter(), 0);
1784 if (debug > 3)
1785 {
1786 Pout<< " Map cell: " << cIter()
1787 << " for cell: " << cIter.key()
1788 << " on proc: " << procIndices_[pI]
1789 << endl;
1790 }
1792 // Add this cell to the map.
1793 map.addCell(cIter(), mapCells);
1795 // Set basic mapping for this cell
1796 setCellMapping(cIter(), mapCells);
1797 }
1799 // Push operation for the slave into coupleMap
1800 cMap.pushOperation
1801 (
1802 sIndex,
1803 coupleMap::REMOVE_CELL
1804 );
1805 }
1807 // Build a list of edges that need to be converted to interior.
1808 // - Do this by looking at edges of master face conversion candidates.
1809 // - Some edges may not need conversion, but deal with this later.
1810 forAll(facesToConvert, pI)
1811 {
1812 // Fetch reference to subMesh
1813 const coupleMap& cMap = recvMeshes_[pI].map();
1814 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
1816 const Map<label>& procFaceMap = facesToConvert[pI];
1818 forAllConstIter(Map<label>, procFaceMap, fIter)
1819 {
1820 const labelList& mfEdges = faceEdges_[fIter()];
1821 const labelList& sfEdges = mesh.faceEdges_[fIter.key()];
1823 forAll(sfEdges, edgeI)
1824 {
1825 if (edgesToInsert[pI][sfEdges[edgeI]] != -1)
1826 {
1827 // Already mapped this edge. Move on.
1828 continue;
1829 }
1831 // Configure the comparison edge
1832 const edge& sEdge = mesh.edges_[sfEdges[edgeI]];
1834 label cMs = cMap.findMaster(coupleMap::POINT, sEdge[0]);
1835 label cMe = cMap.findMaster(coupleMap::POINT, sEdge[1]);
1837 edge cEdge(cMs, cMe);
1839 forAll(mfEdges, edgeJ)
1840 {
1841 const edge& mEdge = edges_[mfEdges[edgeJ]];
1843 if (mEdge == cEdge)
1844 {
1845 edgesToInsert[pI][sfEdges[edgeI]] = mfEdges[edgeJ];
1846 edgesToConvert[pI][sfEdges[edgeI]] = mfEdges[edgeJ];
1847 break;
1848 }
1849 }
1850 }
1851 }
1852 }
1854 // Loop through all insertion points, and merge if necessary
1855 scalar mergeTol =
1856 (
1857 twoDMesh_ ? 0.0 : 1e-4 * mag(edges_[mIndex].vec(points_))
1858 );
1860 // Insert all points
1861 forAll(pointsToInsert, pI)
1862 {
1863 Map<label>& procPointMapI = pointsToInsert[pI];
1865 // Fetch reference to subMesh
1866 const coupleMap& cMapI = recvMeshes_[pI].map();
1867 const dynamicTopoFvMesh& meshI = recvMeshes_[pI].subMesh();
1869 forAllIter(Map<label>, procPointMapI, pItI)
1870 {
1871 if (pItI() != -1)
1872 {
1873 continue;
1874 }
1876 const point& pointI = meshI.points_[pItI.key()];
1878 bool foundMerge = false;
1880 forAll(pointsToInsert, pJ)
1881 {
1882 if (pJ == pI)
1883 {
1884 continue;
1885 }
1887 Map<label>& procPointMapJ = pointsToInsert[pJ];
1889 // Fetch reference to subMesh
1890 const coupleMap& cMapJ = recvMeshes_[pJ].map();
1891 const dynamicTopoFvMesh& meshJ = recvMeshes_[pJ].subMesh();
1893 // Compare points with this processor
1894 forAllIter(Map<label>, procPointMapJ, pItJ)
1895 {
1896 const point& pointJ = meshJ.points_[pItJ.key()];
1898 if (mag(pointI - pointJ) < mergeTol)
1899 {
1900 if (pItJ() == -1)
1901 {
1902 // Make a merge entry
1903 label mergePointIndex =
1904 (
1905 insertPoint
1906 (
1907 pointJ,
1908 meshJ.oldPoints_[pItJ.key()],
1909 labelList(1, -1)
1910 )
1911 );
1913 pItI() = mergePointIndex;
1914 pItJ() = mergePointIndex;
1916 // Update maps for the new point
1917 cMapI.mapSlave
1918 (
1919 coupleMap::POINT,
1920 pItI(), pItI.key()
1921 );
1923 cMapI.mapMaster
1924 (
1925 coupleMap::POINT,
1926 pItI.key(), pItI()
1927 );
1929 cMapJ.mapSlave
1930 (
1931 coupleMap::POINT,
1932 pItJ(), pItJ.key()
1933 );
1935 cMapJ.mapMaster
1936 (
1937 coupleMap::POINT,
1938 pItJ.key(), pItJ()
1939 );
1941 if (debug > 2)
1942 {
1943 Pout<< " Map point: " << mergePointIndex
1944 << " pointI: " << pItI.key()
1945 << " procI: " << procIndices_[pI]
1946 << " pointJ: " << pItJ.key()
1947 << " procJ: " << procIndices_[pJ]
1948 << endl;
1949 }
1951 // Add this point to the map.
1952 map.addPoint(mergePointIndex);
1953 }
1954 else
1955 {
1956 // Point appears to have been inserted
1957 // by a previous operation. Map to it.
1958 if (debug > 2)
1959 {
1960 Pout<< " Inserted point: " << pItJ()
1961 << " pointI: " << pItI.key()
1962 << " procI: " << procIndices_[pI]
1963 << " pointJ: " << pItJ.key()
1964 << " procJ: " << procIndices_[pJ]
1965 << endl;
1966 }
1968 // Set the entry
1969 pItI() = pItJ();
1971 // Update maps for the point
1972 cMapI.mapSlave
1973 (
1974 coupleMap::POINT,
1975 pItI(), pItI.key()
1976 );
1978 cMapI.mapMaster
1979 (
1980 coupleMap::POINT,
1981 pItI.key(), pItI()
1982 );
1983 }
1985 foundMerge = true;
1986 break;
1987 }
1988 }
1989 }
1991 if (foundMerge)
1992 {
1993 continue;
1994 }
1996 // Add a new unique point
1997 label newPointIndex =
1998 (
1999 insertPoint
2000 (
2001 pointI,
2002 meshI.oldPoints_[pItI.key()],
2003 labelList(1, -1)
2004 )
2005 );
2007 // Set the entry
2008 pItI() = newPointIndex;
2010 // Update maps for the new point
2011 cMapI.mapSlave(coupleMap::POINT, pItI(), pItI.key());
2012 cMapI.mapMaster(coupleMap::POINT, pItI.key(), pItI());
2014 if (debug > 2)
2015 {
2016 Pout<< " Map point: " << newPointIndex
2017 << " for point: " << pItI.key()
2018 << " on proc: " << procIndices_[pI]
2019 << endl;
2020 }
2022 // Add this point to the map.
2023 map.addPoint(newPointIndex);
2024 }
2025 }
2027 // Occassionally, inserted edges may already be present.
2028 // Ensure that edges are not added twice.
2029 if (!twoDMesh_)
2030 {
2031 forAll(facesToInsert, pI)
2032 {
2033 // Fetch reference to subMesh
2034 const coupleMap& cMap = recvMeshes_[pI].map();
2035 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2037 const Map<label>& procFaceMap = facesToInsert[pI];
2039 forAllConstIter(Map<label>, procFaceMap, fIter)
2040 {
2041 const labelList& sfEdges = mesh.faceEdges_[fIter.key()];
2043 forAll(sfEdges, edgeI)
2044 {
2045 label seIndex = sfEdges[edgeI];
2047 if (edgesToInsert[pI][seIndex] != -1)
2048 {
2049 // Already mapped this edge. Move on.
2050 continue;
2051 }
2053 label cMe = cMap.findMaster(coupleMap::EDGE, seIndex);
2055 if (cMe > -1)
2056 {
2057 if (debug > 2)
2058 {
2059 Pout<< " Found master edge: " << cMe
2060 << " :: " << edges_[cMe]
2061 << " for slave edge: " << seIndex
2062 << " :: " << mesh.edges_[seIndex]
2063 << " on proc: " << procIndices_[pI]
2064 << endl;
2065 }
2067 edgesToInsert[pI][seIndex] = cMe;
2068 edgesToConvert[pI][seIndex] = cMe;
2070 continue;
2071 }
2073 // Check for point-only coupling
2074 const edge& sEdge = mesh.edges_[seIndex];
2076 edge cEdge
2077 (
2078 cMap.findMaster(coupleMap::POINT, sEdge[0]),
2079 cMap.findMaster(coupleMap::POINT, sEdge[1])
2080 );
2082 if (cEdge[0] > -1 && cEdge[1] > -1)
2083 {
2084 label meIndex = -1;
2086 // Look at pointEdges info for a boundary edge
2087 const labelList& pEdges = pointEdges_[cEdge[0]];
2089 forAll(pEdges, edgeJ)
2090 {
2091 if (edges_[pEdges[edgeJ]] == cEdge)
2092 {
2093 if (whichEdgePatch(pEdges[edgeJ]) > -1)
2094 {
2095 meIndex = pEdges[edgeJ];
2096 break;
2097 }
2098 }
2099 }
2101 if (meIndex > -1)
2102 {
2103 if (debug > 2)
2104 {
2105 Pout<< " Found master edge: " << meIndex
2106 << " :: " << edges_[meIndex]
2107 << " for slave edge: " << seIndex
2108 << " :: " << mesh.edges_[seIndex]
2109 << " on proc: " << procIndices_[pI]
2110 << endl;
2111 }
2113 edgesToInsert[pI][seIndex] = meIndex;
2114 edgesToConvert[pI][seIndex] = meIndex;
2115 }
2116 }
2117 }
2118 }
2119 }
2120 }
2122 // Write out some debug info
2123 if (debug > 2)
2124 {
2125 forAll(cellsToInsert, pI)
2126 {
2127 label procIndex = procIndices_[pI];
2129 // Fetch reference to subMesh
2130 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2132 // Define a convenience output macro
2133 # define writeOutputVTK(fName, eMap, eType, insert) \
2134 { \
2135 labelList entities(eMap.size()); \
2136 \
2137 label nEnt = 0; \
2138 \
2139 forAllConstIter(Map<label>, eMap, eIter) \
2140 { \
2141 if (eIter() == -1 && insert) \
2142 { \
2143 entities[nEnt++] = eIter.key(); \
2144 } \
2145 else \
2146 if (eIter() > -1 && !insert) \
2147 { \
2148 entities[nEnt++] = eIter.key(); \
2149 } \
2150 } \
2151 \
2152 entities.setSize(nEnt); \
2153 \
2154 if (entities.size()) \
2155 { \
2156 mesh.writeVTK \
2157 ( \
2158 word(fName) + '_' \
2159 + Foam::name(mIndex) + '_' \
2160 + Foam::name(procIndex), \
2161 entities, \
2162 eType, false, true \
2163 ); \
2164 } \
2165 }
2167 writeOutputVTK("edgesToConvert", edgesToConvert[pI], 1, false)
2168 writeOutputVTK("facesToConvert", facesToConvert[pI], 2, false)
2169 writeOutputVTK("pointsToConvert", pointsToInsert[pI], 0, false)
2170 writeOutputVTK("pointsToInsert", pointsToInsert[pI], 0, true)
2171 writeOutputVTK("edgesToInsert", edgesToInsert[pI], 1, true)
2172 writeOutputVTK("facesToInsert", facesToInsert[pI], 2, true)
2173 writeOutputVTK("cellsToInsert", cellsToInsert[pI], 3, false)
2174 }
2175 }
2177 // Add edges to the mesh, noting that all
2178 // required points have already been added.
2179 forAll(edgesToInsert, pI)
2180 {
2181 // Fetch reference to subMesh
2182 const coupleMap& cMap = recvMeshes_[pI].map();
2183 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2185 Map<label>& procEdgeMap = edgesToInsert[pI];
2187 forAllIter(Map<label>, procEdgeMap, eIter)
2188 {
2189 // Skip if the edge is a conversion candidate
2190 if (eIter() != -1)
2191 {
2192 continue;
2193 }
2195 const edge& sEdge = mesh.edges_[eIter.key()];
2197 edge newEdge
2198 (
2199 cMap.findMaster(coupleMap::POINT, sEdge[0]),
2200 cMap.findMaster(coupleMap::POINT, sEdge[1])
2201 );
2203 // Ensure that this edge hasn't been added before
2204 label neIndex = -1;
2205 bool foundDuplicate = false;
2207 const List<objectMap>& addedEdges = map.addedEdgeList();
2209 forAll(addedEdges, edgeI)
2210 {
2211 neIndex = addedEdges[edgeI].index();
2213 if (edges_[neIndex] == newEdge)
2214 {
2215 foundDuplicate = true;
2216 break;
2217 }
2218 }
2220 if (foundDuplicate)
2221 {
2222 // Note the duplicate index for later
2223 eIter() = neIndex;
2225 continue;
2226 }
2228 // Insert edge with null edgeFaces for now.
2229 // This can be corrected later.
2230 eIter() =
2231 (
2232 insertEdge
2233 (
2234 masterConvertEdgePatch[pI][eIter.key()],
2235 newEdge,
2236 labelList(0)
2237 )
2238 );
2240 if (debug > 2)
2241 {
2242 Pout<< " Map edge: " << eIter() << "::" << newEdge
2243 << " for edge: " << eIter.key() << "::" << sEdge
2244 << endl;
2245 }
2247 // Add this edge to the map.
2248 map.addEdge(eIter());
2249 }
2250 }
2252 // Add faces to the mesh, noting that all required
2253 // points and edges have already been added.
2254 forAll(facesToInsert, pI)
2255 {
2256 // Fetch reference to subMesh and coupleMap
2257 const coupleMap& cMapI = recvMeshes_[pI].map();
2258 const dynamicTopoFvMesh& meshI = recvMeshes_[pI].subMesh();
2260 Map<label>& procFaceMap = facesToInsert[pI];
2262 forAllIter(Map<label>, procFaceMap, fI)
2263 {
2264 // Skip if the face is a conversion candidate
2265 if (fI() != -1)
2266 {
2267 continue;
2268 }
2270 const face& sFaceI = meshI.faces_[fI.key()];
2271 const labelList& sfEdges = meshI.faceEdges_[fI.key()];
2273 // Configure new / comparison faces
2274 face nF(sFaceI.size()), cF(sFaceI.size());
2275 labelList nFaceEdges(sfEdges.size());
2277 // Configure points from map
2278 forAll(nF, pointI)
2279 {
2280 nF[pointI] =
2281 (
2282 cMapI.findMaster
2283 (
2284 coupleMap::POINT,
2285 sFaceI[pointI]
2286 )
2287 );
2288 }
2290 // This may be a face shared by two other processors.
2291 // Ensure that duplicates are added only once.
2292 label dupeOwnCell = -1;
2293 bool foundDuplicate = false, foundInsertedDuplicate = false;
2295 forAll(facesToInsert, pJ)
2296 {
2297 if (pJ == pI)
2298 {
2299 continue;
2300 }
2302 // Fetch reference to subMesh and coupleMap
2303 const coupleMap& cMapJ = recvMeshes_[pJ].map();
2304 const dynamicTopoFvMesh& meshJ = recvMeshes_[pJ].subMesh();
2306 const Map<label>& procFaceMapJ = facesToInsert[pJ];
2308 forAllConstIter(Map<label>, procFaceMapJ, fJ)
2309 {
2310 const face& sFaceJ = meshJ.faces_[fJ.key()];
2312 // Discard dissimilar face sizes
2313 if (sFaceJ.size() != cF.size())
2314 {
2315 continue;
2316 }
2318 // Configure points from map
2319 forAll(cF, pointI)
2320 {
2321 cF[pointI] =
2322 (
2323 cMapJ.findMaster
2324 (
2325 coupleMap::POINT,
2326 sFaceJ[pointI]
2327 )
2328 );
2329 }
2331 if (cF.size() == 3)
2332 {
2333 // Optimized triangular face comparison
2334 if
2335 (
2336 triFace::compare
2337 (
2338 triFace(nF[0], nF[1], nF[2]),
2339 triFace(cF[0], cF[1], cF[2])
2340 )
2341 )
2342 {
2343 foundDuplicate = true;
2344 }
2345 }
2346 else
2347 {
2348 // Regular face compare
2349 if (face::compare(nF, cF))
2350 {
2351 foundDuplicate = true;
2352 }
2353 }
2355 if (foundDuplicate)
2356 {
2357 // Record the owner for posterity
2358 dupeOwnCell =
2359 (
2360 cellsToInsert[pJ][meshJ.owner_[fJ.key()]]
2361 );
2363 // Was the duplicate face inserted before this?
2364 if (fJ() != -1)
2365 {
2366 // Note the duplicate index for later
2367 fI() = fJ();
2369 // If patch conversion entries were made,
2370 // remove them as well
2371 if
2372 (
2373 convertPatchPoints[pI].found(fJ.key())
2374 && convertPatchPoints[pJ].found(fI.key())
2375 )
2376 {
2377 convertPatchPoints[pI].erase(fJ.key());
2378 convertPatchPoints[pJ].erase(fI.key());
2379 }
2381 foundInsertedDuplicate = true;
2382 }
2384 break;
2385 }
2386 }
2387 }
2389 // Face was inserted before, so don't insert again
2390 if (foundInsertedDuplicate)
2391 {
2392 continue;
2393 }
2395 // Configure edges from edgesToInsert
2396 forAll(sfEdges, edgeI)
2397 {
2398 label meIndex = -1;
2400 // Configure with the appropriate edge
2401 if (edgesToInsert[pI].found(sfEdges[edgeI]))
2402 {
2403 meIndex = edgesToInsert[pI][sfEdges[edgeI]];
2404 }
2406 if (meIndex == -1)
2407 {
2408 // Something is wrong here.
2409 Pout<< " Could not find correspondence for slave edge: "
2410 << sfEdges[edgeI]
2411 << ":: " << meshI.edges_[sfEdges[edgeI]]
2412 << nl << " mIndex: " << mIndex
2413 << abort(FatalError);
2414 }
2416 nFaceEdges[edgeI] = meIndex;
2417 }
2419 // Determine patch, owner and neighbour for this face
2420 label nPatch = -1, nOwner = -1, nNeighbour = -1;
2422 const polyBoundaryMesh& slaveBoundary = meshI.boundaryMesh();
2424 label sfPatch = meshI.whichPatch(fI.key());
2425 label sFaceOwn = meshI.owner_[fI.key()];
2426 label sFaceNei = meshI.neighbour_[fI.key()];
2428 label mFaceOwn =
2429 (
2430 cellsToInsert[pI].found(sFaceOwn) ?
2431 cellsToInsert[pI][sFaceOwn] : -1
2432 );
2434 label mFaceNei =
2435 (
2436 cellsToInsert[pI].found(sFaceNei) ?
2437 cellsToInsert[pI][sFaceNei] : -1
2438 );
2440 // If a duplicate face was found, over-ride neighbour.
2441 // Face-flipping will be taken care of automatically.
2442 if (foundDuplicate)
2443 {
2444 mFaceNei = dupeOwnCell;
2445 }
2447 if (mFaceOwn != -1 && mFaceNei == -1)
2448 {
2449 // Boundary face already has correct orientation
2450 nOwner = mFaceOwn;
2451 nNeighbour = -1;
2453 // Determine patch
2454 nPatch = masterConvertFacePatch[pI][fI.key()];
2455 }
2456 else
2457 if (mFaceOwn == -1 && mFaceNei != -1)
2458 {
2459 // Boundary face is inverted. Flip it
2460 nF = nF.reverseFace();
2461 nOwner = mFaceNei;
2462 nNeighbour = -1;
2464 // Determine patch
2465 nPatch = masterConvertFacePatch[pI][fI.key()];
2466 }
2467 else
2468 if (mFaceOwn != -1 && mFaceNei != -1)
2469 {
2470 // Interior face. Check if a flip is necessary.
2471 if (mFaceNei < mFaceOwn)
2472 {
2473 nF = nF.reverseFace();
2474 nOwner = mFaceNei;
2475 nNeighbour = mFaceOwn;
2476 }
2477 else
2478 {
2479 nOwner = mFaceOwn;
2480 nNeighbour = mFaceNei;
2481 }
2483 nPatch = -1;
2484 }
2485 else
2486 if (mFaceOwn == -1 && mFaceNei == -1)
2487 {
2488 // Something is wrong here.
2489 Pout<< "Could not find correct owner / neighbour info: " << nl
2490 << " Face: " << nF << nl
2491 << " Owner: " << mFaceOwn << nl
2492 << " Neighbour: " << mFaceNei << nl
2493 << " - Slave Face: " << sFaceI << nl
2494 << " - Slave Patch: " << slaveBoundary[sfPatch].name() << nl
2495 << " - Slave Owner: " << sFaceOwn << nl
2496 << " - Slave Neighbour: " << sFaceNei << nl
2497 << abort(FatalError);
2498 }
2500 // Insert the new face
2501 fI() =
2502 (
2503 insertFace
2504 (
2505 nPatch,
2506 nF,
2507 nOwner,
2508 nNeighbour
2509 )
2510 );
2512 // Add the faceEdges entry as well
2513 faceEdges_.append(nFaceEdges);
2515 // Size up edgeFaces for each edge
2516 forAll(nFaceEdges, edgeI)
2517 {
2518 meshOps::sizeUpList
2519 (
2520 fI(),
2521 edgeFaces_[nFaceEdges[edgeI]]
2522 );
2523 }
2525 if (debug > 3)
2526 {
2527 Pout<< " Map face: " << fI() << "::" << nF
2528 << " Own: " << nOwner << " Nei: " << nNeighbour
2529 << " fE: " << nFaceEdges << nl
2530 << " for face: " << fI.key() << "::" << sFaceI
2531 << " Own: " << sFaceOwn << " Nei: " << sFaceNei
2532 << " fE: " << sfEdges
2533 << endl;
2534 }
2536 // Add this face to the map.
2537 map.addFace(fI());
2539 if (nPatch > -1 && getNeighbourProcessor(nPatch) == -1)
2540 {
2541 // Physical patch on subMesh.
2542 // - Set an invalid number so that
2543 // an entry is made in facesFromFaces,
2544 // while faceParents is empty.
2545 setFaceMapping(fI(), labelList(1, -1));
2546 }
2547 else
2548 {
2549 // Interior / processor face
2550 setFaceMapping(fI());
2551 }
2553 // Update cells
2554 cells_[nOwner][nCellFaces[nOwner]++] = fI();
2556 if (nNeighbour > -1)
2557 {
2558 cells_[nNeighbour][nCellFaces[nNeighbour]++] = fI();
2559 }
2560 }
2561 }
2563 // Loop through conversion faces
2564 forAll(facesToConvert, pI)
2565 {
2566 // Fetch reference to subMesh
2567 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2569 const Map<label>& procFaceMap = facesToConvert[pI];
2571 forAllConstIter(Map<label>, procFaceMap, fIter)
2572 {
2573 label fIndex = fIter();
2575 // Fetch the owner information
2576 label mFaceOwn = owner_[fIndex];
2577 label sFaceOwn = mesh.owner_[fIter.key()];
2578 label newNeighbour = cellsToInsert[pI][sFaceOwn];
2580 // Insert a new internal face.
2581 // Orientation should be correct, because this is a boundary
2582 // face converted to an interior, and adjacent to an added cell.
2583 label newFaceIndex =
2584 (
2585 insertFace
2586 (
2587 -1,
2588 face(faces_[fIndex]),
2589 mFaceOwn,
2590 newNeighbour
2591 )
2592 );
2594 // Add the new faceEdges from the existing face. This may contain
2595 // edges that need to be converted, but that will be done later.
2596 faceEdges_.append(labelList(faceEdges_[fIndex]));
2598 // Update map
2599 map.addFace(newFaceIndex, labelList(1, fIndex));
2601 // Set mapping information
2602 setFaceMapping(newFaceIndex);
2604 // Update the owner cell
2605 meshOps::replaceLabel
2606 (
2607 fIndex,
2608 newFaceIndex,
2609 cells_[mFaceOwn]
2610 );
2612 // Update the neighbour cell
2613 cells_[newNeighbour][nCellFaces[newNeighbour]++] = newFaceIndex;
2615 // Replace edgeFaces with the new face index
2616 const labelList& fEdges = faceEdges_[newFaceIndex];
2618 forAll(fEdges, edgeI)
2619 {
2620 meshOps::replaceLabel
2621 (
2622 fIndex,
2623 newFaceIndex,
2624 edgeFaces_[fEdges[edgeI]]
2625 );
2626 }
2628 // Remove the old boundary face
2629 removeFace(fIndex);
2631 // Update map
2632 map.removeFace(fIndex);
2634 // For 2D meshes, the boundary face gets converted
2635 // to an interior one. Note the index for further operations.
2636 if ((mIndex == fIndex) && twoDMesh_)
2637 {
2638 map.index() = newFaceIndex;
2639 }
2640 }
2641 }
2643 // Sequentially add any convert-patch operations
2644 forAll(convertPatchPoints, pI)
2645 {
2646 // Fetch reference to maps / mesh
2647 const coupleMap& cMap = recvMeshes_[pI].map();
2648 dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2650 if (convertPatchPoints[pI].empty())
2651 {
2652 continue;
2653 }
2655 // Find the appropriate processor patch
2656 // for patch conversion operations
2657 label procPatch = -1;
2659 for (label patchI = 0; patchI < mesh.nPatches_; patchI++)
2660 {
2661 label neiProc = mesh.getNeighbourProcessor(patchI);
2663 if (neiProc == Pstream::myProcNo())
2664 {
2665 procPatch = patchI;
2666 break;
2667 }
2668 }
2670 if (procPatch == -1)
2671 {
2672 Pout<< " * * * insertCells() * * * " << nl
2673 << " Could not find patch on slave processor: "
2674 << procIndices_[pI] << nl
2675 << " convertPatchPoints: " << convertPatchPoints[pI]
2676 << abort(FatalError);
2677 }
2679 forAllConstIter(Map<Pair<point> >, convertPatchPoints[pI], fIter)
2680 {
2681 // Search slave mesh for faces
2682 const point& newCentre = fIter().first();
2683 const point& oldCentre = fIter().second();
2685 label replaceFace = -1;
2687 // Accumulate stats in case of failure
2688 scalar minDist = GREAT;
2689 vector minPoint = vector::zero;
2690 DynamicList<label> checkedFaces;
2692 if (debug > 2)
2693 {
2694 // Reserve for append
2695 checkedFaces.setCapacity(50);
2696 }
2698 // Loop through all boundary faces,
2699 // and compute / compare face centres
2700 label sTot = mesh.faces_.size(), sInt = mesh.nOldInternalFaces_;
2702 for (label faceI = sInt; faceI < sTot; faceI++)
2703 {
2704 const face& fCheck = mesh.faces_[faceI];
2706 if (fCheck.empty())
2707 {
2708 continue;
2709 }
2711 label pIndex = mesh.whichPatch(faceI);
2713 if (mesh.getNeighbourProcessor(pIndex) == -1)
2714 {
2715 continue;
2716 }
2718 // Compute face-centre
2719 vector fC = fCheck.centre(mesh.points_);
2721 // Compute tolerance
2722 scalar tol = mag(mesh.points_[fCheck[0]] - fC);
2723 scalar dist = mag(fC - newCentre);
2725 if (dist < (1e-4 * tol))
2726 {
2727 replaceFace = faceI;
2728 break;
2729 }
2730 else
2731 if (dist < minDist)
2732 {
2733 minPoint = fC;
2734 minDist = dist;
2736 if (debug > 2)
2737 {
2738 checkedFaces.append(faceI);
2739 }
2740 }
2741 }
2743 // Ensure that the face was found
2744 if (replaceFace == -1)
2745 {
2746 mesh.writeVTK
2747 (
2748 "checkedFaces_"
2749 + Foam::name(fIter.key()),
2750 checkedFaces,
2751 2, false, true
2752 );
2754 Pout<< " * * * insertCells() * * * " << nl
2755 << " Convert patch Op failed." << nl
2756 << " Face: " << fIter.key() << nl
2757 << " minPoint: " << minPoint << nl
2758 << " minDistance: " << minDist << nl
2759 << " newCentre: " << newCentre << nl
2760 << " oldCentre: " << oldCentre << nl
2761 << abort(FatalError);
2762 }
2764 // Obtain a copy before adding the new face,
2765 // since the reference might become
2766 // invalid during list resizing.
2767 face newFace = mesh.faces_[replaceFace];
2768 label newOwn = mesh.owner_[replaceFace];
2769 labelList newFaceEdges = mesh.faceEdges_[replaceFace];
2771 label newFaceIndex =
2772 (
2773 mesh.insertFace
2774 (
2775 procPatch,
2776 newFace,
2777 newOwn,
2778 -1
2779 )
2780 );
2782 // Add a faceEdges entry as well.
2783 // Edges don't have to change, since they're
2784 // all on the boundary anyway.
2785 mesh.faceEdges_.append(newFaceEdges);
2787 // changeMap update for the new face is not necessary,
2788 // since it is on the slave subMesh
2790 // Set mapping information
2791 mesh.setFaceMapping(newFaceIndex);
2793 meshOps::replaceLabel
2794 (
2795 replaceFace,
2796 newFaceIndex,
2797 mesh.cells_[newOwn]
2798 );
2800 // Correct edgeFaces with the new face label.
2801 forAll(newFaceEdges, edgeI)
2802 {
2803 meshOps::replaceLabel
2804 (
2805 replaceFace,
2806 newFaceIndex,
2807 mesh.edgeFaces_[newFaceEdges[edgeI]]
2808 );
2809 }
2811 if (debug > 3)
2812 {
2813 Pout<< " Pushing CONVERT_PATCH for face: " << replaceFace
2814 << " :: " << mesh.faces_[replaceFace] << nl
2815 << " with new point: " << fIter().first()
2816 << " and old point: " << fIter().second()
2817 << " on proc: " << procIndices_[pI]
2818 << endl;
2819 }
2821 // Finally remove the face
2822 mesh.removeFace(replaceFace);
2824 cMap.pushOperation
2825 (
2826 replaceFace,
2827 coupleMap::CONVERT_PATCH,
2828 fIter().first(),
2829 fIter().second()
2830 );
2831 }
2832 }
2834 // Loop through conversion edges
2835 forAll(edgesToConvert, pI)
2836 {
2837 // Fetch references
2838 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2840 const Map<label>& procEdgeMap = edgesToConvert[pI];
2842 // Loop through conversion edges
2843 forAllConstIter(Map<label>, procEdgeMap, eIter)
2844 {
2845 label eIndex = eIter(), physPatch = -1;
2846 bool allInterior = true, keepEdge = true;
2848 if (eIndex < 0)
2849 {
2850 // Not a conversion edge. Search added edge list for index
2851 eIndex = edgesToInsert[pI][eIter.key()];
2853 if (eIndex < 0)
2854 {
2855 // Something's wrong
2856 Pout<< " Edge: " << eIter.key()
2857 << " :: " << mesh.edges_[eIter.key()]
2858 << " was never inserted."
2859 << abort(FatalError);
2860 }
2861 }
2863 const labelList& eFaces = edgeFaces_[eIndex];
2865 if (eFaces.size())
2866 {
2867 forAll(eFaces, faceI)
2868 {
2869 if (whichPatch(eFaces[faceI]) > -1)
2870 {
2871 allInterior = false;
2872 break;
2873 }
2874 }
2875 }
2876 else
2877 {
2878 // Edge was deleted before, so skip it
2879 allInterior = false;
2880 }
2882 if (allInterior)
2883 {
2884 physPatch = -1;
2885 keepEdge = false;
2886 }
2887 else
2888 if (eFaces.size())
2889 {
2890 // Check if patches need to be switched
2891 label edgePatch = whichEdgePatch(eIndex);
2893 // Is this edge on a processor patch?
2894 bool edgeOnProc = (getNeighbourProcessor(edgePatch) > -1);
2896 if (edgeOnProc)
2897 {
2898 bool foundProc = false;
2900 forAll(eFaces, faceI)
2901 {
2902 label fPatch = whichPatch(eFaces[faceI]);
2904 if (fPatch == -1)
2905 {
2906 continue;
2907 }
2909 if (getNeighbourProcessor(fPatch) == -1)
2910 {
2911 // Note physical patch for later
2912 physPatch = fPatch;
2913 }
2914 else
2915 {
2916 // Still on a processor patch.
2917 foundProc = true;
2918 break;
2919 }
2920 }
2922 if (foundProc)
2923 {
2924 // Reset physPatch
2925 physPatch = -1;
2926 }
2927 }
2929 if (edgeOnProc && physPatch > -1)
2930 {
2931 // Edge needs to be moved to another patch
2932 keepEdge = false;
2933 }
2934 else
2935 if ((mIndex == eIndex) && !twoDMesh_ && map.index() == -1)
2936 {
2937 // Keep the edge, and note index for later
2938 keepEdge = true;
2939 map.index() = eIndex;
2940 }
2941 }
2943 if (keepEdge)
2944 {
2945 continue;
2946 }
2948 // This edge needs to be converted to interior / other patch
2949 label newEdgeIndex =
2950 (
2951 insertEdge
2952 (
2953 physPatch,
2954 edge(edges_[eIndex]),
2955 labelList(edgeFaces_[eIndex])
2956 )
2957 );
2959 // Update map
2960 map.addEdge(newEdgeIndex, labelList(1, eIndex));
2962 // Update faceEdges information for all connected faces
2963 const labelList& neFaces = edgeFaces_[newEdgeIndex];
2965 forAll(neFaces, faceI)
2966 {
2967 meshOps::replaceLabel
2968 (
2969 eIndex,
2970 newEdgeIndex,
2971 faceEdges_[neFaces[faceI]]
2972 );
2973 }
2975 // Remove the old boundary edge
2976 removeEdge(eIndex);
2978 // Update map
2979 map.removeEdge(eIndex);
2981 // For 3D meshes, the boundary edge gets converted
2982 // to an interior one. Note the index for further operations.
2983 if ((mIndex == eIndex) && !twoDMesh_)
2984 {
2985 map.index() = newEdgeIndex;
2986 }
2988 // Replace the entry in edgesToInsert
2989 edgesToInsert[pI][eIter.key()] = newEdgeIndex;
2990 }
2991 }
2993 List<changeMap> slaveMaps(procIndices_.size());
2995 // Loop through all processors, and remove cells
2996 forAll(cellsToInsert, pI)
2997 {
2998 const Map<label>& procCellMap = cellsToInsert[pI];
3000 if (procCellMap.empty())
3001 {
3002 // Set type to something recognizable
3003 slaveMaps[pI].type() = -7;
3005 continue;
3006 }
3008 // Prepare a list of cells
3009 const labelList cellsToRemove = procCellMap.toc();
3011 // Fetch non-const reference to subMesh
3012 dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
3014 // Now remove cells for this processor
3015 slaveMaps[pI] =
3016 (
3017 mesh.removeCells
3018 (
3019 cellsToRemove,
3020 slaveConvertPatch[pI],
3021 "rc_" + Foam::name(mIndex)
3022 )
3023 );
3024 }
3026 // Now map entities from the removeCells operation
3027 forAll(slaveMaps, pI)
3028 {
3029 const changeMap& slaveMap = slaveMaps[pI];
3031 // Skip empty entities
3032 if (slaveMap.type() == -7)
3033 {
3034 continue;
3035 }
3037 // Fetch references
3038 const coupleMap& cMap = recvMeshes_[pI].map();
3039 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
3041 // Map faces from patches.
3042 const Map<label>& procFaceMap = facesToInsert[pI];
3044 // Check removed faces
3045 const labelList& rsfList = slaveMap.removedFaceList();
3047 forAll(rsfList, indexI)
3048 {
3049 label sfIndex = rsfList[indexI], mfIndex = -1;
3051 // Does an entry exist?
3052 mfIndex = cMap.findMaster(coupleMap::FACE, sfIndex);
3054 if (mfIndex > -1)
3055 {
3056 if (debug > 2)
3057 {
3058 Pout<< " Removing map for master face: "
3059 << mfIndex << " :: " << faces_[mfIndex]
3060 << " Master map: "
3061 << cMap.findSlave(coupleMap::FACE, mfIndex)
3062 << " and slave face: " << sfIndex
3063 << " on proc: " << procIndices_[pI]
3064 << endl;
3065 }
3067 cMap.removeSlave(coupleMap::FACE, sfIndex);
3068 cMap.removeMaster(coupleMap::FACE, mfIndex);
3069 }
3070 }
3072 // Check added faces
3073 const List<objectMap>& asfList = slaveMap.addedFaceList();
3075 forAll(asfList, indexI)
3076 {
3077 label sfIndex = asfList[indexI].index(), mfIndex = -1;
3079 if (mesh.whichPatch(sfIndex) == slaveConvertPatch[pI])
3080 {
3081 // Configure a comparison face
3082 face cFace(mesh.faces_[sfIndex]);
3084 forAll(cFace, pointI)
3085 {
3086 cFace[pointI] =
3087 (
3088 cMap.findMaster
3089 (
3090 coupleMap::POINT,
3091 cFace[pointI]
3092 )
3093 );
3094 }
3096 bool foundMatch = false;
3098 forAllConstIter(Map<label>, procFaceMap, fIter)
3099 {
3100 const face& mFace = faces_[fIter()];
3102 // Discard dissimilar face sizes
3103 if (mFace.size() != cFace.size())
3104 {
3105 continue;
3106 }
3108 if (cFace.size() == 3)
3109 {
3110 // Optimized triangular face comparison
3111 if
3112 (
3113 triFace::compare
3114 (
3115 triFace(mFace[0], mFace[1], mFace[2]),
3116 triFace(cFace[0], cFace[1], cFace[2])
3117 )
3118 )
3119 {
3120 mfIndex = fIter();
3121 foundMatch = true;
3122 break;
3123 }
3124 }
3125 else
3126 {
3127 // Regular face compare
3128 if (face::compare(mFace, cFace))
3129 {
3130 mfIndex = fIter();
3131 foundMatch = true;
3132 break;
3133 }
3134 }
3135 }
3137 if (foundMatch)
3138 {
3139 if (debug > 2)
3140 {
3141 Pout<< " Found master face: " << mfIndex
3142 << " :: " << faces_[mfIndex]
3143 << " for slave face: " << sfIndex
3144 << " :: " << mesh.faces_[sfIndex]
3145 << " on proc: " << procIndices_[pI]
3146 << endl;
3147 }
3149 // Update maps for the new face
3150 cMap.mapSlave(coupleMap::FACE, mfIndex, sfIndex);
3151 cMap.mapMaster(coupleMap::FACE, sfIndex, mfIndex);
3152 }
3153 else
3154 {
3155 // Something is wrong here.
3156 Pout<< " Could not find master face for: " << nl
3157 << " Slave face: " << sfIndex
3158 << " :: " << mesh.faces_[sfIndex] << nl
3159 << " cFace: " << cFace << nl
3160 << " on proc: " << procIndices_[pI] << nl
3161 << abort(FatalError);
3162 }
3163 }
3164 }
3166 // Map edges for 3D meshes
3167 if (!twoDMesh_)
3168 {
3169 // Map edges from patches.
3170 const Map<label>& procEdgeMap = edgesToInsert[pI];
3172 // Check removed edges
3173 const labelList& rseList = slaveMap.removedEdgeList();
3175 forAll(rseList, indexI)
3176 {
3177 label seIndex = rseList[indexI], meIndex = -1;
3179 // Does an entry exist?
3180 meIndex = cMap.findMaster(coupleMap::EDGE, seIndex);
3182 if (meIndex > -1)
3183 {
3184 if (debug > 2)
3185 {
3186 Pout<< " Removing map for master edge: "
3187 << meIndex << " :: " << edges_[meIndex]
3188 << " Master map: "
3189 << cMap.findSlave(coupleMap::EDGE, meIndex)
3190 << " and slave edge: " << seIndex
3191 << " on proc: " << procIndices_[pI]
3192 << endl;
3193 }
3195 cMap.removeSlave(coupleMap::EDGE, seIndex);
3196 cMap.removeMaster(coupleMap::EDGE, meIndex);
3197 }
3198 }
3200 // Check added edges
3201 const List<objectMap>& aseList = slaveMap.addedEdgeList();
3203 forAll(aseList, indexI)
3204 {
3205 label seIndex = aseList[indexI].index(), meIndex = -1;
3207 if (mesh.whichEdgePatch(seIndex) == slaveConvertPatch[pI])
3208 {
3209 // Configure a comparison edge
3210 edge cEdge(mesh.edges_[seIndex]);
3212 cEdge[0] = cMap.findMaster(coupleMap::POINT, cEdge[0]);
3213 cEdge[1] = cMap.findMaster(coupleMap::POINT, cEdge[1]);
3215 bool foundMatch = false;
3217 forAllConstIter(Map<label>, procEdgeMap, eIter)
3218 {
3219 const edge& mEdge = edges_[eIter()];
3221 if (mEdge == cEdge)
3222 {
3223 meIndex = eIter();
3224 foundMatch = true;
3225 break;
3226 }
3227 }
3229 if (foundMatch)
3230 {
3231 if (debug > 2)
3232 {
3233 Pout<< " Found master edge: " << meIndex
3234 << " :: " << edges_[meIndex]
3235 << " for slave edge: " << seIndex
3236 << " :: " << mesh.edges_[seIndex]
3237 << " on proc: " << procIndices_[pI]
3238 << endl;
3239 }
3241 // Update maps for the new edge
3242 cMap.mapSlave(coupleMap::EDGE, meIndex, seIndex);
3243 cMap.mapMaster(coupleMap::EDGE, seIndex, meIndex);
3244 }
3245 else
3246 {
3247 // Something is wrong here.
3248 Pout<< " Could not find master edge for: " << nl
3249 << " Slave edge: " << seIndex
3250 << " :: " << mesh.edges_[seIndex] << nl
3251 << " cEdge: " << cEdge << nl
3252 << " on proc: " << procIndices_[pI] << nl
3253 << abort(FatalError);
3254 }
3255 }
3256 }
3257 }
3259 // Check removed points
3260 const labelList& rspList = slaveMap.removedPointList();
3262 forAll(rspList, indexI)
3263 {
3264 label spIndex = rspList[indexI], mpIndex = -1;
3266 // Does an entry exist?
3267 mpIndex = cMap.findMaster(coupleMap::POINT, spIndex);
3269 if (mpIndex > -1)
3270 {
3271 if (debug > 2)
3272 {
3273 Pout<< " Removing map for master point: " << mpIndex
3274 << " :: new: " << points_[mpIndex]
3275 << " :: old: " << oldPoints_[mpIndex] << nl
3276 << " Master point map: "
3277 << cMap.findSlave(coupleMap::POINT, mpIndex)
3278 << " and slave point: " << spIndex
3279 << " on proc: " << procIndices_[pI]
3280 << endl;
3281 }
3283 cMap.removeSlave(coupleMap::POINT, spIndex);
3284 cMap.removeMaster(coupleMap::POINT, mpIndex);
3285 }
3286 }
3287 }
3289 // Write out cells for debug, if necessary
3290 if (debug > 2 || map.index() < 0)
3291 {
3292 const List<objectMap>& acList = map.addedCellList();
3294 labelList addedCells(acList.size(), -1);
3296 forAll(acList, indexI)
3297 {
3298 addedCells[indexI] = acList[indexI].index();
3299 }
3301 // Write it out
3302 writeVTK
3303 (
3304 "insertCells(" + Foam::name(mIndex) + ')',
3305 addedCells, 3, false, true
3306 );
3307 }
3309 // Set mapping information for any added faces
3310 const List<objectMap>& afList = map.addedFaceList();
3312 forAll(afList, indexI)
3313 {
3314 label mfIndex = afList[indexI].index();
3315 label patch = whichPatch(mfIndex);
3316 label neiProc = getNeighbourProcessor(patch);
3318 // Disregard non-processor faces for coupled mapping
3319 if (neiProc == -1)
3320 {
3321 continue;
3322 }
3324 // Update couple maps, if necessary
3325 const face& mFace = faces_[mfIndex];
3327 forAll(procIndices_, pI)
3328 {
3329 // Skip face if not on this processor
3330 if (neiProc != procIndices_[pI])
3331 {
3332 continue;
3333 }
3335 const coupleMap& cMap = recvMeshes_[pI].map();
3336 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
3338 // Configure a comparison face
3339 face cFace(mFace);
3341 forAll(cFace, pointI)
3342 {
3343 cFace[pointI] =
3344 (
3345 cMap.findSlave
3346 (
3347 coupleMap::POINT,
3348 cFace[pointI]
3349 )
3350 );
3351 }
3353 // Ensure all points are found
3354 if (findIndex(cFace, -1) > -1)
3355 {
3356 continue;
3357 }
3359 // Loop through boundary faces on slave and look for a match
3360 label sTot = mesh.faces_.size(), sInt = mesh.nOldInternalFaces_;
3362 label sfIndex = -1;
3364 for (label sfI = sInt; sfI < sTot; sfI++)
3365 {
3366 const face& sFace = mesh.faces_[sfI];
3368 if (sFace.empty())
3369 {
3370 continue;
3371 }
3373 label pIndex = mesh.whichPatch(sfI);
3375 if (mesh.getNeighbourProcessor(pIndex) == -1)
3376 {
3377 continue;
3378 }
3380 if (twoDMesh_)
3381 {
3382 if (face::compare(cFace, sFace))
3383 {
3384 sfIndex = sfI;
3385 }
3386 }
3387 else
3388 {
3389 if
3390 (
3391 triFace::compare
3392 (
3393 triFace(cFace[0], cFace[1], cFace[2]),
3394 triFace(sFace[0], sFace[1], sFace[2])
3395 )
3396 )
3397 {
3398 sfIndex = sfI;
3399 }
3400 }
3402 // Break out if we're done
3403 if (sfIndex > -1)
3404 {
3405 if (debug > 2)
3406 {
3407 Pout<< " Matched master face: " << mfIndex
3408 << " :: " << mFace
3409 << " with slave: " << sfIndex
3410 << " :: " << sFace
3411 << " on proc: " << procIndices_[pI]
3412 << endl;
3413 }
3415 // Found the slave. Add a map entry
3416 cMap.mapSlave
3417 (
3418 coupleMap::FACE,
3419 mfIndex,
3420 sfIndex
3421 );
3423 cMap.mapMaster
3424 (
3425 coupleMap::FACE,
3426 sfIndex,
3427 mfIndex
3428 );
3430 break;
3431 }
3432 }
3434 if (sfIndex == -1)
3435 {
3436 Pout<< " Failed match for master face: " << mfIndex
3437 << " :: " << mFace
3438 << " on proc: " << procIndices_[pI]
3439 << abort(FatalError);
3440 }
3442 if (twoDMesh_)
3443 {
3444 continue;
3445 }
3447 // Map edges on this face as well
3448 const labelList& mfEdges = faceEdges_[mfIndex];
3449 const labelList& sfEdges = mesh.faceEdges_[sfIndex];
3451 forAll(mfEdges, edgeI)
3452 {
3453 label meIndex = mfEdges[edgeI];
3454 const edge& mEdge = edges_[meIndex];
3456 // Configure a comparison edge
3457 edge cEdge
3458 (
3459 cMap.findSlave(coupleMap::POINT, mEdge[0]),
3460 cMap.findSlave(coupleMap::POINT, mEdge[1])
3461 );
3463 forAll(sfEdges, edgeJ)
3464 {
3465 label seIndex = sfEdges[edgeJ];
3466 const edge& sEdge = mesh.edges_[seIndex];
3468 if (cEdge == sEdge)
3469 {
3470 if (debug > 2)
3471 {
3472 Pout<< " Matched master edge: " << meIndex
3473 << " :: " << mEdge
3474 << " with slave: " << seIndex
3475 << " :: " << sEdge
3476 << " on proc: " << procIndices_[pI]
3477 << endl;
3478 }
3480 // Found the slave. Add a map entry
3481 cMap.mapSlave
3482 (
3483 coupleMap::EDGE,
3484 meIndex,
3485 seIndex
3486 );
3488 cMap.mapMaster
3489 (
3490 coupleMap::EDGE,
3491 seIndex,
3492 meIndex
3493 );
3495 break;
3496 }
3497 }
3498 }
3499 }
3500 }
3502 // Specify that the operation was successful
3503 map.type() = 1;
3505 // Return the changeMap
3506 return map;
3507 }
3510 // Remove the specified cells from the mesh,
3511 // and add internal faces/edges to the specified patch
3512 // - Returns a changeMap with a type specifying:
3513 // 1: Operation was successful
3514 // -1: Operation failed
3515 // - checkOnly performs a feasibility check and returns without modifications.
3516 const changeMap dynamicTopoFvMesh::removeCells
3517 (
3518 const labelList& cList,
3519 const label patch,
3520 const word& rcN,
3521 bool checkOnly
3522 )
3523 {
3524 if (cList.empty() || patch < 0)
3525 {
3526 if (cList.size())
3527 {
3528 writeVTK("removeCellsFailed_" + rcN, cList, 3, false, true);
3529 }
3531 FatalErrorIn
3532 (
3533 "const changeMap dynamicTopoFvMesh::removeCells"
3534 "(const labelList&, const label, const word&, bool)"
3535 )
3536 << " Wrong arguments. " << nl
3537 << " cList: " << cList << nl
3538 << " patch: " << patch << nl
3539 << abort(FatalError);
3540 }
3542 changeMap map;
3544 labelHashSet pointsToRemove, edgesToRemove, facesToRemove;
3545 Map<label> facesToConvert, edgesToConvert;
3547 // First loop through all cells and accumulate
3548 // a set of faces to be removed/converted.
3549 forAll(cList, cellI)
3550 {
3551 const cell& cellToCheck = cells_[cList[cellI]];
3553 forAll(cellToCheck, faceI)
3554 {
3555 label own = owner_[cellToCheck[faceI]];
3556 label nei = neighbour_[cellToCheck[faceI]];
3558 if (nei == -1)
3559 {
3560 if (!facesToRemove.found(cellToCheck[faceI]))
3561 {
3562 facesToRemove.insert(cellToCheck[faceI]);
3563 }
3564 }
3565 else
3566 if
3567 (
3568 (findIndex(cList, own) != -1) &&
3569 (findIndex(cList, nei) != -1)
3570 )
3571 {
3572 if (!facesToRemove.found(cellToCheck[faceI]))
3573 {
3574 facesToRemove.insert(cellToCheck[faceI]);
3575 }
3576 }
3577 else
3578 {
3579 facesToConvert.set(cellToCheck[faceI], -1);
3580 }
3581 }
3582 }
3584 // Add all edges as candidates for conversion.
3585 // Some of these will be removed altogether.
3586 forAllConstIter(labelHashSet, facesToRemove, fIter)
3587 {
3588 const labelList& fEdges = faceEdges_[fIter.key()];
3590 forAll(fEdges, edgeI)
3591 {
3592 if (whichEdgePatch(fEdges[edgeI]) == patch)
3593 {
3594 // Make an identical map
3595 edgesToConvert.set(fEdges[edgeI], fEdges[edgeI]);
3596 }
3597 else
3598 {
3599 edgesToConvert.set(fEdges[edgeI], -1);
3600 }
3602 if (twoDMesh_)
3603 {
3604 // Add all points as candidates for removal.
3605 // Some (or all) of these will be weeded out.
3606 const edge& eCheck = edges_[fEdges[edgeI]];
3608 pointsToRemove.set(eCheck[0]);
3609 pointsToRemove.set(eCheck[1]);
3610 }
3611 }
3612 }
3614 forAllConstIter(Map<label>, facesToConvert, fIter)
3615 {
3616 const labelList& fEdges = faceEdges_[fIter.key()];
3618 forAll(fEdges, edgeI)
3619 {
3620 if (whichEdgePatch(fEdges[edgeI]) == patch)
3621 {
3622 // Make an identical map
3623 edgesToConvert.set(fEdges[edgeI], fEdges[edgeI]);
3624 }
3625 else
3626 {
3627 edgesToConvert.set(fEdges[edgeI], -1);
3628 }
3629 }
3630 }
3632 // Build a list of edges to be removed.
3633 forAllConstIter(Map<label>, edgesToConvert, eIter)
3634 {
3635 const labelList& eFaces = edgeFaces_[eIter.key()];
3637 bool allRemove = true;
3639 forAll(eFaces, faceI)
3640 {
3641 if (!facesToRemove.found(eFaces[faceI]))
3642 {
3643 allRemove = false;
3644 break;
3645 }
3646 }
3648 if (allRemove)
3649 {
3650 if (!edgesToRemove.found(eIter.key()))
3651 {
3652 edgesToRemove.insert(eIter.key());
3653 }
3654 }
3655 }
3657 // Weed-out the conversion list.
3658 forAllConstIter(labelHashSet, edgesToRemove, eIter)
3659 {
3660 edgesToConvert.erase(eIter.key());
3661 }
3663 // Build a set of points to be removed.
3664 if (twoDMesh_)
3665 {
3666 forAllConstIter(Map<label>, edgesToConvert, eIter)
3667 {
3668 const edge& eCheck = edges_[eIter.key()];
3670 if (pointsToRemove.found(eCheck[0]))
3671 {
3672 pointsToRemove.erase(eCheck[0]);
3673 }
3675 if (pointsToRemove.found(eCheck[1]))
3676 {
3677 pointsToRemove.erase(eCheck[1]);
3678 }
3679 }
3680 }
3681 else
3682 {
3683 forAllConstIter(labelHashSet, edgesToRemove, eIter)
3684 {
3685 const edge& edgeToCheck = edges_[eIter.key()];
3687 forAll(edgeToCheck, pointI)
3688 {
3689 const labelList& pEdges = pointEdges_[edgeToCheck[pointI]];
3691 bool allRemove = true;
3693 forAll(pEdges, edgeI)
3694 {
3695 if (!edgesToRemove.found(pEdges[edgeI]))
3696 {
3697 allRemove = false;
3698 break;
3699 }
3700 }
3702 if (allRemove)
3703 {
3704 if (!pointsToRemove.found(edgeToCheck[pointI]))
3705 {
3706 pointsToRemove.insert(edgeToCheck[pointI]);
3707 }
3708 }
3709 }
3710 }
3711 }
3713 forAllIter(Map<label>, edgesToConvert, eIter)
3714 {
3715 const labelList& eFaces = edgeFaces_[eIter.key()];
3717 label nConvFaces = 0;
3719 forAll(eFaces, faceI)
3720 {
3721 if (facesToConvert.found(eFaces[faceI]))
3722 {
3723 nConvFaces++;
3724 }
3725 }
3727 if (nConvFaces > 2)
3728 {
3729 Pout<< "Invalid conversion. Bailing out." << endl;
3730 return map;
3731 }
3732 }
3734 // Are we only performing checks?
3735 if (checkOnly)
3736 {
3737 // Make necessary map entries
3738 forAllConstIter(Map<label>, facesToConvert, fIter)
3739 {
3740 map.removeFace(fIter.key());
3741 }
3743 forAllConstIter(Map<label>, edgesToConvert, eIter)
3744 {
3745 map.removeEdge(eIter.key());
3746 }
3748 forAllConstIter(labelHashSet, facesToRemove, fIter)
3749 {
3750 map.removeFace(fIter.key());
3751 }
3753 forAllConstIter(labelHashSet, edgesToRemove, eIter)
3754 {
3755 map.removeEdge(eIter.key());
3756 }
3758 forAllConstIter(labelHashSet, pointsToRemove, pIter)
3759 {
3760 map.removePoint(pIter.key());
3761 }
3763 forAll(cList, cellI)
3764 {
3765 map.removeCell(cList[cellI]);
3766 }
3768 // Specify that the operation was successful
3769 map.type() = 1;
3771 return map;
3772 }
3774 // Perform a few debug calls before modifications
3775 if (debug > 1)
3776 {
3777 Pout<< nl << nl
3778 << " Removing cell(s): " << cList
3779 << " and adding internal faces / edges to patch: "
3780 << boundaryMesh()[patch].name()
3781 << endl;
3782 }
3784 // Write out candidates for post-processing
3785 if (debug > 2)
3786 {
3787 writeVTK("pointsToRemove_" + rcN, pointsToRemove.toc(), 0, false, true);
3788 writeVTK("edgesToRemove_" + rcN, edgesToRemove.toc(), 1, false, true);
3789 writeVTK("facesToRemove_" + rcN, facesToRemove.toc(), 2, false, true);
3790 writeVTK("cellsToRemove_" + rcN, cList, 3, false, true);
3791 writeVTK("edgesToConvert_" + rcN, edgesToConvert.toc(), 1, false, true);
3792 writeVTK("facesToConvert_" + rcN, facesToConvert.toc(), 2, false, true);
3793 }
3795 // Loop through all faces for conversion, check orientation
3796 // and create new faces in their place.
3797 forAllIter(Map<label>, facesToConvert, fIter)
3798 {
3799 // Check if this internal face is oriented properly.
3800 face newFace;
3801 label newOwner = -1;
3802 labelList fEdges = faceEdges_[fIter.key()];
3804 if (findIndex(cList, neighbour_[fIter.key()]) != -1)
3805 {
3806 // Orientation is correct
3807 newFace = faces_[fIter.key()];
3808 newOwner = owner_[fIter.key()];
3809 }
3810 else
3811 if (findIndex(cList, owner_[fIter.key()]) != -1)
3812 {
3813 // Face is to be reversed.
3814 newFace = faces_[fIter.key()].reverseFace();
3815 newOwner = neighbour_[fIter.key()];
3817 setFlip(fIter.key());
3818 }
3819 else
3820 {
3821 // Something's terribly wrong
3822 FatalErrorIn
3823 (
3824 "\n"
3825 "const changeMap dynamicTopoFvMesh::removeCells\n"
3826 "(\n"
3827 " const labelList& cList,\n"
3828 " const label patch,\n"
3829 " const word& rcN,\n"
3830 " bool checkOnly\n"
3831 ")\n"
3832 )
3833 << nl << " Invalid mesh. "
3834 << abort(FatalError);
3835 }
3837 // Insert the reconfigured face at the boundary.
3838 fIter() =
3839 (
3840 insertFace
3841 (
3842 patch,
3843 newFace,
3844 newOwner,
3845 -1
3846 )
3847 );
3849 // Add the faceEdges entry.
3850 // Edges will be corrected later.
3851 faceEdges_.append(fEdges);
3853 // Add this face to the map.
3854 map.addFace(fIter());
3856 // Set mapping information
3857 // - But where to map from?
3858 setFaceMapping(fIter());
3860 // Replace cell with the new face label
3861 meshOps::replaceLabel
3862 (
3863 fIter.key(),
3864 fIter(),
3865 cells_[newOwner]
3866 );
3868 // Remove the internal face.
3869 removeFace(fIter.key());
3871 // Update map
3872 map.removeFace(fIter.key());
3873 }
3875 // Create a new edge for each converted edge
3876 forAllIter(Map<label>, edgesToConvert, eIter)
3877 {
3878 if (eIter() == -1)
3879 {
3880 // Create copies before appending.
3881 edge newEdge = edges_[eIter.key()];
3882 labelList eFaces = edgeFaces_[eIter.key()];
3884 eIter() =
3885 (
3886 insertEdge
3887 (
3888 patch,
3889 newEdge,
3890 eFaces
3891 )
3892 );
3894 // Add this edge to the map.
3895 map.addEdge(eIter());
3897 // Remove the edge
3898 removeEdge(eIter.key());
3900 // Update map
3901 map.removeEdge(eIter.key());
3902 }
3903 }
3905 // Loop through all faces for conversion, and replace edgeFaces.
3906 forAllConstIter(Map<label>, facesToConvert, fIter)
3907 {
3908 // Make a copy, because this list is going to
3909 // be modified within this loop.
3910 labelList fEdges = faceEdges_[fIter()];
3912 forAll(fEdges, edgeI)
3913 {
3914 if (edgesToConvert.found(fEdges[edgeI]))
3915 {
3916 meshOps::replaceLabel
3917 (
3918 fIter.key(),
3919 fIter(),
3920 edgeFaces_[edgesToConvert[fEdges[edgeI]]]
3921 );
3923 meshOps::replaceLabel
3924 (
3925 fEdges[edgeI],
3926 edgesToConvert[fEdges[edgeI]],
3927 faceEdges_[fIter()]
3928 );
3929 }
3930 }
3931 }
3933 // Loop through all edges for conversion, and size-down edgeFaces.
3934 forAllConstIter(Map<label>, edgesToConvert, eIter)
3935 {
3936 // Make a copy, because this list is going to
3937 // be modified within this loop.
3938 labelList eFaces = edgeFaces_[eIter()];
3940 forAll(eFaces, faceI)
3941 {
3942 if (facesToRemove.found(eFaces[faceI]))
3943 {
3944 meshOps::sizeDownList
3945 (
3946 eFaces[faceI],
3947 edgeFaces_[eIter()]
3948 );
3949 }
3951 // Replace old edges with new ones.
3952 labelList& fEdges = faceEdges_[eFaces[faceI]];
3954 forAll(fEdges, edgeI)
3955 {
3956 if (edgesToConvert.found(fEdges[edgeI]))
3957 {
3958 fEdges[edgeI] = edgesToConvert[fEdges[edgeI]];
3959 }
3960 }
3961 }
3962 }
3964 // Remove unwanted faces
3965 forAllConstIter(labelHashSet, facesToRemove, fIter)
3966 {
3967 removeFace(fIter.key());
3969 // Update map
3970 map.removeFace(fIter.key());
3971 }
3973 // Remove unwanted edges
3974 forAllConstIter(labelHashSet, edgesToRemove, eIter)
3975 {
3976 removeEdge(eIter.key());
3978 // Update map
3979 map.removeEdge(eIter.key());
3980 }
3982 // Remove unwanted points
3983 forAllConstIter(labelHashSet, pointsToRemove, pIter)
3984 {
3985 removePoint(pIter.key());
3987 // Update map
3988 map.removePoint(pIter.key());
3989 }
3991 // Remove all cells
3992 forAll(cList, cellI)
3993 {
3994 removeCell(cList[cellI]);
3996 // Update map
3997 map.removeCell(cList[cellI]);
3998 }
4000 // Set the flag
4001 topoChangeFlag_ = true;
4003 // Specify that the operation was successful
4004 map.type() = 1;
4006 return map;
4007 }
4010 // Handle topology changes for coupled patches
4011 void dynamicTopoFvMesh::handleCoupledPatches
4012 (
4013 labelHashSet& entities
4014 )
4015 {
4016 // Initialize coupled patch connectivity for topology modifications.
4017 initCoupledConnectivity(this);
4019 if (patchCoupling_.empty() && procIndices_.empty())
4020 {
4021 return;
4022 }
4024 // Move coupled subMeshes
4025 moveCoupledSubMeshes();
4027 // Exchange length-scale buffers across processors.
4028 exchangeLengthBuffers();
4030 if (debug)
4031 {
4032 // Check coupled-patch sizes first.
4033 forAll(patchCoupling_, patchI)
4034 {
4035 if (patchCoupling_(patchI))
4036 {
4037 const coupleMap& cMap = patchCoupling_[patchI].map();
4039 label mSize = patchSizes_[cMap.masterIndex()];
4040 label sSize = patchSizes_[cMap.slaveIndex()];
4042 if (mSize != sSize)
4043 {
4044 Pout<< "Coupled patch-count is inconsistent." << nl
4045 << " Master Patch: " << cMap.masterIndex()
4046 << " Count: " << mSize << nl
4047 << " Slave Patch: " << cMap.slaveIndex()
4048 << " Count: " << sSize
4049 << endl;
4051 FatalErrorIn
4052 (
4053 "void dynamicTopoFvMesh::handleCoupledPatches()"
4054 )
4055 << " Failures were found in connectivity"
4056 << " prior to coupled topo-changes."
4057 << abort(FatalError);
4058 }
4059 }
4060 }
4062 Info<< "Handling coupled patches...";
4063 }
4065 if (Pstream::parRun())
4066 {
4067 // Calculate mapping offsets
4068 const polyBoundaryMesh& boundary = boundaryMesh();
4070 // Determine number of physical (non-processor) patches
4071 label nPhysical = 0;
4073 forAll(boundary, patchI)
4074 {
4075 if (isA<processorPolyPatch>(boundary[patchI]))
4076 {
4077 continue;
4078 }
4080 nPhysical++;
4081 }
4083 // Fetch number of processors
4084 label nProcs = procIndices_.size();
4086 // Allocate cell, face and patch offsets
4087 labelList cellSizes(nProcs, 0), cellStarts(nProcs, 0);
4088 labelList faceSizes(nProcs, 0), faceStarts(nProcs, 0);
4089 labelListList patchSizes(nProcs, labelList(nPhysical, 0));
4090 labelListList patchStarts(nProcs, labelList(nPhysical, 0));
4092 label nTotalCells = nOldCells_, nTotalIntFaces = nOldInternalFaces_;
4093 labelList nTotalPatchFaces(SubList<label>(oldPatchSizes_, nPhysical));
4095 forAll(procIndices_, pI)
4096 {
4097 const coupleMap& cMap = recvMeshes_[pI].map();
4099 // Fetch size from subMesh
4100 label nCells = cMap.nEntities(coupleMap::CELL);
4101 label nIntFaces = cMap.nEntities(coupleMap::INTERNAL_FACE);
4103 // Set size / offset for this processor
4104 cellSizes[pI] = nCells;
4105 cellStarts[pI] = nTotalCells;
4107 faceSizes[pI] = nIntFaces;
4108 faceStarts[pI] = nTotalIntFaces;
4110 // Update count
4111 nTotalCells += nCells;
4112 nTotalIntFaces += nIntFaces;
4114 // Fetch patch sizes from subMesh
4115 const labelList& nPatchFaces =
4116 (
4117 cMap.entityBuffer(coupleMap::FACE_SIZES)
4118 );
4120 // Loop over physical patches
4121 forAll(nTotalPatchFaces, patchI)
4122 {
4123 // Fetch patch size from subMesh
4124 label nFaces = nPatchFaces[patchI];
4126 // Set patch size / offset for this processor
4127 patchSizes[pI][patchI] = nFaces;
4128 patchStarts[pI][patchI] = nTotalPatchFaces[patchI];
4130 // Update patch count
4131 nTotalPatchFaces[patchI] += nFaces;
4132 }
4133 }
4135 // Set sizes / starts in mapper
4136 mapper_->setOffsets
4137 (
4138 cellSizes,
4139 cellStarts,
4140 faceSizes,
4141 faceStarts,
4142 patchSizes,
4143 patchStarts
4144 );
4146 if (debug > 3)
4147 {
4148 SubList<label> physicalPatches(oldPatchSizes_, nPhysical);
4150 Pout<< " procIndices: " << procIndices_ << nl
4151 << " nCells: " << nOldCells_ << nl
4152 << " proc cellSizes: " << cellSizes << nl
4153 << " cellStarts: " << cellStarts << nl
4154 << " proc faceSizes: " << faceSizes << nl
4155 << " faceStarts: " << faceStarts << nl
4156 << " patchSizes: " << physicalPatches << nl
4157 << " proc patchSizes: " << patchSizes << nl
4158 << " patchStarts: " << patchStarts << endl;
4159 }
4160 }
4162 // Optionally switch off topo-changes
4163 // on processor patches at run-time
4164 bool procTopoChanges = true;
4166 const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
4168 if (meshSubDict.found("procTopoChanges") || mandatory_)
4169 {
4170 procTopoChanges = readBool(meshSubDict.lookup("procTopoChanges"));
4171 }
4173 // Set coupled modifications.
4174 setCoupledModification();
4176 // Loop through the coupled stack and perform changes.
4177 if (edgeRefinement_)
4178 {
4179 // Initialize the coupled stack
4180 if (procTopoChanges)
4181 {
4182 initCoupledStack(entities, false);
4183 }
4185 edgeRefinementEngine(&(handlerPtr_[0]));
4186 }
4188 // Re-Initialize the stack
4189 if (procTopoChanges)
4190 {
4191 initCoupledStack(entities, false);
4192 }
4194 if (twoDMesh_)
4195 {
4196 swap2DEdges(&(handlerPtr_[0]));
4197 }
4198 else
4199 {
4200 swap3DEdges(&(handlerPtr_[0]));
4201 }
4203 // Build a list of entities that need to be avoided
4204 // by regular topo-changes.
4205 buildEntitiesToAvoid(entities, true);
4207 // Reset coupled modifications.
4208 unsetCoupledModification();
4210 if (debug)
4211 {
4212 Info<< "Done." << endl;
4214 // Check coupled-patch sizes after changes.
4215 forAll(patchCoupling_, patchI)
4216 {
4217 if (patchCoupling_(patchI))
4218 {
4219 const coupleMap& cMap = patchCoupling_[patchI].map();
4221 label mSize = patchSizes_[cMap.masterIndex()];
4222 label sSize = patchSizes_[cMap.slaveIndex()];
4224 if (mSize != sSize)
4225 {
4226 Pout<< "Coupled patch-count is inconsistent." << nl
4227 << " Master Patch: " << cMap.masterIndex()
4228 << " Count: " << mSize << nl
4229 << " Slave Patch: " << cMap.slaveIndex()
4230 << " Count: " << sSize
4231 << endl;
4233 FatalErrorIn("dynamicTopoFvMesh::handleCoupledPatches()")
4234 << " Failures were found in connectivity"
4235 << " after coupled topo-changes."
4236 << abort(FatalError);
4237 }
4238 }
4239 }
4240 }
4242 // Schedule transfer of topology operations across processors
4243 forAll(procIndices_, pI)
4244 {
4245 label proc = procIndices_[pI];
4247 coupledInfo& sPM = sendMeshes_[pI];
4248 coupledInfo& rPM = recvMeshes_[pI];
4250 if (proc < Pstream::myProcNo())
4251 {
4252 const coupleMap& cMap = sPM.map();
4254 // How many entities am I receiving..
4255 label nEntities = -1, nMovePoints = -1;
4257 meshOps::pRead(proc, nEntities);
4258 meshOps::pRead(proc, nMovePoints);
4260 if (debug > 3)
4261 {
4262 Pout<< " Op tranfer:"
4263 << " Receiving from [" << proc << "]::"
4264 << " nEntities: " << nEntities
4265 << " nMovePoints: " << nMovePoints
4266 << endl;
4267 }
4269 if (nEntities)
4270 {
4271 // Size up the receipt buffers
4272 labelList& indices = cMap.entityIndices();
4273 List<coupleMap::opType>& operations = cMap.entityOperations();
4275 indices.setSize(nEntities);
4276 operations.setSize(nEntities);
4278 // Schedule indices and operations for receipt
4279 meshOps::pRead(proc, indices);
4280 meshOps::pRead(proc, operations);
4281 }
4283 if (nMovePoints)
4284 {
4285 // Size up the receipt buffers
4286 pointField& newPoints = cMap.moveNewPoints();
4287 pointField& oldPoints = cMap.moveOldPoints();
4289 newPoints.setSize(nMovePoints);
4290 oldPoints.setSize(nMovePoints);
4292 // Schedule old / new points for receipt
4293 meshOps::pRead(proc, newPoints);
4294 meshOps::pRead(proc, oldPoints);
4295 }
4296 }
4297 else
4298 {
4299 const coupleMap& cMap = rPM.map();
4301 label nEntities = cMap.entityIndices().size();
4302 label nMovePoints = cMap.moveNewPoints().size();
4304 if (debug > 3)
4305 {
4306 Pout<< " Op tranfer:"
4307 << " Sending to [" << proc << "]:: "
4308 << " nEntities: " << nEntities << nl
4309 << " entityIndices: " << cMap.entityIndices() << nl
4310 << " entityOperations: " << cMap.entityOperations() << nl
4311 << " nMovePoints: " << nMovePoints << nl
4312 << " moveNewPoints: " << cMap.moveNewPoints() << nl
4313 << " moveOldPoints: " << cMap.moveOldPoints() << nl
4314 << endl;
4315 }
4317 meshOps::pWrite(proc, nEntities);
4318 meshOps::pWrite(proc, nMovePoints);
4320 if (nEntities)
4321 {
4322 // Schedule transfer to processor
4323 meshOps::pWrite(proc, cMap.entityIndices());
4324 meshOps::pWrite(proc, cMap.entityOperations());
4325 }
4327 if (nMovePoints)
4328 {
4329 // Schedule transfer to processor
4330 meshOps::pWrite(proc, cMap.moveNewPoints());
4331 meshOps::pWrite(proc, cMap.moveOldPoints());
4332 }
4333 }
4334 }
4336 // We won't wait for transfers to complete for the moment,
4337 // and will deal with operations once the internal mesh
4338 // has been dealt with.
4339 }
4342 // Synchronize topology operations across processors
4343 void dynamicTopoFvMesh::syncCoupledPatches(labelHashSet& entities)
4344 {
4345 if (!Pstream::parRun())
4346 {
4347 return;
4348 }
4350 const polyBoundaryMesh& boundary = boundaryMesh();
4352 labelList createPatchOrders(Pstream::nProcs(), 0);
4354 forAll(procIndices_, pI)
4355 {
4356 const coupleMap& cMap = sendMeshes_[pI].map();
4358 if (cMap.patchIndex() >= boundary.size())
4359 {
4360 // This processor needs a new patch talking to me
4361 createPatchOrders[procIndices_[pI]] = 1;
4362 }
4363 }
4365 // Send / recv create patch orders, if any
4366 for (label proc = 0; proc < Pstream::nProcs(); proc++)
4367 {
4368 if (proc == Pstream::myProcNo())
4369 {
4370 continue;
4371 }
4373 if (proc > Pstream::myProcNo())
4374 {
4375 meshOps::pWrite(proc, createPatchOrders[proc]);
4376 }
4377 else
4378 {
4379 meshOps::pRead(proc, createPatchOrders[proc]);
4380 }
4381 }
4383 Map<label> addedProcPatches;
4385 for (label proc = 0; proc < Pstream::myProcNo(); proc++)
4386 {
4387 if (createPatchOrders[proc])
4388 {
4389 // Create a new processor patch
4390 addedProcPatches.insert
4391 (
4392 proc,
4393 createProcessorPatch(proc)
4394 );
4395 }
4396 }
4398 // Wait for all transfers to complete.
4399 meshOps::waitForBuffers();
4401 // Buffer for cell-removal
4402 DynamicList<label> rCellList(10);
4404 forAll(procIndices_, pI)
4405 {
4406 label proc = procIndices_[pI];
4408 const coupledInfo& sPM = sendMeshes_[pI];
4410 if (proc < Pstream::myProcNo())
4411 {
4412 const coupleMap& cMap = sPM.map();
4413 const labelList& indices = cMap.entityIndices();
4414 const List<coupleMap::opType>& operations = cMap.entityOperations();
4415 const pointField& newPoints = cMap.moveNewPoints();
4416 const pointField& oldPoints = cMap.moveOldPoints();
4418 // Find the appropriate processor patch
4419 // for cell-removal operations
4420 label procPatch = -1;
4422 forAll(boundary, pI)
4423 {
4424 if (!isA<processorPolyPatch>(boundary[pI]))
4425 {
4426 continue;
4427 }
4429 const processorPolyPatch& pp =
4430 (
4431 refCast<const processorPolyPatch>(boundary[pI])
4432 );
4434 if (pp.neighbProcNo() == proc)
4435 {
4436 procPatch = pI;
4437 break;
4438 }
4439 }
4441 if (procPatch == -1)
4442 {
4443 // Check if this was a newly added patch
4444 if (addedProcPatches.found(proc))
4445 {
4446 procPatch = addedProcPatches[proc];
4447 }
4448 else
4449 if (operations.size())
4450 {
4451 // If we actually have operations from
4452 // this processor, this is a problem.
4453 Pout<< " * * * Sync Operations * * * " << nl
4454 << " Could not find patch for proc: " << proc << nl
4455 << abort(FatalError);
4456 }
4457 else
4458 {
4459 // Move on to the next processor
4460 continue;
4461 }
4462 }
4464 label pointCounter = 0;
4466 // Sequentially execute operations
4467 forAll(indices, indexI)
4468 {
4469 const label index = indices[indexI];
4470 const coupleMap::opType op = operations[indexI];
4472 if (debug > 3)
4473 {
4474 Pout<< " Recv Op index: " << index << endl;
4475 }
4477 // Determine the appropriate local index
4478 label localIndex = -1;
4480 if (op == coupleMap::MOVE_POINT)
4481 {
4482 localIndex = cMap.entityMap(coupleMap::POINT)[index];
4483 }
4484 else
4485 if (op == coupleMap::CONVERT_PATCH)
4486 {
4487 localIndex =
4488 (
4489 cMap.entityMap(coupleMap::FACE).found(index) ?
4490 cMap.entityMap(coupleMap::FACE)[index] : -1
4491 );
4493 if (debug > 3)
4494 {
4495 const point& newCentre = newPoints[pointCounter];
4496 const point& oldCentre = oldPoints[pointCounter];
4498 Pout<< " Convert patch: " << index << nl
4499 << " localIndex: " << localIndex << nl
4500 << " newCentre: " << newCentre << nl
4501 << " oldCentre: " << oldCentre << nl
4502 << endl;
4503 }
4504 }
4505 else
4506 {
4507 // Pick localIndex based on entity type
4508 if (twoDMesh_)
4509 {
4510 const Map<label>& fM = cMap.entityMap(coupleMap::FACE);
4512 Map<label>::const_iterator it = fM.find(index);
4514 if (it != fM.end())
4515 {
4516 localIndex = it();
4517 }
4518 else
4519 {
4520 Pout<< " * * * Sync Operations * * * " << nl
4521 << " Could not find index: " << index
4522 << " in faceMap for proc: " << proc << nl
4523 << " nSentFaces: "
4524 << cMap.nEntities(coupleMap::FACE)
4525 << " opType: " << coupleMap::asText(op)
4526 << abort(FatalError);
4527 }
4528 }
4529 else
4530 {
4531 const Map<label>& eM = cMap.entityMap(coupleMap::EDGE);
4533 Map<label>::const_iterator it = eM.find(index);
4535 if (it != eM.end())
4536 {
4537 localIndex = it();
4538 }
4539 else
4540 {
4541 Pout<< " * * * Sync Operations * * * " << nl
4542 << " Could not find index: " << index
4543 << " in edgeMap for proc: " << proc << nl
4544 << " nSentEdges: "
4545 << cMap.nEntities(coupleMap::EDGE)
4546 << " opType: " << coupleMap::asText(op)
4547 << abort(FatalError);
4548 }
4549 }
4550 }
4552 changeMap opMap;
4554 switch (op)
4555 {
4556 case coupleMap::BISECTION:
4557 {
4558 opMap = bisectEdge(localIndex);
4559 break;
4560 }
4562 case coupleMap::COLLAPSE_FIRST:
4563 {
4564 opMap = collapseEdge(localIndex, 1);
4565 break;
4566 }
4568 case coupleMap::COLLAPSE_SECOND:
4569 {
4570 opMap = collapseEdge(localIndex, 2);
4571 break;
4572 }
4574 case coupleMap::COLLAPSE_MIDPOINT:
4575 {
4576 opMap = collapseEdge(localIndex, 3);
4577 break;
4578 }
4580 case coupleMap::REMOVE_CELL:
4581 {
4582 // Clear existing list
4583 rCellList.clear();
4585 if (twoDMesh_)
4586 {
4587 // Insert the owner cell
4588 rCellList.append(owner_[localIndex]);
4589 }
4590 else
4591 {
4592 // Insert all cells connected to this edge
4593 const labelList& eFaces = edgeFaces_[localIndex];
4595 forAll(eFaces, faceI)
4596 {
4597 const label own = owner_[eFaces[faceI]];
4598 const label nei = neighbour_[eFaces[faceI]];
4600 if (findIndex(rCellList, own) == -1)
4601 {
4602 rCellList.append(own);
4603 }
4605 if (nei != -1)
4606 {
4607 if (findIndex(rCellList, nei) == -1)
4608 {
4609 rCellList.append(nei);
4610 }
4611 }
4612 }
4613 }
4615 opMap =
4616 (
4617 removeCells
4618 (
4619 rCellList,
4620 procPatch,
4621 "rcs_" + Foam::name(localIndex)
4622 )
4623 );
4625 break;
4626 }
4628 case coupleMap::MOVE_POINT:
4629 {
4630 const point& newPoint = newPoints[pointCounter];
4631 const point& oldPoint = oldPoints[pointCounter];
4633 // Move old / new points
4634 points_[localIndex] = newPoint;
4635 oldPoints_[localIndex] = oldPoint;
4637 // Clear the existing map
4638 opMap.clear();
4640 // Force a successful operation
4641 opMap.type() = 1;
4643 pointCounter++;
4645 break;
4646 }
4648 case coupleMap::CONVERT_PATCH:
4649 {
4650 const point& newCentre = newPoints[pointCounter];
4651 const point& oldCentre = oldPoints[pointCounter];
4653 if (localIndex == -1)
4654 {
4655 // Accumulate stats in case of failure
4656 scalar minDist = GREAT;
4657 vector minPoint = vector::zero;
4658 DynamicList<label> checkedFaces;
4660 if (debug > 2)
4661 {
4662 // Reserve for append
4663 checkedFaces.setCapacity(50);
4664 }
4666 // New face. Check all boundary faces
4667 // and match up centre
4668 label sTot = faces_.size();
4669 label sInt = nOldInternalFaces_;
4671 for (label faceI = sInt; faceI < sTot; faceI++)
4672 {
4673 const face& fCheck = faces_[faceI];
4675 if (fCheck.empty())
4676 {
4677 continue;
4678 }
4680 label pIndex = whichPatch(faceI);
4682 if (getNeighbourProcessor(pIndex) == -1)
4683 {
4684 continue;
4685 }
4687 // Compute face-centre
4688 vector fC = fCheck.centre(points_);
4690 // Compute tolerance
4691 scalar tol = mag(points_[fCheck[0]] - fC);
4692 scalar dist = mag(fC - newCentre);
4694 if (dist < (1e-4 * tol))
4695 {
4696 localIndex = faceI;
4697 break;
4698 }
4699 else
4700 if (dist < minDist)
4701 {
4702 minPoint = fC;
4703 minDist = dist;
4705 if (debug > 2)
4706 {
4707 checkedFaces.append(faceI);
4708 }
4709 }
4710 }
4712 // Ensure that the face was found
4713 if (localIndex == -1)
4714 {
4715 writeVTK
4716 (
4717 "checkedFaces_"
4718 + Foam::name(index),
4719 checkedFaces,
4720 2, false, true
4721 );
4723 Pout<< " * * * syncCoupledPatches * * * " << nl
4724 << " Convert patch Op failed." << nl
4725 << " Face: " << index << nl
4726 << " minPoint: " << minPoint << nl
4727 << " minDistance: " << minDist << nl
4728 << " newCentre: " << newCentre << nl
4729 << " oldCentre: " << oldCentre << nl
4730 << abort(FatalError);
4731 }
4732 }
4734 // Fetch reference to face
4735 const face& fCheck = faces_[localIndex];
4737 point fC = fCheck.centre(points_);
4739 // Specify a tolerance
4740 scalar tol = mag(points_[fCheck[0]] - fC);
4741 scalar dist = mag(fC - newCentre);
4743 // Ensure a face-match
4744 if (dist > (1e-4 * tol))
4745 {
4746 Pout<< " * * * Sync Operations * * * " << nl
4747 << " Convert patch Op failed." << nl
4748 << " Index: " << index << nl
4749 << " localIndex: " << localIndex << nl
4750 << " face: " << fCheck << nl
4751 << " faceCentre: " << fC << nl
4752 << " Master processor: " << proc << nl
4753 << " procPatch: " << procPatch << nl
4754 << " tolerance: " << tol << nl
4755 << " distance: " << dist << nl
4756 << " pointCounter: " << pointCounter << nl
4757 << " newCentre: " << newCentre << nl
4758 << " oldCentre: " << oldCentre << nl
4759 << abort(FatalError);
4760 }
4762 // Obtain a copy before adding the new face,
4763 // since the reference might become
4764 // invalid during list resizing.
4765 face newFace = faces_[localIndex];
4766 label newOwn = owner_[localIndex];
4767 labelList newFaceEdges = faceEdges_[localIndex];
4769 label newFaceIndex =
4770 (
4771 insertFace
4772 (
4773 procPatch,
4774 newFace,
4775 newOwn,
4776 -1
4777 )
4778 );
4780 // Add a faceEdges entry as well.
4781 // Edges don't have to change, since they're
4782 // all on the boundary anyway.
4783 faceEdges_.append(newFaceEdges);
4785 meshOps::replaceLabel
4786 (
4787 localIndex,
4788 newFaceIndex,
4789 cells_[newOwn]
4790 );
4792 // Correct edgeFaces with the new face label.
4793 forAll(newFaceEdges, edgeI)
4794 {
4795 meshOps::replaceLabel
4796 (
4797 localIndex,
4798 newFaceIndex,
4799 edgeFaces_[newFaceEdges[edgeI]]
4800 );
4801 }
4803 // Finally remove the face
4804 removeFace(localIndex);
4806 // Clear the existing map
4807 opMap.clear();
4809 // Update map
4810 opMap.addFace
4811 (
4812 newFaceIndex,
4813 labelList(1, localIndex)
4814 );
4816 // Force a successful operation
4817 opMap.type() = 1;
4819 pointCounter++;
4821 break;
4822 }
4824 case coupleMap::INVALID:
4825 {
4826 Pout<< " * * * Sync Operations * * * " << nl
4827 << " Invalid operation." << nl
4828 << " Index: " << index << nl
4829 << " localIndex: " << localIndex << nl
4830 << " operation: " << coupleMap::asText(op) << nl
4831 << abort(FatalError);
4833 break;
4834 }
4835 }
4837 if (opMap.type() <= 0)
4838 {
4839 Pout<< " * * * Sync Operations * * * " << nl
4840 << " Operation failed." << nl
4841 << " Index: " << index << nl
4842 << " localIndex: " << localIndex << nl
4843 << " operation: " << coupleMap::asText(op) << nl
4844 << " opMap.type: " << opMap.type() << nl
4845 << abort(FatalError);
4846 }
4847 }
4848 }
4849 }
4851 // Re-Initialize the stack with avoided entities from subMeshes
4852 // and leave those on processor patches untouched
4853 labelHashSet procEntities;
4855 buildEntitiesToAvoid(procEntities, false);
4857 // First remove processor entries
4858 forAllConstIter(labelHashSet, procEntities, pIter)
4859 {
4860 if (entities.found(pIter.key()))
4861 {
4862 entities.erase(pIter.key());
4863 }
4864 }
4866 // Initialize the coupled stack, using supplied entities
4867 initCoupledStack(entities, true);
4869 // Loop through the coupled stack and perform changes.
4870 if (edgeRefinement_)
4871 {
4872 edgeRefinementEngine(&(handlerPtr_[0]));
4873 }
4875 // Re-Initialize the stack, using supplied entities
4876 initCoupledStack(entities, true);
4878 if (twoDMesh_)
4879 {
4880 swap2DEdges(&(handlerPtr_[0]));
4881 }
4882 else
4883 {
4884 swap3DEdges(&(handlerPtr_[0]));
4885 }
4886 }
4889 // Check the state of coupled boundaries
4890 // - Return true if errors are present
4891 bool dynamicTopoFvMesh::checkCoupledBoundaries(bool report) const
4892 {
4893 const polyBoundaryMesh& boundary = boundaryMesh();
4895 bool sizeError = false, misMatchError = false;
4897 // Check if a geometric tolerance has been specified.
4898 scalar gTol = debug::tolerances("processorMatchTol", 1e-4);
4900 // Maintain a list of master / neighbour anchors
4901 List<vectorField> mAnchors(boundary.size());
4902 List<vectorField> nAnchors(boundary.size());
4904 forAll(boundary, pI)
4905 {
4906 if (isA<processorPolyPatch>(boundary[pI]))
4907 {
4908 const processorPolyPatch& pp =
4909 (
4910 refCast<const processorPolyPatch>(boundary[pI])
4911 );
4913 label neiProcNo = pp.neighbProcNo();
4915 // Send point / face sizes
4916 meshOps::pWrite(neiProcNo, boundary[pI].nPoints());
4917 meshOps::pWrite(neiProcNo, boundary[pI].size());
4919 // Send points / centres / areas
4920 meshOps::pWrite(neiProcNo, boundary[pI].faceAreas());
4921 meshOps::pWrite(neiProcNo, boundary[pI].faceCentres());
4922 meshOps::pWrite(neiProcNo, boundary[pI].localPoints());
4924 // Prepare and send anchor points
4925 mAnchors[pI].setSize(boundary[pI].size());
4927 const faceList& lFaces = boundary[pI].localFaces();
4928 const pointField& lPoints = boundary[pI].localPoints();
4930 forAll(lFaces, faceI)
4931 {
4932 mAnchors[pI][faceI] = lPoints[lFaces[faceI][0]];
4933 }
4935 meshOps::pWrite(neiProcNo, mAnchors[pI]);
4936 }
4938 if (isA<cyclicPolyPatch>(boundary[pI]))
4939 {
4940 const cyclicPolyPatch& cp =
4941 (
4942 refCast<const cyclicPolyPatch>(boundary[pI])
4943 );
4945 if (cp.size() & 1)
4946 {
4947 sizeError = true;
4949 Pout<< " Incorrect size for cyclic patch: "
4950 << cp.size() << endl;
4951 }
4953 if (!sizeError)
4954 {
4955 // Compute halfSize
4956 label halfSize = cp.size() / 2;
4958 vectorField half0Areas
4959 (
4960 SubField<vector>
4961 (
4962 cp.faceAreas(),
4963 halfSize
4964 )
4965 );
4967 vectorField half0Centres
4968 (
4969 SubField<vector>
4970 (
4971 cp.faceCentres(),
4972 halfSize
4973 )
4974 );
4976 vectorField half1Areas
4977 (
4978 SubField<vector>
4979 (
4980 cp.faceAreas(),
4981 halfSize,
4982 halfSize
4983 )
4984 );
4986 vectorField half1Centres
4987 (
4988 SubField<vector>
4989 (
4990 cp.faceCentres(),
4991 halfSize,
4992 halfSize
4993 )
4994 );
4996 const faceList& lF = boundary[pI].localFaces();
4997 const pointField& lP = boundary[pI].localPoints();
4998 const vectorField& lC = boundary[pI].faceCentres();
5000 // Prepare anchor points
5001 mAnchors[pI].setSize(boundary[pI].size());
5003 forAll(lF, faceI)
5004 {
5005 mAnchors[pI][faceI] = lP[lF[faceI][0]];
5006 }
5008 // Transform first-half points and check
5009 if (cp.transform() == cyclicPolyPatch::ROTATIONAL)
5010 {
5011 forAll(half0Centres, faceI)
5012 {
5013 half0Centres[faceI] =
5014 (
5015 Foam::transform
5016 (
5017 cp.transformT(0),
5018 half0Centres[faceI]
5019 )
5020 );
5022 mAnchors[pI][faceI] =
5023 (
5024 Foam::transform
5025 (
5026 cp.transformT(0),
5027 mAnchors[pI][faceI]
5028 )
5029 );
5030 }
5031 }
5032 else
5033 if (cp.transform() == cyclicPolyPatch::TRANSLATIONAL)
5034 {
5035 forAll(half0Centres, faceI)
5036 {
5037 half0Centres[faceI] += cp.separationVector();
5038 mAnchors[pI][faceI] += cp.separationVector();
5039 }
5040 }
5041 else
5042 {
5043 Pout<< " Cyclic check: Unknown transform."
5044 << abort(FatalError);
5045 }
5047 // Calculate a point-match tolerance per patch
5048 scalar pTol = -GREAT;
5050 // Check areas / compute tolerance
5051 forAll(half0Areas, faceI)
5052 {
5053 scalar fMagSf = mag(half0Areas[faceI]);
5054 scalar rMagSf = mag(half1Areas[faceI]);
5055 scalar avSf = 0.5 * (fMagSf + rMagSf);
5057 if (mag(fMagSf - rMagSf)/avSf > gTol)
5058 {
5059 misMatchError = true;
5061 Pout<< " Face: " << faceI
5062 << " area does not match neighbour by: "
5063 << 100 * mag(fMagSf - rMagSf)/avSf
5064 << "% - possible patch ordering problem. "
5065 << " Front area:" << fMagSf
5066 << " Rear area: " << rMagSf
5067 << endl;
5068 }
5070 pTol =
5071 (
5072 Foam::max
5073 (
5074 pTol,
5075 gTol * mag(lP[lF[faceI][0]] - lC[faceI])
5076 )
5077 );
5078 }
5080 // Check centres / anchor points
5081 forAll(half0Centres, faceI)
5082 {
5083 scalar distA =
5084 (
5085 mag
5086 (
5087 mAnchors[pI][faceI]
5088 - mAnchors[pI][faceI + halfSize]
5089 )
5090 );
5092 scalar distC =
5093 (
5094 mag
5095 (
5096 half0Centres[faceI]
5097 - half1Centres[faceI]
5098 )
5099 );
5101 if (distA > pTol || distC > pTol)
5102 {
5103 misMatchError = true;
5105 UIndirectList<point> f1(lP, lF[faceI]);
5106 UIndirectList<point> f2(lP, lF[faceI + halfSize]);
5108 Pout<< " Face: " << faceI << nl
5109 << " Points: " << nl << f1 << nl << f2 << nl
5110 << " Anchors ::" << nl
5111 << mAnchors[pI][faceI] << nl
5112 << mAnchors[pI][faceI + halfSize] << nl
5113 << " Centres ::" << nl
5114 << half0Centres[faceI] << nl
5115 << half1Centres[faceI] << nl
5116 << " Tolerance: " << pTol << nl
5117 << " Measured Anchor distance: " << distA << nl
5118 << " Measured Centre distance: " << distC << nl
5119 << endl;
5120 }
5121 }
5123 if (misMatchError)
5124 {
5125 // Write out to disk
5126 meshOps::writeVTK
5127 (
5128 (*this),
5129 cp.name(),
5130 identity(cp.size()),
5131 2,
5132 cp.localPoints(),
5133 List<edge>(0),
5134 cp.localFaces()
5135 );
5136 }
5137 }
5138 }
5139 }
5141 // Maintain a list of points / areas / centres
5142 List<vectorField> fAreas(boundary.size());
5143 List<vectorField> fCentres(boundary.size());
5144 List<vectorField> neiPoints(boundary.size());
5146 forAll(boundary, pI)
5147 {
5148 if (isA<processorPolyPatch>(boundary[pI]))
5149 {
5150 const processorPolyPatch& pp =
5151 (
5152 refCast<const processorPolyPatch>(boundary[pI])
5153 );
5155 label neiProcNo = pp.neighbProcNo();
5157 // Receive point / face sizes
5158 label neiPSize = -1, neiSize = -1;
5160 meshOps::pRead(neiProcNo, neiPSize);
5161 meshOps::pRead(neiProcNo, neiSize);
5163 // Size up lists
5164 fAreas[pI].setSize(neiSize);
5165 fCentres[pI].setSize(neiSize);
5166 neiPoints[pI].setSize(neiPSize);
5167 nAnchors[pI].setSize(neiSize);
5169 // Receive points / centres / areas / anchors
5170 meshOps::pRead(neiProcNo, fAreas[pI]);
5171 meshOps::pRead(neiProcNo, fCentres[pI]);
5172 meshOps::pRead(neiProcNo, neiPoints[pI]);
5173 meshOps::pRead(neiProcNo, nAnchors[pI]);
5175 if
5176 (
5177 neiSize != boundary[pI].size() ||
5178 neiPSize != boundary[pI].nPoints()
5179 )
5180 {
5181 sizeError = true;
5183 Pout<< "Incorrect send / recv sizes: " << nl
5184 << " Patch: " << boundary[pI].name() << nl
5185 << " My nPoints: " << boundary[pI].nPoints() << nl
5186 << " Nei nPoints: " << neiPSize << nl
5187 << " My nFaces: " << boundary[pI].size() << nl
5188 << " Nei nFaces: " << neiSize << nl
5189 << endl;
5190 }
5191 }
5192 }
5194 // Wait for transfers to complete
5195 meshOps::waitForBuffers();
5197 // Check centres / areas
5198 forAll(boundary, pI)
5199 {
5200 if (!isA<processorPolyPatch>(boundary[pI]))
5201 {
5202 continue;
5203 }
5205 const vectorField& myAreas = boundary[pI].faceAreas();
5206 const vectorField& myCentres = boundary[pI].faceCentres();
5208 // Fetch local connectivity
5209 const faceList& myFaces = boundary[pI].localFaces();
5210 const vectorField& myPoints = boundary[pI].localPoints();
5212 // Calculate a point-match tolerance per patch
5213 scalar pTol = -GREAT;
5215 forAll(myAreas, faceI)
5216 {
5217 scalar magSf = mag(myAreas[faceI]);
5218 scalar nbrMagSf = mag(fAreas[pI][faceI]);
5219 scalar avSf = 0.5 * (magSf + nbrMagSf);
5221 if (mag(magSf - nbrMagSf)/avSf > gTol)
5222 {
5223 misMatchError = true;
5225 Pout<< " Face: " << faceI
5226 << " area does not match neighbour by: "
5227 << 100 * mag(magSf - nbrMagSf)/avSf
5228 << "% - possible patch ordering problem. "
5229 << " My area:" << magSf
5230 << " Neighbour area: " << nbrMagSf
5231 << " My centre: " << myCentres[faceI]
5232 << endl;
5233 }
5235 pTol =
5236 (
5237 Foam::max
5238 (
5239 pTol,
5240 gTol *
5241 mag
5242 (
5243 myPoints[myFaces[faceI][0]]
5244 - myCentres[faceI]
5245 )
5246 )
5247 );
5248 }
5250 const processorPolyPatch& pp =
5251 (
5252 refCast<const processorPolyPatch>(boundary[pI])
5253 );
5255 // Check anchor points
5256 forAll(mAnchors[pI], faceI)
5257 {
5258 scalar dist = mag(mAnchors[pI][faceI] - nAnchors[pI][faceI]);
5260 if (dist > pTol)
5261 {
5262 misMatchError = true;
5264 Pout<< " Face: " << faceI << "::" << mAnchors[pI][faceI] << nl
5265 << " Mismatch for processor: " << pp.neighbProcNo() << nl
5266 << " Neighbour Anchor: " << nAnchors[pI][faceI] << nl
5267 << " Tolerance: " << pTol << nl
5268 << " Measured distance: " << dist << nl
5269 << endl;
5270 }
5271 }
5273 // Check neighbPoints addressing
5274 const labelList& nPts = pp.neighbPoints();
5276 forAll(myPoints, pointI)
5277 {
5278 if (nPts[pointI] < 0)
5279 {
5280 Pout<< " Point: " << pointI << "::" << myPoints[pointI]
5281 << " reported to be multiply connected." << nl
5282 << " neighbPoint: " << nPts[pointI] << nl
5283 << endl;
5285 misMatchError = true;
5287 continue;
5288 }
5290 scalar dist = mag(myPoints[pointI] - neiPoints[pI][nPts[pointI]]);
5292 if (dist > pTol)
5293 {
5294 misMatchError = true;
5296 Pout<< " Point: " << pointI << "::" << myPoints[pointI] << nl
5297 << " Mismatch for processor: " << pp.neighbProcNo() << nl
5298 << " Neighbour Pt: " << neiPoints[pI][nPts[pointI]] << nl
5299 << " Tolerance: " << pTol << nl
5300 << " Measured distance: " << dist << nl
5301 << " Neighbour Index:" << nPts[pointI] << nl
5302 << endl;
5303 }
5304 }
5305 }
5307 if (!sizeError && !misMatchError && report)
5308 {
5309 Pout<< " Coupled boundary check: No problems." << endl;
5310 }
5312 return (sizeError || misMatchError);
5313 }
5316 // Build patch sub-meshes for processor patches
5317 void dynamicTopoFvMesh::buildProcessorPatchMeshes()
5318 {
5319 if (procIndices_.empty())
5320 {
5321 return;
5322 }
5324 // Maintain a list of cells common to multiple processors.
5325 Map<labelList> commonCells;
5327 forAll(procIndices_, pI)
5328 {
5329 label proc = procIndices_[pI];
5331 coupledInfo& sPM = sendMeshes_[pI];
5333 // Build the subMesh.
5334 buildProcessorPatchMesh(sPM, commonCells);
5336 const coupleMap& scMap = sPM.map();
5338 // Send my sub-mesh to the neighbour.
5339 meshOps::pWrite(proc, scMap.nEntities());
5341 if (debug > 3)
5342 {
5343 Pout<< "Sending to [" << proc << "]:: nEntities: "
5344 << scMap.nEntities() << endl;
5345 }
5347 // Send the pointBuffers
5348 meshOps::pWrite(proc, scMap.pointBuffer());
5349 meshOps::pWrite(proc, scMap.oldPointBuffer());
5351 // Send connectivity (points, edges, faces, cells, etc)
5352 forAll(scMap.entityBuffer(), bufferI)
5353 {
5354 if (scMap.entityBuffer(bufferI).size())
5355 {
5356 meshOps::pWrite(proc, scMap.entityBuffer(bufferI));
5357 }
5358 }
5360 // Obtain references
5361 const coupledInfo& rPM = recvMeshes_[pI];
5362 const coupleMap& rcMap = rPM.map();
5364 // First read entity sizes.
5365 meshOps::pRead(proc, rcMap.nEntities());
5367 if (debug > 3)
5368 {
5369 Pout<< "Receiving from [" << proc << "]:: nEntities: "
5370 << rcMap.nEntities() << endl;
5371 }
5373 // Size the buffers.
5374 rcMap.allocateBuffers();
5376 // Receive the pointBuffers
5377 meshOps::pRead(proc, rcMap.pointBuffer());
5378 meshOps::pRead(proc, rcMap.oldPointBuffer());
5380 // Receive connectivity (points, edges, faces, cells, etc)
5381 forAll(rcMap.entityBuffer(), bufferI)
5382 {
5383 if (rcMap.entityBuffer(bufferI).size())
5384 {
5385 meshOps::pRead(proc, rcMap.entityBuffer(bufferI));
5386 }
5387 }
5388 }
5390 // We won't wait for all transfers to complete for the moment.
5391 // Meanwhile, do some other useful work, if possible.
5392 }
5395 // Build patch sub-mesh for a specified processor
5396 // - At this point, procIndices is available as a sorted list
5397 // of neighbouring processors.
5398 void dynamicTopoFvMesh::buildProcessorPatchMesh
5399 (
5400 coupledInfo& subMesh,
5401 Map<labelList>& commonCells
5402 )
5403 {
5404 label nP = 0, nE = 0, nF = 0, nC = 0;
5406 // Obtain references
5407 const coupleMap& cMap = subMesh.map();
5408 const labelList& subMeshPoints = cMap.subMeshPoints();
5409 const List<labelPair>& globalProcPoints = cMap.globalProcPoints();
5411 Map<label>& rPointMap = cMap.reverseEntityMap(coupleMap::POINT);
5412 Map<label>& rEdgeMap = cMap.reverseEntityMap(coupleMap::EDGE);
5413 Map<label>& rFaceMap = cMap.reverseEntityMap(coupleMap::FACE);
5414 Map<label>& rCellMap = cMap.reverseEntityMap(coupleMap::CELL);
5416 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
5417 Map<label>& edgeMap = cMap.entityMap(coupleMap::EDGE);
5418 Map<label>& faceMap = cMap.entityMap(coupleMap::FACE);
5419 Map<label>& cellMap = cMap.entityMap(coupleMap::CELL);
5421 // Add all cells connected to points on the subMeshPoints list
5422 label proc = -1;
5424 if (cMap.masterIndex() == Pstream::myProcNo())
5425 {
5426 proc = cMap.slaveIndex();
5427 }
5428 else
5429 {
5430 proc = cMap.masterIndex();
5431 }
5433 // Check if this is a direct neighbour
5434 const polyBoundaryMesh& boundary = boundaryMesh();
5436 /*
5437 // unused code section - to be removed
5438 // HR: 30/11/2013
5440 label patchIndex = -1;
5442 forAll(boundary, patchI)
5443 {
5444 if (isA<processorPolyPatch>(boundary[patchI]))
5445 {
5446 const processorPolyPatch& pp =
5447 (
5448 refCast<const processorPolyPatch>(boundary[patchI])
5449 );
5451 if (pp.neighbProcNo() == proc)
5452 {
5453 patchIndex = patchI;
5454 break;
5455 }
5456 }
5457 }
5458 */
5460 // Add sub-mesh points first.
5461 // Additional halo points will be added later.
5462 forAll(subMeshPoints, pointI)
5463 {
5464 pointMap.insert(nP, subMeshPoints[pointI]);
5465 rPointMap.insert(subMeshPoints[pointI], nP);
5466 nP++;
5467 }
5469 // Set the number of points (shared) at this point.
5470 cMap.nEntities(coupleMap::SHARED_POINT) = nP;
5472 // Size up the point-buffer with global index information.
5473 // Direct neighbours do not require any addressing.
5474 labelList& gpBuffer = cMap.entityBuffer(coupleMap::POINT);
5475 gpBuffer.setSize(globalProcPoints.size(), -1);
5477 forAll(globalProcPoints, pointI)
5478 {
5479 pointMap.insert(nP, globalProcPoints[pointI].first());
5480 rPointMap.insert(globalProcPoints[pointI].first(), nP);
5481 nP++;
5483 // Fill in buffer with global point index
5484 gpBuffer[pointI] = globalProcPoints[pointI].second();
5485 }
5487 // Set the number of points (shared + global) at this point.
5488 cMap.nEntities(coupleMap::GLOBAL_POINT) = nP;
5490 labelHashSet localCommonCells;
5492 // Detect all cells surrounding shared points.
5493 if (twoDMesh_)
5494 {
5495 // No pointEdges structure, so loop through all cells
5496 // to check for those connected to subMeshPoints
5497 forAll(cells_, cellI)
5498 {
5499 const cell& cellToCheck = cells_[cellI];
5500 const labelList cellPoints = cellToCheck.labels(faces_);
5502 forAll(cellPoints, pointI)
5503 {
5504 if (rPointMap.found(cellPoints[pointI]))
5505 {
5506 if (commonCells.found(cellI))
5507 {
5508 // Add locally common cells at the end.
5509 if (!localCommonCells.found(cellI))
5510 {
5511 localCommonCells.insert(cellI);
5512 }
5514 // Check if the processor exists on the list
5515 // and if not, add it.
5516 labelList& procList = commonCells[cellI];
5518 if (findIndex(procList, proc) == -1)
5519 {
5520 meshOps::sizeUpList(proc, procList);
5521 }
5522 }
5523 else
5524 {
5525 cellMap.insert(nC, cellI);
5526 rCellMap.insert(cellI, nC);
5527 nC++;
5529 commonCells.insert(cellI, labelList(1, proc));
5530 }
5532 break;
5533 }
5534 }
5535 }
5536 }
5537 else
5538 {
5539 forAllConstIter(Map<label>, rPointMap, pIter)
5540 {
5541 // Loop through pointEdges for this point.
5542 const labelList& pEdges = pointEdges_[pIter.key()];
5544 forAll(pEdges, edgeI)
5545 {
5546 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
5548 forAll(eFaces, faceI)
5549 {
5550 label own = owner_[eFaces[faceI]];
5551 label nei = neighbour_[eFaces[faceI]];
5553 // Check owner cell
5554 if (!rCellMap.found(own))
5555 {
5556 if (commonCells.found(own))
5557 {
5558 // Add locally common cells at the end.
5559 if (!localCommonCells.found(own))
5560 {
5561 localCommonCells.insert(own);
5562 }
5564 // Check if the processor exists on the list
5565 // and if not, add it.
5566 labelList& procList = commonCells[own];
5568 if (findIndex(procList, proc) == -1)
5569 {
5570 meshOps::sizeUpList(proc, procList);
5571 }
5572 }
5573 else
5574 {
5575 cellMap.insert(nC, own);
5576 rCellMap.insert(own, nC);
5577 nC++;
5579 commonCells.insert(own, labelList(1, proc));
5580 }
5581 }
5583 // Check neighbour cell
5584 if (!rCellMap.found(nei) && nei != -1)
5585 {
5586 if (commonCells.found(nei))
5587 {
5588 // Add locally common cells at the end.
5589 if (!localCommonCells.found(nei))
5590 {
5591 localCommonCells.insert(nei);
5592 }
5594 // Check if the processor exists on the list
5595 // and if not, add it.
5596 labelList& procList = commonCells[nei];
5598 if (findIndex(procList, proc) == -1)
5599 {
5600 meshOps::sizeUpList(proc, procList);
5601 }
5602 }
5603 else
5604 {
5605 cellMap.insert(nC, nei);
5606 rCellMap.insert(nei, nC);
5607 nC++;
5609 commonCells.insert(nei, labelList(1, proc));
5610 }
5611 }
5612 }
5613 }
5614 }
5615 }
5617 // Now add locally common cells.
5618 forAllConstIter(labelHashSet, localCommonCells, cIter)
5619 {
5620 cellMap.insert(nC, cIter.key());
5621 rCellMap.insert(cIter.key(), nC);
5622 nC++;
5623 }
5625 // Keep track of inserted boundary face indices
5626 // - Add exposed internal faces to a 'default' patch
5627 // at the end of the list.
5628 labelList bdyFaceSizes(boundary.size() + 1, 0);
5629 labelList bdyFaceStarts(boundary.size() + 1, 0);
5630 List<Map<label> > bdyFaceIndices(boundary.size() + 1);
5632 // Allocate the faceMap. Since the number of internal faces is unknown,
5633 // detect internal ones first and update boundaries later.
5634 label sumNFE = 0;
5636 forAllConstIter(Map<label>, rCellMap, cIter)
5637 {
5638 const cell& thisCell = cells_[cIter.key()];
5640 forAll(thisCell, faceI)
5641 {
5642 label fIndex = thisCell[faceI];
5644 if (!rFaceMap.found(fIndex))
5645 {
5646 // Determine the patch index
5647 label patchID = whichPatch(fIndex);
5649 if (patchID == -1)
5650 {
5651 // Internal face. Check if this needs to
5652 // be added to the 'default' patch.
5653 label own = owner_[fIndex];
5654 label nei = neighbour_[fIndex];
5656 if (rCellMap.found(own) && rCellMap.found(nei))
5657 {
5658 faceMap.insert(nF, fIndex);
5659 rFaceMap.insert(fIndex, nF);
5660 nF++;
5661 }
5662 else
5663 {
5664 // Update boundary maps.
5665 label bfI = bdyFaceSizes[boundary.size()]++;
5667 // Skip faceMap and update the reverseMap for now.
5668 // faceMap will be updated once all
5669 // internal faces have been detected.
5670 rFaceMap.insert(fIndex, bfI);
5671 bdyFaceIndices[boundary.size()].insert(bfI, fIndex);
5672 }
5673 }
5674 else
5675 {
5676 // Update boundary maps.
5677 label bfI = bdyFaceSizes[patchID]++;
5679 // Skip faceMap and update the reverseMap for now.
5680 // faceMap will be updated once all
5681 // internal faces have been detected.
5682 rFaceMap.insert(fIndex, bfI);
5683 bdyFaceIndices[patchID].insert(bfI, fIndex);
5684 }
5686 // Accumulate face sizes
5687 sumNFE += faces_[fIndex].size();
5688 }
5689 }
5690 }
5692 // Set the number of internal faces at this point
5693 cMap.nEntities(coupleMap::INTERNAL_FACE) = nF;
5695 // Set patch starts
5696 bdyFaceStarts[0] = nF;
5698 for (label i = 1; i < bdyFaceStarts.size(); i++)
5699 {
5700 bdyFaceStarts[i] = bdyFaceStarts[i-1] + bdyFaceSizes[i-1];
5701 }
5703 // Update faceMap and reverseFaceMap for boundaries
5704 forAll(bdyFaceIndices, patchI)
5705 {
5706 label pStart = bdyFaceStarts[patchI];
5708 forAllConstIter(Map<label>, bdyFaceIndices[patchI], fIter)
5709 {
5710 faceMap.insert(fIter.key() + pStart, fIter());
5711 rFaceMap[fIter()] = fIter.key() + pStart;
5712 }
5714 // Update face-count
5715 nF += bdyFaceSizes[patchI];
5716 }
5718 // Keep track of inserted boundary edge indices
5719 // - Add exposed internal edges to a 'default' patch
5720 // at the end of the list.
5721 labelList bdyEdgeSizes(boundary.size() + 1, 0);
5722 labelList bdyEdgeStarts(boundary.size() + 1, 0);
5723 List<Map<label> > bdyEdgeIndices(boundary.size() + 1);
5725 // Allocate the edgeMap. Since the number of internal edges is unknown,
5726 // detect internal ones first and update boundaries later.
5727 forAllConstIter(Map<label>, rFaceMap, fIter)
5728 {
5729 const labelList& fEdges = faceEdges_[fIter.key()];
5731 forAll(fEdges, edgeI)
5732 {
5733 label eIndex = fEdges[edgeI];
5735 if (!rEdgeMap.found(eIndex))
5736 {
5737 // Determine the patch index
5738 label patchID = whichEdgePatch(eIndex);
5740 if (patchID == -1)
5741 {
5742 bool boundaryEdge = false;
5744 // Check if any cells touching edgeFaces
5745 // do not belong to the cellMap.
5746 const labelList& eFaces = edgeFaces_[eIndex];
5748 forAll(eFaces, faceI)
5749 {
5750 label fIndex = eFaces[faceI];
5752 label own = owner_[fIndex];
5753 label nei = neighbour_[fIndex];
5755 if (!rCellMap.found(own) || !rCellMap.found(nei))
5756 {
5757 boundaryEdge = true;
5758 break;
5759 }
5760 }
5762 if (boundaryEdge)
5763 {
5764 // Update boundary maps.
5765 label beI = bdyEdgeSizes[boundary.size()]++;
5767 // Skip edgeMap and update the reverseMap for now.
5768 // edgeMap will be updated once all
5769 // internal edges have been detected.
5770 rEdgeMap.insert(eIndex, beI);
5771 bdyEdgeIndices[boundary.size()].insert(beI, eIndex);
5772 }
5773 else
5774 {
5775 // Internal edge
5776 edgeMap.insert(nE, eIndex);
5777 rEdgeMap.insert(eIndex, nE);
5778 nE++;
5779 }
5780 }
5781 else
5782 {
5783 // Update boundary maps.
5784 label beI = bdyEdgeSizes[patchID]++;
5786 // Skip edgeMap and update the reverseMap for now.
5787 // edgeMap will be updated once all
5788 // internal edges have been detected.
5789 rEdgeMap.insert(eIndex, beI);
5790 bdyEdgeIndices[patchID].insert(beI, eIndex);
5791 }
5792 }
5793 }
5794 }
5796 // Set the number of internal edges at this point
5797 cMap.nEntities(coupleMap::INTERNAL_EDGE) = nE;
5799 // Set patch starts
5800 bdyEdgeStarts[0] = nE;
5802 for (label i = 1; i < bdyEdgeStarts.size(); i++)
5803 {
5804 bdyEdgeStarts[i] = bdyEdgeStarts[i-1] + bdyEdgeSizes[i-1];
5805 }
5807 // Update edgeMap and reverseEdgeMap for boundaries
5808 forAll(bdyEdgeIndices, patchI)
5809 {
5810 label pStart = bdyEdgeStarts[patchI];
5812 forAllConstIter(Map<label>, bdyEdgeIndices[patchI], eIter)
5813 {
5814 edgeMap.insert(eIter.key() + pStart, eIter());
5815 rEdgeMap[eIter()] = eIter.key() + pStart;
5816 }
5818 // Update edge-count
5819 nE += bdyEdgeSizes[patchI];
5820 }
5822 // Set additional points in the pointMap
5823 forAllConstIter(Map<label>, rEdgeMap, eIter)
5824 {
5825 const edge& thisEdge = edges_[eIter.key()];
5827 if (!rPointMap.found(thisEdge[0]))
5828 {
5829 pointMap.insert(nP, thisEdge[0]);
5830 rPointMap.insert(thisEdge[0], nP);
5831 nP++;
5832 }
5834 if (!rPointMap.found(thisEdge[1]))
5835 {
5836 pointMap.insert(nP, thisEdge[1]);
5837 rPointMap.insert(thisEdge[1], nP);
5838 nP++;
5839 }
5840 }
5842 // Assign sizes to the mesh
5843 cMap.nEntities(coupleMap::POINT) = nP;
5844 cMap.nEntities(coupleMap::EDGE) = nE;
5845 cMap.nEntities(coupleMap::FACE) = nF;
5846 cMap.nEntities(coupleMap::CELL) = nC;
5847 cMap.nEntities(coupleMap::NFE_SIZE) = sumNFE;
5848 cMap.nEntities(coupleMap::NBDY) = boundary.size() + 1;
5850 // Size up buffers and fill them
5851 cMap.allocateBuffers();
5853 pointField& pBuffer = cMap.pointBuffer();
5854 pointField& opBuffer = cMap.oldPointBuffer();
5856 forAllConstIter(Map<label>, pointMap, pIter)
5857 {
5858 pBuffer[pIter.key()] = points_[pIter()];
5859 opBuffer[pIter.key()] = oldPoints_[pIter()];
5860 }
5862 label index = 0;
5864 // Edge buffer size: 2 points for every edge
5865 labelList& eBuffer = cMap.entityBuffer(coupleMap::EDGE);
5867 for (label i = 0; i < nE; i++)
5868 {
5869 label eIndex = edgeMap[i];
5870 const edge& edgeToCheck = edges_[eIndex];
5872 eBuffer[index++] = rPointMap[edgeToCheck[0]];
5873 eBuffer[index++] = rPointMap[edgeToCheck[1]];
5874 }
5876 index = 0;
5877 face thisFace;
5878 labelList& fBuffer = cMap.entityBuffer(coupleMap::FACE);
5879 labelList& fOwner = cMap.entityBuffer(coupleMap::OWNER);
5880 labelList& fNeighbour = cMap.entityBuffer(coupleMap::NEIGHBOUR);
5881 labelList& feBuffer = cMap.entityBuffer(coupleMap::FACE_EDGE);
5883 for (label i = 0; i < nF; i++)
5884 {
5885 label fIndex = faceMap[i];
5886 label own = owner_[fIndex];
5887 label nei = neighbour_[fIndex];
5889 // Fetch mapped addressing
5890 label fOwn = rCellMap.found(own) ? rCellMap[own] : -1;
5891 label fNei = rCellMap.found(nei) ? rCellMap[nei] : -1;
5893 if (fOwn > -1)
5894 {
5895 // Check if this face is pointed the right way
5896 if ((fNei > -1) && (fNei < fOwn))
5897 {
5898 thisFace = faces_[fIndex].reverseFace();
5899 fOwner[i] = fNei;
5900 fNeighbour[i] = fOwn;
5901 }
5902 else
5903 {
5904 thisFace = faces_[fIndex];
5905 fOwner[i] = fOwn;
5907 if (fNei > -1)
5908 {
5909 fNeighbour[i] = fNei;
5910 }
5911 }
5912 }
5913 else
5914 {
5915 // This face is pointed the wrong way.
5916 thisFace = faces_[fIndex].reverseFace();
5917 fOwner[i] = fNei;
5918 }
5920 const labelList& fEdges = faceEdges_[fIndex];
5922 forAll(fEdges, indexI)
5923 {
5924 fBuffer[index] = rPointMap[thisFace[indexI]];
5925 feBuffer[index] = rEdgeMap[fEdges[indexI]];
5927 index++;
5928 }
5929 }
5931 // Fill variable size face-list sizes for 2D
5932 if (twoDMesh_)
5933 {
5934 labelList& nfeBuffer = cMap.entityBuffer(coupleMap::NFE_BUFFER);
5936 forAllConstIter(Map<label>, faceMap, fIter)
5937 {
5938 nfeBuffer[fIter.key()] = faces_[fIter()].size();
5939 }
5940 }
5942 // Fill in boundary information
5943 cMap.entityBuffer(coupleMap::FACE_STARTS) = bdyFaceStarts;
5944 cMap.entityBuffer(coupleMap::FACE_SIZES) = bdyFaceSizes;
5945 cMap.entityBuffer(coupleMap::EDGE_STARTS) = bdyEdgeStarts;
5946 cMap.entityBuffer(coupleMap::EDGE_SIZES) = bdyEdgeSizes;
5948 labelList& ptBuffer = cMap.entityBuffer(coupleMap::PATCH_ID);
5950 // Fill types for all but the last one (which is default).
5951 forAll(boundary, patchI)
5952 {
5953 if (isA<processorPolyPatch>(boundary[patchI]))
5954 {
5955 // Fetch the neighbouring processor ID
5956 const processorPolyPatch& pp =
5957 (
5958 refCast<const processorPolyPatch>(boundary[patchI])
5959 );
5961 // Set processor patches to a special type
5962 ptBuffer[patchI] = -2 - pp.neighbProcNo();
5963 }
5964 else
5965 {
5966 // Conventional physical patch. Make an identical
5967 // map, since physical boundaries are present on
5968 // all processors.
5969 ptBuffer[patchI] = patchI;
5970 }
5971 }
5973 // Fill the default patch as 'internal'
5974 ptBuffer[boundary.size()] = -1;
5976 // Make a temporary dictionary for patch construction
5977 dictionary patchDict;
5979 // Specify the list of patch names and types
5980 wordList patchNames(ptBuffer.size());
5982 forAll(patchNames, patchI)
5983 {
5984 // Add a new subDictionary
5985 dictionary patchSubDict;
5987 if (patchI == patchNames.size() - 1)
5988 {
5989 // Artificially set the last patch
5991 // Set name
5992 patchNames[patchI] = "defaultPatch";
5994 // Add type
5995 patchSubDict.add("type", "empty");
5997 // Add start / size
5998 patchSubDict.add("startFace", bdyFaceStarts[patchI]);
5999 patchSubDict.add("nFaces", bdyFaceSizes[patchI]);
6000 }
6001 else
6002 if (ptBuffer[patchI] <= -2)
6003 {
6004 // Back out the neighbouring processor ID
6005 label neiProcNo = Foam::mag(ptBuffer[patchI] + 2);
6007 // Set name
6008 patchNames[patchI] =
6009 (
6010 "procBoundary"
6011 + Foam::name(Pstream::myProcNo())
6012 + "to"
6013 + Foam::name(neiProcNo)
6014 );
6016 // Add type
6017 patchSubDict.add("type", "processor");
6019 // Add start / size
6020 patchSubDict.add("startFace", bdyFaceStarts[patchI]);
6021 patchSubDict.add("nFaces", bdyFaceSizes[patchI]);
6023 // Add processor-specific info
6024 patchSubDict.add("myProcNo", Pstream::myProcNo());
6025 patchSubDict.add("neighbProcNo", neiProcNo);
6026 }
6027 else
6028 {
6029 // Set name
6030 patchNames[patchI] = boundary[ptBuffer[patchI]].name();
6032 // Add type
6033 patchSubDict.add("type", boundary[ptBuffer[patchI]].type());
6035 // Add start / size
6036 patchSubDict.add("startFace", bdyFaceStarts[patchI]);
6037 patchSubDict.add("nFaces", bdyFaceSizes[patchI]);
6038 }
6040 // Add subdictionary
6041 patchDict.add(patchNames[patchI], patchSubDict);
6042 }
6044 // Set the autoPtr.
6045 subMesh.setMesh
6046 (
6047 proc,
6048 new dynamicTopoFvMesh
6049 (
6050 (*this),
6051 IOobject
6052 (
6053 fvMesh::defaultRegion,
6054 time().timeName(),
6055 time()
6056 ),
6057 xferCopy(cMap.pointBuffer()),
6058 xferCopy(cMap.oldPointBuffer()),
6059 xferCopy(cMap.edges()),
6060 xferCopy(cMap.faces()),
6061 xferCopy(cMap.faceEdges()),
6062 xferCopy(cMap.owner()),
6063 xferCopy(cMap.neighbour()),
6064 bdyFaceStarts,
6065 bdyFaceSizes,
6066 bdyEdgeStarts,
6067 bdyEdgeSizes,
6068 patchNames,
6069 patchDict
6070 )
6071 );
6073 // Set maps as built.
6074 subMesh.setBuiltMaps();
6076 // For debugging purposes...
6077 if (debug > 3)
6078 {
6079 Pout<< "Writing out sent subMesh for processor: "
6080 << proc << endl;
6082 writeVTK
6083 (
6084 "psMesh_"
6085 + Foam::name(Pstream::myProcNo())
6086 + "to"
6087 + Foam::name(proc),
6088 rCellMap.toc()
6089 );
6091 // Write out patch information
6092 Pout<< " faceStarts: " << bdyFaceStarts << nl
6093 << " faceSizes: " << bdyFaceSizes << nl
6094 << " edgeStarts: " << bdyEdgeStarts << nl
6095 << " edgeSizes: " << bdyEdgeSizes << nl
6096 << " patchTypes: " << ptBuffer << endl;
6097 }
6098 }
6101 // Build coupled maps for locally coupled patches.
6102 // - Performs a geometric match initially, since the mesh provides
6103 // no explicit information for topological coupling.
6104 // - For each subsequent topology change, maps are updated during
6105 // the re-ordering stage.
6106 void dynamicTopoFvMesh::buildLocalCoupledMaps()
6107 {
6108 if (patchCoupling_.empty())
6109 {
6110 return;
6111 }
6113 if (debug)
6114 {
6115 Info<< "Building local coupled maps...";
6116 }
6118 // Check if a geometric tolerance has been specified.
6119 const boundBox& box = polyMesh::bounds();
6120 scalar relTol = debug::tolerances("patchFaceMatchTol", 1e-4);
6122 // Compute tolerance
6123 scalar tol = relTol * box.mag();
6125 const polyBoundaryMesh& boundary = boundaryMesh();
6127 forAll(patchCoupling_, patchI)
6128 {
6129 if (!patchCoupling_(patchI))
6130 {
6131 continue;
6132 }
6134 // Build maps only for the first time.
6135 // Information is subsequently mapped for each topo-change.
6136 if (patchCoupling_[patchI].builtMaps())
6137 {
6138 continue;
6139 }
6141 // Deal with cyclics later
6142 if (isA<cyclicPolyPatch>(boundary[patchI]))
6143 {
6144 continue;
6145 }
6147 const coupleMap& cMap = patchCoupling_[patchI].map();
6149 // Map points on each patch.
6150 const labelList& mP = boundary[cMap.masterIndex()].meshPoints();
6151 const labelList& sP = boundary[cMap.slaveIndex()].meshPoints();
6153 // Sanity check: Do patches have equal number of entities?
6154 if (mP.size() != sP.size())
6155 {
6156 FatalErrorIn("void dynamicTopoFvMesh::buildLocalCoupledMaps()")
6157 << "Patch point sizes are not consistent."
6158 << abort(FatalError);
6159 }
6161 // Check if maps were read from disk.
6162 // If so, geometric checking is unnecessary.
6163 if (cMap.entityMap(coupleMap::POINT).size() == 0)
6164 {
6165 // Match points geometrically
6166 labelList pointMap(mP.size(), -1);
6168 bool matchedAll =
6169 (
6170 matchPoints
6171 (
6172 boundary[cMap.masterIndex()].localPoints(),
6173 boundary[cMap.slaveIndex()].localPoints(),
6174 scalarField(mP.size(), tol),
6175 false,
6176 pointMap
6177 )
6178 );
6180 if (!matchedAll)
6181 {
6182 FatalErrorIn("void dynamicTopoFvMesh::buildLocalCoupledMaps()")
6183 << " Failed to match all points"
6184 << " within a tolerance of: " << tol << nl
6185 << " relTol: " << relTol << nl
6186 << abort(FatalError);
6187 }
6189 // Now build maps
6190 forAll(mP, indexI)
6191 {
6192 label masterIndex = mP[indexI];
6193 label slaveIndex = sP[pointMap[indexI]];
6195 cMap.mapSlave
6196 (
6197 coupleMap::POINT,
6198 masterIndex,
6199 slaveIndex
6200 );
6202 cMap.mapMaster
6203 (
6204 coupleMap::POINT,
6205 slaveIndex,
6206 masterIndex
6207 );
6208 }
6209 }
6211 const labelListList& mpF = boundary[cMap.masterIndex()].pointFaces();
6212 const labelListList& spF = boundary[cMap.slaveIndex()].pointFaces();
6214 // Fetch reference to maps
6215 const Map<label>& pMap = cMap.entityMap(coupleMap::POINT);
6216 const Map<label>& eMap = cMap.entityMap(coupleMap::EDGE);
6217 const Map<label>& fMap = cMap.entityMap(coupleMap::FACE);
6219 // Now that all points are matched,
6220 // perform topological matching for higher entities.
6221 forAll(mP, indexI)
6222 {
6223 // Fetch global point indices
6224 label mp = mP[indexI];
6225 label sp = pMap[mp];
6227 // Fetch local point indices
6228 label lmp = boundary[cMap.masterIndex()].whichPoint(mp);
6229 label lsp = boundary[cMap.slaveIndex()].whichPoint(sp);
6231 // Fetch patch starts
6232 label mStart = boundary[cMap.masterIndex()].start();
6233 label sStart = boundary[cMap.slaveIndex()].start();
6235 // Match faces for both 2D and 3D.
6236 const labelList& mpFaces = mpF[lmp];
6237 const labelList& spFaces = spF[lsp];
6239 forAll(mpFaces, faceI)
6240 {
6241 // Fetch the global face index
6242 label mfIndex = (mStart + mpFaces[faceI]);
6244 if (fMap.found(mfIndex))
6245 {
6246 continue;
6247 }
6249 const face& mFace = faces_[mfIndex];
6251 if (twoDMesh_)
6252 {
6253 // Set up a comparison face.
6254 face cFace(4);
6256 // Configure the face for comparison.
6257 forAll(mFace, pointI)
6258 {
6259 cFace[pointI] = pMap[mFace[pointI]];
6260 }
6262 bool matched = false;
6264 forAll(spFaces, faceJ)
6265 {
6266 // Fetch the global face index
6267 label sfIndex = (sStart + spFaces[faceJ]);
6269 const face& sFace = faces_[sfIndex];
6271 if (face::compare(cFace, sFace))
6272 {
6273 // Found the slave. Add a map entry
6274 cMap.mapSlave
6275 (
6276 coupleMap::FACE,
6277 mfIndex,
6278 sfIndex
6279 );
6281 cMap.mapMaster
6282 (
6283 coupleMap::FACE,
6284 sfIndex,
6285 mfIndex
6286 );
6288 matched = true;
6290 break;
6291 }
6292 }
6294 if (!matched)
6295 {
6296 FatalErrorIn
6297 (
6298 "void dynamicTopoFvMesh::buildLocalCoupledMaps()"
6299 )
6300 << " Failed to match face: "
6301 << mfIndex << ": " << mFace << nl
6302 << " Comparison face: " << cFace
6303 << abort(FatalError);
6304 }
6305 }
6306 else
6307 {
6308 label slaveFaceIndex = -1;
6310 // Set up a comparison face.
6311 triFace cFace
6312 (
6313 pMap[mFace[0]],
6314 pMap[mFace[1]],
6315 pMap[mFace[2]]
6316 );
6318 forAll(spFaces, faceJ)
6319 {
6320 // Fetch the global face index
6321 label sfIndex = (sStart + spFaces[faceJ]);
6323 const face& sFace = faces_[sfIndex];
6325 if (triFace::compare(cFace, triFace(sFace)))
6326 {
6327 // Found the slave. Add a map entry
6328 cMap.mapSlave
6329 (
6330 coupleMap::FACE,
6331 mfIndex,
6332 sfIndex
6333 );
6335 cMap.mapMaster
6336 (
6337 coupleMap::FACE,
6338 sfIndex,
6339 mfIndex
6340 );
6342 slaveFaceIndex = sfIndex;
6344 break;
6345 }
6346 }
6348 if (slaveFaceIndex == -1)
6349 {
6350 FatalErrorIn
6351 (
6352 "void dynamicTopoFvMesh::buildLocalCoupledMaps()"
6353 )
6354 << " Failed to match face: "
6355 << mfIndex << ": " << mFace << nl
6356 << " Comparison face: " << cFace
6357 << abort(FatalError);
6358 }
6360 // Match all edges on this face as well.
6361 const labelList& mfEdges = faceEdges_[mfIndex];
6362 const labelList& sfEdges = faceEdges_[slaveFaceIndex];
6364 forAll(mfEdges, edgeI)
6365 {
6366 if (eMap.found(mfEdges[edgeI]))
6367 {
6368 continue;
6369 }
6371 const edge& mEdge = edges_[mfEdges[edgeI]];
6373 // Configure a comparison edge.
6374 edge cEdge(pMap[mEdge[0]], pMap[mEdge[1]]);
6376 bool matchedEdge = false;
6378 forAll(sfEdges, edgeJ)
6379 {
6380 const edge& sEdge = edges_[sfEdges[edgeJ]];
6382 if (cEdge == sEdge)
6383 {
6384 // Found the slave. Add a map entry
6385 cMap.mapSlave
6386 (
6387 coupleMap::EDGE,
6388 mfEdges[edgeI],
6389 sfEdges[edgeJ]
6390 );
6392 cMap.mapMaster
6393 (
6394 coupleMap::EDGE,
6395 sfEdges[edgeJ],
6396 mfEdges[edgeI]
6397 );
6399 matchedEdge = true;
6401 break;
6402 }
6403 }
6405 if (!matchedEdge)
6406 {
6407 FatalErrorIn
6408 (
6409 "void dynamicTopoFvMesh::"
6410 "buildLocalCoupledMaps()"
6411 )
6412 << " Failed to match edge: "
6413 << mfEdges[edgeI] << ": "
6414 << mEdge << nl
6415 << " cEdge: " << cEdge
6416 << abort(FatalError);
6417 }
6418 }
6419 }
6420 }
6421 }
6423 // Set maps as built
6424 patchCoupling_[patchI].setBuiltMaps();
6425 }
6427 // Build maps for cyclics
6428 forAll(patchCoupling_, patchI)
6429 {
6430 if (!patchCoupling_(patchI))
6431 {
6432 continue;
6433 }
6435 // Build maps only for the first time.
6436 // Information is subsequently mapped for each topo-change.
6437 if (patchCoupling_[patchI].builtMaps())
6438 {
6439 continue;
6440 }
6442 // Skip if not cyclic
6443 if (!isA<cyclicPolyPatch>(boundary[patchI]))
6444 {
6445 continue;
6446 }
6448 const coupleMap& cMap = patchCoupling_[patchI].map();
6450 // Match faces and points in one go
6451 label halfSize = boundary[patchI].size() / 2;
6452 label patchStart = boundary[patchI].start();
6454 // Fetch reference to map
6455 const Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
6457 for (label i = 0; i < halfSize; i++)
6458 {
6459 label half0Index = (i + patchStart);
6460 label half1Index = (i + halfSize + patchStart);
6462 // Map the face
6463 cMap.mapSlave
6464 (
6465 coupleMap::FACE,
6466 half0Index,
6467 half1Index
6468 );
6470 cMap.mapMaster
6471 (
6472 coupleMap::FACE,
6473 half1Index,
6474 half0Index
6475 );
6477 const face& half0Face = faces_[half0Index];
6478 const face& half1Face = faces_[half1Index];
6480 label fS = half0Face.size();
6482 forAll(half0Face, pointI)
6483 {
6484 if (pointMap.found(half0Face[pointI]))
6485 {
6486 continue;
6487 }
6489 label masterIndex = half0Face[pointI];
6490 label slaveIndex = half1Face[(fS - pointI) % fS];
6492 // Map the point
6493 cMap.mapSlave
6494 (
6495 coupleMap::POINT,
6496 masterIndex,
6497 slaveIndex
6498 );
6500 cMap.mapMaster
6501 (
6502 coupleMap::POINT,
6503 slaveIndex,
6504 masterIndex
6505 );
6506 }
6508 if (twoDMesh_)
6509 {
6510 continue;
6511 }
6513 // Fetch reference to map
6514 const Map<label>& edgeMap = cMap.entityMap(coupleMap::EDGE);
6516 const labelList& f0Edges = faceEdges_[half0Index];
6517 const labelList& f1Edges = faceEdges_[half1Index];
6519 forAll(f0Edges, edgeI)
6520 {
6521 if (edgeMap.found(f0Edges[edgeI]))
6522 {
6523 continue;
6524 }
6526 const edge& mEdge = edges_[f0Edges[edgeI]];
6528 // Build a comparison edge
6529 edge cEdge(pointMap[mEdge[0]], pointMap[mEdge[1]]);
6531 forAll(f1Edges, edgeJ)
6532 {
6533 const edge& sEdge = edges_[f1Edges[edgeJ]];
6535 if (sEdge == cEdge)
6536 {
6537 // Map the edge
6538 cMap.mapSlave
6539 (
6540 coupleMap::EDGE,
6541 f0Edges[edgeI],
6542 f1Edges[edgeJ]
6543 );
6545 cMap.mapMaster
6546 (
6547 coupleMap::EDGE,
6548 f1Edges[edgeJ],
6549 f0Edges[edgeI]
6550 );
6552 break;
6553 }
6554 }
6555 }
6556 }
6557 }
6559 if (debug)
6560 {
6561 Info<< "Done." << endl;
6562 }
6563 }
6566 // Build coupled maps for coupled processor patches
6567 void dynamicTopoFvMesh::buildProcessorCoupledMaps()
6568 {
6569 if (procIndices_.empty())
6570 {
6571 return;
6572 }
6574 Map<label> nPrc;
6576 const polyBoundaryMesh& boundary = boundaryMesh();
6578 // Build a list of nearest neighbours.
6579 forAll(boundary, patchI)
6580 {
6581 if (isA<processorPolyPatch>(boundary[patchI]))
6582 {
6583 const processorPolyPatch& pp =
6584 (
6585 refCast<const processorPolyPatch>(boundary[patchI])
6586 );
6588 nPrc.insert(pp.neighbProcNo(), patchI);
6589 }
6590 }
6592 // Wait for all transfers to complete.
6593 meshOps::waitForBuffers();
6595 // Put un-matched faces in a list.
6596 labelHashSet unMatchedFaces;
6598 forAll(procIndices_, pI)
6599 {
6600 label proc = procIndices_[pI];
6602 coupledInfo& rPM = recvMeshes_[pI];
6603 const coupleMap& cMap = rPM.map();
6604 const labelList& ptBuffer = cMap.entityBuffer(coupleMap::PATCH_ID);
6606 // Fetch reference to patch starts / sizes
6607 const labelList& faceStarts = cMap.entityBuffer(coupleMap::FACE_STARTS);
6608 const labelList& faceSizes = cMap.entityBuffer(coupleMap::FACE_SIZES);
6609 const labelList& edgeStarts = cMap.entityBuffer(coupleMap::EDGE_STARTS);
6610 const labelList& edgeSizes = cMap.entityBuffer(coupleMap::EDGE_SIZES);
6612 // Make a temporary dictionary for patch construction
6613 dictionary patchDict;
6615 // Specify the list of patch names and types
6616 wordList patchNames(ptBuffer.size());
6618 forAll(patchNames, patchI)
6619 {
6620 // Add a new subDictionary
6621 dictionary patchSubDict;
6623 if (patchI == patchNames.size() - 1)
6624 {
6625 // Artificially set the last patch
6627 // Set name
6628 patchNames[patchI] = "defaultPatch";
6630 // Add type
6631 patchSubDict.add("type", "empty");
6633 // Add start / size
6634 patchSubDict.add("startFace", faceStarts[patchI]);
6635 patchSubDict.add("nFaces", faceSizes[patchI]);
6636 }
6637 else
6638 if (ptBuffer[patchI] <= -2)
6639 {
6640 // Back out the neighbouring processor ID
6641 label neiProcNo = Foam::mag(ptBuffer[patchI] + 2);
6643 // Set name
6644 patchNames[patchI] =
6645 (
6646 "procBoundary"
6647 + Foam::name(proc)
6648 + "to"
6649 + Foam::name(neiProcNo)
6650 );
6652 // Add type
6653 patchSubDict.add("type", "processor");
6655 // Add start / size
6656 patchSubDict.add("startFace", faceStarts[patchI]);
6657 patchSubDict.add("nFaces", faceSizes[patchI]);
6659 // Add processor-specific info
6660 patchSubDict.add("myProcNo", proc);
6661 patchSubDict.add("neighbProcNo", neiProcNo);
6662 }
6663 else
6664 {
6665 // Set name
6666 patchNames[patchI] = boundary[ptBuffer[patchI]].name();
6668 // Add type
6669 patchSubDict.add("type", boundary[ptBuffer[patchI]].type());
6671 // Add start / size
6672 patchSubDict.add("startFace", faceStarts[patchI]);
6673 patchSubDict.add("nFaces", faceSizes[patchI]);
6674 }
6676 // Add subdictionary
6677 patchDict.add(patchNames[patchI], patchSubDict);
6678 }
6680 // Set the autoPtr.
6681 rPM.setMesh
6682 (
6683 proc,
6684 new dynamicTopoFvMesh
6685 (
6686 (*this),
6687 IOobject
6688 (
6689 fvMesh::defaultRegion,
6690 time().timeName(),
6691 time()
6692 ),
6693 xferCopy(cMap.pointBuffer()),
6694 xferCopy(cMap.oldPointBuffer()),
6695 xferCopy(cMap.edges()),
6696 xferCopy(cMap.faces()),
6697 xferCopy(cMap.faceEdges()),
6698 xferCopy(cMap.owner()),
6699 xferCopy(cMap.neighbour()),
6700 faceStarts,
6701 faceSizes,
6702 edgeStarts,
6703 edgeSizes,
6704 patchNames,
6705 patchDict
6706 )
6707 );
6709 if (debug > 3)
6710 {
6711 Pout<< "Writing out received subMesh for processor: "
6712 << proc << endl;
6714 rPM.subMesh().writeVTK
6715 (
6716 "rPatchMesh_"
6717 + Foam::name(Pstream::myProcNo())
6718 + "to"
6719 + Foam::name(proc),
6720 identity(cMap.nEntities(coupleMap::CELL))
6721 );
6722 }
6724 // First, topologically map points based on subMeshPoints.
6725 const labelList& mP = cMap.subMeshPoints();
6727 label nShared = cMap.nEntities(coupleMap::SHARED_POINT);
6728 label nGlobal = cMap.nEntities(coupleMap::GLOBAL_POINT);
6730 // Sanity check: Do sub-mesh point sizes match?
6731 if (mP.size() != nShared)
6732 {
6733 FatalErrorIn("void dynamicTopoFvMesh::buildProcessorCoupledMaps()")
6734 << " Sub-mesh point sizes don't match." << nl
6735 << " My procID: " << Pstream::myProcNo() << nl
6736 << " Neighbour processor: " << proc << nl
6737 << " mP size: " << mP.size() << nl
6738 << " nShared: " << nShared << nl
6739 << abort(FatalError);
6740 }
6742 if (nPrc.found(proc))
6743 {
6744 // This is an immediate neighbour.
6745 // All patch points must be matched.
6747 // Fetch addressing for this patch.
6748 const processorPolyPatch& pp =
6749 (
6750 refCast<const processorPolyPatch>(boundary[nPrc[proc]])
6751 );
6753 const labelList& neiPoints = pp.neighbPoints();
6755 if (debug)
6756 {
6757 // Find all occurrences of multiply connected points
6758 labelList mIdx = findIndices(neiPoints, -1);
6760 if (mIdx.size())
6761 {
6762 // Write out for post-processing
6763 UIndirectList<label> myIdx(mP, mIdx);
6764 writeVTK("mcPoints", labelList(myIdx), 0);
6766 FatalErrorIn
6767 (
6768 "void dynamicTopoFvMesh::buildProcessorCoupledMaps()"
6769 )
6770 << " Multiply connected point." << nl
6771 << " My procID: " << Pstream::myProcNo() << nl
6772 << " Neighbour processor: " << proc << nl
6773 << " Neighbour Indices: " << mIdx << nl
6774 << " My indices: " << myIdx << nl
6775 << abort(FatalError);
6776 }
6777 }
6779 forAll(mP, pointI)
6780 {
6781 // Add a map entry
6782 cMap.mapSlave
6783 (
6784 coupleMap::POINT,
6785 mP[pointI],
6786 neiPoints[pointI]
6787 );
6789 cMap.mapMaster
6790 (
6791 coupleMap::POINT,
6792 neiPoints[pointI],
6793 mP[pointI]
6794 );
6795 }
6796 }
6798 // Prepare point maps for globally shared points
6799 const labelList& gpB = cMap.entityBuffer(coupleMap::POINT);
6801 // Sanity check: Do global point buffer sizes match?
6802 if ((mP.size() + gpB.size()) != nGlobal)
6803 {
6804 FatalErrorIn("void dynamicTopoFvMesh::buildProcessorCoupledMaps()")
6805 << " Global point sizes don't match"
6806 << " for processor: " << proc << nl
6807 << " mP size: " << mP.size() << nl
6808 << " gpB size: " << gpB.size() << nl
6809 << " Total: " << (mP.size() + gpB.size()) << nl
6810 << " nGlobal: " << nGlobal << nl
6811 << abort(FatalError);
6812 }
6814 // Look at globalPoint addressing for its information.
6815 const globalMeshData& gD = polyMesh::globalData();
6816 const labelList& spA = gD.sharedPointAddr();
6817 const labelList& spL = gD.sharedPointLabels();
6819 forAll(gpB, pointI)
6820 {
6821 // Find my equivalent global point
6822 label maIndex = findIndex(spA, gpB[pointI]);
6824 // Add a map entry
6825 cMap.mapSlave
6826 (
6827 coupleMap::POINT,
6828 spL[maIndex],
6829 (pointI + nShared)
6830 );
6832 cMap.mapMaster
6833 (
6834 coupleMap::POINT,
6835 (pointI + nShared),
6836 spL[maIndex]
6837 );
6838 }
6840 if (debug > 1)
6841 {
6842 const boundBox& box = polyMesh::bounds();
6843 const pointField& sPoints = rPM.subMesh().points_;
6844 const Map<label>& pMap = cMap.entityMap(coupleMap::POINT);
6846 // Fetch relative tolerance
6847 scalar relTol = debug::tolerances("processorMatchTol", 1e-4);
6849 // Compute tolerance
6850 scalar tol = relTol * box.mag();
6852 forAllConstIter(Map<label>, pMap, pIter)
6853 {
6854 scalar dist = mag(points_[pIter.key()] - sPoints[pIter()]);
6856 if (dist > tol)
6857 {
6858 FatalErrorIn
6859 (
6860 "void dynamicTopoFvMesh::buildProcessorCoupledMaps()"
6861 )
6862 << " Failed to match point: " << pIter.key()
6863 << ": " << points_[pIter.key()]
6864 << " with point: " << pIter()
6865 << ": " << sPoints[pIter()] << nl
6866 << " Missed by: " << dist << nl
6867 << " Tolerance: " << tol << nl
6868 << abort(FatalError);
6869 }
6870 }
6871 }
6873 // Set up a comparison face.
6874 face cFace(4);
6876 // Now match all faces connected to master points.
6877 if (nPrc.found(proc))
6878 {
6879 // Fetch a global faceList for the slave subMesh
6880 const faceList& slaveFaces = rPM.subMesh().faces();
6882 // This is an immediate neighbour.
6883 label mStart = boundary[nPrc[proc]].start();
6884 label mSize = boundary[nPrc[proc]].size();
6886 // Abbreviate for convenience
6887 const dynamicTopoFvMesh& sMesh = rPM.subMesh();
6888 const Map<label>& pMap = cMap.entityMap(coupleMap::POINT);
6889 const Map<label>& eMap = cMap.entityMap(coupleMap::EDGE);
6890 const Map<label>& fMap = cMap.entityMap(coupleMap::FACE);
6892 // Fetch global pointFaces for the slave.
6893 const labelListList& spF = rPM.subMesh().pointFaces();
6895 // Match patch faces for both 2D and 3D.
6896 for(label i = 0; i < mSize; i++)
6897 {
6898 // Fetch the global face index
6899 label mfIndex = (mStart + i);
6901 if (fMap.found(mfIndex))
6902 {
6903 continue;
6904 }
6906 const face& mFace = faces_[mfIndex];
6908 // Configure the face for comparison.
6909 forAll(mFace, pointI)
6910 {
6911 cFace[pointI] = pMap[mFace[pointI]];
6912 }
6914 // Fetch pointFaces for the zeroth point.
6915 const labelList& spFaces = spF[cFace[0]];
6917 if (twoDMesh_)
6918 {
6919 bool matched = false;
6921 forAll(spFaces, faceJ)
6922 {
6923 // Fetch the global face index
6924 label sfIndex = spFaces[faceJ];
6926 const face& sFace = slaveFaces[sfIndex];
6928 if (face::compare(cFace, sFace))
6929 {
6930 // Found the slave. Add a map entry
6931 cMap.mapSlave
6932 (
6933 coupleMap::FACE,
6934 mfIndex,
6935 sfIndex
6936 );
6938 cMap.mapMaster
6939 (
6940 coupleMap::FACE,
6941 sfIndex,
6942 mfIndex
6943 );
6945 matched = true;
6947 break;
6948 }
6949 }
6951 if (!matched)
6952 {
6953 unMatchedFaces.insert(mfIndex);
6955 continue;
6956 }
6957 }
6958 else
6959 {
6960 label sFaceIndex = -1;
6962 forAll(spFaces, faceJ)
6963 {
6964 // Fetch the global face index
6965 label sfIndex = spFaces[faceJ];
6967 const face& sFace = slaveFaces[sfIndex];
6969 if
6970 (
6971 triFace::compare
6972 (
6973 triFace(cFace[0], cFace[1], cFace[2]),
6974 triFace(sFace[0], sFace[1], sFace[2])
6975 )
6976 )
6977 {
6978 // Found the slave. Add a map entry
6979 cMap.mapSlave
6980 (
6981 coupleMap::FACE,
6982 mfIndex,
6983 sfIndex
6984 );
6986 cMap.mapMaster
6987 (
6988 coupleMap::FACE,
6989 sfIndex,
6990 mfIndex
6991 );
6993 sFaceIndex = sfIndex;
6995 break;
6996 }
6997 }
6999 if (sFaceIndex == -1)
7000 {
7001 unMatchedFaces.insert(mfIndex);
7003 continue;
7004 }
7006 // Match all edges on this face as well.
7007 const labelList& mfEdges = faceEdges_[mfIndex];
7008 const labelList& sfEdges = sMesh.faceEdges_[sFaceIndex];
7010 forAll(mfEdges, edgeI)
7011 {
7012 if (eMap.found(mfEdges[edgeI]))
7013 {
7014 continue;
7015 }
7017 const edge& mEdge = edges_[mfEdges[edgeI]];
7019 // Configure a comparison edge.
7020 edge cEdge(pMap[mEdge[0]], pMap[mEdge[1]]);
7022 bool matchedEdge = false;
7024 forAll(sfEdges, edgeJ)
7025 {
7026 const edge& sEdge =
7027 (
7028 sMesh.edges_[sfEdges[edgeJ]]
7029 );
7031 if (cEdge == sEdge)
7032 {
7033 // Found the slave. Add a map entry
7034 cMap.mapSlave
7035 (
7036 coupleMap::EDGE,
7037 mfEdges[edgeI],
7038 sfEdges[edgeJ]
7039 );
7041 cMap.mapMaster
7042 (
7043 coupleMap::EDGE,
7044 sfEdges[edgeJ],
7045 mfEdges[edgeI]
7046 );
7048 matchedEdge = true;
7050 break;
7051 }
7052 }
7054 if (!matchedEdge)
7055 {
7056 // Write out the edge
7057 writeVTK("mEdge", mfEdges[edgeI], 1);
7059 FatalErrorIn
7060 (
7061 "void dynamicTopoFvMesh::"
7062 "buildProcessorCoupledMaps()"
7063 )
7064 << " Failed to match edge: "
7065 << mfEdges[edgeI] << ": "
7066 << mEdge << nl
7067 << " cEdge: " << cEdge
7068 << " for processor: " << proc
7069 << abort(FatalError);
7070 }
7071 }
7072 }
7073 }
7074 }
7076 // Attempt to match other edges in 3D, provided any common ones exist.
7077 // - This will handle point-only coupling situations for edges.
7078 if (!twoDMesh_)
7079 {
7080 // Fetch Maps for this processor
7081 const Map<label>& edgeMap = cMap.entityMap(coupleMap::EDGE);
7082 const Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
7084 forAllConstIter(Map<label>, pointMap, pIter)
7085 {
7086 const labelList& pEdges = pointEdges_[pIter.key()];
7088 forAll(pEdges, edgeI)
7089 {
7090 const label eIndex = pEdges[edgeI];
7092 // Only pick boundary edges
7093 if (whichEdgePatch(eIndex) == -1)
7094 {
7095 continue;
7096 }
7098 // Skip mapped edges
7099 if (edgeMap.found(eIndex))
7100 {
7101 continue;
7102 }
7104 const edge& eCheck = edges_[eIndex];
7106 // Check if a map exists for the other point
7107 label sIndex =
7108 (
7109 cMap.findSlave
7110 (
7111 coupleMap::POINT,
7112 eCheck.otherVertex(pIter.key())
7113 )
7114 );
7116 if (sIndex > -1)
7117 {
7118 const dynamicTopoFvMesh& mesh = rPM.subMesh();
7120 // Configure a check edge
7121 edge cEdge(sIndex, pIter());
7123 const labelList& spEdges = mesh.pointEdges_[sIndex];
7125 forAll(spEdges, edgeJ)
7126 {
7127 const edge& sEdge = mesh.edges_[spEdges[edgeJ]];
7129 if (sEdge == cEdge)
7130 {
7131 // Found the slave. Add a map entry
7132 cMap.mapSlave
7133 (
7134 coupleMap::EDGE,
7135 eIndex,
7136 spEdges[edgeJ]
7137 );
7139 cMap.mapMaster
7140 (
7141 coupleMap::EDGE,
7142 spEdges[edgeJ],
7143 eIndex
7144 );
7146 break;
7147 }
7148 }
7149 }
7150 }
7151 }
7152 }
7154 if (unMatchedFaces.size())
7155 {
7156 // Write out the face
7157 writeVTK("mFaces", unMatchedFaces.toc(), 2);
7159 FatalErrorIn
7160 (
7161 "void dynamicTopoFvMesh::buildProcessorCoupledMaps()"
7162 )
7163 << " Unmatched faces were found for processor: " << proc
7164 << abort(FatalError);
7165 }
7166 }
7167 }
7170 // Introduce a new processor patch to the mesh
7171 label dynamicTopoFvMesh::createProcessorPatch(const label proc)
7172 {
7173 if (!Pstream::parRun())
7174 {
7175 return -2;
7176 }
7178 // Get the new patch index,
7179 // and increment the number of patches
7180 label patchID = nPatches_++;
7182 // Find index in list of processors
7183 label pI = findIndex(procIndices_, proc);
7185 // Set the new patch index in patchMaps
7186 if (isSubMesh_)
7187 {
7188 if (pI == -1)
7189 {
7190 pI = procIndices_.size();
7192 procIndices_.setSize(pI + 1);
7193 sendMeshes_.setSize(pI + 1);
7194 }
7196 // Create a basic entry
7197 procIndices_[pI] = proc;
7199 sendMeshes_.set
7200 (
7201 pI,
7202 new coupledInfo
7203 (
7204 *this, // Reference to this mesh
7205 twoDMesh_, // 2D or 3D
7206 false, // Not local
7207 true, // Sent to neighbour
7208 patchID, // Patch index
7209 proc, // Master index
7210 -1 // Slave index
7211 )
7212 );
7213 }
7214 else
7215 {
7216 if (pI == -1)
7217 {
7218 Pout<< " Could not find index for processor: " << proc
7219 << " in indices: " << procIndices_
7220 << abort(FatalError);
7221 }
7223 sendMeshes_[pI].map().patchIndex() = patchID;
7224 recvMeshes_[pI].map().patchIndex() = patchID;
7225 }
7227 // Size up patches, and copy old information
7228 label prevPatchID = patchID - 1;
7230 patchSizes_.setSize(nPatches_, 0);
7231 oldPatchSizes_.setSize(nPatches_, 0);
7233 patchStarts_.setSize
7234 (
7235 nPatches_,
7236 patchStarts_[prevPatchID] + patchSizes_[prevPatchID]
7237 );
7238 oldPatchStarts_.setSize
7239 (
7240 nPatches_,
7241 oldPatchStarts_[prevPatchID] + oldPatchSizes_[prevPatchID]
7242 );
7244 edgePatchSizes_.setSize(nPatches_, 0);
7245 oldEdgePatchSizes_.setSize(nPatches_, 0);
7247 edgePatchStarts_.setSize
7248 (
7249 nPatches_,
7250 edgePatchStarts_[prevPatchID] + edgePatchSizes_[prevPatchID]
7251 );
7252 oldEdgePatchStarts_.setSize
7253 (
7254 nPatches_,
7255 oldEdgePatchStarts_[prevPatchID] + oldEdgePatchSizes_[prevPatchID]
7256 );
7258 patchNMeshPoints_.setSize(nPatches_, 0);
7259 oldPatchNMeshPoints_.setSize(nPatches_, 0);
7261 if (debug)
7262 {
7263 Pout<< " dynamicTopoFvMesh::createProcessorPatch :"
7264 << " Created new patch for processor: " << proc << nl
7265 << " On subMesh: " << Switch::asText(isSubMesh_) << nl
7266 << " pI: " << pI << nl
7267 << " patchID: " << patchID << nl
7268 << " patchStarts: " << patchStarts_ << nl
7269 << " patchSizes: " << patchSizes_ << nl
7270 << endl;
7271 }
7273 // Return the new patch index
7274 return patchID;
7275 }
7278 // If a patch is of processor type, get the neighbouring processor ID
7279 label dynamicTopoFvMesh::getNeighbourProcessor(const label patch) const
7280 {
7281 if (!Pstream::parRun())
7282 {
7283 return -1;
7284 }
7286 label neiProcNo = -1;
7288 const polyBoundaryMesh& boundary = boundaryMesh();
7290 if (patch == -1)
7291 {
7292 return -1;
7293 }
7294 else
7295 if (patch < boundary.size())
7296 {
7297 if (isA<processorPolyPatch>(boundary[patch]))
7298 {
7299 const processorPolyPatch& pp =
7300 (
7301 refCast<const processorPolyPatch>(boundary[patch])
7302 );
7304 // Set the neighbour processor ID
7305 neiProcNo = pp.neighbProcNo();
7306 }
7307 }
7308 else
7309 {
7310 // New processor patch
7312 // Find the neighbour processor ID
7313 forAll(procIndices_, pI)
7314 {
7315 label patchID = sendMeshes_[pI].map().patchIndex();
7317 if (patch == patchID)
7318 {
7319 neiProcNo = procIndices_[pI];
7320 break;
7321 }
7322 }
7324 if (neiProcNo == -1)
7325 {
7326 // An index should have been defined by now
7327 Pout<< " Patch: " << patch
7328 << " was not defined on patch subMeshes."
7329 << abort(FatalError);
7330 }
7331 }
7333 return neiProcNo;
7334 }
7337 // If the number of patches have changed during run-time,
7338 // reset boundaries with new processor patches
7339 void dynamicTopoFvMesh::resetBoundaries()
7340 {
7341 // Prepare a new list of patches
7342 List<polyPatch*> patches(nPatches_);
7344 // Fetch reference to existing boundary
7345 // - The removeBoundary member merely resets
7346 // boundary size, so this reference is safe
7347 const polyBoundaryMesh& boundary = boundaryMesh();
7349 // Copy all existing patches first
7350 for (label patchI = 0; patchI < boundaryMesh().size(); patchI++)
7351 {
7352 // Clone the patch
7353 patches[patchI] = boundary[patchI].clone(boundary).ptr();
7354 }
7356 // Create new processor patches
7357 for (label patchI = boundaryMesh().size(); patchI < nPatches_; patchI++)
7358 {
7359 // Make a temporary dictionary for patch construction
7360 dictionary patchDict;
7362 // Back out the neighbour processor ID
7363 label neiProcNo = getNeighbourProcessor(patchI);
7365 // Add relevant info
7366 patchDict.add("type", "processor");
7367 patchDict.add("startFace", patchStarts_[patchI]);
7368 patchDict.add("nFaces", patchSizes_[patchI]);
7369 patchDict.add("myProcNo", Pstream::myProcNo());
7370 patchDict.add("neighbProcNo", neiProcNo);
7372 // Set the pointer
7373 patches[patchI] =
7374 (
7375 polyPatch::New
7376 (
7377 "procBoundary"
7378 + Foam::name(Pstream::myProcNo())
7379 + "to"
7380 + Foam::name(neiProcNo),
7381 patchDict,
7382 patchI,
7383 boundary
7384 ).ptr()
7385 );
7386 }
7388 // Remove the old boundary
7389 fvMesh::removeFvBoundary();
7391 // Add patches, but don't calculate geometry, etc
7392 fvMesh::addFvPatches(patches, false);
7393 }
7396 // Initialize subMesh field transfers for mapping
7397 void dynamicTopoFvMesh::initFieldTransfers
7398 (
7399 wordList& types,
7400 List<wordList>& names,
7401 List<List<char> >& sendBuffer,
7402 List<List<char> >& recvBuffer
7403 )
7404 {
7405 if (!Pstream::parRun())
7406 {
7407 return;
7408 }
7410 if (debug)
7411 {
7412 Info<< " void dynamicTopoFvMesh::initFieldTransfers() :"
7413 << " Initializing subMesh field transfers."
7414 << endl;
7415 }
7417 // Clear out mesh geometry, since those
7418 // are to be wiped out after topo-changes anyway.
7419 fvMesh::clearOut();
7420 polyMesh::resetMotion();
7422 // Size up wordLists
7423 // - Five templated volFields
7424 // - Five templated surfaceFields
7425 // - One scalar volume gradient (conservative mapping)
7426 // - One vector volume gradient (conservative mapping)
7427 names.setSize(12);
7428 types.setSize(12);
7430 // Fill in field-types
7431 types[0] = volScalarField::typeName;
7432 types[1] = volVectorField::typeName;
7433 types[2] = volSphericalTensorField::typeName;
7434 types[3] = volSymmTensorField::typeName;
7435 types[4] = volTensorField::typeName;
7437 types[5] = surfaceScalarField::typeName;
7438 types[6] = surfaceVectorField::typeName;
7439 types[7] = surfaceSphericalTensorField::typeName;
7440 types[8] = surfaceSymmTensorField::typeName;
7441 types[9] = surfaceTensorField::typeName;
7443 types[10] = "grad(" + volScalarField::typeName + ')';
7444 types[11] = "grad(" + volVectorField::typeName + ')';
7446 // Send / recv buffers for field names
7447 List<char> fieldNameSendBuffer, fieldNameRecvBuffer;
7449 // Fetch reference to mapper
7450 const topoMapper& fieldMapper = mapper_();
7452 if (Pstream::master())
7453 {
7454 PtrList<OStringStream> fieldNameStream(1);
7456 // Set the stream
7457 fieldNameStream.set(0, new OStringStream(IOstream::BINARY));
7459 // Alias for convenience
7460 OStringStream& fNStream = fieldNameStream[0];
7462 // Fetch field-names by type
7463 for (label typeI = 0; typeI < 10; typeI++)
7464 {
7465 // Get all fields of type
7466 names[typeI] = objectRegistry::names(types[typeI]);
7467 }
7469 // Fetch scalar gradient names
7470 names[10] = fieldMapper.scalarGrads();
7472 // Fetch vector gradient names
7473 names[11] = fieldMapper.vectorGrads();
7475 // Send field names to Ostream
7476 fNStream << names;
7478 // Size up buffers and fill contents
7479 string contents = fNStream.str();
7480 const char* ptr = contents.data();
7482 fieldNameSendBuffer.setSize(contents.size());
7484 forAll(fieldNameSendBuffer, i)
7485 {
7486 fieldNameSendBuffer[i] = *ptr++;
7487 }
7489 // Clear the stream
7490 fieldNameStream.set(0, NULL);
7492 if (debug > 4)
7493 {
7494 Info<< " Registered fields: " << names << endl;
7495 }
7497 // Send names to slaves
7498 for (label proc = 1; proc < Pstream::nProcs(); proc++)
7499 {
7500 // Send buffer to processor
7501 meshOps::pWrite(proc, fieldNameSendBuffer.size());
7502 meshOps::pWrite(proc, fieldNameSendBuffer);
7503 }
7504 }
7505 else
7506 {
7507 // Receive names from master
7508 label recvBufferSize = -1;
7509 meshOps::pRead(Pstream::masterNo(), recvBufferSize);
7511 // Size up buffer and schedule receive
7512 fieldNameRecvBuffer.setSize(recvBufferSize);
7513 meshOps::pRead(Pstream::masterNo(), fieldNameRecvBuffer);
7514 }
7516 // Wait for transfers to complete
7517 meshOps::waitForBuffers();
7519 // Convert names from stringStream
7520 if (!Pstream::master())
7521 {
7522 string contents
7523 (
7524 fieldNameRecvBuffer.begin(),
7525 fieldNameRecvBuffer.size()
7526 );
7528 // Set the stream
7529 IStringStream fieldNameStream(contents, IOstream::BINARY);
7531 // Get names from stream
7532 fieldNameStream >> names;
7533 }
7535 label nProcs = procIndices_.size();
7537 // Size up buffers
7538 sendBuffer.setSize(nProcs);
7539 recvBuffer.setSize(nProcs);
7541 // Size up the send stringStream
7542 PtrList<OStringStream> stream(nProcs);
7544 // Now fill in subMesh fields
7545 forAll(procIndices_, pI)
7546 {
7547 label proc = procIndices_[pI];
7549 const coupledInfo& cInfo = sendMeshes_[pI];
7551 // Initialize stream
7552 stream.set(pI, new OStringStream(IOstream::BINARY));
7554 // Subset volFields to stream
7555 cInfo.mapVolField<scalar>(names[0], types[0], stream[pI]);
7556 cInfo.mapVolField<vector>(names[1], types[1], stream[pI]);
7557 cInfo.mapVolField<sphericalTensor>(names[2], types[2], stream[pI]);
7558 cInfo.mapVolField<symmTensor>(names[3], types[3], stream[pI]);
7559 cInfo.mapVolField<tensor>(names[4], types[4], stream[pI]);
7561 // Subset surfaceFields to stream
7562 cInfo.mapSurfaceField<scalar>(names[5], types[5], stream[pI]);
7563 cInfo.mapSurfaceField<vector>(names[6], types[6], stream[pI]);
7564 cInfo.mapSurfaceField<sphericalTensor>(names[7], types[7], stream[pI]);
7565 cInfo.mapSurfaceField<symmTensor>(names[8], types[8], stream[pI]);
7566 cInfo.mapSurfaceField<tensor>(names[9], types[9], stream[pI]);
7568 // Subset scalar gradients to stream
7569 stream[pI]
7570 << types[10] << token::NL
7571 << token::BEGIN_BLOCK << token::NL;
7573 forAll(names[10], i)
7574 {
7575 tmp<volVectorField> tvvfFld =
7576 (
7577 cInfo.subSetVolField
7578 (
7579 fieldMapper.gradient<volVectorField>
7580 (
7581 names[10][i]
7582 )
7583 )
7584 );
7586 // Send field through stream
7587 stream[pI]
7588 << names[10][i]
7589 << token::NL << token::BEGIN_BLOCK
7590 << tvvfFld
7591 << token::NL << token::END_BLOCK
7592 << token::NL;
7593 }
7595 stream[pI]
7596 << token::END_BLOCK << token::NL;
7598 // Subset vector gradients to stream
7599 stream[pI]
7600 << types[11] << token::NL
7601 << token::BEGIN_BLOCK << token::NL;
7603 forAll(names[11], i)
7604 {
7605 tmp<volTensorField> tvtfFld =
7606 (
7607 cInfo.subSetVolField
7608 (
7609 fieldMapper.gradient<volTensorField>
7610 (
7611 names[11][i]
7612 )
7613 )
7614 );
7616 // Send field through stream
7617 stream[pI]
7618 << names[11][i]
7619 << token::NL << token::BEGIN_BLOCK
7620 << tvtfFld
7621 << token::NL << token::END_BLOCK
7622 << token::NL;
7623 }
7625 stream[pI]
7626 << token::END_BLOCK << token::NL;
7628 // Size up buffers and fill contents
7629 string contents = stream[pI].str();
7630 const char* ptr = contents.data();
7632 sendBuffer[pI].setSize(contents.size());
7634 forAll(sendBuffer[pI], i)
7635 {
7636 sendBuffer[pI][i] = *ptr++;
7637 }
7639 // Clear the stream
7640 stream.set(pI, NULL);
7642 // Send buffer to processor
7643 meshOps::pWrite(proc, sendBuffer[pI].size());
7644 meshOps::pWrite(proc, sendBuffer[pI]);
7646 // Receive buffer from processor
7647 label recvBufferSize = -1;
7648 meshOps::pRead(proc, recvBufferSize);
7650 // Size up buffer and schedule receive
7651 recvBuffer[pI].setSize(recvBufferSize);
7652 meshOps::pRead(proc, recvBuffer[pI]);
7653 }
7655 // We won't wait for all transfers to complete for the moment.
7656 // Meanwhile, do some other useful work, if possible.
7657 }
7660 // Synchronize field transfers for mapping
7661 void dynamicTopoFvMesh::syncFieldTransfers
7662 (
7663 wordList& types,
7664 List<wordList>& names,
7665 List<List<char> >& recvBuffer
7666 )
7667 {
7668 if (!Pstream::parRun())
7669 {
7670 return;
7671 }
7673 if (debug)
7674 {
7675 Info<< " void dynamicTopoFvMesh::syncFieldTransfers() :"
7676 << " Synchronizing subMesh field transfers."
7677 << endl;
7678 }
7680 // Wait for all transfers to complete.
7681 meshOps::waitForBuffers();
7683 label nProcs = procIndices_.size();
7685 // Size up stringStream
7686 PtrList<IStringStream> stream(nProcs);
7688 // Size up field PtrLists
7690 // Volume Fields
7691 List<PtrList<volScalarField> > vsF(nProcs);
7692 List<PtrList<volVectorField> > vvF(nProcs);
7693 List<PtrList<volSphericalTensorField> > vsptF(nProcs);
7694 List<PtrList<volSymmTensorField> > vsytF(nProcs);
7695 List<PtrList<volTensorField> > vtF(nProcs);
7697 // Surface fields
7698 List<PtrList<surfaceScalarField> > ssF(nProcs);
7699 List<PtrList<surfaceVectorField> > svF(nProcs);
7700 List<PtrList<surfaceSphericalTensorField> > ssptF(nProcs);
7701 List<PtrList<surfaceSymmTensorField> > ssytF(nProcs);
7702 List<PtrList<surfaceTensorField> > stF(nProcs);
7704 // Gradients
7705 List<PtrList<volVectorField> > vgsF(nProcs);
7706 List<PtrList<volTensorField> > vgvF(nProcs);
7708 // Geometry
7709 List<PtrList<volVectorField> > vC(nProcs);
7710 List<PtrList<surfaceScalarField> > sSf(nProcs);
7712 const polyBoundaryMesh& boundary = boundaryMesh();
7714 // Determine number of physical (non-processor) patches
7715 label nPhysical = 0;
7717 forAll(boundary, patchI)
7718 {
7719 if (isA<processorPolyPatch>(boundary[patchI]))
7720 {
7721 continue;
7722 }
7724 nPhysical++;
7725 }
7727 // Keep track of extra / total entities
7728 label nTotalCells = nOldCells_, nTotalIntFaces = nOldInternalFaces_;
7729 labelList nTotalPatchFaces(SubList<label>(oldPatchSizes_, nPhysical));
7731 // Allocate reverse maps
7732 List<labelList> irvMaps(nProcs), irsMaps(nProcs);
7733 List<labelListList> brMaps(nProcs, labelListList(nPhysical));
7735 forAll(procIndices_, pI)
7736 {
7737 const coupledInfo& cInfo = recvMeshes_[pI];
7739 // Convert buffer to string
7740 string contents(recvBuffer[pI].begin(), recvBuffer[pI].size());
7742 // Initialize stream
7743 stream.set(pI, new IStringStream(contents, IOstream::BINARY));
7745 // Construct dictionary from stream
7746 dictionary dict(stream[pI]);
7748 // Set field pointers
7749 cInfo.setField(names[0], dict.subDict(types[0]), vsF[pI]);
7750 cInfo.setField(names[1], dict.subDict(types[1]), vvF[pI]);
7751 cInfo.setField(names[2], dict.subDict(types[2]), vsptF[pI]);
7752 cInfo.setField(names[3], dict.subDict(types[3]), vsytF[pI]);
7753 cInfo.setField(names[4], dict.subDict(types[4]), vtF[pI]);
7755 cInfo.setField(names[5], dict.subDict(types[5]), ssF[pI]);
7756 cInfo.setField(names[6], dict.subDict(types[6]), svF[pI]);
7757 cInfo.setField(names[7], dict.subDict(types[7]), ssptF[pI]);
7758 cInfo.setField(names[8], dict.subDict(types[8]), ssytF[pI]);
7759 cInfo.setField(names[9], dict.subDict(types[9]), stF[pI]);
7761 cInfo.setField(names[10], dict.subDict(types[10]), vgsF[pI]);
7762 cInfo.setField(names[11], dict.subDict(types[11]), vgvF[pI]);
7764 // Set centres
7765 cInfo.setCentres(vC[pI]);
7767 // Count the number of additional entities
7768 const coupleMap& cMap = cInfo.map();
7770 // Fetch size from subMesh
7771 label nCells = cMap.nEntities(coupleMap::CELL);
7772 label nIntFaces = cMap.nEntities(coupleMap::INTERNAL_FACE);
7774 // Set rmap for this processor
7775 irvMaps[pI] = (labelField(identity(nCells)) + nTotalCells);
7776 irsMaps[pI] = (labelField(identity(nIntFaces)) + nTotalIntFaces);
7778 // Update count
7779 nTotalCells += nCells;
7780 nTotalIntFaces += nIntFaces;
7782 // Fetch patch sizes from subMesh
7783 const labelList& nPatchFaces =
7784 (
7785 cMap.entityBuffer(coupleMap::FACE_SIZES)
7786 );
7788 // Loop over physical patches
7789 forAll(nTotalPatchFaces, patchI)
7790 {
7791 // Fetch patch size from subMesh
7792 label nFaces = nPatchFaces[patchI];
7794 // Set rmap for this patch
7795 brMaps[pI][patchI] =
7796 (
7797 labelField(identity(nFaces))
7798 + nTotalPatchFaces[patchI]
7799 );
7801 // Update patch-face count
7802 nTotalPatchFaces[patchI] += nFaces;
7803 }
7804 }
7806 // Prepare internal mappers
7807 labelList cellAddressing(nTotalCells, 0);
7808 labelList faceAddressing(nTotalIntFaces, 0);
7810 // Set identity map for first nCells,
7811 // and map from cell[0] for the rest
7812 for (label i = 0; i < nOldCells_; i++)
7813 {
7814 cellAddressing[i] = i;
7815 }
7817 // Set identity map for first nIntFaces,
7818 // and map from face[0] for the rest
7819 for (label i = 0; i < nOldInternalFaces_; i++)
7820 {
7821 faceAddressing[i] = i;
7822 }
7824 coupledInfo::subMeshMapper vMap
7825 (
7826 nOldCells_,
7827 cellAddressing
7828 );
7830 coupledInfo::subMeshMapper sMap
7831 (
7832 nOldInternalFaces_,
7833 faceAddressing
7834 );
7836 // Prepare boundary mappers
7837 labelListList patchAddressing(nPhysical);
7838 PtrList<coupledInfo::subMeshMapper> bMap(nPhysical);
7840 forAll(bMap, patchI)
7841 {
7842 // Prepare patch mappers
7843 patchAddressing[patchI].setSize(nTotalPatchFaces[patchI], 0);
7845 // Set identity map for first nPatchFaces,
7846 // and map from patch-face[0] for the rest
7847 for (label i = 0; i < oldPatchSizes_[patchI]; i++)
7848 {
7849 patchAddressing[patchI][i] = i;
7850 }
7852 // Set the boundary mapper pointer
7853 bMap.set
7854 (
7855 patchI,
7856 new coupledInfo::subMeshMapper
7857 (
7858 oldPatchSizes_[patchI],
7859 patchAddressing[patchI]
7860 )
7861 );
7862 }
7864 // Loop through all volFields and re-size
7865 // to accomodate additional cells / faces
7866 coupledInfo::resizeMap(names[0], *this, vMap, irvMaps, bMap, brMaps, vsF);
7867 coupledInfo::resizeMap(names[1], *this, vMap, irvMaps, bMap, brMaps, vvF);
7868 coupledInfo::resizeMap(names[2], *this, vMap, irvMaps, bMap, brMaps, vsptF);
7869 coupledInfo::resizeMap(names[3], *this, vMap, irvMaps, bMap, brMaps, vsytF);
7870 coupledInfo::resizeMap(names[4], *this, vMap, irvMaps, bMap, brMaps, vtF);
7872 coupledInfo::resizeMap(names[5], *this, sMap, irsMaps, bMap, brMaps, ssF);
7873 coupledInfo::resizeMap(names[6], *this, sMap, irsMaps, bMap, brMaps, svF);
7874 coupledInfo::resizeMap(names[7], *this, sMap, irsMaps, bMap, brMaps, ssptF);
7875 coupledInfo::resizeMap(names[8], *this, sMap, irsMaps, bMap, brMaps, ssytF);
7876 coupledInfo::resizeMap(names[9], *this, sMap, irsMaps, bMap, brMaps, stF);
7878 // Fetch reference to mapper
7879 const topoMapper& fieldMapper = mapper_();
7881 // Map gradient fields
7882 forAll(names[10], i)
7883 {
7884 volVectorField& sGrad =
7885 (
7886 fieldMapper.gradient<volVectorField>(names[10][i])
7887 );
7889 // Map the field
7890 coupledInfo::resizeMap(i, vMap, irvMaps, bMap, brMaps, vgsF, sGrad);
7891 }
7893 forAll(names[11], i)
7894 {
7895 volTensorField& vGrad =
7896 (
7897 fieldMapper.gradient<volTensorField>(names[11][i])
7898 );
7900 // Map the field
7901 coupledInfo::resizeMap(i, vMap, irvMaps, bMap, brMaps, vgvF, vGrad);
7902 }
7904 // Map stored geometry
7905 volVectorField& mapC = fieldMapper.volCentres();
7907 // Map geometry fields
7908 coupledInfo::resizeMap(0, vMap, irvMaps, bMap, brMaps, vC, mapC);
7909 }
7912 // Initialize coupled boundary ordering
7913 // - Assumes that faces_ and points_ are consistent
7914 // - Assumes that patchStarts_ and patchSizes_ are consistent
7915 void dynamicTopoFvMesh::initCoupledBoundaryOrdering
7916 (
7917 List<pointField>& centres,
7918 List<pointField>& anchors
7919 ) const
7920 {
7921 if (!Pstream::parRun())
7922 {
7923 return;
7924 }
7926 for (label pI = 0; pI < nPatches_; pI++)
7927 {
7928 label neiProcNo = getNeighbourProcessor(pI);
7930 if (neiProcNo == -1)
7931 {
7932 continue;
7933 }
7935 // Check if this is a master processor patch.
7936 label start = patchStarts_[pI];
7937 label size = patchSizes_[pI];
7939 // Prepare centres and anchors
7940 centres[pI].setSize(size, vector::zero);
7941 anchors[pI].setSize(size, vector::zero);
7943 if (Pstream::myProcNo() < neiProcNo)
7944 {
7945 forAll(centres[pI], fI)
7946 {
7947 centres[pI][fI] = faces_[fI + start].centre(points_);
7948 anchors[pI][fI] = points_[faces_[fI + start][0]];
7949 }
7951 {
7952 if (debug > 3)
7953 {
7954 Pout<< "Sending to: " << neiProcNo
7955 << " nCentres: " << size << endl;
7957 // Write out my centres to disk
7958 meshOps::writeVTK
7959 (
7960 (*this),
7961 "centres_"
7962 + Foam::name(Pstream::myProcNo())
7963 + "to"
7964 + Foam::name(neiProcNo),
7965 size, size, size,
7966 centres[pI]
7967 );
7968 }
7970 // Ensure that we're sending the right size
7971 meshOps::pWrite(neiProcNo, size);
7972 }
7974 // Send information to neighbour
7975 meshOps::pWrite(neiProcNo, centres[pI]);
7976 meshOps::pWrite(neiProcNo, anchors[pI]);
7977 }
7978 else
7979 {
7980 {
7981 label nEntities = -1;
7983 // Ensure that we're receiving the right size
7984 meshOps::pRead(neiProcNo, nEntities);
7986 if (debug > 3)
7987 {
7988 Pout<< " Recving from: " << neiProcNo
7989 << " nCentres: " << size << endl;
7990 }
7992 if (nEntities != size)
7993 {
7994 // Set the size to what we're receiving
7995 centres[pI].setSize(nEntities, vector::zero);
7996 anchors[pI].setSize(nEntities, vector::zero);
7998 // Issue a warning now, and let
7999 // syncCoupledBoundaryOrdering fail later
8000 WarningIn
8001 (
8002 "void dynamicTopoFvMesh::"
8003 "initCoupledBoundaryOrdering() const"
8004 )
8005 << "Incorrect send / recv sizes: " << nl
8006 << " nEntities: " << nEntities << nl
8007 << " size: " << size << nl
8008 << endl;
8009 }
8010 }
8012 // Schedule receive from neighbour
8013 meshOps::pRead(neiProcNo, centres[pI]);
8014 meshOps::pRead(neiProcNo, anchors[pI]);
8015 }
8016 }
8017 }
8020 // Synchronize coupled boundary ordering
8021 bool dynamicTopoFvMesh::syncCoupledBoundaryOrdering
8022 (
8023 List<pointField>& centres,
8024 List<pointField>& anchors,
8025 labelListList& patchMaps,
8026 labelListList& rotations
8027 ) const
8028 {
8029 bool anyChange = false, failedPatchMatch = false;
8031 // Fetch tolerance
8032 scalar matchTol = Foam::debug::tolerances("patchFaceMatchTol", 1e-4);
8034 // Calculate centres and tolerances for any slave patches
8035 List<scalarField> slaveTols(nPatches_);
8036 List<pointField> slaveCentres(nPatches_);
8038 // Handle locally coupled patches
8039 forAll(patchCoupling_, pI)
8040 {
8041 if (patchCoupling_(pI))
8042 {
8043 const coupleMap& cMap = patchCoupling_[pI].map();
8045 label masterPatch = cMap.masterIndex();
8046 label slavePatch = cMap.slaveIndex();
8048 label mSize = patchSizes_[masterPatch];
8049 label sSize = patchSizes_[slavePatch];
8051 label mStart = patchStarts_[masterPatch];
8052 label sStart = patchStarts_[slavePatch];
8054 if (mSize != sSize)
8055 {
8056 Pout<< " Patch size mismatch: " << nl
8057 << " mSize: " << mSize << nl
8058 << " sSize: " << sSize << nl
8059 << " mStart: " << mStart << nl
8060 << " sStart: " << sStart << nl
8061 << " master: " << masterPatch << nl
8062 << " slave: " << slavePatch << nl
8063 << abort(FatalError);
8064 }
8066 // Calculate centres / tolerances
8067 anchors[masterPatch].setSize(mSize, vector::zero);
8068 centres[masterPatch].setSize(mSize, vector::zero);
8070 slaveTols[slavePatch].setSize(sSize, 0.0);
8071 slaveCentres[slavePatch].setSize(sSize, vector::zero);
8073 for (label fI = 0; fI < mSize; fI++)
8074 {
8075 point& mfa = anchors[masterPatch][fI];
8076 point& mfc = centres[masterPatch][fI];
8077 point& sfc = slaveCentres[slavePatch][fI];
8079 const face& mFace = faces_[fI + mStart];
8080 const face& sFace = faces_[fI + sStart];
8082 // Calculate anchor / centres
8083 mfa = points_[mFace[0]];
8084 mfc = mFace.centre(points_);
8085 sfc = sFace.centre(points_);
8087 scalar maxLen = -GREAT;
8089 forAll(sFace, fpI)
8090 {
8091 maxLen = max(maxLen, mag(points_[sFace[fpI]] - sfc));
8092 }
8094 slaveTols[slavePatch][fI] = matchTol*maxLen;
8095 }
8097 // For cyclics, additionally test for halves,
8098 // and transform points as necessary.
8099 labelList halfMap;
8101 if (masterPatch == slavePatch)
8102 {
8103 scalar featureCos = 0.9;
8104 label halfSize = (mSize / 2);
8105 label half0Index = 0, half1Index = 0;
8107 // Size up halfMap and slaveAnchors
8108 halfMap.setSize(mSize, -1);
8110 scalarField half1Tols(halfSize, 0.0);
8111 vectorField half1Anchors(halfSize, vector::zero);
8113 // Fetch face-normal for the first face
8114 vector fhN = boundaryMesh()[masterPatch].faceAreas()[0];
8116 // Normalize
8117 fhN /= mag(fhN) + VSMALL;
8119 // Fetch transform type
8120 const cyclicPolyPatch& cyclicPatch =
8121 (
8122 refCast<const cyclicPolyPatch>
8123 (
8124 boundaryMesh()[masterPatch]
8125 )
8126 );
8128 bool translate =
8129 (
8130 cyclicPatch.transform()
8131 == cyclicPolyPatch::TRANSLATIONAL
8132 );
8134 for (label fI = 0; fI < mSize; fI++)
8135 {
8136 const face& mFace = faces_[fI + mStart];
8138 // Compute normal
8139 vector fN = mFace.normal(points_);
8141 // Normalize
8142 fN /= mag(fN) + VSMALL;
8144 // Check which half this face belongs to...
8145 bool isFirstHalf = ((fhN & fN) > featureCos);
8147 if (isFirstHalf)
8148 {
8149 vector mA = anchors[masterPatch][fI];
8150 vector mC = centres[masterPatch][fI];
8152 // Set master anchor
8153 anchors[masterPatch][half0Index] = mA;
8155 // Transform master points
8156 if (translate)
8157 {
8158 mA += cyclicPatch.separationVector();
8159 mC += cyclicPatch.separationVector();
8160 }
8161 else
8162 {
8163 // Assume constant transform tensor
8164 mA = Foam::transform(cyclicPatch.transformT(0), mA);
8165 mC = Foam::transform(cyclicPatch.transformT(0), mC);
8166 }
8168 // Set the transformed anchor / centre
8169 half1Anchors[half0Index] = mA;
8170 centres[masterPatch][half0Index] = mC;
8171 half1Tols[half0Index] = slaveTols[masterPatch][fI];
8173 // Set halfMap
8174 halfMap[fI] = half0Index;
8176 half0Index++;
8177 }
8178 else
8179 {
8180 // Set the slave centre
8181 const scalar& sT = slaveTols[slavePatch][fI];
8182 const point& sC = slaveCentres[slavePatch][fI];
8184 // Duplicate slaveTols for first / second half
8185 slaveTols[slavePatch][half1Index] = sT;
8186 slaveCentres[slavePatch][half1Index] = sC;
8188 // Set halfMap
8189 halfMap[fI] = half1Index + halfSize;
8191 half1Index++;
8192 }
8193 }
8195 // Sanity check for halfSize
8196 if (half0Index != halfSize || half1Index != halfSize)
8197 {
8198 Pout<< " Master: " << masterPatch << nl
8199 << " Slave: " << slavePatch << nl
8200 << " half0Index: " << half0Index << nl
8201 << " half1Index: " << half1Index << nl
8202 << " halfSize: " << halfSize << nl
8203 << " failed to divide into halves."
8204 << abort(FatalError);
8205 }
8207 // Resize lists
8208 centres[masterPatch].setSize(halfSize);
8209 slaveCentres[slavePatch].setSize(halfSize);
8211 // Set slave tols / anchors
8212 forAll(half1Anchors, fI)
8213 {
8214 slaveTols[masterPatch][fI + halfSize] = half1Tols[fI];
8215 anchors[masterPatch][fI + halfSize] = half1Anchors[fI];
8216 }
8217 }
8219 // Fetch reference to slave map
8220 labelList& patchMap = patchMaps[slavePatch];
8222 // Try zero separation automatic matching
8223 bool matchedAll =
8224 (
8225 matchPoints
8226 (
8227 slaveCentres[slavePatch],
8228 centres[masterPatch],
8229 slaveTols[slavePatch],
8230 false,
8231 patchMap
8232 )
8233 );
8235 // Write out master centres to disk
8236 if (debug > 3 || !matchedAll)
8237 {
8238 meshOps::writeVTK
8239 (
8240 (*this),
8241 "localMasterCentres_"
8242 + Foam::name(masterPatch)
8243 + '_'
8244 + Foam::name(slavePatch),
8245 mSize, mSize, mSize,
8246 centres[masterPatch]
8247 );
8249 meshOps::writeVTK
8250 (
8251 (*this),
8252 "localSlaveCentres_"
8253 + Foam::name(masterPatch)
8254 + '_'
8255 + Foam::name(slavePatch),
8256 sSize, sSize, sSize,
8257 slaveCentres[slavePatch]
8258 );
8260 writeVTK
8261 (
8262 "localMasterFaces_"
8263 + Foam::name(masterPatch)
8264 + '_'
8265 + Foam::name(slavePatch),
8266 identity(mSize) + mStart,
8267 2, false, true
8268 );
8270 writeVTK
8271 (
8272 "localSlaveFaces_"
8273 + Foam::name(masterPatch)
8274 + '_'
8275 + Foam::name(slavePatch),
8276 identity(sSize) + sStart,
8277 2, false, true
8278 );
8280 if (!matchedAll)
8281 {
8282 Pout<< " Master: " << masterPatch
8283 << " Slave: " << slavePatch
8284 << " mSize: " << mSize << " sSize: " << sSize
8285 << " failed on match for face centres."
8286 << endl;
8288 // Failed on match-for-all, so patchMaps
8289 // will be invalid. Bail out for now.
8290 failedPatchMatch = true;
8292 continue;
8293 }
8294 }
8296 // Renumber patchMap for cyclics
8297 if (masterPatch == slavePatch)
8298 {
8299 label halfSize = (mSize / 2);
8301 forAll(halfMap, fI)
8302 {
8303 if (halfMap[fI] >= halfSize)
8304 {
8305 halfMap[fI] =
8306 (
8307 patchMap[halfMap[fI] - halfSize]
8308 + halfSize
8309 );
8310 }
8311 }
8313 // Transfer map
8314 patchMap.transfer(halfMap);
8315 }
8317 // Fetch reference to rotation map
8318 labelList& rotation = rotations[slavePatch];
8320 // Initialize rotation to zero
8321 rotation.setSize(sSize, 0);
8323 // Set rotation.
8324 forAll(patchMap, oldFaceI)
8325 {
8326 label newFaceI = patchMap[oldFaceI];
8328 const point& anchor = anchors[slavePatch][newFaceI];
8329 const scalar& faceTol = slaveTols[slavePatch][oldFaceI];
8330 const face& checkFace = faces_[sStart + oldFaceI];
8332 label anchorFp = -1;
8333 scalar minDSqr = GREAT;
8335 forAll(checkFace, fpI)
8336 {
8337 scalar dSqr = magSqr(anchor - points_[checkFace[fpI]]);
8339 if (dSqr < minDSqr)
8340 {
8341 minDSqr = dSqr;
8342 anchorFp = fpI;
8343 }
8344 }
8346 if (anchorFp == -1 || mag(minDSqr) > faceTol)
8347 {
8348 FatalErrorIn
8349 (
8350 "\n"
8351 "void dynamicTopoFvMesh::"
8352 "syncCoupledBoundaryOrdering\n"
8353 "(\n"
8354 " List<pointField>& centres,\n"
8355 " List<pointField>& anchors,\n"
8356 " labelListList& patchMaps,\n"
8357 " labelListList& rotations\n"
8358 ") const\n"
8359 )
8360 << "Cannot find anchor: " << anchor << nl
8361 << " Face: " << checkFace << nl
8362 << " Vertices: "
8363 << UIndirectList<point>(points_, checkFace) << nl
8364 << " on patch: " << slavePatch << nl
8365 << " faceTol: " << faceTol << nl
8366 << " newFaceI: " << newFaceI << nl
8367 << " oldFaceI: " << oldFaceI << nl
8368 << " mSize: " << mSize << nl
8369 << " sSize: " << sSize << nl
8370 << abort(FatalError);
8371 }
8372 else
8373 {
8374 // Positive rotation
8375 // - Set for old face. Will be rotated later
8376 // during the shuffling stage
8377 rotation[oldFaceI] =
8378 (
8379 (checkFace.size() - anchorFp) % checkFace.size()
8380 );
8381 }
8382 }
8384 // Set the flag
8385 anyChange = true;
8386 }
8387 }
8389 if (failedPatchMatch)
8390 {
8391 FatalErrorIn
8392 (
8393 "\n"
8394 "void dynamicTopoFvMesh::"
8395 "syncCoupledBoundaryOrdering\n"
8396 "(\n"
8397 " List<pointField>& centres,\n"
8398 " List<pointField>& anchors,\n"
8399 " labelListList& patchMaps,\n"
8400 " labelListList& rotations\n"
8401 ") const\n"
8402 )
8403 << " Matching for local patches failed. " << nl
8404 << abort(FatalError);
8405 }
8407 if (!Pstream::parRun())
8408 {
8409 return anyChange;
8410 }
8412 for (label pI = 0; pI < nPatches_; pI++)
8413 {
8414 label neiProcNo = getNeighbourProcessor(pI);
8416 if (neiProcNo == -1)
8417 {
8418 continue;
8419 }
8421 // Check if this is a slave processor patch.
8422 label start = patchStarts_[pI];
8423 label size = patchSizes_[pI];
8425 if (Pstream::myProcNo() > neiProcNo)
8426 {
8427 slaveTols[pI].setSize(size, 0.0);
8428 slaveCentres[pI].setSize(size, vector::zero);
8430 forAll(slaveCentres[pI], fI)
8431 {
8432 point& fc = slaveCentres[pI][fI];
8434 const face& checkFace = faces_[fI + start];
8436 // Calculate centre
8437 fc = checkFace.centre(points_);
8439 scalar maxLen = -GREAT;
8441 forAll(checkFace, fpI)
8442 {
8443 maxLen = max(maxLen, mag(points_[checkFace[fpI]] - fc));
8444 }
8446 slaveTols[pI][fI] = matchTol*maxLen;
8447 }
8449 // Write out my centres to disk
8450 if (debug > 3)
8451 {
8452 meshOps::writeVTK
8453 (
8454 (*this),
8455 "slaveCentres_"
8456 + Foam::name(Pstream::myProcNo())
8457 + '_'
8458 + Foam::name(neiProcNo),
8459 size, size, size,
8460 slaveCentres[pI]
8461 );
8462 }
8463 }
8464 }
8466 // Wait for transfers before continuing.
8467 meshOps::waitForBuffers();
8469 for (label pI = 0; pI < nPatches_; pI++)
8470 {
8471 label neiProcNo = getNeighbourProcessor(pI);
8473 if (neiProcNo == -1)
8474 {
8475 continue;
8476 }
8478 // Check if this is a master processor patch.
8479 labelList& patchMap = patchMaps[pI];
8480 labelList& rotation = rotations[pI];
8482 // Initialize map and rotation
8483 patchMap.setSize(patchSizes_[pI], -1);
8484 rotation.setSize(patchSizes_[pI], 0);
8486 if (Pstream::myProcNo() < neiProcNo)
8487 {
8488 // Do nothing (i.e. identical mapping, zero rotation).
8489 forAll(patchMap, pfI)
8490 {
8491 patchMap[pfI] = pfI;
8492 }
8493 }
8494 else
8495 {
8496 // Try zero separation automatic matching
8497 bool matchedAll =
8498 (
8499 matchPoints
8500 (
8501 slaveCentres[pI],
8502 centres[pI],
8503 slaveTols[pI],
8504 false,
8505 patchMap
8506 )
8507 );
8509 // Write out centres to disk
8510 if (debug > 3 || !matchedAll)
8511 {
8512 label mSize = centres[pI].size();
8513 label sSize = slaveCentres[pI].size();
8515 meshOps::writeVTK
8516 (
8517 (*this),
8518 "masterCentres_"
8519 + Foam::name(Pstream::myProcNo())
8520 + '_'
8521 + Foam::name(neiProcNo),
8522 mSize, mSize, mSize,
8523 centres[pI]
8524 );
8526 meshOps::writeVTK
8527 (
8528 (*this),
8529 "slaveCentres_"
8530 + Foam::name(Pstream::myProcNo())
8531 + '_'
8532 + Foam::name(neiProcNo),
8533 sSize, sSize, sSize,
8534 slaveCentres[pI]
8535 );
8537 writeVTK
8538 (
8539 "procPatchFaces_"
8540 + Foam::name(Pstream::myProcNo())
8541 + '_'
8542 + Foam::name(neiProcNo),
8543 identity(mSize) + patchStarts_[pI],
8544 2, false, true
8545 );
8547 if (!matchedAll)
8548 {
8549 Pout<< " Patch: " << pI
8550 << " Processor: " << neiProcNo
8551 << " mSize: " << mSize << " sSize: " << sSize
8552 << " failed on match for face centres."
8553 << endl;
8555 // Failed on match-for-all, so patchMaps
8556 // will be invalid. Bail out for now.
8557 failedPatchMatch = true;
8559 continue;
8560 }
8561 }
8563 label start = patchStarts_[pI];
8565 // Set rotation.
8566 forAll(patchMap, oldFaceI)
8567 {
8568 label newFaceI = patchMap[oldFaceI];
8570 const point& anchor = anchors[pI][newFaceI];
8571 const scalar& faceTol = slaveTols[pI][oldFaceI];
8572 const face& checkFace = faces_[start + oldFaceI];
8574 label anchorFp = -1;
8575 scalar minDSqr = GREAT;
8577 forAll(checkFace, fpI)
8578 {
8579 scalar dSqr = magSqr(anchor - points_[checkFace[fpI]]);
8581 if (dSqr < minDSqr)
8582 {
8583 minDSqr = dSqr;
8584 anchorFp = fpI;
8585 }
8586 }
8588 if (anchorFp == -1 || mag(minDSqr) > faceTol)
8589 {
8590 FatalErrorIn
8591 (
8592 "\n"
8593 "void dynamicTopoFvMesh::"
8594 "syncCoupledBoundaryOrdering\n"
8595 "(\n"
8596 " List<pointField>& centres,\n"
8597 " List<pointField>& anchors,\n"
8598 " labelListList& patchMaps,\n"
8599 " labelListList& rotations\n"
8600 ") const\n"
8601 )
8602 << "Cannot find anchor: " << anchor << nl
8603 << " Face: " << checkFace << nl
8604 << " Vertices: "
8605 << UIndirectList<point>(points_, checkFace) << nl
8606 << " on patch: " << pI
8607 << " for processor: " << neiProcNo
8608 << abort(FatalError);
8609 }
8610 else
8611 {
8612 // Positive rotation
8613 // - Set for old face. Will be rotated later
8614 // during the shuffling stage
8615 rotation[oldFaceI] =
8616 (
8617 (checkFace.size() - anchorFp) % checkFace.size()
8618 );
8619 }
8620 }
8622 // Set the flag
8623 anyChange = true;
8624 }
8625 }
8627 // Reduce failures across processors
8628 reduce(failedPatchMatch, orOp<bool>());
8630 // Write out processor patch faces if a failure was encountered
8631 if (failedPatchMatch)
8632 {
8633 for (label pI = 0; pI < nPatches_; pI++)
8634 {
8635 label neiProcNo = getNeighbourProcessor(pI);
8637 if (neiProcNo == -1)
8638 {
8639 continue;
8640 }
8642 writeVTK
8643 (
8644 "patchFaces_"
8645 + Foam::name(Pstream::myProcNo())
8646 + '_'
8647 + Foam::name(neiProcNo),
8648 identity(patchSizes_[pI]) + patchStarts_[pI],
8649 2
8650 );
8651 }
8653 FatalErrorIn
8654 (
8655 "\n"
8656 "void dynamicTopoFvMesh::"
8657 "syncCoupledBoundaryOrdering\n"
8658 "(\n"
8659 " List<pointField>& centres,\n"
8660 " List<pointField>& anchors,\n"
8661 " labelListList& patchMaps,\n"
8662 " labelListList& rotations\n"
8663 ") const\n"
8664 )
8665 << " Matching for processor patches failed. " << nl
8666 << " Patch faces written out to disk." << nl
8667 << abort(FatalError);
8668 }
8670 return anyChange;
8671 }
8674 // Fill buffers with length-scale info
8675 // and exchange across processors.
8676 void dynamicTopoFvMesh::exchangeLengthBuffers()
8677 {
8678 if (!Pstream::parRun())
8679 {
8680 return;
8681 }
8683 if (!edgeRefinement_)
8684 {
8685 return;
8686 }
8688 forAll(procIndices_, pI)
8689 {
8690 coupledInfo& sPM = sendMeshes_[pI];
8691 coupledInfo& rPM = recvMeshes_[pI];
8693 const coupleMap& scMap = sPM.map();
8694 const coupleMap& rcMap = rPM.map();
8696 if (scMap.slaveIndex() == Pstream::myProcNo())
8697 {
8698 // Clear existing buffer
8699 sPM.subMesh().lengthScale_.clear();
8701 // Fill in buffers to send.
8702 sPM.subMesh().lengthScale_.setSize
8703 (
8704 scMap.nEntities(coupleMap::CELL),
8705 -1.0
8706 );
8708 Map<label>& cellMap = scMap.entityMap(coupleMap::CELL);
8710 forAllIter(Map<label>, cellMap, cIter)
8711 {
8712 sPM.subMesh().lengthScale_[cIter.key()] = lengthScale_[cIter()];
8713 }
8715 meshOps::pWrite
8716 (
8717 scMap.masterIndex(),
8718 sPM.subMesh().lengthScale_
8719 );
8721 if (debug > 4)
8722 {
8723 Pout<< "Sending to: "
8724 << scMap.masterIndex()
8725 << " nCells: "
8726 << scMap.nEntities(coupleMap::CELL)
8727 << endl;
8728 }
8729 }
8731 if (rcMap.masterIndex() == Pstream::myProcNo())
8732 {
8733 // Clear existing buffer
8734 rPM.subMesh().lengthScale_.clear();
8736 // Schedule receipt from neighbour
8737 rPM.subMesh().lengthScale_.setSize
8738 (
8739 rcMap.nEntities(coupleMap::CELL),
8740 -1.0
8741 );
8743 // Schedule for receipt
8744 meshOps::pRead
8745 (
8746 rcMap.slaveIndex(),
8747 rPM.subMesh().lengthScale_
8748 );
8750 if (debug > 4)
8751 {
8752 Pout<< "Receiving from: "
8753 << rcMap.slaveIndex()
8754 << " nCells: "
8755 << rcMap.nEntities(coupleMap::CELL)
8756 << endl;
8757 }
8758 }
8759 }
8761 // Wait for transfers before continuing.
8762 meshOps::waitForBuffers();
8764 if (debug > 4)
8765 {
8766 forAll(procIndices_, pI)
8767 {
8768 const coupledInfo& rPM = recvMeshes_[pI];
8769 const coupleMap& rcMap = rPM.map();
8771 if (rcMap.masterIndex() == Pstream::myProcNo())
8772 {
8773 rPM.subMesh().writeVTK
8774 (
8775 "lengthScale_" + Foam::name(rcMap.slaveIndex()),
8776 identity(rPM.subMesh().nCells()),
8777 3,
8778 false,
8779 false,
8780 rPM.subMesh().lengthScale_
8781 );
8782 }
8783 }
8784 }
8785 }
8788 // Implementing the fillTables operation for coupled edges
8789 bool dynamicTopoFvMesh::coupledFillTables
8790 (
8791 const label eIndex,
8792 scalar& minQuality,
8793 labelList& m,
8794 PtrList<scalarListList>& Q,
8795 PtrList<labelListList>& K,
8796 PtrList<labelListList>& triangulations
8797 ) const
8798 {
8799 bool success = false;
8801 if (locallyCoupledEntity(eIndex))
8802 {
8803 // Fill tables for the slave edge.
8804 label sI = -1;
8806 // Determine the slave index.
8807 forAll(patchCoupling_, patchI)
8808 {
8809 if (patchCoupling_(patchI))
8810 {
8811 const label edgeEnum = coupleMap::EDGE;
8812 const coupleMap& cMap = patchCoupling_[patchI].map();
8814 if ((sI = cMap.findSlave(edgeEnum, eIndex)) > -1)
8815 {
8816 break;
8817 }
8818 }
8819 }
8821 if (sI == -1)
8822 {
8823 FatalErrorIn
8824 (
8825 "\n"
8826 "bool dynamicTopoFvMesh::coupledFillTables\n"
8827 "(\n"
8828 " const label eIndex,\n"
8829 " const scalar minQuality,\n"
8830 " labelList& m,\n"
8831 " PtrList<scalarListList>& Q,\n"
8832 " PtrList<labelListList>& K,\n"
8833 " PtrList<labelListList>& triangulations\n"
8834 ") const\n"
8835 )
8836 << "Coupled maps were improperly specified." << nl
8837 << " Slave index not found for: " << nl
8838 << " Edge: " << eIndex << nl
8839 << abort(FatalError);
8840 }
8842 // Turn off switch temporarily.
8843 unsetCoupledModification();
8845 labelList hullV(10);
8847 // Call fillTables for the slave edge.
8848 success = fillTables(sI, minQuality, m, hullV, Q, K, triangulations, 1);
8850 // Turn it back on.
8851 setCoupledModification();
8852 }
8853 else
8854 if (processorCoupledEntity(eIndex))
8855 {
8856 const edge& checkEdge = edges_[eIndex];
8857 const labelList& eFaces = edgeFaces_[eIndex];
8859 // Reset minQuality
8860 minQuality = GREAT;
8862 // Need to build alternate addressing / point-list
8863 // for swaps across processors.
8864 DynamicList<point> parPts(10);
8865 DynamicList<label> parVtx(10);
8867 bool closed = true;
8868 label nPoints = 0, nProcs = 0;
8869 label otherPoint = -1, nextPoint = -1;
8871 // Define a line for this edge
8872 linePointRef lpr
8873 (
8874 points_[checkEdge.start()],
8875 points_[checkEdge.end()]
8876 );
8878 // Define tangent-to-edge / edge-centre
8879 vector te = -lpr.vec(), xe = lpr.centre();
8881 // Normalize tangent
8882 te /= mag(te) + VSMALL;
8884 // Fill-in vertices for this processor
8885 forAll(eFaces, faceI)
8886 {
8887 label fPatch = whichPatch(eFaces[faceI]);
8888 label neiProc = getNeighbourProcessor(fPatch);
8890 if (neiProc > -1 && neiProc < Pstream::myProcNo())
8891 {
8892 // This edge should not be here
8893 Pout<< " Edge: " << eIndex
8894 << " is talking to processor: " << neiProc
8895 << abort(FatalError);
8896 }
8898 // This face is either internal or on a physical boundary
8899 const face& checkFace = faces_[eFaces[faceI]];
8901 // Find the isolated point
8902 meshOps::findIsolatedPoint
8903 (
8904 checkFace,
8905 checkEdge,
8906 otherPoint,
8907 nextPoint
8908 );
8910 // Insert point and index
8911 parPts.append(points_[otherPoint]);
8912 parVtx.append(nPoints);
8914 // Physical patch: Is this an appropriate start face?
8915 // - If yes, swap with first index
8916 if (fPatch > -1 && neiProc == -1)
8917 {
8918 closed = false;
8920 if (nextPoint == checkEdge[0])
8921 {
8922 Foam::Swap(parPts[0], parPts[nPoints]);
8923 }
8924 }
8926 // Increment point index
8927 nPoints++;
8928 }
8930 // Now look through processors, and add their points
8931 forAll(procIndices_, pI)
8932 {
8933 label proc = procIndices_[pI];
8935 // Fetch reference to subMesh
8936 const coupleMap& cMap = recvMeshes_[pI].map();
8937 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
8939 label sI = -1, sP = -1;
8941 if ((sI = cMap.findSlave(coupleMap::EDGE, eIndex)) == -1)
8942 {
8943 continue;
8944 }
8946 // If this is a new edge, bail out for now.
8947 // This will be handled at the next time-step.
8948 if (sI >= mesh.nOldEdges_)
8949 {
8950 return false;
8951 }
8953 const edge& slaveEdge = mesh.edges_[sI];
8954 const labelList& seFaces = mesh.edgeFaces_[sI];
8956 // Determine the point index that corresponds to checkEdge[0]
8957 edge cE
8958 (
8959 cMap.findMaster(coupleMap::POINT, slaveEdge[0]),
8960 cMap.findMaster(coupleMap::POINT, slaveEdge[1])
8961 );
8963 if (cE[0] == checkEdge[0])
8964 {
8965 sP = slaveEdge[0];
8966 }
8967 else
8968 if (cE[1] == checkEdge[0])
8969 {
8970 sP = slaveEdge[1];
8971 }
8972 else
8973 {
8974 label meP = whichEdgePatch(eIndex);
8975 word mN(meP < 0 ? "Internal" : boundaryMesh()[meP].name());
8977 label seP = mesh.whichEdgePatch(sI);
8978 word sN(seP < 0 ? "Internal" : mesh.boundaryMesh()[seP].name());
8980 Pout<< " Can't find master point: " << checkEdge[0] << nl
8981 << " on master: " << eIndex << "::" << checkEdge
8982 << " Patch: " << mN << nl
8983 << " with slave: " << sI << "::" << slaveEdge
8984 << " Patch: " << sN << nl
8985 << " cE: " << cE << " on proc: " << proc << nl
8986 << abort(FatalError);
8987 }
8989 forAll(seFaces, faceI)
8990 {
8991 label slavePatch = mesh.whichPatch(seFaces[faceI]);
8993 // If this is talking to a lower-ranked processor,
8994 // skip the insertion step.
8995 label neiProc = mesh.getNeighbourProcessor(slavePatch);
8997 if (neiProc > -1 && neiProc < proc)
8998 {
8999 continue;
9000 }
9002 // This face is either internal or on a physical boundary
9003 const face& slaveFace = mesh.faces_[seFaces[faceI]];
9005 // Find the isolated point
9006 meshOps::findIsolatedPoint
9007 (
9008 slaveFace,
9009 slaveEdge,
9010 otherPoint,
9011 nextPoint
9012 );
9014 // Insert point and index
9015 parPts.append(mesh.points_[otherPoint]);
9016 parVtx.append(nPoints);
9018 // Physical patch: Is this an appropriate start face?
9019 // - If yes, swap with first index
9020 if (slavePatch > -1 && neiProc == -1)
9021 {
9022 closed = false;
9024 if (nextPoint == sP)
9025 {
9026 Foam::Swap(parPts[0], parPts[nPoints]);
9027 }
9028 }
9030 // Increment point index
9031 nPoints++;
9032 }
9034 nProcs++;
9035 }
9037 // Sort points / indices in counter-clockwise order
9038 SortableList<scalar> angles(nPoints, 0.0);
9040 // Fetch 2 * pi
9041 scalar twoPi = mathematicalConstant::twoPi;
9043 // Define a base direction
9044 // from the start point
9045 vector dir = (parPts[0] - xe);
9046 dir -= ((dir & te) * te);
9047 dir /= mag(dir) + VSMALL;
9049 // Local coordinate system
9050 coordinateSystem cs("cs", xe, te, dir);
9052 for (label i = 1; i < nPoints; i++)
9053 {
9054 // Convert to local csys and determine angle
9055 vector local = cs.localPosition(parPts[i]);
9056 scalar angle = atan2(local.y(), local.x());
9058 // Account for 3rd and 4th quadrants
9059 angles[i] = (angle < 0.0 ? angle + twoPi : angle);
9060 }
9062 // Sort by angle
9063 angles.sort();
9065 // Reorder points and transfer
9066 List<point> sortedParPts(nPoints);
9068 const labelList& indices = angles.indices();
9070 forAll(sortedParPts, pointI)
9071 {
9072 sortedParPts[pointI] = parPts[indices[pointI]];
9073 }
9075 parPts.transfer(sortedParPts);
9077 // Fill the last two points for the edge
9078 edge parEdge(-1, -1);
9080 parPts.append(lpr.start());
9081 parEdge[0] = nPoints++;
9083 parPts.append(lpr.end());
9084 parEdge[1] = nPoints++;
9086 // Compute minQuality with this loop
9087 minQuality = computeMinQuality(parEdge, parVtx, parPts, closed);
9089 if (debug > 4 || minQuality < 0.0)
9090 {
9091 // Write out edge connectivity
9092 writeEdgeConnectivity(eIndex);
9094 meshOps::writeVTK
9095 (
9096 (*this),
9097 "parPts_" + Foam::name(eIndex),
9098 parPts.size(),
9099 parPts.size(),
9100 parPts.size(),
9101 pointField(parPts)
9102 );
9104 if (minQuality < 0.0)
9105 {
9106 Pout<< " * * * Error in fillTables * * * " << nl
9107 << " Edge: " << eIndex << " :: " << checkEdge << nl
9108 << " minQuality: " << minQuality << nl
9109 << " Closed: " << Switch::asText(closed) << nl
9110 << abort(FatalError);
9111 }
9112 }
9114 // Fill in the size
9115 m[0] = parVtx.size();
9117 // Check if a table-resize is necessary
9118 if (m[0] > maxTetsPerEdge_)
9119 {
9120 if (allowTableResize_)
9121 {
9122 // Resize the tables to account for
9123 // more tets per edge
9124 label& mtpe = const_cast<label&>(maxTetsPerEdge_);
9126 mtpe = m[0];
9128 // Clear tables for this index.
9129 Q[0].clear();
9130 K[0].clear();
9131 triangulations[0].clear();
9133 // Resize for this index.
9134 initTables(m, Q, K, triangulations);
9135 }
9136 else
9137 {
9138 // Can't resize. Bail out.
9139 return false;
9140 }
9141 }
9143 // Fill dynamic programming tables
9144 fillTables
9145 (
9146 parEdge,
9147 minQuality,
9148 m[0],
9149 parVtx,
9150 parPts,
9151 Q[0],
9152 K[0],
9153 triangulations[0]
9154 );
9156 success = true;
9157 }
9159 return success;
9160 }
9163 // Additional mapping contributions for coupled entities
9164 void dynamicTopoFvMesh::computeCoupledWeights
9165 (
9166 const label index,
9167 const label dimension,
9168 labelList& parents,
9169 scalarField& weights,
9170 vectorField& centres,
9171 bool output
9172 )
9173 {
9174 if (!Pstream::parRun())
9175 {
9176 return;
9177 }
9179 // Fetch offsets from mapper
9180 const labelList& cStarts = mapper_->cellStarts();
9181 const labelListList& pSizes = mapper_->patchSizes();
9183 if (dimension == 2)
9184 {
9185 DynamicList<label> faceParents(10);
9187 label patchIndex = whichPatch(index);
9189 forAll(procIndices_, pI)
9190 {
9191 // Ensure that the patch is physical
9192 if (patchIndex < 0 || patchIndex >= pSizes[pI].size())
9193 {
9194 Pout<< " Face: " << index
9195 << " Patch: " << patchIndex
9196 << " does not belong to a physical patch." << nl
9197 << " nPhysicalPatches: " << pSizes[pI].size()
9198 << abort(FatalError);
9199 }
9201 // Fetch reference to subMesh
9202 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
9204 // Prepare lists
9205 labelList coupleObjects;
9206 scalarField coupleWeights;
9207 vectorField coupleCentres;
9209 // Convex-set algorithm for faces
9210 faceSetAlgorithm faceAlgorithm
9211 (
9212 mesh,
9213 oldPoints_,
9214 edges_,
9215 faces_,
9216 cells_,
9217 owner_,
9218 neighbour_
9219 );
9221 // Initialize the bounding box
9222 faceAlgorithm.computeNormFactor(index);
9224 // Loop through all subMesh faces, and check for bounds
9225 const polyBoundaryMesh& boundary = mesh.boundaryMesh();
9227 label pSize = boundary[patchIndex].size();
9228 label pStart = boundary[patchIndex].start();
9230 for (label faceI = pStart; faceI < (pStart + pSize); faceI++)
9231 {
9232 if (faceAlgorithm.contains(faceI))
9233 {
9234 faceParents.append(faceI);
9235 }
9236 }
9238 // Obtain weighting factors for this face.
9239 faceAlgorithm.computeWeights
9240 (
9241 index,
9242 pStart,
9243 faceParents,
9244 boundary[patchIndex].faceFaces(),
9245 coupleObjects,
9246 coupleWeights,
9247 coupleCentres,
9248 output
9249 );
9251 // Add contributions with offsets
9252 if (coupleObjects.size())
9253 {
9254 // Fetch patch size on master
9255 label patchSize = boundaryMesh()[patchIndex].size();
9257 // Resize lists
9258 label oldSize = parents.size();
9260 parents.setSize(oldSize + coupleObjects.size());
9261 weights.setSize(oldSize + coupleObjects.size());
9262 centres.setSize(oldSize + coupleObjects.size());
9264 forAll(coupleObjects, indexI)
9265 {
9266 parents[indexI + oldSize] =
9267 (
9268 patchSize + coupleObjects[indexI]
9269 );
9271 weights[indexI + oldSize] = coupleWeights[indexI];
9272 centres[indexI + oldSize] = coupleCentres[indexI];
9273 }
9274 }
9276 // Clear list
9277 faceParents.clear();
9278 }
9279 }
9280 else
9281 if (dimension == 3)
9282 {
9283 DynamicList<label> cellParents(10);
9285 forAll(procIndices_, pI)
9286 {
9287 // Fetch reference to subMesh
9288 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
9290 // Prepare lists
9291 labelList coupleObjects;
9292 scalarField coupleWeights;
9293 vectorField coupleCentres;
9295 // Convex-set algorithm for cells
9296 cellSetAlgorithm cellAlgorithm
9297 (
9298 mesh,
9299 oldPoints_,
9300 edges_,
9301 faces_,
9302 cells_,
9303 owner_,
9304 neighbour_
9305 );
9307 // Initialize the bounding box
9308 cellAlgorithm.computeNormFactor(index);
9310 // Loop through all subMesh cells, and check for bounds
9311 const cellList& meshCells = mesh.cells();
9313 forAll(meshCells, cellI)
9314 {
9315 if (cellAlgorithm.contains(cellI))
9316 {
9317 cellParents.append(cellI);
9318 }
9319 }
9321 // Obtain weighting factors for this cell.
9322 cellAlgorithm.computeWeights
9323 (
9324 index,
9325 0,
9326 cellParents,
9327 mesh.polyMesh::cellCells(),
9328 coupleObjects,
9329 coupleWeights,
9330 coupleCentres,
9331 output
9332 );
9334 // Add contributions with offsets
9335 if (coupleObjects.size())
9336 {
9337 label cellStart = cStarts[pI];
9339 // Resize lists
9340 label oldSize = parents.size();
9342 parents.setSize(oldSize + coupleObjects.size());
9343 weights.setSize(oldSize + coupleObjects.size());
9344 centres.setSize(oldSize + coupleObjects.size());
9346 forAll(coupleObjects, indexI)
9347 {
9348 parents[indexI + oldSize] =
9349 (
9350 cellStart + coupleObjects[indexI]
9351 );
9353 weights[indexI + oldSize] = coupleWeights[indexI];
9354 centres[indexI + oldSize] = coupleCentres[indexI];
9355 }
9356 }
9358 // Clear list
9359 cellParents.clear();
9360 }
9361 }
9362 else
9363 {
9364 FatalErrorIn("scalar dynamicTopoFvMesh::computeCoupledWeights()")
9365 << " Incorrect dimension: " << dimension << nl
9366 << abort(FatalError);
9367 }
9368 }
9371 // Fetch length-scale info for processor entities
9372 scalar dynamicTopoFvMesh::processorLengthScale(const label index) const
9373 {
9374 scalar procScale = 0.0;
9376 if (twoDMesh_)
9377 {
9378 // First check the master processor
9379 procScale += lengthScale_[owner_[index]];
9381 // Next, check the slave processor
9382 bool foundSlave = false;
9384 forAll(procIndices_, pI)
9385 {
9386 // Fetch non-const reference to subMeshes
9387 const label faceEnum = coupleMap::FACE;
9388 const coupleMap& cMap = recvMeshes_[pI].map();
9389 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
9391 label sIndex = -1;
9393 if ((sIndex = cMap.findSlave(faceEnum, index)) > -1)
9394 {
9395 procScale += mesh.lengthScale_[mesh.owner_[sIndex]];
9397 foundSlave = true;
9398 break;
9399 }
9400 }
9402 // Should have found at least one slave
9403 if (!foundSlave)
9404 {
9405 FatalErrorIn
9406 (
9407 "scalar dynamicTopoFvMesh::processorLengthScale"
9408 "(const label index) const"
9409 )
9410 << "Processor lengthScale lookup failed: " << nl
9411 << " Master face: " << index
9412 << " :: " << faces_[index] << nl
9413 << abort(FatalError);
9414 }
9416 // Average the scale
9417 procScale *= 0.5;
9418 }
9419 else
9420 {
9421 // Check if this is a 'pure' processor edge
9422 bool pure = processorCoupledEntity(index, false, true, true);
9424 const label edgeEnum = coupleMap::EDGE;
9425 const labelList& eFaces = edgeFaces_[index];
9427 bool foundSlave = false;
9429 if (pure)
9430 {
9431 // First check the master processor
9432 label nC = 0;
9434 forAll(eFaces, faceI)
9435 {
9436 label own = owner_[eFaces[faceI]];
9437 label nei = neighbour_[eFaces[faceI]];
9439 procScale += lengthScale_[own];
9440 nC++;
9442 if (nei > -1)
9443 {
9444 procScale += lengthScale_[nei];
9445 nC++;
9446 }
9447 }
9449 // Next check slaves
9450 forAll(procIndices_, pI)
9451 {
9452 const coupleMap& cMap = recvMeshes_[pI].map();
9453 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
9455 label sIndex = -1;
9457 if ((sIndex = cMap.findSlave(edgeEnum, index)) > -1)
9458 {
9459 // Fetch connectivity from patchSubMesh
9460 const labelList& peFaces = mesh.edgeFaces_[sIndex];
9462 foundSlave = true;
9464 forAll(peFaces, faceI)
9465 {
9466 label own = mesh.owner_[peFaces[faceI]];
9467 label nei = mesh.neighbour_[peFaces[faceI]];
9469 procScale += mesh.lengthScale_[own];
9470 nC++;
9472 if (nei > -1)
9473 {
9474 procScale += mesh.lengthScale_[nei];
9475 nC++;
9476 }
9477 }
9478 }
9479 }
9481 // Average the final scale
9482 procScale /= nC;
9483 }
9484 else
9485 {
9486 // Processor is adjacent to physical patch types.
9487 // Search for boundary faces, and average their scale
9489 // First check the master processor
9490 label nBoundary = 0;
9492 forAll(eFaces, faceI)
9493 {
9494 label facePatch = whichPatch(eFaces[faceI]);
9496 // Skip internal faces
9497 if (facePatch == -1)
9498 {
9499 continue;
9500 }
9502 // Skip processor patches
9503 if (getNeighbourProcessor(facePatch) > -1)
9504 {
9505 continue;
9506 }
9508 // If this is a floating face, pick the owner length-scale
9509 if (lengthEstimator().isFreePatch(facePatch))
9510 {
9511 procScale += lengthScale_[owner_[eFaces[faceI]]];
9512 }
9513 else
9514 {
9515 // Fetch fixed length-scale
9516 procScale +=
9517 (
9518 lengthEstimator().fixedLengthScale
9519 (
9520 eFaces[faceI],
9521 facePatch
9522 )
9523 );
9524 }
9526 nBoundary++;
9527 }
9529 forAll(procIndices_, pI)
9530 {
9531 const coupleMap& cMap = recvMeshes_[pI].map();
9532 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
9534 label sIndex = -1;
9536 if ((sIndex = cMap.findSlave(edgeEnum, index)) > -1)
9537 {
9538 // Fetch connectivity from patchSubMesh
9539 const labelList& peFaces = mesh.edgeFaces_[sIndex];
9541 foundSlave = true;
9543 forAll(peFaces, faceI)
9544 {
9545 label facePatch = mesh.whichPatch(peFaces[faceI]);
9547 // Skip internal faces
9548 if (facePatch == -1)
9549 {
9550 continue;
9551 }
9553 // Skip processor patches
9554 if (mesh.getNeighbourProcessor(facePatch) > -1)
9555 {
9556 continue;
9557 }
9559 // If this is a floating face,
9560 // pick the owner length-scale
9561 if (lengthEstimator().isFreePatch(facePatch))
9562 {
9563 procScale +=
9564 (
9565 mesh.lengthScale_[mesh.owner_[peFaces[faceI]]]
9566 );
9567 }
9568 else
9569 {
9570 // Fetch fixed length-scale
9571 procScale +=
9572 (
9573 lengthEstimator().fixedLengthScale
9574 (
9575 peFaces[faceI],
9576 facePatch
9577 )
9578 );
9579 }
9581 nBoundary++;
9582 }
9583 }
9584 }
9586 if (nBoundary != 2)
9587 {
9588 // Write out for post-processing
9589 writeEdgeConnectivity(index);
9591 FatalErrorIn
9592 (
9593 "scalar dynamicTopoFvMesh::processorLengthScale"
9594 "(const label index) const"
9595 )
9596 << " Expected two physical boundary patches: " << nl
9597 << " nBoundary: " << nBoundary
9598 << " Master edge: " << index
9599 << " :: " << edges_[index]
9600 << " Patch: "
9601 << boundaryMesh()[whichEdgePatch(index)].name() << nl
9602 << abort(FatalError);
9603 }
9605 procScale *= 0.5;
9606 }
9608 // Should have found at least one slave
9609 if (!foundSlave)
9610 {
9611 FatalErrorIn
9612 (
9613 "scalar dynamicTopoFvMesh::processorLengthScale"
9614 "(const label index) const"
9615 )
9616 << "Processor lengthScale lookup failed: " << nl
9617 << " Master edge: " << index
9618 << " :: " << edges_[index] << nl
9619 << abort(FatalError);
9620 }
9621 }
9623 return procScale;
9624 }
9627 // Method to determine whether the master face is locally coupled
9628 bool dynamicTopoFvMesh::locallyCoupledEntity
9629 (
9630 const label index,
9631 bool checkSlaves,
9632 bool checkProcs,
9633 bool checkFace
9634 ) const
9635 {
9636 // Bail out if no patchCoupling is present
9637 if (patchCoupling_.empty())
9638 {
9639 return false;
9640 }
9642 if (twoDMesh_ || checkFace)
9643 {
9644 label patch = whichPatch(index);
9646 if (patch == -1)
9647 {
9648 return false;
9649 }
9651 // Processor checks receive priority.
9652 if (checkProcs)
9653 {
9654 if (getNeighbourProcessor(patch) > -1)
9655 {
9656 return false;
9657 }
9658 }
9660 // Check coupled master patches.
9661 if (patchCoupling_(patch))
9662 {
9663 return true;
9664 }
9665 else
9666 if (checkSlaves)
9667 {
9668 // Check on slave patches as well.
9669 forAll(patchCoupling_, pI)
9670 {
9671 if (patchCoupling_(pI))
9672 {
9673 const coupleMap& cMap = patchCoupling_[pI].map();
9675 if (cMap.slaveIndex() == patch)
9676 {
9677 return true;
9678 }
9679 }
9680 }
9681 }
9682 }
9683 else
9684 {
9685 const labelList& eFaces = edgeFaces_[index];
9687 // Search for boundary faces, and determine boundary type.
9688 forAll(eFaces, faceI)
9689 {
9690 if (neighbour_[eFaces[faceI]] == -1)
9691 {
9692 label patch = whichPatch(eFaces[faceI]);
9694 // Processor checks receive priority.
9695 if (checkProcs)
9696 {
9697 if (getNeighbourProcessor(patch) > -1)
9698 {
9699 return false;
9700 }
9701 }
9703 // Check coupled master patches.
9704 if (patchCoupling_(patch))
9705 {
9706 return true;
9707 }
9709 if (checkSlaves)
9710 {
9711 // Check on slave patches as well.
9712 forAll(patchCoupling_, pI)
9713 {
9714 if (patchCoupling_(pI))
9715 {
9716 const coupleMap& cMap =
9717 (
9718 patchCoupling_[pI].map()
9719 );
9721 if (cMap.slaveIndex() == patch)
9722 {
9723 return true;
9724 }
9725 }
9726 }
9727 }
9728 }
9729 }
9730 }
9732 // Could not find any faces on locally coupled patches.
9733 return false;
9734 }
9737 // Method to determine the locally coupled patch index
9738 label dynamicTopoFvMesh::locallyCoupledEdgePatch(const label eIndex) const
9739 {
9740 const labelList& eFaces = edgeFaces_[eIndex];
9742 // Search for boundary faces, and determine boundary type.
9743 forAll(eFaces, faceI)
9744 {
9745 if (neighbour_[eFaces[faceI]] == -1)
9746 {
9747 label patch = whichPatch(eFaces[faceI]);
9749 // Check coupled master patches.
9750 if (patchCoupling_(patch))
9751 {
9752 return patch;
9753 }
9755 // Check on slave patches as well.
9756 forAll(patchCoupling_, pI)
9757 {
9758 if (patchCoupling_(pI))
9759 {
9760 const coupleMap& cMap = patchCoupling_[pI].map();
9762 if (cMap.slaveIndex() == patch)
9763 {
9764 return patch;
9765 }
9766 }
9767 }
9768 }
9769 }
9771 // Could not find any faces on locally coupled patches.
9772 FatalErrorIn
9773 (
9774 "label dynamicTopoFvMesh::locallyCoupledEdgePatch"
9775 "(const label cIndex) const"
9776 )
9777 << "Edge: " << eIndex << ":: " << edges_[eIndex]
9778 << " does not lie on any coupled patches."
9779 << abort(FatalError);
9781 return -1;
9782 }
9785 // Method to determine if the entity is on a processor boundary
9786 // - Also provides an additional check for 'pure' processor edges
9787 // i.e., edges that do not abut a physical patch. This is necessary
9788 // while deciding on collapse cases towards bounding curves.
9789 bool dynamicTopoFvMesh::processorCoupledEntity
9790 (
9791 const label index,
9792 bool checkFace,
9793 bool checkEdge,
9794 bool checkPure,
9795 FixedList<label, 2>* patchLabels,
9796 FixedList<vector, 2>* patchNormals
9797 ) const
9798 {
9799 // Skip check for serial runs
9800 if (!Pstream::parRun())
9801 {
9802 return false;
9803 }
9805 label patch = -2;
9807 if ((twoDMesh_ || checkFace) && !checkEdge)
9808 {
9809 patch = whichPatch(index);
9811 if (patch == -1)
9812 {
9813 return false;
9814 }
9816 if (getNeighbourProcessor(patch) > -1)
9817 {
9818 return true;
9819 }
9820 }
9821 else
9822 {
9823 const labelList& eFaces = edgeFaces_[index];
9825 label nPhysical = 0, nProcessor = 0;
9827 // Search for boundary faces, and determine boundary type.
9828 forAll(eFaces, faceI)
9829 {
9830 label patch = whichPatch(eFaces[faceI]);
9832 if (patch == -1)
9833 {
9834 continue;
9835 }
9837 if (getNeighbourProcessor(patch) > -1)
9838 {
9839 // Increment the processor patch count
9840 nProcessor++;
9842 if (!checkPure)
9843 {
9844 // We don't have to validate if this
9845 // is a 'pure' processor edge, so bail out.
9846 return true;
9847 }
9848 }
9849 else
9850 {
9851 // Physical patch.
9852 if (patchLabels)
9853 {
9854 (*patchLabels)[nPhysical] = patch;
9855 }
9857 if (patchNormals)
9858 {
9859 (*patchNormals)[nPhysical] =
9860 (
9861 faces_[eFaces[faceI]].normal(points_)
9862 );
9863 }
9865 nPhysical++;
9866 }
9867 }
9869 if (patchLabels && patchNormals)
9870 {
9871 // Check other coupled-edges as well
9872 forAll(procIndices_, pI)
9873 {
9874 // Fetch reference to subMesh
9875 const coupleMap& cMap = recvMeshes_[pI].map();
9876 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
9878 label sI = -1;
9880 if ((sI = cMap.findSlave(coupleMap::EDGE, index)) == -1)
9881 {
9882 continue;
9883 }
9885 const labelList& seFaces = mesh.edgeFaces_[sI];
9887 forAll(seFaces, faceI)
9888 {
9889 label sPatch = mesh.whichPatch(seFaces[faceI]);
9890 label neiProc = mesh.getNeighbourProcessor(sPatch);
9892 // Skip internal / processor faces
9893 if (sPatch == -1 || neiProc > -1)
9894 {
9895 continue;
9896 }
9898 const face& sFace = mesh.faces_[seFaces[faceI]];
9900 // Fill patch / normal info
9901 (*patchLabels)[nPhysical] = sPatch;
9902 (*patchNormals)[nPhysical] = sFace.normal(mesh.points_);
9904 nPhysical++;
9905 }
9906 }
9907 }
9909 // Purity check
9910 if (checkPure)
9911 {
9912 if (nProcessor && !nPhysical)
9913 {
9914 return true;
9915 }
9916 }
9917 }
9919 // Could not find any faces on processor patches.
9920 return false;
9921 }
9924 // Build a list of entities that need to be avoided
9925 // by regular topo-changes.
9926 void dynamicTopoFvMesh::buildEntitiesToAvoid
9927 (
9928 labelHashSet& entities,
9929 bool checkSubMesh
9930 )
9931 {
9932 entities.clear();
9934 // Build a set of entities to avoid during regular modifications,
9935 // and build a master stack for coupled modifications.
9937 // Determine locally coupled slave patches.
9938 labelHashSet localMasterPatches, localSlavePatches;
9940 forAll(patchCoupling_, patchI)
9941 {
9942 if (patchCoupling_(patchI))
9943 {
9944 const coupleMap& cMap = patchCoupling_[patchI].map();
9946 localMasterPatches.insert(cMap.masterIndex());
9947 localSlavePatches.insert(cMap.slaveIndex());
9948 }
9949 }
9951 // Loop through boundary faces and check whether
9952 // they belong to master/slave coupled patches.
9953 for (label faceI = nOldInternalFaces_; faceI < faces_.size(); faceI++)
9954 {
9955 // Add only valid faces
9956 if (faces_[faceI].empty())
9957 {
9958 continue;
9959 }
9961 label pIndex = whichPatch(faceI);
9963 if (pIndex == -1)
9964 {
9965 continue;
9966 }
9968 // Check if this is a coupled face.
9969 if
9970 (
9971 localMasterPatches.found(pIndex) ||
9972 localSlavePatches.found(pIndex) ||
9973 getNeighbourProcessor(pIndex) > -1
9974 )
9975 {
9976 if (twoDMesh_)
9977 {
9978 // Avoid this face during regular modification.
9979 if (!entities.found(faceI))
9980 {
9981 entities.insert(faceI);
9982 }
9983 }
9984 else
9985 {
9986 const labelList& fEdges = faceEdges_[faceI];
9988 forAll(fEdges, edgeI)
9989 {
9990 // Avoid this edge during regular modification.
9991 if (!entities.found(fEdges[edgeI]))
9992 {
9993 entities.insert(fEdges[edgeI]);
9994 }
9995 }
9996 }
9997 }
9998 }
10000 // Loop through entities contained in patchSubMeshes, if requested
10001 if (checkSubMesh)
10002 {
10003 forAll(procIndices_, pI)
10004 {
10005 const coupleMap& cMap = sendMeshes_[pI].map();
10006 const Map<label> rEdgeMap = cMap.reverseEntityMap(coupleMap::EDGE);
10008 if (cMap.slaveIndex() == Pstream::myProcNo())
10009 {
10010 forAllConstIter(Map<label>, rEdgeMap, eIter)
10011 {
10012 if (twoDMesh_)
10013 {
10014 const labelList& eFaces = edgeFaces_[eIter.key()];
10016 forAll(eFaces, faceI)
10017 {
10018 if (faces_[eFaces[faceI]].size() == 4)
10019 {
10020 if (!entities.found(eFaces[faceI]))
10021 {
10022 entities.insert(eFaces[faceI]);
10023 }
10024 }
10025 }
10026 }
10027 else
10028 {
10029 const edge& check = edges_[eIter.key()];
10030 const labelList& eFaces = edgeFaces_[eIter.key()];
10032 // Skip deleted edges
10033 if (eFaces.size())
10034 {
10035 forAll(check, pI)
10036 {
10037 const labelList& pE = pointEdges_[check[pI]];
10039 forAll(pE, edgeI)
10040 {
10041 if (!entities.found(pE[edgeI]))
10042 {
10043 entities.insert(pE[edgeI]);
10044 }
10045 }
10046 }
10047 }
10048 }
10049 }
10050 }
10051 }
10052 }
10054 if (debug > 3)
10055 {
10056 Pout<< nl << "nEntitiesToAvoid: " << entities.size() << endl;
10058 if (debug > 4)
10059 {
10060 // Write out entities
10061 label elemType = twoDMesh_ ? 2 : 1;
10063 writeVTK
10064 (
10065 "entitiesToAvoid_"
10066 + Foam::name(Pstream::myProcNo()),
10067 entities.toc(),
10068 elemType
10069 );
10070 }
10071 }
10072 }
10075 // Check whether the specified edge is a coupled master edge.
10076 bool dynamicTopoFvMesh::isCoupledMaster
10077 (
10078 const label eIndex
10079 ) const
10080 {
10081 if (!coupledModification_)
10082 {
10083 return true;
10084 }
10086 return locallyCoupledEntity(eIndex);
10087 }
10090 } // End namespace Foam
10092 // ************************************************************************* //