1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "slidingInterface.H"
28 #include "polyTopoChange.H"
30 #include "primitiveMesh.H"
31 #include "enrichedPatch.H"
32 #include "DynamicList.H"
34 #include "triPointRef.H"
36 #include "polyTopoChanger.H"
37 #include "polyAddPoint.H"
38 #include "polyRemovePoint.H"
39 #include "polyAddFace.H"
40 #include "polyModifyPoint.H"
41 #include "polyModifyFace.H"
42 #include "polyRemoveFace.H"
44 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
46 const Foam::scalar Foam::slidingInterface::edgeCoPlanarTol_ = 0.8;
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
50 // Index of debug signs:
51 // Index of debug signs:
52 // . - loop of the edge-to-face interaction detection
53 // x - reversal of direction in edge-to-face interaction detection
54 // + - complete edge-to-face interaction detection
55 // z - incomplete edge-to-face interaction detection. This may be OK for edges
56 // crossing from one to the other side of multiply connected master patch
58 // e - adding a point on edge
59 // f - adding a point on face
60 // . - collecting edges off another face for edge-to-edge cut
61 // + - finished collection of edges
62 // * - cut both master and slave
63 // n - cutting new edge
64 // t - intersection exists but it is outside of tolerance
65 // x - missed slave edge in cut
66 // - - missed master edge in cut
67 // u - edge already used in cutting
69 void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
73 Pout<< "void slidingInterface::coupleInterface"
74 << "(polyTopoChange& ref) : "
75 << "Coupling sliding interface " << name() << endl;
78 const polyMesh& mesh = topoChanger().mesh();
80 const pointField& points = mesh.points();
81 const faceList& faces = mesh.faces();
83 const labelList& own = mesh.faceOwner();
84 const labelList& nei = mesh.faceNeighbour();
85 const faceZoneMesh& faceZones = mesh.faceZones();
87 const primitiveFacePatch& masterPatch =
88 faceZones[masterFaceZoneID_.index()]();
90 const labelList& masterPatchAddr = faceZones[masterFaceZoneID_.index()];
92 const boolList& masterPatchFlip =
93 faceZones[masterFaceZoneID_.index()].flipMap();
95 const primitiveFacePatch& slavePatch =
96 faceZones[slaveFaceZoneID_.index()]();
98 const labelList& slavePatchAddr = faceZones[slaveFaceZoneID_.index()];
100 const boolList& slavePatchFlip =
101 faceZones[slaveFaceZoneID_.index()].flipMap();
103 const edgeList& masterEdges = masterPatch.edges();
104 const labelListList& masterPointEdges = masterPatch.pointEdges();
105 const labelList& masterMeshPoints = masterPatch.meshPoints();
106 const pointField& masterLocalPoints = masterPatch.localPoints();
107 const labelListList& masterFaceFaces = masterPatch.faceFaces();
108 const labelListList& masterFaceEdges = masterPatch.faceEdges();
109 const Map<label>& masterMeshPointMap = masterPatch.meshPointMap();
111 const edgeList& slaveEdges = slavePatch.edges();
112 const labelListList& slavePointEdges = slavePatch.pointEdges();
113 const labelList& slaveMeshPoints = slavePatch.meshPoints();
114 const pointField& slaveLocalPoints = slavePatch.localPoints();
115 const Map<label>& slaveMeshPointMap = slavePatch.meshPointMap();
116 const vectorField& slavePointNormals = slavePatch.pointNormals();
118 // Collect projection addressing
122 slavePointPointHitsPtr_
123 && slavePointEdgeHitsPtr_
124 && slavePointFaceHitsPtr_
125 && masterPointEdgeHitsPtr_
131 "void slidingInterface::coupleInterface("
132 "polyTopoChange& ref) const"
133 ) << "Point projection addressing not available."
134 << abort(FatalError);
137 const labelList& slavePointPointHits = *slavePointPointHitsPtr_;
138 const labelList& slavePointEdgeHits = *slavePointEdgeHitsPtr_;
139 const List<objectHit>& slavePointFaceHits = *slavePointFaceHitsPtr_;
140 const labelList& masterPointEdgeHits = *masterPointEdgeHitsPtr_;
141 const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
143 // Create enriched patch
144 enrichedPatch cutPatch
153 // Get reference to list of added point for the enriched patch.
154 // This will be used for point addition
155 Map<point>& pointMap = cutPatch.pointMap();
157 // Get reference to the list of merged points
158 Map<label>& pointMergeMap = cutPatch.pointMergeMap();
160 // Create mapping for every merged point of the slave patch
161 forAll (slavePointPointHits, pointI)
163 if (slavePointPointHits[pointI] >= 0)
165 // Pout << "Inserting point merge pair: " << slaveMeshPoints[pointI] << " : " << masterMeshPoints[slavePointPointHits[pointI]] << endl;
168 slaveMeshPoints[pointI],
169 masterMeshPoints[slavePointPointHits[pointI]]
174 // Collect the list of used edges for every slave edge
176 List<labelHashSet> usedMasterEdges(slaveEdges.size());
178 // Collect of slave point hits
179 forAll (slavePointPointHits, pointI)
181 // For point hits, add all point-edges from master side to all point
182 const labelList& curSlaveEdges = slavePointEdges[pointI];
184 if (slavePointPointHits[pointI] > -1)
186 const labelList& curMasterEdges =
187 masterPointEdges[slavePointPointHits[pointI]];
189 // Mark all current master edges as used for all the current slave
191 forAll (curSlaveEdges, slaveEdgeI)
193 labelHashSet& sm = usedMasterEdges[curSlaveEdges[slaveEdgeI]];
195 forAll (curMasterEdges, masterEdgeI)
197 sm.insert(curMasterEdges[masterEdgeI]);
201 else if (slavePointEdgeHits[pointI] > -1)
203 // For edge hits, add the master edge
204 forAll (curSlaveEdges, slaveEdgeI)
206 usedMasterEdges[curSlaveEdges[slaveEdgeI]].insert
208 slavePointEdgeHits[pointI]
214 // Collect off master point hits
215 // For every master point that has hit an edge, add all edges coming from
216 // that point to the slave edge that has been hit
217 forAll (masterPointEdgeHits, masterPointI)
219 if (masterPointEdgeHits[masterPointI] > -1)
221 const labelList& curMasterEdges = masterPointEdges[masterPointI];
224 usedMasterEdges[masterPointEdgeHits[masterPointI]];
226 forAll (curMasterEdges, masterEdgeI)
228 sm.insert(curMasterEdges[masterEdgeI]);
233 // Pout << "used edges: " << endl;
234 // forAll (usedMasterEdges, edgeI)
236 // Pout << "edge: " << edgeI << " used: " << usedMasterEdges[edgeI].toc() << endl;
239 // For every master and slave edge make a list of points to be added into
241 List<DynamicList<label> > pointsIntoMasterEdges(masterEdges.size());
242 List<DynamicList<label> > pointsIntoSlaveEdges(slaveEdges.size());
244 // Add all points from the slave patch that have hit the edge
245 forAll (slavePointEdgeHits, pointI)
247 if (slavePointEdgeHits[pointI] > -1)
249 // Create a new point on the master edge
252 masterEdges[slavePointEdgeHits[pointI]].line
255 ).nearestDist(projectedSlavePoints[pointI]).hitPoint();
262 edgeCutPoint, // point
263 slaveMeshPoints[pointI], // master point
264 cutPointZoneID_.index(), // zone for point
265 true // supports a cell
268 // Pout << "Inserting merge pair off edge: " << slaveMeshPoints[pointI] << " " << newPoint << " cut point: " << edgeCutPoint << " orig: " << slaveLocalPoints[pointI] << " proj: " << projectedSlavePoints[pointI] << endl;
269 // Add the new edge point into the merge map
270 pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);
272 pointsIntoMasterEdges[slavePointEdgeHits[pointI]].append
277 // Add the point into the enriched patch map
287 // Pout<< newPoint << " = " << edgeCutPoint << endl;
297 // Add all points from the slave patch that have hit a face
298 forAll (slavePointFaceHits, pointI)
302 slavePointPointHits[pointI] < 0
303 && slavePointEdgeHits[pointI] < 0
304 && slavePointFaceHits[pointI].hit()
312 projectedSlavePoints[pointI], // point
313 slaveMeshPoints[pointI], // master point
314 cutPointZoneID_.index(), // zone for point
315 true // supports a cell
318 // Pout << "Inserting merge pair off face: " << slaveMeshPoints[pointI] << " " << newPoint << endl;
319 // Add the new edge point into the merge map
320 pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);
322 // Add the point into the enriched patch map
326 projectedSlavePoints[pointI]
331 Pout<< "f: " << newPoint << " = "
332 << projectedSlavePoints[pointI] << endl;
337 forAll (masterPointEdgeHits, pointI)
339 if (masterPointEdgeHits[pointI] > -1)
341 pointsIntoSlaveEdges[masterPointEdgeHits[pointI]].append
343 masterMeshPoints[pointI]
348 // Cut all slave edges.
349 // Collect all master edges the slave edge interacts with. Skip
350 // all the edges already marked as used. For every unused edge,
351 // calculate the cut and insert the new point into the master and
353 // For the edge selection algorithm, see, comment in
354 // slidingInterfaceProjectPoints.C.
355 // Edge cutting algoritm:
356 // As the master patch defines the cutting surface, the newly
357 // inserted point needs to be on the master edge. Also, in 3-D
358 // the pair of edges generally misses each other rather than
359 // intersect. Therefore, the intersection is calculated using the
360 // plane the slave edge defines during projection. The plane is
361 // defined by the centrepoint of the edge in the original
362 // configuration and the projected end points. In case of flat
363 // geometries (when the three points are colinear), the plane is
364 // defined by the two projected end-points and the slave edge
365 // normal used as the in-plane vector. When the intersection
366 // point is created, it is added as a new point for both the
367 // master and the slave edge.
368 // In order to be able to re-create the points on edges in mesh
369 // motion without the topology change, the edge pair used to
370 // create the cut point needs to be recorded in
371 // cutPointEdgePairMap
373 // Note. "Processing slave edges" code is repeated twice in the
374 // sliding intergace coupling in order to allow the point
375 // projection to be done separately from the actual cutting.
376 // Please change consistently with slidingInterfaceProjectPoints.C
380 Pout << "Processing slave edges " << endl;
383 if (!cutPointEdgePairMapPtr_)
387 "void slidingInterface::coupleInterface("
388 "polyTopoChange& ref) const"
389 ) << "Cut point edge pair map pointer not set."
390 << abort(FatalError);
393 Map<Pair<edge> >& addToCpepm = *cutPointEdgePairMapPtr_;
398 // Create a map of faces the edge can interact with
399 labelHashSet curFaceMap
401 nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
404 labelHashSet addedFaces(2*primitiveMesh::edgesPerFace_);
406 forAll (slaveEdges, edgeI)
408 const edge& curEdge = slaveEdges[edgeI];
412 slavePointFaceHits[curEdge.start()].hit()
413 || slavePointFaceHits[curEdge.end()].hit()
416 labelHashSet& curUme = usedMasterEdges[edgeI];
417 // Pout<< "Doing edge " << edgeI << " curEdge: " << curEdge << " curUme: " << curUme << endl;
422 // Grab the faces for start and end points.
423 const label startFace =
424 slavePointFaceHits[curEdge.start()].hitObject();
425 const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
426 // Pout << "startFace: " << slavePointFaceHits[curEdge.start()] << " endFace: " << slavePointFaceHits[curEdge.end()] << endl;
427 // Insert the start face into the list
428 curFaceMap.insert(startFace);
429 addedFaces.insert(startFace);
430 // Pout << "curFaceMap: " << curFaceMap.toc() << endl;
432 bool completed = false;
434 while (nSweeps < edgeFaceEscapeLimit_)
438 if (addedFaces.found(endFace))
443 // Add all face neighbours of face in the map
444 const labelList cf = addedFaces.toc();
449 const labelList& curNbrs = masterFaceFaces[cf[cfI]];
451 forAll (curNbrs, nbrI)
453 if (!curFaceMap.found(curNbrs[nbrI]))
455 curFaceMap.insert(curNbrs[nbrI]);
456 addedFaces.insert(curNbrs[nbrI]);
461 if (completed) break;
476 // It is impossible to reach the end from the start, probably
477 // due to disconnected domain. Do search in opposite direction
479 label nReverseSweeps = 0;
482 addedFaces.insert(endFace);
484 while (nReverseSweeps < edgeFaceEscapeLimit_)
488 if (addedFaces.found(startFace))
493 // Add all face neighbours of face in the map
494 const labelList cf = addedFaces.toc();
499 const labelList& curNbrs = masterFaceFaces[cf[cfI]];
501 forAll (curNbrs, nbrI)
503 if (!curFaceMap.found(curNbrs[nbrI]))
505 curFaceMap.insert(curNbrs[nbrI]);
506 addedFaces.insert(curNbrs[nbrI]);
511 if (completed) break;
537 // Create a map of edges the edge can interact with
538 labelHashSet curMasterEdgesMap
540 nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
543 const labelList curFaces = curFaceMap.toc();
544 // Pout << "curFaces: " << curFaces << endl;
545 forAll (curFaces, faceI)
547 // Pout<< "face: " << curFaces[faceI] << " "
548 // << masterPatch[curFaces[faceI]]
550 // << masterPatch.localFaces()[curFaces[faceI]]
552 const labelList& me = masterFaceEdges[curFaces[faceI]];
556 curMasterEdgesMap.insert(me[meI]);
560 const labelList curMasterEdges = curMasterEdgesMap.toc();
562 // For all master edges to intersect, skip the ones
563 // already used and cut the rest with a cutting plane. If
564 // the intersection point, falls inside of both edges, it
568 // The edge cutting code is repeated in
569 // slidingInterface::modifyMotionPoints. This is done for
570 // efficiency reasons and avoids multiple creation of cutting
571 // planes. Please update both simultaneously.
573 const point& a = projectedSlavePoints[curEdge.start()];
574 const point& b = projectedSlavePoints[curEdge.end()];
579 slaveLocalPoints[curEdge.start()]
580 + slavePointNormals[curEdge.start()] // Add start normal
581 + slaveLocalPoints[curEdge.end()]
582 + slavePointNormals[curEdge.end()] // Add end normal
586 plane cutPlane(a, b, c);
587 // Pout << "a: " << a << " b: " << b << " c: " << c << " plane: " << cutPlane << endl;
589 linePointRef curSlaveLine = curEdge.line(projectedSlavePoints);
590 const scalar curSlaveLineMag = curSlaveLine.mag();
591 // Pout << "curSlaveLine: " << curSlaveLine << endl;
592 forAll (curMasterEdges, masterEdgeI)
594 if (!curUme.found(curMasterEdges[masterEdgeI]))
602 const label cmeIndex = curMasterEdges[masterEdgeI];
603 const edge& cme = masterEdges[cmeIndex];
604 // Pout<< "Edge " << cmeIndex << " cme: " << cme << " line: " << cme.line(masterLocalPoints) << endl;
606 cutPlane.lineIntersect
608 cme.line(masterLocalPoints)
613 cutOnMaster > edgeEndCutoffTol_
614 && cutOnMaster < 1.0 - edgeEndCutoffTol_
617 // Master is cut, check the slave
618 point masterCutPoint =
619 masterLocalPoints[cme.start()]
620 + cutOnMaster*cme.vec(masterLocalPoints);
623 curSlaveLine.nearestDist(masterCutPoint);
627 // Strict checking of slave cut to avoid capturing
633 - curSlaveLine.start()
634 ) & curSlaveLine.vec()
635 )/sqr(curSlaveLineMag);
637 // Calculate merge tolerance from the
638 // target edge length
640 edgeCoPlanarTol_*mag(b - a);
641 // 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;
644 cutOnSlave > edgeEndCutoffTol_
645 && cutOnSlave < 1.0 - edgeEndCutoffTol_
646 && slaveCut.distance() < mergeTol
649 // Cut both master and slave. Add point
650 // to edge points The point is nominally
651 // added from the start of the master edge
652 // and added to the cut point zone
658 masterCutPoint, // point
659 masterMeshPoints[cme.start()],// m p
660 cutPointZoneID_.index(), // zone
664 // Pout << "Inserting point: " << newPoint << " as edge to edge intersection. Slave edge: " << edgeI << " " << curEdge << " master edge: " << cmeIndex << " " << cme << endl;
665 pointsIntoSlaveEdges[edgeI].append(newPoint);
666 pointsIntoMasterEdges[cmeIndex].append(newPoint);
668 // Add the point into the enriched patch map
675 // Record which two edges intersect to
679 newPoint, // Cut point index
684 masterMeshPoints[cme.start()],
685 masterMeshPoints[cme.end()]
689 slaveMeshPoints[curEdge.start()],
690 slaveMeshPoints[curEdge.end()]
697 Pout<< " " << newPoint << " = "
698 << masterCutPoint << " ";
705 // Intersection exists but it is too far
723 // Missed master edge
741 } // End if both ends missing
742 } // End for all slave edges
744 // Re-pack the points into edges lists
745 labelListList pime(pointsIntoMasterEdges.size());
747 forAll (pointsIntoMasterEdges, i)
749 pime[i].transfer(pointsIntoMasterEdges[i].shrink());
752 labelListList pise(pointsIntoSlaveEdges.size());
754 forAll (pointsIntoSlaveEdges, i)
756 pise[i].transfer(pointsIntoSlaveEdges[i].shrink());
759 // Prepare the enriched faces
760 cutPatch.calcEnrichedFaces
767 // Demand driven calculate the cut faces. Apart from the
768 // cutFaces/cutFaceMaster/cutFaceSlave no information from the cutPatch
770 const faceList& cutFaces = cutPatch.cutFaces();
771 const labelList& cutFaceMaster = cutPatch.cutFaceMaster();
772 const labelList& cutFaceSlave = cutPatch.cutFaceSlave();
774 const labelList& masterFc = masterFaceCells();
775 const labelList& slaveFc = slaveFaceCells();
777 // Couple the interface
780 // 1) Go through the cut faces and check if the cut face is the same as the
781 // defining master or slave. If so, modify the appropriate
782 // face and mark the other for relegation into the face zone.
783 // 2) If different, mark both sides for relegation and insert the new face
786 boolList orphanedMaster(masterPatch.size(), false);
787 boolList orphanedSlave(slavePatch.size(), false);
789 forAll (cutFaces, faceI)
791 const face& curCutFace = cutFaces[faceI];
792 const label curMaster = cutFaceMaster[faceI];
793 const label curSlave = cutFaceSlave[faceI];
795 // Pout << "Doing insertion of face " << faceI << ": ";
797 // Check if the face has changed topologically
798 bool insertedFace = false;
802 // Face has got a master
803 if (curCutFace == masterPatch[curMaster])
805 // Face is equal to master. Modify master face.
806 // Pout << "Face is equal to master and is ";
808 // If the face has got both master and slave, it is an
809 // internal face; otherwise it is a patch face in the
810 // master patch. Keep it in the master face zone.
814 // Pout << "internal" << endl;
815 if (masterFc[curMaster] < slaveFc[curSlave])
817 // Cut face should point into slave.
818 // Be careful about flips in zone!
823 curCutFace, // new face
824 masterPatchAddr[curMaster], // master face id
825 masterFc[curMaster], // owner
826 slaveFc[curSlave], // neighbour
829 false, // remove from zone
830 masterFaceZoneID_.index(), // zone ID
831 masterPatchFlip[curMaster] // zone flip
834 // Pout << "modifying master face. Old master: " << masterPatch[curMaster] << " new face: " << curCutFace.reverseFace() << " own: " << masterFc[curMaster] << " nei: " << slaveFc[curSlave] << endl;
838 // Cut face should point into master. Flip required.
839 // Be careful about flips in zone!
844 curCutFace.reverseFace(), // new face
845 masterPatchAddr[curMaster], // master face id
846 slaveFc[curSlave], // owner
847 masterFc[curMaster], // neighbour
850 false, // remove from zone
851 masterFaceZoneID_.index(), // zone ID
852 !masterPatchFlip[curMaster] // zone flip
858 orphanedSlave[curSlave] = true;
862 // Pout << "master boundary" << endl;
867 curCutFace, // new face
868 masterPatchAddr[curMaster], // master face index
869 masterFc[curMaster], // owner
872 masterPatchID_.index(), // patch ID
873 false, // remove from zone
874 masterFaceZoneID_.index(), // zone ID
875 masterPatchFlip[curMaster] // zone flip
883 else if (curSlave >= 0)
885 // Face has got a slave
887 // Renumber the slave face into merged labels
888 face rsf(slavePatch[curSlave]);
892 Map<label>::const_iterator mpIter = pointMergeMap.find(rsf[i]);
894 if (mpIter != pointMergeMap.end())
900 if (curCutFace == rsf)
902 // Face is equal to slave. Modify slave face.
903 // Pout << "Face is equal to slave and is ";
907 // Pout << "regular internal" << endl;
908 if (masterFc[curMaster] < slaveFc[curSlave])
914 curCutFace, // new face
915 slavePatchAddr[curSlave], // master face id
916 masterFc[curMaster], // owner
917 slaveFc[curSlave], // neighbour
920 false, // remove from zone
921 slaveFaceZoneID_.index(), // zone ID
922 !slavePatchFlip[curMaster] // zone flip
928 // Cut face should point into master.
929 // Be careful about flips in zone!
930 // Pout << "flipped internal" << endl;
935 curCutFace.reverseFace(), // new face
936 slavePatchAddr[curSlave], // master face id
937 slaveFc[curSlave], // owner
938 masterFc[curMaster], // neighbour
941 false, // remove from zone
942 slaveFaceZoneID_.index(), // zone ID
943 slavePatchFlip[curSlave] // zone flip
949 orphanedMaster[curMaster] = true;
953 // Pout << "slave boundary" << endl;
958 curCutFace, // new face
959 slavePatchAddr[curSlave], // master face index
960 slaveFc[curSlave], // owner
963 slavePatchID_.index(), // patch ID
964 false, // remove from zone
965 slaveFaceZoneID_.index(), // zone ID
966 slavePatchFlip[curSlave] // zone flip
978 "void slidingInterface::coupleInterface("
979 "polyTopoChange& ref) const"
980 ) << "Face " << faceI << " in cut faces has neither a master "
981 << "nor a slave. Error in the cutting algorithm on modify."
982 << abort(FatalError);
987 // Face is different from both master and slave
988 // Pout << "Face different from both master and slave" << endl;
995 if (masterFc[curMaster] < slaveFc[curSlave])
997 // Pout << "Adding internal face " << curCutFace << " owner: " << masterFc[curMaster] << " slave: " << slaveFc[curSlave] << " master face: " << masterPatchAddr[curMaster] << endl;
998 // Cut face should point into slave.
1003 curCutFace, // new face
1004 masterFc[curMaster], // owner
1005 slaveFc[curSlave], // neighbour
1008 masterPatchAddr[curMaster], // master face id
1011 cutFaceZoneID_.index(), // zone ID
1018 // Cut face should point into master. Flip required.
1023 curCutFace.reverseFace(), // new face
1024 slaveFc[curSlave], // owner
1025 masterFc[curMaster], // neighbour
1028 masterPatchAddr[curMaster], // master face id
1031 cutFaceZoneID_.index(), // zone ID
1038 orphanedSlave[curSlave] = true;
1042 // Pout << "Adding solo master face " << curCutFace << " owner: " << masterFc[curMaster] << " master face: " << masterPatchAddr[curMaster] << endl;
1043 // Add master patch face
1048 curCutFace, // new face
1049 masterFc[curMaster], // owner
1053 masterPatchAddr[curMaster], // master face index
1055 masterPatchID_.index(), // patch ID
1056 cutFaceZoneID_.index(), // zone ID
1063 orphanedMaster[curMaster] = true;
1065 else if (curSlave >= 0)
1067 // Pout << "Adding solo slave face " << curCutFace << " owner: " << slaveFc[curSlave] << " master face: " << slavePatchAddr[curSlave] << endl;
1068 // Add slave patch face
1073 curCutFace, // new face
1074 slaveFc[curSlave], // owner
1078 slavePatchAddr[curSlave], // master face index
1080 slavePatchID_.index(), // patch ID
1081 cutFaceZoneID_.index(), // zone ID
1087 orphanedSlave[curSlave] = true;
1093 "void slidingInterface::coupleInterface("
1094 "polyTopoChange& ref) const"
1095 ) << "Face " << faceI << " in cut faces has neither a master "
1096 << "nor a slave. Error in the cutting algorithm on add."
1097 << abort(FatalError);
1102 // Move the orphaned faces into the face zone
1103 // Pout << "Orphaned master faces: " << orphanedMaster << endl;
1104 // Pout << "Orphaned slave faces: " << orphanedSlave << endl;
1106 label nOrphanedMasters = 0;
1108 forAll (orphanedMaster, faceI)
1110 if (orphanedMaster[faceI])
1114 //// Recover original orientation
1119 // masterPatch[faceI], // new face
1120 // masterPatchAddr[faceI], // master face index
1123 // false, // flux flip
1125 // false, // remove from zone
1126 // masterFaceZoneID_.index(), // zone ID
1127 // false // zone flip
1131 //Pout<< "**MJ:deleting master face " << masterPatchAddr[faceI]
1132 // << " old verts:" << masterPatch[faceI] << endl;
1133 ref.setAction(polyRemoveFace(masterPatchAddr[faceI]));
1137 label nOrphanedSlaves = 0;
1139 forAll (orphanedSlave, faceI)
1141 if (orphanedSlave[faceI])
1145 //// Recover original orientation
1150 // slavePatch[faceI], // new face
1151 // slavePatchAddr[faceI], // slave face index
1154 // false, // flux flip
1156 // false, // remove from zone
1157 // slaveFaceZoneID_.index(), // zone ID
1158 // false // zone flip
1162 //Pout<< "**MJ:deleting slave face " << slavePatchAddr[faceI]
1163 // << " old verts:" << slavePatch[faceI] << endl;
1164 ref.setAction(polyRemoveFace(slavePatchAddr[faceI]));
1170 Pout<< "Number of orphaned faces: "
1171 << "master = " << nOrphanedMasters << " out of "
1172 << orphanedMaster.size()
1173 << " slave = " << nOrphanedSlaves << " out of "
1174 << orphanedSlave.size() << endl;
1177 // Finished coupling the plane of sliding interface.
1179 // Modify faces influenced by the sliding interface
1180 // These are the faces belonging to the master and slave cell
1181 // layer which have not already been modified.
1182 // Modification comes in three types:
1183 // 1) eliminate the points that have been removed when the old sliding
1184 // interface has been removed
1185 // 2) Merge the slave points that have hit points on the master patch
1186 // 3) Introduce new points resulting from the new sliding interface cut
1188 // Collect all affected faces
1192 // Grab the list of faces in the layer
1193 const labelList& masterStickOuts = masterStickOutFaces();
1195 // Pout << "masterStickOuts: " << masterStickOuts << endl;
1197 // Re-create the master stick-out faces
1198 forAll (masterStickOuts, faceI)
1200 // Renumber the face and remove additional points
1202 const label curFaceID = masterStickOuts[faceI];
1204 const face& oldRichFace = faces[curFaceID];
1206 bool changed = false;
1208 // Remove removed points from the face
1209 face oldFace(oldRichFace.size());
1212 forAll (oldRichFace, pointI)
1214 if (ref.pointRemoved(oldRichFace[pointI]))
1221 oldFace[nOldFace] = oldRichFace[pointI];
1226 oldFace.setSize(nOldFace);
1228 // Pout << "old rich master face: " << oldRichFace << " old face: " << oldFace << endl;
1229 DynamicList<label> newFaceLabels(2*oldFace.size());
1231 forAll (oldFace, pointI)
1233 if (masterMeshPointMap.found(oldFace[pointI]))
1235 // Point is in master patch. Add it
1237 // If the point is a direct hit, grab its label; otherwise
1238 // keep the original
1239 if (pointMergeMap.found(oldFace[pointI]))
1242 newFaceLabels.append
1244 pointMergeMap.find(oldFace[pointI])()
1249 newFaceLabels.append(oldFace[pointI]);
1252 // Find if there are additional points inserted onto the edge
1253 // between current point and the next point
1255 // 1) Find all the edges in the master patch coming
1256 // out of the current point.
1257 // 2) If the next point in the face to pick the right edge
1258 const label localFirstLabel =
1259 masterMeshPointMap.find(oldFace[pointI])();
1261 const labelList& curEdges = masterPointEdges[localFirstLabel];
1263 const label nextLabel = oldFace.nextLabel(pointI);
1265 Map<label>::const_iterator mmpmIter =
1266 masterMeshPointMap.find(nextLabel);
1268 if (mmpmIter != masterMeshPointMap.end())
1270 // Pout << "found label pair " << oldFace[pointI] << " and " << nextLabel;
1271 // Find the points on the edge between them
1272 const label localNextLabel = mmpmIter();
1274 forAll (curEdges, curEdgeI)
1278 masterEdges[curEdges[curEdgeI]].otherVertex
1285 // Pout << " found edge: " << curEdges[curEdgeI] << endl;
1287 // Get points on current edge
1288 const labelList& curPime = pime[curEdges[curEdgeI]];
1290 if (curPime.size() > 0)
1293 // Pout << "curPime: " << curPime << endl;
1294 // Insert the edge points into the face
1295 // in the correct order
1296 const point& startPoint =
1297 masterLocalPoints[localFirstLabel];
1300 masterLocalPoints[localNextLabel]
1305 scalarField edgePointWeights(curPime.size());
1307 forAll (curPime, curPimeI)
1309 edgePointWeights[curPimeI] =
1313 pointMap.find(curPime[curPimeI])()
1323 min(edgePointWeights) < 0
1324 || max(edgePointWeights) > 1
1329 "void slidingInterface::"
1331 "polyTopoChange& ref) const"
1332 ) << "Error in master stick-out edge "
1333 << "point collection."
1334 << abort(FatalError);
1338 // Go through the points and collect
1339 // them based on weights from lower to
1340 // higher. This gives the correct
1341 // order of points along the edge.
1345 passI < edgePointWeights.size();
1349 // Max weight can only be one, so
1350 // the sorting is done by
1352 label nextPoint = -1;
1355 forAll (edgePointWeights, wI)
1357 if (edgePointWeights[wI] < dist)
1359 dist = edgePointWeights[wI];
1364 // Insert the next point and reset
1365 // its weight to exclude it from
1367 newFaceLabels.append(curPime[nextPoint]);
1368 edgePointWeights[nextPoint] = GREAT;
1373 } // End of found edge
1374 } // End of looking through current edges
1379 newFaceLabels.append(oldFace[pointI]);
1385 if (newFaceLabels.size() < 3)
1389 "void slidingInterface::coupleInterface("
1390 "polyTopoChange& ref) const"
1391 ) << "Face " << curFaceID << " reduced to less than "
1392 << "3 points. Topological/cutting error A." << nl
1393 << "Old face: " << oldFace << " new face: " << newFaceLabels
1394 << abort(FatalError);
1397 // Get face zone and its flip
1398 label modifiedFaceZone = faceZones.whichZone(curFaceID);
1399 bool modifiedFaceZoneFlip = false;
1401 if (modifiedFaceZone >= 0)
1403 modifiedFaceZoneFlip =
1404 faceZones[modifiedFaceZone].flipMap()
1406 faceZones[modifiedFaceZone].whichFace(curFaceID)
1411 newFace.transfer(newFaceLabels.shrink());
1413 // Pout << "Modifying master stick-out face " << curFaceID << " old face: " << oldFace << " new face: " << newFace << endl;
1416 if (mesh.isInternalFace(curFaceID))
1422 newFace, // modified face
1423 curFaceID, // label of face being modified
1424 own[curFaceID], // owner
1425 nei[curFaceID], // neighbour
1427 -1, // patch for face
1428 false, // remove from zone
1429 modifiedFaceZone, // zone for face
1430 modifiedFaceZoneFlip // face flip in zone
1440 newFace, // modified face
1441 curFaceID, // label of face being modified
1442 own[curFaceID], // owner
1445 mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
1446 false, // remove from zone
1447 modifiedFaceZone, // zone for face
1448 modifiedFaceZoneFlip // face flip in zone
1455 // Pout << "Finished master side" << endl;
1459 // Grab the list of faces in the layer
1460 const labelList& slaveStickOuts = slaveStickOutFaces();
1462 // Pout << "slaveStickOuts: " << slaveStickOuts << endl;
1464 const Map<label>& rpm = retiredPointMap();
1466 // Re-create the slave stick-out faces
1468 forAll (slaveStickOuts, faceI)
1470 // Renumber the face and remove additional points
1471 const label curFaceID = slaveStickOuts[faceI];
1473 const face& oldRichFace = faces[curFaceID];
1475 bool changed = false;
1477 // Remove removed points from the face
1478 face oldFace(oldRichFace.size());
1481 forAll (oldRichFace, pointI)
1485 rpm.found(oldRichFace[pointI])
1486 || slaveMeshPointMap.found(oldRichFace[pointI])
1489 // Point definitely live. Add it
1490 oldFace[nOldFace] = oldRichFace[pointI];
1495 ref.pointRemoved(oldRichFace[pointI])
1496 || masterMeshPointMap.found(oldRichFace[pointI])
1499 // Point removed and not on slave patch
1500 // (first if takes care of that!) or
1501 // point belonging to master patch
1507 oldFace[nOldFace] = oldRichFace[pointI];
1512 oldFace.setSize(nOldFace);
1514 DynamicList<label> newFaceLabels(2*oldFace.size());
1516 // Pout << "old rich slave face: " << oldRichFace << " old face: " << oldFace << endl;
1517 forAll (oldFace, pointI)
1519 // Try to find the point in retired points
1520 label curP = oldFace[pointI];
1522 Map<label>::const_iterator rpmIter = rpm.find(oldFace[pointI]);
1524 if (rpmIter != rpm.end())
1530 if (slaveMeshPointMap.found(curP))
1532 // Point is in slave patch. Add it
1534 // If the point is a direct hit, grab its label; otherwise
1535 // keep the original
1536 if (pointMergeMap.found(curP))
1539 newFaceLabels.append
1541 pointMergeMap.find(curP)()
1546 newFaceLabels.append(curP);
1549 // Find if there are additional points inserted onto the edge
1550 // between current point and the next point
1552 // 1) Find all the edges in the slave patch coming
1553 // out of the current point.
1554 // 2) Use the next point in the face to pick the right edge
1556 const label localFirstLabel =
1557 slaveMeshPointMap.find(curP)();
1559 const labelList& curEdges = slavePointEdges[localFirstLabel];
1561 label nextLabel = oldFace.nextLabel(pointI);
1563 Map<label>::const_iterator rpmNextIter =
1564 rpm.find(nextLabel);
1566 if (rpmNextIter != rpm.end())
1568 nextLabel = rpmNextIter();
1571 Map<label>::const_iterator mmpmIter =
1572 slaveMeshPointMap.find(nextLabel);
1574 if (mmpmIter != slaveMeshPointMap.end())
1576 // Both points on the slave patch.
1577 // Find the points on the edge between them
1578 const label localNextLabel = mmpmIter();
1580 forAll (curEdges, curEdgeI)
1584 slaveEdges[curEdges[curEdgeI]].otherVertex
1591 // Pout << " found edge: " << curEdges[curEdgeI] << endl;
1593 // Get points on current edge
1594 const labelList& curPise = pise[curEdges[curEdgeI]];
1596 if (curPise.size() > 0)
1599 // Pout << "curPise: " << curPise << endl;
1600 // Insert the edge points into the face
1601 // in the correct order
1602 const point& startPoint =
1603 projectedSlavePoints[localFirstLabel];
1606 projectedSlavePoints[localNextLabel]
1611 scalarField edgePointWeights(curPise.size());
1613 forAll (curPise, curPiseI)
1615 edgePointWeights[curPiseI] =
1619 pointMap.find(curPise[curPiseI])()
1629 min(edgePointWeights) < 0
1630 || max(edgePointWeights) > 1
1635 "void slidingInterface::"
1637 "polyTopoChange& ref) const"
1638 ) << "Error in slave stick-out edge "
1639 << "point collection."
1640 << abort(FatalError);
1644 // Go through the points and collect
1645 // them based on weights from lower to
1646 // higher. This gives the correct
1647 // order of points along the edge.
1651 passI < edgePointWeights.size();
1655 // Max weight can only be one, so
1656 // the sorting is done by
1658 label nextPoint = -1;
1661 forAll (edgePointWeights, wI)
1663 if (edgePointWeights[wI] < dist)
1665 dist = edgePointWeights[wI];
1670 // Insert the next point and reset
1671 // its weight to exclude it from
1673 newFaceLabels.append(curPise[nextPoint]);
1674 edgePointWeights[nextPoint] = GREAT;
1680 } // End of found edge
1681 } // End of looking through current edges
1685 newFaceLabels.append(oldFace[pointI]);
1691 if (newFaceLabels.size() < 3)
1695 "void slidingInterface::coupleInterface("
1696 "polyTopoChange& ref) const"
1697 ) << "Face " << curFaceID << " reduced to less than "
1698 << "3 points. Topological/cutting error B." << nl
1699 << "Old face: " << oldFace << " new face: " << newFaceLabels
1700 << abort(FatalError);
1703 // Get face zone and its flip
1704 label modifiedFaceZone = faceZones.whichZone(curFaceID);
1705 bool modifiedFaceZoneFlip = false;
1707 if (modifiedFaceZone >= 0)
1709 modifiedFaceZoneFlip =
1710 faceZones[modifiedFaceZone].flipMap()
1712 faceZones[modifiedFaceZone].whichFace(curFaceID)
1717 newFace.transfer(newFaceLabels.shrink());
1719 // Pout << "Modifying slave stick-out face " << curFaceID << " old face: " << oldFace << " new face: " << newFace << endl;
1722 if (mesh.isInternalFace(curFaceID))
1728 newFace, // modified face
1729 curFaceID, // label of face being modified
1730 own[curFaceID], // owner
1731 nei[curFaceID], // neighbour
1733 -1, // patch for face
1734 false, // remove from zone
1735 modifiedFaceZone, // zone for face
1736 modifiedFaceZoneFlip // face flip in zone
1746 newFace, // modified face
1747 curFaceID, // label of face being modified
1748 own[curFaceID], // owner
1751 mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
1752 false, // remove from zone
1753 modifiedFaceZone, // zone for face
1754 modifiedFaceZoneFlip // face flip in zone
1761 // Activate and retire slave patch points
1762 // This needs to be done last, so that the map of removed points
1763 // does not get damaged by point modifications
1765 if (!retiredPointMapPtr_)
1769 "void slidingInterface::coupleInterface("
1770 "polyTopoChange& ref) const"
1771 ) << "Retired point map pointer not set."
1772 << abort(FatalError);
1775 Map<label>& addToRpm = *retiredPointMapPtr_;
1777 // Clear the old map
1780 label nRetiredPoints = 0;
1782 forAll (slaveMeshPoints, pointI)
1784 if (pointMergeMap.found(slaveMeshPoints[pointI]))
1786 // Retire the point - only used for supporting the detached
1794 // slaveMeshPoints[pointI], // point ID
1795 // points[slaveMeshPoints[pointI]], // point
1796 // false, // remove from zone
1797 // mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
1798 // false // in a cell
1801 //Pout<< "MJ retire slave point " << slaveMeshPoints[pointI]
1802 // << " coord " << points[slaveMeshPoints[pointI]]
1808 slaveMeshPoints[pointI]
1814 pointMergeMap.find(slaveMeshPoints[pointI])(),
1815 slaveMeshPoints[pointI]
1824 slaveMeshPoints[pointI], // point ID
1825 points[slaveMeshPoints[pointI]], // point
1826 false, // remove from zone
1827 mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
1836 Pout << "Retired " << nRetiredPoints << " out of "
1837 << slaveMeshPoints.size() << " points." << endl;
1840 // Grab cut face master and slave addressing
1841 if (cutFaceMasterPtr_) deleteDemandDrivenData(cutFaceMasterPtr_);
1842 cutFaceMasterPtr_ = new labelList(cutPatch.cutFaceMaster());
1844 if (cutFaceSlavePtr_) deleteDemandDrivenData(cutFaceSlavePtr_);
1845 cutFaceSlavePtr_ = new labelList(cutPatch.cutFaceSlave());
1847 // Finished coupling
1852 Pout<< "void slidingInterface::coupleInterface("
1853 << "polyTopoChange& ref) : "
1854 << "Finished coupling sliding interface " << name() << endl;
1859 // ************************************************************************* //