1 /*---------------------------------------------------------------------------* \
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*----------------------------------------------------------------------------*/
27 #include "faceCoupleInfo.H"
29 #include "matchPoints.H"
30 #include "indirectPrimitivePatch.H"
31 #include "meshTools.H"
32 #include "octreeDataFace.H"
35 #include "IndirectList.H"
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 defineTypeNameAndDebug(faceCoupleInfo, 0);
46 const Foam::scalar Foam::faceCoupleInfo::angleTol_ = 1E-3;
49 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
52 void Foam::faceCoupleInfo::writeOBJ
54 const fileName& fName,
55 const edgeList& edges,
56 const pointField& points,
62 labelList pointMap(points.size(), -1);
70 const edge& e = edges[edgeI];
76 if (pointMap[pointI] == -1)
78 pointMap[pointI] = newPointI++;
80 meshTools::writeOBJ(str, points[pointI]);
87 forAll(points, pointI)
89 meshTools::writeOBJ(str, points[pointI]);
92 pointMap = identity(points.size());
97 const edge& e = edges[edgeI];
99 str<< "l " << pointMap[e[0]]+1 << ' ' << pointMap[e[1]]+1 << nl;
105 void Foam::faceCoupleInfo::writeOBJ
107 const fileName& fName,
108 const pointField& points0,
109 const pointField& points1
112 Pout<< "Writing connections as edges to " << fName << endl;
120 meshTools::writeOBJ(str, points0[i]);
122 meshTools::writeOBJ(str, points1[i]);
124 str << "l " << vertI-1 << ' ' << vertI << nl;
129 //- Writes face and point connectivity as .obj files.
130 void Foam::faceCoupleInfo::writePointsFaces() const
132 const indirectPrimitivePatch& m = masterPatch();
133 const indirectPrimitivePatch& s = slavePatch();
134 const primitiveFacePatch& c = cutFaces();
138 OFstream str("masterPatch.obj");
139 Pout<< "Writing masterPatch to " << str.name() << endl;
140 meshTools::writeOBJ(str, m.localFaces(), m.localPoints());
143 OFstream str("slavePatch.obj");
144 Pout<< "Writing slavePatch to " << str.name() << endl;
145 meshTools::writeOBJ(str, s.localFaces(), s.localPoints());
148 OFstream str("cutFaces.obj");
149 Pout<< "Writing cutFaces to " << str.name() << endl;
150 meshTools::writeOBJ(str, c.localFaces(), c.localPoints());
153 // Point connectivity
155 Pout<< "Writing cutToMasterPoints to cutToMasterPoints.obj" << endl;
159 "cutToMasterPoints.obj",
163 IndirectList<point>(c.localPoints(), masterToCutPoints_)()
168 Pout<< "Writing cutToSlavePoints to cutToSlavePoints.obj" << endl;
172 "cutToSlavePoints.obj",
176 IndirectList<point>(c.localPoints(), slaveToCutPoints_)()
183 Pout<< "Writing cutToMasterFaces to cutToMasterFaces.obj" << endl;
185 pointField equivMasterFaces(c.size());
187 forAll(cutToMasterFaces(), cutFaceI)
189 label masterFaceI = cutToMasterFaces()[cutFaceI];
191 if (masterFaceI != -1)
193 equivMasterFaces[cutFaceI] = m[masterFaceI].centre(m.points());
197 WarningIn("writePointsFaces()")
198 << "No master face for cut face " << cutFaceI
199 << " at position " << c[cutFaceI].centre(c.points())
202 equivMasterFaces[cutFaceI] = vector::zero;
208 "cutToMasterFaces.obj",
209 calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
215 Pout<< "Writing cutToSlaveFaces to cutToSlaveFaces.obj" << endl;
217 pointField equivSlaveFaces(c.size());
219 forAll(cutToSlaveFaces(), cutFaceI)
221 label slaveFaceI = cutToSlaveFaces()[cutFaceI];
223 equivSlaveFaces[cutFaceI] = s[slaveFaceI].centre(s.points());
228 "cutToSlaveFaces.obj",
229 calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
238 void Foam::faceCoupleInfo::writeEdges
240 const labelList& cutToMasterEdges,
241 const labelList& cutToSlaveEdges
244 const indirectPrimitivePatch& m = masterPatch();
245 const indirectPrimitivePatch& s = slavePatch();
246 const primitiveFacePatch& c = cutFaces();
250 OFstream str("cutToMasterEdges.obj");
251 Pout<< "Writing cutToMasterEdges to " << str.name() << endl;
255 forAll(cutToMasterEdges, cutEdgeI)
257 if (cutToMasterEdges[cutEdgeI] != -1)
259 const edge& masterEdge =
260 m.edges()[cutToMasterEdges[cutEdgeI]];
261 const edge& cutEdge = c.edges()[cutEdgeI];
263 meshTools::writeOBJ(str, m.localPoints()[masterEdge[0]]);
265 meshTools::writeOBJ(str, m.localPoints()[masterEdge[1]]);
267 meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
269 meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
271 str << "l " << vertI-3 << ' ' << vertI-2 << nl;
272 str << "l " << vertI-3 << ' ' << vertI-1 << nl;
273 str << "l " << vertI-3 << ' ' << vertI << nl;
274 str << "l " << vertI-2 << ' ' << vertI-1 << nl;
275 str << "l " << vertI-2 << ' ' << vertI << nl;
276 str << "l " << vertI-1 << ' ' << vertI << nl;
281 OFstream str("cutToSlaveEdges.obj");
282 Pout<< "Writing cutToSlaveEdges to " << str.name() << endl;
286 labelList slaveToCut(invert(s.nEdges(), cutToSlaveEdges));
288 forAll(slaveToCut, edgeI)
290 if (slaveToCut[edgeI] != -1)
292 const edge& slaveEdge = s.edges()[edgeI];
293 const edge& cutEdge = c.edges()[slaveToCut[edgeI]];
295 meshTools::writeOBJ(str, s.localPoints()[slaveEdge[0]]);
297 meshTools::writeOBJ(str, s.localPoints()[slaveEdge[1]]);
299 meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
301 meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
303 str << "l " << vertI-3 << ' ' << vertI-2 << nl;
304 str << "l " << vertI-3 << ' ' << vertI-1 << nl;
305 str << "l " << vertI-3 << ' ' << vertI << nl;
306 str << "l " << vertI-2 << ' ' << vertI-1 << nl;
307 str << "l " << vertI-2 << ' ' << vertI << nl;
308 str << "l " << vertI-1 << ' ' << vertI << nl;
317 // Given an edgelist and a map for the points on the edges it tries to find
318 // the corresponding patch edges.
319 Foam::labelList Foam::faceCoupleInfo::findMappedEdges
321 const edgeList& edges,
322 const labelList& pointMap,
323 const indirectPrimitivePatch& patch
326 labelList toPatchEdges(edges.size());
328 forAll(toPatchEdges, edgeI)
330 const edge& e = edges[edgeI];
332 label v0 = pointMap[e[0]];
333 label v1 = pointMap[e[1]];
335 toPatchEdges[edgeI] =
339 patch.pointEdges()[v0],
348 // Detect a cut edge which originates from two boundary faces having different
350 bool Foam::faceCoupleInfo::regionEdge
352 const polyMesh& slaveMesh,
353 const label slaveEdgeI
356 const labelList& eFaces = slavePatch().edgeFaces()[slaveEdgeI];
358 if (eFaces.size() == 1)
364 // Count how many different patches connected to this edge.
370 label faceI = eFaces[i];
372 label meshFaceI = slavePatch().addressing()[faceI];
374 label patchI = slaveMesh.boundaryMesh().whichPatch(meshFaceI);
380 else if (patchI != patch0)
382 // Found two different patches connected to this edge.
391 // Find edge using pointI that is most aligned with vector between
392 // master points. Patchdivision tells us whether or not to use
393 // patch information to match edges.
394 Foam::label Foam::faceCoupleInfo::mostAlignedCutEdge
397 const polyMesh& slaveMesh,
398 const bool patchDivision,
399 const labelList& cutToMasterEdges,
400 const labelList& cutToSlaveEdges,
402 const label edgeStart,
406 const pointField& localPoints = cutFaces().localPoints();
408 const labelList& pEdges = cutFaces().pointEdges()[pointI];
412 Pout<< "mostAlignedEdge : finding nearest edge among "
413 << IndirectList<edge>(cutFaces().edges(), pEdges)()
414 << " connected to point " << pointI
415 << " coord:" << localPoints[pointI]
416 << " running between " << edgeStart << " coord:"
417 << localPoints[edgeStart]
418 << " and " << edgeEnd << " coord:"
419 << localPoints[edgeEnd]
423 // Find the edge that gets us nearest end.
426 scalar maxCos = -GREAT;
430 label edgeI = pEdges[i];
436 && cutToMasterEdges[edgeI] == -1
440 && regionEdge(slaveMesh, cutToSlaveEdges[edgeI])
444 const edge& e = cutFaces().edges()[edgeI];
446 label otherPointI = e.otherVertex(pointI);
448 if (otherPointI == edgeEnd)
450 // Shortcut: found edge end point.
453 Pout<< " mostAlignedEdge : found end point " << edgeEnd
459 // Get angle between edge and edge to masterEnd
461 vector eVec(localPoints[otherPointI] - localPoints[pointI]);
463 scalar magEVec = mag(eVec);
465 if (magEVec < VSMALL)
467 WarningIn("faceCoupleInfo::mostAlignedEdge")
468 << "Crossing zero sized edge " << edgeI
469 << " coords:" << localPoints[otherPointI]
470 << localPoints[pointI]
471 << " when walking from " << localPoints[edgeStart]
472 << " to " << localPoints[edgeEnd]
479 vector eToEndPoint(localPoints[edgeEnd] - localPoints[otherPointI]);
480 eToEndPoint /= mag(eToEndPoint);
482 scalar cosAngle = eVec & eToEndPoint;
486 Pout<< " edge:" << e << " points:" << localPoints[pointI]
487 << localPoints[otherPointI]
489 << " vecToEnd:" << eToEndPoint
490 << " cosAngle:" << cosAngle
494 if (cosAngle > maxCos)
502 if (maxCos > 1 - angleTol_)
513 // Construct points to split points map (in cut addressing)
514 void Foam::faceCoupleInfo::setCutEdgeToPoints(const labelList& cutToMasterEdges)
516 labelListList masterToCutEdges
520 masterPatch().nEdges(),
525 const edgeList& cutEdges = cutFaces().edges();
527 // Size extra big so searching is faster
528 cutEdgeToPoints_.resize
530 masterPatch().nEdges()
531 + slavePatch().nEdges()
535 forAll(masterToCutEdges, masterEdgeI)
537 const edge& masterE = masterPatch().edges()[masterEdgeI];
539 //Pout<< "Master:" << masterPatch().localPoints()[masterE[0]] << ' '
540 // << masterPatch().localPoints()[masterE[1]] << endl;
542 const labelList& stringedEdges = masterToCutEdges[masterEdgeI];
544 if (stringedEdges.size() == 0)
548 "faceCoupleInfo::setCutEdgeToPoints"
550 ) << "Did not match all of master edges to cutFace edges"
552 << "First unmatched edge:" << masterEdgeI << " endPoints:"
553 << masterPatch().localPoints()[masterE[0]]
554 << masterPatch().localPoints()[masterE[1]]
556 << "This usually means that the slave patch is not a"
557 << " subdivision of the master patch"
558 << abort(FatalError);
560 else if (stringedEdges.size() > 1)
562 // String up the edges between e[0] and e[1]. Store the points
563 // inbetween e[0] and e[1] (all in cutFaces() labels)
565 DynamicList<label> splitPoints(stringedEdges.size()-1);
567 // Unsplit edge endpoints
568 const edge unsplitEdge
570 masterToCutPoints_[masterE[0]],
571 masterToCutPoints_[masterE[1]]
574 label startVertI = unsplitEdge[0];
575 label startEdgeI = -1;
577 while (startVertI != unsplitEdge[1])
579 // Loop over all string of edges. Update
580 // - startVertI : previous vertex
581 // - startEdgeI : previous edge
582 // and insert any points into splitPoints
585 label oldStart = startVertI;
587 forAll(stringedEdges, i)
589 label edgeI = stringedEdges[i];
591 if (edgeI != startEdgeI)
593 const edge& e = cutEdges[edgeI];
595 //Pout<< " cut:" << e << " at:"
596 // << cutFaces().localPoints()[e[0]]
597 // << ' ' << cutFaces().localPoints()[e[1]] << endl;
599 if (e[0] == startVertI)
603 if (e[1] != unsplitEdge[1])
605 splitPoints.append(e[1]);
609 else if (e[1] == startVertI)
613 if (e[0] != unsplitEdge[1])
615 splitPoints.append(e[0]);
623 if (oldStart == startVertI)
627 "faceCoupleInfo::setCutEdgeToPoints"
629 ) << " unsplitEdge:" << unsplitEdge
630 << " does not correspond to split edges "
631 << IndirectList<edge>(cutEdges, stringedEdges)()
632 << abort(FatalError);
636 //Pout<< "For master edge:"
638 // << " Found stringed points "
639 // << IndirectList<point>
641 // cutFaces().localPoints(),
642 // splitPoints.shrink()
646 cutEdgeToPoints_.insert(unsplitEdge, splitPoints.shrink());
652 // Determines rotation for f1 to match up with f0, i.e. the index in f0 of
653 // the first point of f1.
654 Foam::label Foam::faceCoupleInfo::matchFaces
657 const pointField& points0,
659 const pointField& points1,
661 const bool sameOrientation
664 if (f0.size() != f1.size())
668 "faceCoupleInfo::matchFaces"
669 "(const scalar, const face&, const pointField&"
670 ", const face&, const pointField&)"
671 ) << "Different sizes for supposedly matching faces." << nl
672 << "f0:" << f0 << " coords:" << IndirectList<point>(points0, f0)()
674 << "f1:" << f1 << " coords:" << IndirectList<point>(points1, f1)()
675 << abort(FatalError);
678 const scalar absTolSqr = sqr(absTol);
685 // See -if starting from startFp on f0- the two faces match.
686 bool fullMatch = true;
693 scalar distSqr = Foam::magSqr(points0[f0[fp0]] - points1[f1[fp1]]);
695 if (distSqr > absTolSqr)
701 fp0 = f0.fcIndex(fp0); // walk forward
705 fp1 = f1.fcIndex(fp1);
709 fp1 = f1.rcIndex(fp1);
724 "faceCoupleInfo::matchFaces"
725 "(const scalar, const face&, const pointField&"
726 ", const face&, const pointField&)"
727 ) << "No unique match between two faces" << nl
728 << "Face " << f0 << " coords "
729 << IndirectList<point>(points0, f0)()
731 << "Face " << f1 << " coords "
732 << IndirectList<point>(points1, f1)()
733 << "when using tolerance " << absTol
734 << " and forwardMatching:" << sameOrientation
735 << abort(FatalError);
742 // Find correspondence from patch points to cut points. This might
743 // detect shared points so the output is a patch-to-cut point list
744 // and a compaction list for the cut points (which will always be equal or more
745 // connected than the patch).
746 // Returns true if there are any duplicates.
747 bool Foam::faceCoupleInfo::matchPointsThroughFaces
750 const pointField& cutPoints,
751 const faceList& cutFaces,
752 const pointField& patchPoints,
753 const faceList& patchFaces,
754 const bool sameOrientation,
756 labelList& patchToCutPoints, // patch to (uncompacted) cut points
757 labelList& cutToCompact, // compaction list for cut points
758 labelList& compactToCut // inverse ,,
762 // From slave to cut point
763 patchToCutPoints.setSize(patchPoints.size());
764 patchToCutPoints = -1;
766 // Compaction list for cut points: either -1 or index into master which
767 // gives the point to compact to.
768 labelList cutPointRegion(cutPoints.size(), -1);
769 DynamicList<label> cutPointRegionMaster;
771 forAll(patchFaces, patchFaceI)
773 const face& patchF = patchFaces[patchFaceI];
775 //const face& cutF = cutFaces[patchToCutFaces[patchFaceI]];
776 const face& cutF = cutFaces[patchFaceI];
778 // Do geometric matching to get position of cutF[0] in patchF
779 label patchFp = matchFaces
786 sameOrientation // orientation
791 label cutPointI = cutF[cutFp];
792 label patchPointI = patchF[patchFp];
794 //const point& cutPt = cutPoints[cutPointI];
795 //const point& patchPt = patchPoints[patchPointI];
796 //if (mag(cutPt - patchPt) > SMALL)
798 // FatalErrorIn("matchPointsThroughFaces")
799 // << "cutP:" << cutPt
800 // << " patchP:" << patchPt
801 // << abort(FatalError);
804 if (patchToCutPoints[patchPointI] == -1)
806 patchToCutPoints[patchPointI] = cutPointI;
808 else if (patchToCutPoints[patchPointI] != cutPointI)
810 // Multiple cut points connecting to same patch.
811 // Check if already have region & region master for this set
812 label otherCutPointI = patchToCutPoints[patchPointI];
814 //Pout<< "PatchPoint:" << patchPt
815 // << " matches to:" << cutPointI
816 // << " coord:" << cutPoints[cutPointI]
817 // << " and to:" << otherCutPointI
818 // << " coord:" << cutPoints[otherCutPointI]
821 if (cutPointRegion[otherCutPointI] != -1)
823 // Have region for this set. Copy.
824 label region = cutPointRegion[otherCutPointI];
825 cutPointRegion[cutPointI] = region;
827 // Update region master with min point label
828 cutPointRegionMaster[region] = min
830 cutPointRegionMaster[region],
836 // Create new region.
837 label region = cutPointRegionMaster.size();
838 cutPointRegionMaster.append
840 min(cutPointI, otherCutPointI)
842 cutPointRegion[cutPointI] = region;
843 cutPointRegion[otherCutPointI] = region;
849 patchFp = patchF.fcIndex(patchFp);
853 patchFp = patchF.rcIndex(patchFp);
858 // Rework region&master into compaction array
859 compactToCut.setSize(cutPointRegion.size());
860 cutToCompact.setSize(cutPointRegion.size());
862 label compactPointI = 0;
864 forAll(cutPointRegion, i)
866 if (cutPointRegion[i] == -1)
868 // Unduplicated point. Allocate new compacted point.
869 cutToCompact[i] = compactPointI;
870 compactToCut[compactPointI] = i;
875 // Duplicate point. Get master.
877 label masterPointI = cutPointRegionMaster[cutPointRegion[i]];
879 if (cutToCompact[masterPointI] == -1)
881 cutToCompact[masterPointI] = compactPointI;
882 compactToCut[compactPointI] = masterPointI;
885 cutToCompact[i] = cutToCompact[masterPointI];
888 compactToCut.setSize(compactPointI);
890 return compactToCut.size() != cutToCompact.size();
894 // Return max distance from any point on cutF to masterF
895 Foam::scalar Foam::faceCoupleInfo::maxDistance
898 const pointField& cutPoints,
900 const pointField& masterPoints
903 scalar maxDist = -GREAT;
907 const point& cutPt = cutPoints[cutF[fp]];
909 pointHit pHit = masterF.nearestPoint(cutPt, masterPoints);
911 maxDist = max(maxDist, pHit.distance());
917 void Foam::faceCoupleInfo::findPerfectMatchingFaces
919 const primitiveMesh& mesh0,
920 const primitiveMesh& mesh1,
923 labelList& mesh0Faces,
924 labelList& mesh1Faces
927 // Face centres of external faces (without invoking
928 // mesh.faceCentres since mesh might have been clearedOut)
932 calcFaceCentres<List>
936 mesh0.nInternalFaces(),
937 mesh0.nFaces() - mesh0.nInternalFaces()
943 calcFaceCentres<List>
947 mesh1.nInternalFaces(),
948 mesh1.nFaces() - mesh1.nInternalFaces()
955 Pout<< "Face matching tolerance : " << absTol << endl;
959 // Match geometrically
961 bool matchedAllFaces = matchPoints
965 scalarField(fc1.size(), absTol),
973 << "faceCoupleInfo::faceCoupleInfo : "
974 << "Matched ALL " << fc1.size()
975 << " boundary faces of mesh0 to boundary faces of mesh1." << endl
976 << "This is only valid if the mesh to add is fully"
977 << " enclosed by the mesh it is added to." << endl;
984 mesh0Faces.setSize(fc0.size());
985 mesh1Faces.setSize(fc1.size());
989 if (from1To0[i] != -1)
991 mesh1Faces[nMatched] = i + mesh1.nInternalFaces();
992 mesh0Faces[nMatched] = from1To0[i] + mesh0.nInternalFaces();
998 mesh0Faces.setSize(nMatched);
999 mesh1Faces.setSize(nMatched);
1003 void Foam::faceCoupleInfo::findSlavesCoveringMaster
1005 const primitiveMesh& mesh0,
1006 const primitiveMesh& mesh1,
1007 const scalar absTol,
1009 labelList& mesh0Faces,
1010 labelList& mesh1Faces
1013 // Construct octree from all mesh0 boundary faces
1014 octreeDataFace shapes(mesh0);
1016 treeBoundBox overallBb(mesh0.points());
1018 octree<octreeDataFace> tree
1020 overallBb, // overall search domain
1021 shapes, // all information needed to do checks on cells
1023 20.0, // maximum ratio of cubes v.s. cells
1030 Pout<< "findSlavesCoveringMaster :"
1031 << " constructed octree for mesh0 boundary faces" << endl;
1034 // Search nearest face for every mesh1 boundary face.
1036 labelHashSet mesh0Set(mesh0.nFaces() - mesh0.nInternalFaces());
1037 labelHashSet mesh1Set(mesh1.nFaces() - mesh1.nInternalFaces());
1041 label mesh1FaceI = mesh1.nInternalFaces();
1042 mesh1FaceI < mesh1.nFaces();
1046 const face& f1 = mesh1.faces()[mesh1FaceI];
1048 // Generate face centre (prevent cellCentres() reconstruction)
1049 point fc(f1.centre(mesh1.points()));
1051 // Search in bounding box of face only.
1052 treeBoundBox tightest(static_cast<const pointField&>(f1.points(mesh1.points())));
1054 scalar tightestDist = GREAT;
1056 label index = tree.findNearest(fc, tightest, tightestDist);
1061 label mesh0FaceI = shapes.meshFaces()[index];
1063 // Check if points of f1 actually lie on top of mesh0 face
1064 // This is the bit that might fail since if f0 is severely warped
1065 // and f1's points are not present on f0 (since subdivision)
1066 // then the distance of the points to f0 might be quite large
1067 // and the test will fail. This all should in fact be some kind
1068 // of walk from known corresponding points/edges/faces.
1074 mesh0.faces()[mesh0FaceI],
1080 mesh0Set.insert(mesh0FaceI);
1081 mesh1Set.insert(mesh1FaceI);
1088 Pout<< "findSlavesCoveringMaster :"
1089 << " matched " << mesh1Set.size() << " mesh1 faces to "
1090 << mesh0Set.size() << " mesh0 faces" << endl;
1093 mesh0Faces = mesh0Set.toc();
1094 mesh1Faces = mesh1Set.toc();
1098 // Grow cutToMasterFace across 'internal' edges.
1099 Foam::label Foam::faceCoupleInfo::growCutFaces
1101 const labelList& cutToMasterEdges,
1102 Map<labelList>& candidates
1107 Pout<< "growCutFaces :"
1108 << " growing cut faces to masterPatch" << endl;
1111 label nTotChanged = 0;
1115 const labelListList& cutFaceEdges = cutFaces().faceEdges();
1116 const labelListList& cutEdgeFaces = cutFaces().edgeFaces();
1120 forAll(cutToMasterFaces_, cutFaceI)
1122 const label masterFaceI = cutToMasterFaces_[cutFaceI];
1124 if (masterFaceI != -1)
1126 // We now have a cutFace for which we have already found a
1127 // master face. Grow this masterface across any internal edge
1128 // (internal: no corresponding master edge)
1130 const labelList& fEdges = cutFaceEdges[cutFaceI];
1134 const label cutEdgeI = fEdges[i];
1136 if (cutToMasterEdges[cutEdgeI] == -1)
1139 // - internal to the cutPatch (since all region edges
1141 // - on cutFaceI which corresponds to masterFace.
1142 // Mark all connected faces with this masterFace.
1144 const labelList& eFaces = cutEdgeFaces[cutEdgeI];
1148 const label faceI = eFaces[j];
1150 if (cutToMasterFaces_[faceI] == -1)
1152 cutToMasterFaces_[faceI] = masterFaceI;
1153 candidates.erase(faceI);
1156 else if (cutToMasterFaces_[faceI] != masterFaceI)
1158 const pointField& cutPoints =
1159 cutFaces().points();
1160 const pointField& masterPoints =
1161 masterPatch().points();
1163 const edge& e = cutFaces().edges()[cutEdgeI];
1165 label myMaster = cutToMasterFaces_[faceI];
1166 const face& myF = masterPatch()[myMaster];
1168 const face& nbrF = masterPatch()[masterFaceI];
1172 "faceCoupleInfo::growCutFaces"
1173 "(const labelList&, Map<labelList>&)"
1175 << cutFaces()[faceI].points(cutPoints)
1178 << " but also connects to nbr face "
1179 << cutFaces()[cutFaceI].points(cutPoints)
1180 << " with master " << masterFaceI
1183 << myF.points(masterPoints)
1184 << " nbrMasterFace:"
1185 << nbrF.points(masterPoints) << nl
1186 << "Across cut edge " << cutEdgeI
1188 << cutFaces().localPoints()[e[0]]
1189 << cutFaces().localPoints()[e[1]]
1190 << abort(FatalError);
1200 Pout<< "growCutFaces : Grown an additional " << nChanged
1201 << " cut to master face" << " correspondence" << endl;
1204 nTotChanged += nChanged;
1216 void Foam::faceCoupleInfo::checkMatch(const labelList& cutToMasterEdges) const
1218 const pointField& cutLocalPoints = cutFaces().localPoints();
1220 const pointField& masterLocalPoints = masterPatch().localPoints();
1221 const faceList& masterLocalFaces = masterPatch().localFaces();
1223 forAll(cutToMasterEdges, cutEdgeI)
1225 const edge& e = cutFaces().edges()[cutEdgeI];
1227 if (cutToMasterEdges[cutEdgeI] == -1)
1229 // Internal edge. Check that master face is same on both sides.
1230 const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1232 label masterFaceI = -1;
1234 forAll(cutEFaces, i)
1236 label cutFaceI = cutEFaces[i];
1238 if (cutToMasterFaces_[cutFaceI] != -1)
1240 if (masterFaceI == -1)
1242 masterFaceI = cutToMasterFaces_[cutFaceI];
1244 else if (masterFaceI != cutToMasterFaces_[cutFaceI])
1246 label myMaster = cutToMasterFaces_[cutFaceI];
1247 const face& myF = masterLocalFaces[myMaster];
1249 const face& nbrF = masterLocalFaces[masterFaceI];
1253 "faceCoupleInfo::checkMatch(const labelList&) const"
1255 << "Internal CutEdge " << e
1257 << cutLocalPoints[e[0]]
1258 << cutLocalPoints[e[1]]
1259 << " connects to master " << myMaster
1260 << " and to master " << masterFaceI << nl
1262 << myF.points(masterLocalPoints)
1263 << " nbrMasterFace:"
1264 << nbrF.points(masterLocalPoints)
1265 << abort(FatalError);
1274 // Extends matching information by elimination across cutFaces using more
1275 // than one region edge. Updates cutToMasterFaces_ and sets candidates
1276 // which is for every cutface on a region edge the possible master faces.
1277 Foam::label Foam::faceCoupleInfo::matchEdgeFaces
1279 const labelList& cutToMasterEdges,
1280 Map<labelList>& candidates
1283 // For every unassigned cutFaceI the possible list of master faces.
1285 candidates.resize(cutFaces().size());
1289 forAll(cutToMasterEdges, cutEdgeI)
1291 label masterEdgeI = cutToMasterEdges[cutEdgeI];
1293 if (masterEdgeI != -1)
1295 // So cutEdgeI is matched to masterEdgeI. For all cut faces using
1296 // the cut edge store the master face as a candidate match.
1298 const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1299 const labelList& masterEFaces =
1300 masterPatch().edgeFaces()[masterEdgeI];
1302 forAll(cutEFaces, i)
1304 label cutFaceI = cutEFaces[i];
1306 if (cutToMasterFaces_[cutFaceI] == -1)
1308 // So this face (cutFaceI) is on an edge (cutEdgeI) that
1309 // is used by some master faces (masterEFaces). Check if
1310 // the cutFaceI and some of these masterEFaces share more
1311 // than one edge (which uniquely defines face)
1313 // Combine master faces with current set of candidate
1315 Map<labelList>::iterator fnd = candidates.find(cutFaceI);
1317 if (fnd == candidates.end())
1319 // No info yet for cutFaceI. Add all master faces as
1321 candidates.insert(cutFaceI, masterEFaces);
1325 // From some other cutEdgeI there are already some
1326 // candidate master faces. Check the overlap with
1327 // the current set of master faces.
1328 const labelList& masterFaces = fnd();
1330 DynamicList<label> newCandidates(masterFaces.size());
1332 forAll(masterEFaces, j)
1334 if (findIndex(masterFaces, masterEFaces[j]) != -1)
1336 newCandidates.append(masterEFaces[j]);
1340 if (newCandidates.size() == 1)
1342 // We found a perfect match. Delete entry from
1344 cutToMasterFaces_[cutFaceI] = newCandidates[0];
1345 candidates.erase(cutFaceI);
1350 // Should not happen?
1351 //Pout<< "On edge:" << cutEdgeI
1352 // << " have connected masterFaces:"
1354 // << " and from previous edge we have"
1355 // << " connected masterFaces" << masterFaces
1356 // << " . Overlap:" << newCandidates << endl;
1358 fnd() = newCandidates.shrink();
1369 Pout<< "matchEdgeFaces : Found " << nChanged
1370 << " faces where there was"
1371 << " only one remaining choice for cut-master correspondence"
1379 // Gets a list of cutFaces (that use a master edge) and the candidate
1381 // Finds most aligned master face.
1382 Foam::label Foam::faceCoupleInfo::geometricMatchEdgeFaces
1384 Map<labelList>& candidates
1387 const pointField& cutPoints = cutFaces().points();
1391 // Mark all master faces that have been matched so far.
1393 labelListList masterToCutFaces
1397 masterPatch().size(),
1402 forAllConstIter(Map<labelList>, candidates, iter)
1404 label cutFaceI = iter.key();
1406 const face& cutF = cutFaces()[cutFaceI];
1408 if (cutToMasterFaces_[cutFaceI] == -1)
1410 const labelList& masterFaces = iter();
1412 // Find the best matching master face.
1413 scalar minDist = GREAT;
1414 label minMasterFaceI = -1;
1416 forAll(masterFaces, i)
1418 label masterFaceI = masterFaces[i];
1420 if (masterToCutFaces[masterFaces[i]].size() == 0)
1427 masterPatch()[masterFaceI],
1428 masterPatch().points()
1434 minMasterFaceI = masterFaceI;
1439 if (minMasterFaceI != -1)
1441 cutToMasterFaces_[cutFaceI] = minMasterFaceI;
1442 masterToCutFaces[minMasterFaceI] = cutFaceI;
1448 // (inefficiently) force candidates to be uptodate.
1449 forAll(cutToMasterFaces_, cutFaceI)
1451 if (cutToMasterFaces_[cutFaceI] != -1)
1453 candidates.erase(cutFaceI);
1460 Pout<< "geometricMatchEdgeFaces : Found " << nChanged
1461 << " faces where there was"
1462 << " only one remaining choice for cut-master correspondence"
1470 // Calculate the set of cut faces inbetween master and slave patch
1471 // assuming perfect match (and optional face ordering on slave)
1472 void Foam::faceCoupleInfo::perfectPointMatch
1474 const scalar absTol,
1475 const bool slaveFacesOrdered
1480 Pout<< "perfectPointMatch :"
1481 << " Matching master and slave to cut."
1482 << " Master and slave faces are identical" << nl;
1484 if (slaveFacesOrdered)
1486 Pout<< "and master and slave faces are ordered"
1487 << " (on coupled patches)" << endl;
1491 Pout<< "and master and slave faces are not ordered"
1492 << " (on coupled patches)" << endl;
1496 cutToMasterFaces_ = identity(masterPatch().size());
1497 cutPoints_ = masterPatch().localPoints();
1500 new primitiveFacePatch
1502 masterPatch().localFaces(),
1506 masterToCutPoints_ = identity(cutPoints_.size());
1509 // Cut faces to slave patch.
1510 bool matchedAllFaces = false;
1512 if (slaveFacesOrdered)
1514 cutToSlaveFaces_ = identity(cutFaces().size());
1515 matchedAllFaces = (cutFaces().size() == slavePatch().size());
1519 // Faces do not have to be ordered (but all have
1520 // to match). Note: Faces will be already ordered if we enter here from
1521 // construct from meshes.
1522 matchedAllFaces = matchPoints
1524 calcFaceCentres<List>
1531 calcFaceCentres<IndirectList>
1534 slavePatch().points(),
1538 scalarField(slavePatch().size(), absTol),
1545 if (!matchedAllFaces)
1549 "faceCoupleInfo::perfectPointMatch"
1550 "(const scalar&, const bool)"
1551 ) << "Did not match all of the master faces to the slave faces"
1553 << "This usually means that the slave patch and master patch"
1554 << " do not align to within " << absTol << " meter."
1555 << abort(FatalError);
1559 // Find correspondence from slave points to cut points. This might
1560 // detect shared points so the output is a slave-to-cut point list
1561 // and a compaction list.
1563 labelList cutToCompact, compactToCut;
1564 matchPointsThroughFaces
1567 cutFaces().localPoints(),
1568 reorder(cutToSlaveFaces_, cutFaces().localFaces()),
1569 slavePatch().localPoints(),
1570 slavePatch().localFaces(),
1571 false, // slave and cut have opposite orientation
1573 slaveToCutPoints_, // slave to (uncompacted) cut points
1574 cutToCompact, // compaction map: from cut to compacted
1575 compactToCut // compaction map: from compacted to cut
1579 // Use compaction lists to renumber cutPoints.
1580 cutPoints_ = IndirectList<point>(cutPoints_, compactToCut)();
1582 const faceList& cutLocalFaces = cutFaces().localFaces();
1584 faceList compactFaces(cutLocalFaces.size());
1585 forAll(cutLocalFaces, i)
1587 compactFaces[i] = renumber(cutToCompact, cutLocalFaces[i]);
1591 new primitiveFacePatch
1598 inplaceRenumber(cutToCompact, slaveToCutPoints_);
1599 inplaceRenumber(cutToCompact, masterToCutPoints_);
1603 // Calculate the set of cut faces inbetween master and slave patch
1604 // assuming that slave patch is subdivision of masterPatch.
1605 void Foam::faceCoupleInfo::subDivisionMatch
1607 const polyMesh& slaveMesh,
1608 const bool patchDivision,
1614 Pout<< "subDivisionMatch :"
1615 << " Matching master and slave to cut."
1616 << " Slave can be subdivision of master but all master edges"
1617 << " have to be covered by slave edges." << endl;
1620 // Assume slave patch is subdivision of the master patch so cutFaces is a
1621 // copy of the slave (but reversed since orientation has to be like master).
1622 cutPoints_ = slavePatch().localPoints();
1624 faceList cutFaces(slavePatch().size());
1628 cutFaces[i] = slavePatch().localFaces()[i].reverseFace();
1630 cutFacesPtr_.reset(new primitiveFacePatch(cutFaces, cutPoints_));
1633 // Cut is copy of slave so addressing to slave is trivial.
1634 cutToSlaveFaces_ = identity(cutFaces().size());
1635 slaveToCutPoints_ = identity(slavePatch().nPoints());
1637 // Cut edges to slave patch
1638 labelList cutToSlaveEdges
1643 slaveToCutPoints_, //note:should be cutToSlavePoints but since iden
1651 OFstream str("regionEdges.obj");
1653 Pout<< "subDivisionMatch :"
1654 << " addressing for slave patch fully done."
1655 << " Dumping region edges to " << str.name() << endl;
1659 forAll(slavePatch().edges(), slaveEdgeI)
1661 if (regionEdge(slaveMesh, slaveEdgeI))
1663 const edge& e = slavePatch().edges()[slaveEdgeI];
1665 meshTools::writeOBJ(str, slavePatch().localPoints()[e[0]]);
1667 meshTools::writeOBJ(str, slavePatch().localPoints()[e[1]]);
1669 str<< "l " << vertI-1 << ' ' << vertI << nl;
1675 // Addressing from cut to master.
1677 // Cut points to master patch
1678 // - determine master points to cut points. Has to be full!
1679 // - invert to get cut to master
1682 Pout<< "subDivisionMatch :"
1683 << " matching master points to cut points" << endl;
1686 bool matchedAllPoints = matchPoints
1688 masterPatch().localPoints(),
1690 scalarField(masterPatch().nPoints(), absTol),
1695 if (!matchedAllPoints)
1699 "faceCoupleInfo::subDivisionMatch"
1700 "(const polyMesh&, const bool, const scalar)"
1701 ) << "Did not match all of the master points to the slave points"
1703 << "This usually means that the slave patch is not a"
1704 << " subdivision of the master patch"
1705 << abort(FatalError);
1709 // Do masterEdges to cutEdges :
1710 // - mark all edges between two masterEdge endpoints. (geometric test since
1711 // this is the only distinction)
1712 // - this gives the borders inbetween which all cutfaces come from
1713 // a single master face.
1716 Pout<< "subDivisionMatch :"
1717 << " matching cut edges to master edges" << endl;
1720 const edgeList& masterEdges = masterPatch().edges();
1721 const pointField& masterPoints = masterPatch().localPoints();
1723 const edgeList& cutEdges = cutFaces().edges();
1725 labelList cutToMasterEdges(cutFaces().nEdges(), -1);
1727 forAll(masterEdges, masterEdgeI)
1729 const edge& masterEdge = masterEdges[masterEdgeI];
1731 label cutPoint0 = masterToCutPoints_[masterEdge[0]];
1732 label cutPoint1 = masterToCutPoints_[masterEdge[1]];
1734 // Find edges between cutPoint0 and cutPoint1.
1736 label cutPointI = cutPoint0;
1739 // Find edge (starting at pointI on cut), aligned with master
1756 // Problem. Report while matching to produce nice error message
1769 Pout<< "Dumping unmatched pointEdges to errorEdges.obj"
1778 cutFaces().pointEdges()[cutPointI]
1780 cutFaces().localPoints(),
1786 "faceCoupleInfo::subDivisionMatch"
1787 "(const polyMesh&, const bool, const scalar)"
1788 ) << "Problem in finding cut edges corresponding to"
1789 << " master edge " << masterEdge
1790 << " points:" << masterPoints[masterEdge[0]]
1791 << ' ' << masterPoints[masterEdge[1]]
1792 << " corresponding cut edge: (" << cutPoint0
1794 << abort(FatalError);
1797 cutToMasterEdges[cutEdgeI] = masterEdgeI;
1799 cutPointI = cutEdges[cutEdgeI].otherVertex(cutPointI);
1801 } while(cutPointI != cutPoint1);
1806 // Write all matched edges.
1807 writeEdges(cutToMasterEdges, cutToSlaveEdges);
1810 // Rework cutToMasterEdges into list of points inbetween two endpoints
1811 // (cutEdgeToPoints_)
1812 setCutEdgeToPoints(cutToMasterEdges);
1815 // Now that we have marked the cut edges that lie on top of master edges
1816 // we can find cut faces that have two (or more) of these edges.
1817 // Note that there can be multiple faces having two or more matched edges
1818 // since the cut faces can form a non-manifold surface(?)
1819 // So the matching is done as an elimination process where for every
1820 // cutFace on a matched edge we store the possible master faces and
1821 // eliminate more and more until we only have one possible master face
1826 Pout<< "subDivisionMatch :"
1827 << " matching (topological) cut faces to masterPatch" << endl;
1830 // For every unassigned cutFaceI the possible list of master faces.
1831 Map<labelList> candidates(cutFaces().size());
1833 cutToMasterFaces_.setSize(cutFaces().size());
1834 cutToMasterFaces_ = -1;
1838 // See if there are any edges left where there is only one remaining
1840 label nChanged = matchEdgeFaces(cutToMasterEdges, candidates);
1842 checkMatch(cutToMasterEdges);
1845 // Now we should have found a face correspondence for every cutFace
1846 // that uses two or more cutEdges. Floodfill from these across
1847 // 'internal' cutedges (i.e. edges that do not have a master
1848 // equivalent) to determine some more cutToMasterFaces
1849 nChanged += growCutFaces(cutToMasterEdges, candidates);
1851 checkMatch(cutToMasterEdges);
1861 Pout<< "subDivisionMatch :"
1862 << " matching (geometric) cut faces to masterPatch" << endl;
1865 // Above should have matched all cases fully. Cannot prove this yet so
1866 // do any remaining unmatched edges by a geometric test.
1869 // See if there are any edges left where there is only one remaining
1871 label nChanged = geometricMatchEdgeFaces(candidates);
1873 checkMatch(cutToMasterEdges);
1875 nChanged += growCutFaces(cutToMasterEdges, candidates);
1877 checkMatch(cutToMasterEdges);
1886 // All cut faces matched?
1887 forAll(cutToMasterFaces_, cutFaceI)
1889 if (cutToMasterFaces_[cutFaceI] == -1)
1891 const face& cutF = cutFaces()[cutFaceI];
1895 "faceCoupleInfo::subDivisionMatch"
1896 "(const polyMesh&, const bool, const scalar)"
1897 ) << "Did not match all of cutFaces to a master face" << nl
1898 << "First unmatched cut face:" << cutFaceI << " with points:"
1899 << IndirectList<point>(cutFaces().points(), cutF)()
1901 << "This usually means that the slave patch is not a"
1902 << " subdivision of the master patch"
1903 << abort(FatalError);
1909 Pout<< "subDivisionMatch :"
1910 << " finished matching master and slave to cut" << endl;
1920 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1922 // Construct from mesh data
1923 Foam::faceCoupleInfo::faceCoupleInfo
1925 const polyMesh& masterMesh,
1926 const polyMesh& slaveMesh,
1927 const scalar absTol,
1928 const bool perfectMatch
1931 masterPatchPtr_(NULL),
1932 slavePatchPtr_(NULL),
1935 cutToMasterFaces_(0),
1936 masterToCutPoints_(0),
1937 cutToSlaveFaces_(0),
1938 slaveToCutPoints_(0),
1941 // Get faces on both meshes that are aligned.
1942 // (not ordered i.e. masterToMesh[0] does
1943 // not couple to slaveToMesh[0]. This ordering is done later on)
1944 labelList masterToMesh;
1945 labelList slaveToMesh;
1949 // Identical faces so use tight face-centre comparison.
1950 findPerfectMatchingFaces
1961 // Slave subdivision of master so use 'nearest'. Bit dodgy, especially
1962 // with using absTol (which is quite small)
1963 findSlavesCoveringMaster
1973 // Construct addressing engines for both sides
1974 masterPatchPtr_.reset
1976 new indirectPrimitivePatch
1978 IndirectList<face>(masterMesh.faces(), masterToMesh),
1983 slavePatchPtr_.reset
1985 new indirectPrimitivePatch
1987 IndirectList<face>(slaveMesh.faces(), slaveToMesh),
1995 // Faces are perfectly aligned but probably not ordered.
1996 perfectPointMatch(absTol, false);
2000 // Slave faces are subdivision of master face. Faces not ordered.
2001 subDivisionMatch(slaveMesh, false, absTol);
2011 // Slave is subdivision of master patch.
2012 // (so -both cover the same area -all of master points are present in slave)
2013 Foam::faceCoupleInfo::faceCoupleInfo
2015 const polyMesh& masterMesh,
2016 const labelList& masterAddressing,
2017 const polyMesh& slaveMesh,
2018 const labelList& slaveAddressing,
2019 const scalar absTol,
2020 const bool perfectMatch,
2021 const bool orderedFaces,
2022 const bool patchDivision
2027 new indirectPrimitivePatch
2029 IndirectList<face>(masterMesh.faces(), masterAddressing),
2035 new indirectPrimitivePatch
2037 IndirectList<face>(slaveMesh.faces(), slaveAddressing),
2043 cutToMasterFaces_(0),
2044 masterToCutPoints_(0),
2045 cutToSlaveFaces_(0),
2046 slaveToCutPoints_(0),
2049 if (perfectMatch && (masterAddressing.size() != slaveAddressing.size()))
2053 "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2054 ", const primitiveMesh&, const scalar, const bool"
2055 ) << "Perfect match specified but number of master and slave faces"
2056 << " differ." << endl
2057 << "master:" << masterAddressing.size()
2058 << " slave:" << slaveAddressing.size()
2059 << abort(FatalError);
2064 masterAddressing.size() > 0
2065 && min(masterAddressing) < masterMesh.nInternalFaces()
2070 "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2071 ", const primitiveMesh&, const scalar, const bool"
2072 ) << "Supplied internal face on master mesh to couple." << nl
2073 << "Faces to be coupled have to be boundary faces."
2074 << abort(FatalError);
2078 slaveAddressing.size() > 0
2079 && min(slaveAddressing) < slaveMesh.nInternalFaces()
2084 "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2085 ", const primitiveMesh&, const scalar, const bool"
2086 ) << "Supplied internal face on slave mesh to couple." << nl
2087 << "Faces to be coupled have to be boundary faces."
2088 << abort(FatalError);
2094 perfectPointMatch(absTol, orderedFaces);
2098 // Slave faces are subdivision of master face. Faces not ordered.
2099 subDivisionMatch(slaveMesh, patchDivision, absTol);
2109 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
2111 Foam::faceCoupleInfo::~faceCoupleInfo()
2115 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2117 Foam::labelList Foam::faceCoupleInfo::faceLabels(const polyPatch& pp)
2119 labelList faces(pp.size());
2121 label faceI = pp.start();
2131 Foam::Map<Foam::label> Foam::faceCoupleInfo::makeMap(const labelList& lst)
2133 Map<label> map(lst.size());
2139 map.insert(i, lst[i]);
2146 Foam::Map<Foam::labelList> Foam::faceCoupleInfo::makeMap
2148 const labelListList& lst
2151 Map<labelList> map(lst.size());
2155 if (lst[i].size() > 0)
2157 map.insert(i, lst[i]);
2164 // ************************************************************************* //