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