coupleSlidingInterface.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 "slidingInterface.H"
27 #include "polyTopoChange.H"
28 #include "polyMesh.H"
29 #include "primitiveMesh.H"
30 #include "enrichedPatch.H"
31 #include "DynamicList.H"
32 #include "pointHit.H"
33 #include "triPointRef.H"
34 #include "plane.H"
35 #include "polyTopoChanger.H"
36 #include "polyAddPoint.H"
37 #include "polyRemovePoint.H"
38 #include "polyAddFace.H"
39 #include "polyModifyPoint.H"
40 #include "polyModifyFace.H"
41 #include "polyRemoveFace.H"
42 
43 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44 
45 const Foam::scalar Foam::slidingInterface::edgeCoPlanarTolDefault_ = 0.8;
46 
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
49 // Index of debug signs:
50 // Index of debug signs:
51 // . - loop of the edge-to-face interaction detection
52 // x - reversal of direction in edge-to-face interaction detection
53 // + - complete edge-to-face interaction detection
54 // z - incomplete edge-to-face interaction detection. This may be OK for edges
55 // crossing from one to the other side of multiply connected master patch
56 
57 // e - adding a point on edge
58 // f - adding a point on face
59 // . - collecting edges off another face for edge-to-edge cut
60 // + - finished collection of edges
61 // * - cut both master and slave
62 // n - cutting new edge
63 // t - intersection exists but it is outside of tolerance
64 // x - missed slave edge in cut
65 // - - missed master edge in cut
66 // u - edge already used in cutting
67 
68 void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
69 {
70  if (debug)
71  {
72  Pout<< "void slidingInterface::coupleInterface"
73  << "(polyTopoChange& ref) : "
74  << "Coupling sliding interface " << name() << endl;
75  }
76 
77  const polyMesh& mesh = topoChanger().mesh();
78 
79  const pointField& points = mesh.points();
80  const faceList& faces = mesh.faces();
81 
82  const labelList& own = mesh.faceOwner();
83  const labelList& nei = mesh.faceNeighbour();
84  const faceZoneMesh& faceZones = mesh.faceZones();
85 
86  const primitiveFacePatch& masterPatch =
87  faceZones[masterFaceZoneID_.index()]();
88 
89  const labelList& masterPatchAddr = faceZones[masterFaceZoneID_.index()];
90 
91  const boolList& masterPatchFlip =
92  faceZones[masterFaceZoneID_.index()].flipMap();
93 
94  const primitiveFacePatch& slavePatch =
95  faceZones[slaveFaceZoneID_.index()]();
96 
97  const labelList& slavePatchAddr = faceZones[slaveFaceZoneID_.index()];
98 
99  const boolList& slavePatchFlip =
100  faceZones[slaveFaceZoneID_.index()].flipMap();
101 
102  const edgeList& masterEdges = masterPatch.edges();
103  const labelListList& masterPointEdges = masterPatch.pointEdges();
104  const labelList& masterMeshPoints = masterPatch.meshPoints();
105  const pointField& masterLocalPoints = masterPatch.localPoints();
106  const labelListList& masterFaceFaces = masterPatch.faceFaces();
107  const labelListList& masterFaceEdges = masterPatch.faceEdges();
108  const Map<label>& masterMeshPointMap = masterPatch.meshPointMap();
109 
110  const edgeList& slaveEdges = slavePatch.edges();
111  const labelListList& slavePointEdges = slavePatch.pointEdges();
112  const labelList& slaveMeshPoints = slavePatch.meshPoints();
113  const pointField& slaveLocalPoints = slavePatch.localPoints();
114  const Map<label>& slaveMeshPointMap = slavePatch.meshPointMap();
115  const vectorField& slavePointNormals = slavePatch.pointNormals();
116 
117  // Collect projection addressing
118  if
119  (
120  !(
121  slavePointPointHitsPtr_
122  && slavePointEdgeHitsPtr_
123  && slavePointFaceHitsPtr_
124  && masterPointEdgeHitsPtr_
125  )
126  )
127  {
129  << "Point projection addressing not available."
130  << abort(FatalError);
131  }
132 
133  const labelList& slavePointPointHits = *slavePointPointHitsPtr_;
134  const labelList& slavePointEdgeHits = *slavePointEdgeHitsPtr_;
135  const List<objectHit>& slavePointFaceHits = *slavePointFaceHitsPtr_;
136  const labelList& masterPointEdgeHits = *masterPointEdgeHitsPtr_;
137  const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
138 
139  // Create enriched patch
140  enrichedPatch cutPatch
141  (
142  masterPatch,
143  slavePatch,
144  slavePointPointHits,
145  slavePointEdgeHits,
146  slavePointFaceHits
147  );
148 
149  // Get reference to list of added point for the enriched patch.
150  // This will be used for point addition
151  Map<point>& pointMap = cutPatch.pointMap();
152 
153  // Get reference to the list of merged points
154  Map<label>& pointMergeMap = cutPatch.pointMergeMap();
155 
156  // Create mapping for every merged point of the slave patch
157  forAll(slavePointPointHits, pointi)
158  {
159  if (slavePointPointHits[pointi] >= 0)
160  {
161  // Pout<< "Inserting point merge pair: " << slaveMeshPoints[pointi]
162  // << " : " << masterMeshPoints[slavePointPointHits[pointi]]
163  // << endl;
164 
165  pointMergeMap.insert
166  (
167  slaveMeshPoints[pointi],
168  masterMeshPoints[slavePointPointHits[pointi]]
169  );
170  }
171  }
172 
173  // Collect the list of used edges for every slave edge
174 
175  List<labelHashSet> usedMasterEdges(slaveEdges.size());
176 
177  // Collect of slave point hits
178  forAll(slavePointPointHits, pointi)
179  {
180  // For point hits, add all point-edges from master side to all point
181  const labelList& curSlaveEdges = slavePointEdges[pointi];
182 
183  if (slavePointPointHits[pointi] > -1)
184  {
185  const labelList& curMasterEdges =
186  masterPointEdges[slavePointPointHits[pointi]];
187 
188  // Mark all current master edges as used for all the current slave
189  // edges
190  forAll(curSlaveEdges, slaveEdgeI)
191  {
192  labelHashSet& sm = usedMasterEdges[curSlaveEdges[slaveEdgeI]];
193 
194  forAll(curMasterEdges, masterEdgeI)
195  {
196  sm.insert(curMasterEdges[masterEdgeI]);
197  }
198  }
199  }
200  else if (slavePointEdgeHits[pointi] > -1)
201  {
202  // For edge hits, add the master edge
203  forAll(curSlaveEdges, slaveEdgeI)
204  {
205  usedMasterEdges[curSlaveEdges[slaveEdgeI]].insert
206  (
207  slavePointEdgeHits[pointi]
208  );
209  }
210  }
211  }
212 
213  // Collect off master point hits
214  // For every master point that has hit an edge, add all edges coming from
215  // that point to the slave edge that has been hit
216  forAll(masterPointEdgeHits, masterPointi)
217  {
218  if (masterPointEdgeHits[masterPointi] > -1)
219  {
220  const labelList& curMasterEdges = masterPointEdges[masterPointi];
221 
222  labelHashSet& sm =
223  usedMasterEdges[masterPointEdgeHits[masterPointi]];
224 
225  forAll(curMasterEdges, masterEdgeI)
226  {
227  sm.insert(curMasterEdges[masterEdgeI]);
228  }
229  }
230  }
231 
232  // Pout<< "used edges: " << endl;
233  // forAll(usedMasterEdges, edgeI)
234  // {
235  // Pout<< "edge: " << edgeI
236  // << " used: " << usedMasterEdges[edgeI].toc()
237  // << endl;
238  // }
239 
240  // For every master and slave edge make a list of points to be added into
241  // that edge.
242  List<DynamicList<label>> pointsIntoMasterEdges(masterEdges.size());
243  List<DynamicList<label>> pointsIntoSlaveEdges(slaveEdges.size());
244 
245  // Add all points from the slave patch that have hit the edge
246  forAll(slavePointEdgeHits, pointi)
247  {
248  if (slavePointEdgeHits[pointi] > -1)
249  {
250  // Create a new point on the master edge
251 
252  point edgeCutPoint =
253  masterEdges[slavePointEdgeHits[pointi]].line
254  (
255  masterLocalPoints
256  ).nearestDist(projectedSlavePoints[pointi]).hitPoint();
257 
258  label newPoint =
259  ref.setAction
260  (
261  polyAddPoint
262  (
263  edgeCutPoint, // point
264  slaveMeshPoints[pointi], // master point
265  cutPointZoneID_.index(), // zone for point
266  true // supports a cell
267  )
268  );
269 
270  // Pout<< "Inserting merge pair off edge: "
271  // << slaveMeshPoints[pointi] << " " << newPoint
272  // << " cut point: " << edgeCutPoint
273  // << " orig: " << slaveLocalPoints[pointi]
274  // << " proj: " << projectedSlavePoints[pointi]
275  // << endl;
276 
277  // Add the new edge point into the merge map
278  pointMergeMap.insert(slaveMeshPoints[pointi], newPoint);
279 
280  pointsIntoMasterEdges[slavePointEdgeHits[pointi]].append
281  (
282  newPoint
283  );
284 
285  // Add the point into the enriched patch map
286  pointMap.insert
287  (
288  newPoint,
289  edgeCutPoint
290  );
291 
292  if (debug)
293  {
294  Pout<< "e";
295  // Pout<< newPoint << " = " << edgeCutPoint << endl;
296  }
297  }
298  }
299 
300  if (debug)
301  {
302  Pout<< endl;
303  }
304 
305  // Add all points from the slave patch that have hit a face
306  forAll(slavePointFaceHits, pointi)
307  {
308  if
309  (
310  slavePointPointHits[pointi] < 0
311  && slavePointEdgeHits[pointi] < 0
312  && slavePointFaceHits[pointi].hit()
313  )
314  {
315  label newPoint =
316  ref.setAction
317  (
318  polyAddPoint
319  (
320  projectedSlavePoints[pointi], // point
321  slaveMeshPoints[pointi], // master point
322  cutPointZoneID_.index(), // zone for point
323  true // supports a cell
324  )
325  );
326 
327  // Pout<< "Inserting merge pair off face: "
328  // << slaveMeshPoints[pointi]
329  // << " " << newPoint
330  // << endl;
331 
332  // Add the new edge point into the merge map
333  pointMergeMap.insert(slaveMeshPoints[pointi], newPoint);
334 
335  // Add the point into the enriched patch map
336  pointMap.insert
337  (
338  newPoint,
339  projectedSlavePoints[pointi]
340  );
341 
342  if (debug)
343  {
344  Pout<< "f: " << newPoint << " = "
345  << projectedSlavePoints[pointi] << endl;
346  }
347  }
348  }
349 
350  forAll(masterPointEdgeHits, pointi)
351  {
352  if (masterPointEdgeHits[pointi] > -1)
353  {
354  pointsIntoSlaveEdges[masterPointEdgeHits[pointi]].append
355  (
356  masterMeshPoints[pointi]
357  );
358  }
359  }
360 
361  // Cut all slave edges.
362  // Collect all master edges the slave edge interacts with. Skip
363  // all the edges already marked as used. For every unused edge,
364  // calculate the cut and insert the new point into the master and
365  // slave edge.
366  // For the edge selection algorithm, see, comment in
367  // slidingInterfaceProjectPoints.C.
368  // Edge cutting algorithm:
369  // As the master patch defines the cutting surface, the newly
370  // inserted point needs to be on the master edge. Also, in 3-D
371  // the pair of edges generally misses each other rather than
372  // intersect. Therefore, the intersection is calculated using the
373  // plane the slave edge defines during projection. The plane is
374  // defined by the centrepoint of the edge in the original
375  // configuration and the projected end points. In case of flat
376  // geometries (when the three points are colinear), the plane is
377  // defined by the two projected end-points and the slave edge
378  // normal used as the in-plane vector. When the intersection
379  // point is created, it is added as a new point for both the
380  // master and the slave edge.
381  // In order to be able to re-create the points on edges in mesh
382  // motion without the topology change, the edge pair used to
383  // create the cut point needs to be recorded in
384  // cutPointEdgePairMap
385 
386  // Note. "Processing slave edges" code is repeated twice in the
387  // sliding intergace coupling in order to allow the point
388  // projection to be done separately from the actual cutting.
389  // Please change consistently with slidingInterfaceProjectPoints.C
390  //
391  if (debug)
392  {
393  Pout<< "Processing slave edges " << endl;
394  }
395 
396  if (!cutPointEdgePairMapPtr_)
397  {
399  << "Cut point edge pair map pointer not set."
400  << abort(FatalError);
401  }
402 
403  Map<Pair<edge>>& addToCpepm = *cutPointEdgePairMapPtr_;
404 
405  // Clear the old map
406  addToCpepm.clear();
407 
408  // Create a map of faces the edge can interact with
409  labelHashSet curFaceMap
410  (
411  nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
412  );
413 
415 
416  forAll(slaveEdges, edgeI)
417  {
418  const edge& curEdge = slaveEdges[edgeI];
419 
420  if
421  (
422  slavePointFaceHits[curEdge.start()].hit()
423  || slavePointFaceHits[curEdge.end()].hit()
424  )
425  {
426  labelHashSet& curUme = usedMasterEdges[edgeI];
427 
428  // Pout<< "Doing edge " << edgeI
429  // << " curEdge: " << curEdge
430  // << " curUme: " << curUme
431  // << endl;
432 
433  // Clear the maps
434  curFaceMap.clear();
435  addedFaces.clear();
436 
437  // Grab the faces for start and end points.
438  const label startFace =
439  slavePointFaceHits[curEdge.start()].hitObject();
440  const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
441 
442  // Pout<< "startFace: " << slavePointFaceHits[curEdge.start()]
443  // << " endFace: " << slavePointFaceHits[curEdge.end()]
444  // << endl;
445 
446  // Insert the start face into the list
447  curFaceMap.insert(startFace);
448  addedFaces.insert(startFace);
449 
450  // Pout<< "curFaceMap: " << curFaceMap.toc() << endl;
451 
452  label nSweeps = 0;
453  bool completed = false;
454 
455  while (nSweeps < edgeFaceEscapeLimit_)
456  {
457  nSweeps++;
458 
459  if (addedFaces.found(endFace))
460  {
461  completed = true;
462  }
463 
464  // Add all face neighbours of face in the map
465  const labelList cf = addedFaces.toc();
466  addedFaces.clear();
467 
468  forAll(cf, cfI)
469  {
470  const labelList& curNbrs = masterFaceFaces[cf[cfI]];
471 
472  forAll(curNbrs, nbrI)
473  {
474  if (!curFaceMap.found(curNbrs[nbrI]))
475  {
476  curFaceMap.insert(curNbrs[nbrI]);
477  addedFaces.insert(curNbrs[nbrI]);
478  }
479  }
480  }
481 
482  if (completed) break;
483 
484  if (debug)
485  {
486  Pout<< ".";
487  }
488  }
489 
490  if (!completed)
491  {
492  if (debug)
493  {
494  Pout<< "x";
495  }
496 
497  // It is impossible to reach the end from the start, probably
498  // due to disconnected domain. Do search in opposite direction
499 
500  label nReverseSweeps = 0;
501 
502  addedFaces.clear();
503  addedFaces.insert(endFace);
504 
505  while (nReverseSweeps < edgeFaceEscapeLimit_)
506  {
507  nReverseSweeps++;
508 
509  if (addedFaces.found(startFace))
510  {
511  completed = true;
512  }
513 
514  // Add all face neighbours of face in the map
515  const labelList cf = addedFaces.toc();
516  addedFaces.clear();
517 
518  forAll(cf, cfI)
519  {
520  const labelList& curNbrs = masterFaceFaces[cf[cfI]];
521 
522  forAll(curNbrs, nbrI)
523  {
524  if (!curFaceMap.found(curNbrs[nbrI]))
525  {
526  curFaceMap.insert(curNbrs[nbrI]);
527  addedFaces.insert(curNbrs[nbrI]);
528  }
529  }
530  }
531 
532  if (completed) break;
533 
534  if (debug)
535  {
536  Pout<< ".";
537  }
538  }
539  }
540 
541  if (completed)
542  {
543  if (debug)
544  {
545  Pout<< "+ ";
546  }
547  }
548  else
549  {
550  if (debug)
551  {
552  Pout<< "z ";
553  }
554  }
555 
556  // Collect the edges
557 
558  // Create a map of edges the edge can interact with
559  labelHashSet curMasterEdgesMap
560  (
561  nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
562  );
563 
564  const labelList curFaces = curFaceMap.toc();
565 
566  // Pout<< "curFaces: " << curFaces << endl;
567 
568  forAll(curFaces, facei)
569  {
570  // Pout<< "face: " << curFaces[facei] << " "
571  // << masterPatch[curFaces[facei]]
572  // << " local: "
573  // << masterPatch.localFaces()[curFaces[facei]]
574  // << endl;
575 
576  const labelList& me = masterFaceEdges[curFaces[facei]];
577 
578  forAll(me, meI)
579  {
580  curMasterEdgesMap.insert(me[meI]);
581  }
582  }
583 
584  const labelList curMasterEdges = curMasterEdgesMap.toc();
585 
586  // For all master edges to intersect, skip the ones
587  // already used and cut the rest with a cutting plane. If
588  // the intersection point, falls inside of both edges, it
589  // is valid.
590 
591  // Note.
592  // The edge cutting code is repeated in
593  // slidingInterface::modifyMotionPoints. This is done for
594  // efficiency reasons and avoids multiple creation of cutting
595  // planes. Please update both simultaneously.
596 
597  const point& a = projectedSlavePoints[curEdge.start()];
598  const point& b = projectedSlavePoints[curEdge.end()];
599 
600  point c =
601  0.5*
602  (
603  slaveLocalPoints[curEdge.start()]
604  + slavePointNormals[curEdge.start()] // Add start normal
605  + slaveLocalPoints[curEdge.end()]
606  + slavePointNormals[curEdge.end()] // Add end normal
607  );
608 
609  // Create the plane
610  plane cutPlane(a, b, c);
611 
612  // Pout<< "a: " << a
613  // << " b: " << b
614  // << " c: " << c
615  // << " plane: " << cutPlane
616  // << endl;
617 
618  linePointRef curSlaveLine = curEdge.line(projectedSlavePoints);
619  const scalar curSlaveLineMag = curSlaveLine.mag();
620 
621  // Pout<< "curSlaveLine: " << curSlaveLine << endl;
622 
623  forAll(curMasterEdges, masterEdgeI)
624  {
625  if (!curUme.found(curMasterEdges[masterEdgeI]))
626  {
627  // New edge
628  if (debug)
629  {
630  Pout<< "n";
631  }
632 
633  const label cmeIndex = curMasterEdges[masterEdgeI];
634  const edge& cme = masterEdges[cmeIndex];
635 
636  // Pout<< "Edge " << cmeIndex
637  // << " cme: " << cme
638  // << " line: " << cme.line(masterLocalPoints)
639  // << endl;
640 
641  scalar cutOnMaster =
642  cutPlane.lineIntersect
643  (
644  cme.line(masterLocalPoints)
645  );
646 
647  if
648  (
649  cutOnMaster > edgeEndCutoffTol_
650  && cutOnMaster < 1.0 - edgeEndCutoffTol_
651  )
652  {
653  // Master is cut, check the slave
654  point masterCutPoint =
655  masterLocalPoints[cme.start()]
656  + cutOnMaster*cme.vec(masterLocalPoints);
657 
658  pointHit slaveCut =
659  curSlaveLine.nearestDist(masterCutPoint);
660 
661  if (slaveCut.hit())
662  {
663  // Strict checking of slave cut to avoid capturing
664  // end points.
665  scalar cutOnSlave =
666  (
667  (
668  slaveCut.hitPoint()
669  - curSlaveLine.start()
670  ) & curSlaveLine.vec()
671  )/sqr(curSlaveLineMag);
672 
673  // Calculate merge tolerance from the
674  // target edge length
675  scalar mergeTol = edgeCoPlanarTol_*mag(b - a);
676 
677  // Pout<< "cutOnMaster: " << cutOnMaster
678  // << " masterCutPoint: " << masterCutPoint
679  // << " slaveCutPoint: " << slaveCut.hitPoint()
680  // << " slaveCut.distance(): "
681  // << slaveCut.distance()
682  // << " slave length: " << mag(b - a)
683  // << " mergeTol: " << mergeTol
684  // << " 1: " << mag(b - a)
685  // << " 2: "
686  // << cme.line(masterLocalPoints).mag()
687  // << endl;
688 
689  if
690  (
691  cutOnSlave > edgeEndCutoffTol_
692  && cutOnSlave < 1.0 - edgeEndCutoffTol_
693  && slaveCut.distance() < mergeTol
694  )
695  {
696  // Cut both master and slave. Add point
697  // to edge points The point is nominally
698  // added from the start of the master edge
699  // and added to the cut point zone
700  label newPoint =
701  ref.setAction
702  (
703  polyAddPoint
704  (
705  masterCutPoint, // point
706  masterMeshPoints[cme.start()],// m p
707  cutPointZoneID_.index(), // zone
708  true // active
709  )
710  );
711 
712  // Pout<< "Inserting point: " << newPoint
713  // << " as edge to edge intersection. "
714  // << "Slave edge: "
715  // << edgeI << " " << curEdge
716  // << " master edge: "
717  // << cmeIndex << " " << cme
718  // << endl;
719 
720  pointsIntoSlaveEdges[edgeI].append(newPoint);
721  pointsIntoMasterEdges[cmeIndex].append
722  (
723  newPoint
724  );
725 
726  // Add the point into the enriched patch map
727  pointMap.insert
728  (
729  newPoint,
730  masterCutPoint
731  );
732 
733  // Record which two edges intersect to
734  // create cut point
735  addToCpepm.insert
736  (
737  newPoint, // Cut point index
738  Pair<edge>
739  (
740  edge
741  (
742  masterMeshPoints[cme.start()],
743  masterMeshPoints[cme.end()]
744  ), // Master edge
745  edge
746  (
747  slaveMeshPoints[curEdge.start()],
748  slaveMeshPoints[curEdge.end()]
749  )// Slave edge
750  )
751  );
752 
753  if (debug)
754  {
755  Pout<< " " << newPoint << " = "
756  << masterCutPoint << " ";
757  }
758  }
759  else
760  {
761  if (debug)
762  {
763  // Intersection exists but it is too far
764  Pout<< "t";
765  }
766  }
767  }
768  else
769  {
770  if (debug)
771  {
772  // Missed slave edge
773  Pout<< "x";
774  }
775  }
776  }
777  else
778  {
779  if (debug)
780  {
781  // Missed master edge
782  Pout<< "-";
783  }
784  }
785  }
786  else
787  {
788  if (debug)
789  {
790  Pout<< "u";
791  }
792  }
793  }
794 
795  if (debug)
796  {
797  Pout<< endl;
798  }
799  } // End if both ends missing
800  } // End for all slave edges
801 
802 // Pout<< "pointsIntoMasterEdges: " << pointsIntoMasterEdges << endl;
803 // Pout<< "pointsIntoSlaveEdges: " << pointsIntoSlaveEdges << endl;
804 
805  // Re-pack the points into edges lists
806  labelListList pime(pointsIntoMasterEdges.size());
807 
808  forAll(pointsIntoMasterEdges, i)
809  {
810  pime[i].transfer(pointsIntoMasterEdges[i]);
811  }
812 
813  labelListList pise(pointsIntoSlaveEdges.size());
814 
815  forAll(pointsIntoSlaveEdges, i)
816  {
817  pise[i].transfer(pointsIntoSlaveEdges[i]);
818  }
819 
820  // Prepare the enriched faces
821  cutPatch.calcEnrichedFaces
822  (
823  pime,
824  pise,
825  projectedSlavePoints
826  );
827 
828  // Demand driven calculate the cut faces. Apart from the
829  // cutFaces/cutFaceMaster/cutFaceSlave no information from the cutPatch
830  // is used anymore!
831  const faceList& cutFaces = cutPatch.cutFaces();
832  const labelList& cutFaceMaster = cutPatch.cutFaceMaster();
833  const labelList& cutFaceSlave = cutPatch.cutFaceSlave();
834 
835  const labelList& masterFc = masterFaceCells();
836  const labelList& slaveFc = slaveFaceCells();
837 
838  // Couple the interface
839 
840  // Algorithm:
841  // 1) Go through the cut faces and check if the cut face is the same as the
842  // defining master or slave. If so, modify the appropriate
843  // face and mark the other for relegation into the face zone.
844  // 2) If different, mark both sides for relegation and insert the new face
845 
846 
847  boolList orphanedMaster(masterPatch.size(), false);
848  boolList orphanedSlave(slavePatch.size(), false);
849 
850  forAll(cutFaces, facei)
851  {
852  const face& curCutFace = cutFaces[facei];
853  const label curMaster = cutFaceMaster[facei];
854  const label curSlave = cutFaceSlave[facei];
855 
856 // Pout<< "Doing insertion of face " << facei << ": ";
857 
858  // Check if the face has changed topologically
859  bool insertedFace = false;
860 
861  if (curMaster >= 0)
862  {
863  // Face has got a master
864  if (curCutFace == masterPatch[curMaster])
865  {
866  // Face is equal to master. Modify master face.
867 // Pout<< "Face is equal to master and is ";
868 
869  // If the face has got both master and slave, it is an
870  // internal face; otherwise it is a patch face in the
871  // master patch. Keep it in the master face zone.
872 
873  if (curSlave >= 0)
874  {
875 // Pout<< "internal" << endl;
876  if (masterFc[curMaster] < slaveFc[curSlave])
877  {
878  // Cut face should point into slave.
879  // Be careful about flips in zone!
880  ref.setAction
881  (
882  polyModifyFace
883  (
884  curCutFace, // new face
885  masterPatchAddr[curMaster], // master face id
886  masterFc[curMaster], // owner
887  slaveFc[curSlave], // neighbour
888  false, // flux flip
889  -1, // patch ID
890  false, // remove from zone
891  masterFaceZoneID_.index(), // zone ID
892  masterPatchFlip[curMaster] // zone flip
893  )
894  );
895 
896  // Pout<< "modifying master face. Old master: "
897  // << masterPatch[curMaster]
898  // << " new face: " << curCutFace.reverseFace()
899  // << " own: " << masterFc[curMaster]
900  // << " nei: " << slaveFc[curSlave] << endl;
901  }
902  else
903  {
904  // Cut face should point into master. Flip required.
905  // Be careful about flips in zone!
906  ref.setAction
907  (
908  polyModifyFace
909  (
910  curCutFace.reverseFace(), // new face
911  masterPatchAddr[curMaster], // master face id
912  slaveFc[curSlave], // owner
913  masterFc[curMaster], // neighbour
914  true, // flux flip
915  -1, // patch ID
916  false, // remove from zone
917  masterFaceZoneID_.index(), // zone ID
918  !masterPatchFlip[curMaster] // zone flip
919  )
920  );
921  }
922 
923  // Orphan the slave
924  orphanedSlave[curSlave] = true;
925  }
926  else
927  {
928 // Pout<< "master boundary" << endl;
929  ref.setAction
930  (
931  polyModifyFace
932  (
933  curCutFace, // new face
934  masterPatchAddr[curMaster], // master face index
935  masterFc[curMaster], // owner
936  -1, // neighbour
937  false, // flux flip
938  masterPatchID_.index(), // patch ID
939  false, // remove from zone
940  masterFaceZoneID_.index(), // zone ID
941  masterPatchFlip[curMaster] // zone flip
942  )
943  );
944  }
945 
946  insertedFace = true;
947  }
948  }
949  else if (curSlave >= 0)
950  {
951  // Face has got a slave
952 
953  // Renumber the slave face into merged labels
954  face rsf(slavePatch[curSlave]);
955 
956  forAll(rsf, i)
957  {
958  Map<label>::const_iterator mpIter = pointMergeMap.find(rsf[i]);
959 
960  if (mpIter != pointMergeMap.end())
961  {
962  rsf[i] = mpIter();
963  }
964  }
965 
966  if (curCutFace == rsf)
967  {
968  // Face is equal to slave. Modify slave face.
969  // Pout<< "Face is equal to slave and is ";
970 
971  if (curMaster >= 0)
972  {
973  // Pout<< "regular internal" << endl;
974  if (masterFc[curMaster] < slaveFc[curSlave])
975  {
976  ref.setAction
977  (
978  polyModifyFace
979  (
980  curCutFace, // new face
981  slavePatchAddr[curSlave], // master face id
982  masterFc[curMaster], // owner
983  slaveFc[curSlave], // neighbour
984  true, // flux flip
985  -1, // patch ID
986  false, // remove from zone
987  slaveFaceZoneID_.index(), // zone ID
988  !slavePatchFlip[curMaster] // zone flip
989  )
990  );
991  }
992  else
993  {
994  // Cut face should point into master.
995  // Be careful about flips in zone!
996  // Pout<< "flipped internal" << endl;
997  ref.setAction
998  (
999  polyModifyFace
1000  (
1001  curCutFace.reverseFace(), // new face
1002  slavePatchAddr[curSlave], // master face id
1003  slaveFc[curSlave], // owner
1004  masterFc[curMaster], // neighbour
1005  false, // flux flip
1006  -1, // patch ID
1007  false, // remove from zone
1008  slaveFaceZoneID_.index(), // zone ID
1009  slavePatchFlip[curSlave] // zone flip
1010  )
1011  );
1012  }
1013 
1014  // Orphan the master
1015  orphanedMaster[curMaster] = true;
1016  }
1017  else
1018  {
1019  // Pout<< "slave boundary" << endl;
1020  ref.setAction
1021  (
1022  polyModifyFace
1023  (
1024  curCutFace, // new face
1025  slavePatchAddr[curSlave], // master face index
1026  slaveFc[curSlave], // owner
1027  -1, // neighbour
1028  false, // flux flip
1029  slavePatchID_.index(), // patch ID
1030  false, // remove from zone
1031  slaveFaceZoneID_.index(), // zone ID
1032  slavePatchFlip[curSlave] // zone flip
1033  )
1034  );
1035  }
1036 
1037  insertedFace = true;
1038  }
1039  }
1040  else
1041  {
1043  << "Face " << facei << " in cut faces has neither a master "
1044  << "nor a slave. Error in the cutting algorithm on modify."
1045  << abort(FatalError);
1046  }
1047 
1048  if (!insertedFace)
1049  {
1050  // Face is different from both master and slave
1051  // Pout<< "Face different from both master and slave" << endl;
1052 
1053  if (curMaster >= 0)
1054  {
1055  if (curSlave >= 0)
1056  {
1057  // Add internal face
1058  if (masterFc[curMaster] < slaveFc[curSlave])
1059  {
1060  // Pout<< "Adding internal face " << curCutFace
1061  // << " owner: " << masterFc[curMaster]
1062  // << " slave: " << slaveFc[curSlave]
1063  // << " master face: " << masterPatchAddr[curMaster]
1064  // << endl;
1065 
1066  // Cut face should point into slave.
1067  ref.setAction
1068  (
1069  polyAddFace
1070  (
1071  curCutFace, // new face
1072  masterFc[curMaster], // owner
1073  slaveFc[curSlave], // neighbour
1074  -1, // master point
1075  -1, // master edge
1076  masterPatchAddr[curMaster], // master face id
1077  false, // flux flip
1078  -1, // patch ID
1079  cutFaceZoneID_.index(), // zone ID
1080  false // zone flip
1081  )
1082  );
1083  }
1084  else
1085  {
1086  // Cut face should point into master. Flip required.
1087  ref.setAction
1088  (
1089  polyAddFace
1090  (
1091  curCutFace.reverseFace(), // new face
1092  slaveFc[curSlave], // owner
1093  masterFc[curMaster], // neighbour
1094  -1, // master point
1095  -1, // master edge
1096  masterPatchAddr[curMaster], // master face id
1097  true, // flux flip
1098  -1, // patch ID
1099  cutFaceZoneID_.index(), // zone ID
1100  true // zone flip
1101  )
1102  );
1103  }
1104 
1105  // Orphan slave
1106  orphanedSlave[curSlave] = true;
1107  }
1108  else
1109  {
1110  // Pout<< "Adding solo master face " << curCutFace
1111  // << " owner: " << masterFc[curMaster]
1112  // << " master face: " << masterPatchAddr[curMaster]
1113  // << endl;
1114 
1115  // Add master patch face
1116  ref.setAction
1117  (
1118  polyAddFace
1119  (
1120  curCutFace, // new face
1121  masterFc[curMaster], // owner
1122  -1, // neighbour
1123  -1, // master point
1124  -1, // master edge
1125  masterPatchAddr[curMaster], // master face index
1126  false, // flux flip
1127  masterPatchID_.index(), // patch ID
1128  cutFaceZoneID_.index(), // zone ID
1129  false // zone flip
1130  )
1131  );
1132  }
1133 
1134  // Orphan master
1135  orphanedMaster[curMaster] = true;
1136  }
1137  else if (curSlave >= 0)
1138  {
1139  // Pout<< "Adding solo slave face " << curCutFace
1140  // << " owner: " << slaveFc[curSlave]
1141  // << " master face: " << slavePatchAddr[curSlave]
1142  // << endl;
1143 
1144  // Add slave patch face
1145  ref.setAction
1146  (
1147  polyAddFace
1148  (
1149  curCutFace, // new face
1150  slaveFc[curSlave], // owner
1151  -1, // neighbour
1152  -1, // master point
1153  -1, // master edge
1154  slavePatchAddr[curSlave], // master face index
1155  false, // flux flip
1156  slavePatchID_.index(), // patch ID
1157  cutFaceZoneID_.index(), // zone ID
1158  false // zone flip
1159  )
1160  );
1161 
1162  // Orphan slave
1163  orphanedSlave[curSlave] = true;
1164  }
1165  else
1166  {
1168  << "Face " << facei << " in cut faces has neither a master "
1169  << "nor a slave. Error in the cutting algorithm on add."
1170  << abort(FatalError);
1171  }
1172  }
1173  }
1174 
1175  // Move the orphaned faces into the face zone
1176  // Pout<< "Orphaned master faces: " << orphanedMaster << endl;
1177  // Pout<< "Orphaned slave faces: " << orphanedSlave << endl;
1178 
1179  label nOrphanedMasters = 0;
1180 
1181  forAll(orphanedMaster, facei)
1182  {
1183  if (orphanedMaster[facei])
1184  {
1185  nOrphanedMasters++;
1186 
1188  // ref.setAction
1189  //(
1190  // polyModifyFace
1191  // (
1192  // masterPatch[facei], // new face
1193  // masterPatchAddr[facei], // master face index
1194  // -1, // owner
1195  // -1, // neighbour
1196  // false, // flux flip
1197  // -1, // patch ID
1198  // false, // remove from zone
1199  // masterFaceZoneID_.index(), // zone ID
1200  // false // zone flip
1201  // )
1202  //);
1203 
1204  // Pout<< "**MJ:deleting master face " << masterPatchAddr[facei]
1205  // << " old verts:" << masterPatch[facei] << endl;
1206  ref.setAction(polyRemoveFace(masterPatchAddr[facei]));
1207  }
1208  }
1209 
1210  label nOrphanedSlaves = 0;
1211 
1212  forAll(orphanedSlave, facei)
1213  {
1214  if (orphanedSlave[facei])
1215  {
1216  nOrphanedSlaves++;
1217 
1219  // ref.setAction
1220  //(
1221  // polyModifyFace
1222  // (
1223  // slavePatch[facei], // new face
1224  // slavePatchAddr[facei], // slave face index
1225  // -1, // owner
1226  // -1, // neighbour
1227  // false, // flux flip
1228  // -1, // patch ID
1229  // false, // remove from zone
1230  // slaveFaceZoneID_.index(), // zone ID
1231  // false // zone flip
1232  // )
1233  //);
1234 
1235  // Pout<< "**MJ:deleting slave face " << slavePatchAddr[facei]
1236  // << " old verts:" << slavePatch[facei] << endl;
1237  ref.setAction(polyRemoveFace(slavePatchAddr[facei]));
1238  }
1239  }
1240 
1241  if (debug)
1242  {
1243  Pout<< "Number of orphaned faces: "
1244  << "master = " << nOrphanedMasters << " out of "
1245  << orphanedMaster.size()
1246  << " slave = " << nOrphanedSlaves << " out of "
1247  << orphanedSlave.size() << endl;
1248  }
1249 
1250  // Finished coupling the plane of sliding interface.
1251 
1252  // Modify faces influenced by the sliding interface
1253  // These are the faces belonging to the master and slave cell
1254  // layer which have not already been modified.
1255  // Modification comes in three types:
1256  // 1) eliminate the points that have been removed when the old sliding
1257  // interface has been removed
1258  // 2) Merge the slave points that have hit points on the master patch
1259  // 3) Introduce new points resulting from the new sliding interface cut
1260 
1261  // Collect all affected faces
1262 
1263  // Master side
1264 
1265  // Grab the list of faces in the layer
1266  const labelList& masterStickOuts = masterStickOutFaces();
1267 
1268  // Pout<< "masterStickOuts: " << masterStickOuts << endl;
1269 
1270  // Re-create the master stick-out faces
1271  forAll(masterStickOuts, facei)
1272  {
1273  // Renumber the face and remove additional points
1274 
1275  const label curFaceID = masterStickOuts[facei];
1276 
1277  const face& oldRichFace = faces[curFaceID];
1278 
1279  bool changed = false;
1280 
1281  // Remove removed points from the face
1282  face oldFace(oldRichFace.size());
1283  label nOldFace = 0;
1284 
1285  forAll(oldRichFace, pointi)
1286  {
1287  if (ref.pointRemoved(oldRichFace[pointi]))
1288  {
1289  changed = true;
1290  }
1291  else
1292  {
1293  // Point off patch
1294  oldFace[nOldFace] = oldRichFace[pointi];
1295  nOldFace++;
1296  }
1297  }
1298 
1299  oldFace.setSize(nOldFace);
1300 
1301  // Pout<< "old rich master face: " << oldRichFace
1302  // << " old face: " << oldFace
1303  // << endl;
1304 
1305  DynamicList<label> newFaceLabels(2*oldFace.size());
1306 
1307  forAll(oldFace, pointi)
1308  {
1309  if (masterMeshPointMap.found(oldFace[pointi]))
1310  {
1311  // Point is in master patch. Add it
1312 
1313  // If the point is a direct hit, grab its label; otherwise
1314  // keep the original
1315  if (pointMergeMap.found(oldFace[pointi]))
1316  {
1317  changed = true;
1318  newFaceLabels.append
1319  (
1320  pointMergeMap.find(oldFace[pointi])()
1321  );
1322  }
1323  else
1324  {
1325  newFaceLabels.append(oldFace[pointi]);
1326  }
1327 
1328  // Find if there are additional points inserted onto the edge
1329  // between current point and the next point
1330  // Algorithm:
1331  // 1) Find all the edges in the master patch coming
1332  // out of the current point.
1333  // 2) If the next point in the face to pick the right edge
1334  const label localFirstLabel =
1335  masterMeshPointMap.find(oldFace[pointi])();
1336 
1337  const labelList& curEdges = masterPointEdges[localFirstLabel];
1338 
1339  const label nextLabel = oldFace.nextLabel(pointi);
1340 
1341  Map<label>::const_iterator mmpmIter =
1342  masterMeshPointMap.find(nextLabel);
1343 
1344  if (mmpmIter != masterMeshPointMap.end())
1345  {
1346  // Pout<< "found label pair " << oldFace[pointi]
1347  // << " and " << nextLabel;
1348  // Find the points on the edge between them
1349  const label localNextLabel = mmpmIter();
1350 
1351  forAll(curEdges, curEdgeI)
1352  {
1353  if
1354  (
1355  masterEdges[curEdges[curEdgeI]].otherVertex
1356  (
1357  localFirstLabel
1358  )
1359  == localNextLabel
1360  )
1361  {
1362  // Pout<< " found edge: " << curEdges[curEdgeI]
1363  // << endl;
1364 
1365  // Get points on current edge
1366  const labelList& curPime = pime[curEdges[curEdgeI]];
1367 
1368  if (curPime.size())
1369  {
1370  changed = true;
1371  // Pout<< "curPime: " << curPime << endl;
1372  // Insert the edge points into the face
1373  // in the correct order
1374  const point& startPoint =
1375  masterLocalPoints[localFirstLabel];
1376 
1377  vector e =
1378  masterLocalPoints[localNextLabel]
1379  - startPoint;
1380 
1381  e /= magSqr(e);
1382 
1383  scalarField edgePointWeights(curPime.size());
1384 
1385  forAll(curPime, curPimeI)
1386  {
1387  edgePointWeights[curPimeI] =
1388  (
1389  e
1390  & (
1391  pointMap.find(curPime[curPimeI])()
1392  - startPoint
1393  )
1394  );
1395  }
1396 
1397  if (debug)
1398  {
1399  if
1400  (
1401  min(edgePointWeights) < 0
1402  || max(edgePointWeights) > 1
1403  )
1404  {
1406  << "Error in master stick-out edge "
1407  << "point collection."
1408  << abort(FatalError);
1409  }
1410  }
1411 
1412  // Go through the points and collect
1413  // them based on weights from lower to
1414  // higher. This gives the correct
1415  // order of points along the edge.
1416  for
1417  (
1418  label passI = 0;
1419  passI < edgePointWeights.size();
1420  passI++
1421  )
1422  {
1423  // Max weight can only be one, so
1424  // the sorting is done by
1425  // elimination.
1426  label nextPoint = -1;
1427  scalar dist = 2;
1428 
1429  forAll(edgePointWeights, wI)
1430  {
1431  if (edgePointWeights[wI] < dist)
1432  {
1433  dist = edgePointWeights[wI];
1434  nextPoint = wI;
1435  }
1436  }
1437 
1438  // Insert the next point and reset
1439  // its weight to exclude it from
1440  // future picks
1441  newFaceLabels.append(curPime[nextPoint]);
1442  edgePointWeights[nextPoint] = great;
1443  }
1444  }
1445 
1446  break;
1447  } // End of found edge
1448  } // End of looking through current edges
1449  }
1450  }
1451  else
1452  {
1453  newFaceLabels.append(oldFace[pointi]);
1454  }
1455  }
1456 
1457  if (changed)
1458  {
1459  if (newFaceLabels.size() < 3)
1460  {
1462  << "Face " << curFaceID << " reduced to less than "
1463  << "3 points. Topological/cutting error A." << nl
1464  << "Old face: " << oldFace << " new face: " << newFaceLabels
1465  << abort(FatalError);
1466  }
1467 
1468  // Get face zone and its flip
1469  label modifiedFaceZone = faceZones.whichZone(curFaceID);
1470  bool modifiedFaceZoneFlip = false;
1471 
1472  if (modifiedFaceZone >= 0)
1473  {
1474  modifiedFaceZoneFlip =
1475  faceZones[modifiedFaceZone].flipMap()
1476  [
1477  faceZones[modifiedFaceZone].whichFace(curFaceID)
1478  ];
1479  }
1480 
1481  face newFace;
1482  newFace.transfer(newFaceLabels);
1483 
1484  // Pout<< "Modifying master stick-out face " << curFaceID
1485  // << " old face: " << oldFace
1486  // << " new face: " << newFace
1487  // << endl;
1488 
1489  // Modify the face
1490  if (mesh.isInternalFace(curFaceID))
1491  {
1492  ref.setAction
1493  (
1494  polyModifyFace
1495  (
1496  newFace, // modified face
1497  curFaceID, // label of face being modified
1498  own[curFaceID], // owner
1499  nei[curFaceID], // neighbour
1500  false, // face flip
1501  -1, // patch for face
1502  false, // remove from zone
1503  modifiedFaceZone, // zone for face
1504  modifiedFaceZoneFlip // face flip in zone
1505  )
1506  );
1507  }
1508  else
1509  {
1510  ref.setAction
1511  (
1512  polyModifyFace
1513  (
1514  newFace, // modified face
1515  curFaceID, // label of face being modified
1516  own[curFaceID], // owner
1517  -1, // neighbour
1518  false, // face flip
1519  mesh.boundaryMesh().whichPatch(curFaceID),
1520  // patch for face
1521  false, // remove from zone
1522  modifiedFaceZone, // zone for face
1523  modifiedFaceZoneFlip // face flip in zone
1524  )
1525  );
1526  }
1527  }
1528  }
1529 
1530  // Pout<< "Finished master side" << endl;
1531 
1532  // Slave side
1533 
1534  // Grab the list of faces in the layer
1535  const labelList& slaveStickOuts = slaveStickOutFaces();
1536 
1537  // Pout<< "slaveStickOuts: " << slaveStickOuts << endl;
1538 
1539  const Map<label>& rpm = retiredPointMap();
1540 
1541  // Re-create the slave stick-out faces
1542 
1543  forAll(slaveStickOuts, facei)
1544  {
1545  // Renumber the face and remove additional points
1546  const label curFaceID = slaveStickOuts[facei];
1547 
1548  const face& oldRichFace = faces[curFaceID];
1549 
1550  bool changed = false;
1551 
1552  // Remove removed points from the face
1553  face oldFace(oldRichFace.size());
1554  label nOldFace = 0;
1555 
1556  forAll(oldRichFace, pointi)
1557  {
1558  if
1559  (
1560  rpm.found(oldRichFace[pointi])
1561  || slaveMeshPointMap.found(oldRichFace[pointi])
1562  )
1563  {
1564  // Point definitely live. Add it
1565  oldFace[nOldFace] = oldRichFace[pointi];
1566  nOldFace++;
1567  }
1568  else if
1569  (
1570  ref.pointRemoved(oldRichFace[pointi])
1571  || masterMeshPointMap.found(oldRichFace[pointi])
1572  )
1573  {
1574  // Point removed and not on slave patch
1575  // (first if takes care of that!) or
1576  // point belonging to master patch
1577  changed = true;
1578  }
1579  else
1580  {
1581  // Point off patch
1582  oldFace[nOldFace] = oldRichFace[pointi];
1583  nOldFace++;
1584  }
1585  }
1586 
1587  oldFace.setSize(nOldFace);
1588 
1589  DynamicList<label> newFaceLabels(2*oldFace.size());
1590 
1591  // Pout<< "old rich slave face: " << oldRichFace
1592  // << " old face: " << oldFace
1593  // << endl;
1594 
1595  forAll(oldFace, pointi)
1596  {
1597  // Try to find the point in retired points
1598  label curP = oldFace[pointi];
1599 
1600  Map<label>::const_iterator rpmIter = rpm.find(oldFace[pointi]);
1601 
1602  if (rpmIter != rpm.end())
1603  {
1604  changed = true;
1605  curP = rpmIter();
1606  }
1607 
1608  if (slaveMeshPointMap.found(curP))
1609  {
1610  // Point is in slave patch. Add it
1611 
1612  // If the point is a direct hit, grab its label; otherwise
1613  // keep the original
1614  if (pointMergeMap.found(curP))
1615  {
1616  changed = true;
1617  newFaceLabels.append
1618  (
1619  pointMergeMap.find(curP)()
1620  );
1621  }
1622  else
1623  {
1624  newFaceLabels.append(curP);
1625  }
1626 
1627  // Find if there are additional points inserted onto the edge
1628  // between current point and the next point
1629  // Algorithm:
1630  // 1) Find all the edges in the slave patch coming
1631  // out of the current point.
1632  // 2) Use the next point in the face to pick the right edge
1633 
1634  const label localFirstLabel =
1635  slaveMeshPointMap.find(curP)();
1636 
1637  const labelList& curEdges = slavePointEdges[localFirstLabel];
1638 
1639  label nextLabel = oldFace.nextLabel(pointi);
1640 
1641  Map<label>::const_iterator rpmNextIter =
1642  rpm.find(nextLabel);
1643 
1644  if (rpmNextIter != rpm.end())
1645  {
1646  nextLabel = rpmNextIter();
1647  }
1648 
1649  Map<label>::const_iterator mmpmIter =
1650  slaveMeshPointMap.find(nextLabel);
1651 
1652  if (mmpmIter != slaveMeshPointMap.end())
1653  {
1654  // Both points on the slave patch.
1655  // Find the points on the edge between them
1656  const label localNextLabel = mmpmIter();
1657 
1658  forAll(curEdges, curEdgeI)
1659  {
1660  if
1661  (
1662  slaveEdges[curEdges[curEdgeI]].otherVertex
1663  (
1664  localFirstLabel
1665  )
1666  == localNextLabel
1667  )
1668  {
1669  // Pout<< " found edge: " << curEdges[curEdgeI]
1670  // << endl;
1671 
1672  // Get points on current edge
1673  const labelList& curPise = pise[curEdges[curEdgeI]];
1674 
1675  if (curPise.size())
1676  {
1677  changed = true;
1678  // Pout<< "curPise: " << curPise << endl;
1679  // Insert the edge points into the face
1680  // in the correct order
1681  const point& startPoint =
1682  projectedSlavePoints[localFirstLabel];
1683 
1684  vector e =
1685  projectedSlavePoints[localNextLabel]
1686  - startPoint;
1687 
1688  e /= magSqr(e);
1689 
1690  scalarField edgePointWeights(curPise.size());
1691 
1692  forAll(curPise, curPiseI)
1693  {
1694  edgePointWeights[curPiseI] =
1695  (
1696  e
1697  & (
1698  pointMap.find(curPise[curPiseI])()
1699  - startPoint
1700  )
1701  );
1702  }
1703 
1704  if (debug)
1705  {
1706  if
1707  (
1708  min(edgePointWeights) < 0
1709  || max(edgePointWeights) > 1
1710  )
1711  {
1713  << "Error in slave stick-out edge "
1714  << "point collection."
1715  << abort(FatalError);
1716  }
1717  }
1718 
1719  // Go through the points and collect
1720  // them based on weights from lower to
1721  // higher. This gives the correct
1722  // order of points along the edge.
1723  for
1724  (
1725  label passI = 0;
1726  passI < edgePointWeights.size();
1727  passI++
1728  )
1729  {
1730  // Max weight can only be one, so
1731  // the sorting is done by
1732  // elimination.
1733  label nextPoint = -1;
1734  scalar dist = 2;
1735 
1736  forAll(edgePointWeights, wI)
1737  {
1738  if (edgePointWeights[wI] < dist)
1739  {
1740  dist = edgePointWeights[wI];
1741  nextPoint = wI;
1742  }
1743  }
1744 
1745  // Insert the next point and reset
1746  // its weight to exclude it from
1747  // future picks
1748  newFaceLabels.append(curPise[nextPoint]);
1749  edgePointWeights[nextPoint] = great;
1750  }
1751  }
1752 
1753  break;
1754  }
1755  } // End of found edge
1756  } // End of looking through current edges
1757  }
1758  else
1759  {
1760  newFaceLabels.append(oldFace[pointi]);
1761  }
1762  }
1763 
1764  if (changed)
1765  {
1766  if (newFaceLabels.size() < 3)
1767  {
1769  << "Face " << curFaceID << " reduced to less than "
1770  << "3 points. Topological/cutting error B." << nl
1771  << "Old face: " << oldFace << " new face: " << newFaceLabels
1772  << abort(FatalError);
1773  }
1774 
1775  // Get face zone and its flip
1776  label modifiedFaceZone = faceZones.whichZone(curFaceID);
1777  bool modifiedFaceZoneFlip = false;
1778 
1779  if (modifiedFaceZone >= 0)
1780  {
1781  modifiedFaceZoneFlip =
1782  faceZones[modifiedFaceZone].flipMap()
1783  [
1784  faceZones[modifiedFaceZone].whichFace(curFaceID)
1785  ];
1786  }
1787 
1788  face newFace;
1789  newFace.transfer(newFaceLabels);
1790 
1791  // Pout<< "Modifying slave stick-out face " << curFaceID
1792  // << " old face: " << oldFace
1793  // << " new face: " << newFace
1794  // << endl;
1795 
1796  // Modify the face
1797  if (mesh.isInternalFace(curFaceID))
1798  {
1799  ref.setAction
1800  (
1801  polyModifyFace
1802  (
1803  newFace, // modified face
1804  curFaceID, // label of face being modified
1805  own[curFaceID], // owner
1806  nei[curFaceID], // neighbour
1807  false, // face flip
1808  -1, // patch for face
1809  false, // remove from zone
1810  modifiedFaceZone, // zone for face
1811  modifiedFaceZoneFlip // face flip in zone
1812  )
1813  );
1814  }
1815  else
1816  {
1817  ref.setAction
1818  (
1819  polyModifyFace
1820  (
1821  newFace, // modified face
1822  curFaceID, // label of face being modified
1823  own[curFaceID], // owner
1824  -1, // neighbour
1825  false, // face flip
1826  mesh.boundaryMesh().whichPatch(curFaceID),
1827  // patch for face
1828  false, // remove from zone
1829  modifiedFaceZone, // zone for face
1830  modifiedFaceZoneFlip // face flip in zone
1831  )
1832  );
1833  }
1834  }
1835  }
1836 
1837  // Activate and retire slave patch points
1838  // This needs to be done last, so that the map of removed points
1839  // does not get damaged by point modifications
1840 
1841  if (!retiredPointMapPtr_)
1842  {
1844  << "Retired point map pointer not set."
1845  << abort(FatalError);
1846  }
1847 
1848  Map<label>& addToRpm = *retiredPointMapPtr_;
1849 
1850  // Clear the old map
1851  addToRpm.clear();
1852 
1853  label nRetiredPoints = 0;
1854 
1855  forAll(slaveMeshPoints, pointi)
1856  {
1857  if (pointMergeMap.found(slaveMeshPoints[pointi]))
1858  {
1859  // Retire the point - only used for supporting the detached
1860  // slave patch
1861  nRetiredPoints++;
1862 
1863  // ref.setAction
1864  // (
1865  // polyModifyPoint
1866  // (
1867  // slaveMeshPoints[pointi], // point ID
1868  // points[slaveMeshPoints[pointi]], // point
1869  // false, // remove from zone
1870  // mesh.pointZones().whichZone(slaveMeshPoints[pointi]),
1871  // // zone
1872  // false // in a cell
1873  // )
1874  // );
1875  // Pout<< "MJ retire slave point " << slaveMeshPoints[pointi]
1876  // << " coord " << points[slaveMeshPoints[pointi]]
1877  // << endl;
1878  ref.setAction
1879  (
1880  polyRemovePoint
1881  (
1882  slaveMeshPoints[pointi]
1883  )
1884  );
1885 
1886  addToRpm.insert
1887  (
1888  pointMergeMap.find(slaveMeshPoints[pointi])(),
1889  slaveMeshPoints[pointi]
1890  );
1891  }
1892  else
1893  {
1894  ref.setAction
1895  (
1896  polyModifyPoint
1897  (
1898  slaveMeshPoints[pointi], // point ID
1899  points[slaveMeshPoints[pointi]], // point
1900  false, // remove from zone
1901  mesh.pointZones().whichZone(slaveMeshPoints[pointi]),// zone
1902  true // in a cell
1903  )
1904  );
1905  }
1906  }
1907 
1908  if (debug)
1909  {
1910  Pout<< "Retired " << nRetiredPoints << " out of "
1911  << slaveMeshPoints.size() << " points." << endl;
1912  }
1913 
1914  // Grab cut face master and slave addressing
1915  if (cutFaceMasterPtr_) deleteDemandDrivenData(cutFaceMasterPtr_);
1916  cutFaceMasterPtr_ = new labelList(cutPatch.cutFaceMaster());
1917 
1918  if (cutFaceSlavePtr_) deleteDemandDrivenData(cutFaceSlavePtr_);
1919  cutFaceSlavePtr_ = new labelList(cutPatch.cutFaceSlave());
1920 
1921  // Finished coupling
1922  attached_ = true;
1923 
1924  if (debug)
1925  {
1926  Pout<< "void slidingInterface::coupleInterface("
1927  << "polyTopoChange& ref) : "
1928  << "Finished coupling sliding interface " << name() << endl;
1929  }
1930 }
1931 
1932 
1933 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
dimensionedSymmTensor sqr(const dimensionedVector &dv)
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
HashTable< label, label, Hash< label > >::const_iterator const_iterator
Definition: Map.H:59
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static const unsigned edgesPerFace_
Estimated number of edges per cell.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
List< edge > edgeList
Definition: edgeList.H:38
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:124
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void clear()
Clear all entries from table.
Definition: HashTable.C:468
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label index() const
Return index of first matching zone.
Definition: DynamicID.H:107
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static const char nl
Definition: Ostream.H:265
PrimitivePatch< face, List, const pointField & > primitiveFacePatch
Foam::primitiveFacePatch.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const polyTopoChanger & topoChanger() const
Return reference to morph engine.
vector point
Point is a vector.
Definition: point.H:41
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const word & name() const
Return name of this modifier.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
PointHit< point > pointHit
Definition: pointHit.H:41
const polyMesh & mesh() const
Return the mesh reference.
void deleteDemandDrivenData(DataPtr &dataPtr)