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