1 /*---------------------------------------------------------------------------* \
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 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"
28 #include <OpenFOAM/polyMesh.H>
29 #include <OpenFOAM/matchPoints.H>
30 #include <OpenFOAM/indirectPrimitivePatch.H>
31 #include <meshTools/meshTools.H>
32 #include <meshTools/treeDataFace.H>
33 #include <meshTools/indexedOctree.H>
34 #include <OpenFOAM/OFstream.H>
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 defineTypeNameAndDebug(Foam::faceCoupleInfo, 0);
40 const Foam::scalar Foam::faceCoupleInfo::angleTol_ = 1E-3;
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 void Foam::faceCoupleInfo::writeOBJ
48 const fileName& fName,
49 const edgeList& edges,
50 const pointField& points,
56 labelList pointMap(points.size(), -1);
64 const edge& e = edges[edgeI];
70 if (pointMap[pointI] == -1)
72 pointMap[pointI] = newPointI++;
74 meshTools::writeOBJ(str, points[pointI]);
81 forAll(points, pointI)
83 meshTools::writeOBJ(str, points[pointI]);
86 pointMap = identity(points.size());
91 const edge& e = edges[edgeI];
93 str<< "l " << pointMap[e[0]]+1 << ' ' << pointMap[e[1]]+1 << nl;
99 void Foam::faceCoupleInfo::writeOBJ
101 const fileName& fName,
102 const pointField& points0,
103 const pointField& points1
106 Pout<< "Writing connections as edges to " << fName << endl;
114 meshTools::writeOBJ(str, points0[i]);
116 meshTools::writeOBJ(str, points1[i]);
118 str << "l " << vertI-1 << ' ' << vertI << nl;
123 //- Writes face and point connectivity as .obj files.
124 void Foam::faceCoupleInfo::writePointsFaces() const
126 const indirectPrimitivePatch& m = masterPatch();
127 const indirectPrimitivePatch& s = slavePatch();
128 const primitiveFacePatch& c = cutFaces();
132 OFstream str("masterPatch.obj");
133 Pout<< "Writing masterPatch to " << str.name() << endl;
134 meshTools::writeOBJ(str, m.localFaces(), m.localPoints());
137 OFstream str("slavePatch.obj");
138 Pout<< "Writing slavePatch to " << str.name() << endl;
139 meshTools::writeOBJ(str, s.localFaces(), s.localPoints());
142 OFstream str("cutFaces.obj");
143 Pout<< "Writing cutFaces to " << str.name() << endl;
144 meshTools::writeOBJ(str, c.localFaces(), c.localPoints());
147 // Point connectivity
149 Pout<< "Writing cutToMasterPoints to cutToMasterPoints.obj" << endl;
153 "cutToMasterPoints.obj",
155 pointField(c.localPoints(), masterToCutPoints_));
158 Pout<< "Writing cutToSlavePoints to cutToSlavePoints.obj" << endl;
162 "cutToSlavePoints.obj",
164 pointField(c.localPoints(), slaveToCutPoints_)
170 Pout<< "Writing cutToMasterFaces to cutToMasterFaces.obj" << endl;
172 pointField equivMasterFaces(c.size());
174 forAll(cutToMasterFaces(), cutFaceI)
176 label masterFaceI = cutToMasterFaces()[cutFaceI];
178 if (masterFaceI != -1)
180 equivMasterFaces[cutFaceI] = m[masterFaceI].centre(m.points());
184 WarningIn("writePointsFaces()")
185 << "No master face for cut face " << cutFaceI
186 << " at position " << c[cutFaceI].centre(c.points())
189 equivMasterFaces[cutFaceI] = vector::zero;
195 "cutToMasterFaces.obj",
196 calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
202 Pout<< "Writing cutToSlaveFaces to cutToSlaveFaces.obj" << endl;
204 pointField equivSlaveFaces(c.size());
206 forAll(cutToSlaveFaces(), cutFaceI)
208 label slaveFaceI = cutToSlaveFaces()[cutFaceI];
210 equivSlaveFaces[cutFaceI] = s[slaveFaceI].centre(s.points());
215 "cutToSlaveFaces.obj",
216 calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
225 void Foam::faceCoupleInfo::writeEdges
227 const labelList& cutToMasterEdges,
228 const labelList& cutToSlaveEdges
231 const indirectPrimitivePatch& m = masterPatch();
232 const indirectPrimitivePatch& s = slavePatch();
233 const primitiveFacePatch& c = cutFaces();
237 OFstream str("cutToMasterEdges.obj");
238 Pout<< "Writing cutToMasterEdges to " << str.name() << endl;
242 forAll(cutToMasterEdges, cutEdgeI)
244 if (cutToMasterEdges[cutEdgeI] != -1)
246 const edge& masterEdge =
247 m.edges()[cutToMasterEdges[cutEdgeI]];
248 const edge& cutEdge = c.edges()[cutEdgeI];
250 meshTools::writeOBJ(str, m.localPoints()[masterEdge[0]]);
252 meshTools::writeOBJ(str, m.localPoints()[masterEdge[1]]);
254 meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
256 meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
258 str << "l " << vertI-3 << ' ' << vertI-2 << nl;
259 str << "l " << vertI-3 << ' ' << vertI-1 << nl;
260 str << "l " << vertI-3 << ' ' << vertI << nl;
261 str << "l " << vertI-2 << ' ' << vertI-1 << nl;
262 str << "l " << vertI-2 << ' ' << vertI << nl;
263 str << "l " << vertI-1 << ' ' << vertI << nl;
268 OFstream str("cutToSlaveEdges.obj");
269 Pout<< "Writing cutToSlaveEdges to " << str.name() << endl;
273 labelList slaveToCut(invert(s.nEdges(), cutToSlaveEdges));
275 forAll(slaveToCut, edgeI)
277 if (slaveToCut[edgeI] != -1)
279 const edge& slaveEdge = s.edges()[edgeI];
280 const edge& cutEdge = c.edges()[slaveToCut[edgeI]];
282 meshTools::writeOBJ(str, s.localPoints()[slaveEdge[0]]);
284 meshTools::writeOBJ(str, s.localPoints()[slaveEdge[1]]);
286 meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
288 meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
290 str << "l " << vertI-3 << ' ' << vertI-2 << nl;
291 str << "l " << vertI-3 << ' ' << vertI-1 << nl;
292 str << "l " << vertI-3 << ' ' << vertI << nl;
293 str << "l " << vertI-2 << ' ' << vertI-1 << nl;
294 str << "l " << vertI-2 << ' ' << vertI << nl;
295 str << "l " << vertI-1 << ' ' << vertI << nl;
304 // Given an edgelist and a map for the points on the edges it tries to find
305 // the corresponding patch edges.
306 Foam::labelList Foam::faceCoupleInfo::findMappedEdges
308 const edgeList& edges,
309 const labelList& pointMap,
310 const indirectPrimitivePatch& patch
313 labelList toPatchEdges(edges.size());
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 // Detect a cut edge which originates from two boundary faces having different
337 bool Foam::faceCoupleInfo::regionEdge
339 const polyMesh& slaveMesh,
340 const label slaveEdgeI
343 const labelList& eFaces = slavePatch().edgeFaces()[slaveEdgeI];
345 if (eFaces.size() == 1)
351 // Count how many different patches connected to this edge.
357 label faceI = eFaces[i];
359 label meshFaceI = slavePatch().addressing()[faceI];
361 label patchI = slaveMesh.boundaryMesh().whichPatch(meshFaceI);
367 else if (patchI != patch0)
369 // Found two different patches connected to this edge.
378 // Find edge using pointI that is most aligned with vector between
379 // master points. Patchdivision tells us whether or not to use
380 // patch information to match edges.
381 Foam::label Foam::faceCoupleInfo::mostAlignedCutEdge
384 const polyMesh& slaveMesh,
385 const bool patchDivision,
386 const labelList& cutToMasterEdges,
387 const labelList& cutToSlaveEdges,
389 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]
410 // Find the edge that gets us nearest end.
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)
437 // Shortcut: found edge end point.
440 Pout<< " mostAlignedEdge : found end point " << edgeEnd
446 // Get angle between edge and edge to masterEnd
448 vector eVec(localPoints[otherPointI] - localPoints[pointI]);
450 scalar magEVec = mag(eVec);
452 if (magEVec < VSMALL)
454 WarningIn("faceCoupleInfo::mostAlignedEdge")
455 << "Crossing zero sized edge " << edgeI
456 << " coords:" << localPoints[otherPointI]
457 << localPoints[pointI]
458 << " when walking from " << localPoints[edgeStart]
459 << " to " << localPoints[edgeEnd]
466 vector eToEndPoint(localPoints[edgeEnd] - localPoints[otherPointI]);
467 eToEndPoint /= mag(eToEndPoint);
469 scalar cosAngle = eVec & eToEndPoint;
473 Pout<< " edge:" << e << " points:" << localPoints[pointI]
474 << localPoints[otherPointI]
476 << " vecToEnd:" << eToEndPoint
477 << " cosAngle:" << cosAngle
481 if (cosAngle > maxCos)
489 if (maxCos > 1 - angleTol_)
500 // Construct points to split points map (in cut addressing)
501 void Foam::faceCoupleInfo::setCutEdgeToPoints(const labelList& cutToMasterEdges)
503 labelListList masterToCutEdges
507 masterPatch().nEdges(),
512 const edgeList& cutEdges = cutFaces().edges();
514 // Size extra big so searching is faster
515 cutEdgeToPoints_.resize
517 masterPatch().nEdges()
518 + slavePatch().nEdges()
522 forAll(masterToCutEdges, masterEdgeI)
524 const edge& masterE = masterPatch().edges()[masterEdgeI];
526 //Pout<< "Master:" << masterPatch().localPoints()[masterE[0]] << ' '
527 // << masterPatch().localPoints()[masterE[1]] << endl;
529 const labelList& stringedEdges = masterToCutEdges[masterEdgeI];
531 if (stringedEdges.empty())
535 "faceCoupleInfo::setCutEdgeToPoints"
537 ) << "Did not match all of master edges to cutFace edges"
539 << "First unmatched edge:" << masterEdgeI << " endPoints:"
540 << masterPatch().localPoints()[masterE[0]]
541 << masterPatch().localPoints()[masterE[1]]
543 << "This usually means that the slave patch is not a"
544 << " subdivision of the master patch"
545 << abort(FatalError);
547 else if (stringedEdges.size() > 1)
549 // String up the edges between e[0] and e[1]. Store the points
550 // inbetween e[0] and e[1] (all in cutFaces() labels)
552 DynamicList<label> splitPoints(stringedEdges.size()-1);
554 // Unsplit edge endpoints
555 const edge unsplitEdge
557 masterToCutPoints_[masterE[0]],
558 masterToCutPoints_[masterE[1]]
561 label startVertI = unsplitEdge[0];
562 label startEdgeI = -1;
564 while (startVertI != unsplitEdge[1])
566 // Loop over all string of edges. Update
567 // - startVertI : previous vertex
568 // - startEdgeI : previous edge
569 // and insert any points into splitPoints
572 label oldStart = startVertI;
574 forAll(stringedEdges, i)
576 label edgeI = stringedEdges[i];
578 if (edgeI != startEdgeI)
580 const edge& e = cutEdges[edgeI];
582 //Pout<< " cut:" << e << " at:"
583 // << cutFaces().localPoints()[e[0]]
584 // << ' ' << cutFaces().localPoints()[e[1]] << endl;
586 if (e[0] == startVertI)
590 if (e[1] != unsplitEdge[1])
592 splitPoints.append(e[1]);
596 else if (e[1] == startVertI)
600 if (e[0] != unsplitEdge[1])
602 splitPoints.append(e[0]);
610 if (oldStart == startVertI)
614 "faceCoupleInfo::setCutEdgeToPoints"
616 ) << " unsplitEdge:" << unsplitEdge
617 << " does not correspond to split edges "
618 << UIndirectList<edge>(cutEdges, stringedEdges)()
619 << abort(FatalError);
623 //Pout<< "For master edge:"
625 // << " Found stringed points "
626 // << UIndirectList<point>
628 // cutFaces().localPoints(),
629 // splitPoints.shrink()
633 cutEdgeToPoints_.insert(unsplitEdge, splitPoints.shrink());
639 // Determines rotation for f1 to match up with f0, i.e. the index in f0 of
640 // the first point of f1.
641 Foam::label Foam::faceCoupleInfo::matchFaces
644 const pointField& points0,
646 const pointField& points1,
648 const bool sameOrientation
651 if (f0.size() != f1.size())
655 "faceCoupleInfo::matchFaces"
656 "(const scalar, const face&, const pointField&"
657 ", const face&, const pointField&)"
658 ) << "Different sizes for supposedly matching faces." << nl
659 << "f0:" << f0 << " coords:" << UIndirectList<point>(points0, f0)()
661 << "f1:" << f1 << " coords:" << UIndirectList<point>(points1, f1)()
662 << abort(FatalError);
665 const scalar absTolSqr = sqr(absTol);
672 // See -if starting from startFp on f0- the two faces match.
673 bool fullMatch = true;
680 scalar distSqr = Foam::magSqr(points0[f0[fp0]] - points1[f1[fp1]]);
682 if (distSqr > absTolSqr)
688 fp0 = f0.fcIndex(fp0); // walk forward
692 fp1 = f1.fcIndex(fp1);
696 fp1 = f1.rcIndex(fp1);
711 "faceCoupleInfo::matchFaces"
712 "(const scalar, const face&, const pointField&"
713 ", const face&, const pointField&)"
714 ) << "No unique match between two faces" << nl
715 << "Face " << f0 << " coords "
716 << UIndirectList<point>(points0, f0)() << nl
717 << "Face " << f1 << " coords "
718 << UIndirectList<point>(points1, f1)()
719 << "when using tolerance " << absTol
720 << " and forwardMatching:" << sameOrientation
721 << abort(FatalError);
728 // Find correspondence from patch points to cut points. This might
729 // detect shared points so the output is a patch-to-cut point list
730 // and a compaction list for the cut points (which will always be equal or more
731 // connected than the patch).
732 // Returns true if there are any duplicates.
733 bool Foam::faceCoupleInfo::matchPointsThroughFaces
736 const pointField& cutPoints,
737 const faceList& cutFaces,
738 const pointField& patchPoints,
739 const faceList& patchFaces,
740 const bool sameOrientation,
742 labelList& patchToCutPoints, // patch to (uncompacted) cut points
743 labelList& cutToCompact, // compaction list for cut points
744 labelList& compactToCut // inverse ,,
748 // From slave to cut point
749 patchToCutPoints.setSize(patchPoints.size());
750 patchToCutPoints = -1;
752 // Compaction list for cut points: either -1 or index into master which
753 // gives the point to compact to.
754 labelList cutPointRegion(cutPoints.size(), -1);
755 DynamicList<label> cutPointRegionMaster;
757 forAll(patchFaces, patchFaceI)
759 const face& patchF = patchFaces[patchFaceI];
761 //const face& cutF = cutFaces[patchToCutFaces[patchFaceI]];
762 const face& cutF = cutFaces[patchFaceI];
764 // Do geometric matching to get position of cutF[0] in patchF
765 label patchFp = matchFaces
772 sameOrientation // orientation
777 label cutPointI = cutF[cutFp];
778 label patchPointI = patchF[patchFp];
780 //const point& cutPt = cutPoints[cutPointI];
781 //const point& patchPt = patchPoints[patchPointI];
782 //if (mag(cutPt - patchPt) > SMALL)
784 // FatalErrorIn("matchPointsThroughFaces")
785 // << "cutP:" << cutPt
786 // << " patchP:" << patchPt
787 // << abort(FatalError);
790 if (patchToCutPoints[patchPointI] == -1)
792 patchToCutPoints[patchPointI] = cutPointI;
794 else if (patchToCutPoints[patchPointI] != cutPointI)
796 // Multiple cut points connecting to same patch.
797 // Check if already have region & region master for this set
798 label otherCutPointI = patchToCutPoints[patchPointI];
800 //Pout<< "PatchPoint:" << patchPt
801 // << " matches to:" << cutPointI
802 // << " coord:" << cutPoints[cutPointI]
803 // << " and to:" << otherCutPointI
804 // << " coord:" << cutPoints[otherCutPointI]
807 if (cutPointRegion[otherCutPointI] != -1)
809 // Have region for this set. Copy.
810 label region = cutPointRegion[otherCutPointI];
811 cutPointRegion[cutPointI] = region;
813 // Update region master with min point label
814 cutPointRegionMaster[region] = min
816 cutPointRegionMaster[region],
822 // Create new region.
823 label region = cutPointRegionMaster.size();
824 cutPointRegionMaster.append
826 min(cutPointI, otherCutPointI)
828 cutPointRegion[cutPointI] = region;
829 cutPointRegion[otherCutPointI] = region;
835 patchFp = patchF.fcIndex(patchFp);
839 patchFp = patchF.rcIndex(patchFp);
844 // Rework region&master into compaction array
845 compactToCut.setSize(cutPointRegion.size());
846 cutToCompact.setSize(cutPointRegion.size());
848 label compactPointI = 0;
850 forAll(cutPointRegion, i)
852 if (cutPointRegion[i] == -1)
854 // Unduplicated point. Allocate new compacted point.
855 cutToCompact[i] = compactPointI;
856 compactToCut[compactPointI] = i;
861 // Duplicate point. Get master.
863 label masterPointI = cutPointRegionMaster[cutPointRegion[i]];
865 if (cutToCompact[masterPointI] == -1)
867 cutToCompact[masterPointI] = compactPointI;
868 compactToCut[compactPointI] = masterPointI;
871 cutToCompact[i] = cutToCompact[masterPointI];
874 compactToCut.setSize(compactPointI);
876 return compactToCut.size() != cutToCompact.size();
880 // Return max distance from any point on cutF to masterF
881 Foam::scalar Foam::faceCoupleInfo::maxDistance
884 const pointField& cutPoints,
886 const pointField& masterPoints
889 scalar maxDist = -GREAT;
893 const point& cutPt = cutPoints[cutF[fp]];
895 pointHit pHit = masterF.nearestPoint(cutPt, masterPoints);
897 maxDist = max(maxDist, pHit.distance());
903 void Foam::faceCoupleInfo::findPerfectMatchingFaces
905 const primitiveMesh& mesh0,
906 const primitiveMesh& mesh1,
909 labelList& mesh0Faces,
910 labelList& mesh1Faces
913 // Face centres of external faces (without invoking
914 // mesh.faceCentres since mesh might have been clearedOut)
918 calcFaceCentres<List>
922 mesh0.nInternalFaces(),
923 mesh0.nFaces() - mesh0.nInternalFaces()
929 calcFaceCentres<List>
933 mesh1.nInternalFaces(),
934 mesh1.nFaces() - mesh1.nInternalFaces()
941 Pout<< "Face matching tolerance : " << absTol << endl;
945 // Match geometrically
947 bool matchedAllFaces = matchPoints
951 scalarField(fc1.size(), absTol),
959 << "faceCoupleInfo::faceCoupleInfo : "
960 << "Matched ALL " << fc1.size()
961 << " boundary faces of mesh0 to boundary faces of mesh1." << endl
962 << "This is only valid if the mesh to add is fully"
963 << " enclosed by the mesh it is added to." << endl;
970 mesh0Faces.setSize(fc0.size());
971 mesh1Faces.setSize(fc1.size());
975 if (from1To0[i] != -1)
977 mesh1Faces[nMatched] = i + mesh1.nInternalFaces();
978 mesh0Faces[nMatched] = from1To0[i] + mesh0.nInternalFaces();
984 mesh0Faces.setSize(nMatched);
985 mesh1Faces.setSize(nMatched);
989 void Foam::faceCoupleInfo::findSlavesCoveringMaster
991 const primitiveMesh& mesh0,
992 const primitiveMesh& mesh1,
995 labelList& mesh0Faces,
996 labelList& mesh1Faces
999 // Construct octree from all mesh0 boundary faces
1000 labelList bndFaces(mesh0.nFaces()-mesh0.nInternalFaces());
1003 bndFaces[i] = mesh0.nInternalFaces() + i;
1006 treeBoundBox overallBb(mesh0.points());
1008 Random rndGen(123456);
1010 indexedOctree<treeDataFace> tree
1012 treeDataFace // all information needed to search faces
1014 false, // do not cache bb
1016 bndFaces // boundary faces only
1018 overallBb.extend(rndGen, 1E-4), // overall search domain
1026 Pout<< "findSlavesCoveringMaster :"
1027 << " constructed octree for mesh0 boundary faces" << endl;
1030 // Search nearest face for every mesh1 boundary face.
1032 labelHashSet mesh0Set(mesh0.nFaces() - mesh0.nInternalFaces());
1033 labelHashSet mesh1Set(mesh1.nFaces() - mesh1.nInternalFaces());
1037 label mesh1FaceI = mesh1.nInternalFaces();
1038 mesh1FaceI < mesh1.nFaces();
1042 const face& f1 = mesh1.faces()[mesh1FaceI];
1044 // Generate face centre (prevent cellCentres() reconstruction)
1045 point fc(f1.centre(mesh1.points()));
1047 pointIndexHit nearInfo = tree.findNearest(fc, Foam::sqr(absTol));
1051 label mesh0FaceI = tree.shapes().faceLabels()[nearInfo.index()];
1053 // Check if points of f1 actually lie on top of mesh0 face
1054 // This is the bit that might fail since if f0 is severely warped
1055 // and f1's points are not present on f0 (since subdivision)
1056 // then the distance of the points to f0 might be quite large
1057 // and the test will fail. This all should in fact be some kind
1058 // of walk from known corresponding points/edges/faces.
1064 mesh0.faces()[mesh0FaceI],
1070 mesh0Set.insert(mesh0FaceI);
1071 mesh1Set.insert(mesh1FaceI);
1078 Pout<< "findSlavesCoveringMaster :"
1079 << " matched " << mesh1Set.size() << " mesh1 faces to "
1080 << mesh0Set.size() << " mesh0 faces" << endl;
1083 mesh0Faces = mesh0Set.toc();
1084 mesh1Faces = mesh1Set.toc();
1088 // Grow cutToMasterFace across 'internal' edges.
1089 Foam::label Foam::faceCoupleInfo::growCutFaces
1091 const labelList& cutToMasterEdges,
1092 Map<labelList>& candidates
1097 Pout<< "growCutFaces :"
1098 << " growing cut faces to masterPatch" << endl;
1101 label nTotChanged = 0;
1105 const labelListList& cutFaceEdges = cutFaces().faceEdges();
1106 const labelListList& cutEdgeFaces = cutFaces().edgeFaces();
1110 forAll(cutToMasterFaces_, cutFaceI)
1112 const label masterFaceI = cutToMasterFaces_[cutFaceI];
1114 if (masterFaceI != -1)
1116 // We now have a cutFace for which we have already found a
1117 // master face. Grow this masterface across any internal edge
1118 // (internal: no corresponding master edge)
1120 const labelList& fEdges = cutFaceEdges[cutFaceI];
1124 const label cutEdgeI = fEdges[i];
1126 if (cutToMasterEdges[cutEdgeI] == -1)
1129 // - internal to the cutPatch (since all region edges
1131 // - on cutFaceI which corresponds to masterFace.
1132 // Mark all connected faces with this masterFace.
1134 const labelList& eFaces = cutEdgeFaces[cutEdgeI];
1138 const label faceI = eFaces[j];
1140 if (cutToMasterFaces_[faceI] == -1)
1142 cutToMasterFaces_[faceI] = masterFaceI;
1143 candidates.erase(faceI);
1146 else if (cutToMasterFaces_[faceI] != masterFaceI)
1148 const pointField& cutPoints =
1149 cutFaces().points();
1150 const pointField& masterPoints =
1151 masterPatch().points();
1153 const edge& e = cutFaces().edges()[cutEdgeI];
1155 label myMaster = cutToMasterFaces_[faceI];
1156 const face& myF = masterPatch()[myMaster];
1158 const face& nbrF = masterPatch()[masterFaceI];
1162 "faceCoupleInfo::growCutFaces"
1163 "(const labelList&, Map<labelList>&)"
1165 << cutFaces()[faceI].points(cutPoints)
1168 << " but also connects to nbr face "
1169 << cutFaces()[cutFaceI].points(cutPoints)
1170 << " with master " << masterFaceI
1173 << myF.points(masterPoints)
1174 << " nbrMasterFace:"
1175 << nbrF.points(masterPoints) << nl
1176 << "Across cut edge " << cutEdgeI
1178 << cutFaces().localPoints()[e[0]]
1179 << cutFaces().localPoints()[e[1]]
1180 << abort(FatalError);
1190 Pout<< "growCutFaces : Grown an additional " << nChanged
1191 << " cut to master face" << " correspondence" << endl;
1194 nTotChanged += nChanged;
1206 void Foam::faceCoupleInfo::checkMatch(const labelList& cutToMasterEdges) const
1208 const pointField& cutLocalPoints = cutFaces().localPoints();
1210 const pointField& masterLocalPoints = masterPatch().localPoints();
1211 const faceList& masterLocalFaces = masterPatch().localFaces();
1213 forAll(cutToMasterEdges, cutEdgeI)
1215 const edge& e = cutFaces().edges()[cutEdgeI];
1217 if (cutToMasterEdges[cutEdgeI] == -1)
1219 // Internal edge. Check that master face is same on both sides.
1220 const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1222 label masterFaceI = -1;
1224 forAll(cutEFaces, i)
1226 label cutFaceI = cutEFaces[i];
1228 if (cutToMasterFaces_[cutFaceI] != -1)
1230 if (masterFaceI == -1)
1232 masterFaceI = cutToMasterFaces_[cutFaceI];
1234 else if (masterFaceI != cutToMasterFaces_[cutFaceI])
1236 label myMaster = cutToMasterFaces_[cutFaceI];
1237 const face& myF = masterLocalFaces[myMaster];
1239 const face& nbrF = masterLocalFaces[masterFaceI];
1243 "faceCoupleInfo::checkMatch(const labelList&) const"
1245 << "Internal CutEdge " << e
1247 << cutLocalPoints[e[0]]
1248 << cutLocalPoints[e[1]]
1249 << " connects to master " << myMaster
1250 << " and to master " << masterFaceI << nl
1252 << myF.points(masterLocalPoints)
1253 << " nbrMasterFace:"
1254 << nbrF.points(masterLocalPoints)
1255 << abort(FatalError);
1264 // Extends matching information by elimination across cutFaces using more
1265 // than one region edge. Updates cutToMasterFaces_ and sets candidates
1266 // which is for every cutface on a region edge the possible master faces.
1267 Foam::label Foam::faceCoupleInfo::matchEdgeFaces
1269 const labelList& cutToMasterEdges,
1270 Map<labelList>& candidates
1273 // For every unassigned cutFaceI the possible list of master faces.
1275 candidates.resize(cutFaces().size());
1279 forAll(cutToMasterEdges, cutEdgeI)
1281 label masterEdgeI = cutToMasterEdges[cutEdgeI];
1283 if (masterEdgeI != -1)
1285 // So cutEdgeI is matched to masterEdgeI. For all cut faces using
1286 // the cut edge store the master face as a candidate match.
1288 const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1289 const labelList& masterEFaces =
1290 masterPatch().edgeFaces()[masterEdgeI];
1292 forAll(cutEFaces, i)
1294 label cutFaceI = cutEFaces[i];
1296 if (cutToMasterFaces_[cutFaceI] == -1)
1298 // So this face (cutFaceI) is on an edge (cutEdgeI) that
1299 // is used by some master faces (masterEFaces). Check if
1300 // the cutFaceI and some of these masterEFaces share more
1301 // than one edge (which uniquely defines face)
1303 // Combine master faces with current set of candidate
1305 Map<labelList>::iterator fnd = candidates.find(cutFaceI);
1307 if (fnd == candidates.end())
1309 // No info yet for cutFaceI. Add all master faces as
1311 candidates.insert(cutFaceI, masterEFaces);
1315 // From some other cutEdgeI there are already some
1316 // candidate master faces. Check the overlap with
1317 // the current set of master faces.
1318 const labelList& masterFaces = fnd();
1320 DynamicList<label> newCandidates(masterFaces.size());
1322 forAll(masterEFaces, j)
1324 if (findIndex(masterFaces, masterEFaces[j]) != -1)
1326 newCandidates.append(masterEFaces[j]);
1330 if (newCandidates.size() == 1)
1332 // We found a perfect match. Delete entry from
1334 cutToMasterFaces_[cutFaceI] = newCandidates[0];
1335 candidates.erase(cutFaceI);
1340 // Should not happen?
1341 //Pout<< "On edge:" << cutEdgeI
1342 // << " have connected masterFaces:"
1344 // << " and from previous edge we have"
1345 // << " connected masterFaces" << masterFaces
1346 // << " . Overlap:" << newCandidates << endl;
1348 fnd() = newCandidates.shrink();
1359 Pout<< "matchEdgeFaces : Found " << nChanged
1360 << " faces where there was"
1361 << " only one remaining choice for cut-master correspondence"
1369 // Gets a list of cutFaces (that use a master edge) and the candidate
1371 // Finds most aligned master face.
1372 Foam::label Foam::faceCoupleInfo::geometricMatchEdgeFaces
1374 Map<labelList>& candidates
1377 const pointField& cutPoints = cutFaces().points();
1381 // Mark all master faces that have been matched so far.
1383 labelListList masterToCutFaces
1387 masterPatch().size(),
1392 forAllConstIter(Map<labelList>, candidates, iter)
1394 label cutFaceI = iter.key();
1396 const face& cutF = cutFaces()[cutFaceI];
1398 if (cutToMasterFaces_[cutFaceI] == -1)
1400 const labelList& masterFaces = iter();
1402 // Find the best matching master face.
1403 scalar minDist = GREAT;
1404 label minMasterFaceI = -1;
1406 forAll(masterFaces, i)
1408 label masterFaceI = masterFaces[i];
1410 if (masterToCutFaces[masterFaces[i]].empty())
1412 scalar dist = maxDistance
1416 masterPatch()[masterFaceI],
1417 masterPatch().points()
1423 minMasterFaceI = masterFaceI;
1428 if (minMasterFaceI != -1)
1430 cutToMasterFaces_[cutFaceI] = minMasterFaceI;
1431 masterToCutFaces[minMasterFaceI] = cutFaceI;
1437 // (inefficiently) force candidates to be uptodate.
1438 forAll(cutToMasterFaces_, cutFaceI)
1440 if (cutToMasterFaces_[cutFaceI] != -1)
1442 candidates.erase(cutFaceI);
1449 Pout<< "geometricMatchEdgeFaces : Found " << nChanged
1450 << " faces where there was"
1451 << " only one remaining choice for cut-master correspondence"
1459 // Calculate the set of cut faces inbetween master and slave patch
1460 // assuming perfect match (and optional face ordering on slave)
1461 void Foam::faceCoupleInfo::perfectPointMatch
1463 const scalar absTol,
1464 const bool slaveFacesOrdered
1469 Pout<< "perfectPointMatch :"
1470 << " Matching master and slave to cut."
1471 << " Master and slave faces are identical" << nl;
1473 if (slaveFacesOrdered)
1475 Pout<< "and master and slave faces are ordered"
1476 << " (on coupled patches)" << endl;
1480 Pout<< "and master and slave faces are not ordered"
1481 << " (on coupled patches)" << endl;
1485 cutToMasterFaces_ = identity(masterPatch().size());
1486 cutPoints_ = masterPatch().localPoints();
1489 new primitiveFacePatch
1491 masterPatch().localFaces(),
1495 masterToCutPoints_ = identity(cutPoints_.size());
1498 // Cut faces to slave patch.
1499 bool matchedAllFaces = false;
1501 if (slaveFacesOrdered)
1503 cutToSlaveFaces_ = identity(cutFaces().size());
1504 matchedAllFaces = (cutFaces().size() == slavePatch().size());
1508 // Faces do not have to be ordered (but all have
1509 // to match). Note: Faces will be already ordered if we enter here from
1510 // construct from meshes.
1511 matchedAllFaces = matchPoints
1513 calcFaceCentres<List>
1520 calcFaceCentres<IndirectList>
1523 slavePatch().points(),
1527 scalarField(slavePatch().size(), absTol),
1534 if (!matchedAllFaces)
1538 "faceCoupleInfo::perfectPointMatch"
1539 "(const scalar&, const bool)"
1540 ) << "Did not match all of the master faces to the slave faces"
1542 << "This usually means that the slave patch and master patch"
1543 << " do not align to within " << absTol << " meter."
1544 << abort(FatalError);
1548 // Find correspondence from slave points to cut points. This might
1549 // detect shared points so the output is a slave-to-cut point list
1550 // and a compaction list.
1552 labelList cutToCompact, compactToCut;
1553 matchPointsThroughFaces
1556 cutFaces().localPoints(),
1557 reorder(cutToSlaveFaces_, cutFaces().localFaces()),
1558 slavePatch().localPoints(),
1559 slavePatch().localFaces(),
1560 false, // slave and cut have opposite orientation
1562 slaveToCutPoints_, // slave to (uncompacted) cut points
1563 cutToCompact, // compaction map: from cut to compacted
1564 compactToCut // compaction map: from compacted to cut
1568 // Use compaction lists to renumber cutPoints.
1569 cutPoints_ = UIndirectList<point>(cutPoints_, compactToCut)();
1571 const faceList& cutLocalFaces = cutFaces().localFaces();
1573 faceList compactFaces(cutLocalFaces.size());
1574 forAll(cutLocalFaces, i)
1576 compactFaces[i] = renumber(cutToCompact, cutLocalFaces[i]);
1580 new primitiveFacePatch
1587 inplaceRenumber(cutToCompact, slaveToCutPoints_);
1588 inplaceRenumber(cutToCompact, masterToCutPoints_);
1592 // Calculate the set of cut faces inbetween master and slave patch
1593 // assuming that slave patch is subdivision of masterPatch.
1594 void Foam::faceCoupleInfo::subDivisionMatch
1596 const polyMesh& slaveMesh,
1597 const bool patchDivision,
1603 Pout<< "subDivisionMatch :"
1604 << " Matching master and slave to cut."
1605 << " Slave can be subdivision of master but all master edges"
1606 << " have to be covered by slave edges." << endl;
1609 // Assume slave patch is subdivision of the master patch so cutFaces is a
1610 // copy of the slave (but reversed since orientation has to be like master).
1611 cutPoints_ = slavePatch().localPoints();
1613 faceList cutFaces(slavePatch().size());
1617 cutFaces[i] = slavePatch().localFaces()[i].reverseFace();
1619 cutFacesPtr_.reset(new primitiveFacePatch(cutFaces, cutPoints_));
1622 // Cut is copy of slave so addressing to slave is trivial.
1623 cutToSlaveFaces_ = identity(cutFaces().size());
1624 slaveToCutPoints_ = identity(slavePatch().nPoints());
1626 // Cut edges to slave patch
1627 labelList cutToSlaveEdges
1632 slaveToCutPoints_, //note:should be cutToSlavePoints but since iden
1640 OFstream str("regionEdges.obj");
1642 Pout<< "subDivisionMatch :"
1643 << " addressing for slave patch fully done."
1644 << " Dumping region edges to " << str.name() << endl;
1648 forAll(slavePatch().edges(), slaveEdgeI)
1650 if (regionEdge(slaveMesh, slaveEdgeI))
1652 const edge& e = slavePatch().edges()[slaveEdgeI];
1654 meshTools::writeOBJ(str, slavePatch().localPoints()[e[0]]);
1656 meshTools::writeOBJ(str, slavePatch().localPoints()[e[1]]);
1658 str<< "l " << vertI-1 << ' ' << vertI << nl;
1664 // Addressing from cut to master.
1666 // Cut points to master patch
1667 // - determine master points to cut points. Has to be full!
1668 // - invert to get cut to master
1671 Pout<< "subDivisionMatch :"
1672 << " matching master points to cut points" << endl;
1675 bool matchedAllPoints = matchPoints
1677 masterPatch().localPoints(),
1679 scalarField(masterPatch().nPoints(), absTol),
1684 if (!matchedAllPoints)
1688 "faceCoupleInfo::subDivisionMatch"
1689 "(const polyMesh&, const bool, const scalar)"
1690 ) << "Did not match all of the master points to the slave points"
1692 << "This usually means that the slave patch is not a"
1693 << " subdivision of the master patch"
1694 << abort(FatalError);
1698 // Do masterEdges to cutEdges :
1699 // - mark all edges between two masterEdge endpoints. (geometric test since
1700 // this is the only distinction)
1701 // - this gives the borders inbetween which all cutfaces come from
1702 // a single master face.
1705 Pout<< "subDivisionMatch :"
1706 << " matching cut edges to master edges" << endl;
1709 const edgeList& masterEdges = masterPatch().edges();
1710 const pointField& masterPoints = masterPatch().localPoints();
1712 const edgeList& cutEdges = cutFaces().edges();
1714 labelList cutToMasterEdges(cutFaces().nEdges(), -1);
1716 forAll(masterEdges, masterEdgeI)
1718 const edge& masterEdge = masterEdges[masterEdgeI];
1720 label cutPoint0 = masterToCutPoints_[masterEdge[0]];
1721 label cutPoint1 = masterToCutPoints_[masterEdge[1]];
1723 // Find edges between cutPoint0 and cutPoint1.
1725 label cutPointI = cutPoint0;
1728 // Find edge (starting at pointI on cut), aligned with master
1745 // Problem. Report while matching to produce nice error message
1758 Pout<< "Dumping unmatched pointEdges to errorEdges.obj"
1767 cutFaces().pointEdges()[cutPointI]
1769 cutFaces().localPoints(),
1775 "faceCoupleInfo::subDivisionMatch"
1776 "(const polyMesh&, const bool, const scalar)"
1777 ) << "Problem in finding cut edges corresponding to"
1778 << " master edge " << masterEdge
1779 << " points:" << masterPoints[masterEdge[0]]
1780 << ' ' << masterPoints[masterEdge[1]]
1781 << " corresponding cut edge: (" << cutPoint0
1783 << abort(FatalError);
1786 cutToMasterEdges[cutEdgeI] = masterEdgeI;
1788 cutPointI = cutEdges[cutEdgeI].otherVertex(cutPointI);
1790 } while(cutPointI != cutPoint1);
1795 // Write all matched edges.
1796 writeEdges(cutToMasterEdges, cutToSlaveEdges);
1799 // Rework cutToMasterEdges into list of points inbetween two endpoints
1800 // (cutEdgeToPoints_)
1801 setCutEdgeToPoints(cutToMasterEdges);
1804 // Now that we have marked the cut edges that lie on top of master edges
1805 // we can find cut faces that have two (or more) of these edges.
1806 // Note that there can be multiple faces having two or more matched edges
1807 // since the cut faces can form a non-manifold surface(?)
1808 // So the matching is done as an elimination process where for every
1809 // cutFace on a matched edge we store the possible master faces and
1810 // eliminate more and more until we only have one possible master face
1815 Pout<< "subDivisionMatch :"
1816 << " matching (topological) cut faces to masterPatch" << endl;
1819 // For every unassigned cutFaceI the possible list of master faces.
1820 Map<labelList> candidates(cutFaces().size());
1822 cutToMasterFaces_.setSize(cutFaces().size());
1823 cutToMasterFaces_ = -1;
1827 // See if there are any edges left where there is only one remaining
1829 label nChanged = matchEdgeFaces(cutToMasterEdges, candidates);
1831 checkMatch(cutToMasterEdges);
1834 // Now we should have found a face correspondence for every cutFace
1835 // that uses two or more cutEdges. Floodfill from these across
1836 // 'internal' cutedges (i.e. edges that do not have a master
1837 // equivalent) to determine some more cutToMasterFaces
1838 nChanged += growCutFaces(cutToMasterEdges, candidates);
1840 checkMatch(cutToMasterEdges);
1850 Pout<< "subDivisionMatch :"
1851 << " matching (geometric) cut faces to masterPatch" << endl;
1854 // Above should have matched all cases fully. Cannot prove this yet so
1855 // do any remaining unmatched edges by a geometric test.
1858 // See if there are any edges left where there is only one remaining
1860 label nChanged = geometricMatchEdgeFaces(candidates);
1862 checkMatch(cutToMasterEdges);
1864 nChanged += growCutFaces(cutToMasterEdges, candidates);
1866 checkMatch(cutToMasterEdges);
1875 // All cut faces matched?
1876 forAll(cutToMasterFaces_, cutFaceI)
1878 if (cutToMasterFaces_[cutFaceI] == -1)
1880 const face& cutF = cutFaces()[cutFaceI];
1884 "faceCoupleInfo::subDivisionMatch"
1885 "(const polyMesh&, const bool, const scalar)"
1886 ) << "Did not match all of cutFaces to a master face" << nl
1887 << "First unmatched cut face:" << cutFaceI << " with points:"
1888 << UIndirectList<point>(cutFaces().points(), cutF)()
1890 << "This usually means that the slave patch is not a"
1891 << " subdivision of the master patch"
1892 << abort(FatalError);
1898 Pout<< "subDivisionMatch :"
1899 << " finished matching master and slave to cut" << endl;
1909 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1911 // Construct from mesh data
1912 Foam::faceCoupleInfo::faceCoupleInfo
1914 const polyMesh& masterMesh,
1915 const polyMesh& slaveMesh,
1916 const scalar absTol,
1917 const bool perfectMatch
1920 masterPatchPtr_(NULL),
1921 slavePatchPtr_(NULL),
1924 cutToMasterFaces_(0),
1925 masterToCutPoints_(0),
1926 cutToSlaveFaces_(0),
1927 slaveToCutPoints_(0),
1930 // Get faces on both meshes that are aligned.
1931 // (not ordered i.e. masterToMesh[0] does
1932 // not couple to slaveToMesh[0]. This ordering is done later on)
1933 labelList masterToMesh;
1934 labelList slaveToMesh;
1938 // Identical faces so use tight face-centre comparison.
1939 findPerfectMatchingFaces
1950 // Slave subdivision of master so use 'nearest'. Bit dodgy, especially
1951 // with using absTol (which is quite small)
1952 findSlavesCoveringMaster
1962 // Construct addressing engines for both sides
1963 masterPatchPtr_.reset
1965 new indirectPrimitivePatch
1967 IndirectList<face>(masterMesh.faces(), masterToMesh),
1972 slavePatchPtr_.reset
1974 new indirectPrimitivePatch
1976 IndirectList<face>(slaveMesh.faces(), slaveToMesh),
1984 // Faces are perfectly aligned but probably not ordered.
1985 perfectPointMatch(absTol, false);
1989 // Slave faces are subdivision of master face. Faces not ordered.
1990 subDivisionMatch(slaveMesh, false, absTol);
2000 // Slave is subdivision of master patch.
2001 // (so -both cover the same area -all of master points are present in slave)
2002 Foam::faceCoupleInfo::faceCoupleInfo
2004 const polyMesh& masterMesh,
2005 const labelList& masterAddressing,
2006 const polyMesh& slaveMesh,
2007 const labelList& slaveAddressing,
2008 const scalar absTol,
2009 const bool perfectMatch,
2010 const bool orderedFaces,
2011 const bool patchDivision
2016 new indirectPrimitivePatch
2018 IndirectList<face>(masterMesh.faces(), masterAddressing),
2024 new indirectPrimitivePatch
2026 IndirectList<face>(slaveMesh.faces(), slaveAddressing),
2032 cutToMasterFaces_(0),
2033 masterToCutPoints_(0),
2034 cutToSlaveFaces_(0),
2035 slaveToCutPoints_(0),
2038 if (perfectMatch && (masterAddressing.size() != slaveAddressing.size()))
2042 "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2043 ", const primitiveMesh&, const scalar, const bool"
2044 ) << "Perfect match specified but number of master and slave faces"
2045 << " differ." << endl
2046 << "master:" << masterAddressing.size()
2047 << " slave:" << slaveAddressing.size()
2048 << abort(FatalError);
2053 masterAddressing.size()
2054 && min(masterAddressing) < masterMesh.nInternalFaces()
2059 "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2060 ", const primitiveMesh&, const scalar, const bool"
2061 ) << "Supplied internal face on master mesh to couple." << nl
2062 << "Faces to be coupled have to be boundary faces."
2063 << abort(FatalError);
2067 slaveAddressing.size()
2068 && min(slaveAddressing) < slaveMesh.nInternalFaces()
2073 "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2074 ", const primitiveMesh&, const scalar, const bool"
2075 ) << "Supplied internal face on slave mesh to couple." << nl
2076 << "Faces to be coupled have to be boundary faces."
2077 << abort(FatalError);
2083 perfectPointMatch(absTol, orderedFaces);
2087 // Slave faces are subdivision of master face. Faces not ordered.
2088 subDivisionMatch(slaveMesh, patchDivision, absTol);
2098 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
2100 Foam::faceCoupleInfo::~faceCoupleInfo()
2104 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2106 Foam::labelList Foam::faceCoupleInfo::faceLabels(const polyPatch& pp)
2108 labelList faces(pp.size());
2110 label faceI = pp.start();
2120 Foam::Map<Foam::label> Foam::faceCoupleInfo::makeMap(const labelList& lst)
2122 Map<label> map(lst.size());
2128 map.insert(i, lst[i]);
2135 Foam::Map<Foam::labelList> Foam::faceCoupleInfo::makeMap
2137 const labelListList& lst
2140 Map<labelList> map(lst.size());
2146 map.insert(i, lst[i]);
2153 // ************************ vim: set sw=4 sts=4 et: ************************ //