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