faceCoupleInfo.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "faceCoupleInfo.H"
27 #include "polyMesh.H"
28 #include "matchPoints.H"
29 #include "indirectPrimitivePatch.H"
30 #include "meshTools.H"
31 #include "treeDataFace.H"
32 #include "indexedOctree.H"
33 #include "OFstream.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 defineTypeNameAndDebug(faceCoupleInfo, 0);
40 
41 const scalar faceCoupleInfo::angleTol_ = 1e-3;
42 }
43 
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 //- Write edges
48 void Foam::faceCoupleInfo::writeOBJ
49 (
50  const fileName& fName,
51  const edgeList& edges,
52  const pointField& points,
53  const bool compact
54 )
55 {
56  OFstream str(fName);
57 
58  labelList pointMap(points.size(), -1);
59 
60  if (compact)
61  {
62  label newPointI = 0;
63 
64  forAll(edges, edgeI)
65  {
66  const edge& e = edges[edgeI];
67 
68  forAll(e, eI)
69  {
70  label pointI = e[eI];
71 
72  if (pointMap[pointI] == -1)
73  {
74  pointMap[pointI] = newPointI++;
75 
76  meshTools::writeOBJ(str, points[pointI]);
77  }
78  }
79  }
80  }
81  else
82  {
83  forAll(points, pointI)
84  {
85  meshTools::writeOBJ(str, points[pointI]);
86  }
87 
88  pointMap = identity(points.size());
89  }
90 
91  forAll(edges, edgeI)
92  {
93  const edge& e = edges[edgeI];
94 
95  str<< "l " << pointMap[e[0]]+1 << ' ' << pointMap[e[1]]+1 << nl;
96  }
97 }
98 
99 
100 //- Writes edges.
101 void Foam::faceCoupleInfo::writeOBJ
102 (
103  const fileName& fName,
104  const pointField& points0,
105  const pointField& points1
106 )
107 {
108  Pout<< "Writing connections as edges to " << fName << endl;
109 
110  OFstream str(fName);
111 
112  label vertI = 0;
113 
114  forAll(points0, i)
115  {
116  meshTools::writeOBJ(str, points0[i]);
117  vertI++;
118  meshTools::writeOBJ(str, points1[i]);
119  vertI++;
120  str << "l " << vertI-1 << ' ' << vertI << nl;
121  }
122 }
123 
124 
125 //- Writes face and point connectivity as .obj files.
126 void Foam::faceCoupleInfo::writePointsFaces() const
127 {
128  const indirectPrimitivePatch& m = masterPatch();
129  const indirectPrimitivePatch& s = slavePatch();
130  const primitiveFacePatch& c = cutFaces();
131 
132  // Patches
133  {
134  OFstream str("masterPatch.obj");
135  Pout<< "Writing masterPatch to " << str.name() << endl;
136  meshTools::writeOBJ(str, m.localFaces(), m.localPoints());
137  }
138  {
139  OFstream str("slavePatch.obj");
140  Pout<< "Writing slavePatch to " << str.name() << endl;
141  meshTools::writeOBJ(str, s.localFaces(), s.localPoints());
142  }
143  {
144  OFstream str("cutFaces.obj");
145  Pout<< "Writing cutFaces to " << str.name() << endl;
146  meshTools::writeOBJ(str, c.localFaces(), c.localPoints());
147  }
148 
149  // Point connectivity
150  {
151  Pout<< "Writing cutToMasterPoints to cutToMasterPoints.obj" << endl;
152 
153  writeOBJ
154  (
155  "cutToMasterPoints.obj",
156  m.localPoints(),
157  pointField(c.localPoints(), masterToCutPoints_));
158  }
159  {
160  Pout<< "Writing cutToSlavePoints to cutToSlavePoints.obj" << endl;
161 
162  writeOBJ
163  (
164  "cutToSlavePoints.obj",
165  s.localPoints(),
166  pointField(c.localPoints(), slaveToCutPoints_)
167  );
168  }
169 
170  // Face connectivity
171  {
172  Pout<< "Writing cutToMasterFaces to cutToMasterFaces.obj" << endl;
173 
174  pointField equivMasterFaces(c.size());
175 
176  forAll(cutToMasterFaces(), cutFaceI)
177  {
178  label masterFaceI = cutToMasterFaces()[cutFaceI];
179 
180  if (masterFaceI != -1)
181  {
182  equivMasterFaces[cutFaceI] = m[masterFaceI].centre(m.points());
183  }
184  else
185  {
186  WarningIn("writePointsFaces()")
187  << "No master face for cut face " << cutFaceI
188  << " at position " << c[cutFaceI].centre(c.points())
189  << endl;
190 
191  equivMasterFaces[cutFaceI] = vector::zero;
192  }
193  }
194 
195  writeOBJ
196  (
197  "cutToMasterFaces.obj",
198  calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
199  equivMasterFaces
200  );
201  }
202 
203  {
204  Pout<< "Writing cutToSlaveFaces to cutToSlaveFaces.obj" << endl;
205 
206  pointField equivSlaveFaces(c.size());
207 
208  forAll(cutToSlaveFaces(), cutFaceI)
209  {
210  label slaveFaceI = cutToSlaveFaces()[cutFaceI];
211 
212  equivSlaveFaces[cutFaceI] = s[slaveFaceI].centre(s.points());
213  }
214 
215  writeOBJ
216  (
217  "cutToSlaveFaces.obj",
218  calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
219  equivSlaveFaces
220  );
221  }
222 
223  Pout<< endl;
224 }
225 
226 
227 void Foam::faceCoupleInfo::writeEdges
228 (
229  const labelList& cutToMasterEdges,
230  const labelList& cutToSlaveEdges
231 ) const
232 {
233  const indirectPrimitivePatch& m = masterPatch();
234  const indirectPrimitivePatch& s = slavePatch();
235  const primitiveFacePatch& c = cutFaces();
236 
237  // Edge connectivity
238  {
239  OFstream str("cutToMasterEdges.obj");
240  Pout<< "Writing cutToMasterEdges to " << str.name() << endl;
241 
242  label vertI = 0;
243 
244  forAll(cutToMasterEdges, cutEdgeI)
245  {
246  if (cutToMasterEdges[cutEdgeI] != -1)
247  {
248  const edge& masterEdge =
249  m.edges()[cutToMasterEdges[cutEdgeI]];
250  const edge& cutEdge = c.edges()[cutEdgeI];
251 
252  meshTools::writeOBJ(str, m.localPoints()[masterEdge[0]]);
253  vertI++;
254  meshTools::writeOBJ(str, m.localPoints()[masterEdge[1]]);
255  vertI++;
256  meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
257  vertI++;
258  meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
259  vertI++;
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;
266  }
267  }
268  }
269  {
270  OFstream str("cutToSlaveEdges.obj");
271  Pout<< "Writing cutToSlaveEdges to " << str.name() << endl;
272 
273  label vertI = 0;
274 
275  labelList slaveToCut(invert(s.nEdges(), cutToSlaveEdges));
276 
277  forAll(slaveToCut, edgeI)
278  {
279  if (slaveToCut[edgeI] != -1)
280  {
281  const edge& slaveEdge = s.edges()[edgeI];
282  const edge& cutEdge = c.edges()[slaveToCut[edgeI]];
283 
284  meshTools::writeOBJ(str, s.localPoints()[slaveEdge[0]]);
285  vertI++;
286  meshTools::writeOBJ(str, s.localPoints()[slaveEdge[1]]);
287  vertI++;
288  meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
289  vertI++;
290  meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
291  vertI++;
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;
298  }
299  }
300  }
301 
302  Pout<< endl;
303 }
304 
305 
306 // Given an edgelist and a map for the points on the edges it tries to find
307 // the corresponding patch edges.
308 Foam::labelList Foam::faceCoupleInfo::findMappedEdges
309 (
310  const edgeList& edges,
311  const labelList& pointMap,
312  const indirectPrimitivePatch& patch
313 )
314 {
315  labelList toPatchEdges(edges.size());
316 
317  forAll(toPatchEdges, edgeI)
318  {
319  const edge& e = edges[edgeI];
320 
321  label v0 = pointMap[e[0]];
322  label v1 = pointMap[e[1]];
323 
324  toPatchEdges[edgeI] =
326  (
327  patch.edges(),
328  patch.pointEdges()[v0],
329  v0,
330  v1
331  );
332  }
333  return toPatchEdges;
334 }
335 
336 
337 // Detect a cut edge which originates from two boundary faces having different
338 // polyPatches.
339 bool Foam::faceCoupleInfo::regionEdge
340 (
341  const polyMesh& slaveMesh,
342  const label slaveEdgeI
343 ) const
344 {
345  const labelList& eFaces = slavePatch().edgeFaces()[slaveEdgeI];
346 
347  if (eFaces.size() == 1)
348  {
349  return true;
350  }
351  else
352  {
353  // Count how many different patches connected to this edge.
354 
355  label patch0 = -1;
356 
357  forAll(eFaces, i)
358  {
359  label faceI = eFaces[i];
360 
361  label meshFaceI = slavePatch().addressing()[faceI];
362 
363  label patchI = slaveMesh.boundaryMesh().whichPatch(meshFaceI);
364 
365  if (patch0 == -1)
366  {
367  patch0 = patchI;
368  }
369  else if (patchI != patch0)
370  {
371  // Found two different patches connected to this edge.
372  return true;
373  }
374  }
375  return false;
376  }
377 }
378 
379 
380 // Find edge using pointI that is most aligned with vector between
381 // master points. Patchdivision tells us whether or not to use
382 // patch information to match edges.
383 Foam::label Foam::faceCoupleInfo::mostAlignedCutEdge
384 (
385  const bool report,
386  const polyMesh& slaveMesh,
387  const bool patchDivision,
388  const labelList& cutToMasterEdges,
389  const labelList& cutToSlaveEdges,
390  const label pointI,
391  const label edgeStart,
392  const label edgeEnd
393 ) const
394 {
395  const pointField& localPoints = cutFaces().localPoints();
396 
397  const labelList& pEdges = cutFaces().pointEdges()[pointI];
398 
399  if (report)
400  {
401  Pout<< "mostAlignedEdge : finding nearest edge among "
402  << UIndirectList<edge>(cutFaces().edges(), pEdges)()
403  << " connected to point " << pointI
404  << " coord:" << localPoints[pointI]
405  << " running between " << edgeStart << " coord:"
406  << localPoints[edgeStart]
407  << " and " << edgeEnd << " coord:"
408  << localPoints[edgeEnd]
409  << endl;
410  }
411 
412  // Find the edge that gets us nearest end.
413 
414  label maxEdgeI = -1;
415  scalar maxCos = -GREAT;
416 
417  forAll(pEdges, i)
418  {
419  label edgeI = pEdges[i];
420 
421  if
422  (
423  !(
424  patchDivision
425  && cutToMasterEdges[edgeI] == -1
426  )
427  || (
428  patchDivision
429  && regionEdge(slaveMesh, cutToSlaveEdges[edgeI])
430  )
431  )
432  {
433  const edge& e = cutFaces().edges()[edgeI];
434 
435  label otherPointI = e.otherVertex(pointI);
436 
437  if (otherPointI == edgeEnd)
438  {
439  // Shortcut: found edge end point.
440  if (report)
441  {
442  Pout<< " mostAlignedEdge : found end point " << edgeEnd
443  << endl;
444  }
445  return edgeI;
446  }
447 
448  // Get angle between edge and edge to masterEnd
449 
450  vector eVec(localPoints[otherPointI] - localPoints[pointI]);
451 
452  scalar magEVec = mag(eVec);
453 
454  if (magEVec < VSMALL)
455  {
456  WarningIn("faceCoupleInfo::mostAlignedEdge")
457  << "Crossing zero sized edge " << edgeI
458  << " coords:" << localPoints[otherPointI]
459  << localPoints[pointI]
460  << " when walking from " << localPoints[edgeStart]
461  << " to " << localPoints[edgeEnd]
462  << endl;
463  return edgeI;
464  }
465 
466  eVec /= magEVec;
467 
468  vector eToEndPoint(localPoints[edgeEnd] - localPoints[otherPointI]);
469  eToEndPoint /= mag(eToEndPoint);
470 
471  scalar cosAngle = eVec & eToEndPoint;
472 
473  if (report)
474  {
475  Pout<< " edge:" << e << " points:" << localPoints[pointI]
476  << localPoints[otherPointI]
477  << " vec:" << eVec
478  << " vecToEnd:" << eToEndPoint
479  << " cosAngle:" << cosAngle
480  << endl;
481  }
482 
483  if (cosAngle > maxCos)
484  {
485  maxCos = cosAngle;
486  maxEdgeI = edgeI;
487  }
488  }
489  }
490 
491  if (maxCos > 1 - angleTol_)
492  {
493  return maxEdgeI;
494  }
495  else
496  {
497  return -1;
498  }
499 }
500 
501 
502 // Construct points to split points map (in cut addressing)
503 void Foam::faceCoupleInfo::setCutEdgeToPoints(const labelList& cutToMasterEdges)
504 {
505  labelListList masterToCutEdges
506  (
508  (
509  masterPatch().nEdges(),
510  cutToMasterEdges
511  )
512  );
513 
514  const edgeList& cutEdges = cutFaces().edges();
515 
516  // Size extra big so searching is faster
517  cutEdgeToPoints_.resize
518  (
519  masterPatch().nEdges()
520  + slavePatch().nEdges()
521  + cutEdges.size()
522  );
523 
524  forAll(masterToCutEdges, masterEdgeI)
525  {
526  const edge& masterE = masterPatch().edges()[masterEdgeI];
527 
528  //Pout<< "Master:" << masterPatch().localPoints()[masterE[0]] << ' '
529  // << masterPatch().localPoints()[masterE[1]] << endl;
530 
531  const labelList& stringedEdges = masterToCutEdges[masterEdgeI];
532 
533  if (stringedEdges.empty())
534  {
536  (
537  "faceCoupleInfo::setCutEdgeToPoints"
538  "(const labelList&)"
539  ) << "Did not match all of master edges to cutFace edges"
540  << nl
541  << "First unmatched edge:" << masterEdgeI << " endPoints:"
542  << masterPatch().localPoints()[masterE[0]]
543  << masterPatch().localPoints()[masterE[1]]
544  << endl
545  << "This usually means that the slave patch is not a"
546  << " subdivision of the master patch"
547  << abort(FatalError);
548  }
549  else if (stringedEdges.size() > 1)
550  {
551  // String up the edges between e[0] and e[1]. Store the points
552  // inbetween e[0] and e[1] (all in cutFaces() labels)
553 
554  DynamicList<label> splitPoints(stringedEdges.size()-1);
555 
556  // Unsplit edge endpoints
557  const edge unsplitEdge
558  (
559  masterToCutPoints_[masterE[0]],
560  masterToCutPoints_[masterE[1]]
561  );
562 
563  label startVertI = unsplitEdge[0];
564  label startEdgeI = -1;
565 
566  while (startVertI != unsplitEdge[1])
567  {
568  // Loop over all string of edges. Update
569  // - startVertI : previous vertex
570  // - startEdgeI : previous edge
571  // and insert any points into splitPoints
572 
573  // For checking
574  label oldStart = startVertI;
575 
576  forAll(stringedEdges, i)
577  {
578  label edgeI = stringedEdges[i];
579 
580  if (edgeI != startEdgeI)
581  {
582  const edge& e = cutEdges[edgeI];
583 
584  //Pout<< " cut:" << e << " at:"
585  // << cutFaces().localPoints()[e[0]]
586  // << ' ' << cutFaces().localPoints()[e[1]] << endl;
587 
588  if (e[0] == startVertI)
589  {
590  startEdgeI = edgeI;
591  startVertI = e[1];
592  if (e[1] != unsplitEdge[1])
593  {
594  splitPoints.append(e[1]);
595  }
596  break;
597  }
598  else if (e[1] == startVertI)
599  {
600  startEdgeI = edgeI;
601  startVertI = e[0];
602  if (e[0] != unsplitEdge[1])
603  {
604  splitPoints.append(e[0]);
605  }
606  break;
607  }
608  }
609  }
610 
611  // Check
612  if (oldStart == startVertI)
613  {
615  (
616  "faceCoupleInfo::setCutEdgeToPoints"
617  "(const labelList&)"
618  ) << " unsplitEdge:" << unsplitEdge
619  << " does not correspond to split edges "
620  << UIndirectList<edge>(cutEdges, stringedEdges)()
621  << abort(FatalError);
622  }
623  }
624 
625  //Pout<< "For master edge:"
626  // << unsplitEdge
627  // << " Found stringed points "
628  // << UIndirectList<point>
629  // (
630  // cutFaces().localPoints(),
631  // splitPoints.shrink()
632  // )()
633  // << endl;
634 
635  cutEdgeToPoints_.insert(unsplitEdge, splitPoints.shrink());
636  }
637  }
638 }
639 
640 
641 // Determines rotation for f1 to match up with f0, i.e. the index in f0 of
642 // the first point of f1.
643 Foam::label Foam::faceCoupleInfo::matchFaces
644 (
645  const scalar absTol,
646  const pointField& points0,
647  const face& f0,
648  const pointField& points1,
649  const face& f1,
650  const bool sameOrientation
651 )
652 {
653  if (f0.size() != f1.size())
654  {
656  (
657  "faceCoupleInfo::matchFaces"
658  "(const scalar, const face&, const pointField&"
659  ", const face&, const pointField&)"
660  ) << "Different sizes for supposedly matching faces." << nl
661  << "f0:" << f0 << " coords:" << UIndirectList<point>(points0, f0)()
662  << nl
663  << "f1:" << f1 << " coords:" << UIndirectList<point>(points1, f1)()
664  << abort(FatalError);
665  }
666 
667  const scalar absTolSqr = sqr(absTol);
668 
669 
670  label matchFp = -1;
671 
672  forAll(f0, startFp)
673  {
674  // See -if starting from startFp on f0- the two faces match.
675  bool fullMatch = true;
676 
677  label fp0 = startFp;
678  label fp1 = 0;
679 
680  forAll(f1, i)
681  {
682  scalar distSqr = Foam::magSqr(points0[f0[fp0]] - points1[f1[fp1]]);
683 
684  if (distSqr > absTolSqr)
685  {
686  fullMatch = false;
687  break;
688  }
689 
690  fp0 = f0.fcIndex(fp0); // walk forward
691 
692  if (sameOrientation)
693  {
694  fp1 = f1.fcIndex(fp1);
695  }
696  else
697  {
698  fp1 = f1.rcIndex(fp1);
699  }
700  }
701 
702  if (fullMatch)
703  {
704  matchFp = startFp;
705  break;
706  }
707  }
708 
709  if (matchFp == -1)
710  {
712  (
713  "faceCoupleInfo::matchFaces"
714  "(const scalar, const face&, const pointField&"
715  ", const face&, const pointField&)"
716  ) << "No unique match between two faces" << nl
717  << "Face " << f0 << " coords "
718  << UIndirectList<point>(points0, f0)() << nl
719  << "Face " << f1 << " coords "
720  << UIndirectList<point>(points1, f1)()
721  << "when using tolerance " << absTol
722  << " and forwardMatching:" << sameOrientation
723  << abort(FatalError);
724  }
725 
726  return matchFp;
727 }
728 
729 
730 // Find correspondence from patch points to cut points. This might
731 // detect shared points so the output is a patch-to-cut point list
732 // and a compaction list for the cut points (which will always be equal or more
733 // connected than the patch).
734 // Returns true if there are any duplicates.
735 bool Foam::faceCoupleInfo::matchPointsThroughFaces
736 (
737  const scalar absTol,
738  const pointField& cutPoints,
739  const faceList& cutFaces,
740  const pointField& patchPoints,
741  const faceList& patchFaces,
742  const bool sameOrientation,
743 
744  labelList& patchToCutPoints, // patch to (uncompacted) cut points
745  labelList& cutToCompact, // compaction list for cut points
746  labelList& compactToCut // inverse ,,
747 )
748 {
749 
750  // From slave to cut point
751  patchToCutPoints.setSize(patchPoints.size());
752  patchToCutPoints = -1;
753 
754  // Compaction list for cut points: either -1 or index into master which
755  // gives the point to compact to.
756  labelList cutPointRegion(cutPoints.size(), -1);
757  DynamicList<label> cutPointRegionMaster;
758 
759  forAll(patchFaces, patchFaceI)
760  {
761  const face& patchF = patchFaces[patchFaceI];
762 
763  //const face& cutF = cutFaces[patchToCutFaces[patchFaceI]];
764  const face& cutF = cutFaces[patchFaceI];
765 
766  // Do geometric matching to get position of cutF[0] in patchF
767  label patchFp = matchFaces
768  (
769  absTol,
770  patchPoints,
771  patchF,
772  cutPoints,
773  cutF,
774  sameOrientation // orientation
775  );
776 
777  forAll(cutF, cutFp)
778  {
779  label cutPointI = cutF[cutFp];
780  label patchPointI = patchF[patchFp];
781 
782  //const point& cutPt = cutPoints[cutPointI];
783  //const point& patchPt = patchPoints[patchPointI];
784  //if (mag(cutPt - patchPt) > SMALL)
785  //{
786  // FatalErrorIn("matchPointsThroughFaces")
787  // << "cutP:" << cutPt
788  // << " patchP:" << patchPt
789  // << abort(FatalError);
790  //}
791 
792  if (patchToCutPoints[patchPointI] == -1)
793  {
794  patchToCutPoints[patchPointI] = cutPointI;
795  }
796  else if (patchToCutPoints[patchPointI] != cutPointI)
797  {
798  // Multiple cut points connecting to same patch.
799  // Check if already have region & region master for this set
800  label otherCutPointI = patchToCutPoints[patchPointI];
801 
802  //Pout<< "PatchPoint:" << patchPt
803  // << " matches to:" << cutPointI
804  // << " coord:" << cutPoints[cutPointI]
805  // << " and to:" << otherCutPointI
806  // << " coord:" << cutPoints[otherCutPointI]
807  // << endl;
808 
809  if (cutPointRegion[otherCutPointI] != -1)
810  {
811  // Have region for this set. Copy.
812  label region = cutPointRegion[otherCutPointI];
813  cutPointRegion[cutPointI] = region;
814 
815  // Update region master with min point label
816  cutPointRegionMaster[region] = min
817  (
818  cutPointRegionMaster[region],
819  cutPointI
820  );
821  }
822  else
823  {
824  // Create new region.
825  label region = cutPointRegionMaster.size();
826  cutPointRegionMaster.append
827  (
828  min(cutPointI, otherCutPointI)
829  );
830  cutPointRegion[cutPointI] = region;
831  cutPointRegion[otherCutPointI] = region;
832  }
833  }
834 
835  if (sameOrientation)
836  {
837  patchFp = patchF.fcIndex(patchFp);
838  }
839  else
840  {
841  patchFp = patchF.rcIndex(patchFp);
842  }
843  }
844  }
845 
846  // Rework region&master into compaction array
847  compactToCut.setSize(cutPointRegion.size());
848  cutToCompact.setSize(cutPointRegion.size());
849  cutToCompact = -1;
850  label compactPointI = 0;
851 
852  forAll(cutPointRegion, i)
853  {
854  if (cutPointRegion[i] == -1)
855  {
856  // Unduplicated point. Allocate new compacted point.
857  cutToCompact[i] = compactPointI;
858  compactToCut[compactPointI] = i;
859  compactPointI++;
860  }
861  else
862  {
863  // Duplicate point. Get master.
864 
865  label masterPointI = cutPointRegionMaster[cutPointRegion[i]];
866 
867  if (cutToCompact[masterPointI] == -1)
868  {
869  cutToCompact[masterPointI] = compactPointI;
870  compactToCut[compactPointI] = masterPointI;
871  compactPointI++;
872  }
873  cutToCompact[i] = cutToCompact[masterPointI];
874  }
875  }
876  compactToCut.setSize(compactPointI);
877 
878  return compactToCut.size() != cutToCompact.size();
879 }
880 
881 
882 // Return max distance from any point on cutF to masterF
883 Foam::scalar Foam::faceCoupleInfo::maxDistance
884 (
885  const face& cutF,
886  const pointField& cutPoints,
887  const face& masterF,
888  const pointField& masterPoints
889 )
890 {
891  scalar maxDist = -GREAT;
892 
893  forAll(cutF, fp)
894  {
895  const point& cutPt = cutPoints[cutF[fp]];
896 
897  pointHit pHit = masterF.nearestPoint(cutPt, masterPoints);
898 
899  maxDist = max(maxDist, pHit.distance());
900  }
901  return maxDist;
902 }
903 
904 
905 void Foam::faceCoupleInfo::findPerfectMatchingFaces
906 (
907  const primitiveMesh& mesh0,
908  const primitiveMesh& mesh1,
909  const scalar absTol,
910 
911  labelList& mesh0Faces,
912  labelList& mesh1Faces
913 )
914 {
915  // Face centres of external faces (without invoking
916  // mesh.faceCentres since mesh might have been clearedOut)
917 
918  pointField fc0
919  (
920  calcFaceCentres<List>
921  (
922  mesh0.faces(),
923  mesh0.points(),
924  mesh0.nInternalFaces(),
925  mesh0.nFaces() - mesh0.nInternalFaces()
926  )
927  );
928 
929  pointField fc1
930  (
931  calcFaceCentres<List>
932  (
933  mesh1.faces(),
934  mesh1.points(),
935  mesh1.nInternalFaces(),
936  mesh1.nFaces() - mesh1.nInternalFaces()
937  )
938  );
939 
940 
941  if (debug)
942  {
943  Pout<< "Face matching tolerance : " << absTol << endl;
944  }
945 
946 
947  // Match geometrically
948  labelList from1To0;
949  bool matchedAllFaces = matchPoints
950  (
951  fc1,
952  fc0,
953  scalarField(fc1.size(), absTol),
954  false,
955  from1To0
956  );
957 
958  if (matchedAllFaces)
959  {
960  Warning
961  << "faceCoupleInfo::faceCoupleInfo : "
962  << "Matched ALL " << fc1.size()
963  << " boundary faces of mesh0 to boundary faces of mesh1." << endl
964  << "This is only valid if the mesh to add is fully"
965  << " enclosed by the mesh it is added to." << endl;
966  }
967 
968 
969  // Collect matches.
970  label nMatched = 0;
971 
972  mesh0Faces.setSize(fc0.size());
973  mesh1Faces.setSize(fc1.size());
974 
975  forAll(from1To0, i)
976  {
977  if (from1To0[i] != -1)
978  {
979  mesh1Faces[nMatched] = i + mesh1.nInternalFaces();
980  mesh0Faces[nMatched] = from1To0[i] + mesh0.nInternalFaces();
981 
982  nMatched++;
983  }
984  }
985 
986  mesh0Faces.setSize(nMatched);
987  mesh1Faces.setSize(nMatched);
988 }
989 
990 
991 void Foam::faceCoupleInfo::findSlavesCoveringMaster
992 (
993  const primitiveMesh& mesh0,
994  const primitiveMesh& mesh1,
995  const scalar absTol,
996 
997  labelList& mesh0Faces,
998  labelList& mesh1Faces
999 )
1000 {
1001  // Construct octree from all mesh0 boundary faces
1002  labelList bndFaces(mesh0.nFaces()-mesh0.nInternalFaces());
1003  forAll(bndFaces, i)
1004  {
1005  bndFaces[i] = mesh0.nInternalFaces() + i;
1006  }
1007 
1008  treeBoundBox overallBb(mesh0.points());
1009 
1010  Random rndGen(123456);
1011 
1012  indexedOctree<treeDataFace> tree
1013  (
1014  treeDataFace // all information needed to search faces
1015  (
1016  false, // do not cache bb
1017  mesh0,
1018  bndFaces // boundary faces only
1019  ),
1020  overallBb.extend(rndGen, 1e-4), // overall search domain
1021  8, // maxLevel
1022  10, // leafsize
1023  3.0 // duplicity
1024  );
1025 
1026  if (debug)
1027  {
1028  Pout<< "findSlavesCoveringMaster :"
1029  << " constructed octree for mesh0 boundary faces" << endl;
1030  }
1031 
1032  // Search nearest face for every mesh1 boundary face.
1033 
1034  labelHashSet mesh0Set(mesh0.nFaces() - mesh0.nInternalFaces());
1035  labelHashSet mesh1Set(mesh1.nFaces() - mesh1.nInternalFaces());
1036 
1037  for
1038  (
1039  label mesh1FaceI = mesh1.nInternalFaces();
1040  mesh1FaceI < mesh1.nFaces();
1041  mesh1FaceI++
1042  )
1043  {
1044  const face& f1 = mesh1.faces()[mesh1FaceI];
1045 
1046  // Generate face centre (prevent cellCentres() reconstruction)
1047  point fc(f1.centre(mesh1.points()));
1048 
1049  pointIndexHit nearInfo = tree.findNearest(fc, Foam::sqr(absTol));
1050 
1051  if (nearInfo.hit())
1052  {
1053  label mesh0FaceI = tree.shapes().faceLabels()[nearInfo.index()];
1054 
1055  // Check if points of f1 actually lie on top of mesh0 face
1056  // This is the bit that might fail since if f0 is severely warped
1057  // and f1's points are not present on f0 (since subdivision)
1058  // then the distance of the points to f0 might be quite large
1059  // and the test will fail. This all should in fact be some kind
1060  // of walk from known corresponding points/edges/faces.
1061  scalar dist =
1062  maxDistance
1063  (
1064  f1,
1065  mesh1.points(),
1066  mesh0.faces()[mesh0FaceI],
1067  mesh0.points()
1068  );
1069 
1070  if (dist < absTol)
1071  {
1072  mesh0Set.insert(mesh0FaceI);
1073  mesh1Set.insert(mesh1FaceI);
1074  }
1075  }
1076  }
1077 
1078  if (debug)
1079  {
1080  Pout<< "findSlavesCoveringMaster :"
1081  << " matched " << mesh1Set.size() << " mesh1 faces to "
1082  << mesh0Set.size() << " mesh0 faces" << endl;
1083  }
1084 
1085  mesh0Faces = mesh0Set.toc();
1086  mesh1Faces = mesh1Set.toc();
1087 }
1088 
1089 
1090 // Grow cutToMasterFace across 'internal' edges.
1091 Foam::label Foam::faceCoupleInfo::growCutFaces
1092 (
1093  const labelList& cutToMasterEdges,
1094  Map<labelList>& candidates
1095 )
1096 {
1097  if (debug)
1098  {
1099  Pout<< "growCutFaces :"
1100  << " growing cut faces to masterPatch" << endl;
1101  }
1102 
1103  label nTotChanged = 0;
1104 
1105  while (true)
1106  {
1107  const labelListList& cutFaceEdges = cutFaces().faceEdges();
1108  const labelListList& cutEdgeFaces = cutFaces().edgeFaces();
1109 
1110  label nChanged = 0;
1111 
1112  forAll(cutToMasterFaces_, cutFaceI)
1113  {
1114  const label masterFaceI = cutToMasterFaces_[cutFaceI];
1115 
1116  if (masterFaceI != -1)
1117  {
1118  // We now have a cutFace for which we have already found a
1119  // master face. Grow this masterface across any internal edge
1120  // (internal: no corresponding master edge)
1121 
1122  const labelList& fEdges = cutFaceEdges[cutFaceI];
1123 
1124  forAll(fEdges, i)
1125  {
1126  const label cutEdgeI = fEdges[i];
1127 
1128  if (cutToMasterEdges[cutEdgeI] == -1)
1129  {
1130  // So this edge:
1131  // - internal to the cutPatch (since all region edges
1132  // marked before)
1133  // - on cutFaceI which corresponds to masterFace.
1134  // Mark all connected faces with this masterFace.
1135 
1136  const labelList& eFaces = cutEdgeFaces[cutEdgeI];
1137 
1138  forAll(eFaces, j)
1139  {
1140  const label faceI = eFaces[j];
1141 
1142  if (cutToMasterFaces_[faceI] == -1)
1143  {
1144  cutToMasterFaces_[faceI] = masterFaceI;
1145  candidates.erase(faceI);
1146  nChanged++;
1147  }
1148  else if (cutToMasterFaces_[faceI] != masterFaceI)
1149  {
1150  const pointField& cutPoints =
1151  cutFaces().points();
1152  const pointField& masterPoints =
1153  masterPatch().points();
1154 
1155  const edge& e = cutFaces().edges()[cutEdgeI];
1156 
1157  label myMaster = cutToMasterFaces_[faceI];
1158  const face& myF = masterPatch()[myMaster];
1159 
1160  const face& nbrF = masterPatch()[masterFaceI];
1161 
1162  FatalErrorIn
1163  (
1164  "faceCoupleInfo::growCutFaces"
1165  "(const labelList&, Map<labelList>&)"
1166  ) << "Cut face "
1167  << cutFaces()[faceI].points(cutPoints)
1168  << " has master "
1169  << myMaster
1170  << " but also connects to nbr face "
1171  << cutFaces()[cutFaceI].points(cutPoints)
1172  << " with master " << masterFaceI
1173  << nl
1174  << "myMasterFace:"
1175  << myF.points(masterPoints)
1176  << " nbrMasterFace:"
1177  << nbrF.points(masterPoints) << nl
1178  << "Across cut edge " << cutEdgeI
1179  << " coord:"
1180  << cutFaces().localPoints()[e[0]]
1181  << cutFaces().localPoints()[e[1]]
1182  << abort(FatalError);
1183  }
1184  }
1185  }
1186  }
1187  }
1188  }
1189 
1190  if (debug)
1191  {
1192  Pout<< "growCutFaces : Grown an additional " << nChanged
1193  << " cut to master face" << " correspondence" << endl;
1194  }
1195 
1196  nTotChanged += nChanged;
1197 
1198  if (nChanged == 0)
1199  {
1200  break;
1201  }
1202  }
1203 
1204  return nTotChanged;
1205 }
1206 
1207 
1208 void Foam::faceCoupleInfo::checkMatch(const labelList& cutToMasterEdges) const
1209 {
1210  const pointField& cutLocalPoints = cutFaces().localPoints();
1211 
1212  const pointField& masterLocalPoints = masterPatch().localPoints();
1213  const faceList& masterLocalFaces = masterPatch().localFaces();
1214 
1215  forAll(cutToMasterEdges, cutEdgeI)
1216  {
1217  const edge& e = cutFaces().edges()[cutEdgeI];
1218 
1219  if (cutToMasterEdges[cutEdgeI] == -1)
1220  {
1221  // Internal edge. Check that master face is same on both sides.
1222  const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1223 
1224  label masterFaceI = -1;
1225 
1226  forAll(cutEFaces, i)
1227  {
1228  label cutFaceI = cutEFaces[i];
1229 
1230  if (cutToMasterFaces_[cutFaceI] != -1)
1231  {
1232  if (masterFaceI == -1)
1233  {
1234  masterFaceI = cutToMasterFaces_[cutFaceI];
1235  }
1236  else if (masterFaceI != cutToMasterFaces_[cutFaceI])
1237  {
1238  label myMaster = cutToMasterFaces_[cutFaceI];
1239  const face& myF = masterLocalFaces[myMaster];
1240 
1241  const face& nbrF = masterLocalFaces[masterFaceI];
1242 
1243  FatalErrorIn
1244  (
1245  "faceCoupleInfo::checkMatch(const labelList&) const"
1246  )
1247  << "Internal CutEdge " << e
1248  << " coord:"
1249  << cutLocalPoints[e[0]]
1250  << cutLocalPoints[e[1]]
1251  << " connects to master " << myMaster
1252  << " and to master " << masterFaceI << nl
1253  << "myMasterFace:"
1254  << myF.points(masterLocalPoints)
1255  << " nbrMasterFace:"
1256  << nbrF.points(masterLocalPoints)
1257  << abort(FatalError);
1258  }
1259  }
1260  }
1261  }
1262  }
1263 }
1264 
1265 
1266 // Extends matching information by elimination across cutFaces using more
1267 // than one region edge. Updates cutToMasterFaces_ and sets candidates
1268 // which is for every cutface on a region edge the possible master faces.
1269 Foam::label Foam::faceCoupleInfo::matchEdgeFaces
1270 (
1271  const labelList& cutToMasterEdges,
1272  Map<labelList>& candidates
1273 )
1274 {
1275  // For every unassigned cutFaceI the possible list of master faces.
1276  candidates.clear();
1277  candidates.resize(cutFaces().size());
1278 
1279  label nChanged = 0;
1280 
1281  forAll(cutToMasterEdges, cutEdgeI)
1282  {
1283  label masterEdgeI = cutToMasterEdges[cutEdgeI];
1284 
1285  if (masterEdgeI != -1)
1286  {
1287  // So cutEdgeI is matched to masterEdgeI. For all cut faces using
1288  // the cut edge store the master face as a candidate match.
1289 
1290  const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1291  const labelList& masterEFaces =
1292  masterPatch().edgeFaces()[masterEdgeI];
1293 
1294  forAll(cutEFaces, i)
1295  {
1296  label cutFaceI = cutEFaces[i];
1297 
1298  if (cutToMasterFaces_[cutFaceI] == -1)
1299  {
1300  // So this face (cutFaceI) is on an edge (cutEdgeI) that
1301  // is used by some master faces (masterEFaces). Check if
1302  // the cutFaceI and some of these masterEFaces share more
1303  // than one edge (which uniquely defines face)
1304 
1305  // Combine master faces with current set of candidate
1306  // master faces.
1307  Map<labelList>::iterator fnd = candidates.find(cutFaceI);
1308 
1309  if (fnd == candidates.end())
1310  {
1311  // No info yet for cutFaceI. Add all master faces as
1312  // candidates
1313  candidates.insert(cutFaceI, masterEFaces);
1314  }
1315  else
1316  {
1317  // From some other cutEdgeI there are already some
1318  // candidate master faces. Check the overlap with
1319  // the current set of master faces.
1320  const labelList& masterFaces = fnd();
1321 
1322  DynamicList<label> newCandidates(masterFaces.size());
1323 
1324  forAll(masterEFaces, j)
1325  {
1326  if (findIndex(masterFaces, masterEFaces[j]) != -1)
1327  {
1328  newCandidates.append(masterEFaces[j]);
1329  }
1330  }
1331 
1332  if (newCandidates.size() == 1)
1333  {
1334  // We found a perfect match. Delete entry from
1335  // candidates map.
1336  cutToMasterFaces_[cutFaceI] = newCandidates[0];
1337  candidates.erase(cutFaceI);
1338  nChanged++;
1339  }
1340  else
1341  {
1342  // Should not happen?
1343  //Pout<< "On edge:" << cutEdgeI
1344  // << " have connected masterFaces:"
1345  // << masterEFaces
1346  // << " and from previous edge we have"
1347  // << " connected masterFaces" << masterFaces
1348  // << " . Overlap:" << newCandidates << endl;
1349 
1350  fnd() = newCandidates.shrink();
1351  }
1352  }
1353  }
1354 
1355  }
1356  }
1357  }
1358 
1359  if (debug)
1360  {
1361  Pout<< "matchEdgeFaces : Found " << nChanged
1362  << " faces where there was"
1363  << " only one remaining choice for cut-master correspondence"
1364  << endl;
1365  }
1366 
1367  return nChanged;
1368 }
1369 
1370 
1371 // Gets a list of cutFaces (that use a master edge) and the candidate
1372 // master faces.
1373 // Finds most aligned master face.
1374 Foam::label Foam::faceCoupleInfo::geometricMatchEdgeFaces
1375 (
1376  Map<labelList>& candidates
1377 )
1378 {
1379  const pointField& cutPoints = cutFaces().points();
1380 
1381  label nChanged = 0;
1382 
1383  // Mark all master faces that have been matched so far.
1384 
1385  labelListList masterToCutFaces
1386  (
1388  (
1389  masterPatch().size(),
1390  cutToMasterFaces_
1391  )
1392  );
1393 
1394  forAllConstIter(Map<labelList>, candidates, iter)
1395  {
1396  label cutFaceI = iter.key();
1397 
1398  const face& cutF = cutFaces()[cutFaceI];
1399 
1400  if (cutToMasterFaces_[cutFaceI] == -1)
1401  {
1402  const labelList& masterFaces = iter();
1403 
1404  // Find the best matching master face.
1405  scalar minDist = GREAT;
1406  label minMasterFaceI = -1;
1407 
1408  forAll(masterFaces, i)
1409  {
1410  label masterFaceI = masterFaces[i];
1411 
1412  if (masterToCutFaces[masterFaces[i]].empty())
1413  {
1414  scalar dist = maxDistance
1415  (
1416  cutF,
1417  cutPoints,
1418  masterPatch()[masterFaceI],
1419  masterPatch().points()
1420  );
1421 
1422  if (dist < minDist)
1423  {
1424  minDist = dist;
1425  minMasterFaceI = masterFaceI;
1426  }
1427  }
1428  }
1429 
1430  if (minMasterFaceI != -1)
1431  {
1432  cutToMasterFaces_[cutFaceI] = minMasterFaceI;
1433  masterToCutFaces[minMasterFaceI] = cutFaceI;
1434  nChanged++;
1435  }
1436  }
1437  }
1438 
1439  // (inefficiently) force candidates to be uptodate.
1440  forAll(cutToMasterFaces_, cutFaceI)
1441  {
1442  if (cutToMasterFaces_[cutFaceI] != -1)
1443  {
1444  candidates.erase(cutFaceI);
1445  }
1446  }
1447 
1448 
1449  if (debug)
1450  {
1451  Pout<< "geometricMatchEdgeFaces : Found " << nChanged
1452  << " faces where there was"
1453  << " only one remaining choice for cut-master correspondence"
1454  << endl;
1455  }
1456 
1457  return nChanged;
1458 }
1459 
1460 
1461 // Calculate the set of cut faces inbetween master and slave patch
1462 // assuming perfect match (and optional face ordering on slave)
1463 void Foam::faceCoupleInfo::perfectPointMatch
1464 (
1465  const scalar absTol,
1466  const bool slaveFacesOrdered
1467 )
1468 {
1469  if (debug)
1470  {
1471  Pout<< "perfectPointMatch :"
1472  << " Matching master and slave to cut."
1473  << " Master and slave faces are identical" << nl;
1474 
1475  if (slaveFacesOrdered)
1476  {
1477  Pout<< "and master and slave faces are ordered"
1478  << " (on coupled patches)" << endl;
1479  }
1480  else
1481  {
1482  Pout<< "and master and slave faces are not ordered"
1483  << " (on coupled patches)" << endl;
1484  }
1485  }
1486 
1487  cutToMasterFaces_ = identity(masterPatch().size());
1488  cutPoints_ = masterPatch().localPoints();
1489  cutFacesPtr_.reset
1490  (
1491  new primitiveFacePatch
1492  (
1493  masterPatch().localFaces(),
1494  cutPoints_
1495  )
1496  );
1497  masterToCutPoints_ = identity(cutPoints_.size());
1498 
1499 
1500  // Cut faces to slave patch.
1501  bool matchedAllFaces = false;
1502 
1503  if (slaveFacesOrdered)
1504  {
1505  cutToSlaveFaces_ = identity(cutFaces().size());
1506  matchedAllFaces = (cutFaces().size() == slavePatch().size());
1507  }
1508  else
1509  {
1510  // Faces do not have to be ordered (but all have
1511  // to match). Note: Faces will be already ordered if we enter here from
1512  // construct from meshes.
1513  matchedAllFaces = matchPoints
1514  (
1515  calcFaceCentres<List>
1516  (
1517  cutFaces(),
1518  cutPoints_,
1519  0,
1520  cutFaces().size()
1521  ),
1522  calcFaceCentres<IndirectList>
1523  (
1524  slavePatch(),
1525  slavePatch().points(),
1526  0,
1527  slavePatch().size()
1528  ),
1529  scalarField(slavePatch().size(), absTol),
1530  true,
1531  cutToSlaveFaces_
1532  );
1533  }
1534 
1535 
1536  if (!matchedAllFaces)
1537  {
1538  FatalErrorIn
1539  (
1540  "faceCoupleInfo::perfectPointMatch"
1541  "(const scalar, const bool)"
1542  ) << "Did not match all of the master faces to the slave faces"
1543  << endl
1544  << "This usually means that the slave patch and master patch"
1545  << " do not align to within " << absTol << " metre."
1546  << abort(FatalError);
1547  }
1548 
1549 
1550  // Find correspondence from slave points to cut points. This might
1551  // detect shared points so the output is a slave-to-cut point list
1552  // and a compaction list.
1553 
1554  labelList cutToCompact, compactToCut;
1555  matchPointsThroughFaces
1556  (
1557  absTol,
1558  cutFaces().localPoints(),
1559  reorder(cutToSlaveFaces_, cutFaces().localFaces()),
1560  slavePatch().localPoints(),
1561  slavePatch().localFaces(),
1562  false, // slave and cut have opposite orientation
1563 
1564  slaveToCutPoints_, // slave to (uncompacted) cut points
1565  cutToCompact, // compaction map: from cut to compacted
1566  compactToCut // compaction map: from compacted to cut
1567  );
1568 
1569 
1570  // Use compaction lists to renumber cutPoints.
1571  cutPoints_ = UIndirectList<point>(cutPoints_, compactToCut)();
1572  {
1573  const faceList& cutLocalFaces = cutFaces().localFaces();
1574 
1575  faceList compactFaces(cutLocalFaces.size());
1576  forAll(cutLocalFaces, i)
1577  {
1578  compactFaces[i] = renumber(cutToCompact, cutLocalFaces[i]);
1579  }
1580  cutFacesPtr_.reset
1581  (
1582  new primitiveFacePatch
1583  (
1584  compactFaces,
1585  cutPoints_
1586  )
1587  );
1588  }
1589  inplaceRenumber(cutToCompact, slaveToCutPoints_);
1590  inplaceRenumber(cutToCompact, masterToCutPoints_);
1591 }
1592 
1593 
1594 // Calculate the set of cut faces inbetween master and slave patch
1595 // assuming that slave patch is subdivision of masterPatch.
1596 void Foam::faceCoupleInfo::subDivisionMatch
1597 (
1598  const polyMesh& slaveMesh,
1599  const bool patchDivision,
1600  const scalar absTol
1601 )
1602 {
1603  if (debug)
1604  {
1605  Pout<< "subDivisionMatch :"
1606  << " Matching master and slave to cut."
1607  << " Slave can be subdivision of master but all master edges"
1608  << " have to be covered by slave edges." << endl;
1609  }
1610 
1611  // Assume slave patch is subdivision of the master patch so cutFaces is a
1612  // copy of the slave (but reversed since orientation has to be like master).
1613  cutPoints_ = slavePatch().localPoints();
1614  {
1615  faceList cutFaces(slavePatch().size());
1616 
1617  forAll(cutFaces, i)
1618  {
1619  cutFaces[i] = slavePatch().localFaces()[i].reverseFace();
1620  }
1621  cutFacesPtr_.reset(new primitiveFacePatch(cutFaces, cutPoints_));
1622  }
1623 
1624  // Cut is copy of slave so addressing to slave is trivial.
1625  cutToSlaveFaces_ = identity(cutFaces().size());
1626  slaveToCutPoints_ = identity(slavePatch().nPoints());
1627 
1628  // Cut edges to slave patch
1629  labelList cutToSlaveEdges
1630  (
1631  findMappedEdges
1632  (
1633  cutFaces().edges(),
1634  slaveToCutPoints_, //note:should be cutToSlavePoints but since iden
1635  slavePatch()
1636  )
1637  );
1638 
1639 
1640  if (debug)
1641  {
1642  OFstream str("regionEdges.obj");
1643 
1644  Pout<< "subDivisionMatch :"
1645  << " addressing for slave patch fully done."
1646  << " Dumping region edges to " << str.name() << endl;
1647 
1648  label vertI = 0;
1649 
1650  forAll(slavePatch().edges(), slaveEdgeI)
1651  {
1652  if (regionEdge(slaveMesh, slaveEdgeI))
1653  {
1654  const edge& e = slavePatch().edges()[slaveEdgeI];
1655 
1656  meshTools::writeOBJ(str, slavePatch().localPoints()[e[0]]);
1657  vertI++;
1658  meshTools::writeOBJ(str, slavePatch().localPoints()[e[1]]);
1659  vertI++;
1660  str<< "l " << vertI-1 << ' ' << vertI << nl;
1661  }
1662  }
1663  }
1664 
1665 
1666  // Addressing from cut to master.
1667 
1668  // Cut points to master patch
1669  // - determine master points to cut points. Has to be full!
1670  // - invert to get cut to master
1671  if (debug)
1672  {
1673  Pout<< "subDivisionMatch :"
1674  << " matching master points to cut points" << endl;
1675  }
1676 
1677  bool matchedAllPoints = matchPoints
1678  (
1679  masterPatch().localPoints(),
1680  cutPoints_,
1681  scalarField(masterPatch().nPoints(), absTol),
1682  true,
1683  masterToCutPoints_
1684  );
1685 
1686  if (!matchedAllPoints)
1687  {
1688  FatalErrorIn
1689  (
1690  "faceCoupleInfo::subDivisionMatch"
1691  "(const polyMesh&, const bool, const scalar)"
1692  ) << "Did not match all of the master points to the slave points"
1693  << endl
1694  << "This usually means that the slave patch is not a"
1695  << " subdivision of the master patch"
1696  << abort(FatalError);
1697  }
1698 
1699 
1700  // Do masterEdges to cutEdges :
1701  // - mark all edges between two masterEdge endpoints. (geometric test since
1702  // this is the only distinction)
1703  // - this gives the borders inbetween which all cutfaces come from
1704  // a single master face.
1705  if (debug)
1706  {
1707  Pout<< "subDivisionMatch :"
1708  << " matching cut edges to master edges" << endl;
1709  }
1710 
1711  const edgeList& masterEdges = masterPatch().edges();
1712  const pointField& masterPoints = masterPatch().localPoints();
1713 
1714  const edgeList& cutEdges = cutFaces().edges();
1715 
1716  labelList cutToMasterEdges(cutFaces().nEdges(), -1);
1717 
1718  forAll(masterEdges, masterEdgeI)
1719  {
1720  const edge& masterEdge = masterEdges[masterEdgeI];
1721 
1722  label cutPoint0 = masterToCutPoints_[masterEdge[0]];
1723  label cutPoint1 = masterToCutPoints_[masterEdge[1]];
1724 
1725  // Find edges between cutPoint0 and cutPoint1.
1726 
1727  label cutPointI = cutPoint0;
1728  do
1729  {
1730  // Find edge (starting at pointI on cut), aligned with master
1731  // edge.
1732  label cutEdgeI =
1733  mostAlignedCutEdge
1734  (
1735  false,
1736  slaveMesh,
1737  patchDivision,
1738  cutToMasterEdges,
1739  cutToSlaveEdges,
1740  cutPointI,
1741  cutPoint0,
1742  cutPoint1
1743  );
1744 
1745  if (cutEdgeI == -1)
1746  {
1747  // Problem. Report while matching to produce nice error message
1748  mostAlignedCutEdge
1749  (
1750  true,
1751  slaveMesh,
1752  patchDivision,
1753  cutToMasterEdges,
1754  cutToSlaveEdges,
1755  cutPointI,
1756  cutPoint0,
1757  cutPoint1
1758  );
1759 
1760  Pout<< "Dumping unmatched pointEdges to errorEdges.obj"
1761  << endl;
1762 
1763  writeOBJ
1764  (
1765  "errorEdges.obj",
1766  edgeList
1767  (
1768  UIndirectList<edge>
1769  (
1770  cutFaces().edges(),
1771  cutFaces().pointEdges()[cutPointI]
1772  )
1773  ),
1774  cutFaces().localPoints(),
1775  false
1776  );
1777 
1778  FatalErrorIn
1779  (
1780  "faceCoupleInfo::subDivisionMatch"
1781  "(const polyMesh&, const bool, const scalar)"
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
1788  << abort(FatalError);
1789  }
1790 
1791  cutToMasterEdges[cutEdgeI] = masterEdgeI;
1792 
1793  cutPointI = cutEdges[cutEdgeI].otherVertex(cutPointI);
1794 
1795  } while (cutPointI != cutPoint1);
1796  }
1797 
1798  if (debug)
1799  {
1800  // Write all matched edges.
1801  writeEdges(cutToMasterEdges, cutToSlaveEdges);
1802  }
1803 
1804  // Rework cutToMasterEdges into list of points inbetween two endpoints
1805  // (cutEdgeToPoints_)
1806  setCutEdgeToPoints(cutToMasterEdges);
1807 
1808 
1809  // Now that we have marked the cut edges that lie on top of master edges
1810  // we can find cut faces that have two (or more) of these edges.
1811  // Note that there can be multiple faces having two or more matched edges
1812  // since the cut faces can form a non-manifold surface(?)
1813  // So the matching is done as an elimination process where for every
1814  // cutFace on a matched edge we store the possible master faces and
1815  // eliminate more and more until we only have one possible master face
1816  // left.
1817 
1818  if (debug)
1819  {
1820  Pout<< "subDivisionMatch :"
1821  << " matching (topological) cut faces to masterPatch" << endl;
1822  }
1823 
1824  // For every unassigned cutFaceI the possible list of master faces.
1825  Map<labelList> candidates(cutFaces().size());
1826 
1827  cutToMasterFaces_.setSize(cutFaces().size());
1828  cutToMasterFaces_ = -1;
1829 
1830  while (true)
1831  {
1832  // See if there are any edges left where there is only one remaining
1833  // candidate.
1834  label nChanged = matchEdgeFaces(cutToMasterEdges, candidates);
1835 
1836  checkMatch(cutToMasterEdges);
1837 
1838 
1839  // Now we should have found a face correspondence for every cutFace
1840  // that uses two or more cutEdges. Floodfill from these across
1841  // 'internal' cutedges (i.e. edges that do not have a master
1842  // equivalent) to determine some more cutToMasterFaces
1843  nChanged += growCutFaces(cutToMasterEdges, candidates);
1844 
1845  checkMatch(cutToMasterEdges);
1846 
1847  if (nChanged == 0)
1848  {
1849  break;
1850  }
1851  }
1852 
1853  if (debug)
1854  {
1855  Pout<< "subDivisionMatch :"
1856  << " matching (geometric) cut faces to masterPatch" << endl;
1857  }
1858 
1859  // Above should have matched all cases fully. Cannot prove this yet so
1860  // do any remaining unmatched edges by a geometric test.
1861  while (true)
1862  {
1863  // See if there are any edges left where there is only one remaining
1864  // candidate.
1865  label nChanged = geometricMatchEdgeFaces(candidates);
1866 
1867  checkMatch(cutToMasterEdges);
1868 
1869  nChanged += growCutFaces(cutToMasterEdges, candidates);
1870 
1871  checkMatch(cutToMasterEdges);
1872 
1873  if (nChanged == 0)
1874  {
1875  break;
1876  }
1877  }
1878 
1879 
1880  // All cut faces matched?
1881  forAll(cutToMasterFaces_, cutFaceI)
1882  {
1883  if (cutToMasterFaces_[cutFaceI] == -1)
1884  {
1885  const face& cutF = cutFaces()[cutFaceI];
1886 
1887  FatalErrorIn
1888  (
1889  "faceCoupleInfo::subDivisionMatch"
1890  "(const polyMesh&, const bool, const scalar)"
1891  ) << "Did not match all of cutFaces to a master face" << nl
1892  << "First unmatched cut face:" << cutFaceI << " with points:"
1893  << UIndirectList<point>(cutFaces().points(), cutF)()
1894  << nl
1895  << "This usually means that the slave patch is not a"
1896  << " subdivision of the master patch"
1897  << abort(FatalError);
1898  }
1899  }
1900 
1901  if (debug)
1902  {
1903  Pout<< "subDivisionMatch :"
1904  << " finished matching master and slave to cut" << endl;
1905  }
1906 
1907  if (debug)
1908  {
1909  writePointsFaces();
1910  }
1911 }
1912 
1913 
1914 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1915 
1916 // Construct from mesh data
1919  const polyMesh& masterMesh,
1920  const polyMesh& slaveMesh,
1921  const scalar absTol,
1922  const bool perfectMatch
1923 )
1924 :
1925  masterPatchPtr_(NULL),
1926  slavePatchPtr_(NULL),
1927  cutPoints_(0),
1928  cutFacesPtr_(NULL),
1929  cutToMasterFaces_(0),
1930  masterToCutPoints_(0),
1931  cutToSlaveFaces_(0),
1932  slaveToCutPoints_(0),
1933  cutEdgeToPoints_(0)
1934 {
1935  // Get faces on both meshes that are aligned.
1936  // (not ordered i.e. masterToMesh[0] does
1937  // not couple to slaveToMesh[0]. This ordering is done later on)
1938  labelList masterToMesh;
1939  labelList slaveToMesh;
1940 
1941  if (perfectMatch)
1942  {
1943  // Identical faces so use tight face-centre comparison.
1944  findPerfectMatchingFaces
1945  (
1946  masterMesh,
1947  slaveMesh,
1948  absTol,
1949  masterToMesh,
1950  slaveToMesh
1951  );
1952  }
1953  else
1954  {
1955  // Slave subdivision of master so use 'nearest'. Bit dodgy, especially
1956  // with using absTol (which is quite small)
1957  findSlavesCoveringMaster
1958  (
1959  masterMesh,
1960  slaveMesh,
1961  absTol,
1962  masterToMesh,
1963  slaveToMesh
1964  );
1965  }
1966 
1967  // Construct addressing engines for both sides
1968  masterPatchPtr_.reset
1969  (
1971  (
1972  IndirectList<face>(masterMesh.faces(), masterToMesh),
1973  masterMesh.points()
1974  )
1975  );
1976 
1977  slavePatchPtr_.reset
1978  (
1980  (
1981  IndirectList<face>(slaveMesh.faces(), slaveToMesh),
1982  slaveMesh.points()
1983  )
1984  );
1985 
1986 
1987  if (perfectMatch)
1988  {
1989  // Faces are perfectly aligned but probably not ordered.
1990  perfectPointMatch(absTol, false);
1991  }
1992  else
1993  {
1994  // Slave faces are subdivision of master face. Faces not ordered.
1995  subDivisionMatch(slaveMesh, false, absTol);
1996  }
1997 
1998  if (debug)
1999  {
2000  writePointsFaces();
2001  }
2002 }
2003 
2004 
2005 // Slave is subdivision of master patch.
2006 // (so -both cover the same area -all of master points are present in slave)
2009  const polyMesh& masterMesh,
2010  const labelList& masterAddressing,
2011  const polyMesh& slaveMesh,
2012  const labelList& slaveAddressing,
2013  const scalar absTol,
2014  const bool perfectMatch,
2015  const bool orderedFaces,
2016  const bool patchDivision
2017 )
2018 :
2019  masterPatchPtr_
2020  (
2022  (
2023  IndirectList<face>(masterMesh.faces(), masterAddressing),
2024  masterMesh.points()
2025  )
2026  ),
2027  slavePatchPtr_
2028  (
2030  (
2031  IndirectList<face>(slaveMesh.faces(), slaveAddressing),
2032  slaveMesh.points()
2033  )
2034  ),
2035  cutPoints_(0),
2036  cutFacesPtr_(NULL),
2037  cutToMasterFaces_(0),
2038  masterToCutPoints_(0),
2039  cutToSlaveFaces_(0),
2040  slaveToCutPoints_(0),
2041  cutEdgeToPoints_(0)
2042 {
2043  if (perfectMatch && (masterAddressing.size() != slaveAddressing.size()))
2044  {
2045  FatalErrorIn
2046  (
2047  "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2048  ", const primitiveMesh&, const scalar, const bool"
2049  ) << "Perfect match specified but number of master and slave faces"
2050  << " differ." << endl
2051  << "master:" << masterAddressing.size()
2052  << " slave:" << slaveAddressing.size()
2053  << abort(FatalError);
2054  }
2055 
2056  if
2057  (
2058  masterAddressing.size()
2059  && min(masterAddressing) < masterMesh.nInternalFaces()
2060  )
2061  {
2062  FatalErrorIn
2063  (
2064  "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2065  ", const primitiveMesh&, const scalar, const bool"
2066  ) << "Supplied internal face on master mesh to couple." << nl
2067  << "Faces to be coupled have to be boundary faces."
2068  << abort(FatalError);
2069  }
2070  if
2071  (
2072  slaveAddressing.size()
2073  && min(slaveAddressing) < slaveMesh.nInternalFaces()
2074  )
2075  {
2076  FatalErrorIn
2077  (
2078  "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2079  ", const primitiveMesh&, const scalar, const bool"
2080  ) << "Supplied internal face on slave mesh to couple." << nl
2081  << "Faces to be coupled have to be boundary faces."
2082  << abort(FatalError);
2083  }
2084 
2085 
2086  if (perfectMatch)
2087  {
2088  perfectPointMatch(absTol, orderedFaces);
2089  }
2090  else
2091  {
2092  // Slave faces are subdivision of master face. Faces not ordered.
2093  subDivisionMatch(slaveMesh, patchDivision, absTol);
2094  }
2095 
2096  if (debug)
2097  {
2098  writePointsFaces();
2099  }
2100 }
2101 
2102 
2103 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
2104 
2106 {}
2107 
2108 
2109 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2110 
2112 {
2113  labelList faces(pp.size());
2114 
2115  label faceI = pp.start();
2116 
2117  forAll(pp, i)
2118  {
2119  faces[i] = faceI++;
2120  }
2121  return faces;
2122 }
2123 
2124 
2126 {
2127  Map<label> map(lst.size());
2128 
2129  forAll(lst, i)
2130  {
2131  if (lst[i] != -1)
2132  {
2133  map.insert(i, lst[i]);
2134  }
2135  }
2136  return map;
2137 }
2138 
2139 
2142  const labelListList& lst
2143 )
2144 {
2145  Map<labelList> map(lst.size());
2146 
2147  forAll(lst, i)
2148  {
2149  if (lst[i].size())
2150  {
2151  map.insert(i, lst[i]);
2152  }
2153  }
2154  return map;
2155 }
2156 
2157 
2158 // ************************************************************************* //
A List with indirect addressing.
Definition: IndirectList.H:102
const pointField & points
static labelList faceLabels(const polyPatch &)
Utility functions.
cachedRandom rndGen(label(0),-1)
vector point
Point is a vector.
Definition: point.H:41
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
dimensioned< scalar > mag(const dimensioned< Type > &)
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
dimensioned< scalar > magSqr(const dimensioned< Type > &)
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49
scalar f1
Definition: createFields.H:28
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
List< face > faceList
Definition: faceListFwd.H:43
PrimitivePatch< face, List, const pointField & > primitiveFacePatch
Foam::primitiveFacePatch.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
PointHit< point > pointHit
Definition: pointHit.H:41
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
messageStream Warning
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
ListType renumber(const labelUList &oldToNew, const ListType &)
Renumber the values (not the indices) of a list.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
faceCoupleInfo(const polyMesh &mesh0, const polyMesh &mesh1, const scalar absTol, const bool perfectMatch)
Construct from two meshes and absolute tolerance.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:88
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
List< edge > edgeList
Definition: edgeList.H:38
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
A list of faces which address into the list of points.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
label nInternalFaces() const
~faceCoupleInfo()
Destructor.
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:50
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Determine correspondence between points. See below.
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const Vector zero
Definition: Vector.H:80
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.
Definition: matchPoints.C:33
HashTable< T, label, Hash< label > >::iterator iterator
Definition: Map.H:56
static Map< label > makeMap(const labelList &)
Create Map from List.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1060
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:139
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
label nPoints
label findEdge(const edgeList &edges, const labelList &candidates, const label v0, const label v1)
Return edge among candidates that uses the two vertices.
Definition: meshTools.C:343
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:209
ListType reorder(const labelUList &oldToNew, const ListType &)
Reorder the elements (indices, not values) of a list.
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53