43 const scalar faceCoupleInfo::angleTol_ = 1
e-3;
49 void Foam::faceCoupleInfo::writeOBJ
51 const fileName& fName,
69 const edge&
e = edges[edgeI];
73 const label pointi =
e[eI];
75 if (pointMap[pointi] == -1)
96 const edge&
e = edges[edgeI];
98 str<<
"l " << pointMap[
e[0]]+1 <<
' ' << pointMap[
e[1]]+1 <<
nl;
103 void Foam::faceCoupleInfo::writeOBJ
105 const fileName& fName,
110 Pout<<
"Writing connections as edges to " << fName <<
endl;
122 str <<
"l " << vertI-1 <<
' ' << vertI <<
nl;
127 void Foam::faceCoupleInfo::writePointsFaces()
const 135 OFstream str(
"masterPatch.obj");
140 OFstream str(
"slavePatch.obj");
145 OFstream str(
"cutFaces.obj");
152 Pout<<
"Writing cutToMasterPoints to cutToMasterPoints.obj" <<
endl;
156 "cutToMasterPoints.obj",
161 Pout<<
"Writing cutToSlavePoints to cutToSlavePoints.obj" <<
endl;
165 "cutToSlavePoints.obj",
173 Pout<<
"Writing cutToMasterFaces to cutToMasterFaces.obj" <<
endl;
177 forAll(cutToMasterFaces(), cutFacei)
179 label masterFacei = cutToMasterFaces()[cutFacei];
181 if (masterFacei != -1)
183 equivMasterFaces[cutFacei] = m[masterFacei].centre(m.points());
188 <<
"No master face for cut face " << cutFacei
189 <<
" at position " <<
c[cutFacei].centre(
c.points())
192 equivMasterFaces[cutFacei] =
Zero;
198 "cutToMasterFaces.obj",
199 calcFaceCentres<List>(
c, cutPoints(), 0,
c.size()),
205 Pout<<
"Writing cutToSlaveFaces to cutToSlaveFaces.obj" <<
endl;
209 forAll(cutToSlaveFaces(), cutFacei)
211 label slaveFacei = cutToSlaveFaces()[cutFacei];
213 equivSlaveFaces[cutFacei] =
s[slaveFacei].centre(
s.points());
218 "cutToSlaveFaces.obj",
219 calcFaceCentres<List>(
c, cutPoints(), 0,
c.size()),
228 void Foam::faceCoupleInfo::writeEdges
240 OFstream str(
"cutToMasterEdges.obj");
241 Pout<<
"Writing cutToMasterEdges to " << str.
name() <<
endl;
245 forAll(cutToMasterEdges, cutEdgeI)
247 if (cutToMasterEdges[cutEdgeI] != -1)
249 const edge& masterEdge = m.edges()[cutToMasterEdges[cutEdgeI]];
250 const edge& cutEdge =
c.edges()[cutEdgeI];
260 str <<
"l " << vertI-3 <<
' ' << vertI-2 <<
nl;
261 str <<
"l " << vertI-3 <<
' ' << vertI-1 <<
nl;
262 str <<
"l " << vertI-3 <<
' ' << vertI <<
nl;
263 str <<
"l " << vertI-2 <<
' ' << vertI-1 <<
nl;
264 str <<
"l " << vertI-2 <<
' ' << vertI <<
nl;
265 str <<
"l " << vertI-1 <<
' ' << vertI <<
nl;
270 OFstream str(
"cutToSlaveEdges.obj");
271 Pout<<
"Writing cutToSlaveEdges to " << str.
name() <<
endl;
279 if (slaveToCut[edgeI] != -1)
281 const edge& slaveEdge =
s.edges()[edgeI];
282 const edge& cutEdge =
c.edges()[slaveToCut[edgeI]];
292 str <<
"l " << vertI-3 <<
' ' << vertI-2 <<
nl;
293 str <<
"l " << vertI-3 <<
' ' << vertI-1 <<
nl;
294 str <<
"l " << vertI-3 <<
' ' << vertI <<
nl;
295 str <<
"l " << vertI-2 <<
' ' << vertI-1 <<
nl;
296 str <<
"l " << vertI-2 <<
' ' << vertI <<
nl;
297 str <<
"l " << vertI-1 <<
' ' << vertI <<
nl;
315 forAll(toPatchEdges, edgeI)
317 const edge&
e = edges[edgeI];
319 label v0 = pointMap[
e[0]];
320 label v1 = pointMap[
e[1]];
322 toPatchEdges[edgeI] =
326 patch.pointEdges()[v0],
335 bool Foam::faceCoupleInfo::regionEdge
337 const polyMesh& slaveMesh,
338 const label slaveEdgeI
341 const labelList& eFaces = slavePatch().edgeFaces()[slaveEdgeI];
343 if (eFaces.size() == 1)
355 label facei = eFaces[i];
357 label meshFacei = slavePatch().addressing()[facei];
359 label patchi = slaveMesh.boundaryMesh().whichPatch(meshFacei);
365 else if (patchi != patch0)
377 Foam::label Foam::faceCoupleInfo::mostAlignedCutEdge
380 const polyMesh& slaveMesh,
381 const bool patchDivision,
385 const label edgeStart,
393 const pointField& localPoints = cutFaces().localPoints();
395 const labelList& pEdges = cutFaces().pointEdges()[pointi];
399 Pout<<
"mostAlignedEdge : finding nearest edge among " 400 << UIndirectList<edge>(cutFaces().edges(), pEdges)
401 <<
" connected to point " << pointi
402 <<
" coord:" << localPoints[pointi]
403 <<
" running between " << edgeStart <<
" coord:" 404 << localPoints[edgeStart]
405 <<
" and " << edgeEnd <<
" coord:" 406 << localPoints[edgeEnd]
413 scalar maxCos = -GREAT;
417 label edgeI = pEdges[i];
423 && cutToMasterEdges[edgeI] == -1
427 && regionEdge(slaveMesh, cutToSlaveEdges[edgeI])
431 const edge&
e = cutFaces().edges()[edgeI];
433 label otherPointi =
e.otherVertex(pointi);
435 if (otherPointi == edgeEnd)
440 Pout<<
" mostAlignedEdge : found end point " << edgeEnd
448 vector eVec(localPoints[otherPointi] - localPoints[pointi]);
450 scalar magEVec =
mag(eVec);
452 if (magEVec < VSMALL)
455 <<
"Crossing zero sized edge " << edgeI
456 <<
" coords:" << localPoints[otherPointi]
457 << localPoints[pointi]
458 <<
" when walking from " << localPoints[edgeStart]
459 <<
" to " << localPoints[edgeEnd]
466 const vector eToEndPoint =
469 localPoints[edgeEnd] - localPoints[otherPointi]
472 scalar cosAngle = eVec & eToEndPoint;
476 Pout<<
" edge:" <<
e <<
" points:" << localPoints[pointi]
477 << localPoints[otherPointi]
479 <<
" vecToEnd:" << eToEndPoint
480 <<
" cosAngle:" << cosAngle
484 if (cosAngle > maxCos)
492 if (maxCos > 1 - angleTol_)
503 void Foam::faceCoupleInfo::setCutEdgeToPoints(
const labelList& cutToMasterEdges)
509 masterPatch().nEdges(),
514 const edgeList& cutEdges = cutFaces().edges();
517 cutEdgeToPoints_.reserve
519 masterPatch().nEdges()
520 + slavePatch().nEdges()
524 forAll(masterToCutEdges, masterEdgeI)
526 const edge& masterE = masterPatch().edges()[masterEdgeI];
531 const labelList& stringedEdges = masterToCutEdges[masterEdgeI];
533 if (stringedEdges.empty())
536 <<
"Did not match all of master edges to cutFace edges" 538 <<
"First unmatched edge:" << masterEdgeI <<
" endPoints:" 539 << masterPatch().localPoints()[masterE[0]]
540 << masterPatch().localPoints()[masterE[1]]
542 <<
"This usually means that the slave patch is not a" 543 <<
" subdivision of the master patch" 546 else if (stringedEdges.size() > 1)
551 DynamicList<label> splitPoints(stringedEdges.size()-1);
554 const edge unsplitEdge
556 masterToCutPoints_[masterE[0]],
557 masterToCutPoints_[masterE[1]]
560 label startVertI = unsplitEdge[0];
561 label startEdgeI = -1;
563 while (startVertI != unsplitEdge[1])
571 label oldStart = startVertI;
575 label edgeI = stringedEdges[i];
577 if (edgeI != startEdgeI)
579 const edge&
e = cutEdges[edgeI];
585 if (
e[0] == startVertI)
589 if (
e[1] != unsplitEdge[1])
591 splitPoints.append(
e[1]);
595 else if (
e[1] == startVertI)
599 if (
e[0] != unsplitEdge[1])
601 splitPoints.append(
e[0]);
609 if (oldStart == startVertI)
612 <<
" unsplitEdge:" << unsplitEdge
613 <<
" does not correspond to split edges " 614 << UIndirectList<edge>(cutEdges, stringedEdges)
629 cutEdgeToPoints_.insert(unsplitEdge, splitPoints.shrink());
635 Foam::label Foam::faceCoupleInfo::matchFaces
642 const bool sameOrientation
645 if (f0.size() != f1.size())
648 <<
"Different sizes for supposedly matching faces." <<
nl 649 <<
"f0:" << f0 <<
" coords:" << UIndirectList<point>(
points0, f0)
651 <<
"f1:" << f1 <<
" coords:" << UIndirectList<point>(points1, f1)
655 const scalar absTolSqr =
sqr(absTol);
663 bool fullMatch =
true;
672 if (distSqr > absTolSqr)
678 fp0 = f0.fcIndex(fp0);
682 fp1 = f1.fcIndex(fp1);
686 fp1 = f1.rcIndex(fp1);
700 <<
"No unique match between two faces" <<
nl 701 <<
"Face " << f0 <<
" coords " 702 << UIndirectList<point>(
points0, f0) <<
nl 703 <<
"Face " << f1 <<
" coords " 704 << UIndirectList<point>(points1, f1)
705 <<
"when using tolerance " << absTol
706 <<
" and forwardMatching:" << sameOrientation
714 bool Foam::faceCoupleInfo::matchPointsThroughFaces
721 const bool sameOrientation,
734 patchToCutPoints.setSize(patchPoints.size());
735 patchToCutPoints = -1;
739 labelList cutPointRegion(cutPoints.size(), -1);
740 DynamicList<label> cutPointRegionMaster;
742 forAll(patchFaces, patchFacei)
744 const face& patchF = patchFaces[patchFacei];
747 const face& cutF = cutFaces[patchFacei];
750 label patchFp = matchFaces
762 label cutPointi = cutF[cutFp];
763 label patchPointi = patchF[patchFp];
775 if (patchToCutPoints[patchPointi] == -1)
777 patchToCutPoints[patchPointi] = cutPointi;
779 else if (patchToCutPoints[patchPointi] != cutPointi)
783 label otherCutPointi = patchToCutPoints[patchPointi];
792 if (cutPointRegion[otherCutPointi] != -1)
795 label region = cutPointRegion[otherCutPointi];
796 cutPointRegion[cutPointi] = region;
799 cutPointRegionMaster[region] =
min 801 cutPointRegionMaster[region],
808 label region = cutPointRegionMaster.size();
809 cutPointRegionMaster.append
811 min(cutPointi, otherCutPointi)
813 cutPointRegion[cutPointi] = region;
814 cutPointRegion[otherCutPointi] = region;
820 patchFp = patchF.fcIndex(patchFp);
824 patchFp = patchF.rcIndex(patchFp);
830 compactToCut.setSize(cutPointRegion.size());
831 cutToCompact.setSize(cutPointRegion.size());
833 label compactPointi = 0;
837 if (cutPointRegion[i] == -1)
840 cutToCompact[i] = compactPointi;
841 compactToCut[compactPointi] = i;
848 label masterPointi = cutPointRegionMaster[cutPointRegion[i]];
850 if (cutToCompact[masterPointi] == -1)
852 cutToCompact[masterPointi] = compactPointi;
853 compactToCut[compactPointi] = masterPointi;
856 cutToCompact[i] = cutToCompact[masterPointi];
859 compactToCut.setSize(compactPointi);
861 return compactToCut.size() != cutToCompact.size();
865 Foam::scalar Foam::faceCoupleInfo::maxDistance
873 scalar maxDist = -GREAT;
877 const point& cutPt = cutPoints[cutF[fp]];
879 pointHit pHit = masterF.nearestPoint(cutPt, masterPoints);
881 maxDist =
max(maxDist, pHit.distance());
887 void Foam::faceCoupleInfo::findPerfectMatchingFaces
889 const primitiveMesh& mesh0,
890 const primitiveMesh& mesh1,
898 if (!mesh0.nFaces() || !mesh1.nFaces())
911 calcFaceCentres<List>
915 mesh0.nInternalFaces(),
916 mesh0.nBoundaryFaces()
922 calcFaceCentres<List>
926 mesh1.nInternalFaces(),
927 mesh1.nBoundaryFaces()
934 Pout<<
"Face matching tolerance : " << absTol <<
endl;
952 <<
"Matched ALL " << fc1.size()
953 <<
" boundary faces of mesh0 to boundary faces of mesh1." <<
endl 954 <<
"This is only valid if the mesh to add is fully" 955 <<
" enclosed by the mesh it is added to." <<
endl;
962 mesh0Faces.setSize(fc0.size());
963 mesh1Faces.setSize(fc1.size());
967 if (from1To0[i] != -1)
969 mesh1Faces[nMatched] = i + mesh1.nInternalFaces();
970 mesh0Faces[nMatched] = from1To0[i] + mesh0.nInternalFaces();
976 mesh0Faces.setSize(nMatched);
977 mesh1Faces.setSize(nMatched);
981 void Foam::faceCoupleInfo::findSlavesCoveringMaster
983 const primitiveMesh& mesh0,
984 const primitiveMesh& mesh1,
995 treeBoundBox overallBb(mesh0.points());
997 indexedOctree<treeDataFace>
tree 1003 identity(mesh0.nBoundaryFaces(), mesh0.nInternalFaces())
1006 overallBb.extend(
rndGen, 1
e-4),
1014 Pout<<
"findSlavesCoveringMaster :" 1015 <<
" constructed octree for mesh0 boundary faces" <<
endl;
1025 label mesh1Facei = mesh1.nInternalFaces();
1026 mesh1Facei < mesh1.nFaces();
1030 const face& f1 = mesh1.faces()[mesh1Facei];
1033 point fc(f1.centre(mesh1.points()));
1039 label mesh0Facei =
tree.shapes().objectIndex(nearInfo.index());
1052 mesh0.faces()[mesh0Facei],
1058 mesh0Set.insert(mesh0Facei);
1059 mesh1Set.insert(mesh1Facei);
1066 Pout<<
"findSlavesCoveringMaster :" 1067 <<
" matched " << mesh1Set.size() <<
" mesh1 faces to " 1068 << mesh0Set.size() <<
" mesh0 faces" <<
endl;
1071 mesh0Faces = mesh0Set.toc();
1072 mesh1Faces = mesh1Set.toc();
1076 Foam::label Foam::faceCoupleInfo::growCutFaces
1079 Map<labelList>& candidates
1084 Pout<<
"growCutFaces :" 1085 <<
" growing cut faces to masterPatch" <<
endl;
1088 label nTotChanged = 0;
1097 forAll(cutToMasterFaces_, cutFacei)
1099 const label masterFacei = cutToMasterFaces_[cutFacei];
1101 if (masterFacei != -1)
1107 const labelList& fEdges = cutFaceEdges[cutFacei];
1111 const label cutEdgeI = fEdges[i];
1113 if (cutToMasterEdges[cutEdgeI] == -1)
1121 const labelList& eFaces = cutEdgeFaces[cutEdgeI];
1125 const label facei = eFaces[j];
1127 if (cutToMasterFaces_[facei] == -1)
1129 cutToMasterFaces_[facei] = masterFacei;
1130 candidates.erase(facei);
1133 else if (cutToMasterFaces_[facei] != masterFacei)
1136 cutFaces().points();
1138 masterPatch().points();
1140 const edge&
e = cutFaces().edges()[cutEdgeI];
1142 label myMaster = cutToMasterFaces_[facei];
1143 const face& myF = masterPatch()[myMaster];
1145 const face& nbrF = masterPatch()[masterFacei];
1149 << cutFaces()[facei].points(cutPoints)
1152 <<
" but also connects to nbr face " 1153 << cutFaces()[cutFacei].points(cutPoints)
1154 <<
" with master " << masterFacei
1157 << myF.points(masterPoints)
1158 <<
" nbrMasterFace:" 1159 << nbrF.points(masterPoints) <<
nl 1160 <<
"Across cut edge " << cutEdgeI
1162 << cutFaces().localPoints()[
e[0]]
1163 << cutFaces().localPoints()[
e[1]]
1174 Pout<<
"growCutFaces : Grown an additional " << nChanged
1175 <<
" cut to master face" <<
" correspondence" <<
endl;
1178 nTotChanged += nChanged;
1190 void Foam::faceCoupleInfo::checkMatch(
const labelList& cutToMasterEdges)
const 1192 const pointField& cutLocalPoints = cutFaces().localPoints();
1194 const pointField& masterLocalPoints = masterPatch().localPoints();
1195 const faceList& masterLocalFaces = masterPatch().localFaces();
1197 forAll(cutToMasterEdges, cutEdgeI)
1199 const edge&
e = cutFaces().edges()[cutEdgeI];
1201 if (cutToMasterEdges[cutEdgeI] == -1)
1204 const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1206 label masterFacei = -1;
1210 label cutFacei = cutEFaces[i];
1212 if (cutToMasterFaces_[cutFacei] != -1)
1214 if (masterFacei == -1)
1216 masterFacei = cutToMasterFaces_[cutFacei];
1218 else if (masterFacei != cutToMasterFaces_[cutFacei])
1220 label myMaster = cutToMasterFaces_[cutFacei];
1221 const face& myF = masterLocalFaces[myMaster];
1223 const face& nbrF = masterLocalFaces[masterFacei];
1226 <<
"Internal CutEdge " <<
e 1228 << cutLocalPoints[
e[0]]
1229 << cutLocalPoints[
e[1]]
1230 <<
" connects to master " << myMaster
1231 <<
" and to master " << masterFacei <<
nl 1233 << myF.points(masterLocalPoints)
1234 <<
" nbrMasterFace:" 1235 << nbrF.points(masterLocalPoints)
1245 Foam::label Foam::faceCoupleInfo::matchEdgeFaces
1248 Map<labelList>& candidates
1257 candidates.reserve(cutFaces().size());
1261 forAll(cutToMasterEdges, cutEdgeI)
1263 label masterEdgeI = cutToMasterEdges[cutEdgeI];
1265 if (masterEdgeI != -1)
1270 const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1272 masterPatch().edgeFaces()[masterEdgeI];
1276 label cutFacei = cutEFaces[i];
1278 if (cutToMasterFaces_[cutFacei] == -1)
1287 auto fnd = candidates.find(cutFacei);
1293 candidates.insert(cutFacei, masterEFaces);
1300 const labelList& masterFaces = fnd.val();
1302 DynamicList<label> newCandidates(masterFaces.size());
1306 if (masterFaces.found(masterEFaces[j]))
1308 newCandidates.append(masterEFaces[j]);
1312 if (newCandidates.size() == 1)
1315 cutToMasterFaces_[cutFacei] = newCandidates[0];
1316 candidates.erase(cutFacei);
1329 fnd() = newCandidates.shrink();
1340 Pout<<
"matchEdgeFaces : Found " << nChanged
1341 <<
" faces where there was" 1342 <<
" only one remaining choice for cut-master correspondence" 1350 Foam::label Foam::faceCoupleInfo::geometricMatchEdgeFaces
1352 Map<labelList>& candidates
1355 const pointField& cutPoints = cutFaces().points();
1365 masterPatch().size(),
1372 label cutFacei = iter.key();
1373 const labelList& masterFaces = iter.val();
1375 const face& cutF = cutFaces()[cutFacei];
1377 if (cutToMasterFaces_[cutFacei] == -1)
1380 scalar minDist = GREAT;
1381 label minMasterFacei = -1;
1385 label masterFacei = masterFaces[i];
1387 if (masterToCutFaces[masterFaces[i]].empty())
1389 scalar dist = maxDistance
1393 masterPatch()[masterFacei],
1400 minMasterFacei = masterFacei;
1405 if (minMasterFacei != -1)
1407 cutToMasterFaces_[cutFacei] = minMasterFacei;
1408 masterToCutFaces[minMasterFacei] = cutFacei;
1415 forAll(cutToMasterFaces_, cutFacei)
1417 if (cutToMasterFaces_[cutFacei] != -1)
1419 candidates.erase(cutFacei);
1426 Pout<<
"geometricMatchEdgeFaces : Found " << nChanged
1427 <<
" faces where there was" 1428 <<
" only one remaining choice for cut-master correspondence" 1436 void Foam::faceCoupleInfo::perfectPointMatch
1438 const scalar absTol,
1439 const bool slaveFacesOrdered
1447 Pout<<
"perfectPointMatch :" 1448 <<
" Matching master and slave to cut." 1449 <<
" Master and slave faces are identical" <<
nl;
1451 if (slaveFacesOrdered)
1453 Pout<<
"and master and slave faces are ordered" 1454 <<
" (on coupled patches)" <<
endl;
1458 Pout<<
"and master and slave faces are not ordered" 1459 <<
" (on coupled patches)" <<
endl;
1463 cutToMasterFaces_ =
identity(masterPatch().size());
1464 cutPoints_ = masterPatch().localPoints();
1469 masterPatch().localFaces(),
1473 masterToCutPoints_ =
identity(cutPoints_.size());
1477 bool matchedAllFaces =
false;
1479 if (slaveFacesOrdered)
1481 cutToSlaveFaces_ =
identity(cutFaces().size());
1482 matchedAllFaces = (cutFaces().size() == slavePatch().size());
1492 calcFaceCentres<List>
1499 calcFaceCentres<IndirectList>
1515 if (!matchedAllFaces)
1517 labelList cutToSlaveFacesTemp(cutToSlaveFaces_.size(), -1);
1521 calcFacePointAverages<List>
1528 calcFacePointAverages<IndirectList>
1540 cutToSlaveFaces_ =
max(cutToSlaveFaces_, cutToSlaveFacesTemp);
1542 matchedAllFaces =
min(cutToSlaveFaces_) != -1;
1547 if (!matchedAllFaces)
1550 <<
"Did not match all of the master faces to the slave faces" 1552 <<
"This usually means that the slave patch and master patch" 1553 <<
" do not align to within " << absTol <<
" metre." 1563 matchPointsThroughFaces
1566 cutFaces().localPoints(),
1567 reorder(cutToSlaveFaces_, cutFaces().localFaces()),
1568 slavePatch().localPoints(),
1569 slavePatch().localFaces(),
1579 cutPoints_ = UIndirectList<point>(cutPoints_, compactToCut)();
1581 const faceList& cutLocalFaces = cutFaces().localFaces();
1583 faceList compactFaces(cutLocalFaces.size());
1586 compactFaces[i] =
renumber(cutToCompact, cutLocalFaces[i]);
1602 void Foam::faceCoupleInfo::subDivisionMatch
1604 const polyMesh& slaveMesh,
1605 const bool patchDivision,
1611 Pout<<
"subDivisionMatch :" 1612 <<
" Matching master and slave to cut." 1613 <<
" Slave can be subdivision of master but all master edges" 1614 <<
" have to be covered by slave edges." <<
endl;
1619 cutPoints_ = slavePatch().localPoints();
1621 faceList cutFaces(slavePatch().size());
1625 cutFaces[i] = slavePatch().localFaces()[i].reverseFace();
1631 cutToSlaveFaces_ =
identity(cutFaces().size());
1648 OFstream str(
"regionEdges.obj");
1650 Pout<<
"subDivisionMatch :" 1651 <<
" addressing for slave patch fully done." 1652 <<
" Dumping region edges to " << str.
name() <<
endl;
1656 forAll(slavePatch().edges(), slaveEdgeI)
1658 if (regionEdge(slaveMesh, slaveEdgeI))
1660 const edge&
e = slavePatch().edges()[slaveEdgeI];
1666 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
1679 Pout<<
"subDivisionMatch :" 1680 <<
" matching master points to cut points" <<
endl;
1685 masterPatch().localPoints(),
1692 if (!matchedAllPoints)
1695 <<
"Did not match all of the master points to the slave points" 1697 <<
"This usually means that the slave patch is not a" 1698 <<
" subdivision of the master patch" 1710 Pout<<
"subDivisionMatch :" 1711 <<
" matching cut edges to master edges" <<
endl;
1714 const edgeList& masterEdges = masterPatch().edges();
1715 const pointField& masterPoints = masterPatch().localPoints();
1717 const edgeList& cutEdges = cutFaces().edges();
1719 labelList cutToMasterEdges(cutFaces().nEdges(), -1);
1721 forAll(masterEdges, masterEdgeI)
1723 const edge& masterEdge = masterEdges[masterEdgeI];
1725 label cutPoint0 = masterToCutPoints_[masterEdge[0]];
1726 label cutPoint1 = masterToCutPoints_[masterEdge[1]];
1730 label cutPointi = cutPoint0;
1763 Pout<<
"Dumping unmatched pointEdges to errorEdges.obj" 1774 cutFaces().pointEdges()[cutPointi]
1777 cutFaces().localPoints(),
1782 <<
"Problem in finding cut edges corresponding to" 1783 <<
" master edge " << masterEdge
1784 <<
" points:" << masterPoints[masterEdge[0]]
1785 <<
' ' << masterPoints[masterEdge[1]]
1786 <<
" corresponding cut edge: (" << cutPoint0
1787 <<
") " << cutPoint1
1791 cutToMasterEdges[cutEdgeI] = masterEdgeI;
1793 cutPointi = cutEdges[cutEdgeI].otherVertex(cutPointi);
1795 }
while (cutPointi != cutPoint1);
1801 writeEdges(cutToMasterEdges, cutToSlaveEdges);
1806 setCutEdgeToPoints(cutToMasterEdges);
1820 Pout<<
"subDivisionMatch :" 1821 <<
" matching (topological) cut faces to masterPatch" <<
endl;
1825 Map<labelList> candidates(cutFaces().size());
1827 cutToMasterFaces_.setSize(cutFaces().size());
1828 cutToMasterFaces_ = -1;
1834 label nChanged = matchEdgeFaces(cutToMasterEdges, candidates);
1836 checkMatch(cutToMasterEdges);
1843 nChanged += growCutFaces(cutToMasterEdges, candidates);
1845 checkMatch(cutToMasterEdges);
1855 Pout<<
"subDivisionMatch :" 1856 <<
" matching (geometric) cut faces to masterPatch" <<
endl;
1865 label nChanged = geometricMatchEdgeFaces(candidates);
1867 checkMatch(cutToMasterEdges);
1869 nChanged += growCutFaces(cutToMasterEdges, candidates);
1871 checkMatch(cutToMasterEdges);
1881 forAll(cutToMasterFaces_, cutFacei)
1883 if (cutToMasterFaces_[cutFacei] == -1)
1885 const face& cutF = cutFaces()[cutFacei];
1888 <<
"Did not match all of cutFaces to a master face" <<
nl 1889 <<
"First unmatched cut face:" << cutFacei <<
" with points:" 1890 << UIndirectList<point>(cutFaces().points(), cutF)
1892 <<
"This usually means that the slave patch is not a" 1893 <<
" subdivision of the master patch" 1900 Pout<<
"subDivisionMatch :" 1901 <<
" finished matching master and slave to cut" <<
endl;
1915 const polyMesh& masterMesh,
1916 const polyMesh& slaveMesh,
1917 const scalar absTol,
1918 const bool perfectMatch
1921 masterPatchPtr_(nullptr),
1922 slavePatchPtr_(nullptr),
1924 cutFacesPtr_(nullptr),
1925 cutToMasterFaces_(0),
1926 masterToCutPoints_(0),
1927 cutToSlaveFaces_(0),
1928 slaveToCutPoints_(0),
1940 findPerfectMatchingFaces
1953 findSlavesCoveringMaster
1964 masterPatchPtr_.reset
1968 IndirectList<face>(masterMesh.faces(), masterToMesh),
1973 slavePatchPtr_.reset
1977 IndirectList<face>(slaveMesh.faces(), slaveToMesh),
1986 perfectPointMatch(absTol,
false);
1991 subDivisionMatch(slaveMesh,
false, absTol);
2003 const polyMesh& masterMesh,
2005 const polyMesh& slaveMesh,
2007 const scalar absTol,
2008 const bool perfectMatch,
2009 const bool orderedFaces,
2010 const bool patchDivision
2017 IndirectList<face>(masterMesh.faces(), masterAddressing),
2025 IndirectList<face>(slaveMesh.faces(), slaveAddressing),
2030 cutFacesPtr_(nullptr),
2031 cutToMasterFaces_(0),
2032 masterToCutPoints_(0),
2033 cutToSlaveFaces_(0),
2034 slaveToCutPoints_(0),
2037 if (perfectMatch && (masterAddressing.size() != slaveAddressing.size()))
2040 <<
"Perfect match specified but number of master and slave faces" 2041 <<
" differ." <<
endl 2042 <<
"master:" << masterAddressing.size()
2043 <<
" slave:" << slaveAddressing.size()
2049 masterAddressing.size()
2050 &&
min(masterAddressing) < masterMesh.nInternalFaces()
2054 <<
"Supplied internal face on master mesh to couple." <<
nl 2055 <<
"Faces to be coupled have to be boundary faces." 2060 slaveAddressing.size()
2061 &&
min(slaveAddressing) < slaveMesh.nInternalFaces()
2065 <<
"Supplied internal face on slave mesh to couple." <<
nl 2066 <<
"Faces to be coupled have to be boundary faces." 2073 perfectPointMatch(absTol, orderedFaces);
2078 subDivisionMatch(slaveMesh, patchDivision, absTol);
2094 label facei =
pp.start();
2106 Map<label> map(lst.size());
2112 map.insert(i, lst[i]);
2124 Map<labelList> map(lst.size());
2130 map.insert(i, lst[i]);
static labelList faceLabels(const polyPatch &)
Utility functions.
faceCoupleInfo(const polyMesh &mesh0, const polyMesh &mesh1, const scalar absTol, const bool perfectMatch)
Construct from two meshes and absolute tolerance.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void resize(const label len)
Adjust allocated size of list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
List< edge > edgeList
List of edge.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
IntListType renumber(const labelUList &oldToNew, const IntListType &input)
Renumber the values within a list.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field...
virtual const fileName & name() const override
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
List< labelList > labelListList
List of labelList.
Determine correspondence between points. See below.
#define forAll(list, i)
Loop across all elements in list.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
List< face > faceList
List of faces.
bool matchPoints(const UList< point > &pts0, const UList< point > &pts1, const UList< scalar > &matchDistance, const bool verbose, labelList &from0To1, const point &origin=point::zero)
Determine correspondence between pointFields. Gets passed.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Tree tree(triangles.begin(), triangles.end())
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
static Map< label > makeMap(const labelList &)
Create Map from List.
PrimitivePatch< List< face >, const pointField & > primitiveFacePatch
A PrimitivePatch with List storage for the faces, const reference for the point field.
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
pointField points0(pointIOField(IOobject("points", mesh.time().constant(), polyMesh::meshSubDir, mesh, IOobject::MUST_READ, IOobject::NO_WRITE, IOobject::NO_REGISTER)))
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values within a list.
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar c
Speed of light in a vacuum.
const std::string patch
OpenFOAM patch number as a std::string.
List< label > labelList
A List of labels.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
ListType reorder(const labelUList &oldToNew, const ListType &input, const bool prune=false)
Reorder the elements of a list.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
PointHit< point > pointHit
A PointHit with a 3D point.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
forAllConstIters(mixture.phases(), phase)
A HashTable to objects of type <T> with a label key.
static constexpr const zero Zero
Global zero (0)