1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | foam-extend: Open Source CFD
4 \\ / O peration | Version: 3.2
5 \\ / A nd | Web: http://www.foam-extend.org
6 \\/ M anipulation | For copyright notice see file Copyright
7 -------------------------------------------------------------------------------
9 This file is part of foam-extend.
11 foam-extend is free software: you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation, either version 3 of the License, or (at your
14 option) any later version.
16 foam-extend is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 General Public License for more details.
21 You should have received a copy of the GNU General Public License
22 along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
28 Functions specific to coupled connectivity
32 University of Massachusetts Amherst
35 \*---------------------------------------------------------------------------*/
37 #include "dynamicTopoFvMesh.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"
64 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
66 defineTypeNameAndDebug(coupleMap, 0);
69 // Geometric relative match tolerance
70 static const Foam::debug::tolerancesSwitch geomMatchTol_
76 // Priority scheme enumerants
86 // Priority scheme names
87 static const char* prioritySchemeNames_[MAX_PRIORITY_SCHEMES + 1] =
98 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
100 // Set coupled modification
101 void dynamicTopoFvMesh::setCoupledModification() const
103 coupledModification_ = true;
107 // Unset coupled modification
108 void dynamicTopoFvMesh::unsetCoupledModification() const
110 coupledModification_ = false;
114 // Initialize the coupled stack
115 void dynamicTopoFvMesh::initCoupledStack
117 const labelHashSet& entities,
121 // Clear existing lists/stacks.
128 // Initialize the stack with entries
129 // in the labelHashSet and return
130 forAllConstIter(labelHashSet, entities, eIter)
132 // Add only valid entities
136 faces_[eIter.key()].empty() :
137 edgeFaces_[eIter.key()].empty()
145 stack(0).insert(eIter.key());
148 if (debug > 3 && Pstream::parRun())
150 Pout<< nl << "Entity stack size: " << stack(0).size() << endl;
154 // Write out stack entities
155 labelList stackElements(stack(0).size(), -1);
157 forAll(stackElements, elemI)
159 stackElements[elemI] = stack(0)[elemI];
162 label elemType = is2D() ? 2 : 1;
167 + Foam::name(Pstream::myProcNo()),
177 // Loop though boundary faces and check whether
178 // they belong to master/slave coupled patches.
179 for (label faceI = nOldInternalFaces_; faceI < faces_.size(); faceI++)
181 // Add only valid faces
182 if (faces_[faceI].empty())
187 label pIndex = whichPatch(faceI);
194 // Check if this is a locally coupled master face.
195 if (patchCoupling_.size())
197 if (patchCoupling_(pIndex))
200 bool addIndex = true;
202 if (isA<cyclicPolyPatch>(boundaryMesh()[pIndex]))
204 // Check if this is a cyclic master face
205 const coupleMap& cMap = patchCoupling_[pIndex].map();
206 const Map<label>& fMap = cMap.entityMap(coupleMap::FACE);
208 if (!fMap.found(faceI))
214 // Add this to the coupled modification stack.
219 stack(0).push(faceI);
223 const labelList& mfEdges = faceEdges_[faceI];
225 forAll(mfEdges, edgeI)
227 // Add this to the coupled modification stack.
228 stack(0).push(mfEdges[edgeI]);
235 // Check if this is a processor patch.
236 label myProcID = Pstream::myProcNo();
237 label neiProcID = getNeighbourProcessor(pIndex);
239 // Check if this is a master processor patch.
243 && priority(neiProcID, greaterOp<label>(), myProcID)
246 // Add this to the coupled modification stack.
249 stack(0).push(faceI);
253 const label edgeEnum = coupleMap::EDGE;
254 const label pointEnum = coupleMap::POINT;
256 const labelList& mfEdges = faceEdges_[faceI];
258 forAll(mfEdges, edgeI)
260 label eIndex = mfEdges[edgeI];
262 const edge& checkEdge = edges_[eIndex];
264 // Need to avoid this edge if it is
265 // talking to higher-priority processors.
266 bool permissible = true;
268 forAll(procIndices_, pI)
270 label nProc = procIndices_[pI];
272 // List is specified in ascending order
273 // of neighbour processors, so break out.
274 if (priority(nProc, greaterOp<label>(), myProcID))
279 // Fetch reference to subMeshes
280 const coupledMesh& recvMesh = recvMeshes_[pI];
281 const coupleMap& rcMap = recvMesh.map();
283 if (rcMap.findSlave(edgeEnum, eIndex) > -1)
289 // Check points as well, since there might be
290 // situations where both points lie on a lower
291 // ranked processor, but not the edge itself.
294 (rcMap.findSlave(pointEnum, checkEdge[0]) > -1)
295 && (rcMap.findSlave(pointEnum, checkEdge[1]) > -1)
305 stack(0).push(eIndex);
312 if (debug > 3 && Pstream::parRun())
314 Pout<< nl << "Coupled stack size: " << stack(0).size() << endl;
318 // Write out stack entities
319 labelList stackElements(stack(0).size(), -1);
321 forAll(stackElements, elemI)
323 stackElements[elemI] = stack(0)[elemI];
326 label elemType = is2D() ? 2 : 1;
331 + Foam::name(Pstream::myProcNo()),
340 // Execute load balancing operations on the mesh
341 void dynamicTopoFvMesh::executeLoadBalancing
343 dictionary& balanceDict
346 if (!Pstream::parRun())
351 bool enabled = readBool(balanceDict.lookup("enabled"));
359 label interval = readLabel(balanceDict.lookup("interval"));
361 // Are we at the interval?
362 if ((interval < 0) || (time().timeIndex() % interval != 0))
369 Info<< " void dynamicTopoFvMesh::executeLoadBalancing() :"
370 << " Executing load-balancing operations on the mesh."
373 // Write out patch boundaries before re-distribution
374 for (label pI = 0; pI < nPatches_; pI++)
376 label neiProcNo = getNeighbourProcessor(pI);
386 + Foam::name(Pstream::myProcNo())
388 + Foam::name(neiProcNo),
389 identity(patchSizes_[pI]) + patchStarts_[pI],
396 + Foam::name(Pstream::myProcNo())
398 + Foam::name(neiProcNo),
399 identity(patchSizes_[pI]) + patchStarts_[pI],
405 // Ensure that the number of processors is consistent,
406 // and silently modify dictionary entry if not
407 balanceDict.set("numberOfSubdomains", Pstream::nProcs());
409 // Read decomposition options and initialize
410 autoPtr<decompositionMethod> decomposerPtr
412 decompositionMethod::New
419 // Alias for convenience...
420 decompositionMethod& decomposer = decomposerPtr();
422 if (!decomposer.parallelAware())
424 FatalErrorIn("void dynamicTopoFvMesh::executeLoadBalancing()")
425 << "You have selected decomposition method "
426 << decomposer.typeName
427 << " which is not parallel aware."
431 // Calculate a merge-distance
432 const boundBox& box = polyMesh::bounds();
433 scalar mergeTol = readScalar(balanceDict.lookup("mergeTol"));
436 scalar mergeDist = mergeTol * box.mag();
438 // Compute write tolerance
444 -scalar(IOstream::defaultPrecision())
449 << " Overall mesh bounding box : " << box << nl
450 << " Relative tolerance : " << mergeTol << nl
451 << " Absolute matching distance : " << mergeDist << nl
454 // Check for compatibility
455 if (time().writeFormat() == IOstream::ASCII && mergeTol < writeTol)
457 FatalErrorIn("void dynamicTopoFvMesh::executeLoadBalancing()")
458 << " Your current settings specify ASCII writing with "
459 << IOstream::defaultPrecision() << " digits precision." << nl
460 << " Your mergeTol (" << mergeTol << ") is finer than this." << nl
461 << " Remedies: " << nl
462 << " - Change writeFormat to binary" << nl
463 << " - Increase writePrecision" << nl
464 << " - Adjust the merge tolerance (mergeTol)" << endl
468 // Decompose the mesh and obtain distribution
469 labelList cellDistribution
473 primitiveMesh::cellCentres(),
474 scalarField(nCells_, 1)
478 // Set the coupledModification switch
479 setCoupledModification();
481 // Initialize the mesh distribution engine
482 fvMeshDistribute distributor(*this, mergeDist);
484 // Re-distribute mesh according to new decomposition
485 distributor.distribute(cellDistribution);
487 // Reset the coupledModification switch
488 unsetCoupledModification();
490 // Set old / new sizes
491 nPoints_ = nOldPoints_ = primitiveMesh::nPoints();
492 nEdges_ = nOldEdges_ = 0;
493 nFaces_ = nOldFaces_ = primitiveMesh::nFaces();
494 nCells_ = nOldCells_ = primitiveMesh::nCells();
495 nInternalFaces_ = nOldInternalFaces_ = primitiveMesh::nInternalFaces();
496 nInternalEdges_ = nOldInternalEdges_ = 0;
498 // Now re-initialize all connectivity structures
499 oldPoints_ = polyMesh::points();
500 points_ = polyMesh::points();
501 faces_ = polyMesh::faces();
502 owner_ = polyMesh::faceOwner();
503 neighbour_ = polyMesh::faceNeighbour();
504 cells_ = primitiveMesh::cells();
506 // Check the size of owner / neighbour
507 if (owner_.size() != neighbour_.size())
509 // Size up to number of faces
510 neighbour_.setSize(nFaces_, -1);
513 const polyBoundaryMesh& boundary = polyMesh::boundaryMesh();
515 // Initialize patch-size information
516 nPatches_ = boundary.size();
518 oldPatchSizes_.setSize(nPatches_, 0);
519 oldPatchStarts_.setSize(nPatches_, -1);
520 oldEdgePatchSizes_.setSize(nPatches_, 0);
521 oldEdgePatchStarts_.setSize(nPatches_, -1);
522 oldPatchNMeshPoints_.setSize(nPatches_, -1);
524 patchSizes_.setSize(nPatches_, 0);
525 patchStarts_.setSize(nPatches_, -1);
526 edgePatchSizes_.setSize(nPatches_, 0);
527 edgePatchStarts_.setSize(nPatches_, -1);
528 patchNMeshPoints_.setSize(nPatches_, -1);
530 for (label i = 0; i < nPatches_; i++)
532 patchNMeshPoints_[i] = boundary[i].meshPoints().size();
533 oldPatchSizes_[i] = patchSizes_[i] = boundary[i].size();
534 oldPatchStarts_[i] = patchStarts_[i] = boundary[i].start();
535 oldPatchNMeshPoints_[i] = patchNMeshPoints_[i];
538 // Clear pointers and re-initialize
541 motionSolver_.clear();
542 lengthEstimator_.clear();
544 // Initialize edge-related connectivity structures
547 // Load the mesh-motion solver
550 // Load the field-mapper
553 // Load the length-scale estimator,
554 // and read refinement options
555 loadLengthScaleEstimator();
557 // Clear parallel structures
558 procIndices_.clear();
559 procPriority_.clear();
565 // Set processor rank priority
566 void dynamicTopoFvMesh::initProcessorPriority()
568 if (!Pstream::parRun())
573 if (Pstream::master())
578 const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
580 if (meshSubDict.found("priorityScheme") || mandatory_)
582 word schemeType(meshSubDict.lookup("priorityScheme"));
584 type = MAX_PRIORITY_SCHEMES;
586 for (label i = 0; i < MAX_PRIORITY_SCHEMES; i++)
588 if (prioritySchemeNames_[i] == schemeType)
600 // Initialize to identity map
601 procPriority_ = identity(Pstream::nProcs());
607 Random randomizer(std::time(NULL));
609 // Initialize to identity map
610 procPriority_ = identity(Pstream::nProcs());
612 for (label i = 0; i < Pstream::nProcs(); i++)
614 // Specify a random index to shuffle
615 label j = randomizer.integer(0, Pstream::nProcs() - 1);
617 Foam::Swap(procPriority_[i], procPriority_[j]);
625 const label listSize = Pstream::nProcs();
627 // Initialize to identity map
628 procPriority_ = identity(listSize);
630 // In-place reverse macro
631 # define inPlaceReverse(list) \
633 const label listSize = list.size(); \
634 const label lastIndex = listSize - 1; \
635 const label nIterations = listSize >> 1; \
638 while (elemI < nIterations) \
640 Swap(list[elemI], list[lastIndex - elemI]); \
645 // Rotate by time index
646 label n = time().timeIndex() % listSize;
653 // Rotate list in place
654 SubList<label> firstHalf(procPriority_, n, 0);
655 SubList<label> secondHalf(procPriority_, listSize - n, n);
657 inPlaceReverse(firstHalf);
658 inPlaceReverse(secondHalf);
660 inPlaceReverse(procPriority_);
667 procPriority_ = labelList(meshSubDict.lookup("priorityList"));
669 if (procPriority_.size() != Pstream::nProcs())
673 "void dynamicTopoFvMesh::initProcessorPriority()"
675 << " Priority scheme size mismatch." << nl
676 << " List size: " << procPriority_.size() << nl
677 << " Processor count: " << Pstream::nProcs() << nl
678 << abort(FatalError);
686 Info<< "Available schemes: " << endl;
688 for (label i = 0; i < MAX_PRIORITY_SCHEMES; i++)
690 Info<< ' ' << prioritySchemeNames_[i] << endl;
695 "void dynamicTopoFvMesh::initProcessorPriority()"
697 << " Unknown priority scheme." << nl
698 << abort(FatalError);
702 // Send priority list to others
703 for (label proc = 1; proc < Pstream::nProcs(); proc++)
705 meshOps::pWrite(proc, procPriority_);
710 // Receive priority from master
711 procPriority_.setSize(Pstream::nProcs());
712 meshOps::pRead(Pstream::masterNo(), procPriority_);
715 // Wait for transfers to complete
716 meshOps::waitForBuffers();
720 labelList ranks(Pstream::nProcs(), -1);
724 ranks[procPriority_[procI]] = procI;
727 Info<< " Processor priority: " << procPriority_ << nl
728 << " Processor rankings: " << ranks
734 // Identify coupled patches.
735 // - Also builds global shared point information.
736 // - Returns true if no coupled patches were found.
737 bool dynamicTopoFvMesh::identifyCoupledPatches()
739 bool coupledPatchesAbsent = true;
741 // Check if patches are explicitly coupled
742 if (patchCoupling_.size())
744 coupledPatchesAbsent = false;
747 // Maintain a separate list of processor IDs in procIndices.
748 // This is done because this sub-domain may talk to processors
749 // that share only edges/points.
750 if (!Pstream::parRun() || procIndices_.size())
752 return coupledPatchesAbsent;
755 // Prepare a list of points for sub-mesh creation.
756 // - Obtain global shared-points information, if necessary.
757 const polyBoundaryMesh& boundary = boundaryMesh();
759 // Set flag as default for parallel runs
760 coupledPatchesAbsent = false;
762 // Fetch the list of global points from polyMesh.
763 // - This should be the first evaluation of globalData,
764 // so this involves global communication
765 const globalMeshData& gData = polyMesh::globalData();
767 const labelList& spA = gData.sharedPointAddr();
768 const labelList& spL = gData.sharedPointLabels();
770 labelListList spB(Pstream::nProcs(), labelList(0));
772 if (gData.nGlobalPoints())
776 Info<< " dynamicTopoFvMesh::identifyCoupledPatches :"
777 << " Found " << gData.nGlobalPoints()
778 << " global points." << endl;
781 // Send others my addressing.
782 for (label proc = 0; proc < Pstream::nProcs(); proc++)
784 if (proc != Pstream::myProcNo())
786 // Send number of entities first.
787 meshOps::pWrite(proc, spA.size());
792 meshOps::pWrite(proc, spA);
797 // Receive addressing from others
798 for (label proc = 0; proc < Pstream::nProcs(); proc++)
800 if (proc != Pstream::myProcNo())
802 label procInfoSize = -1;
804 // How many entities am I going to be receiving?
805 meshOps::pRead(proc, procInfoSize);
809 // Size the receive buffer.
810 spB[proc].setSize(procInfoSize, -1);
812 // Schedule for receipt.
813 meshOps::pRead(proc, spB[proc]);
821 Info<< " dynamicTopoFvMesh::identifyCoupledPatches :"
822 << " Did not find any global points." << endl;
825 Map<label> immNeighbours;
826 labelListList procPoints(Pstream::nProcs());
828 // Track globally shared points
829 List<DynamicList<labelPair> > globalProcPoints
832 DynamicList<labelPair>(5)
835 // Insert my immediate neighbours into the list.
838 if (isA<processorPolyPatch>(boundary[pI]))
840 const processorPolyPatch& pp =
842 refCast<const processorPolyPatch>(boundary[pI])
845 label neiProcNo = pp.neighbProcNo();
847 // Insert all boundary points.
848 procPoints[neiProcNo] = pp.meshPoints();
850 // Keep track of immediate neighbours.
851 immNeighbours.insert(neiProcNo, pI);
855 if (gData.nGlobalPoints())
857 // Wait for all transfers to complete.
858 meshOps::waitForBuffers();
860 // Now loop through all processor addressing, and check if
861 // any labels coincide with my global shared points.
862 // If this is true, we need to be talking to that neighbour
863 // as well (if not already).
864 for (label proc = 0; proc < Pstream::nProcs(); proc++)
866 if (proc == Pstream::myProcNo())
871 bool foundGlobalMatch = false;
873 // Fetch reference to buffer
874 const labelList& procBuffer = spB[proc];
876 forAll(procBuffer, pointI)
880 if (spA[pointJ] == procBuffer[pointI])
882 // Make an entry, if one wasn't made already
883 if (findIndex(procPoints[proc], spL[pointJ]) == -1)
885 globalProcPoints[proc].append
887 labelPair(spL[pointJ], spA[pointJ])
891 foundGlobalMatch = true;
898 if (!immNeighbours.found(proc))
900 if (foundGlobalMatch && debug)
902 Pout<< " dynamicTopoFvMesh::identifyCoupledPatches :"
903 << " Additionally talking to processor: "
910 // Estimate an initial size
911 procIndices_.setSize(Pstream::nProcs());
913 label nTotalProcs = 0;
915 forAll(procPoints, procI)
917 if (procPoints[procI].size())
919 procIndices_[nTotalProcs++] = procI;
922 // Check for point / edge coupling
923 if (globalProcPoints[procI].size() && !procPoints[procI].size())
925 procIndices_[nTotalProcs++] = procI;
929 // Shrink to actual size
930 procIndices_.setSize(nTotalProcs);
932 // Initialize the processor priority list
933 initProcessorPriority();
935 // Sort by order of neighbouring
936 // processor priority - highest to lowest
940 UList<label>::less(procPriority_)
943 // Size the PtrLists.
944 sendMeshes_.setSize(nTotalProcs);
945 recvMeshes_.setSize(nTotalProcs);
947 // Create send/recv patch meshes, and copy
948 // the list of points for each processor.
949 forAll(procIndices_, pI)
951 label proc = procIndices_[pI], master = -1, slave = -1;
953 // For processors that have only point / edge coupling
954 // specify an invalid patch ID for now
955 label patchID = immNeighbours.found(proc) ? immNeighbours[proc] : -1;
957 // Check if this processor is of higher priority
958 if (priority(proc, lessOp<label>(), Pstream::myProcNo()))
961 slave = Pstream::myProcNo();
965 master = Pstream::myProcNo();
974 *this, // Reference to this mesh
977 true, // Sent to neighbour
978 patchID, // Patch index
979 master, // Master index
984 sendMeshes_[pI].map().subMeshPoints() =
986 procPoints[procIndices_[pI]]
989 sendMeshes_[pI].map().globalProcPoints() =
991 globalProcPoints[procIndices_[pI]]
999 *this, // Reference to this mesh
1002 false, // Not sent to neighbour
1003 patchID, // Patch index
1004 master, // Master index
1005 slave // Slave index
1009 recvMeshes_[pI].map().subMeshPoints() =
1011 procPoints[procIndices_[pI]]
1014 recvMeshes_[pI].map().globalProcPoints() =
1016 globalProcPoints[procIndices_[pI]]
1022 Pout<< " identifyCoupledPatches :"
1023 << " Talking to processors: "
1024 << procIndices_ << endl;
1026 forAll(procIndices_, pI)
1028 label proc = procIndices_[pI];
1030 // Write out subMeshPoints as a VTK
1031 const coupleMap& cMap = sendMeshes_[pI].map();
1036 Foam::name(Pstream::myProcNo()) + "to" +
1038 cMap.subMeshPoints(),
1042 // Write out globalProcPoints as a VTK
1043 const List<labelPair>& gpp = cMap.globalProcPoints();
1047 labelList gpPoints(gpp.size(), -1);
1051 gpPoints[pointI] = gpp[pointI].first();
1056 "globalProcPoints_" +
1057 Foam::name(Pstream::myProcNo()) + "to" +
1066 return coupledPatchesAbsent;
1070 // Read coupled patch information from dictionary.
1071 void dynamicTopoFvMesh::readCoupledPatches()
1073 const polyBoundaryMesh& boundary = boundaryMesh();
1076 patchCoupling_.clear();
1078 if (dict_.found("coupledPatches") || mandatory_)
1080 // Size the initial list
1081 patchCoupling_.setSize(boundary.size());
1083 const dictionary& coupledPatches = dict_.subDict("coupledPatches");
1085 // Determine master and slave patches
1086 forAllConstIter(dictionary, coupledPatches, dIter)
1088 const dictionary& dictI = dIter().dict();
1090 // Lookup the master / slave patches
1091 word masterPatch = dictI.lookup("master");
1092 word slavePatch = dictI.lookup("slave");
1094 // Determine patch indices
1095 label mPatch = boundary.findPatchID(masterPatch);
1096 label sPatch = boundary.findPatchID(slavePatch);
1100 Info<< " Adding master: " << masterPatch << " : " << mPatch
1101 << " with slave: " << slavePatch << " : " << sPatch
1105 // Add to the list if entries are legitimate
1109 boundary[mPatch].size() == boundary[sPatch].size()
1112 // Check whether patches are associated with zones.
1115 dictI.lookup("specifyZones")
1118 label mZone = -1, sZone = -1;
1122 const faceZoneMesh& faceZones = polyMesh::faceZones();
1124 mZone = faceZones.findZoneID
1126 dictI.lookup("masterZone")
1129 sZone = faceZones.findZoneID
1131 dictI.lookup("slaveZone")
1135 // Configure with regIOobject for check-in
1141 + Foam::name(mPatch)
1143 + Foam::name(sPatch)
1147 IOobject::READ_IF_PRESENT,
1148 IOobject::AUTO_WRITE,
1174 FatalErrorIn("void dynamicTopoFvMesh::readCoupledPatches()")
1175 << " Coupled patches are either wrongly specified,"
1176 << " or the sizes don't match." << nl
1177 << " Master: " << mPatch << ":" << masterPatch
1178 << " Size: " << boundary[mPatch].size() << nl
1179 << " Slave: " << sPatch << ":" << slavePatch
1180 << " Size: " << boundary[sPatch].size() << nl
1181 << abort(FatalError);
1186 // Loop through boundaries and add any cyclic patches
1187 forAll(boundary, patchI)
1189 if (!isA<cyclicPolyPatch>(boundary[patchI]))
1194 // Size up patchCoupling, if necessary
1195 if (patchCoupling_.empty())
1197 patchCoupling_.setSize(boundary.size());
1202 Info<< " Adding cyclic: " << patchI
1203 << " : " << boundary[patchI].name()
1207 // Configure with regIOobject for check-in
1213 + Foam::name(patchI)
1215 + Foam::name(patchI)
1219 IOobject::READ_IF_PRESENT,
1220 IOobject::AUTO_WRITE,
1235 new coupledMesh(*this, cMap, -1, -1)
1241 // Initialize coupled patch connectivity for topology modifications.
1242 // - Send and receive sub meshes for processor patches.
1243 // - Made static because this function may be run in a separate thread.
1244 void dynamicTopoFvMesh::initCoupledConnectivity
1249 // Recast the argument
1250 dynamicTopoFvMesh *mesh = static_cast<dynamicTopoFvMesh*>(argument);
1252 // Identify coupled patches.
1253 if (mesh->identifyCoupledPatches())
1258 // Build and send patch sub-meshes.
1259 mesh->buildProcessorPatchMeshes();
1261 // Build inital maps for locally coupled patches.
1262 mesh->buildLocalCoupledMaps();
1264 // Build maps for coupled processor patches.
1265 mesh->buildProcessorCoupledMaps();
1269 // Move coupled subMesh points
1270 void dynamicTopoFvMesh::moveCoupledSubMeshes()
1272 if (!Pstream::parRun())
1279 Info<< " void dynamicTopoFvMesh::moveCoupledSubMeshes() :"
1280 << " Moving points for coupled subMeshes."
1284 forAll(procIndices_, pI)
1286 label proc = procIndices_[pI];
1288 const coupledMesh& sPM = sendMeshes_[pI];
1289 const coupledMesh& rPM = recvMeshes_[pI];
1291 // Fetch the coupleMap
1292 const coupleMap& scMap = sPM.map();
1293 const coupleMap& rcMap = rPM.map();
1295 Map<label>& pointMap = scMap.entityMap(coupleMap::POINT);
1297 // Fill point buffers
1298 pointField& pBuffer = scMap.pointBuffer();
1299 pointField& opBuffer = scMap.oldPointBuffer();
1301 forAllConstIter(Map<label>, pointMap, pIter)
1303 pBuffer[pIter.key()] = points_[pIter()];
1304 opBuffer[pIter.key()] = oldPoints_[pIter()];
1307 // Buffers have already been allocated
1308 // to the right size, so just transfer points
1310 // Send point buffers to neighbour
1311 meshOps::pWrite(proc, scMap.pointBuffer());
1312 meshOps::pWrite(proc, scMap.oldPointBuffer());
1314 // Receive point buffers from neighbour
1315 meshOps::pRead(proc, rcMap.pointBuffer());
1316 meshOps::pRead(proc, rcMap.oldPointBuffer());
1319 // Wait for transfers to complete before moving on
1320 meshOps::waitForBuffers();
1322 // Set points in mesh
1323 forAll(procIndices_, pI)
1325 // Fetch non-const reference to patchSubMesh
1326 coupledMesh& rPM = recvMeshes_[pI];
1328 // Fetch the coupleMap
1329 const coupleMap& rcMap = rPM.map();
1331 dynamicTopoFvMesh& mesh = rPM.subMesh();
1332 const polyBoundaryMesh& boundary = mesh.boundaryMesh();
1334 // Set points / oldPoints
1335 mesh.points_ = rcMap.pointBuffer();
1336 mesh.oldPoints_ = rcMap.oldPointBuffer();
1338 labelList patchSizes(boundary.size(), -1);
1339 labelList patchStarts(boundary.size(), -1);
1341 forAll(boundary, patchI)
1343 patchSizes[patchI] = boundary[patchI].size();
1344 patchStarts[patchI] = boundary[patchI].start();
1347 // Reset underlying mesh.
1348 // - Use null lists for addressing to avoid over-writes
1349 // - Specify non-valid boundary to avoid globalData creation
1350 mesh.resetPrimitives
1352 xferCopy(rcMap.oldPointBuffer()),
1353 (*reinterpret_cast<Xfer<faceList>*>(0)),
1354 (*reinterpret_cast<Xfer<labelList>*>(0)),
1355 (*reinterpret_cast<Xfer<labelList>*>(0)),
1364 // Insert the cells around the coupled master entity to the mesh
1365 // - Returns a changeMap with a type specifying:
1366 // 1: Insertion was successful
1367 // -1: Insertion failed
1368 // -2: Failed because entity was being handled elsewhere
1369 // - The changeMap index specifies the converted mIndex.
1370 const changeMap dynamicTopoFvMesh::insertCells(const label mIndex)
1372 // Prepare the changeMaps
1375 // Maintain face counts for each inserted cell
1376 Map<label> nCellFaces;
1378 // Track inserted entities
1379 List<Map<label> > pointsToInsert(procIndices_.size());
1380 List<Map<label> > edgesToInsert(procIndices_.size());
1381 List<Map<label> > facesToInsert(procIndices_.size());
1382 List<Map<label> > cellsToInsert(procIndices_.size());
1384 // Track converted entities
1385 List<Map<label> > edgesToConvert(procIndices_.size());
1386 List<Map<label> > facesToConvert(procIndices_.size());
1388 // Track edge / face patch information
1389 List<Map<label> > masterConvertEdgePatch(procIndices_.size());
1390 List<Map<label> > masterConvertFacePatch(procIndices_.size());
1392 Map<label> createPatch;
1393 labelList slaveConvertPatch(procIndices_.size(), -1);
1394 List<Map<Pair<point> > > convertPatchPoints(procIndices_.size());
1396 // First check to ensure that this case can be handled
1397 forAll(procIndices_, pI)
1399 const label cplEnum = is2D() ? coupleMap::FACE : coupleMap::EDGE;
1401 // Fetch reference to maps
1402 const coupleMap& cMap = recvMeshes_[pI].map();
1404 // Does a coupling exist?
1405 if (cMap.findSlave(cplEnum, mIndex) == -1)
1410 // If this entity is being handled elsewhere, bail out
1411 if (priority(procIndices_[pI], lessOp<label>(), Pstream::myProcNo()))
1419 const polyBoundaryMesh& boundary = boundaryMesh();
1421 // Agglomerate cells from surrounding subMeshes
1422 forAll(procIndices_, pI)
1424 const label cplEnum = is2D() ? coupleMap::FACE : coupleMap::EDGE;
1426 // Fetch reference to maps
1427 const coupleMap& cMap = recvMeshes_[pI].map();
1428 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
1432 if ((sIndex = cMap.findSlave(cplEnum, mIndex)) == -1)
1438 Map<label>& procCellMap = cellsToInsert[pI];
1442 // Insert the owner cell
1443 procCellMap.insert(mesh.owner_[sIndex], -1);
1447 // Insert all cells connected to this edge
1448 const labelList& eFaces = mesh.edgeFaces_[sIndex];
1450 forAll(eFaces, faceI)
1452 const label own = mesh.owner_[eFaces[faceI]];
1453 const label nei = mesh.neighbour_[eFaces[faceI]];
1455 if (!procCellMap.found(own))
1457 procCellMap.insert(own, -1);
1460 if (!procCellMap.found(nei) && (nei != -1))
1462 procCellMap.insert(nei, -1);
1467 // Find a patch that talks to this processor
1468 label neiProcPatch = -1;
1470 forAll(boundary, patchI)
1472 if (getNeighbourProcessor(patchI) == procIndices_[pI])
1474 neiProcPatch = patchI;
1479 // Prepare information about inserted / converted faces and edges
1480 const polyBoundaryMesh& slaveBoundary = mesh.boundaryMesh();
1482 // Set the slave conversion patch
1483 forAll(slaveBoundary, patchI)
1485 if (mesh.getNeighbourProcessor(patchI) == Pstream::myProcNo())
1487 slaveConvertPatch[pI] = patchI;
1492 forAllIter(Map<label>, procCellMap, cIter)
1494 const cell& checkCell = mesh.cells_[cIter.key()];
1496 forAll(checkCell, fI)
1498 label mfIndex = -1, sfIndex = checkCell[fI];
1499 label sfPatch = mesh.whichPatch(sfIndex);
1501 // Check whether a face mapping exists for this face
1502 if ((mfIndex = cMap.findMaster(coupleMap::FACE, sfIndex)) > -1)
1504 // Mark as a candidate for conversion
1505 facesToInsert[pI].insert(sfIndex, mfIndex);
1506 facesToConvert[pI].insert(sfIndex, mfIndex);
1510 // Mark for insertion. Some faces may be duplicated
1511 // across processors, but this will be dealt
1512 // with at a later stage.
1513 facesToInsert[pI].insert(sfIndex, -1);
1516 // Check face points for coupling
1517 const face& sFace = mesh.faces_[sfIndex];
1519 // Configure points which need to be inserted
1520 forAll(sFace, pointI)
1522 label sFacePoint = sFace[pointI];
1524 // Skip if added already
1525 if (pointsToInsert[pI].found(sFacePoint))
1530 // Check for coupled points
1531 pointsToInsert[pI].insert
1542 // Loop through edges and check whether edge-mapping exists
1543 const labelList& sfEdges = mesh.faceEdges_[sfIndex];
1545 forAll(sfEdges, edgeI)
1547 const label seIndex = sfEdges[edgeI];
1548 const edge& sEdge = mesh.edges_[seIndex];
1550 // Meshes in 2D don't have edge-mapping, so check
1551 // point maps instead. If either point doesn't exist,
1552 // this is an edge that needs to be inserted.
1553 label cMs = cMap.findMaster(coupleMap::POINT, sEdge[0]);
1554 label cMe = cMap.findMaster(coupleMap::POINT, sEdge[1]);
1556 if (!edgesToInsert[pI].found(seIndex))
1558 edgesToInsert[pI].insert(seIndex, -1);
1561 // If both points have maps, this is a conversion edge
1562 if (cMs > -1 && cMe > -1)
1564 if (!edgesToConvert[pI].found(seIndex))
1566 edgesToConvert[pI].insert(seIndex, -1);
1570 // Add a patch entry as well
1571 if (masterConvertEdgePatch[pI].found(seIndex))
1576 label edgePatch = -1;
1577 label sePatch = mesh.whichEdgePatch(seIndex);
1578 label neiProcNo = mesh.getNeighbourProcessor(sePatch);
1583 // Slave edge was an interior one
1584 if (neiProcPatch > -1)
1586 edgePatch = neiProcPatch;
1590 neiProcNo = procIndices_[pI];
1595 << " No face contact with"
1596 << " processor: " << neiProcNo
1600 label pIdx = sendMeshes_[pI].map().patchIndex();
1602 // Add a patch creation order, if necessary
1603 if (pIdx == -1 && !createPatch.found(neiProcNo))
1605 createPatch.insert(neiProcNo, -1);
1608 // Specify a value for edgePatch:
1609 // - Specify a value that we can use
1610 // to back out the patch after creation
1611 // - Also needs to bypass failure check
1612 edgePatch = (-2 - neiProcNo);
1616 if (sePatch == (slaveBoundary.size() - 1))
1618 // The 'defaultPatch'
1619 edgePatch = neiProcPatch;
1624 if (neiProcNo == Pstream::myProcNo())
1626 // Set the master edge patch
1627 edgePatch = neiProcPatch;
1631 // If this other processor is
1632 // lesser-ranked, bail out.
1651 << " Edge: " << seIndex
1652 << " :: " << mesh.edges_[seIndex]
1653 << " is talking to processor: "
1654 << neiProcNo << endl;
1657 // Find an appropriate boundary patch
1658 forAll(boundary, patchI)
1660 label eP = getNeighbourProcessor(patchI);
1662 if (eP == neiProcNo)
1669 if (edgePatch == -1)
1674 << " No direct contact with"
1675 << " processor: " << neiProcNo
1676 << " so adding to: " << neiProcPatch
1680 // Add this to neiProcPatch instead.
1681 edgePatch = neiProcPatch;
1688 edgePatch = sePatch;
1691 if (edgePatch == -1)
1693 // Write out the edge
1694 mesh.writeVTK("seEdge", seIndex, 1);
1696 Pout<< " Could not find correct patch info: " << nl
1697 << " sEdge: " << sEdge << nl
1698 << " seIndex: " << seIndex << nl
1699 << " sePatch: " << sePatch << nl
1700 << " neiProcPatch: " << neiProcPatch << nl
1701 << " neiProcNo: " << neiProcNo << nl
1702 << " proc: " << procIndices_[pI] << nl
1703 << " cMs: " << cMs << " cMe: " << cMe << nl
1704 << " Patch Name: " <<
1707 slaveBoundary[sePatch].name() :
1710 << abort(FatalError);
1713 // Add the patch entry
1714 masterConvertEdgePatch[pI].insert
1722 if (masterConvertFacePatch[pI].found(sfIndex))
1727 label facePatch = -1;
1728 label neiProcNo = mesh.getNeighbourProcessor(sfPatch);
1732 // Slave face was an interior one
1733 if (neiProcPatch > -1)
1735 facePatch = neiProcPatch;
1739 neiProcNo = procIndices_[pI];
1744 << " No face contact with"
1745 << " processor: " << neiProcNo
1749 label pIdx = sendMeshes_[pI].map().patchIndex();
1751 // Add a patch creation order, if necessary
1752 if (pIdx == -1 && !createPatch.found(neiProcNo))
1754 createPatch.insert(neiProcNo, -1);
1757 // Specify a value for facePatch:
1758 // - Specify a value that we can use
1759 // to back out the patch after creation
1760 // - Also needs to bypass failure check
1761 facePatch = (-2 - neiProcNo);
1765 if (sfPatch == (slaveBoundary.size() - 1))
1767 // The 'defaultPatch'
1768 facePatch = neiProcPatch;
1773 if (neiProcNo == Pstream::myProcNo())
1775 // Check to see if this face was converted before
1780 cMap.entityOperations(),
1781 coupleMap::CONVERT_PATCH
1785 // Loop through points and check for match
1786 const pointField& pts = cMap.moveNewPoints();
1787 const face& fCheck = mesh.faces_[sfIndex];
1788 const point fC = fCheck.centre(mesh.points_);
1790 scalar tol = mag(fC - mesh.points_[fCheck[0]]);
1794 scalar dist = mag(fC - pts[ptI]);
1796 if (dist < (geomMatchTol_() * tol))
1798 // Face was converted before
1802 << " Face: " << sfIndex
1803 << " :: " << mesh.faces_[sfIndex]
1804 << " was converted before."
1815 // Set the master face patch
1816 facePatch = neiProcPatch;
1820 // If this other processor is
1821 // lesser-ranked, bail out.
1840 << " Face: " << sfIndex
1841 << " :: " << mesh.faces_[sfIndex]
1842 << " is talking to processor: "
1843 << neiProcNo << endl;
1846 // Find an appropriate boundary patch
1847 forAll(boundary, patchI)
1849 label fP = getNeighbourProcessor(patchI);
1851 if (fP == neiProcNo)
1858 // Specify a convert-patch operation
1859 // for later, if this operation succeeds
1860 const face& sfCheck = mesh.faces_[sfIndex];
1862 point nfC = sfCheck.centre(mesh.points_);
1863 point ofC = sfCheck.centre(mesh.oldPoints_);
1865 // Are we talking to this processor already?
1866 label prI = findIndex(procIndices_, neiProcNo);
1868 // Check for face-contact
1869 if (facePatch == -1)
1874 << " No face contact with"
1875 << " processor: " << neiProcNo
1881 Pout<< " No contact with: " << neiProcNo
1882 << abort(FatalError);
1885 label pIdx = sendMeshes_[prI].map().patchIndex();
1887 // Add a patch creation order, if necessary
1888 if (pIdx == -1 && !createPatch.found(neiProcNo))
1890 createPatch.insert(neiProcNo, -1);
1893 // Specify a value for facePatch:
1894 // - Specify a value that we can use
1895 // to back out the patch after creation
1896 // - Also needs to bypass failure check
1897 facePatch = (-2 - neiProcNo);
1901 convertPatchPoints[prI].insert
1904 Pair<point>(nfC, ofC)
1911 facePatch = sfPatch;
1914 if (facePatch == -1)
1916 // Write out the face
1917 mesh.writeVTK("sfFace", sfIndex, 2);
1919 Pout<< " Could not find correct patch info: " << nl
1920 << " sfIndex: " << sfIndex << nl
1921 << " sfPatch: " << sfPatch << nl
1922 << " neiProcPatch: " << neiProcPatch << nl
1923 << " slavePatch: " <<
1926 slaveBoundary[sfPatch].name() :
1929 << abort(FatalError);
1932 // Add the patch entry
1933 masterConvertFacePatch[pI].insert
1942 // Specify a merge tolerance for insertion points
1945 is2D() ? 0.0 : geomMatchTol_() * magSqr(edges_[mIndex].vec(points_))
1948 // Check for point / edge processor connections
1951 const label myProcNo = Pstream::myProcNo();
1953 forAll(procIndices_, pI)
1955 const Map<label>& cEdges = edgesToInsert[pI];
1956 const coupleMap& cMap = recvMeshes_[pI].map();
1957 const Map<labelList>& pEdgeMap = cMap.subMeshEdgeMap();
1958 const Map<labelList>& pPointMap = cMap.subMeshPointMap();
1959 const dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
1961 // Some points may have been inserted by a prior
1962 // operation involving disconnected entities.
1963 // Check if points need to be truly inserted.
1964 Map<label>& cPoints = pointsToInsert[pI];
1966 forAllIter(Map<label>, cPoints, pIter)
1973 label cSp = pIter.key();
1975 Map<labelList>::const_iterator pIt = pPointMap.find(cSp);
1977 if (pIt == pPointMap.end())
1982 const point& pointI = sMesh.points_[pIt.key()];
1984 bool foundMaster = false;
1985 const labelList& procs = pIt();
1987 forAll(procs, procJ)
1989 label pJ = -1, nProcNo = procs[procJ];
1993 (nProcNo == procIndices_[pI]) ||
1994 ((pJ = findIndex(procIndices_, nProcNo)) == -1)
2000 // Fetch map from the other processor,
2001 // and check if a mapping already exists
2002 const coupleMap& cMapJ = recvMeshes_[pJ].map();
2003 const Map<labelList>& pMapJ = cMapJ.subMeshPointMap();
2004 const dynamicTopoFvMesh& sMeshJ = recvMeshes_[pJ].subMesh();
2006 forAllConstIter(Map<labelList>, pMapJ, pjIter)
2017 if (pointMaster == -1)
2022 const point& pointJ = sMeshJ.points_[pjIter.key()];
2024 if (magSqr(pointI - pointJ) < mergeTol)
2028 Pout<< " Using mapped point: " << pjIter.key()
2030 << " procs: " << pjIter()
2031 << " master: " << pointMaster
2035 // Assign to existing master point
2036 pIter() = pointMaster;
2038 // Update maps for the point
2042 pIter(), pIter.key()
2048 pIter.key(), pIter()
2063 forAllConstIter(Map<label>, cEdges, eIter)
2065 label cSe = eIter.key();
2067 Map<labelList>::const_iterator eIt = pEdgeMap.find(cSe);
2069 if (eIt != pEdgeMap.end())
2071 const labelList& procs = eIt();
2073 forAll(procs, procI)
2075 label nProcNo = procs[procI];
2077 if (priority(nProcNo, greaterEqOp<label>(), myProcNo))
2086 << " :: " << sMesh.edges_[cSe]
2087 << " is talking to processors: " << procs
2088 << " so bailing out."
2098 const edge& checkEdge = sMesh.edges_[cSe];
2100 forAll(checkEdge, pointI)
2102 label cSp = checkEdge[pointI];
2104 Map<labelList>::const_iterator pIt = pPointMap.find(cSp);
2106 if (pIt == pPointMap.end())
2111 const labelList& procs = pIt();
2113 forAll(procs, procI)
2115 label nProcNo = procs[procI];
2117 if (priority(nProcNo, greaterEqOp<label>(), myProcNo))
2125 << " Points: " << nl
2126 << " Edge[0]: " << checkEdge[0] << nl
2127 << " Edge[1]: " << checkEdge[1] << nl
2128 << " Processor conflict: " << nl
2129 << " Point: " << cSp << nl
2130 << " Processors: " << procs
2143 // Perform a few debug calls before modifications
2147 << "Inserting cell(s) around coupled "
2148 << (is2D() ? "face: " : "edge: ") << mIndex
2152 // Check to see if any new processor
2153 // patches need to be created
2154 forAllIter(Map<label>, createPatch, procIter)
2156 label neiProcNo = procIter.key();
2158 // Create a new processor patch
2159 procIter() = createProcessorPatch(neiProcNo);
2161 // Renumber edge conversion patches
2162 forAll(masterConvertEdgePatch, pI)
2164 Map<label>& convertMap = masterConvertEdgePatch[pI];
2166 forAllIter(Map<label>, convertMap, mIter)
2170 // Back out the neighbouring processor ID
2171 label proc = Foam::mag(mIter() + 2);
2173 if (proc == neiProcNo)
2176 mIter() = procIter();
2182 // Renumber face conversion patches
2183 forAll(masterConvertFacePatch, pI)
2185 Map<label>& convertMap = masterConvertFacePatch[pI];
2187 forAllIter(Map<label>, convertMap, mIter)
2191 // Back out the neighbouring processor ID
2192 label proc = Foam::mag(mIter() + 2);
2194 if (proc == neiProcNo)
2197 mIter() = procIter();
2203 // Find index in processor list
2204 label pI = findIndex(procIndices_, neiProcNo);
2208 Pout<< " Could not find index for processor: " << neiProcNo
2209 << " in indices: " << procIndices_
2210 << abort(FatalError);
2213 // Create patch on subMesh
2214 dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2216 slaveConvertPatch[pI] = mesh.createProcessorPatch(Pstream::myProcNo());
2219 // Build a list of mapping entities on this processor
2220 dynamicLabelList mapCells(10);
2224 label own = owner_[mIndex];
2226 mapCells.append(own);
2230 const labelList& eFaces = edgeFaces_[mIndex];
2232 forAll(eFaces, faceI)
2234 label own = owner_[eFaces[faceI]];
2235 label nei = neighbour_[eFaces[faceI]];
2237 if (findIndex(mapCells, own) == -1)
2239 mapCells.append(own);
2244 if (findIndex(mapCells, nei) == -1)
2246 mapCells.append(nei);
2252 // Loop through insertion cells and
2253 // create an equivalent on this mesh
2254 forAll(procIndices_, pI)
2256 const label cplEnum = is2D() ? coupleMap::FACE : coupleMap::EDGE;
2258 // Fetch reference to maps
2259 const coupleMap& cMap = recvMeshes_[pI].map();
2260 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2264 if ((sIndex = cMap.findSlave(cplEnum, mIndex)) == -1)
2270 Map<label>& procCellMap = cellsToInsert[pI];
2272 forAllIter(Map<label>, procCellMap, cIter)
2274 const cell& checkCell = mesh.cells_[cIter.key()];
2276 scalar sLengthScale = -1.0;
2278 if (edgeRefinement_)
2280 sLengthScale = mesh.lengthScale_[cIter.key()];
2283 // Add an empty cell for now, and update
2284 // with face information at a later stage.
2285 cIter() = insertCell(cell(checkCell.size()), sLengthScale);
2287 // Initialize a face counter
2288 nCellFaces.insert(cIter(), 0);
2292 Pout<< " Map cell: " << cIter()
2293 << " for cell: " << cIter.key()
2294 << " on proc: " << procIndices_[pI]
2298 // Add this cell to the map.
2299 map.addCell(cIter(), mapCells);
2301 // Set basic mapping for this cell
2302 setCellMapping(cIter(), mapCells);
2305 // Push operation for the slave into coupleMap
2309 coupleMap::REMOVE_CELL
2313 // Build a list of edges that need to be converted to interior.
2314 // - Do this by looking at edges of master face conversion candidates.
2315 // - Some edges may not need conversion, but deal with this later.
2316 forAll(facesToConvert, pI)
2318 // Fetch reference to subMesh
2319 const coupleMap& cMap = recvMeshes_[pI].map();
2320 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2322 const Map<label>& procFaceMap = facesToConvert[pI];
2324 forAllConstIter(Map<label>, procFaceMap, fIter)
2326 const labelList& mfEdges = faceEdges_[fIter()];
2327 const labelList& sfEdges = mesh.faceEdges_[fIter.key()];
2329 forAll(sfEdges, edgeI)
2331 if (edgesToInsert[pI][sfEdges[edgeI]] != -1)
2333 // Already mapped this edge. Move on.
2337 // Configure the comparison edge
2338 const edge& sEdge = mesh.edges_[sfEdges[edgeI]];
2340 label cMs = cMap.findMaster(coupleMap::POINT, sEdge[0]);
2341 label cMe = cMap.findMaster(coupleMap::POINT, sEdge[1]);
2343 edge cEdge(cMs, cMe);
2345 forAll(mfEdges, edgeJ)
2347 const edge& mEdge = edges_[mfEdges[edgeJ]];
2351 edgesToInsert[pI][sfEdges[edgeI]] = mfEdges[edgeJ];
2352 edgesToConvert[pI][sfEdges[edgeI]] = mfEdges[edgeJ];
2360 // Insert all points, with merging if necessary
2361 forAll(pointsToInsert, pI)
2363 Map<label>& procPointMapI = pointsToInsert[pI];
2365 // Fetch reference to subMesh
2366 const coupleMap& cMapI = recvMeshes_[pI].map();
2367 const dynamicTopoFvMesh& meshI = recvMeshes_[pI].subMesh();
2369 forAllIter(Map<label>, procPointMapI, pItI)
2376 const point& pointI = meshI.points_[pItI.key()];
2378 bool foundMerge = false;
2380 forAll(pointsToInsert, pJ)
2387 Map<label>& procPointMapJ = pointsToInsert[pJ];
2389 // Fetch reference to subMesh
2390 const coupleMap& cMapJ = recvMeshes_[pJ].map();
2391 const dynamicTopoFvMesh& meshJ = recvMeshes_[pJ].subMesh();
2393 // Compare points with this processor
2394 forAllIter(Map<label>, procPointMapJ, pItJ)
2396 const point& pointJ = meshJ.points_[pItJ.key()];
2398 if (magSqr(pointI - pointJ) < mergeTol)
2402 // Make a merge entry
2403 label mergePointIndex =
2408 meshJ.oldPoints_[pItJ.key()],
2413 pItI() = mergePointIndex;
2414 pItJ() = mergePointIndex;
2416 // Update maps for the new point
2443 Pout<< " Map point: " << mergePointIndex
2444 << " pointI: " << pItI.key()
2445 << " procI: " << procIndices_[pI]
2446 << " pointJ: " << pItJ.key()
2447 << " procJ: " << procIndices_[pJ]
2451 // Add this point to the map.
2452 map.addPoint(mergePointIndex);
2456 // Point appears to have been inserted
2457 // by a previous operation. Map to it.
2460 Pout<< " Inserted point: " << pItJ()
2461 << " pointI: " << pItI.key()
2462 << " procI: " << procIndices_[pI]
2463 << " pointJ: " << pItJ.key()
2464 << " procJ: " << procIndices_[pJ]
2471 // Update maps for the point
2496 // Add a new unique point
2497 label newPointIndex =
2502 meshI.oldPoints_[pItI.key()],
2508 pItI() = newPointIndex;
2510 // Update maps for the new point
2511 cMapI.mapSlave(coupleMap::POINT, pItI(), pItI.key());
2512 cMapI.mapMaster(coupleMap::POINT, pItI.key(), pItI());
2516 Pout<< " Map point: " << newPointIndex
2517 << " for point: " << pItI.key()
2518 << " on proc: " << procIndices_[pI]
2522 // Add this point to the map.
2523 map.addPoint(newPointIndex);
2527 // Occassionally, inserted edges may already be present.
2528 // Ensure that edges are not added twice.
2531 forAll(facesToInsert, pI)
2533 // Fetch reference to subMesh
2534 const coupleMap& cMap = recvMeshes_[pI].map();
2535 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2537 const Map<label>& procFaceMap = facesToInsert[pI];
2539 forAllConstIter(Map<label>, procFaceMap, fIter)
2541 const labelList& sfEdges = mesh.faceEdges_[fIter.key()];
2543 forAll(sfEdges, edgeI)
2545 label seIndex = sfEdges[edgeI];
2547 if (edgesToInsert[pI][seIndex] != -1)
2549 // Already mapped this edge. Move on.
2553 label cMe = cMap.findMaster(coupleMap::EDGE, seIndex);
2559 Pout<< " Found master edge: " << cMe
2560 << " :: " << edges_[cMe]
2561 << " for slave edge: " << seIndex
2562 << " :: " << mesh.edges_[seIndex]
2563 << " on proc: " << procIndices_[pI]
2567 edgesToInsert[pI][seIndex] = cMe;
2568 edgesToConvert[pI][seIndex] = cMe;
2573 // Check for point-only coupling
2574 const edge& sEdge = mesh.edges_[seIndex];
2578 cMap.findMaster(coupleMap::POINT, sEdge[0]),
2579 cMap.findMaster(coupleMap::POINT, sEdge[1])
2582 if (cEdge[0] > -1 && cEdge[1] > -1)
2586 // Look at pointEdges info for a boundary edge
2587 const labelList& pEdges = pointEdges_[cEdge[0]];
2589 forAll(pEdges, edgeJ)
2591 if (edges_[pEdges[edgeJ]] == cEdge)
2593 if (whichEdgePatch(pEdges[edgeJ]) > -1)
2595 meIndex = pEdges[edgeJ];
2605 Pout<< " Found master edge: " << meIndex
2606 << " :: " << edges_[meIndex]
2607 << " for slave edge: " << seIndex
2608 << " :: " << mesh.edges_[seIndex]
2609 << " on proc: " << procIndices_[pI]
2613 edgesToInsert[pI][seIndex] = meIndex;
2615 if (edgesToConvert[pI].found(seIndex))
2617 edgesToConvert[pI][seIndex] = meIndex;
2621 edgesToConvert[pI].insert(seIndex, meIndex);
2630 // Write out some debug info
2633 forAll(cellsToInsert, pI)
2635 label procIndex = procIndices_[pI];
2637 // Fetch reference to subMesh
2638 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2640 // Define a convenience output macro
2641 # define writeOutputVTK(fName, eMap, eType, insert) \
2643 labelList entities(eMap.size()); \
2647 forAllConstIter(Map<label>, eMap, eIter) \
2649 if (eIter() == -1 && insert) \
2651 entities[nEnt++] = eIter.key(); \
2654 if (eIter() > -1 && !insert) \
2656 entities[nEnt++] = eIter.key(); \
2660 entities.setSize(nEnt); \
2662 if (entities.size()) \
2667 + Foam::name(mIndex) + '_' \
2668 + Foam::name(procIndex), \
2670 eType, false, true \
2675 writeOutputVTK("edgesToConvert", edgesToConvert[pI], 1, false)
2676 writeOutputVTK("facesToConvert", facesToConvert[pI], 2, false)
2677 writeOutputVTK("pointsToConvert", pointsToInsert[pI], 0, false)
2678 writeOutputVTK("pointsToInsert", pointsToInsert[pI], 0, true)
2679 writeOutputVTK("edgesToInsert", edgesToInsert[pI], 1, true)
2680 writeOutputVTK("facesToInsert", facesToInsert[pI], 2, true)
2681 writeOutputVTK("cellsToInsert", cellsToInsert[pI], 3, false)
2685 // Add edges to the mesh, noting that all
2686 // required points have already been added.
2687 forAll(edgesToInsert, pI)
2689 // Fetch reference to subMesh
2690 const coupleMap& cMap = recvMeshes_[pI].map();
2691 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2693 Map<label>& procEdgeMap = edgesToInsert[pI];
2695 forAllIter(Map<label>, procEdgeMap, eIter)
2697 // Skip if the edge is a conversion candidate
2703 const edge& sEdge = mesh.edges_[eIter.key()];
2707 cMap.findMaster(coupleMap::POINT, sEdge[0]),
2708 cMap.findMaster(coupleMap::POINT, sEdge[1])
2711 // Ensure that this edge hasn't been added before
2713 bool foundDuplicate = false;
2715 const List<objectMap>& addedEdges = map.addedEdgeList();
2717 forAll(addedEdges, edgeI)
2719 neIndex = addedEdges[edgeI].index();
2721 if (edges_[neIndex] == newEdge)
2723 foundDuplicate = true;
2730 // Note the duplicate index for later
2736 // Insert edge with null edgeFaces for now.
2737 // This can be corrected later.
2742 masterConvertEdgePatch[pI][eIter.key()],
2750 Pout<< " Map edge: " << eIter() << "::" << newEdge
2751 << " for edge: " << eIter.key() << "::" << sEdge
2755 // Add this edge to the map.
2756 map.addEdge(eIter());
2760 // Add faces to the mesh, noting that all required
2761 // points and edges have already been added.
2762 forAll(facesToInsert, pI)
2764 // Fetch reference to subMesh and coupleMap
2765 const coupleMap& cMapI = recvMeshes_[pI].map();
2766 const dynamicTopoFvMesh& meshI = recvMeshes_[pI].subMesh();
2768 Map<label>& procFaceMap = facesToInsert[pI];
2770 forAllIter(Map<label>, procFaceMap, fI)
2772 // Skip if the face is a conversion candidate
2778 const face& sFaceI = meshI.faces_[fI.key()];
2779 const labelList& sfEdges = meshI.faceEdges_[fI.key()];
2781 // Configure new / comparison faces
2782 face nF(sFaceI.size()), cF(sFaceI.size());
2783 labelList nFaceEdges(sfEdges.size());
2785 // Configure points from map
2798 // This may be a face shared by two other processors.
2799 // Ensure that duplicates are added only once.
2800 label dupeOwnCell = -1;
2801 bool foundDuplicate = false, foundInsertedDuplicate = false;
2803 forAll(facesToInsert, pJ)
2810 // Fetch reference to subMesh and coupleMap
2811 const coupleMap& cMapJ = recvMeshes_[pJ].map();
2812 const dynamicTopoFvMesh& meshJ = recvMeshes_[pJ].subMesh();
2814 const Map<label>& procFaceMapJ = facesToInsert[pJ];
2816 forAllConstIter(Map<label>, procFaceMapJ, fJ)
2818 const face& sFaceJ = meshJ.faces_[fJ.key()];
2820 // Discard dissimilar face sizes
2821 if (sFaceJ.size() != cF.size())
2826 // Configure points from map
2841 // Optimized triangular face comparison
2846 triFace(nF[0], nF[1], nF[2]),
2847 triFace(cF[0], cF[1], cF[2])
2851 foundDuplicate = true;
2856 // Regular face compare
2857 if (face::compare(nF, cF))
2859 foundDuplicate = true;
2865 // Record the owner for posterity
2868 cellsToInsert[pJ][meshJ.owner_[fJ.key()]]
2871 // Was the duplicate face inserted before this?
2874 // Note the duplicate index for later
2877 // If patch conversion entries were made,
2878 // remove them as well
2881 convertPatchPoints[pI].found(fJ.key())
2882 && convertPatchPoints[pJ].found(fI.key())
2885 convertPatchPoints[pI].erase(fJ.key());
2886 convertPatchPoints[pJ].erase(fI.key());
2889 foundInsertedDuplicate = true;
2902 // Face was inserted before, so don't insert again
2903 if (foundInsertedDuplicate)
2908 // Configure edges from edgesToInsert
2909 forAll(sfEdges, edgeI)
2913 // Configure with the appropriate edge
2914 if (edgesToInsert[pI].found(sfEdges[edgeI]))
2916 meIndex = edgesToInsert[pI][sfEdges[edgeI]];
2921 // Something is wrong here.
2922 Pout<< " Could not find correspondence for slave edge: "
2924 << ":: " << meshI.edges_[sfEdges[edgeI]]
2925 << nl << " mIndex: " << mIndex
2926 << abort(FatalError);
2929 nFaceEdges[edgeI] = meIndex;
2932 // Determine patch, owner and neighbour for this face
2933 label nPatch = -1, nOwner = -1, nNeighbour = -1;
2935 const polyBoundaryMesh& slaveBoundary = meshI.boundaryMesh();
2937 label sfPatch = meshI.whichPatch(fI.key());
2938 label sFaceOwn = meshI.owner_[fI.key()];
2939 label sFaceNei = meshI.neighbour_[fI.key()];
2943 cellsToInsert[pI].found(sFaceOwn) ?
2944 cellsToInsert[pI][sFaceOwn] : -1
2949 cellsToInsert[pI].found(sFaceNei) ?
2950 cellsToInsert[pI][sFaceNei] : -1
2953 // If a duplicate face was found, over-ride neighbour.
2954 // Face-flipping will be taken care of automatically.
2957 mFaceNei = dupeOwnCell;
2960 if (mFaceOwn != -1 && mFaceNei == -1)
2962 // Boundary face already has correct orientation
2967 nPatch = masterConvertFacePatch[pI][fI.key()];
2970 if (mFaceOwn == -1 && mFaceNei != -1)
2972 // Boundary face is inverted. Flip it
2973 nF = nF.reverseFace();
2978 nPatch = masterConvertFacePatch[pI][fI.key()];
2981 if (mFaceOwn != -1 && mFaceNei != -1)
2983 // Interior face. Check if a flip is necessary.
2984 if (mFaceNei < mFaceOwn)
2986 nF = nF.reverseFace();
2988 nNeighbour = mFaceOwn;
2993 nNeighbour = mFaceNei;
2999 if (mFaceOwn == -1 && mFaceNei == -1)
3001 // Something is wrong here.
3002 Pout<< "Could not find correct owner / neighbour info: " << nl
3003 << " Face: " << nF << nl
3004 << " Owner: " << mFaceOwn << nl
3005 << " Neighbour: " << mFaceNei << nl
3006 << " - Slave Face: " << sFaceI << nl
3007 << " - Slave Patch: " << slaveBoundary[sfPatch].name() << nl
3008 << " - Slave Owner: " << sFaceOwn << nl
3009 << " - Slave Neighbour: " << sFaceNei << nl
3010 << abort(FatalError);
3013 // Insert the new face
3026 // Size up edgeFaces for each edge
3027 forAll(nFaceEdges, edgeI)
3032 edgeFaces_[nFaceEdges[edgeI]]
3038 Pout<< " Map face: " << fI() << "::" << nF
3039 << " Own: " << nOwner << " Nei: " << nNeighbour
3040 << " fE: " << nFaceEdges << nl
3041 << " for face: " << fI.key() << "::" << sFaceI
3042 << " Own: " << sFaceOwn << " Nei: " << sFaceNei
3043 << " fE: " << sfEdges
3047 // Add this face to the map.
3050 if (nPatch > -1 && getNeighbourProcessor(nPatch) == -1)
3052 // Physical patch on subMesh.
3053 // - Set an invalid number so that
3054 // an entry is made in facesFromFaces,
3055 // while faceParents is empty.
3056 setFaceMapping(fI(), labelList(1, -1));
3060 // Interior / processor face
3061 setFaceMapping(fI());
3065 cells_[nOwner][nCellFaces[nOwner]++] = fI();
3067 if (nNeighbour > -1)
3069 cells_[nNeighbour][nCellFaces[nNeighbour]++] = fI();
3074 // Loop through conversion faces
3075 forAll(facesToConvert, pI)
3077 // Fetch reference to subMesh
3078 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
3080 const Map<label>& procFaceMap = facesToConvert[pI];
3082 forAllConstIter(Map<label>, procFaceMap, fIter)
3084 label fIndex = fIter();
3086 // Fetch the owner information
3087 label mFaceOwn = owner_[fIndex];
3088 label sFaceOwn = mesh.owner_[fIter.key()];
3089 label newNeighbour = cellsToInsert[pI][sFaceOwn];
3091 // Insert a new internal face.
3092 // - Orientation should be correct, because this is a boundary
3093 // face converted to an interior, and adjacent to an added cell.
3094 // - Add the new faceEdges from the existing face. This may contain
3095 // edges that need to be converted, but that will be done later.
3096 label newFaceIndex =
3101 face(faces_[fIndex]),
3104 labelList(faceEdges_[fIndex])
3109 map.addFace(newFaceIndex, labelList(1, fIndex));
3111 // Set mapping information
3112 setFaceMapping(newFaceIndex);
3114 // Update the owner cell
3115 meshOps::replaceLabel
3122 // Update the neighbour cell
3123 cells_[newNeighbour][nCellFaces[newNeighbour]++] = newFaceIndex;
3125 // Replace edgeFaces with the new face index
3126 const labelList& fEdges = faceEdges_[newFaceIndex];
3128 forAll(fEdges, edgeI)
3130 meshOps::replaceLabel
3134 edgeFaces_[fEdges[edgeI]]
3138 // Remove the old boundary face
3142 map.removeFace(fIndex);
3144 // For 2D meshes, the boundary face gets converted
3145 // to an interior one. Note the index for further operations.
3146 if ((mIndex == fIndex) && is2D())
3148 map.index() = newFaceIndex;
3153 // Sequentially add any convert-patch operations
3154 forAll(convertPatchPoints, pI)
3156 // Fetch reference to maps / mesh
3157 const coupleMap& cMap = recvMeshes_[pI].map();
3158 dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
3160 if (convertPatchPoints[pI].empty())
3165 // Find the appropriate processor patch
3166 // for patch conversion operations
3167 label procPatch = -1;
3169 for (label patchI = 0; patchI < mesh.nPatches_; patchI++)
3171 label neiProc = mesh.getNeighbourProcessor(patchI);
3173 if (neiProc == Pstream::myProcNo())
3180 if (procPatch == -1)
3182 Pout<< " * * * insertCells() * * * " << nl
3183 << " Could not find patch on slave processor: "
3184 << procIndices_[pI] << nl
3185 << " convertPatchPoints: " << convertPatchPoints[pI]
3186 << abort(FatalError);
3189 forAllConstIter(Map<Pair<point> >, convertPatchPoints[pI], fIter)
3191 // Search slave mesh for faces
3192 const point& newCentre = fIter().first();
3193 const point& oldCentre = fIter().second();
3195 label replaceFace = -1;
3197 // Accumulate stats in case of failure
3198 scalar minDist = GREAT;
3199 vector minPoint = vector::zero;
3200 dynamicLabelList checkedFaces;
3204 // Reserve for append
3205 checkedFaces.setCapacity(50);
3208 // Loop through all boundary faces,
3209 // and compute / compare face centres
3210 label sTot = mesh.faces_.size(), sInt = mesh.nOldInternalFaces_;
3212 for (label faceI = sInt; faceI < sTot; faceI++)
3214 const face& fCheck = mesh.faces_[faceI];
3221 label pIndex = mesh.whichPatch(faceI);
3223 if (mesh.getNeighbourProcessor(pIndex) == -1)
3228 // Compute face-centre
3229 vector fC = fCheck.centre(mesh.points_);
3231 // Compute tolerance
3232 scalar tol = mag(mesh.points_[fCheck[0]] - fC);
3233 scalar dist = mag(fC - newCentre);
3235 if (dist < (geomMatchTol_() * tol))
3237 replaceFace = faceI;
3248 checkedFaces.append(faceI);
3253 // Ensure that the face was found
3254 if (replaceFace == -1)
3259 + Foam::name(fIter.key()),
3264 Pout<< " * * * insertCells() * * * " << nl
3265 << " Convert patch Op failed." << nl
3266 << " Face: " << fIter.key() << nl
3267 << " minPoint: " << minPoint << nl
3268 << " minDistance: " << minDist << nl
3269 << " newCentre: " << newCentre << nl
3270 << " oldCentre: " << oldCentre << nl
3271 << abort(FatalError);
3274 // Obtain a copy before adding the new face,
3275 // since the reference might become
3276 // invalid during list resizing.
3277 // Edges don't have to change, since they're
3278 // all on the boundary anyway.
3279 face newFace = mesh.faces_[replaceFace];
3280 label newOwn = mesh.owner_[replaceFace];
3281 labelList newFaceEdges = mesh.faceEdges_[replaceFace];
3283 label newFaceIndex =
3295 // changeMap update for the new face is not necessary,
3296 // since it is on the slave subMesh
3298 // Set mapping information
3299 mesh.setFaceMapping(newFaceIndex);
3301 meshOps::replaceLabel
3308 // Correct edgeFaces with the new face label.
3309 forAll(newFaceEdges, edgeI)
3311 meshOps::replaceLabel
3315 mesh.edgeFaces_[newFaceEdges[edgeI]]
3321 Pout<< " Pushing CONVERT_PATCH for face: " << replaceFace
3322 << " :: " << mesh.faces_[replaceFace] << nl
3323 << " with new point: " << fIter().first()
3324 << " and old point: " << fIter().second()
3325 << " on proc: " << procIndices_[pI]
3329 // Finally remove the face
3330 mesh.removeFace(replaceFace);
3335 coupleMap::CONVERT_PATCH,
3342 // Loop through conversion edges
3343 forAll(edgesToConvert, pI)
3346 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
3348 const Map<label>& procEdgeMap = edgesToConvert[pI];
3350 // Loop through conversion edges
3351 forAllConstIter(Map<label>, procEdgeMap, eIter)
3353 label eIndex = eIter(), physPatch = -1;
3354 bool allInterior = true, keepEdge = true;
3358 // Not a conversion edge. Search added edge list for index
3359 eIndex = edgesToInsert[pI][eIter.key()];
3363 // Something's wrong
3364 Pout<< " Edge: " << eIter.key()
3365 << " :: " << mesh.edges_[eIter.key()]
3366 << " was never inserted."
3367 << abort(FatalError);
3371 const labelList& eFaces = edgeFaces_[eIndex];
3375 forAll(eFaces, faceI)
3377 if (whichPatch(eFaces[faceI]) > -1)
3379 allInterior = false;
3386 // Edge was deleted before, so skip it
3387 allInterior = false;
3398 // Check if patches need to be switched
3399 label edgePatch = whichEdgePatch(eIndex);
3401 // Is this edge on a processor patch?
3402 bool edgeOnProc = (getNeighbourProcessor(edgePatch) > -1);
3406 bool foundProc = false;
3408 forAll(eFaces, faceI)
3410 label fPatch = whichPatch(eFaces[faceI]);
3417 if (getNeighbourProcessor(fPatch) == -1)
3419 // Note physical patch for later
3424 // Still on a processor patch.
3437 if (edgeOnProc && physPatch > -1)
3439 // Edge needs to be moved to another patch
3443 if ((mIndex == eIndex) && is3D() && map.index() == -1)
3445 // Keep the edge, and note index for later
3447 map.index() = eIndex;
3456 // This edge needs to be converted to interior / other patch
3457 label newEdgeIndex =
3462 edge(edges_[eIndex]),
3463 labelList(edgeFaces_[eIndex])
3468 map.addEdge(newEdgeIndex, labelList(1, eIndex));
3470 // Update faceEdges information for all connected faces
3471 const labelList& neFaces = edgeFaces_[newEdgeIndex];
3473 forAll(neFaces, faceI)
3475 meshOps::replaceLabel
3479 faceEdges_[neFaces[faceI]]
3483 // Remove the old boundary edge
3487 map.removeEdge(eIndex);
3489 // For 3D meshes, the boundary edge gets converted
3490 // to an interior one. Note the index for further operations.
3491 if ((mIndex == eIndex) && is3D())
3493 map.index() = newEdgeIndex;
3496 // Replace the entry in edgesToInsert
3497 edgesToInsert[pI][eIter.key()] = newEdgeIndex;
3501 List<changeMap> slaveMaps(procIndices_.size());
3503 // Loop through all processors, and remove cells
3504 const label emptyMap = -7;
3506 forAll(cellsToInsert, pI)
3508 const Map<label>& procCellMap = cellsToInsert[pI];
3510 if (procCellMap.empty())
3512 // Set type to something recognizable
3513 slaveMaps[pI].type() = emptyMap;
3518 // Prepare a list of cells
3519 const labelList cellsToRemove = procCellMap.toc();
3521 // Fetch non-const reference to subMesh
3522 dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
3524 // Now remove cells for this processor
3530 slaveConvertPatch[pI],
3531 "rc_" + Foam::name(mIndex)
3536 // Now map entities from the removeCells operation
3537 forAll(slaveMaps, pI)
3539 const changeMap& slaveMap = slaveMaps[pI];
3541 // Skip empty entities
3542 if (slaveMap.type() == emptyMap)
3548 const coupleMap& cMap = recvMeshes_[pI].map();
3549 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
3551 // Map faces from patches.
3552 const Map<label>& procFaceMap = facesToInsert[pI];
3554 // Check removed faces
3555 const labelList& rsfList = slaveMap.removedFaceList();
3557 forAll(rsfList, indexI)
3559 label sfIndex = rsfList[indexI], mfIndex = -1;
3561 // Does an entry exist?
3562 mfIndex = cMap.findMaster(coupleMap::FACE, sfIndex);
3568 Pout<< " Removing map for master face: "
3569 << mfIndex << " :: " << faces_[mfIndex]
3571 << cMap.findSlave(coupleMap::FACE, mfIndex)
3572 << " and slave face: " << sfIndex
3573 << " on proc: " << procIndices_[pI]
3577 cMap.removeSlave(coupleMap::FACE, sfIndex);
3578 cMap.removeMaster(coupleMap::FACE, mfIndex);
3582 // Check added faces
3583 const List<objectMap>& asfList = slaveMap.addedFaceList();
3585 forAll(asfList, indexI)
3587 label sfIndex = asfList[indexI].index(), mfIndex = -1;
3589 if (mesh.whichPatch(sfIndex) == slaveConvertPatch[pI])
3591 // Configure a comparison face
3592 face cFace(mesh.faces_[sfIndex]);
3594 forAll(cFace, pointI)
3606 bool foundMatch = false;
3608 forAllConstIter(Map<label>, procFaceMap, fIter)
3610 const face& mFace = faces_[fIter()];
3612 // Discard dissimilar face sizes
3613 if (mFace.size() != cFace.size())
3618 if (cFace.size() == 3)
3620 // Optimized triangular face comparison
3625 triFace(mFace[0], mFace[1], mFace[2]),
3626 triFace(cFace[0], cFace[1], cFace[2])
3637 // Regular face compare
3638 if (face::compare(mFace, cFace))
3651 Pout<< " Found master face: " << mfIndex
3652 << " :: " << faces_[mfIndex]
3653 << " for slave face: " << sfIndex
3654 << " :: " << mesh.faces_[sfIndex]
3655 << " on proc: " << procIndices_[pI]
3659 // Update maps for the new face
3660 cMap.mapSlave(coupleMap::FACE, mfIndex, sfIndex);
3661 cMap.mapMaster(coupleMap::FACE, sfIndex, mfIndex);
3665 // Something is wrong here.
3666 Pout<< " Could not find master face for: " << nl
3667 << " Slave face: " << sfIndex
3668 << " :: " << mesh.faces_[sfIndex] << nl
3669 << " cFace: " << cFace << nl
3670 << " on proc: " << procIndices_[pI] << nl
3671 << abort(FatalError);
3676 // Map edges for 3D meshes
3679 // Map edges from patches.
3680 const Map<label>& procEdgeMap = edgesToInsert[pI];
3682 // Check removed edges
3683 const labelList& rseList = slaveMap.removedEdgeList();
3685 forAll(rseList, indexI)
3687 label seIndex = rseList[indexI], meIndex = -1;
3689 // Does an entry exist?
3690 meIndex = cMap.findMaster(coupleMap::EDGE, seIndex);
3696 Pout<< " Removing map for master edge: "
3697 << meIndex << " :: " << edges_[meIndex]
3699 << cMap.findSlave(coupleMap::EDGE, meIndex)
3700 << " and slave edge: " << seIndex
3701 << " on proc: " << procIndices_[pI]
3705 cMap.removeSlave(coupleMap::EDGE, seIndex);
3706 cMap.removeMaster(coupleMap::EDGE, meIndex);
3710 // Check added edges
3711 const List<objectMap>& aseList = slaveMap.addedEdgeList();
3713 forAll(aseList, indexI)
3715 label seIndex = aseList[indexI].index(), meIndex = -1;
3717 if (mesh.whichEdgePatch(seIndex) == slaveConvertPatch[pI])
3719 // Configure a comparison edge
3720 edge cEdge(mesh.edges_[seIndex]);
3722 cEdge[0] = cMap.findMaster(coupleMap::POINT, cEdge[0]);
3723 cEdge[1] = cMap.findMaster(coupleMap::POINT, cEdge[1]);
3725 bool foundMatch = false;
3727 forAllConstIter(Map<label>, procEdgeMap, eIter)
3729 const edge& mEdge = edges_[eIter()];
3743 Pout<< " Found master edge: " << meIndex
3744 << " :: " << edges_[meIndex]
3745 << " for slave edge: " << seIndex
3746 << " :: " << mesh.edges_[seIndex]
3747 << " on proc: " << procIndices_[pI]
3751 // Update maps for the new edge
3752 cMap.mapSlave(coupleMap::EDGE, meIndex, seIndex);
3753 cMap.mapMaster(coupleMap::EDGE, seIndex, meIndex);
3757 // Something is wrong here.
3758 Pout<< " Could not find master edge for: " << nl
3759 << " Slave edge: " << seIndex
3760 << " :: " << mesh.edges_[seIndex] << nl
3761 << " cEdge: " << cEdge << nl
3762 << " on proc: " << procIndices_[pI] << nl
3763 << abort(FatalError);
3769 // Check removed points
3770 const labelList& rspList = slaveMap.removedPointList();
3772 forAll(rspList, indexI)
3774 label spIndex = rspList[indexI], mpIndex = -1;
3776 // Does an entry exist?
3777 mpIndex = cMap.findMaster(coupleMap::POINT, spIndex);
3783 Pout<< " Removing map for master point: " << mpIndex
3784 << " :: new: " << points_[mpIndex]
3785 << " :: old: " << oldPoints_[mpIndex] << nl
3786 << " Master point map: "
3787 << cMap.findSlave(coupleMap::POINT, mpIndex)
3788 << " and slave point: " << spIndex
3789 << " on proc: " << procIndices_[pI]
3793 cMap.removeSlave(coupleMap::POINT, spIndex);
3794 cMap.removeMaster(coupleMap::POINT, mpIndex);
3799 // Write out cells for debug, if necessary
3800 if (debug > 2 || map.index() < 0)
3802 const List<objectMap>& acList = map.addedCellList();
3804 labelList addedCells(acList.size(), -1);
3806 forAll(acList, indexI)
3808 addedCells[indexI] = acList[indexI].index();
3814 "insertCells(" + Foam::name(mIndex) + ')',
3815 addedCells, 3, false, true
3819 // Set mapping information for any added faces
3820 const List<objectMap>& afList = map.addedFaceList();
3822 forAll(afList, indexI)
3824 label mfIndex = afList[indexI].index();
3825 label patch = whichPatch(mfIndex);
3826 label neiProc = getNeighbourProcessor(patch);
3828 // Disregard non-processor faces for coupled mapping
3834 // Update couple maps, if necessary
3835 const face& mFace = faces_[mfIndex];
3837 forAll(procIndices_, pI)
3839 // Skip face if not on this processor
3840 if (neiProc != procIndices_[pI])
3845 const coupleMap& cMap = recvMeshes_[pI].map();
3846 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
3848 // Configure a comparison face
3851 forAll(cFace, pointI)
3863 // Ensure all points are found
3864 if (findIndex(cFace, -1) > -1)
3869 // Loop through boundary faces on slave and look for a match
3870 label sTot = mesh.faces_.size(), sInt = mesh.nOldInternalFaces_;
3874 for (label sfI = sInt; sfI < sTot; sfI++)
3876 const face& sFace = mesh.faces_[sfI];
3883 label pIndex = mesh.whichPatch(sfI);
3885 if (mesh.getNeighbourProcessor(pIndex) == -1)
3892 if (face::compare(cFace, sFace))
3903 triFace(cFace[0], cFace[1], cFace[2]),
3904 triFace(sFace[0], sFace[1], sFace[2])
3912 // Break out if we're done
3917 Pout<< " Matched master face: " << mfIndex
3919 << " with slave: " << sfIndex
3921 << " on proc: " << procIndices_[pI]
3925 // Found the slave. Add a map entry
3946 Pout<< " Failed match for master face: " << mfIndex
3948 << " on proc: " << procIndices_[pI]
3949 << abort(FatalError);
3957 // Map edges on this face as well
3958 const labelList& mfEdges = faceEdges_[mfIndex];
3959 const labelList& sfEdges = mesh.faceEdges_[sfIndex];
3961 forAll(mfEdges, edgeI)
3963 label meIndex = mfEdges[edgeI];
3964 const edge& mEdge = edges_[meIndex];
3966 // Configure a comparison edge
3969 cMap.findSlave(coupleMap::POINT, mEdge[0]),
3970 cMap.findSlave(coupleMap::POINT, mEdge[1])
3973 forAll(sfEdges, edgeJ)
3975 label seIndex = sfEdges[edgeJ];
3976 const edge& sEdge = mesh.edges_[seIndex];
3982 Pout<< " Matched master edge: " << meIndex
3984 << " with slave: " << seIndex
3986 << " on proc: " << procIndices_[pI]
3990 // Found the slave. Add a map entry
4012 // Specify that the operation was successful
4015 // Return the changeMap
4020 // Remove the specified cells from the mesh,
4021 // and add internal faces/edges to the specified patch
4022 // - Returns a changeMap with a type specifying:
4023 // 1: Operation was successful
4024 // -1: Operation failed
4025 // - checkOnly performs a feasibility check and returns without modifications.
4026 const changeMap dynamicTopoFvMesh::removeCells
4028 const labelList& cList,
4034 if (cList.empty() || patch < 0)
4038 writeVTK("removeCellsFailed_" + rcN, cList, 3, false, true);
4043 "const changeMap dynamicTopoFvMesh::removeCells"
4044 "(const labelList&, const label, const word&, bool)"
4046 << " Wrong arguments. " << nl
4047 << " cList: " << cList << nl
4048 << " patch: " << patch << nl
4049 << abort(FatalError);
4054 labelHashSet pointsToRemove, edgesToRemove, facesToRemove;
4055 Map<label> facesToConvert, edgesToConvert;
4057 // First loop through all cells and accumulate
4058 // a set of faces to be removed/converted.
4059 forAll(cList, cellI)
4061 const cell& cellToCheck = cells_[cList[cellI]];
4063 forAll(cellToCheck, faceI)
4065 label own = owner_[cellToCheck[faceI]];
4066 label nei = neighbour_[cellToCheck[faceI]];
4070 if (!facesToRemove.found(cellToCheck[faceI]))
4072 facesToRemove.insert(cellToCheck[faceI]);
4078 (findIndex(cList, own) != -1) &&
4079 (findIndex(cList, nei) != -1)
4082 if (!facesToRemove.found(cellToCheck[faceI]))
4084 facesToRemove.insert(cellToCheck[faceI]);
4089 facesToConvert.set(cellToCheck[faceI], -1);
4094 // Add all edges as candidates for conversion.
4095 // Some of these will be removed altogether.
4096 forAllConstIter(labelHashSet, facesToRemove, fIter)
4098 const labelList& fEdges = faceEdges_[fIter.key()];
4100 forAll(fEdges, edgeI)
4102 if (whichEdgePatch(fEdges[edgeI]) == patch)
4104 // Make an identical map
4105 edgesToConvert.set(fEdges[edgeI], fEdges[edgeI]);
4109 edgesToConvert.set(fEdges[edgeI], -1);
4114 // Add all points as candidates for removal.
4115 // Some (or all) of these will be weeded out.
4116 const edge& eCheck = edges_[fEdges[edgeI]];
4118 pointsToRemove.set(eCheck[0]);
4119 pointsToRemove.set(eCheck[1]);
4124 forAllConstIter(Map<label>, facesToConvert, fIter)
4126 const labelList& fEdges = faceEdges_[fIter.key()];
4128 forAll(fEdges, edgeI)
4130 if (whichEdgePatch(fEdges[edgeI]) == patch)
4132 // Make an identical map
4133 edgesToConvert.set(fEdges[edgeI], fEdges[edgeI]);
4137 edgesToConvert.set(fEdges[edgeI], -1);
4142 // Build a list of edges to be removed.
4143 forAllConstIter(Map<label>, edgesToConvert, eIter)
4145 const labelList& eFaces = edgeFaces_[eIter.key()];
4147 bool allRemove = true;
4149 forAll(eFaces, faceI)
4151 if (!facesToRemove.found(eFaces[faceI]))
4160 if (!edgesToRemove.found(eIter.key()))
4162 edgesToRemove.insert(eIter.key());
4167 // Weed-out the conversion list.
4168 forAllConstIter(labelHashSet, edgesToRemove, eIter)
4170 edgesToConvert.erase(eIter.key());
4173 // Build a set of points to be removed.
4176 forAllConstIter(Map<label>, edgesToConvert, eIter)
4178 const edge& eCheck = edges_[eIter.key()];
4180 if (pointsToRemove.found(eCheck[0]))
4182 pointsToRemove.erase(eCheck[0]);
4185 if (pointsToRemove.found(eCheck[1]))
4187 pointsToRemove.erase(eCheck[1]);
4193 forAllConstIter(labelHashSet, edgesToRemove, eIter)
4195 const edge& edgeToCheck = edges_[eIter.key()];
4197 forAll(edgeToCheck, pointI)
4199 const labelList& pEdges = pointEdges_[edgeToCheck[pointI]];
4201 bool allRemove = true;
4203 forAll(pEdges, edgeI)
4205 if (!edgesToRemove.found(pEdges[edgeI]))
4214 if (!pointsToRemove.found(edgeToCheck[pointI]))
4216 pointsToRemove.insert(edgeToCheck[pointI]);
4223 forAllIter(Map<label>, edgesToConvert, eIter)
4225 const labelList& eFaces = edgeFaces_[eIter.key()];
4227 label nConvFaces = 0;
4229 forAll(eFaces, faceI)
4231 if (facesToConvert.found(eFaces[faceI]))
4239 Pout<< "Invalid conversion. Bailing out." << endl;
4244 // Are we only performing checks?
4247 // Make necessary map entries
4248 forAllConstIter(Map<label>, facesToConvert, fIter)
4250 map.removeFace(fIter.key());
4253 forAllConstIter(Map<label>, edgesToConvert, eIter)
4255 map.removeEdge(eIter.key());
4258 forAllConstIter(labelHashSet, facesToRemove, fIter)
4260 map.removeFace(fIter.key());
4263 forAllConstIter(labelHashSet, edgesToRemove, eIter)
4265 map.removeEdge(eIter.key());
4268 forAllConstIter(labelHashSet, pointsToRemove, pIter)
4270 map.removePoint(pIter.key());
4273 forAll(cList, cellI)
4275 map.removeCell(cList[cellI]);
4278 // Specify that the operation was successful
4284 // Perform a few debug calls before modifications
4288 << " Removing cell(s): " << cList
4289 << " and adding internal faces / edges to patch: "
4290 << boundaryMesh()[patch].name()
4294 // Write out candidates for post-processing
4297 writeVTK("pointsToRemove_" + rcN, pointsToRemove.toc(), 0, false, true);
4298 writeVTK("edgesToRemove_" + rcN, edgesToRemove.toc(), 1, false, true);
4299 writeVTK("facesToRemove_" + rcN, facesToRemove.toc(), 2, false, true);
4300 writeVTK("cellsToRemove_" + rcN, cList, 3, false, true);
4301 writeVTK("edgesToConvert_" + rcN, edgesToConvert.toc(), 1, false, true);
4302 writeVTK("facesToConvert_" + rcN, facesToConvert.toc(), 2, false, true);
4305 // Loop through all faces for conversion, check orientation
4306 // and create new faces in their place.
4307 forAllIter(Map<label>, facesToConvert, fIter)
4309 // Check if this internal face is oriented properly.
4311 label newOwner = -1;
4312 labelList fEdges = faceEdges_[fIter.key()];
4314 if (findIndex(cList, neighbour_[fIter.key()]) != -1)
4316 // Orientation is correct
4317 newFace = faces_[fIter.key()];
4318 newOwner = owner_[fIter.key()];
4321 if (findIndex(cList, owner_[fIter.key()]) != -1)
4323 // Face is to be reversed.
4324 newFace = faces_[fIter.key()].reverseFace();
4325 newOwner = neighbour_[fIter.key()];
4327 setFlip(fIter.key());
4331 // Something's terribly wrong
4335 "const changeMap dynamicTopoFvMesh::removeCells\n"
4337 " const labelList& cList,\n"
4338 " const label patch,\n"
4339 " const word& rcN,\n"
4343 << nl << " Invalid mesh. "
4344 << abort(FatalError);
4347 // Insert the reconfigured face at the boundary.
4348 // faceEdges will be corrected later.
4361 // Add this face to the map.
4362 map.addFace(fIter());
4364 // Set mapping information
4365 // - But where to map from?
4366 setFaceMapping(fIter());
4368 // Replace cell with the new face label
4369 meshOps::replaceLabel
4376 // Remove the internal face.
4377 removeFace(fIter.key());
4380 map.removeFace(fIter.key());
4383 // Create a new edge for each converted edge
4384 forAllIter(Map<label>, edgesToConvert, eIter)
4388 // Create copies before appending.
4389 edge newEdge = edges_[eIter.key()];
4390 labelList eFaces = edgeFaces_[eIter.key()];
4402 // Add this edge to the map.
4403 map.addEdge(eIter());
4406 removeEdge(eIter.key());
4409 map.removeEdge(eIter.key());
4413 // Loop through all faces for conversion, and replace edgeFaces.
4414 forAllConstIter(Map<label>, facesToConvert, fIter)
4416 // Make a copy, because this list is going to
4417 // be modified within this loop.
4418 labelList fEdges = faceEdges_[fIter()];
4420 forAll(fEdges, edgeI)
4422 if (edgesToConvert.found(fEdges[edgeI]))
4424 meshOps::replaceLabel
4428 edgeFaces_[edgesToConvert[fEdges[edgeI]]]
4431 meshOps::replaceLabel
4434 edgesToConvert[fEdges[edgeI]],
4441 // Loop through all edges for conversion, and size-down edgeFaces.
4442 forAllConstIter(Map<label>, edgesToConvert, eIter)
4444 // Make a copy, because this list is going to
4445 // be modified within this loop.
4446 labelList eFaces = edgeFaces_[eIter()];
4448 forAll(eFaces, faceI)
4450 if (facesToRemove.found(eFaces[faceI]))
4452 meshOps::sizeDownList
4459 // Replace old edges with new ones.
4460 labelList& fEdges = faceEdges_[eFaces[faceI]];
4462 forAll(fEdges, edgeI)
4464 if (edgesToConvert.found(fEdges[edgeI]))
4466 fEdges[edgeI] = edgesToConvert[fEdges[edgeI]];
4472 // Remove unwanted faces
4473 forAllConstIter(labelHashSet, facesToRemove, fIter)
4475 removeFace(fIter.key());
4478 map.removeFace(fIter.key());
4481 // Remove unwanted edges
4482 forAllConstIter(labelHashSet, edgesToRemove, eIter)
4484 removeEdge(eIter.key());
4487 map.removeEdge(eIter.key());
4490 // Remove unwanted points
4491 forAllConstIter(labelHashSet, pointsToRemove, pIter)
4493 removePoint(pIter.key());
4496 map.removePoint(pIter.key());
4500 forAll(cList, cellI)
4502 removeCell(cList[cellI]);
4505 map.removeCell(cList[cellI]);
4509 topoChangeFlag_ = true;
4511 // Specify that the operation was successful
4518 // Handle topology changes for coupled patches
4519 void dynamicTopoFvMesh::handleCoupledPatches
4521 labelHashSet& entities
4524 // Initialize coupled patch connectivity for topology modifications.
4525 initCoupledConnectivity(this);
4527 if (patchCoupling_.empty() && procIndices_.empty())
4532 // Move coupled subMeshes
4533 moveCoupledSubMeshes();
4535 // Exchange length-scale buffers across processors.
4536 exchangeLengthBuffers();
4540 // Check coupled-patch sizes first.
4541 forAll(patchCoupling_, patchI)
4543 if (patchCoupling_(patchI))
4545 const coupleMap& cMap = patchCoupling_[patchI].map();
4547 label mSize = patchSizes_[cMap.masterIndex()];
4548 label sSize = patchSizes_[cMap.slaveIndex()];
4552 Pout<< "Coupled patch-count is inconsistent." << nl
4553 << " Master Patch: " << cMap.masterIndex()
4554 << " Count: " << mSize << nl
4555 << " Slave Patch: " << cMap.slaveIndex()
4556 << " Count: " << sSize
4561 "void dynamicTopoFvMesh::handleCoupledPatches()"
4563 << " Failures were found in connectivity"
4564 << " prior to coupled topo-changes."
4565 << abort(FatalError);
4570 Info<< "Handling coupled patches...";
4573 if (Pstream::parRun())
4575 // Calculate mapping offsets
4576 const polyBoundaryMesh& boundary = boundaryMesh();
4578 // Determine number of physical (non-processor) patches
4579 label nPhysical = 0;
4581 forAll(boundary, patchI)
4583 if (isA<processorPolyPatch>(boundary[patchI]))
4591 // Fetch number of processors
4592 label nProcs = procIndices_.size();
4594 // Allocate cell, face and patch offsets
4595 labelList cellSizes(nProcs, 0), cellStarts(nProcs, 0);
4596 labelList faceSizes(nProcs, 0), faceStarts(nProcs, 0);
4597 labelListList patchSizes(nProcs, labelList(nPhysical, 0));
4598 labelListList patchStarts(nProcs, labelList(nPhysical, 0));
4600 label nTotalCells = nOldCells_, nTotalIntFaces = nOldInternalFaces_;
4601 labelList nTotalPatchFaces(SubList<label>(oldPatchSizes_, nPhysical));
4603 forAll(procIndices_, pI)
4605 const coupleMap& cMap = recvMeshes_[pI].map();
4607 // Fetch size from subMesh
4608 label nCells = cMap.nEntities(coupleMap::CELL);
4609 label nIntFaces = cMap.nEntities(coupleMap::INTERNAL_FACE);
4611 // Set size / offset for this processor
4612 cellSizes[pI] = nCells;
4613 cellStarts[pI] = nTotalCells;
4615 faceSizes[pI] = nIntFaces;
4616 faceStarts[pI] = nTotalIntFaces;
4619 nTotalCells += nCells;
4620 nTotalIntFaces += nIntFaces;
4622 // Fetch patch sizes from subMesh
4623 const labelList& nPatchFaces =
4625 cMap.entityBuffer(coupleMap::FACE_SIZES)
4628 // Loop over physical patches
4629 forAll(nTotalPatchFaces, patchI)
4631 // Fetch patch size from subMesh
4632 label nFaces = nPatchFaces[patchI];
4634 // Set patch size / offset for this processor
4635 patchSizes[pI][patchI] = nFaces;
4636 patchStarts[pI][patchI] = nTotalPatchFaces[patchI];
4638 // Update patch count
4639 nTotalPatchFaces[patchI] += nFaces;
4643 // Set sizes / starts in mapper
4656 SubList<label> physicalPatches(oldPatchSizes_, nPhysical);
4658 Pout<< " procIndices: " << procIndices_ << nl
4659 << " nCells: " << nOldCells_ << nl
4660 << " proc cellSizes: " << cellSizes << nl
4661 << " cellStarts: " << cellStarts << nl
4662 << " proc faceSizes: " << faceSizes << nl
4663 << " faceStarts: " << faceStarts << nl
4664 << " patchSizes: " << physicalPatches << nl
4665 << " proc patchSizes: " << patchSizes << nl
4666 << " patchStarts: " << patchStarts << endl;
4670 // Optionally switch off topo-changes
4671 // on processor patches at run-time
4672 bool procTopoChanges = true;
4674 const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
4676 if (meshSubDict.found("procTopoChanges") || mandatory_)
4678 procTopoChanges = readBool(meshSubDict.lookup("procTopoChanges"));
4681 // Set coupled modifications.
4682 setCoupledModification();
4684 // Loop through the coupled stack and perform changes.
4685 if (edgeRefinement_)
4687 // Initialize the coupled stack
4688 if (procTopoChanges)
4690 initCoupledStack(entities, false);
4693 edgeRefinementEngine(&(handlerPtr_[0]));
4696 // Re-Initialize the stack
4697 if (procTopoChanges)
4699 initCoupledStack(entities, false);
4704 swap2DEdges(&(handlerPtr_[0]));
4708 swap3DEdges(&(handlerPtr_[0]));
4711 // Build a list of entities that need to be avoided
4712 // by regular topo-changes.
4713 buildEntitiesToAvoid(entities, true);
4715 // Reset coupled modifications.
4716 unsetCoupledModification();
4720 Info<< "Done." << endl;
4722 // Check coupled-patch sizes after changes.
4723 forAll(patchCoupling_, patchI)
4725 if (patchCoupling_(patchI))
4727 const coupleMap& cMap = patchCoupling_[patchI].map();
4729 label mSize = patchSizes_[cMap.masterIndex()];
4730 label sSize = patchSizes_[cMap.slaveIndex()];
4734 Pout<< "Coupled patch-count is inconsistent." << nl
4735 << " Master Patch: " << cMap.masterIndex()
4736 << " Count: " << mSize << nl
4737 << " Slave Patch: " << cMap.slaveIndex()
4738 << " Count: " << sSize
4741 FatalErrorIn("dynamicTopoFvMesh::handleCoupledPatches()")
4742 << " Failures were found in connectivity"
4743 << " after coupled topo-changes."
4744 << abort(FatalError);
4750 // Schedule transfer of topology operations across processors
4751 forAll(procIndices_, pI)
4753 label proc = procIndices_[pI];
4755 coupledMesh& sPM = sendMeshes_[pI];
4756 coupledMesh& rPM = recvMeshes_[pI];
4758 if (priority(proc, lessOp<label>(), Pstream::myProcNo()))
4760 const coupleMap& cMap = sPM.map();
4762 // How many entities am I receiving..
4763 FixedList<label, 3> nOps(-1);
4765 meshOps::pRead(proc, nOps);
4767 label nEntities = nOps[0];
4768 label nMovePoints = nOps[1];
4769 label nPatchConvert = nOps[2];
4773 Pout<< " Op transfer:"
4774 << " Receiving from [" << proc << "]::"
4775 << " nEntities: " << nEntities
4776 << " nMovePoints: " << nMovePoints
4777 << " nPatchConvert: " << nPatchConvert
4783 // Size up the receipt buffers
4784 labelList& indices = cMap.entityIndices();
4785 List<coupleMap::opType>& operations = cMap.entityOperations();
4787 indices.setSize(nEntities);
4788 operations.setSize(nEntities);
4790 // Schedule indices and operations for receipt
4791 meshOps::pRead(proc, indices);
4792 meshOps::pRead(proc, operations);
4797 // Size up the receipt buffers
4798 pointField& newPoints = cMap.moveNewPoints();
4799 pointField& oldPoints = cMap.moveOldPoints();
4801 newPoints.setSize(nMovePoints);
4802 oldPoints.setSize(nMovePoints);
4804 // Schedule old / new points for receipt
4805 meshOps::pRead(proc, newPoints);
4806 meshOps::pRead(proc, oldPoints);
4811 // Size up the receipt buffers
4812 labelList& patches = cMap.patchIndices();
4814 patches.setSize(nPatchConvert);
4816 // Schedule patch indices for receipt
4817 meshOps::pRead(proc, patches);
4822 const coupleMap& cMap = rPM.map();
4824 FixedList<label, 3> nOps(-1);
4826 label nEntities = cMap.entityIndices().size();
4827 label nMovePoints = cMap.moveNewPoints().size();
4828 label nPatchConvert = cMap.patchIndices().size();
4832 Pout<< " Op transfer:"
4833 << " Sending to [" << proc << "]:: "
4834 << " nEntities: " << nEntities << nl
4835 << " entityIndices: " << cMap.entityIndices() << nl
4836 << " entityOperations: " << cMap.entityOperations() << nl
4837 << " nMovePoints: " << nMovePoints << nl
4838 << " moveNewPoints: " << cMap.moveNewPoints() << nl
4839 << " moveOldPoints: " << cMap.moveOldPoints() << nl
4840 << " nPatchConvert: " << nPatchConvert << nl
4841 << " patchIndices: " << cMap.patchIndices() << nl
4845 nOps[0] = nEntities;
4846 nOps[1] = nMovePoints;
4847 nOps[2] = nPatchConvert;
4849 meshOps::pWrite(proc, nOps);
4853 // Schedule transfer to processor
4854 meshOps::pWrite(proc, cMap.entityIndices());
4855 meshOps::pWrite(proc, cMap.entityOperations());
4860 // Schedule transfer to processor
4861 meshOps::pWrite(proc, cMap.moveNewPoints());
4862 meshOps::pWrite(proc, cMap.moveOldPoints());
4867 // Schedule transfer to processor
4868 meshOps::pWrite(proc, cMap.patchIndices());
4873 // We won't wait for transfers to complete for the moment,
4874 // and will deal with operations once the internal mesh
4875 // has been dealt with.
4879 // Synchronize topology operations across processors
4880 void dynamicTopoFvMesh::syncCoupledPatches(labelHashSet& entities)
4882 if (!Pstream::parRun())
4887 // Temporarily reset maxModifications to
4888 // ensure that synchronization succeeds
4889 label maxModSave = maxModifications_;
4891 maxModifications_ = -1;
4893 const polyBoundaryMesh& boundary = boundaryMesh();
4895 labelList createPatchOrders(procIndices_.size(), 0);
4897 forAll(procIndices_, pI)
4899 const coupleMap& cMap = recvMeshes_[pI].map();
4901 if (cMap.patchIndex() >= boundary.size())
4903 // This processor needs a new patch talking to me
4904 createPatchOrders[pI] = 1;
4908 // Send / recv create patch orders, if any
4909 forAll(procIndices_, pI)
4911 label proc = procIndices_[pI];
4913 if (priority(proc, greaterOp<label>(), Pstream::myProcNo()))
4915 // Non-blocking send to processor
4916 meshOps::pWrite(proc, createPatchOrders[pI], false);
4920 // Non-blocking receive from processor
4921 meshOps::pRead(proc, createPatchOrders[pI], false);
4925 // Wait for all transfers to complete.
4926 meshOps::waitForBuffers();
4928 Map<label> addedProcPatches;
4930 forAll(procIndices_, pI)
4932 label proc = procIndices_[pI];
4934 if (createPatchOrders[pI])
4936 // Create a new processor patch
4937 addedProcPatches.insert
4940 createProcessorPatch(proc)
4945 // Buffer for cell-removal
4946 dynamicLabelList rCellList(10);
4948 forAll(procIndices_, pI)
4950 label proc = procIndices_[pI];
4952 const coupledMesh& sPM = sendMeshes_[pI];
4954 if (priority(proc, lessOp<label>(), Pstream::myProcNo()))
4956 const coupleMap& cMap = sPM.map();
4957 const labelList& indices = cMap.entityIndices();
4958 const List<coupleMap::opType>& operations = cMap.entityOperations();
4960 const labelList& patches = cMap.patchIndices();
4961 const pointField& newPoints = cMap.moveNewPoints();
4962 const pointField& oldPoints = cMap.moveOldPoints();
4964 // Find the appropriate processor patch
4965 // for cell-removal operations
4966 label procPatch = -1;
4968 forAll(boundary, pI)
4970 if (!isA<processorPolyPatch>(boundary[pI]))
4975 const processorPolyPatch& pp =
4977 refCast<const processorPolyPatch>(boundary[pI])
4980 if (pp.neighbProcNo() == proc)
4987 if (procPatch == -1)
4989 // Check if this was a newly added patch
4990 if (addedProcPatches.found(proc))
4992 procPatch = addedProcPatches[proc];
4995 if (operations.size())
4997 bool required = false;
4999 forAll(operations, opI)
5001 const coupleMap::opType op = operations[opI];
5005 op == coupleMap::REMOVE_CELL ||
5006 op == coupleMap::CONVERT_PATCH
5014 // If we actually have operations from
5015 // this processor that require patch conversion,
5016 // this is a problem.
5019 Pout<< " * * * Sync Operations * * * " << nl
5020 << " Could not find patch for proc: " << proc << nl
5021 << " procIndices: " << procIndices_
5022 << abort(FatalError);
5027 // Move on to the next processor
5032 label pointCounter = 0, patchCounter = 0;
5034 // Sequentially execute operations
5035 forAll(indices, indexI)
5037 const label index = indices[indexI];
5038 const coupleMap::opType op = operations[indexI];
5042 Pout<< " Recv Op index: " << index << endl;
5045 // Determine the appropriate local index
5046 label localIndex = -1;
5048 if (op == coupleMap::MOVE_POINT)
5050 localIndex = cMap.entityMap(coupleMap::POINT)[index];
5055 op == coupleMap::CONVERT_PATCH ||
5056 op == coupleMap::CONVERT_PHYSICAL
5061 cMap.entityMap(coupleMap::FACE).found(index) ?
5062 cMap.entityMap(coupleMap::FACE)[index] : -1
5067 if (op == coupleMap::CONVERT_PATCH)
5069 const point& newCentre = newPoints[pointCounter];
5070 const point& oldCentre = oldPoints[pointCounter];
5072 Pout<< " Convert patch: " << index << nl
5073 << " localIndex: " << localIndex << nl
5074 << " newCentre: " << newCentre << nl
5075 << " oldCentre: " << oldCentre << nl
5079 if (op == coupleMap::CONVERT_PHYSICAL)
5081 Pout<< " Convert patch: " << index << nl
5082 << " localIndex: " << localIndex << nl
5083 << " patch: " << patches[patchCounter] << nl
5090 // Pick localIndex based on entity type
5093 const Map<label>& fM = cMap.entityMap(coupleMap::FACE);
5095 Map<label>::const_iterator it = fM.find(index);
5103 Pout<< " * * * Sync Operations * * * " << nl
5104 << " Could not find index: " << index
5105 << " in faceMap for proc: " << proc << nl
5107 << cMap.nEntities(coupleMap::FACE)
5108 << " opType: " << coupleMap::asText(op)
5109 << abort(FatalError);
5114 const Map<label>& eM = cMap.entityMap(coupleMap::EDGE);
5116 Map<label>::const_iterator it = eM.find(index);
5124 Pout<< " * * * Sync Operations * * * " << nl
5125 << " Could not find index: " << index
5126 << " in edgeMap for proc: " << proc << nl
5128 << cMap.nEntities(coupleMap::EDGE)
5129 << " opType: " << coupleMap::asText(op)
5130 << abort(FatalError);
5139 case coupleMap::BISECTION:
5141 opMap = bisectEdge(localIndex);
5145 case coupleMap::COLLAPSE_FIRST:
5147 opMap = collapseEdge(localIndex, 1);
5151 case coupleMap::COLLAPSE_SECOND:
5153 opMap = collapseEdge(localIndex, 2);
5157 case coupleMap::COLLAPSE_MIDPOINT:
5159 opMap = collapseEdge(localIndex, 3);
5163 case coupleMap::REMOVE_CELL:
5165 // Clear existing list
5170 // Insert the owner cell
5171 rCellList.append(owner_[localIndex]);
5175 // Insert all cells connected to this edge
5176 const labelList& eFaces = edgeFaces_[localIndex];
5178 forAll(eFaces, faceI)
5180 const label own = owner_[eFaces[faceI]];
5181 const label nei = neighbour_[eFaces[faceI]];
5183 if (findIndex(rCellList, own) == -1)
5185 rCellList.append(own);
5190 if (findIndex(rCellList, nei) == -1)
5192 rCellList.append(nei);
5204 "rcs_" + Foam::name(localIndex)
5211 case coupleMap::MOVE_POINT:
5213 const point& newPoint = newPoints[pointCounter];
5214 const point& oldPoint = oldPoints[pointCounter];
5216 // Move old / new points
5217 points_[localIndex] = newPoint;
5218 oldPoints_[localIndex] = oldPoint;
5220 // Clear the existing map
5223 // Force a successful operation
5231 case coupleMap::CONVERT_PATCH:
5233 const point& newCentre = newPoints[pointCounter];
5234 const point& oldCentre = oldPoints[pointCounter];
5236 if (localIndex == -1)
5238 // Accumulate stats in case of failure
5239 scalar minDist = GREAT;
5240 vector minPoint = vector::zero;
5241 dynamicLabelList checkedFaces;
5245 // Reserve for append
5246 checkedFaces.setCapacity(50);
5249 // New face. Check all boundary faces
5250 // and match up centre
5251 label sTot = faces_.size();
5252 label sInt = nOldInternalFaces_;
5254 for (label faceI = sInt; faceI < sTot; faceI++)
5256 const face& fCheck = faces_[faceI];
5263 label pIndex = whichPatch(faceI);
5265 if (getNeighbourProcessor(pIndex) == -1)
5270 // Compute face-centre
5271 vector fC = fCheck.centre(points_);
5273 // Compute tolerance
5274 scalar tol = mag(points_[fCheck[0]] - fC);
5275 scalar dist = mag(fC - newCentre);
5277 if (dist < (geomMatchTol_()*tol))
5290 checkedFaces.append(faceI);
5295 // Ensure that the face was found
5296 if (localIndex == -1)
5301 + Foam::name(index),
5306 Pout<< " * * * syncCoupledPatches * * * " << nl
5307 << " Convert patch Op failed." << nl
5308 << " Face: " << index << nl
5309 << " minPoint: " << minPoint << nl
5310 << " minDistance: " << minDist << nl
5311 << " newCentre: " << newCentre << nl
5312 << " oldCentre: " << oldCentre << nl
5313 << abort(FatalError);
5317 // Fetch reference to face
5318 const face& fCheck = faces_[localIndex];
5320 point fC = fCheck.centre(points_);
5322 // Specify a tolerance
5323 scalar tol = mag(points_[fCheck[0]] - fC);
5324 scalar dist = mag(fC - newCentre);
5326 // Ensure a face-match
5327 if (dist > (geomMatchTol_() * tol))
5329 Pout<< " * * * Sync Operations * * * " << nl
5330 << " Convert patch Op failed." << nl
5331 << " Index: " << index << nl
5332 << " localIndex: " << localIndex << nl
5333 << " face: " << fCheck << nl
5334 << " faceCentre: " << fC << nl
5335 << " Master processor: " << proc << nl
5336 << " procPatch: " << procPatch << nl
5337 << " tolerance: " << (geomMatchTol_() * tol) << nl
5338 << " distance: " << dist << nl
5339 << " pointCounter: " << pointCounter << nl
5340 << " newCentre: " << newCentre << nl
5341 << " oldCentre: " << oldCentre << nl
5342 << abort(FatalError);
5345 // Obtain a copy before adding the new face,
5346 // since the reference might become
5347 // invalid during list resizing.
5348 // Edges don't have to change, since they're
5349 // all on the boundary anyway.
5350 face newFace = faces_[localIndex];
5351 label newOwn = owner_[localIndex];
5352 labelList newFaceEdges = faceEdges_[localIndex];
5354 label newFaceIndex =
5366 meshOps::replaceLabel
5373 // Correct edgeFaces with the new face label.
5374 forAll(newFaceEdges, edgeI)
5376 meshOps::replaceLabel
5380 edgeFaces_[newFaceEdges[edgeI]]
5384 // Finally remove the face
5385 removeFace(localIndex);
5387 // Clear the existing map
5394 labelList(1, localIndex)
5397 // Force a successful operation
5405 case coupleMap::CONVERT_PHYSICAL:
5407 const label patchIndex = patches[patchCounter];
5409 // Obtain a copy before adding the new face,
5410 // since the reference might become
5411 // invalid during list resizing.
5412 // Edges don't have to change, since they're
5413 // all on the boundary anyway.
5414 face newFace = faces_[localIndex];
5415 label newOwn = owner_[localIndex];
5416 labelList newFaceEdges = faceEdges_[localIndex];
5418 label newFaceIndex =
5430 meshOps::replaceLabel
5437 // Correct edgeFaces with the new face label.
5438 forAll(newFaceEdges, edgeI)
5440 meshOps::replaceLabel
5444 edgeFaces_[newFaceEdges[edgeI]]
5448 // Finally remove the face
5449 removeFace(localIndex);
5451 // Clear the existing map
5458 labelList(1, localIndex)
5461 // Force a successful operation
5469 case coupleMap::INVALID:
5471 Pout<< " * * * Sync Operations * * * " << nl
5472 << " Invalid operation." << nl
5473 << " Index: " << index << nl
5474 << " localIndex: " << localIndex << nl
5475 << " operation: " << coupleMap::asText(op) << nl
5476 << abort(FatalError);
5482 if (opMap.type() <= 0)
5484 Pout<< " * * * Sync Operations * * * " << nl
5485 << " Operation failed." << nl
5486 << " Index: " << index << nl
5487 << " localIndex: " << localIndex << nl
5488 << " operation: " << coupleMap::asText(op) << nl
5489 << " opMap.type: " << opMap.type() << nl
5490 << abort(FatalError);
5496 // Reinstate max modifications
5497 maxModifications_ = maxModSave;
5499 // Re-Initialize the stack with avoided entities from subMeshes
5500 // and leave those on processor patches untouched
5501 labelHashSet procEntities;
5503 buildEntitiesToAvoid(procEntities, false);
5505 // First remove processor entries
5506 forAllConstIter(labelHashSet, procEntities, pIter)
5508 if (entities.found(pIter.key()))
5510 entities.erase(pIter.key());
5514 // Initialize the coupled stack, using supplied entities
5515 initCoupledStack(entities, true);
5517 // Loop through the coupled stack and perform changes.
5518 if (edgeRefinement_)
5520 edgeRefinementEngine(&(handlerPtr_[0]));
5523 // Re-Initialize the stack, using supplied entities
5524 initCoupledStack(entities, true);
5528 swap2DEdges(&(handlerPtr_[0]));
5532 swap3DEdges(&(handlerPtr_[0]));
5537 // Check the state of coupled boundaries
5538 // - Return true if errors are present
5539 bool dynamicTopoFvMesh::checkCoupledBoundaries(bool report) const
5541 const polyBoundaryMesh& boundary = boundaryMesh();
5543 bool sizeError = false, misMatchError = false;
5545 // Maintain a list of master / neighbour anchors
5546 List<vectorField> mAnchors(boundary.size());
5547 List<vectorField> nAnchors(boundary.size());
5549 forAll(boundary, pI)
5551 if (isA<processorPolyPatch>(boundary[pI]))
5553 const processorPolyPatch& pp =
5555 refCast<const processorPolyPatch>(boundary[pI])
5558 label neiProcNo = pp.neighbProcNo();
5560 // Send point / face sizes
5561 meshOps::pWrite(neiProcNo, boundary[pI].nPoints());
5562 meshOps::pWrite(neiProcNo, boundary[pI].size());
5564 // Send points / centres / areas
5565 meshOps::pWrite(neiProcNo, boundary[pI].faceAreas());
5566 meshOps::pWrite(neiProcNo, boundary[pI].faceCentres());
5567 meshOps::pWrite(neiProcNo, boundary[pI].localPoints());
5569 // Prepare and send anchor points
5570 mAnchors[pI].setSize(boundary[pI].size());
5572 const faceList& lFaces = boundary[pI].localFaces();
5573 const pointField& lPoints = boundary[pI].localPoints();
5575 forAll(lFaces, faceI)
5577 mAnchors[pI][faceI] = lPoints[lFaces[faceI][0]];
5580 meshOps::pWrite(neiProcNo, mAnchors[pI]);
5583 if (isA<cyclicPolyPatch>(boundary[pI]))
5585 const cyclicPolyPatch& cp =
5587 refCast<const cyclicPolyPatch>(boundary[pI])
5594 Pout<< " Incorrect size for cyclic patch: "
5595 << cp.size() << endl;
5601 label halfSize = cp.size() / 2;
5603 vectorField half0Areas
5612 vectorField half0Centres
5621 vectorField half1Areas
5631 vectorField half1Centres
5641 const faceList& lF = boundary[pI].localFaces();
5642 const pointField& lP = boundary[pI].localPoints();
5643 const vectorField& lC = boundary[pI].faceCentres();
5645 // Prepare anchor points
5646 mAnchors[pI].setSize(boundary[pI].size());
5650 mAnchors[pI][faceI] = lP[lF[faceI][0]];
5653 // Transform first-half points and check
5654 if (cp.transform() == cyclicPolyPatch::ROTATIONAL)
5656 forAll(half0Centres, faceI)
5658 half0Centres[faceI] =
5667 mAnchors[pI][faceI] =
5678 if (cp.transform() == cyclicPolyPatch::TRANSLATIONAL)
5680 forAll(half0Centres, faceI)
5682 half0Centres[faceI] += cp.separationVector();
5683 mAnchors[pI][faceI] += cp.separationVector();
5688 Pout<< " Cyclic check: Unknown transform."
5689 << abort(FatalError);
5692 // Calculate a point-match tolerance per patch
5693 scalar pTol = -GREAT;
5695 // Check areas / compute tolerance
5696 forAll(half0Areas, faceI)
5698 scalar fMagSf = mag(half0Areas[faceI]);
5699 scalar rMagSf = mag(half1Areas[faceI]);
5700 scalar avSf = 0.5 * (fMagSf + rMagSf);
5702 if (mag(fMagSf - rMagSf)/avSf > geomMatchTol_())
5704 misMatchError = true;
5706 Pout<< " Face: " << faceI
5707 << " area does not match neighbour by: "
5708 << 100 * mag(fMagSf - rMagSf)/avSf
5709 << "% - possible patch ordering problem. "
5710 << " Front area:" << fMagSf
5711 << " Rear area: " << rMagSf
5720 geomMatchTol_() * mag(lP[lF[faceI][0]] - lC[faceI])
5725 // Check centres / anchor points
5726 forAll(half0Centres, faceI)
5733 - mAnchors[pI][faceI + halfSize]
5742 - half1Centres[faceI]
5746 if (distA > pTol || distC > pTol)
5748 misMatchError = true;
5750 UIndirectList<point> f1(lP, lF[faceI]);
5751 UIndirectList<point> f2(lP, lF[faceI + halfSize]);
5753 Pout<< " Face: " << faceI << nl
5754 << " Points: " << nl << f1 << nl << f2 << nl
5755 << " Anchors ::" << nl
5756 << mAnchors[pI][faceI] << nl
5757 << mAnchors[pI][faceI + halfSize] << nl
5758 << " Centres ::" << nl
5759 << half0Centres[faceI] << nl
5760 << half1Centres[faceI] << nl
5761 << " Tolerance: " << pTol << nl
5762 << " Measured Anchor distance: " << distA << nl
5763 << " Measured Centre distance: " << distC << nl
5770 // Write out to disk
5775 identity(cp.size()),
5786 // Maintain a list of points / areas / centres
5787 List<vectorField> fAreas(boundary.size());
5788 List<vectorField> fCentres(boundary.size());
5789 List<vectorField> neiPoints(boundary.size());
5791 forAll(boundary, pI)
5793 if (isA<processorPolyPatch>(boundary[pI]))
5795 const processorPolyPatch& pp =
5797 refCast<const processorPolyPatch>(boundary[pI])
5800 label neiProcNo = pp.neighbProcNo();
5802 // Receive point / face sizes
5803 label neiPSize = -1, neiSize = -1;
5805 meshOps::pRead(neiProcNo, neiPSize);
5806 meshOps::pRead(neiProcNo, neiSize);
5809 fAreas[pI].setSize(neiSize);
5810 fCentres[pI].setSize(neiSize);
5811 neiPoints[pI].setSize(neiPSize);
5812 nAnchors[pI].setSize(neiSize);
5814 // Receive points / centres / areas / anchors
5815 meshOps::pRead(neiProcNo, fAreas[pI]);
5816 meshOps::pRead(neiProcNo, fCentres[pI]);
5817 meshOps::pRead(neiProcNo, neiPoints[pI]);
5818 meshOps::pRead(neiProcNo, nAnchors[pI]);
5822 neiSize != boundary[pI].size() ||
5823 neiPSize != boundary[pI].nPoints()
5828 Pout<< "Incorrect send / recv sizes: " << nl
5829 << " Patch: " << boundary[pI].name() << nl
5830 << " My nPoints: " << boundary[pI].nPoints() << nl
5831 << " Nei nPoints: " << neiPSize << nl
5832 << " My nFaces: " << boundary[pI].size() << nl
5833 << " Nei nFaces: " << neiSize << nl
5839 // Wait for transfers to complete
5840 meshOps::waitForBuffers();
5842 // Check centres / areas
5843 forAll(boundary, pI)
5845 if (!isA<processorPolyPatch>(boundary[pI]))
5850 const vectorField& myAreas = boundary[pI].faceAreas();
5851 const vectorField& myCentres = boundary[pI].faceCentres();
5853 // Fetch local connectivity
5854 const faceList& myFaces = boundary[pI].localFaces();
5855 const vectorField& myPoints = boundary[pI].localPoints();
5857 // Calculate a point-match tolerance per patch
5858 scalar pTol = -GREAT;
5860 forAll(myAreas, faceI)
5862 scalar magSf = mag(myAreas[faceI]);
5863 scalar nbrMagSf = mag(fAreas[pI][faceI]);
5864 scalar avSf = 0.5 * (magSf + nbrMagSf);
5866 if (mag(magSf - nbrMagSf)/avSf > geomMatchTol_())
5868 misMatchError = true;
5870 Pout<< " Face: " << faceI
5871 << " area does not match neighbour by: "
5872 << 100 * mag(magSf - nbrMagSf)/avSf
5873 << "% - possible patch ordering problem. "
5874 << " My area:" << magSf
5875 << " Neighbour area: " << nbrMagSf
5876 << " My centre: " << myCentres[faceI]
5888 myPoints[myFaces[faceI][0]]
5895 const processorPolyPatch& pp =
5897 refCast<const processorPolyPatch>(boundary[pI])
5900 // Check anchor points
5901 forAll(mAnchors[pI], faceI)
5903 scalar dist = mag(mAnchors[pI][faceI] - nAnchors[pI][faceI]);
5907 misMatchError = true;
5909 Pout<< " Face: " << faceI << "::" << mAnchors[pI][faceI] << nl
5910 << " Mismatch for processor: " << pp.neighbProcNo() << nl
5911 << " Neighbour Anchor: " << nAnchors[pI][faceI] << nl
5912 << " Tolerance: " << pTol << nl
5913 << " Measured distance: " << dist << nl
5918 // Check neighbPoints addressing
5919 const labelList& nPts = pp.neighbPoints();
5921 forAll(myPoints, pointI)
5923 if (nPts[pointI] < 0)
5925 Pout<< " Point: " << pointI << "::" << myPoints[pointI]
5926 << " reported to be multiply connected." << nl
5927 << " neighbPoint: " << nPts[pointI] << nl
5930 misMatchError = true;
5935 scalar dist = mag(myPoints[pointI] - neiPoints[pI][nPts[pointI]]);
5939 misMatchError = true;
5941 Pout<< " Point: " << pointI << "::" << myPoints[pointI] << nl
5942 << " Mismatch for processor: " << pp.neighbProcNo() << nl
5943 << " Neighbour Pt: " << neiPoints[pI][nPts[pointI]] << nl
5944 << " Tolerance: " << pTol << nl
5945 << " Measured distance: " << dist << nl
5946 << " Neighbour Index:" << nPts[pointI] << nl
5952 if (!sizeError && !misMatchError && report)
5954 Pout<< " Coupled boundary check: No problems." << endl;
5957 return (sizeError || misMatchError);
5961 // Build patch sub-meshes for processor patches
5962 void dynamicTopoFvMesh::buildProcessorPatchMeshes()
5964 if (procIndices_.empty())
5969 // Maintain a list of cells common to multiple processors.
5970 Map<labelList> commonCells;
5972 forAll(procIndices_, pI)
5974 label proc = procIndices_[pI];
5976 coupledMesh& sPM = sendMeshes_[pI];
5978 // Build the subMesh.
5979 buildProcessorPatchMesh(sPM, commonCells);
5981 const coupleMap& scMap = sPM.map();
5983 // Send my sub-mesh to the neighbour.
5984 meshOps::pWrite(proc, scMap.nEntities());
5988 Pout<< "Sending to [" << proc << "]:: nEntities: "
5989 << scMap.nEntities() << endl;
5992 // Send the pointBuffers
5993 meshOps::pWrite(proc, scMap.pointBuffer());
5994 meshOps::pWrite(proc, scMap.oldPointBuffer());
5996 // Send connectivity (points, edges, faces, cells, etc)
5997 forAll(scMap.entityBuffer(), bufferI)
5999 if (scMap.entityBuffer(bufferI).size())
6001 meshOps::pWrite(proc, scMap.entityBuffer(bufferI));
6005 // Obtain references
6006 const coupledMesh& rPM = recvMeshes_[pI];
6007 const coupleMap& rcMap = rPM.map();
6009 // First read entity sizes.
6010 meshOps::pRead(proc, rcMap.nEntities());
6014 Pout<< "Receiving from [" << proc << "]:: nEntities: "
6015 << rcMap.nEntities() << endl;
6018 // Size the buffers.
6019 rcMap.allocateBuffers();
6021 // Receive the pointBuffers
6022 meshOps::pRead(proc, rcMap.pointBuffer());
6023 meshOps::pRead(proc, rcMap.oldPointBuffer());
6025 // Receive connectivity (points, edges, faces, cells, etc)
6026 forAll(rcMap.entityBuffer(), bufferI)
6028 if (rcMap.entityBuffer(bufferI).size())
6030 meshOps::pRead(proc, rcMap.entityBuffer(bufferI));
6035 // We won't wait for all transfers to complete for the moment.
6036 // Meanwhile, do some other useful work, if possible.
6040 // Build patch sub-mesh for a specified processor
6041 // - At this point, procIndices is available as a sorted list
6042 // of neighbouring processors.
6043 void dynamicTopoFvMesh::buildProcessorPatchMesh
6045 coupledMesh& subMesh,
6046 Map<labelList>& commonCells
6049 label nP = 0, nE = 0, nF = 0, nC = 0;
6051 // Obtain references
6052 const coupleMap& cMap = subMesh.map();
6053 const labelList& subMeshPoints = cMap.subMeshPoints();
6054 const List<labelPair>& globalProcPoints = cMap.globalProcPoints();
6056 Map<label>& rPointMap = cMap.reverseEntityMap(coupleMap::POINT);
6057 Map<label>& rEdgeMap = cMap.reverseEntityMap(coupleMap::EDGE);
6058 Map<label>& rFaceMap = cMap.reverseEntityMap(coupleMap::FACE);
6059 Map<label>& rCellMap = cMap.reverseEntityMap(coupleMap::CELL);
6061 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
6062 Map<label>& edgeMap = cMap.entityMap(coupleMap::EDGE);
6063 Map<label>& faceMap = cMap.entityMap(coupleMap::FACE);
6064 Map<label>& cellMap = cMap.entityMap(coupleMap::CELL);
6066 // Add all cells connected to points on the subMeshPoints list
6069 if (cMap.masterIndex() == Pstream::myProcNo())
6071 proc = cMap.slaveIndex();
6075 proc = cMap.masterIndex();
6078 // Check if this is a direct neighbour
6079 const polyBoundaryMesh& boundary = boundaryMesh();
6081 // Add sub-mesh points first.
6082 // Additional halo points will be added later.
6083 forAll(subMeshPoints, pointI)
6085 pointMap.insert(nP, subMeshPoints[pointI]);
6086 rPointMap.insert(subMeshPoints[pointI], nP);
6090 // Set the number of points (shared) at this point.
6091 cMap.nEntities(coupleMap::SHARED_POINT) = nP;
6093 // Size up the point-buffer with global index information.
6094 // Direct neighbours do not require any addressing.
6095 labelList& gpBuffer = cMap.entityBuffer(coupleMap::POINT);
6096 gpBuffer.setSize(globalProcPoints.size(), -1);
6098 forAll(globalProcPoints, pointI)
6100 pointMap.insert(nP, globalProcPoints[pointI].first());
6101 rPointMap.insert(globalProcPoints[pointI].first(), nP);
6104 // Fill in buffer with global point index
6105 gpBuffer[pointI] = globalProcPoints[pointI].second();
6108 // Set the number of points (shared + global) at this point.
6109 cMap.nEntities(coupleMap::GLOBAL_POINT) = nP;
6111 labelHashSet localCommonCells;
6113 // Detect all cells surrounding shared points.
6116 // No pointEdges structure, so loop through all cells
6117 // to check for those connected to subMeshPoints
6118 forAll(cells_, cellI)
6120 const cell& cellToCheck = cells_[cellI];
6122 if (cellToCheck.empty())
6127 const labelList cellPoints = cellToCheck.labels(faces_);
6129 forAll(cellPoints, pointI)
6131 if (rPointMap.found(cellPoints[pointI]))
6133 if (commonCells.found(cellI))
6135 // Add locally common cells at the end.
6136 if (!localCommonCells.found(cellI))
6138 localCommonCells.insert(cellI);
6141 // Check if the processor exists on the list
6142 // and if not, add it.
6143 labelList& procList = commonCells[cellI];
6145 if (findIndex(procList, proc) == -1)
6147 meshOps::sizeUpList(proc, procList);
6152 cellMap.insert(nC, cellI);
6153 rCellMap.insert(cellI, nC);
6156 commonCells.insert(cellI, labelList(1, proc));
6166 forAllConstIter(Map<label>, rPointMap, pIter)
6168 // Loop through pointEdges for this point.
6169 const labelList& pEdges = pointEdges_[pIter.key()];
6171 forAll(pEdges, edgeI)
6173 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
6175 forAll(eFaces, faceI)
6177 label own = owner_[eFaces[faceI]];
6178 label nei = neighbour_[eFaces[faceI]];
6181 if (!rCellMap.found(own))
6183 if (commonCells.found(own))
6185 // Add locally common cells at the end.
6186 if (!localCommonCells.found(own))
6188 localCommonCells.insert(own);
6191 // Check if the processor exists on the list
6192 // and if not, add it.
6193 labelList& procList = commonCells[own];
6195 if (findIndex(procList, proc) == -1)
6197 meshOps::sizeUpList(proc, procList);
6202 cellMap.insert(nC, own);
6203 rCellMap.insert(own, nC);
6206 commonCells.insert(own, labelList(1, proc));
6210 // Check neighbour cell
6211 if (!rCellMap.found(nei) && nei != -1)
6213 if (commonCells.found(nei))
6215 // Add locally common cells at the end.
6216 if (!localCommonCells.found(nei))
6218 localCommonCells.insert(nei);
6221 // Check if the processor exists on the list
6222 // and if not, add it.
6223 labelList& procList = commonCells[nei];
6225 if (findIndex(procList, proc) == -1)
6227 meshOps::sizeUpList(proc, procList);
6232 cellMap.insert(nC, nei);
6233 rCellMap.insert(nei, nC);
6236 commonCells.insert(nei, labelList(1, proc));
6244 // Now add locally common cells.
6245 forAllConstIter(labelHashSet, localCommonCells, cIter)
6247 cellMap.insert(nC, cIter.key());
6248 rCellMap.insert(cIter.key(), nC);
6252 // Keep track of inserted boundary face indices
6253 // - Add exposed internal faces to a 'default' patch
6254 // at the end of the list.
6255 labelList bdyFaceSizes(boundary.size() + 1, 0);
6256 labelList bdyFaceStarts(boundary.size() + 1, 0);
6257 List<Map<label> > bdyFaceIndices(boundary.size() + 1);
6259 // Allocate the faceMap. Since the number of internal faces is unknown,
6260 // detect internal ones first and update boundaries later.
6263 forAllConstIter(Map<label>, rCellMap, cIter)
6265 const cell& thisCell = cells_[cIter.key()];
6267 forAll(thisCell, faceI)
6269 label fIndex = thisCell[faceI];
6271 if (!rFaceMap.found(fIndex))
6273 // Determine the patch index
6274 label patchID = whichPatch(fIndex);
6278 // Internal face. Check if this needs to
6279 // be added to the 'default' patch.
6280 label own = owner_[fIndex];
6281 label nei = neighbour_[fIndex];
6283 if (rCellMap.found(own) && rCellMap.found(nei))
6285 faceMap.insert(nF, fIndex);
6286 rFaceMap.insert(fIndex, nF);
6291 // Update boundary maps.
6292 label bfI = bdyFaceSizes[boundary.size()]++;
6294 // Skip faceMap and update the reverseMap for now.
6295 // faceMap will be updated once all
6296 // internal faces have been detected.
6297 rFaceMap.insert(fIndex, bfI);
6298 bdyFaceIndices[boundary.size()].insert(bfI, fIndex);
6303 // Update boundary maps.
6304 label bfI = bdyFaceSizes[patchID]++;
6306 // Skip faceMap and update the reverseMap for now.
6307 // faceMap will be updated once all
6308 // internal faces have been detected.
6309 rFaceMap.insert(fIndex, bfI);
6310 bdyFaceIndices[patchID].insert(bfI, fIndex);
6313 // Accumulate face sizes
6314 sumNFE += faces_[fIndex].size();
6319 // Set the number of internal faces at this point
6320 cMap.nEntities(coupleMap::INTERNAL_FACE) = nF;
6323 bdyFaceStarts[0] = nF;
6325 for (label i = 1; i < bdyFaceStarts.size(); i++)
6327 bdyFaceStarts[i] = bdyFaceStarts[i-1] + bdyFaceSizes[i-1];
6330 // Update faceMap and reverseFaceMap for boundaries
6331 forAll(bdyFaceIndices, patchI)
6333 label pStart = bdyFaceStarts[patchI];
6335 forAllConstIter(Map<label>, bdyFaceIndices[patchI], fIter)
6337 faceMap.insert(fIter.key() + pStart, fIter());
6338 rFaceMap[fIter()] = fIter.key() + pStart;
6341 // Update face-count
6342 nF += bdyFaceSizes[patchI];
6345 // Keep track of inserted boundary edge indices
6346 // - Add exposed internal edges to a 'default' patch
6347 // at the end of the list.
6348 labelList bdyEdgeSizes(boundary.size() + 1, 0);
6349 labelList bdyEdgeStarts(boundary.size() + 1, 0);
6350 List<Map<label> > bdyEdgeIndices(boundary.size() + 1);
6352 // Allocate the edgeMap. Since the number of internal edges is unknown,
6353 // detect internal ones first and update boundaries later.
6354 forAllConstIter(Map<label>, rFaceMap, fIter)
6356 const labelList& fEdges = faceEdges_[fIter.key()];
6358 forAll(fEdges, edgeI)
6360 label eIndex = fEdges[edgeI];
6362 if (!rEdgeMap.found(eIndex))
6364 // Determine the patch index
6365 label patchID = whichEdgePatch(eIndex);
6369 bool boundaryEdge = false;
6371 // Check if any cells touching edgeFaces
6372 // do not belong to the cellMap.
6373 const labelList& eFaces = edgeFaces_[eIndex];
6375 forAll(eFaces, faceI)
6377 label fIndex = eFaces[faceI];
6379 label own = owner_[fIndex];
6380 label nei = neighbour_[fIndex];
6382 if (!rCellMap.found(own) || !rCellMap.found(nei))
6384 boundaryEdge = true;
6391 // Update boundary maps.
6392 label beI = bdyEdgeSizes[boundary.size()]++;
6394 // Skip edgeMap and update the reverseMap for now.
6395 // edgeMap will be updated once all
6396 // internal edges have been detected.
6397 rEdgeMap.insert(eIndex, beI);
6398 bdyEdgeIndices[boundary.size()].insert(beI, eIndex);
6403 edgeMap.insert(nE, eIndex);
6404 rEdgeMap.insert(eIndex, nE);
6410 // Update boundary maps.
6411 label beI = bdyEdgeSizes[patchID]++;
6413 // Skip edgeMap and update the reverseMap for now.
6414 // edgeMap will be updated once all
6415 // internal edges have been detected.
6416 rEdgeMap.insert(eIndex, beI);
6417 bdyEdgeIndices[patchID].insert(beI, eIndex);
6423 // Set the number of internal edges at this point
6424 cMap.nEntities(coupleMap::INTERNAL_EDGE) = nE;
6427 bdyEdgeStarts[0] = nE;
6429 for (label i = 1; i < bdyEdgeStarts.size(); i++)
6431 bdyEdgeStarts[i] = bdyEdgeStarts[i-1] + bdyEdgeSizes[i-1];
6434 // Update edgeMap and reverseEdgeMap for boundaries
6435 forAll(bdyEdgeIndices, patchI)
6437 label pStart = bdyEdgeStarts[patchI];
6439 forAllConstIter(Map<label>, bdyEdgeIndices[patchI], eIter)
6441 edgeMap.insert(eIter.key() + pStart, eIter());
6442 rEdgeMap[eIter()] = eIter.key() + pStart;
6445 // Update edge-count
6446 nE += bdyEdgeSizes[patchI];
6449 // Set additional points in the pointMap
6450 forAllConstIter(Map<label>, rEdgeMap, eIter)
6452 const edge& thisEdge = edges_[eIter.key()];
6454 if (!rPointMap.found(thisEdge[0]))
6456 pointMap.insert(nP, thisEdge[0]);
6457 rPointMap.insert(thisEdge[0], nP);
6461 if (!rPointMap.found(thisEdge[1]))
6463 pointMap.insert(nP, thisEdge[1]);
6464 rPointMap.insert(thisEdge[1], nP);
6469 // Assign sizes to the mesh
6470 cMap.nEntities(coupleMap::POINT) = nP;
6471 cMap.nEntities(coupleMap::EDGE) = nE;
6472 cMap.nEntities(coupleMap::FACE) = nF;
6473 cMap.nEntities(coupleMap::CELL) = nC;
6474 cMap.nEntities(coupleMap::NFE_SIZE) = sumNFE;
6475 cMap.nEntities(coupleMap::NBDY) = boundary.size() + 1;
6477 // Size up buffers and fill them
6478 cMap.allocateBuffers();
6480 pointField& pBuffer = cMap.pointBuffer();
6481 pointField& opBuffer = cMap.oldPointBuffer();
6483 forAllConstIter(Map<label>, pointMap, pIter)
6485 pBuffer[pIter.key()] = points_[pIter()];
6486 opBuffer[pIter.key()] = oldPoints_[pIter()];
6491 // Edge buffer size: 2 points for every edge
6492 labelList& eBuffer = cMap.entityBuffer(coupleMap::EDGE);
6494 for (label i = 0; i < nE; i++)
6496 label eIndex = edgeMap[i];
6497 const edge& edgeToCheck = edges_[eIndex];
6499 eBuffer[index++] = rPointMap[edgeToCheck[0]];
6500 eBuffer[index++] = rPointMap[edgeToCheck[1]];
6505 labelList& fBuffer = cMap.entityBuffer(coupleMap::FACE);
6506 labelList& fOwner = cMap.entityBuffer(coupleMap::OWNER);
6507 labelList& fNeighbour = cMap.entityBuffer(coupleMap::NEIGHBOUR);
6508 labelList& feBuffer = cMap.entityBuffer(coupleMap::FACE_EDGE);
6510 for (label i = 0; i < nF; i++)
6512 label fIndex = faceMap[i];
6513 label own = owner_[fIndex];
6514 label nei = neighbour_[fIndex];
6516 // Fetch mapped addressing
6517 label fOwn = rCellMap.found(own) ? rCellMap[own] : -1;
6518 label fNei = rCellMap.found(nei) ? rCellMap[nei] : -1;
6522 // Check if this face is pointed the right way
6523 if ((fNei > -1) && (fNei < fOwn))
6525 thisFace = faces_[fIndex].reverseFace();
6527 fNeighbour[i] = fOwn;
6531 thisFace = faces_[fIndex];
6536 fNeighbour[i] = fNei;
6542 // This face is pointed the wrong way.
6543 thisFace = faces_[fIndex].reverseFace();
6547 const labelList& fEdges = faceEdges_[fIndex];
6549 forAll(fEdges, indexI)
6551 fBuffer[index] = rPointMap[thisFace[indexI]];
6552 feBuffer[index] = rEdgeMap[fEdges[indexI]];
6558 // Fill variable size face-list sizes
6559 labelList& nfeBuffer = cMap.entityBuffer(coupleMap::NFE_BUFFER);
6561 forAllConstIter(Map<label>, faceMap, fIter)
6563 nfeBuffer[fIter.key()] = faces_[fIter()].size();
6566 // Fill in boundary information
6567 cMap.entityBuffer(coupleMap::FACE_STARTS) = bdyFaceStarts;
6568 cMap.entityBuffer(coupleMap::FACE_SIZES) = bdyFaceSizes;
6569 cMap.entityBuffer(coupleMap::EDGE_STARTS) = bdyEdgeStarts;
6570 cMap.entityBuffer(coupleMap::EDGE_SIZES) = bdyEdgeSizes;
6572 labelList& ptBuffer = cMap.entityBuffer(coupleMap::PATCH_ID);
6574 // Fill types for all but the last one (which is default).
6575 forAll(boundary, patchI)
6577 if (isA<processorPolyPatch>(boundary[patchI]))
6579 // Fetch the neighbouring processor ID
6580 const processorPolyPatch& pp =
6582 refCast<const processorPolyPatch>(boundary[patchI])
6585 // Set processor patches to a special type
6586 ptBuffer[patchI] = -2 - pp.neighbProcNo();
6590 // Conventional physical patch. Make an identical
6591 // map, since physical boundaries are present on
6593 ptBuffer[patchI] = patchI;
6597 // Fill the default patch as 'internal'
6598 ptBuffer[boundary.size()] = -1;
6600 // Make a temporary dictionary for patch construction
6601 dictionary patchDict;
6603 // Specify the list of patch names and types
6604 wordList patchNames(ptBuffer.size());
6606 forAll(patchNames, patchI)
6608 // Add a new subDictionary
6609 dictionary patchSubDict;
6611 if (patchI == patchNames.size() - 1)
6613 // Artificially set the last patch
6616 patchNames[patchI] = "defaultPatch";
6619 patchSubDict.add("type", "empty");
6622 patchSubDict.add("startFace", bdyFaceStarts[patchI]);
6623 patchSubDict.add("nFaces", bdyFaceSizes[patchI]);
6626 if (ptBuffer[patchI] <= -2)
6628 // Back out the neighbouring processor ID
6629 label neiProcNo = Foam::mag(ptBuffer[patchI] + 2);
6632 patchNames[patchI] =
6635 + Foam::name(Pstream::myProcNo())
6637 + Foam::name(neiProcNo)
6641 patchSubDict.add("type", "subMeshProcessor");
6644 patchSubDict.add("startFace", bdyFaceStarts[patchI]);
6645 patchSubDict.add("nFaces", bdyFaceSizes[patchI]);
6647 // Add processor-specific info
6648 patchSubDict.add("myProcNo", Pstream::myProcNo());
6649 patchSubDict.add("neighbProcNo", neiProcNo);
6654 patchNames[patchI] = boundary[ptBuffer[patchI]].name();
6657 patchSubDict.add("type", boundary[ptBuffer[patchI]].type());
6660 patchSubDict.add("startFace", bdyFaceStarts[patchI]);
6661 patchSubDict.add("nFaces", bdyFaceSizes[patchI]);
6664 // Add subdictionary
6665 patchDict.add(patchNames[patchI], patchSubDict);
6672 new dynamicTopoFvMesh
6677 fvMesh::defaultRegion,
6681 xferCopy(cMap.pointBuffer()),
6682 xferCopy(cMap.oldPointBuffer()),
6683 xferCopy(cMap.edges()),
6684 xferCopy(cMap.faces()),
6685 xferCopy(cMap.faceEdges()),
6686 xferCopy(cMap.owner()),
6687 xferCopy(cMap.neighbour()),
6697 // Set maps as built.
6698 subMesh.setBuiltMaps();
6700 // For debugging purposes...
6703 Pout<< "Writing out sent subMesh for processor: "
6709 + Foam::name(Pstream::myProcNo())
6715 // Write out patch information
6716 Pout<< " faceStarts: " << bdyFaceStarts << nl
6717 << " faceSizes: " << bdyFaceSizes << nl
6718 << " edgeStarts: " << bdyEdgeStarts << nl
6719 << " edgeSizes: " << bdyEdgeSizes << nl
6720 << " patchTypes: " << ptBuffer << endl;
6725 // Build coupled maps for locally coupled patches.
6726 // - Performs a geometric match initially, since the mesh provides
6727 // no explicit information for topological coupling.
6728 // - For each subsequent topology change, maps are updated during
6729 // the re-ordering stage.
6730 void dynamicTopoFvMesh::buildLocalCoupledMaps()
6732 if (patchCoupling_.empty())
6739 Info<< "Building local coupled maps...";
6742 // Check if a geometric tolerance has been specified.
6743 const boundBox& box = polyMesh::bounds();
6745 // Compute tolerance
6746 scalar tol = geomMatchTol_()*box.mag();
6748 const polyBoundaryMesh& boundary = boundaryMesh();
6750 forAll(patchCoupling_, patchI)
6752 if (!patchCoupling_(patchI))
6757 // Build maps only for the first time.
6758 // Information is subsequently mapped for each topo-change.
6759 if (patchCoupling_[patchI].builtMaps())
6764 // Deal with cyclics later
6765 if (isA<cyclicPolyPatch>(boundary[patchI]))
6770 const coupleMap& cMap = patchCoupling_[patchI].map();
6772 // Map points on each patch.
6773 const labelList& mP = boundary[cMap.masterIndex()].meshPoints();
6774 const labelList& sP = boundary[cMap.slaveIndex()].meshPoints();
6776 // Sanity check: Do patches have equal number of entities?
6777 if (mP.size() != sP.size())
6779 FatalErrorIn("void dynamicTopoFvMesh::buildLocalCoupledMaps()")
6780 << "Patch point sizes are not consistent."
6781 << abort(FatalError);
6784 // Check if maps were read from disk.
6785 // If so, geometric checking is unnecessary.
6786 if (cMap.entityMap(coupleMap::POINT).size() == 0)
6788 // Match points geometrically
6789 labelList pointMap(mP.size(), -1);
6795 boundary[cMap.masterIndex()].localPoints(),
6796 boundary[cMap.slaveIndex()].localPoints(),
6797 scalarField(mP.size(), tol),
6805 FatalErrorIn("void dynamicTopoFvMesh::buildLocalCoupledMaps()")
6806 << " Failed to match all points"
6807 << " within a tolerance of: " << tol << nl
6808 << " matchTol: " << geomMatchTol_() << nl
6809 << abort(FatalError);
6815 label masterIndex = mP[indexI];
6816 label slaveIndex = sP[pointMap[indexI]];
6834 const labelListList& mpF = boundary[cMap.masterIndex()].pointFaces();
6835 const labelListList& spF = boundary[cMap.slaveIndex()].pointFaces();
6837 // Fetch reference to maps
6838 const Map<label>& pMap = cMap.entityMap(coupleMap::POINT);
6839 const Map<label>& eMap = cMap.entityMap(coupleMap::EDGE);
6840 const Map<label>& fMap = cMap.entityMap(coupleMap::FACE);
6842 // Now that all points are matched,
6843 // perform topological matching for higher entities.
6846 // Fetch global point indices
6847 label mp = mP[indexI];
6848 label sp = pMap[mp];
6850 // Fetch local point indices
6851 label lmp = boundary[cMap.masterIndex()].whichPoint(mp);
6852 label lsp = boundary[cMap.slaveIndex()].whichPoint(sp);
6854 // Fetch patch starts
6855 label mStart = boundary[cMap.masterIndex()].start();
6856 label sStart = boundary[cMap.slaveIndex()].start();
6858 // Match faces for both 2D and 3D.
6859 const labelList& mpFaces = mpF[lmp];
6860 const labelList& spFaces = spF[lsp];
6862 forAll(mpFaces, faceI)
6864 // Fetch the global face index
6865 label mfIndex = (mStart + mpFaces[faceI]);
6867 if (fMap.found(mfIndex))
6872 const face& mFace = faces_[mfIndex];
6876 // Set up a comparison face.
6879 // Configure the face for comparison.
6880 forAll(mFace, pointI)
6882 cFace[pointI] = pMap[mFace[pointI]];
6885 bool matched = false;
6887 forAll(spFaces, faceJ)
6889 // Fetch the global face index
6890 label sfIndex = (sStart + spFaces[faceJ]);
6892 const face& sFace = faces_[sfIndex];
6894 if (face::compare(cFace, sFace))
6896 // Found the slave. Add a map entry
6921 "void dynamicTopoFvMesh::buildLocalCoupledMaps()"
6923 << " Failed to match face: "
6924 << mfIndex << ": " << mFace << nl
6925 << " Comparison face: " << cFace
6926 << abort(FatalError);
6931 label slaveFaceIndex = -1;
6933 // Set up a comparison face.
6941 forAll(spFaces, faceJ)
6943 // Fetch the global face index
6944 label sfIndex = (sStart + spFaces[faceJ]);
6946 const face& sFace = faces_[sfIndex];
6948 if (triFace::compare(cFace, triFace(sFace)))
6950 // Found the slave. Add a map entry
6965 slaveFaceIndex = sfIndex;
6971 if (slaveFaceIndex == -1)
6975 "void dynamicTopoFvMesh::buildLocalCoupledMaps()"
6977 << " Failed to match face: "
6978 << mfIndex << ": " << mFace << nl
6979 << " Comparison face: " << cFace
6980 << abort(FatalError);
6983 // Match all edges on this face as well.
6984 const labelList& mfEdges = faceEdges_[mfIndex];
6985 const labelList& sfEdges = faceEdges_[slaveFaceIndex];
6987 forAll(mfEdges, edgeI)
6989 if (eMap.found(mfEdges[edgeI]))
6994 const edge& mEdge = edges_[mfEdges[edgeI]];
6996 // Configure a comparison edge.
6997 edge cEdge(pMap[mEdge[0]], pMap[mEdge[1]]);
6999 bool matchedEdge = false;
7001 forAll(sfEdges, edgeJ)
7003 const edge& sEdge = edges_[sfEdges[edgeJ]];
7007 // Found the slave. Add a map entry
7032 "void dynamicTopoFvMesh::"
7033 "buildLocalCoupledMaps()"
7035 << " Failed to match edge: "
7036 << mfEdges[edgeI] << ": "
7038 << " cEdge: " << cEdge
7039 << abort(FatalError);
7046 // Set maps as built
7047 patchCoupling_[patchI].setBuiltMaps();
7050 // Build maps for cyclics
7051 forAll(patchCoupling_, patchI)
7053 if (!patchCoupling_(patchI))
7058 // Build maps only for the first time.
7059 // Information is subsequently mapped for each topo-change.
7060 if (patchCoupling_[patchI].builtMaps())
7065 // Skip if not cyclic
7066 if (!isA<cyclicPolyPatch>(boundary[patchI]))
7071 const coupleMap& cMap = patchCoupling_[patchI].map();
7073 // Match faces and points in one go
7074 label halfSize = boundary[patchI].size() / 2;
7075 label patchStart = boundary[patchI].start();
7077 // Fetch reference to map
7078 const Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
7080 for (label i = 0; i < halfSize; i++)
7082 label half0Index = (i + patchStart);
7083 label half1Index = (i + halfSize + patchStart);
7100 const face& half0Face = faces_[half0Index];
7101 const face& half1Face = faces_[half1Index];
7103 label fS = half0Face.size();
7105 forAll(half0Face, pointI)
7107 if (pointMap.found(half0Face[pointI]))
7112 label masterIndex = half0Face[pointI];
7113 label slaveIndex = half1Face[(fS - pointI) % fS];
7136 // Fetch reference to map
7137 const Map<label>& edgeMap = cMap.entityMap(coupleMap::EDGE);
7139 const labelList& f0Edges = faceEdges_[half0Index];
7140 const labelList& f1Edges = faceEdges_[half1Index];
7142 forAll(f0Edges, edgeI)
7144 if (edgeMap.found(f0Edges[edgeI]))
7149 const edge& mEdge = edges_[f0Edges[edgeI]];
7151 // Build a comparison edge
7152 edge cEdge(pointMap[mEdge[0]], pointMap[mEdge[1]]);
7154 forAll(f1Edges, edgeJ)
7156 const edge& sEdge = edges_[f1Edges[edgeJ]];
7184 Info<< "Done." << endl;
7189 // Build coupled maps for coupled processor patches
7190 void dynamicTopoFvMesh::buildProcessorCoupledMaps()
7192 if (procIndices_.empty())
7199 const polyBoundaryMesh& boundary = boundaryMesh();
7201 // Build a list of nearest neighbours.
7202 forAll(boundary, patchI)
7204 if (isA<processorPolyPatch>(boundary[patchI]))
7206 const processorPolyPatch& pp =
7208 refCast<const processorPolyPatch>(boundary[patchI])
7211 nPrc.insert(pp.neighbProcNo(), patchI);
7215 // Wait for all transfers to complete.
7216 meshOps::waitForBuffers();
7218 // Put un-matched faces in a list.
7219 labelHashSet unMatchedFaces;
7221 forAll(procIndices_, pI)
7223 label proc = procIndices_[pI];
7225 coupledMesh& rPM = recvMeshes_[pI];
7226 const coupleMap& cMap = rPM.map();
7227 const labelList& ptBuffer = cMap.entityBuffer(coupleMap::PATCH_ID);
7229 // Fetch reference to patch starts / sizes
7230 const labelList& faceStarts = cMap.entityBuffer(coupleMap::FACE_STARTS);
7231 const labelList& faceSizes = cMap.entityBuffer(coupleMap::FACE_SIZES);
7232 const labelList& edgeStarts = cMap.entityBuffer(coupleMap::EDGE_STARTS);
7233 const labelList& edgeSizes = cMap.entityBuffer(coupleMap::EDGE_SIZES);
7235 // Make a temporary dictionary for patch construction
7236 dictionary patchDict;
7238 // Specify the list of patch names and types
7239 wordList patchNames(ptBuffer.size());
7241 forAll(patchNames, patchI)
7243 // Add a new subDictionary
7244 dictionary patchSubDict;
7246 if (patchI == patchNames.size() - 1)
7248 // Artificially set the last patch
7251 patchNames[patchI] = "defaultPatch";
7254 patchSubDict.add("type", "empty");
7257 patchSubDict.add("startFace", faceStarts[patchI]);
7258 patchSubDict.add("nFaces", faceSizes[patchI]);
7261 if (ptBuffer[patchI] <= -2)
7263 // Back out the neighbouring processor ID
7264 label neiProcNo = Foam::mag(ptBuffer[patchI] + 2);
7267 patchNames[patchI] =
7272 + Foam::name(neiProcNo)
7276 patchSubDict.add("type", "subMeshProcessor");
7279 patchSubDict.add("startFace", faceStarts[patchI]);
7280 patchSubDict.add("nFaces", faceSizes[patchI]);
7282 // Add processor-specific info
7283 patchSubDict.add("myProcNo", proc);
7284 patchSubDict.add("neighbProcNo", neiProcNo);
7289 patchNames[patchI] = boundary[ptBuffer[patchI]].name();
7292 patchSubDict.add("type", boundary[ptBuffer[patchI]].type());
7295 patchSubDict.add("startFace", faceStarts[patchI]);
7296 patchSubDict.add("nFaces", faceSizes[patchI]);
7299 // Add subdictionary
7300 patchDict.add(patchNames[patchI], patchSubDict);
7307 new dynamicTopoFvMesh
7312 fvMesh::defaultRegion,
7316 xferCopy(cMap.pointBuffer()),
7317 xferCopy(cMap.oldPointBuffer()),
7318 xferCopy(cMap.edges()),
7319 xferCopy(cMap.faces()),
7320 xferCopy(cMap.faceEdges()),
7321 xferCopy(cMap.owner()),
7322 xferCopy(cMap.neighbour()),
7334 Pout<< "Writing out received subMesh for processor: "
7337 rPM.subMesh().writeVTK
7340 + Foam::name(Pstream::myProcNo())
7343 identity(cMap.nEntities(coupleMap::CELL)),
7348 // First, topologically map points based on subMeshPoints.
7349 const labelList& mP = cMap.subMeshPoints();
7351 label nShared = cMap.nEntities(coupleMap::SHARED_POINT);
7352 label nGlobal = cMap.nEntities(coupleMap::GLOBAL_POINT);
7354 // Sanity check: Do sub-mesh point sizes match?
7355 if (mP.size() != nShared)
7357 FatalErrorIn("void dynamicTopoFvMesh::buildProcessorCoupledMaps()")
7358 << " Sub-mesh point sizes don't match." << nl
7359 << " My procID: " << Pstream::myProcNo() << nl
7360 << " Neighbour processor: " << proc << nl
7361 << " mP size: " << mP.size() << nl
7362 << " nShared: " << nShared << nl
7363 << abort(FatalError);
7366 if (nPrc.found(proc))
7368 // This is an immediate neighbour.
7369 // All patch points must be matched.
7371 // Fetch addressing for this patch.
7372 const processorPolyPatch& pp =
7374 refCast<const processorPolyPatch>(boundary[nPrc[proc]])
7377 const labelList& neiPoints = pp.neighbPoints();
7381 // Find all occurrences of multiply connected points
7382 labelList mIdx = findIndices(neiPoints, -1);
7386 // Write out for post-processing
7387 UIndirectList<label> myIdx(mP, mIdx);
7388 writeVTK("mcPoints", labelList(myIdx), 0);
7392 "void dynamicTopoFvMesh::buildProcessorCoupledMaps()"
7394 << " Multiply connected point." << nl
7395 << " My procID: " << Pstream::myProcNo() << nl
7396 << " Neighbour processor: " << proc << nl
7397 << " Neighbour Indices: " << mIdx << nl
7398 << " My indices: " << myIdx << nl
7399 << abort(FatalError);
7422 // Prepare point maps for globally shared points
7423 const labelList& gpB = cMap.entityBuffer(coupleMap::POINT);
7425 // Sanity check: Do global point buffer sizes match?
7426 if ((mP.size() + gpB.size()) != nGlobal)
7428 FatalErrorIn("void dynamicTopoFvMesh::buildProcessorCoupledMaps()")
7429 << " Global point sizes don't match"
7430 << " for processor: " << proc << nl
7431 << " mP size: " << mP.size() << nl
7432 << " gpB size: " << gpB.size() << nl
7433 << " Total: " << (mP.size() + gpB.size()) << nl
7434 << " nGlobal: " << nGlobal << nl
7435 << abort(FatalError);
7438 // Look at globalPoint addressing for its information.
7439 const globalMeshData& gD = polyMesh::globalData();
7440 const labelList& spA = gD.sharedPointAddr();
7441 const labelList& spL = gD.sharedPointLabels();
7445 // Find my equivalent global point
7446 label maIndex = findIndex(spA, gpB[pointI]);
7466 const boundBox& box = polyMesh::bounds();
7467 const pointField& sPoints = rPM.subMesh().points_;
7468 const Map<label>& pMap = cMap.entityMap(coupleMap::POINT);
7470 // Compute tolerance
7471 scalar tol = geomMatchTol_()*box.mag();
7473 forAllConstIter(Map<label>, pMap, pIter)
7475 scalar dist = mag(points_[pIter.key()] - sPoints[pIter()]);
7481 "void dynamicTopoFvMesh::buildProcessorCoupledMaps()"
7483 << " Failed to match point: " << pIter.key()
7484 << ": " << points_[pIter.key()]
7485 << " with point: " << pIter()
7486 << ": " << sPoints[pIter()] << nl
7487 << " Missed by: " << dist << nl
7488 << " Tolerance: " << tol << nl
7489 << abort(FatalError);
7494 // Set up a comparison face.
7497 // Now match all faces connected to master points.
7498 if (nPrc.found(proc))
7500 // This is an immediate neighbour.
7501 label mStart = boundary[nPrc[proc]].start();
7502 label mSize = boundary[nPrc[proc]].size();
7504 // Abbreviate for convenience
7505 const dynamicTopoFvMesh& sMesh = rPM.subMesh();
7506 const Map<label>& pMap = cMap.entityMap(coupleMap::POINT);
7507 const Map<label>& eMap = cMap.entityMap(coupleMap::EDGE);
7508 const Map<label>& fMap = cMap.entityMap(coupleMap::FACE);
7510 // Fetch global pointFaces for the slave.
7511 const faceList& slaveFaces = sMesh.faces();
7512 const labelListList& spF = sMesh.pointFaces();
7514 // Match patch faces for both 2D and 3D.
7515 for (label i = 0; i < mSize; i++)
7517 // Fetch the global face index
7518 label mfIndex = (mStart + i);
7520 if (fMap.found(mfIndex))
7525 const face& mFace = faces_[mfIndex];
7527 // Configure the face for comparison.
7528 forAll(mFace, pointI)
7530 cFace[pointI] = pMap[mFace[pointI]];
7533 // Fetch pointFaces for the zeroth point.
7534 const labelList& spFaces = spF[cFace[0]];
7536 label sFaceIndex = -1, compareVal = -1;
7538 forAll(spFaces, faceJ)
7540 // Fetch the global face index
7541 label sfIndex = spFaces[faceJ];
7543 const face& sFace = slaveFaces[sfIndex];
7545 // Check for triangle face optimization
7546 if (mFace.size() == 3)
7552 triFace(cFace[0], cFace[1], cFace[2]),
7553 triFace(sFace[0], sFace[1], sFace[2])
7559 compareVal = face::compare(cFace, sFace);
7564 // Found the slave. Add a map entry
7579 sFaceIndex = sfIndex;
7585 if (sFaceIndex == -1)
7587 unMatchedFaces.insert(mfIndex);
7592 // No need to match edges in 2D
7598 // Match all edges on this face as well.
7599 const labelList& mfEdges = faceEdges_[mfIndex];
7600 const labelList& sfEdges = sMesh.faceEdges_[sFaceIndex];
7602 forAll(mfEdges, edgeI)
7604 if (eMap.found(mfEdges[edgeI]))
7609 const edge& mEdge = edges_[mfEdges[edgeI]];
7611 // Configure a comparison edge.
7612 edge cEdge(pMap[mEdge[0]], pMap[mEdge[1]]);
7614 bool matchedEdge = false;
7616 forAll(sfEdges, edgeJ)
7618 const edge& sEdge = sMesh.edges_[sfEdges[edgeJ]];
7622 // Found the slave. Add a map entry
7645 // Write out the edge
7646 writeVTK("mEdge", mfEdges[edgeI], 1);
7650 "void dynamicTopoFvMesh::"
7651 "buildProcessorCoupledMaps()"
7653 << " Failed to match edge: "
7654 << mfEdges[edgeI] << ": "
7656 << " cEdge: " << cEdge
7657 << " for processor: " << proc
7658 << abort(FatalError);
7664 // Attempt to match other edges in 3D, provided any common ones exist.
7665 // - This will handle point-only coupling situations for edges.
7668 // Fetch Maps for this processor
7669 const Map<label>& edgeMap = cMap.entityMap(coupleMap::EDGE);
7670 const Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
7672 forAllConstIter(Map<label>, pointMap, pIter)
7674 const labelList& pEdges = pointEdges_[pIter.key()];
7676 forAll(pEdges, edgeI)
7678 const label eIndex = pEdges[edgeI];
7680 // Only pick boundary edges
7681 if (whichEdgePatch(eIndex) == -1)
7686 // Skip mapped edges
7687 if (edgeMap.found(eIndex))
7692 const edge& eCheck = edges_[eIndex];
7694 // Check if a map exists for the other point
7700 eCheck.otherVertex(pIter.key())
7706 const dynamicTopoFvMesh& mesh = rPM.subMesh();
7708 // Configure a check edge
7709 edge cEdge(sIndex, pIter());
7711 const labelList& spEdges = mesh.pointEdges_[sIndex];
7713 forAll(spEdges, edgeJ)
7715 const edge& sEdge = mesh.edges_[spEdges[edgeJ]];
7719 // Found the slave. Add a map entry
7741 // Prepare processor point and edge maps
7742 // - For each point / edge on the subMesh, create
7743 // a list of processors (other than this one)
7744 // that are connected to it
7745 Map<labelList>& pEdgeMap = cMap.subMeshEdgeMap();
7746 Map<labelList>& pPointMap = cMap.subMeshPointMap();
7748 const dynamicTopoFvMesh& sMesh = rPM.subMesh();
7749 const polyBoundaryMesh& bdy = sMesh.boundaryMesh();
7753 if (!isA<processorPolyPatch>(bdy[pI]))
7758 const processorPolyPatch& pp =
7760 refCast<const processorPolyPatch>(bdy[pI])
7763 label neiProcNo = pp.neighbProcNo();
7764 label myProcNo = Pstream::myProcNo();
7766 if (neiProcNo == myProcNo)
7771 label mStart = pp.start();
7772 label mSize = pp.size();
7774 for (label i = 0; i < mSize; i++)
7776 const face& f = sMesh.faces_[i + mStart];
7777 const labelList& fe = sMesh.faceEdges_[i + mStart];
7781 const label pIndex = f[j];
7782 const label eIndex = fe[j];
7784 Map<labelList>::iterator pIt = pPointMap.find(pIndex);
7786 if (pIt == pPointMap.end())
7788 pPointMap.insert(pIndex, labelList(1, neiProcNo));
7792 if (findIndex(pIt(), neiProcNo) == -1)
7794 meshOps::sizeUpList(neiProcNo, pIt());
7798 Map<labelList>::iterator eIt = pEdgeMap.find(eIndex);
7800 if (eIt == pEdgeMap.end())
7802 pEdgeMap.insert(eIndex, labelList(1, neiProcNo));
7806 if (findIndex(eIt(), neiProcNo) == -1)
7808 meshOps::sizeUpList(neiProcNo, eIt());
7816 if (unMatchedFaces.size())
7818 // Write out the face
7819 writeVTK("mFaces", unMatchedFaces.toc(), 2);
7823 "void dynamicTopoFvMesh::buildProcessorCoupledMaps()"
7825 << " Unmatched faces were found for processor: " << proc
7826 << abort(FatalError);
7830 // Use subMesh and global points to identify additional
7831 // processors, and update maps accordingly
7834 forAll(procIndices_, pI)
7836 const label neiProcNo = procIndices_[pI];
7837 const coupleMap& cMi = recvMeshes_[pI].map();
7839 const labelList& smPts = cMi.subMeshPoints();
7840 const List<labelPair>& gpp = cMi.globalProcPoints();
7842 forAll(procIndices_, pJ)
7849 const coupleMap& cMj = recvMeshes_[pJ].map();
7851 Map<labelList>& pPointMap = cMj.subMeshPointMap();
7853 forAll(smPts, pointI)
7855 label mP = smPts[pointI];
7856 label sP = cMj.findSlave(coupleMap::POINT, mP);
7863 Map<labelList>::iterator pIt = pPointMap.find(sP);
7865 if (pIt == pPointMap.end())
7867 pPointMap.insert(sP, labelList(1, neiProcNo));
7871 if (findIndex(pIt(), neiProcNo) == -1)
7873 meshOps::sizeUpList(neiProcNo, pIt());
7880 label mP = gpp[pointI].first();
7881 label sP = cMj.findSlave(coupleMap::POINT, mP);
7888 Map<labelList>::iterator pIt = pPointMap.find(sP);
7890 if (pIt == pPointMap.end())
7892 pPointMap.insert(sP, labelList(1, neiProcNo));
7896 if (findIndex(pIt(), neiProcNo) == -1)
7898 meshOps::sizeUpList(neiProcNo, pIt());
7908 // Introduce a new processor patch to the mesh
7909 label dynamicTopoFvMesh::createProcessorPatch(const label proc)
7911 if (!Pstream::parRun())
7916 // Find index in list of processors
7917 label pI = findIndex(procIndices_, proc);
7919 // Check if an entry already exists
7922 label patchIndex = sendMeshes_[pI].map().patchIndex();
7924 if (patchIndex > -1)
7930 // Get the new patch index,
7931 // and increment the number of patches
7932 label patchID = nPatches_++;
7934 // Set the new patch index in patchMaps
7939 pI = procIndices_.size();
7941 procIndices_.setSize(pI + 1);
7942 sendMeshes_.setSize(pI + 1);
7945 // Create a basic entry on the subMesh
7946 procIndices_[pI] = proc;
7953 *this, // Reference to this mesh
7956 true, // Sent to neighbour
7957 patchID, // Patch index
7958 proc, // Master index
7967 Pout<< " Could not find index for processor: " << proc
7968 << " in indices: " << procIndices_
7969 << abort(FatalError);
7972 sendMeshes_[pI].map().patchIndex() = patchID;
7973 recvMeshes_[pI].map().patchIndex() = patchID;
7976 // Size up patches, and copy old information
7977 label prevPatchID = patchID - 1;
7979 patchSizes_.setSize(nPatches_, 0);
7980 oldPatchSizes_.setSize(nPatches_, 0);
7982 patchStarts_.setSize
7985 patchStarts_[prevPatchID] + patchSizes_[prevPatchID]
7987 oldPatchStarts_.setSize
7990 oldPatchStarts_[prevPatchID] + oldPatchSizes_[prevPatchID]
7993 edgePatchSizes_.setSize(nPatches_, 0);
7994 oldEdgePatchSizes_.setSize(nPatches_, 0);
7996 edgePatchStarts_.setSize
7999 edgePatchStarts_[prevPatchID] + edgePatchSizes_[prevPatchID]
8001 oldEdgePatchStarts_.setSize
8004 oldEdgePatchStarts_[prevPatchID] + oldEdgePatchSizes_[prevPatchID]
8007 patchNMeshPoints_.setSize(nPatches_, 0);
8008 oldPatchNMeshPoints_.setSize(nPatches_, 0);
8012 Pout<< " dynamicTopoFvMesh::createProcessorPatch :"
8013 << " Created new patch for processor: " << proc << nl
8014 << " On subMesh: " << isSubMesh_ << nl
8015 << " pI: " << pI << nl
8016 << " patchID: " << patchID << nl
8017 << " oldPatchStarts: " << oldPatchStarts_ << nl
8018 << " oldPatchSizes: " << oldPatchSizes_ << nl
8019 << " patchStarts: " << patchStarts_ << nl
8020 << " patchSizes: " << patchSizes_ << nl
8024 // Return the new patch index
8029 // If a patch is of processor type, get the neighbouring processor ID
8030 label dynamicTopoFvMesh::getNeighbourProcessor(const label patch) const
8032 if (!Pstream::parRun())
8037 label neiProcNo = -1;
8039 const polyBoundaryMesh& boundary = boundaryMesh();
8046 if (patch < boundary.size())
8048 if (isA<processorPolyPatch>(boundary[patch]))
8050 const processorPolyPatch& pp =
8052 refCast<const processorPolyPatch>(boundary[patch])
8055 // Set the neighbour processor ID
8056 neiProcNo = pp.neighbProcNo();
8061 // New processor patch
8063 // Find the neighbour processor ID
8064 forAll(procIndices_, pI)
8066 label patchID = sendMeshes_[pI].map().patchIndex();
8068 if (patch == patchID)
8070 neiProcNo = procIndices_[pI];
8075 if (neiProcNo == -1)
8077 // An index should have been defined by now
8078 Pout<< " Patch: " << patch
8079 << " was not defined on patch subMeshes."
8080 << abort(FatalError);
8088 // If the number of patches have changed during run-time,
8089 // reset boundaries with new processor patches
8090 void dynamicTopoFvMesh::resetBoundaries()
8092 // Prepare a new list of patches
8093 List<polyPatch*> patches(nPatches_);
8095 // Fetch reference to existing boundary
8096 // - The removeBoundary member merely resets
8097 // boundary size, so this reference is safe
8098 const polyBoundaryMesh& polyBoundary = boundaryMesh();
8100 // Copy all existing patches first
8101 label nOldPatches = polyBoundary.size();
8103 for (label patchI = 0; patchI < nOldPatches; patchI++)
8106 patches[patchI] = polyBoundary[patchI].clone(polyBoundary).ptr();
8109 // Create new processor patches
8110 for (label patchI = nOldPatches; patchI < nPatches_; patchI++)
8112 // Make a temporary dictionary for patch construction
8113 dictionary patchDict;
8115 // Back out the neighbour processor ID
8116 label neiProcNo = getNeighbourProcessor(patchI);
8118 // Add relevant info
8119 patchDict.add("type", "processor");
8120 patchDict.add("startFace", oldPatchStarts_[patchI]);
8121 patchDict.add("nFaces", oldPatchSizes_[patchI]);
8122 patchDict.add("myProcNo", Pstream::myProcNo());
8123 patchDict.add("neighbProcNo", neiProcNo);
8131 + Foam::name(Pstream::myProcNo())
8133 + Foam::name(neiProcNo),
8141 // Remove the old boundary
8142 fvMesh::removeFvBoundary();
8144 // Add patches, but don't calculate geometry, etc
8145 fvMesh::addFvPatches(patches, false);
8147 // Since all fvPatches in fvBoundaryMesh have been reset,
8148 // fvPatch references in all fvPatchFields are no longer
8149 // valid, and must therefore be remapped.
8150 const fvBoundaryMesh& bdy = fvMesh::boundary();
8152 coupledMesh::resizeBoundaries<volScalarField>(*this, bdy);
8153 coupledMesh::resizeBoundaries<volVectorField>(*this, bdy);
8154 coupledMesh::resizeBoundaries<volSphericalTensorField>(*this, bdy);
8155 coupledMesh::resizeBoundaries<volSymmTensorField>(*this, bdy);
8156 coupledMesh::resizeBoundaries<volTensorField>(*this, bdy);
8158 coupledMesh::resizeBoundaries<surfaceScalarField>(*this, bdy);
8159 coupledMesh::resizeBoundaries<surfaceVectorField>(*this, bdy);
8160 coupledMesh::resizeBoundaries<surfaceSphericalTensorField>(*this, bdy);
8161 coupledMesh::resizeBoundaries<surfaceSymmTensorField>(*this, bdy);
8162 coupledMesh::resizeBoundaries<surfaceTensorField>(*this, bdy);
8166 // Convenience macro for field sub-setting
8167 #define sendFieldsOfType(type, info, names, types, offset, map, stream) \
8170 scalar zeroValue = pTraits<scalar>::zero; \
8171 info.send<type##ScalarField> \
8172 (names[0 + offset], types[0 + offset], zeroValue, map, stream); \
8175 vector zeroValue = pTraits<vector>::zero; \
8176 info.send<type##VectorField> \
8177 (names[1 + offset], types[1 + offset], zeroValue, map, stream); \
8180 sphericalTensor zeroValue = pTraits<sphericalTensor>::zero; \
8181 info.send<type##SphericalTensorField> \
8182 (names[2 + offset], types[2 + offset], zeroValue, map, stream); \
8185 symmTensor zeroValue = pTraits<symmTensor>::zero; \
8186 info.send<type##SymmTensorField> \
8187 (names[3 + offset], types[3 + offset], zeroValue, map, stream); \
8190 tensor zeroValue = pTraits<tensor>::zero; \
8191 info.send<type##TensorField> \
8192 (names[4 + offset], types[4 + offset], zeroValue, map, stream); \
8197 // Initialize subMesh field transfers for mapping
8198 void dynamicTopoFvMesh::initFieldTransfers
8201 List<wordList>& names,
8202 List<List<char> >& sendBuffer,
8203 List<List<char> >& recvBuffer
8206 if (!Pstream::parRun())
8213 Info<< " void dynamicTopoFvMesh::initFieldTransfers() :"
8214 << " Initializing subMesh field transfers."
8218 // Clear out mesh geometry, since those
8219 // are to be wiped out after topo-changes anyway.
8221 polyMesh::resetMotion();
8223 // Size up wordLists
8224 // - Five templated volFields
8225 // - Five templated surfaceFields
8229 // Fill in field-types
8230 types[0] = volScalarField::typeName;
8231 types[1] = volVectorField::typeName;
8232 types[2] = volSphericalTensorField::typeName;
8233 types[3] = volSymmTensorField::typeName;
8234 types[4] = volTensorField::typeName;
8236 types[5] = surfaceScalarField::typeName;
8237 types[6] = surfaceVectorField::typeName;
8238 types[7] = surfaceSphericalTensorField::typeName;
8239 types[8] = surfaceSymmTensorField::typeName;
8240 types[9] = surfaceTensorField::typeName;
8242 // Send / recv buffers for field names
8243 List<char> fieldNameSendBuffer, fieldNameRecvBuffer;
8245 if (Pstream::master())
8247 PtrList<OStringStream> fieldNameStream(1);
8250 fieldNameStream.set(0, new OStringStream(IOstream::BINARY));
8252 // Alias for convenience
8253 OStringStream& fNStream = fieldNameStream[0];
8255 // Fetch field-names by type
8256 forAll(types, typeI)
8258 // Get all fields of type
8259 names[typeI] = objectRegistry::names(types[typeI]);
8262 // Send field names to Ostream
8265 // Size up buffers and fill contents
8266 string contents = fNStream.str();
8267 const char* ptr = contents.data();
8269 fieldNameSendBuffer.setSize(contents.size());
8271 forAll(fieldNameSendBuffer, i)
8273 fieldNameSendBuffer[i] = *ptr++;
8277 fieldNameStream.set(0, NULL);
8281 Info<< " Registered fields: " << names << endl;
8284 // Send names to slaves
8285 for (label proc = 1; proc < Pstream::nProcs(); proc++)
8287 // Send buffer to processor
8288 meshOps::pWrite(proc, fieldNameSendBuffer.size());
8289 meshOps::pWrite(proc, fieldNameSendBuffer);
8294 // Receive names from master
8295 label recvBufferSize = -1;
8296 meshOps::pRead(Pstream::masterNo(), recvBufferSize);
8298 // Size up buffer and schedule receive
8299 fieldNameRecvBuffer.setSize(recvBufferSize);
8300 meshOps::pRead(Pstream::masterNo(), fieldNameRecvBuffer);
8303 // Wait for transfers to complete
8304 meshOps::waitForBuffers();
8306 // Convert names from stringStream
8307 if (!Pstream::master())
8311 fieldNameRecvBuffer.begin(),
8312 fieldNameRecvBuffer.size()
8316 IStringStream fieldNameStream(contents, IOstream::BINARY);
8318 // Get names from stream
8319 fieldNameStream >> names;
8322 label nProcs = procIndices_.size();
8325 sendBuffer.setSize(nProcs);
8326 recvBuffer.setSize(nProcs);
8328 // Size up the send stringStream
8329 PtrList<OStringStream> stream(nProcs);
8331 // Now fill in subMesh fields
8332 forAll(procIndices_, pI)
8334 label proc = procIndices_[pI];
8336 const coupledMesh& cInfo = sendMeshes_[pI];
8338 // Initialize stream
8339 stream.set(pI, new OStringStream(IOstream::BINARY));
8341 // Subset and send volFields to stream
8342 const labelList& cMap = cInfo.map().cellMap();
8344 sendFieldsOfType(vol, cInfo, names, types, 0, cMap, stream[pI]);
8346 // Subset and send surfaceFields to stream
8347 const labelList& fMap = cInfo.map().internalFaceMap();
8349 sendFieldsOfType(surface, cInfo, names, types, 5, fMap, stream[pI]);
8351 // Size up buffers and fill contents
8352 string contents = stream[pI].str();
8353 const char* ptr = contents.data();
8355 sendBuffer[pI].setSize(contents.size());
8357 forAll(sendBuffer[pI], i)
8359 sendBuffer[pI][i] = *ptr++;
8363 stream.set(pI, NULL);
8365 // Send buffer to processor
8366 meshOps::pWrite(proc, sendBuffer[pI].size());
8367 meshOps::pWrite(proc, sendBuffer[pI]);
8369 // Receive buffer from processor
8370 label recvBufferSize = -1;
8371 meshOps::pRead(proc, recvBufferSize);
8373 // Size up buffer and schedule receive
8374 recvBuffer[pI].setSize(recvBufferSize);
8375 meshOps::pRead(proc, recvBuffer[pI]);
8378 // We won't wait for all transfers to complete for the moment.
8379 // Meanwhile, do some other useful work, if possible.
8383 // Synchronize field transfers for mapping
8384 void dynamicTopoFvMesh::syncFieldTransfers
8387 List<wordList>& names,
8388 List<List<char> >& recvBuffer
8391 if (!Pstream::parRun())
8398 Info<< " void dynamicTopoFvMesh::syncFieldTransfers() :"
8399 << " Synchronizing subMesh field transfers."
8403 // Wait for all transfers to complete.
8404 meshOps::waitForBuffers();
8406 label nProcs = procIndices_.size();
8408 // Size up stringStream
8409 PtrList<IStringStream> stream(nProcs);
8411 // Size up field PtrLists
8414 List<PtrList<volScalarField> > vsF(nProcs);
8415 List<PtrList<volVectorField> > vvF(nProcs);
8416 List<PtrList<volSphericalTensorField> > vsptF(nProcs);
8417 List<PtrList<volSymmTensorField> > vsytF(nProcs);
8418 List<PtrList<volTensorField> > vtF(nProcs);
8421 List<PtrList<surfaceScalarField> > ssF(nProcs);
8422 List<PtrList<surfaceVectorField> > svF(nProcs);
8423 List<PtrList<surfaceSphericalTensorField> > ssptF(nProcs);
8424 List<PtrList<surfaceSymmTensorField> > ssytF(nProcs);
8425 List<PtrList<surfaceTensorField> > stF(nProcs);
8427 const polyBoundaryMesh& polyBoundary = boundaryMesh();
8429 // Determine number of physical (non-processor) patches
8430 label nPhysical = 0;
8432 forAll(polyBoundary, patchI)
8434 if (isA<processorPolyPatch>(polyBoundary[patchI]))
8442 // Keep track of extra / total entities
8443 label nTotalCells = nOldCells_, nTotalIntFaces = nOldInternalFaces_;
8444 labelList nTotalPatchFaces(SubList<label>(oldPatchSizes_, nPhysical));
8446 // Allocate reverse maps
8447 List<labelList> irvMaps(nProcs), irsMaps(nProcs);
8448 List<labelListList> brMaps(nProcs, labelListList(nPhysical));
8450 forAll(procIndices_, pI)
8452 const coupledMesh& cInfo = recvMeshes_[pI];
8454 // Convert buffer to string
8455 string contents(recvBuffer[pI].begin(), recvBuffer[pI].size());
8457 // Initialize stream
8458 stream.set(pI, new IStringStream(contents, IOstream::BINARY));
8460 // Construct dictionary from stream
8461 dictionary dict(stream[pI]);
8463 // Count the number of additional entities
8464 const coupleMap& cMap = cInfo.map();
8466 // Fetch size from subMesh
8467 label nCells = cMap.nEntities(coupleMap::CELL);
8468 label nIntFaces = cMap.nEntities(coupleMap::INTERNAL_FACE);
8470 // Set field pointers
8471 cInfo.setField(names[0], dict.subDict(types[0]), nCells, vsF[pI]);
8472 cInfo.setField(names[1], dict.subDict(types[1]), nCells, vvF[pI]);
8473 cInfo.setField(names[2], dict.subDict(types[2]), nCells, vsptF[pI]);
8474 cInfo.setField(names[3], dict.subDict(types[3]), nCells, vsytF[pI]);
8475 cInfo.setField(names[4], dict.subDict(types[4]), nCells, vtF[pI]);
8477 cInfo.setField(names[5], dict.subDict(types[5]), nIntFaces, ssF[pI]);
8478 cInfo.setField(names[6], dict.subDict(types[6]), nIntFaces, svF[pI]);
8479 cInfo.setField(names[7], dict.subDict(types[7]), nIntFaces, ssptF[pI]);
8480 cInfo.setField(names[8], dict.subDict(types[8]), nIntFaces, ssytF[pI]);
8481 cInfo.setField(names[9], dict.subDict(types[9]), nIntFaces, stF[pI]);
8483 // Set rmap for this processor
8484 irvMaps[pI] = (labelField(identity(nCells)) + nTotalCells);
8485 irsMaps[pI] = (labelField(identity(nIntFaces)) + nTotalIntFaces);
8488 nTotalCells += nCells;
8489 nTotalIntFaces += nIntFaces;
8491 // Fetch patch sizes from subMesh
8492 const labelList& nPatchFaces =
8494 cMap.entityBuffer(coupleMap::FACE_SIZES)
8497 // Loop over physical patches
8498 forAll(nTotalPatchFaces, patchI)
8500 // Fetch patch size from subMesh
8501 label nFaces = nPatchFaces[patchI];
8503 // Set rmap for this patch
8504 brMaps[pI][patchI] =
8506 labelField(identity(nFaces))
8507 + nTotalPatchFaces[patchI]
8510 // Update patch-face count
8511 nTotalPatchFaces[patchI] += nFaces;
8515 // Prepare internal mappers
8516 labelList cellAddressing(nTotalCells, 0);
8517 labelList faceAddressing(nTotalIntFaces, 0);
8519 // Set identity map for first nCells,
8520 // and map from cell[0] for the rest
8521 for (label i = 0; i < nOldCells_; i++)
8523 cellAddressing[i] = i;
8526 // Set identity map for first nIntFaces,
8527 // and map from face[0] for the rest
8528 for (label i = 0; i < nOldInternalFaces_; i++)
8530 faceAddressing[i] = i;
8533 coupledMesh::subMeshMapper vMap
8539 coupledMesh::subMeshMapper sMap
8545 // Prepare boundary mappers
8546 labelListList patchAddressing(nPhysical);
8547 PtrList<coupledMesh::subMeshMapper> bMap(nPhysical);
8549 forAll(bMap, patchI)
8551 // Prepare patch mappers
8552 patchAddressing[patchI].setSize(nTotalPatchFaces[patchI], 0);
8554 // Set identity map for first nPatchFaces,
8555 // and map from patch-face[0] for the rest
8556 for (label i = 0; i < oldPatchSizes_[patchI]; i++)
8558 patchAddressing[patchI][i] = i;
8561 // Set the boundary mapper pointer
8565 new coupledMesh::subMeshMapper
8567 oldPatchSizes_[patchI],
8568 patchAddressing[patchI]
8573 // Loop through all volFields and re-size
8574 // to accomodate additional cells / faces
8575 coupledMesh::resizeMap(names[0], *this, vMap, irvMaps, bMap, brMaps, vsF);
8576 coupledMesh::resizeMap(names[1], *this, vMap, irvMaps, bMap, brMaps, vvF);
8577 coupledMesh::resizeMap(names[2], *this, vMap, irvMaps, bMap, brMaps, vsptF);
8578 coupledMesh::resizeMap(names[3], *this, vMap, irvMaps, bMap, brMaps, vsytF);
8579 coupledMesh::resizeMap(names[4], *this, vMap, irvMaps, bMap, brMaps, vtF);
8581 coupledMesh::resizeMap(names[5], *this, sMap, irsMaps, bMap, brMaps, ssF);
8582 coupledMesh::resizeMap(names[6], *this, sMap, irsMaps, bMap, brMaps, svF);
8583 coupledMesh::resizeMap(names[7], *this, sMap, irsMaps, bMap, brMaps, ssptF);
8584 coupledMesh::resizeMap(names[8], *this, sMap, irsMaps, bMap, brMaps, ssytF);
8585 coupledMesh::resizeMap(names[9], *this, sMap, irsMaps, bMap, brMaps, stF);
8589 // Initialize coupled boundary ordering
8590 // - Assumes that faces_ and points_ are consistent
8591 // - Assumes that patchStarts_ and patchSizes_ are consistent
8592 void dynamicTopoFvMesh::initCoupledBoundaryOrdering
8594 List<pointField>& centres,
8595 List<pointField>& anchors
8598 if (!Pstream::parRun())
8603 for (label pI = 0; pI < nPatches_; pI++)
8605 label neiProcNo = getNeighbourProcessor(pI);
8607 if (neiProcNo == -1)
8612 // Check if this is a master processor patch.
8613 label start = patchStarts_[pI];
8614 label size = patchSizes_[pI];
8616 // Prepare centres and anchors
8617 centres[pI].setSize(size, vector::zero);
8618 anchors[pI].setSize(size, vector::zero);
8620 if (Pstream::myProcNo() < neiProcNo)
8622 forAll(centres[pI], fI)
8624 centres[pI][fI] = faces_[fI + start].centre(points_);
8625 anchors[pI][fI] = points_[faces_[fI + start][0]];
8631 Pout<< "Sending to: " << neiProcNo
8632 << " nCentres: " << size << endl;
8634 // Write out my centres to disk
8639 + Foam::name(Pstream::myProcNo())
8641 + Foam::name(neiProcNo),
8647 // Ensure that we're sending the right size
8648 meshOps::pWrite(neiProcNo, size);
8651 // Send information to neighbour
8652 meshOps::pWrite(neiProcNo, centres[pI]);
8653 meshOps::pWrite(neiProcNo, anchors[pI]);
8658 label nEntities = -1;
8660 // Ensure that we're receiving the right size
8661 meshOps::pRead(neiProcNo, nEntities);
8665 Pout<< " Recving from: " << neiProcNo
8666 << " nCentres: " << size << endl;
8669 if (nEntities != size)
8671 // Set the size to what we're receiving
8672 centres[pI].setSize(nEntities, vector::zero);
8673 anchors[pI].setSize(nEntities, vector::zero);
8675 // Issue a warning now, and let
8676 // syncCoupledBoundaryOrdering fail later
8679 "void dynamicTopoFvMesh::"
8680 "initCoupledBoundaryOrdering() const"
8682 << "Incorrect send / recv sizes: " << nl
8683 << " nEntities: " << nEntities << nl
8684 << " size: " << size << nl
8689 // Schedule receive from neighbour
8690 meshOps::pRead(neiProcNo, centres[pI]);
8691 meshOps::pRead(neiProcNo, anchors[pI]);
8697 // Synchronize coupled boundary ordering
8698 bool dynamicTopoFvMesh::syncCoupledBoundaryOrdering
8700 List<pointField>& centres,
8701 List<pointField>& anchors,
8702 labelListList& patchMaps,
8703 labelListList& rotations
8706 bool anyChange = false, failedPatchMatch = false;
8708 // Calculate centres and tolerances for any slave patches
8709 List<scalarField> slaveTols(nPatches_);
8710 List<pointField> slaveCentres(nPatches_);
8712 // Handle locally coupled patches
8713 forAll(patchCoupling_, pI)
8715 if (patchCoupling_(pI))
8717 const coupleMap& cMap = patchCoupling_[pI].map();
8719 label masterPatch = cMap.masterIndex();
8720 label slavePatch = cMap.slaveIndex();
8722 label mSize = patchSizes_[masterPatch];
8723 label sSize = patchSizes_[slavePatch];
8725 label mStart = patchStarts_[masterPatch];
8726 label sStart = patchStarts_[slavePatch];
8730 Pout<< " Patch size mismatch: " << nl
8731 << " mSize: " << mSize << nl
8732 << " sSize: " << sSize << nl
8733 << " mStart: " << mStart << nl
8734 << " sStart: " << sStart << nl
8735 << " master: " << masterPatch << nl
8736 << " slave: " << slavePatch << nl
8737 << abort(FatalError);
8740 // Calculate centres / tolerances
8741 anchors[masterPatch].setSize(mSize, vector::zero);
8742 centres[masterPatch].setSize(mSize, vector::zero);
8744 slaveTols[slavePatch].setSize(sSize, 0.0);
8745 slaveCentres[slavePatch].setSize(sSize, vector::zero);
8747 for (label fI = 0; fI < mSize; fI++)
8749 point& mfa = anchors[masterPatch][fI];
8750 point& mfc = centres[masterPatch][fI];
8751 point& sfc = slaveCentres[slavePatch][fI];
8753 const face& mFace = faces_[fI + mStart];
8754 const face& sFace = faces_[fI + sStart];
8756 // Calculate anchor / centres
8757 mfa = points_[mFace[0]];
8758 mfc = mFace.centre(points_);
8759 sfc = sFace.centre(points_);
8761 scalar maxLen = -GREAT;
8765 maxLen = max(maxLen, mag(points_[sFace[fpI]] - sfc));
8768 slaveTols[slavePatch][fI] = geomMatchTol_()*maxLen;
8771 // For cyclics, additionally test for halves,
8772 // and transform points as necessary.
8775 if (masterPatch == slavePatch)
8777 scalar featureCos = 0.9;
8778 label halfSize = (mSize / 2);
8779 label half0Index = 0, half1Index = 0;
8781 // Size up halfMap and slaveAnchors
8782 halfMap.setSize(mSize, -1);
8784 scalarField half1Tols(halfSize, 0.0);
8785 vectorField half1Anchors(halfSize, vector::zero);
8787 // Fetch face-normal for the first face
8788 vector fhN = boundaryMesh()[masterPatch].faceAreas()[0];
8791 fhN /= mag(fhN) + VSMALL;
8793 // Fetch transform type
8794 const cyclicPolyPatch& cyclicPatch =
8796 refCast<const cyclicPolyPatch>
8798 boundaryMesh()[masterPatch]
8804 cyclicPatch.transform()
8805 == cyclicPolyPatch::TRANSLATIONAL
8808 for (label fI = 0; fI < mSize; fI++)
8810 const face& mFace = faces_[fI + mStart];
8813 vector fN = mFace.normal(points_);
8816 fN /= mag(fN) + VSMALL;
8818 // Check which half this face belongs to...
8819 bool isFirstHalf = ((fhN & fN) > featureCos);
8823 vector mA = anchors[masterPatch][fI];
8824 vector mC = centres[masterPatch][fI];
8826 // Set master anchor
8827 anchors[masterPatch][half0Index] = mA;
8829 // Transform master points
8832 mA += cyclicPatch.separationVector();
8833 mC += cyclicPatch.separationVector();
8837 // Assume constant transform tensor
8838 mA = Foam::transform(cyclicPatch.forwardT()[0], mA);
8839 mC = Foam::transform(cyclicPatch.forwardT()[0], mC);
8842 // Set the transformed anchor / centre
8843 half1Anchors[half0Index] = mA;
8844 centres[masterPatch][half0Index] = mC;
8845 half1Tols[half0Index] = slaveTols[masterPatch][fI];
8848 halfMap[fI] = half0Index;
8854 // Set the slave centre
8855 const scalar& sT = slaveTols[slavePatch][fI];
8856 const point& sC = slaveCentres[slavePatch][fI];
8858 // Duplicate slaveTols for first / second half
8859 slaveTols[slavePatch][half1Index] = sT;
8860 slaveCentres[slavePatch][half1Index] = sC;
8863 halfMap[fI] = half1Index + halfSize;
8869 // Sanity check for halfSize
8870 if (half0Index != halfSize || half1Index != halfSize)
8872 Pout<< " Master: " << masterPatch << nl
8873 << " Slave: " << slavePatch << nl
8874 << " half0Index: " << half0Index << nl
8875 << " half1Index: " << half1Index << nl
8876 << " halfSize: " << halfSize << nl
8877 << " failed to divide into halves."
8878 << abort(FatalError);
8882 centres[masterPatch].setSize(halfSize);
8883 slaveCentres[slavePatch].setSize(halfSize);
8885 // Set slave tols / anchors
8886 forAll(half1Anchors, fI)
8888 slaveTols[masterPatch][fI + halfSize] = half1Tols[fI];
8889 anchors[masterPatch][fI + halfSize] = half1Anchors[fI];
8893 // Fetch reference to slave map
8894 labelList& patchMap = patchMaps[slavePatch];
8896 // Try zero separation automatic matching
8903 slaveCentres[slavePatch],
8904 centres[masterPatch],
8905 slaveTols[slavePatch],
8911 // Write out master centres to disk
8912 if (debug > 3 || !matchedAll)
8917 "localMasterCentres_"
8918 + Foam::name(masterPatch)
8920 + Foam::name(slavePatch),
8921 mSize, mSize, mSize,
8922 centres[masterPatch]
8928 "localSlaveCentres_"
8929 + Foam::name(masterPatch)
8931 + Foam::name(slavePatch),
8932 sSize, sSize, sSize,
8933 slaveCentres[slavePatch]
8939 + Foam::name(masterPatch)
8941 + Foam::name(slavePatch),
8942 identity(mSize) + mStart,
8949 + Foam::name(masterPatch)
8951 + Foam::name(slavePatch),
8952 identity(sSize) + sStart,
8958 Pout<< " Master: " << masterPatch
8959 << " Slave: " << slavePatch
8960 << " mSize: " << mSize << " sSize: " << sSize
8961 << " failed on match for face centres."
8964 // Failed on match-for-all, so patchMaps
8965 // will be invalid. Bail out for now.
8966 failedPatchMatch = true;
8972 // Renumber patchMap for cyclics
8973 if (masterPatch == slavePatch)
8975 label halfSize = (mSize / 2);
8979 if (halfMap[fI] >= halfSize)
8983 patchMap[halfMap[fI] - halfSize]
8990 patchMap.transfer(halfMap);
8993 // Fetch reference to rotation map
8994 labelList& rotation = rotations[slavePatch];
8996 // Initialize rotation to zero
8997 rotation.setSize(sSize, 0);
9000 forAll(patchMap, oldFaceI)
9002 label newFaceI = patchMap[oldFaceI];
9004 const point& anchor = anchors[slavePatch][newFaceI];
9005 const scalar& faceTol = slaveTols[slavePatch][oldFaceI];
9006 const face& checkFace = faces_[sStart + oldFaceI];
9008 label anchorFp = -1;
9009 scalar minDSqr = GREAT;
9011 forAll(checkFace, fpI)
9013 scalar dSqr = magSqr(anchor - points_[checkFace[fpI]]);
9022 if (anchorFp == -1 || mag(minDSqr) > faceTol)
9027 "void dynamicTopoFvMesh::"
9028 "syncCoupledBoundaryOrdering\n"
9030 " List<pointField>& centres,\n"
9031 " List<pointField>& anchors,\n"
9032 " labelListList& patchMaps,\n"
9033 " labelListList& rotations\n"
9036 << "Cannot find anchor: " << anchor << nl
9037 << " Face: " << checkFace << nl
9039 << UIndirectList<point>(points_, checkFace) << nl
9040 << " on patch: " << slavePatch << nl
9041 << " faceTol: " << faceTol << nl
9042 << " newFaceI: " << newFaceI << nl
9043 << " oldFaceI: " << oldFaceI << nl
9044 << " mSize: " << mSize << nl
9045 << " sSize: " << sSize << nl
9046 << abort(FatalError);
9050 // Positive rotation
9051 // - Set for old face. Will be rotated later
9052 // during the shuffling stage
9053 rotation[oldFaceI] =
9055 (checkFace.size() - anchorFp) % checkFace.size()
9065 if (failedPatchMatch)
9070 "void dynamicTopoFvMesh::"
9071 "syncCoupledBoundaryOrdering\n"
9073 " List<pointField>& centres,\n"
9074 " List<pointField>& anchors,\n"
9075 " labelListList& patchMaps,\n"
9076 " labelListList& rotations\n"
9079 << " Matching for local patches failed. " << nl
9080 << abort(FatalError);
9083 if (!Pstream::parRun())
9088 for (label pI = 0; pI < nPatches_; pI++)
9090 label neiProcNo = getNeighbourProcessor(pI);
9092 if (neiProcNo == -1)
9097 // Check if this is a slave processor patch.
9098 label start = patchStarts_[pI];
9099 label size = patchSizes_[pI];
9101 if (Pstream::myProcNo() > neiProcNo)
9103 slaveTols[pI].setSize(size, 0.0);
9104 slaveCentres[pI].setSize(size, vector::zero);
9106 forAll(slaveCentres[pI], fI)
9108 point& fc = slaveCentres[pI][fI];
9110 const face& checkFace = faces_[fI + start];
9113 fc = checkFace.centre(points_);
9115 scalar maxLen = -GREAT;
9117 forAll(checkFace, fpI)
9119 maxLen = max(maxLen, mag(points_[checkFace[fpI]] - fc));
9122 slaveTols[pI][fI] = geomMatchTol_()*maxLen;
9125 // Write out my centres to disk
9132 + Foam::name(Pstream::myProcNo())
9134 + Foam::name(neiProcNo),
9142 // Wait for transfers before continuing.
9143 meshOps::waitForBuffers();
9145 for (label pI = 0; pI < nPatches_; pI++)
9147 label neiProcNo = getNeighbourProcessor(pI);
9149 if (neiProcNo == -1)
9154 // Check if this is a master processor patch.
9155 labelList& patchMap = patchMaps[pI];
9156 labelList& rotation = rotations[pI];
9158 // Initialize map and rotation
9159 patchMap.setSize(patchSizes_[pI], -1);
9160 rotation.setSize(patchSizes_[pI], 0);
9162 if (Pstream::myProcNo() < neiProcNo)
9164 // Do nothing (i.e. identical mapping, zero rotation).
9165 forAll(patchMap, pfI)
9167 patchMap[pfI] = pfI;
9172 // Try zero separation automatic matching
9173 label mSize = centres[pI].size();
9174 label sSize = slaveCentres[pI].size();
9190 // Write out centres to disk
9191 if (debug > 3 || !matchedAll)
9197 + Foam::name(Pstream::myProcNo())
9199 + Foam::name(neiProcNo),
9200 mSize, mSize, mSize,
9208 + Foam::name(Pstream::myProcNo())
9210 + Foam::name(neiProcNo),
9211 sSize, sSize, sSize,
9218 + Foam::name(Pstream::myProcNo())
9220 + Foam::name(neiProcNo),
9221 identity(patchSizes_[pI]) + patchStarts_[pI],
9227 Pout<< " Patch: " << pI
9228 << " Processor: " << neiProcNo
9229 << " mSize: " << mSize << " sSize: " << sSize
9230 << " failed on match for face centres."
9233 // Failed on match-for-all, so patchMaps
9234 // will be invalid. Bail out for now.
9235 failedPatchMatch = true;
9241 label start = patchStarts_[pI];
9244 forAll(patchMap, oldFaceI)
9246 label newFaceI = patchMap[oldFaceI];
9248 const point& anchor = anchors[pI][newFaceI];
9249 const scalar& faceTol = slaveTols[pI][oldFaceI];
9250 const face& checkFace = faces_[start + oldFaceI];
9252 label anchorFp = -1;
9253 scalar minDSqr = GREAT;
9255 forAll(checkFace, fpI)
9257 scalar dSqr = magSqr(anchor - points_[checkFace[fpI]]);
9266 if (anchorFp == -1 || mag(minDSqr) > faceTol)
9271 "void dynamicTopoFvMesh::"
9272 "syncCoupledBoundaryOrdering\n"
9274 " List<pointField>& centres,\n"
9275 " List<pointField>& anchors,\n"
9276 " labelListList& patchMaps,\n"
9277 " labelListList& rotations\n"
9280 << "Cannot find anchor: " << anchor << nl
9281 << " Face: " << checkFace << nl
9283 << UIndirectList<point>(points_, checkFace) << nl
9284 << " on patch: " << pI
9285 << " for processor: " << neiProcNo
9286 << abort(FatalError);
9290 // Positive rotation
9291 // - Set for old face. Will be rotated later
9292 // during the shuffling stage
9293 rotation[oldFaceI] =
9295 (checkFace.size() - anchorFp) % checkFace.size()
9305 // Reduce failures across processors
9306 reduce(failedPatchMatch, orOp<bool>());
9308 // Write out processor patch faces if a failure was encountered
9309 if (failedPatchMatch)
9311 for (label pI = 0; pI < nPatches_; pI++)
9313 label neiProcNo = getNeighbourProcessor(pI);
9315 if (neiProcNo == -1)
9323 + Foam::name(Pstream::myProcNo())
9325 + Foam::name(neiProcNo),
9326 identity(patchSizes_[pI]) + patchStarts_[pI],
9334 "void dynamicTopoFvMesh::"
9335 "syncCoupledBoundaryOrdering\n"
9337 " List<pointField>& centres,\n"
9338 " List<pointField>& anchors,\n"
9339 " labelListList& patchMaps,\n"
9340 " labelListList& rotations\n"
9343 << " Matching for processor patches failed. " << nl
9344 << " Patch faces written out to disk." << nl
9345 << abort(FatalError);
9352 // Fill buffers with length-scale info
9353 // and exchange across processors.
9354 void dynamicTopoFvMesh::exchangeLengthBuffers()
9356 if (!Pstream::parRun())
9361 if (!edgeRefinement_)
9366 forAll(procIndices_, pI)
9368 coupledMesh& sPM = sendMeshes_[pI];
9369 coupledMesh& rPM = recvMeshes_[pI];
9371 const coupleMap& scMap = sPM.map();
9372 const coupleMap& rcMap = rPM.map();
9374 if (scMap.slaveIndex() == Pstream::myProcNo())
9376 // Clear existing buffer
9377 sPM.subMesh().lengthScale_.clear();
9379 // Fill in buffers to send.
9380 sPM.subMesh().lengthScale_.setSize
9382 scMap.nEntities(coupleMap::CELL),
9386 Map<label>& cellMap = scMap.entityMap(coupleMap::CELL);
9388 forAllIter(Map<label>, cellMap, cIter)
9390 sPM.subMesh().lengthScale_[cIter.key()] = lengthScale_[cIter()];
9395 scMap.masterIndex(),
9396 sPM.subMesh().lengthScale_
9401 Pout<< "Sending to: "
9402 << scMap.masterIndex()
9404 << scMap.nEntities(coupleMap::CELL)
9409 if (rcMap.masterIndex() == Pstream::myProcNo())
9411 // Clear existing buffer
9412 rPM.subMesh().lengthScale_.clear();
9414 // Schedule receipt from neighbour
9415 rPM.subMesh().lengthScale_.setSize
9417 rcMap.nEntities(coupleMap::CELL),
9421 // Schedule for receipt
9425 rPM.subMesh().lengthScale_
9430 Pout<< "Receiving from: "
9431 << rcMap.slaveIndex()
9433 << rcMap.nEntities(coupleMap::CELL)
9439 // Wait for transfers before continuing.
9440 meshOps::waitForBuffers();
9444 forAll(procIndices_, pI)
9446 const coupledMesh& rPM = recvMeshes_[pI];
9447 const coupleMap& rcMap = rPM.map();
9449 if (rcMap.masterIndex() == Pstream::myProcNo())
9451 rPM.subMesh().writeVTK
9453 "lengthScale_" + Foam::name(rcMap.slaveIndex()),
9454 identity(rPM.subMesh().nCells()),
9458 rPM.subMesh().lengthScale_
9466 // Implementing the fillTables operation for coupled edges
9467 bool dynamicTopoFvMesh::coupledFillTables
9472 PtrList<scalarListList>& Q,
9473 PtrList<labelListList>& K,
9474 PtrList<labelListList>& triangulations
9477 bool success = false;
9479 if (locallyCoupledEntity(eIndex))
9481 // Fill tables for the slave edge.
9484 // Determine the slave index.
9485 forAll(patchCoupling_, patchI)
9487 if (patchCoupling_(patchI))
9489 const label edgeEnum = coupleMap::EDGE;
9490 const coupleMap& cMap = patchCoupling_[patchI].map();
9492 if ((sI = cMap.findSlave(edgeEnum, eIndex)) > -1)
9504 "bool dynamicTopoFvMesh::coupledFillTables\n"
9506 " const label eIndex,\n"
9507 " const scalar minQuality,\n"
9509 " PtrList<scalarListList>& Q,\n"
9510 " PtrList<labelListList>& K,\n"
9511 " PtrList<labelListList>& triangulations\n"
9514 << "Coupled maps were improperly specified." << nl
9515 << " Slave index not found for: " << nl
9516 << " Edge: " << eIndex << nl
9517 << abort(FatalError);
9520 // Turn off switch temporarily.
9521 unsetCoupledModification();
9523 labelList hullV(10);
9525 // Call fillTables for the slave edge.
9526 success = fillTables(sI, minQuality, m, hullV, Q, K, triangulations, 1);
9529 setCoupledModification();
9532 if (processorCoupledEntity(eIndex))
9534 const edge& checkEdge = edges_[eIndex];
9535 const labelList& eFaces = edgeFaces_[eIndex];
9540 // Need to build alternate addressing / point-list
9541 // for swaps across processors.
9542 DynamicList<point> parPts(10);
9543 dynamicLabelList parVtx(10);
9546 label nPoints = 0, nProcs = 0;
9547 label otherPoint = -1, nextPoint = -1;
9549 // Define a line for this edge
9552 points_[checkEdge.start()],
9553 points_[checkEdge.end()]
9556 // Define tangent-to-edge / edge-centre
9557 vector te = -lpr.vec(), xe = lpr.centre();
9559 // Normalize tangent
9560 te /= mag(te) + VSMALL;
9562 // Fill-in vertices for this processor
9563 forAll(eFaces, faceI)
9565 label fPatch = whichPatch(eFaces[faceI]);
9566 label neiProc = getNeighbourProcessor(fPatch);
9571 && priority(neiProc, lessOp<label>(), Pstream::myProcNo())
9574 // This edge should not be here
9575 Pout<< " Edge: " << eIndex
9576 << " is talking to processor: " << neiProc
9577 << abort(FatalError);
9580 // This face is either internal or on a physical boundary
9581 const face& checkFace = faces_[eFaces[faceI]];
9583 // Find the isolated point
9584 meshOps::findIsolatedPoint
9592 // Insert point and index
9593 parPts.append(points_[otherPoint]);
9594 parVtx.append(nPoints);
9596 // Physical patch: Is this an appropriate start face?
9597 // - If yes, swap with first index
9598 if (fPatch > -1 && neiProc == -1)
9602 if (nextPoint == checkEdge[0])
9604 Foam::Swap(parPts[0], parPts[nPoints]);
9608 // Increment point index
9612 // Now look through processors, and add their points
9613 forAll(procIndices_, pI)
9615 label proc = procIndices_[pI];
9617 // Fetch reference to subMesh
9618 const coupleMap& cMap = recvMeshes_[pI].map();
9619 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
9621 label sI = -1, sP = -1;
9623 if ((sI = cMap.findSlave(coupleMap::EDGE, eIndex)) == -1)
9628 // If this is a new edge, bail out for now.
9629 // This will be handled at the next time-step.
9630 if (sI >= mesh.nOldEdges_)
9635 const edge& slaveEdge = mesh.edges_[sI];
9636 const labelList& seFaces = mesh.edgeFaces_[sI];
9638 // Determine the point index that corresponds to checkEdge[0]
9641 cMap.findMaster(coupleMap::POINT, slaveEdge[0]),
9642 cMap.findMaster(coupleMap::POINT, slaveEdge[1])
9645 if (cE[0] == checkEdge[0])
9650 if (cE[1] == checkEdge[0])
9656 label meP = whichEdgePatch(eIndex);
9657 word mN(meP < 0 ? "Internal" : boundaryMesh()[meP].name());
9659 label seP = mesh.whichEdgePatch(sI);
9660 word sN(seP < 0 ? "Internal" : mesh.boundaryMesh()[seP].name());
9662 Pout<< " Can't find master point: " << checkEdge[0] << nl
9663 << " on master: " << eIndex << "::" << checkEdge
9664 << " Patch: " << mN << nl
9665 << " with slave: " << sI << "::" << slaveEdge
9666 << " Patch: " << sN << nl
9667 << " cE: " << cE << " on proc: " << proc << nl
9668 << abort(FatalError);
9671 forAll(seFaces, faceI)
9673 label slavePatch = mesh.whichPatch(seFaces[faceI]);
9675 // If this is talking to a lower-ranked processor,
9676 // skip the insertion step.
9677 label neiProc = mesh.getNeighbourProcessor(slavePatch);
9682 && priority(neiProc, lessOp<label>(), proc)
9688 // This face is either internal or on a physical boundary
9689 const face& slaveFace = mesh.faces_[seFaces[faceI]];
9691 // Find the isolated point
9692 meshOps::findIsolatedPoint
9700 // Insert point and index
9701 parPts.append(mesh.points_[otherPoint]);
9702 parVtx.append(nPoints);
9704 // Physical patch: Is this an appropriate start face?
9705 // - If yes, swap with first index
9706 if (slavePatch > -1 && neiProc == -1)
9710 if (nextPoint == sP)
9712 Foam::Swap(parPts[0], parPts[nPoints]);
9716 // Increment point index
9723 // Sort points / indices in counter-clockwise order
9724 SortableList<scalar> angles(nPoints, 0.0);
9727 scalar twoPi = mathematicalConstant::twoPi;
9729 // Define a base direction
9730 // from the start point
9731 vector dir = (parPts[0] - xe);
9732 dir -= ((dir & te) * te);
9733 dir /= mag(dir) + VSMALL;
9735 // Local coordinate system
9736 coordinateSystem cs("cs", xe, te, dir);
9738 for (label i = 1; i < nPoints; i++)
9740 // Convert to local csys and determine angle
9741 vector local = cs.localPosition(parPts[i]);
9742 scalar angle = atan2(local.y(), local.x());
9744 // Account for 3rd and 4th quadrants
9745 angles[i] = (angle < 0.0 ? angle + twoPi : angle);
9751 // Reorder points and transfer
9752 List<point> sortedParPts(nPoints);
9754 const labelList& indices = angles.indices();
9756 forAll(sortedParPts, pointI)
9758 sortedParPts[pointI] = parPts[indices[pointI]];
9761 parPts.transfer(sortedParPts);
9763 // Fill the last two points for the edge
9764 edge parEdge(-1, -1);
9766 parPts.append(lpr.start());
9767 parEdge[0] = nPoints++;
9769 parPts.append(lpr.end());
9770 parEdge[1] = nPoints++;
9772 // Compute minQuality with this loop
9773 minQuality = computeMinQuality(parEdge, parVtx, parPts, closed);
9775 if (debug > 4 || minQuality < 0.0)
9777 // Write out edge connectivity
9778 writeEdgeConnectivity(eIndex);
9783 "parPts_" + Foam::name(eIndex),
9790 if (minQuality < 0.0)
9792 Switch isClosed(closed);
9794 Pout<< " * * * Error in fillTables * * * " << nl
9795 << " Edge: " << eIndex << " :: " << checkEdge << nl
9796 << " minQuality: " << minQuality << nl
9797 << " Closed: " << isClosed << nl
9798 << abort(FatalError);
9803 m[0] = parVtx.size();
9805 // Check if a table-resize is necessary
9806 if (m[0] > maxTetsPerEdge_)
9808 if (allowTableResize_)
9810 // Resize the tables to account for
9811 // more tets per edge
9812 label& mtpe = const_cast<label&>(maxTetsPerEdge_);
9816 // Clear tables for this index.
9819 triangulations[0].clear();
9821 // Resize for this index.
9822 initTables(m, Q, K, triangulations);
9826 // Can't resize. Bail out.
9831 // Fill dynamic programming tables
9851 // Initialize coupled weights calculation
9852 void dynamicTopoFvMesh::initCoupledWeights()
9854 if (!Pstream::parRun())
9859 forAll(procIndices_, pI)
9861 // Fetch reference to subMesh
9862 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
9868 const polyBoundaryMesh& sBoundary = mesh.boundaryMesh();
9870 forAll(sBoundary, patchI)
9872 sBoundary[patchI].faceFaces();
9878 // Additional mapping contributions for coupled entities
9879 void dynamicTopoFvMesh::computeCoupledWeights
9882 const label dimension,
9884 scalarField& weights,
9885 vectorField& centres,
9889 if (!Pstream::parRun())
9894 // Fetch offsets from mapper
9895 const labelList& cStarts = mapper_->cellStarts();
9896 const labelListList& pSizes = mapper_->patchSizes();
9900 dynamicLabelList faceParents(10);
9902 label patchIndex = whichPatch(index);
9904 forAll(procIndices_, pI)
9906 // Ensure that the patch is physical
9907 if (patchIndex < 0 || patchIndex >= pSizes[pI].size())
9909 Pout<< " Face: " << index
9910 << " Patch: " << patchIndex
9911 << " does not belong to a physical patch." << nl
9912 << " nPhysicalPatches: " << pSizes[pI].size()
9913 << abort(FatalError);
9916 // Fetch reference to subMesh
9917 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
9920 labelList coupleObjects;
9921 scalarField coupleWeights;
9922 vectorField coupleCentres;
9924 // Convex-set algorithm for faces
9925 faceSetAlgorithm faceAlgorithm
9936 // Initialize the bounding box
9937 faceAlgorithm.computeNormFactor(index);
9939 // Loop through all subMesh faces, and check for bounds
9940 const polyBoundaryMesh& boundary = mesh.boundaryMesh();
9942 label pSize = boundary[patchIndex].size();
9943 label pStart = boundary[patchIndex].start();
9945 for (label faceI = pStart; faceI < (pStart + pSize); faceI++)
9947 if (faceAlgorithm.contains(faceI))
9949 faceParents.append(faceI);
9953 // Obtain weighting factors for this face.
9954 faceAlgorithm.computeWeights
9959 boundary[patchIndex].faceFaces(),
9966 // Add contributions with offsets
9967 if (coupleObjects.size())
9969 // Fetch patch size on master
9970 label patchSize = boundaryMesh()[patchIndex].size();
9973 label oldSize = parents.size();
9975 parents.setSize(oldSize + coupleObjects.size());
9976 weights.setSize(oldSize + coupleObjects.size());
9977 centres.setSize(oldSize + coupleObjects.size());
9979 forAll(coupleObjects, indexI)
9981 parents[indexI + oldSize] =
9983 patchSize + coupleObjects[indexI]
9986 weights[indexI + oldSize] = coupleWeights[indexI];
9987 centres[indexI + oldSize] = coupleCentres[indexI];
9992 faceParents.clear();
9998 dynamicLabelList cellParents(10);
10000 forAll(procIndices_, pI)
10002 // Fetch reference to subMesh
10003 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
10006 labelList coupleObjects;
10007 scalarField coupleWeights;
10008 vectorField coupleCentres;
10010 // Convex-set algorithm for cells
10011 cellSetAlgorithm cellAlgorithm
10022 // Initialize the bounding box
10023 cellAlgorithm.computeNormFactor(index);
10025 // Loop through all subMesh cells, and check for bounds
10026 const cellList& meshCells = mesh.cells();
10028 forAll(meshCells, cellI)
10030 if (cellAlgorithm.contains(cellI))
10032 cellParents.append(cellI);
10036 // Obtain weighting factors for this cell.
10037 cellAlgorithm.computeWeights
10042 mesh.polyMesh::cellCells(),
10049 // Add contributions with offsets
10050 if (coupleObjects.size())
10052 label cellStart = cStarts[pI];
10055 label oldSize = parents.size();
10057 parents.setSize(oldSize + coupleObjects.size());
10058 weights.setSize(oldSize + coupleObjects.size());
10059 centres.setSize(oldSize + coupleObjects.size());
10061 forAll(coupleObjects, indexI)
10063 parents[indexI + oldSize] =
10065 cellStart + coupleObjects[indexI]
10068 weights[indexI + oldSize] = coupleWeights[indexI];
10069 centres[indexI + oldSize] = coupleCentres[indexI];
10074 cellParents.clear();
10079 FatalErrorIn("scalar dynamicTopoFvMesh::computeCoupledWeights()")
10080 << " Incorrect dimension: " << dimension << nl
10081 << abort(FatalError);
10086 // Fetch length-scale info for processor entities
10087 scalar dynamicTopoFvMesh::processorLengthScale(const label index) const
10089 scalar procScale = 0.0;
10093 // First check the master processor
10094 procScale += lengthScale_[owner_[index]];
10096 // Next, check the slave processor
10097 bool foundSlave = false;
10099 forAll(procIndices_, pI)
10101 // Fetch non-const reference to subMeshes
10102 const label faceEnum = coupleMap::FACE;
10103 const coupleMap& cMap = recvMeshes_[pI].map();
10104 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
10108 if ((sIndex = cMap.findSlave(faceEnum, index)) > -1)
10110 procScale += mesh.lengthScale_[mesh.owner_[sIndex]];
10117 // Should have found at least one slave
10122 "scalar dynamicTopoFvMesh::processorLengthScale"
10123 "(const label index) const"
10125 << "Processor lengthScale lookup failed: " << nl
10126 << " Master face: " << index
10127 << " :: " << faces_[index] << nl
10128 << abort(FatalError);
10131 // Average the scale
10136 // Check if this is a 'pure' processor edge
10137 bool pure = processorCoupledEntity(index, false, true, true);
10139 const label edgeEnum = coupleMap::EDGE;
10140 const labelList& eFaces = edgeFaces_[index];
10142 bool foundSlave = false;
10146 // First check the master processor
10149 forAll(eFaces, faceI)
10151 label own = owner_[eFaces[faceI]];
10152 label nei = neighbour_[eFaces[faceI]];
10154 procScale += lengthScale_[own];
10159 procScale += lengthScale_[nei];
10164 // Next check slaves
10165 forAll(procIndices_, pI)
10167 const coupleMap& cMap = recvMeshes_[pI].map();
10168 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
10172 if ((sIndex = cMap.findSlave(edgeEnum, index)) > -1)
10174 // Fetch connectivity from patchSubMesh
10175 const labelList& peFaces = mesh.edgeFaces_[sIndex];
10179 forAll(peFaces, faceI)
10181 label own = mesh.owner_[peFaces[faceI]];
10182 label nei = mesh.neighbour_[peFaces[faceI]];
10184 procScale += mesh.lengthScale_[own];
10189 procScale += mesh.lengthScale_[nei];
10196 // Average the final scale
10201 // Processor is adjacent to physical patch types.
10202 // Search for boundary faces, and average their scale
10204 // First check the master processor
10205 label nBoundary = 0;
10207 forAll(eFaces, faceI)
10209 label facePatch = whichPatch(eFaces[faceI]);
10211 // Skip internal faces
10212 if (facePatch == -1)
10217 // Skip processor patches
10218 if (getNeighbourProcessor(facePatch) > -1)
10223 // If this is a floating face, pick the owner length-scale
10224 if (lengthEstimator().isFreePatch(facePatch))
10226 procScale += lengthScale_[owner_[eFaces[faceI]]];
10230 // Fetch fixed length-scale
10233 lengthEstimator().fixedLengthScale
10244 forAll(procIndices_, pI)
10246 const coupleMap& cMap = recvMeshes_[pI].map();
10247 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
10251 if ((sIndex = cMap.findSlave(edgeEnum, index)) > -1)
10253 // Fetch connectivity from patchSubMesh
10254 const labelList& peFaces = mesh.edgeFaces_[sIndex];
10258 forAll(peFaces, faceI)
10260 label facePatch = mesh.whichPatch(peFaces[faceI]);
10262 // Skip internal faces
10263 if (facePatch == -1)
10268 // Skip processor patches
10269 if (mesh.getNeighbourProcessor(facePatch) > -1)
10274 // If this is a floating face,
10275 // pick the owner length-scale
10276 if (lengthEstimator().isFreePatch(facePatch))
10280 mesh.lengthScale_[mesh.owner_[peFaces[faceI]]]
10285 // Fetch fixed length-scale
10288 lengthEstimator().fixedLengthScale
10301 if (nBoundary != 2)
10303 // Write out for post-processing
10304 writeEdgeConnectivity(index);
10308 "scalar dynamicTopoFvMesh::processorLengthScale"
10309 "(const label index) const"
10311 << " Expected two physical boundary patches: " << nl
10312 << " nBoundary: " << nBoundary
10313 << " Master edge: " << index
10314 << " :: " << edges_[index]
10316 << boundaryMesh()[whichEdgePatch(index)].name() << nl
10317 << abort(FatalError);
10323 // Should have found at least one slave
10328 "scalar dynamicTopoFvMesh::processorLengthScale"
10329 "(const label index) const"
10331 << "Processor lengthScale lookup failed: " << nl
10332 << " Master edge: " << index
10333 << " :: " << edges_[index] << nl
10334 << abort(FatalError);
10342 // Method to determine whether the master face is locally coupled
10343 bool dynamicTopoFvMesh::locallyCoupledEntity
10351 // Bail out if no patchCoupling is present
10352 if (patchCoupling_.empty())
10357 if (is2D() || checkFace)
10359 label patch = whichPatch(index);
10366 // Processor checks receive priority.
10369 if (getNeighbourProcessor(patch) > -1)
10375 // Check coupled master patches.
10376 if (patchCoupling_(patch))
10383 // Check on slave patches as well.
10384 forAll(patchCoupling_, pI)
10386 if (patchCoupling_(pI))
10388 const coupleMap& cMap = patchCoupling_[pI].map();
10390 if (cMap.slaveIndex() == patch)
10400 const labelList& eFaces = edgeFaces_[index];
10402 // Search for boundary faces, and determine boundary type.
10403 forAll(eFaces, faceI)
10405 if (neighbour_[eFaces[faceI]] == -1)
10407 label patch = whichPatch(eFaces[faceI]);
10409 // Processor checks receive priority.
10412 if (getNeighbourProcessor(patch) > -1)
10418 // Check coupled master patches.
10419 if (patchCoupling_(patch))
10426 // Check on slave patches as well.
10427 forAll(patchCoupling_, pI)
10429 if (patchCoupling_(pI))
10431 const coupleMap& cMap =
10433 patchCoupling_[pI].map()
10436 if (cMap.slaveIndex() == patch)
10447 // Could not find any faces on locally coupled patches.
10452 // Method to determine the locally coupled patch index
10453 label dynamicTopoFvMesh::locallyCoupledEdgePatch(const label eIndex) const
10455 const labelList& eFaces = edgeFaces_[eIndex];
10457 // Search for boundary faces, and determine boundary type.
10458 forAll(eFaces, faceI)
10460 if (neighbour_[eFaces[faceI]] == -1)
10462 label patch = whichPatch(eFaces[faceI]);
10464 // Check coupled master patches.
10465 if (patchCoupling_(patch))
10470 // Check on slave patches as well.
10471 forAll(patchCoupling_, pI)
10473 if (patchCoupling_(pI))
10475 const coupleMap& cMap = patchCoupling_[pI].map();
10477 if (cMap.slaveIndex() == patch)
10486 // Could not find any faces on locally coupled patches.
10489 "label dynamicTopoFvMesh::locallyCoupledEdgePatch"
10490 "(const label cIndex) const"
10492 << "Edge: " << eIndex << ":: " << edges_[eIndex]
10493 << " does not lie on any coupled patches."
10494 << abort(FatalError);
10500 // Method to determine if the entity is on a processor boundary
10501 // - Also provides an additional check for 'pure' processor edges
10502 // i.e., edges that do not abut a physical patch. This is necessary
10503 // while deciding on collapse cases towards bounding curves.
10504 bool dynamicTopoFvMesh::processorCoupledEntity
10510 FixedList<label, 2>* patchLabels,
10511 FixedList<vector, 2>* patchNormals
10514 // Skip check for serial runs
10515 if (!Pstream::parRun())
10522 if ((is2D() || checkFace) && !checkEdge)
10524 patch = whichPatch(index);
10531 if (getNeighbourProcessor(patch) > -1)
10538 const labelList& eFaces = edgeFaces_[index];
10540 label nPhysical = 0, nProcessor = 0;
10542 // Search for boundary faces, and determine boundary type.
10543 forAll(eFaces, faceI)
10545 label patch = whichPatch(eFaces[faceI]);
10552 if (getNeighbourProcessor(patch) > -1)
10554 // Increment the processor patch count
10559 // We don't have to validate if this
10560 // is a 'pure' processor edge, so bail out.
10569 (*patchLabels)[nPhysical] = patch;
10574 (*patchNormals)[nPhysical] =
10576 faces_[eFaces[faceI]].normal(points_)
10584 if (patchLabels && patchNormals)
10586 // Check other coupled-edges as well
10587 forAll(procIndices_, pI)
10589 // Fetch reference to subMesh
10590 const coupleMap& cMap = recvMeshes_[pI].map();
10591 const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
10595 if ((sI = cMap.findSlave(coupleMap::EDGE, index)) == -1)
10600 const labelList& seFaces = mesh.edgeFaces_[sI];
10602 forAll(seFaces, faceI)
10604 label sPatch = mesh.whichPatch(seFaces[faceI]);
10605 label neiProc = mesh.getNeighbourProcessor(sPatch);
10607 // Skip internal / processor faces
10608 if (sPatch == -1 || neiProc > -1)
10613 const face& sFace = mesh.faces_[seFaces[faceI]];
10615 // Fill patch / normal info
10616 (*patchLabels)[nPhysical] = sPatch;
10617 (*patchNormals)[nPhysical] = sFace.normal(mesh.points_);
10627 if (nProcessor && !nPhysical)
10634 // Could not find any faces on processor patches.
10639 // Build a list of entities that need to be avoided
10640 // by regular topo-changes.
10641 void dynamicTopoFvMesh::buildEntitiesToAvoid
10643 labelHashSet& entities,
10649 // Build a set of entities to avoid during regular modifications,
10650 // and build a master stack for coupled modifications.
10652 // Determine locally coupled slave patches.
10653 labelHashSet localMasterPatches, localSlavePatches;
10655 forAll(patchCoupling_, patchI)
10657 if (patchCoupling_(patchI))
10659 const coupleMap& cMap = patchCoupling_[patchI].map();
10661 localMasterPatches.insert(cMap.masterIndex());
10662 localSlavePatches.insert(cMap.slaveIndex());
10666 // Loop through boundary faces and check whether
10667 // they belong to master/slave coupled patches.
10668 for (label faceI = nOldInternalFaces_; faceI < faces_.size(); faceI++)
10670 // Add only valid faces
10671 if (faces_[faceI].empty())
10676 label pIndex = whichPatch(faceI);
10683 // Check if this is a coupled face.
10686 localMasterPatches.found(pIndex) ||
10687 localSlavePatches.found(pIndex) ||
10688 getNeighbourProcessor(pIndex) > -1
10693 // Avoid this face during regular modification.
10694 if (!entities.found(faceI))
10696 entities.insert(faceI);
10701 const labelList& fEdges = faceEdges_[faceI];
10703 forAll(fEdges, edgeI)
10705 // Avoid this edge during regular modification.
10706 if (!entities.found(fEdges[edgeI]))
10708 entities.insert(fEdges[edgeI]);
10715 // Loop through entities contained in patchSubMeshes, if requested
10718 forAll(procIndices_, pI)
10720 const coupleMap& cMap = sendMeshes_[pI].map();
10721 const Map<label> rEdgeMap = cMap.reverseEntityMap(coupleMap::EDGE);
10723 if (cMap.slaveIndex() == Pstream::myProcNo())
10725 forAllConstIter(Map<label>, rEdgeMap, eIter)
10729 const labelList& eFaces = edgeFaces_[eIter.key()];
10731 forAll(eFaces, faceI)
10733 if (faces_[eFaces[faceI]].size() == 4)
10735 if (!entities.found(eFaces[faceI]))
10737 entities.insert(eFaces[faceI]);
10744 const edge& check = edges_[eIter.key()];
10745 const labelList& eFaces = edgeFaces_[eIter.key()];
10747 // Skip deleted edges
10752 const labelList& pE = pointEdges_[check[pI]];
10756 if (!entities.found(pE[edgeI]))
10758 entities.insert(pE[edgeI]);
10771 Pout<< nl << "nEntitiesToAvoid: " << entities.size() << endl;
10775 // Write out entities
10776 label elemType = is2D() ? 2 : 1;
10781 + Foam::name(Pstream::myProcNo()),
10790 // Check whether the specified edge is a coupled master edge.
10791 bool dynamicTopoFvMesh::isCoupledMaster
10796 if (!coupledModification_)
10801 return locallyCoupledEntity(eIndex);
10805 } // End namespace Foam
10807 // ************************************************************************* //