All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
faceCoupleInfo.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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  // in between 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  indexedOctree<treeDataFace> tree
984  (
985  treeDataFace // all information needed to search faces
986  (
987  false, // do not cache bb
988  mesh0,
989  bndFaces // boundary faces only
990  ),
991  overallBb.extend(1e-4), // overall search domain
992  8, // maxLevel
993  10, // leafsize
994  3.0 // duplicity
995  );
996 
997  if (debug)
998  {
999  Pout<< "findSlavesCoveringMaster :"
1000  << " constructed octree for mesh0 boundary faces" << endl;
1001  }
1002 
1003  // Search nearest face for every mesh1 boundary face.
1004 
1005  labelHashSet mesh0Set(mesh0.nFaces() - mesh0.nInternalFaces());
1006  labelHashSet mesh1Set(mesh1.nFaces() - mesh1.nInternalFaces());
1007 
1008  for
1009  (
1010  label mesh1Facei = mesh1.nInternalFaces();
1011  mesh1Facei < mesh1.nFaces();
1012  mesh1Facei++
1013  )
1014  {
1015  const face& f1 = mesh1.faces()[mesh1Facei];
1016 
1017  // Generate face centre (prevent cellCentres() reconstruction)
1018  point fc(f1.centre(mesh1.points()));
1019 
1020  pointIndexHit nearInfo = tree.findNearest(fc, Foam::sqr(absTol));
1021 
1022  if (nearInfo.hit())
1023  {
1024  label mesh0Facei = tree.shapes().faceLabels()[nearInfo.index()];
1025 
1026  // Check if points of f1 actually lie on top of mesh0 face
1027  // This is the bit that might fail since if f0 is severely warped
1028  // and f1's points are not present on f0 (since subdivision)
1029  // then the distance of the points to f0 might be quite large
1030  // and the test will fail. This all should in fact be some kind
1031  // of walk from known corresponding points/edges/faces.
1032  scalar dist =
1033  maxDistance
1034  (
1035  f1,
1036  mesh1.points(),
1037  mesh0.faces()[mesh0Facei],
1038  mesh0.points()
1039  );
1040 
1041  if (dist < absTol)
1042  {
1043  mesh0Set.insert(mesh0Facei);
1044  mesh1Set.insert(mesh1Facei);
1045  }
1046  }
1047  }
1048 
1049  if (debug)
1050  {
1051  Pout<< "findSlavesCoveringMaster :"
1052  << " matched " << mesh1Set.size() << " mesh1 faces to "
1053  << mesh0Set.size() << " mesh0 faces" << endl;
1054  }
1055 
1056  mesh0Faces = mesh0Set.toc();
1057  mesh1Faces = mesh1Set.toc();
1058 }
1059 
1060 
1061 Foam::label Foam::faceCoupleInfo::growCutFaces
1062 (
1063  const labelList& cutToMasterEdges,
1064  Map<labelList>& candidates
1065 )
1066 {
1067  if (debug)
1068  {
1069  Pout<< "growCutFaces :"
1070  << " growing cut faces to masterPatch" << endl;
1071  }
1072 
1073  label nTotChanged = 0;
1074 
1075  while (true)
1076  {
1077  const labelListList& cutFaceEdges = cutFaces().faceEdges();
1078  const labelListList& cutEdgeFaces = cutFaces().edgeFaces();
1079 
1080  label nChanged = 0;
1081 
1082  forAll(cutToMasterFaces_, cutFacei)
1083  {
1084  const label masterFacei = cutToMasterFaces_[cutFacei];
1085 
1086  if (masterFacei != -1)
1087  {
1088  // We now have a cutFace for which we have already found a
1089  // master face. Grow this masterface across any internal edge
1090  // (internal: no corresponding master edge)
1091 
1092  const labelList& fEdges = cutFaceEdges[cutFacei];
1093 
1094  forAll(fEdges, i)
1095  {
1096  const label cutEdgeI = fEdges[i];
1097 
1098  if (cutToMasterEdges[cutEdgeI] == -1)
1099  {
1100  // So this edge:
1101  // - internal to the cutPatch (since all region edges
1102  // marked before)
1103  // - on cutFacei which corresponds to masterFace.
1104  // Mark all connected faces with this masterFace.
1105 
1106  const labelList& eFaces = cutEdgeFaces[cutEdgeI];
1107 
1108  forAll(eFaces, j)
1109  {
1110  const label facei = eFaces[j];
1111 
1112  if (cutToMasterFaces_[facei] == -1)
1113  {
1114  cutToMasterFaces_[facei] = masterFacei;
1115  candidates.erase(facei);
1116  nChanged++;
1117  }
1118  else if (cutToMasterFaces_[facei] != masterFacei)
1119  {
1120  const pointField& cutPoints =
1121  cutFaces().points();
1122  const pointField& masterPoints =
1123  masterPatch().points();
1124 
1125  const edge& e = cutFaces().edges()[cutEdgeI];
1126 
1127  label myMaster = cutToMasterFaces_[facei];
1128  const face& myF = masterPatch()[myMaster];
1129 
1130  const face& nbrF = masterPatch()[masterFacei];
1131 
1133  << "Cut face "
1134  << cutFaces()[facei].points(cutPoints)
1135  << " has master "
1136  << myMaster
1137  << " but also connects to nbr face "
1138  << cutFaces()[cutFacei].points(cutPoints)
1139  << " with master " << masterFacei
1140  << nl
1141  << "myMasterFace:"
1142  << myF.points(masterPoints)
1143  << " nbrMasterFace:"
1144  << nbrF.points(masterPoints) << nl
1145  << "Across cut edge " << cutEdgeI
1146  << " coord:"
1147  << cutFaces().localPoints()[e[0]]
1148  << cutFaces().localPoints()[e[1]]
1149  << abort(FatalError);
1150  }
1151  }
1152  }
1153  }
1154  }
1155  }
1156 
1157  if (debug)
1158  {
1159  Pout<< "growCutFaces : Grown an additional " << nChanged
1160  << " cut to master face" << " correspondence" << endl;
1161  }
1162 
1163  nTotChanged += nChanged;
1164 
1165  if (nChanged == 0)
1166  {
1167  break;
1168  }
1169  }
1170 
1171  return nTotChanged;
1172 }
1173 
1174 
1175 void Foam::faceCoupleInfo::checkMatch(const labelList& cutToMasterEdges) const
1176 {
1177  const pointField& cutLocalPoints = cutFaces().localPoints();
1178 
1179  const pointField& masterLocalPoints = masterPatch().localPoints();
1180  const faceList& masterLocalFaces = masterPatch().localFaces();
1181 
1182  forAll(cutToMasterEdges, cutEdgeI)
1183  {
1184  const edge& e = cutFaces().edges()[cutEdgeI];
1185 
1186  if (cutToMasterEdges[cutEdgeI] == -1)
1187  {
1188  // Internal edge. Check that master face is same on both sides.
1189  const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1190 
1191  label masterFacei = -1;
1192 
1193  forAll(cutEFaces, i)
1194  {
1195  label cutFacei = cutEFaces[i];
1196 
1197  if (cutToMasterFaces_[cutFacei] != -1)
1198  {
1199  if (masterFacei == -1)
1200  {
1201  masterFacei = cutToMasterFaces_[cutFacei];
1202  }
1203  else if (masterFacei != cutToMasterFaces_[cutFacei])
1204  {
1205  label myMaster = cutToMasterFaces_[cutFacei];
1206  const face& myF = masterLocalFaces[myMaster];
1207 
1208  const face& nbrF = masterLocalFaces[masterFacei];
1209 
1211  << "Internal CutEdge " << e
1212  << " coord:"
1213  << cutLocalPoints[e[0]]
1214  << cutLocalPoints[e[1]]
1215  << " connects to master " << myMaster
1216  << " and to master " << masterFacei << nl
1217  << "myMasterFace:"
1218  << myF.points(masterLocalPoints)
1219  << " nbrMasterFace:"
1220  << nbrF.points(masterLocalPoints)
1221  << abort(FatalError);
1222  }
1223  }
1224  }
1225  }
1226  }
1227 }
1228 
1229 
1230 Foam::label Foam::faceCoupleInfo::matchEdgeFaces
1231 (
1232  const labelList& cutToMasterEdges,
1233  Map<labelList>& candidates
1234 )
1235 {
1236  // Extends matching information by elimination across cutFaces using more
1237  // than one region edge. Updates cutToMasterFaces_ and sets candidates which
1238  // is for every cutface on a region edge the possible master faces.
1239 
1240  // For every unassigned cutFacei the possible list of master faces.
1241  candidates.clear();
1242  candidates.resize(cutFaces().size());
1243 
1244  label nChanged = 0;
1245 
1246  forAll(cutToMasterEdges, cutEdgeI)
1247  {
1248  label masterEdgeI = cutToMasterEdges[cutEdgeI];
1249 
1250  if (masterEdgeI != -1)
1251  {
1252  // So cutEdgeI is matched to masterEdgeI. For all cut faces using
1253  // the cut edge store the master face as a candidate match.
1254 
1255  const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1256  const labelList& masterEFaces =
1257  masterPatch().edgeFaces()[masterEdgeI];
1258 
1259  forAll(cutEFaces, i)
1260  {
1261  label cutFacei = cutEFaces[i];
1262 
1263  if (cutToMasterFaces_[cutFacei] == -1)
1264  {
1265  // So this face (cutFacei) is on an edge (cutEdgeI) that
1266  // is used by some master faces (masterEFaces). Check if
1267  // the cutFacei and some of these masterEFaces share more
1268  // than one edge (which uniquely defines face)
1269 
1270  // Combine master faces with current set of candidate
1271  // master faces.
1272  Map<labelList>::iterator fnd = candidates.find(cutFacei);
1273 
1274  if (fnd == candidates.end())
1275  {
1276  // No info yet for cutFacei. Add all master faces as
1277  // candidates
1278  candidates.insert(cutFacei, masterEFaces);
1279  }
1280  else
1281  {
1282  // From some other cutEdgeI there are already some
1283  // candidate master faces. Check the overlap with
1284  // the current set of master faces.
1285  const labelList& masterFaces = fnd();
1286 
1287  DynamicList<label> newCandidates(masterFaces.size());
1288 
1289  forAll(masterEFaces, j)
1290  {
1291  if (findIndex(masterFaces, masterEFaces[j]) != -1)
1292  {
1293  newCandidates.append(masterEFaces[j]);
1294  }
1295  }
1296 
1297  if (newCandidates.size() == 1)
1298  {
1299  // We found a perfect match. Delete entry from
1300  // candidates map.
1301  cutToMasterFaces_[cutFacei] = newCandidates[0];
1302  candidates.erase(cutFacei);
1303  nChanged++;
1304  }
1305  else
1306  {
1307  // Should not happen?
1308  // Pout<< "On edge:" << cutEdgeI
1309  // << " have connected masterFaces:"
1310  // << masterEFaces
1311  // << " and from previous edge we have"
1312  // << " connected masterFaces" << masterFaces
1313  // << " . Overlap:" << newCandidates << endl;
1314 
1315  fnd() = newCandidates.shrink();
1316  }
1317  }
1318  }
1319 
1320  }
1321  }
1322  }
1323 
1324  if (debug)
1325  {
1326  Pout<< "matchEdgeFaces : Found " << nChanged
1327  << " faces where there was"
1328  << " only one remaining choice for cut-master correspondence"
1329  << endl;
1330  }
1331 
1332  return nChanged;
1333 }
1334 
1335 
1336 Foam::label Foam::faceCoupleInfo::geometricMatchEdgeFaces
1337 (
1338  Map<labelList>& candidates
1339 )
1340 {
1341  const pointField& cutPoints = cutFaces().points();
1342 
1343  label nChanged = 0;
1344 
1345  // Mark all master faces that have been matched so far.
1346 
1347  labelListList masterToCutFaces
1348  (
1350  (
1351  masterPatch().size(),
1352  cutToMasterFaces_
1353  )
1354  );
1355 
1356  forAllConstIter(Map<labelList>, candidates, iter)
1357  {
1358  label cutFacei = iter.key();
1359 
1360  const face& cutF = cutFaces()[cutFacei];
1361 
1362  if (cutToMasterFaces_[cutFacei] == -1)
1363  {
1364  const labelList& masterFaces = iter();
1365 
1366  // Find the best matching master face.
1367  scalar minDist = great;
1368  label minMasterFacei = -1;
1369 
1370  forAll(masterFaces, i)
1371  {
1372  label masterFacei = masterFaces[i];
1373 
1374  if (masterToCutFaces[masterFaces[i]].empty())
1375  {
1376  scalar dist = maxDistance
1377  (
1378  cutF,
1379  cutPoints,
1380  masterPatch()[masterFacei],
1381  masterPatch().points()
1382  );
1383 
1384  if (dist < minDist)
1385  {
1386  minDist = dist;
1387  minMasterFacei = masterFacei;
1388  }
1389  }
1390  }
1391 
1392  if (minMasterFacei != -1)
1393  {
1394  cutToMasterFaces_[cutFacei] = minMasterFacei;
1395  masterToCutFaces[minMasterFacei] = cutFacei;
1396  nChanged++;
1397  }
1398  }
1399  }
1400 
1401  // (inefficiently) force candidates to be uptodate.
1402  forAll(cutToMasterFaces_, cutFacei)
1403  {
1404  if (cutToMasterFaces_[cutFacei] != -1)
1405  {
1406  candidates.erase(cutFacei);
1407  }
1408  }
1409 
1410 
1411  if (debug)
1412  {
1413  Pout<< "geometricMatchEdgeFaces : Found " << nChanged
1414  << " faces where there was"
1415  << " only one remaining choice for cut-master correspondence"
1416  << endl;
1417  }
1418 
1419  return nChanged;
1420 }
1421 
1422 
1423 void Foam::faceCoupleInfo::perfectPointMatch
1424 (
1425  const scalar absTol,
1426  const bool slaveFacesOrdered
1427 )
1428 {
1429  // Calculate the set of cut faces in between master and slave patch assuming
1430  // perfect match (and optional face ordering on slave)
1431 
1432  if (debug)
1433  {
1434  Pout<< "perfectPointMatch :"
1435  << " Matching master and slave to cut."
1436  << " Master and slave faces are identical" << nl;
1437 
1438  if (slaveFacesOrdered)
1439  {
1440  Pout<< "and master and slave faces are ordered"
1441  << " (on coupled patches)" << endl;
1442  }
1443  else
1444  {
1445  Pout<< "and master and slave faces are not ordered"
1446  << " (on coupled patches)" << endl;
1447  }
1448  }
1449 
1450  cutToMasterFaces_ = identity(masterPatch().size());
1451  cutPoints_ = masterPatch().localPoints();
1452  cutFacesPtr_.reset
1453  (
1454  new primitiveFacePatch
1455  (
1456  masterPatch().localFaces(),
1457  cutPoints_
1458  )
1459  );
1460  masterToCutPoints_ = identity(cutPoints_.size());
1461 
1462 
1463  // Cut faces to slave patch.
1464  bool matchedAllFaces = false;
1465 
1466  if (slaveFacesOrdered)
1467  {
1468  cutToSlaveFaces_ = identity(cutFaces().size());
1469  matchedAllFaces = (cutFaces().size() == slavePatch().size());
1470  }
1471  else
1472  {
1473  // Faces do not have to be ordered (but all have
1474  // to match). Note: Faces will be already ordered if we enter here from
1475  // construct from meshes.
1476 
1477  matchedAllFaces = matchPoints
1478  (
1479  calcFaceCentres<List>
1480  (
1481  cutFaces(),
1482  cutFaces().points(),
1483  0,
1484  cutFaces().size()
1485  ),
1486  calcFaceCentres<IndirectList>
1487  (
1488  slavePatch(),
1489  slavePatch().points(),
1490  0,
1491  slavePatch().size()
1492  ),
1493  scalarField(slavePatch().size(), absTol),
1494  false,
1495  cutToSlaveFaces_
1496  );
1497 
1498  // If some of the face centres did not match, then try to match the
1499  // point averages instead. There is no division by the face area in
1500  // calculating the point average, so this is more stable when faces
1501  // collapse onto a line or point.
1502  if (!matchedAllFaces)
1503  {
1504  labelList cutToSlaveFacesTemp(cutToSlaveFaces_.size(), -1);
1505 
1506  matchPoints
1507  (
1508  calcFacePointAverages<List>
1509  (
1510  cutFaces(),
1511  cutFaces().points(),
1512  0,
1513  cutFaces().size()
1514  ),
1515  calcFacePointAverages<IndirectList>
1516  (
1517  slavePatch(),
1518  slavePatch().points(),
1519  0,
1520  slavePatch().size()
1521  ),
1522  scalarField(slavePatch().size(), absTol),
1523  true,
1524  cutToSlaveFacesTemp
1525  );
1526 
1527  cutToSlaveFaces_ = max(cutToSlaveFaces_, cutToSlaveFacesTemp);
1528 
1529  matchedAllFaces = min(cutToSlaveFaces_) != -1;
1530  }
1531  }
1532 
1533 
1534  if (!matchedAllFaces)
1535  {
1537  << "Did not match all of the master faces to the slave faces"
1538  << endl
1539  << "This usually means that the slave patch and master patch"
1540  << " do not align to within " << absTol << " metre."
1541  << abort(FatalError);
1542  }
1543 
1544 
1545  // Find correspondence from slave points to cut points. This might
1546  // detect shared points so the output is a slave-to-cut point list
1547  // and a compaction list.
1548 
1549  labelList cutToCompact, compactToCut;
1550  matchPointsThroughFaces
1551  (
1552  absTol,
1553  cutFaces().localPoints(),
1554  reorder(cutToSlaveFaces_, cutFaces().localFaces()),
1555  slavePatch().localPoints(),
1556  slavePatch().localFaces(),
1557  false, // slave and cut have opposite orientation
1558 
1559  slaveToCutPoints_, // slave to (uncompacted) cut points
1560  cutToCompact, // compaction map: from cut to compacted
1561  compactToCut // compaction map: from compacted to cut
1562  );
1563 
1564 
1565  // Use compaction lists to renumber cutPoints.
1566  cutPoints_ = UIndirectList<point>(cutPoints_, compactToCut)();
1567  {
1568  const faceList& cutLocalFaces = cutFaces().localFaces();
1569 
1570  faceList compactFaces(cutLocalFaces.size());
1571  forAll(cutLocalFaces, i)
1572  {
1573  compactFaces[i] = renumber(cutToCompact, cutLocalFaces[i]);
1574  }
1575  cutFacesPtr_.reset
1576  (
1577  new primitiveFacePatch
1578  (
1579  compactFaces,
1580  cutPoints_
1581  )
1582  );
1583  }
1584  inplaceRenumber(cutToCompact, slaveToCutPoints_);
1585  inplaceRenumber(cutToCompact, masterToCutPoints_);
1586 }
1587 
1588 
1589 void Foam::faceCoupleInfo::subDivisionMatch
1590 (
1591  const polyMesh& slaveMesh,
1592  const bool patchDivision,
1593  const scalar absTol
1594 )
1595 {
1596  if (debug)
1597  {
1598  Pout<< "subDivisionMatch :"
1599  << " Matching master and slave to cut."
1600  << " Slave can be subdivision of master but all master edges"
1601  << " have to be covered by slave edges." << endl;
1602  }
1603 
1604  // Assume slave patch is subdivision of the master patch so cutFaces is a
1605  // copy of the slave (but reversed since orientation has to be like master).
1606  cutPoints_ = slavePatch().localPoints();
1607  {
1608  faceList cutFaces(slavePatch().size());
1609 
1610  forAll(cutFaces, i)
1611  {
1612  cutFaces[i] = slavePatch().localFaces()[i].reverseFace();
1613  }
1614  cutFacesPtr_.reset(new primitiveFacePatch(cutFaces, cutPoints_));
1615  }
1616 
1617  // Cut is copy of slave so addressing to slave is trivial.
1618  cutToSlaveFaces_ = identity(cutFaces().size());
1619  slaveToCutPoints_ = identity(slavePatch().nPoints());
1620 
1621  // Cut edges to slave patch
1622  labelList cutToSlaveEdges
1623  (
1624  findMappedEdges
1625  (
1626  cutFaces().edges(),
1627  slaveToCutPoints_,
1628  slavePatch()
1629  )
1630  );
1631 
1632 
1633  if (debug)
1634  {
1635  OFstream str("regionEdges.obj");
1636 
1637  Pout<< "subDivisionMatch :"
1638  << " addressing for slave patch fully done."
1639  << " Dumping region edges to " << str.name() << endl;
1640 
1641  label vertI = 0;
1642 
1643  forAll(slavePatch().edges(), slaveEdgeI)
1644  {
1645  if (regionEdge(slaveMesh, slaveEdgeI))
1646  {
1647  const edge& e = slavePatch().edges()[slaveEdgeI];
1648 
1649  meshTools::writeOBJ(str, slavePatch().localPoints()[e[0]]);
1650  vertI++;
1651  meshTools::writeOBJ(str, slavePatch().localPoints()[e[1]]);
1652  vertI++;
1653  str<< "l " << vertI-1 << ' ' << vertI << nl;
1654  }
1655  }
1656  }
1657 
1658 
1659  // Addressing from cut to master.
1660 
1661  // Cut points to master patch
1662  // - determine master points to cut points. Has to be full!
1663  // - invert to get cut to master
1664  if (debug)
1665  {
1666  Pout<< "subDivisionMatch :"
1667  << " matching master points to cut points" << endl;
1668  }
1669 
1670  bool matchedAllPoints = matchPoints
1671  (
1672  masterPatch().localPoints(),
1673  cutPoints_,
1674  scalarField(masterPatch().nPoints(), absTol),
1675  true,
1676  masterToCutPoints_
1677  );
1678 
1679  if (!matchedAllPoints)
1680  {
1682  << "Did not match all of the master points to the slave points"
1683  << endl
1684  << "This usually means that the slave patch is not a"
1685  << " subdivision of the master patch"
1686  << abort(FatalError);
1687  }
1688 
1689 
1690  // Do masterEdges to cutEdges :
1691  // - mark all edges between two masterEdge endpoints. (geometric test since
1692  // this is the only distinction)
1693  // - this gives the borders in between which all cutfaces come from
1694  // a single master face.
1695  if (debug)
1696  {
1697  Pout<< "subDivisionMatch :"
1698  << " matching cut edges to master edges" << endl;
1699  }
1700 
1701  const edgeList& masterEdges = masterPatch().edges();
1702  const pointField& masterPoints = masterPatch().localPoints();
1703 
1704  const edgeList& cutEdges = cutFaces().edges();
1705 
1706  labelList cutToMasterEdges(cutFaces().nEdges(), -1);
1707 
1708  forAll(masterEdges, masterEdgeI)
1709  {
1710  const edge& masterEdge = masterEdges[masterEdgeI];
1711 
1712  label cutPoint0 = masterToCutPoints_[masterEdge[0]];
1713  label cutPoint1 = masterToCutPoints_[masterEdge[1]];
1714 
1715  // Find edges between cutPoint0 and cutPoint1.
1716 
1717  label cutPointi = cutPoint0;
1718  do
1719  {
1720  // Find edge (starting at pointi on cut), aligned with master
1721  // edge.
1722  label cutEdgeI =
1723  mostAlignedCutEdge
1724  (
1725  false,
1726  slaveMesh,
1727  patchDivision,
1728  cutToMasterEdges,
1729  cutToSlaveEdges,
1730  cutPointi,
1731  cutPoint0,
1732  cutPoint1
1733  );
1734 
1735  if (cutEdgeI == -1)
1736  {
1737  // Problem. Report while matching to produce nice error message
1738  mostAlignedCutEdge
1739  (
1740  true,
1741  slaveMesh,
1742  patchDivision,
1743  cutToMasterEdges,
1744  cutToSlaveEdges,
1745  cutPointi,
1746  cutPoint0,
1747  cutPoint1
1748  );
1749 
1750  Pout<< "Dumping unmatched pointEdges to errorEdges.obj"
1751  << endl;
1752 
1753  writeOBJ
1754  (
1755  "errorEdges.obj",
1756  edgeList
1757  (
1758  UIndirectList<edge>
1759  (
1760  cutFaces().edges(),
1761  cutFaces().pointEdges()[cutPointi]
1762  )
1763  ),
1764  cutFaces().localPoints(),
1765  false
1766  );
1767 
1769  << "Problem in finding cut edges corresponding to"
1770  << " master edge " << masterEdge
1771  << " points:" << masterPoints[masterEdge[0]]
1772  << ' ' << masterPoints[masterEdge[1]]
1773  << " corresponding cut edge: (" << cutPoint0
1774  << ") " << cutPoint1
1775  << abort(FatalError);
1776  }
1777 
1778  cutToMasterEdges[cutEdgeI] = masterEdgeI;
1779 
1780  cutPointi = cutEdges[cutEdgeI].otherVertex(cutPointi);
1781 
1782  } while (cutPointi != cutPoint1);
1783  }
1784 
1785  if (debug)
1786  {
1787  // Write all matched edges.
1788  writeEdges(cutToMasterEdges, cutToSlaveEdges);
1789  }
1790 
1791  // Rework cutToMasterEdges into list of points in between two endpoints
1792  // (cutEdgeToPoints_)
1793  setCutEdgeToPoints(cutToMasterEdges);
1794 
1795 
1796  // Now that we have marked the cut edges that lie on top of master edges
1797  // we can find cut faces that have two (or more) of these edges.
1798  // Note that there can be multiple faces having two or more matched edges
1799  // since the cut faces can form a non-manifold surface(?)
1800  // So the matching is done as an elimination process where for every
1801  // cutFace on a matched edge we store the possible master faces and
1802  // eliminate more and more until we only have one possible master face
1803  // left.
1804 
1805  if (debug)
1806  {
1807  Pout<< "subDivisionMatch :"
1808  << " matching (topological) cut faces to masterPatch" << endl;
1809  }
1810 
1811  // For every unassigned cutFacei the possible list of master faces.
1812  Map<labelList> candidates(cutFaces().size());
1813 
1814  cutToMasterFaces_.setSize(cutFaces().size());
1815  cutToMasterFaces_ = -1;
1816 
1817  while (true)
1818  {
1819  // See if there are any edges left where there is only one remaining
1820  // candidate.
1821  label nChanged = matchEdgeFaces(cutToMasterEdges, candidates);
1822 
1823  checkMatch(cutToMasterEdges);
1824 
1825 
1826  // Now we should have found a face correspondence for every cutFace
1827  // that uses two or more cutEdges. Floodfill from these across
1828  // 'internal' cutedges (i.e. edges that do not have a master
1829  // equivalent) to determine some more cutToMasterFaces
1830  nChanged += growCutFaces(cutToMasterEdges, candidates);
1831 
1832  checkMatch(cutToMasterEdges);
1833 
1834  if (nChanged == 0)
1835  {
1836  break;
1837  }
1838  }
1839 
1840  if (debug)
1841  {
1842  Pout<< "subDivisionMatch :"
1843  << " matching (geometric) cut faces to masterPatch" << endl;
1844  }
1845 
1846  // Above should have matched all cases fully. Cannot prove this yet so
1847  // do any remaining unmatched edges by a geometric test.
1848  while (true)
1849  {
1850  // See if there are any edges left where there is only one remaining
1851  // candidate.
1852  label nChanged = geometricMatchEdgeFaces(candidates);
1853 
1854  checkMatch(cutToMasterEdges);
1855 
1856  nChanged += growCutFaces(cutToMasterEdges, candidates);
1857 
1858  checkMatch(cutToMasterEdges);
1859 
1860  if (nChanged == 0)
1861  {
1862  break;
1863  }
1864  }
1865 
1866 
1867  // All cut faces matched?
1868  forAll(cutToMasterFaces_, cutFacei)
1869  {
1870  if (cutToMasterFaces_[cutFacei] == -1)
1871  {
1872  const face& cutF = cutFaces()[cutFacei];
1873 
1875  << "Did not match all of cutFaces to a master face" << nl
1876  << "First unmatched cut face:" << cutFacei << " with points:"
1877  << UIndirectList<point>(cutFaces().points(), cutF)()
1878  << nl
1879  << "This usually means that the slave patch is not a"
1880  << " subdivision of the master patch"
1881  << abort(FatalError);
1882  }
1883  }
1884 
1885  if (debug)
1886  {
1887  Pout<< "subDivisionMatch :"
1888  << " finished matching master and slave to cut" << endl;
1889  }
1890 
1891  if (debug)
1892  {
1893  writePointsFaces();
1894  }
1895 }
1896 
1897 
1898 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1899 
1902  const polyMesh& masterMesh,
1903  const polyMesh& slaveMesh,
1904  const scalar absTol,
1905  const bool perfectMatch
1906 )
1907 :
1908  masterPatchPtr_(nullptr),
1909  slavePatchPtr_(nullptr),
1910  cutPoints_(0),
1911  cutFacesPtr_(nullptr),
1912  cutToMasterFaces_(0),
1913  masterToCutPoints_(0),
1914  cutToSlaveFaces_(0),
1915  slaveToCutPoints_(0),
1916  cutEdgeToPoints_(0)
1917 {
1918  // Get faces on both meshes that are aligned.
1919  // (not ordered i.e. masterToMesh[0] does
1920  // not couple to slaveToMesh[0]. This ordering is done later on)
1921  labelList masterToMesh;
1922  labelList slaveToMesh;
1923 
1924  if (perfectMatch)
1925  {
1926  // Identical faces so use tight face-centre comparison.
1927  findPerfectMatchingFaces
1928  (
1929  masterMesh,
1930  slaveMesh,
1931  absTol,
1932  masterToMesh,
1933  slaveToMesh
1934  );
1935  }
1936  else
1937  {
1938  // Slave subdivision of master so use 'nearest'. Bit dodgy, especially
1939  // with using absTol (which is quite small)
1940  findSlavesCoveringMaster
1941  (
1942  masterMesh,
1943  slaveMesh,
1944  absTol,
1945  masterToMesh,
1946  slaveToMesh
1947  );
1948  }
1949 
1950  // Construct addressing engines for both sides
1951  masterPatchPtr_.reset
1952  (
1954  (
1955  IndirectList<face>(masterMesh.faces(), masterToMesh),
1956  masterMesh.points()
1957  )
1958  );
1959 
1960  slavePatchPtr_.reset
1961  (
1963  (
1964  IndirectList<face>(slaveMesh.faces(), slaveToMesh),
1965  slaveMesh.points()
1966  )
1967  );
1968 
1969 
1970  if (perfectMatch)
1971  {
1972  // Faces are perfectly aligned but probably not ordered.
1973  perfectPointMatch(absTol, false);
1974  }
1975  else
1976  {
1977  // Slave faces are subdivision of master face. Faces not ordered.
1978  subDivisionMatch(slaveMesh, false, absTol);
1979  }
1980 
1981  if (debug)
1982  {
1983  writePointsFaces();
1984  }
1985 }
1986 
1987 
1990  const polyMesh& masterMesh,
1991  const labelList& masterAddressing,
1992  const polyMesh& slaveMesh,
1993  const labelList& slaveAddressing,
1994  const scalar absTol,
1995  const bool perfectMatch,
1996  const bool orderedFaces,
1997  const bool patchDivision
1998 )
1999 :
2000  masterPatchPtr_
2001  (
2003  (
2004  IndirectList<face>(masterMesh.faces(), masterAddressing),
2005  masterMesh.points()
2006  )
2007  ),
2008  slavePatchPtr_
2009  (
2011  (
2012  IndirectList<face>(slaveMesh.faces(), slaveAddressing),
2013  slaveMesh.points()
2014  )
2015  ),
2016  cutPoints_(0),
2017  cutFacesPtr_(nullptr),
2018  cutToMasterFaces_(0),
2019  masterToCutPoints_(0),
2020  cutToSlaveFaces_(0),
2021  slaveToCutPoints_(0),
2022  cutEdgeToPoints_(0)
2023 {
2024  if (perfectMatch && (masterAddressing.size() != slaveAddressing.size()))
2025  {
2027  << "Perfect match specified but number of master and slave faces"
2028  << " differ." << endl
2029  << "master:" << masterAddressing.size()
2030  << " slave:" << slaveAddressing.size()
2031  << abort(FatalError);
2032  }
2033 
2034  if
2035  (
2036  masterAddressing.size()
2037  && min(masterAddressing) < masterMesh.nInternalFaces()
2038  )
2039  {
2041  << "Supplied internal face on master mesh to couple." << nl
2042  << "Faces to be coupled have to be boundary faces."
2043  << abort(FatalError);
2044  }
2045  if
2046  (
2047  slaveAddressing.size()
2048  && min(slaveAddressing) < slaveMesh.nInternalFaces()
2049  )
2050  {
2052  << "Supplied internal face on slave mesh to couple." << nl
2053  << "Faces to be coupled have to be boundary faces."
2054  << abort(FatalError);
2055  }
2056 
2057 
2058  if (perfectMatch)
2059  {
2060  perfectPointMatch(absTol, orderedFaces);
2061  }
2062  else
2063  {
2064  // Slave faces are subdivision of master face. Faces not ordered.
2065  subDivisionMatch(slaveMesh, patchDivision, absTol);
2066  }
2067 
2068  if (debug)
2069  {
2070  writePointsFaces();
2071  }
2072 }
2073 
2074 
2075 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
2076 
2078 {}
2079 
2080 
2081 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2082 
2084 {
2085  labelList faces(pp.size());
2086 
2087  label facei = pp.start();
2088 
2089  forAll(pp, i)
2090  {
2091  faces[i] = facei++;
2092  }
2093  return faces;
2094 }
2095 
2096 
2098 {
2099  Map<label> map(lst.size());
2100 
2101  forAll(lst, i)
2102  {
2103  if (lst[i] != -1)
2104  {
2105  map.insert(i, lst[i]);
2106  }
2107  }
2108  return map;
2109 }
2110 
2111 
2114  const labelListList& lst
2115 )
2116 {
2117  Map<labelList> map(lst.size());
2118 
2119  forAll(lst, i)
2120  {
2121  if (lst[i].size())
2122  {
2123  map.insert(i, lst[i]);
2124  }
2125  }
2126  return map;
2127 }
2128 
2129 
2130 // ************************************************************************* //
ListType renumber(const labelUList &oldToNew, const ListType &)
Renumber the values (not the indices) of a list.
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:82
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
label nInternalFaces() const
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:164
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
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:138
PrimitivePatch< IndirectList< face >, 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:1131
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
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
const pointField & points
List< edge > edgeList
Definition: edgeList.H:38
label nPoints
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
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 > &)
~faceCoupleInfo()
Destructor.
static const char nl
Definition: Ostream.H:260
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 occurrence of given element and return index,.
PrimitivePatch< List< face >, 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.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:303
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
A List with indirect addressing.
Definition: IndirectList.H:101
Namespace for OpenFOAM.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49