1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 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 "polyAddFace.H"
39 #include "polyModifyPoint.H"
40 #include "polyModifyFace.H"
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44 const Foam::scalar Foam::slidingInterface::edgeCoPlanarTol_ = 0.8;
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 // Index of debug signs:
49 // Index of debug signs:
50 // . - loop of the edge-to-face interaction detection
51 // x - reversal of direction in edge-to-face interaction detection
52 // + - complete edge-to-face interaction detection
53 // z - incomplete edge-to-face interaction detection. This may be OK for edges
54 // crossing from one to the other side of multiply connected master patch
56 // e - adding a point on edge
57 // f - adding a point on face
58 // . - collecting edges off another face for edge-to-edge cut
59 // + - finished collection of edges
60 // * - cut both master and slave
61 // n - cutting new edge
62 // t - intersection exists but it is outside of tolerance
63 // x - missed slave edge in cut
64 // - - missed master edge in cut
65 // u - edge already used in cutting
67 void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
71 Pout<< "void slidingInterface::coupleInterface"
72 << "(polyTopoChange& ref) : "
73 << "Coupling sliding interface " << name() << endl;
76 const polyMesh& mesh = topoChanger().mesh();
78 const pointField& points = mesh.points();
79 const faceList& faces = mesh.faces();
81 const labelList& own = mesh.faceOwner();
82 const labelList& nei = mesh.faceNeighbour();
83 const faceZoneMesh& faceZones = mesh.faceZones();
85 const primitiveFacePatch& masterPatch =
86 faceZones[masterFaceZoneID_.index()]();
88 const labelList& masterPatchAddr = faceZones[masterFaceZoneID_.index()];
90 const boolList& masterPatchFlip =
91 faceZones[masterFaceZoneID_.index()].flipMap();
93 const primitiveFacePatch& slavePatch =
94 faceZones[slaveFaceZoneID_.index()]();
96 const labelList& slavePatchAddr = faceZones[slaveFaceZoneID_.index()];
98 const boolList& slavePatchFlip =
99 faceZones[slaveFaceZoneID_.index()].flipMap();
101 const edgeList& masterEdges = masterPatch.edges();
102 const labelListList& masterPointEdges = masterPatch.pointEdges();
103 const labelList& masterMeshPoints = masterPatch.meshPoints();
104 const pointField& masterLocalPoints = masterPatch.localPoints();
105 const labelListList& masterFaceFaces = masterPatch.faceFaces();
106 const labelListList& masterFaceEdges = masterPatch.faceEdges();
107 const Map<label>& masterMeshPointMap = masterPatch.meshPointMap();
109 const edgeList& slaveEdges = slavePatch.edges();
110 const labelListList& slavePointEdges = slavePatch.pointEdges();
111 const labelList& slaveMeshPoints = slavePatch.meshPoints();
112 const pointField& slaveLocalPoints = slavePatch.localPoints();
113 const Map<label>& slaveMeshPointMap = slavePatch.meshPointMap();
114 const vectorField& slavePointNormals = slavePatch.pointNormals();
116 // Collect projection addressing
120 slavePointPointHitsPtr_
121 && slavePointEdgeHitsPtr_
122 && slavePointFaceHitsPtr_
123 && masterPointEdgeHitsPtr_
129 "void slidingInterface::coupleInterface("
130 "polyTopoChange& ref) const"
131 ) << "Point projection addressing not available."
132 << abort(FatalError);
135 const labelList& slavePointPointHits = *slavePointPointHitsPtr_;
136 const labelList& slavePointEdgeHits = *slavePointEdgeHitsPtr_;
137 const List<objectHit>& slavePointFaceHits = *slavePointFaceHitsPtr_;
138 const labelList& masterPointEdgeHits = *masterPointEdgeHitsPtr_;
139 const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
141 // Create enriched patch
142 enrichedPatch cutPatch
151 // Get reference to list of added point for the enriched patch.
152 // This will be used for point addition
153 Map<point>& pointMap = cutPatch.pointMap();
155 // Get reference to the list of merged points
156 Map<label>& pointMergeMap = cutPatch.pointMergeMap();
158 // Create mapping for every merged point of the slave patch
159 forAll (slavePointPointHits, pointI)
161 if (slavePointPointHits[pointI] >= 0)
163 // Pout << "Inserting point merge pair: " << slaveMeshPoints[pointI] << " : " << masterMeshPoints[slavePointPointHits[pointI]] << endl;
166 slaveMeshPoints[pointI],
167 masterMeshPoints[slavePointPointHits[pointI]]
172 // Collect the list of used edges for every slave edge
174 List<labelHashSet> usedMasterEdges(slaveEdges.size());
176 // Collect of slave point hits
177 forAll (slavePointPointHits, pointI)
179 // For point hits, add all point-edges from master side to all point
180 const labelList& curSlaveEdges = slavePointEdges[pointI];
182 if (slavePointPointHits[pointI] > -1)
184 const labelList& curMasterEdges =
185 masterPointEdges[slavePointPointHits[pointI]];
187 // Mark all current master edges as used for all the current slave
189 forAll (curSlaveEdges, slaveEdgeI)
191 labelHashSet& sm = usedMasterEdges[curSlaveEdges[slaveEdgeI]];
193 forAll (curMasterEdges, masterEdgeI)
195 sm.insert(curMasterEdges[masterEdgeI]);
199 else if (slavePointEdgeHits[pointI] > -1)
201 // For edge hits, add the master edge
202 forAll (curSlaveEdges, slaveEdgeI)
204 usedMasterEdges[curSlaveEdges[slaveEdgeI]].insert
206 slavePointEdgeHits[pointI]
212 // Collect off master point hits
213 // For every master point that has hit an edge, add all edges coming from
214 // that point to the slave edge that has been hit
215 forAll (masterPointEdgeHits, masterPointI)
217 if (masterPointEdgeHits[masterPointI] > -1)
219 const labelList& curMasterEdges = masterPointEdges[masterPointI];
222 usedMasterEdges[masterPointEdgeHits[masterPointI]];
224 forAll (curMasterEdges, masterEdgeI)
226 sm.insert(curMasterEdges[masterEdgeI]);
231 // Pout << "used edges: " << endl;
232 // forAll (usedMasterEdges, edgeI)
234 // Pout << "edge: " << edgeI << " used: " << usedMasterEdges[edgeI].toc() << endl;
237 // For every master and slave edge make a list of points to be added into
239 List<DynamicList<label> > pointsIntoMasterEdges(masterEdges.size());
240 List<DynamicList<label> > pointsIntoSlaveEdges(slaveEdges.size());
242 // Add all points from the slave patch that have hit the edge
243 forAll (slavePointEdgeHits, pointI)
245 if (slavePointEdgeHits[pointI] > -1)
247 // Create a new point on the master edge
250 masterEdges[slavePointEdgeHits[pointI]].line
253 ).nearestDist(projectedSlavePoints[pointI]).hitPoint();
260 edgeCutPoint, // point
261 slaveMeshPoints[pointI], // master point
262 cutPointZoneID_.index(), // zone for point
263 true // supports a cell
266 // Pout << "Inserting merge pair off edge: " << slaveMeshPoints[pointI] << " " << newPoint << " cut point: " << edgeCutPoint << " orig: " << slaveLocalPoints[pointI] << " proj: " << projectedSlavePoints[pointI] << endl;
267 // Add the new edge point into the merge map
268 pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);
270 pointsIntoMasterEdges[slavePointEdgeHits[pointI]].append
275 // Add the point into the enriched patch map
285 // Pout<< newPoint << " = " << edgeCutPoint << endl;
295 // Add all points from the slave patch that have hit a face
296 forAll (slavePointFaceHits, pointI)
300 slavePointPointHits[pointI] < 0
301 && slavePointEdgeHits[pointI] < 0
302 && slavePointFaceHits[pointI].hit()
310 projectedSlavePoints[pointI], // point
311 slaveMeshPoints[pointI], // master point
312 cutPointZoneID_.index(), // zone for point
313 true // supports a cell
316 // Pout << "Inserting merge pair off face: " << slaveMeshPoints[pointI] << " " << newPoint << endl;
317 // Add the new edge point into the merge map
318 pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);
320 // Add the point into the enriched patch map
324 projectedSlavePoints[pointI]
329 Pout<< "f: " << newPoint << " = "
330 << projectedSlavePoints[pointI] << endl;
335 forAll (masterPointEdgeHits, pointI)
337 if (masterPointEdgeHits[pointI] > -1)
339 pointsIntoSlaveEdges[masterPointEdgeHits[pointI]].append
341 masterMeshPoints[pointI]
346 // Cut all slave edges.
347 // Collect all master edges the slave edge interacts with. Skip
348 // all the edges already marked as used. For every unused edge,
349 // calculate the cut and insert the new point into the master and
351 // For the edge selection algorithm, see, comment in
352 // slidingInterfaceProjectPoints.C.
353 // Edge cutting algoritm:
354 // As the master patch defines the cutting surface, the newly
355 // inserted point needs to be on the master edge. Also, in 3-D
356 // the pair of edges generally misses each other rather than
357 // intersect. Therefore, the intersection is calculated using the
358 // plane the slave edge defines during projection. The plane is
359 // defined by the centrepoint of the edge in the original
360 // configuration and the projected end points. In case of flat
361 // geometries (when the three points are colinear), the plane is
362 // defined by the two projected end-points and the slave edge
363 // normal used as the in-plane vector. When the intersection
364 // point is created, it is added as a new point for both the
365 // master and the slave edge.
366 // In order to be able to re-create the points on edges in mesh
367 // motion without the topology change, the edge pair used to
368 // create the cut point needs to be recorded in
369 // cutPointEdgePairMap
371 // Note. "Processing slave edges" code is repeated twice in the
372 // sliding intergace coupling in order to allow the point
373 // projection to be done separately from the actual cutting.
374 // Please change consistently with slidingInterfaceProjectPoints.C
378 Pout << "Processing slave edges " << endl;
381 if (!cutPointEdgePairMapPtr_)
385 "void slidingInterface::coupleInterface("
386 "polyTopoChange& ref) const"
387 ) << "Cut point edge pair map pointer not set."
388 << abort(FatalError);
391 Map<Pair<edge> >& addToCpepm = *cutPointEdgePairMapPtr_;
396 // Create a map of faces the edge can interact with
397 labelHashSet curFaceMap
399 nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
402 labelHashSet addedFaces(2*primitiveMesh::edgesPerFace_);
404 forAll (slaveEdges, edgeI)
406 const edge& curEdge = slaveEdges[edgeI];
410 slavePointFaceHits[curEdge.start()].hit()
411 || slavePointFaceHits[curEdge.end()].hit()
414 labelHashSet& curUme = usedMasterEdges[edgeI];
415 // Pout<< "Doing edge " << edgeI << " curEdge: " << curEdge << " curUme: " << curUme << endl;
420 // Grab the faces for start and end points.
421 const label startFace =
422 slavePointFaceHits[curEdge.start()].hitObject();
423 const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
424 // Pout << "startFace: " << slavePointFaceHits[curEdge.start()] << " endFace: " << slavePointFaceHits[curEdge.end()] << endl;
425 // Insert the start face into the list
426 curFaceMap.insert(startFace);
427 addedFaces.insert(startFace);
428 // Pout << "curFaceMap: " << curFaceMap.toc() << endl;
430 bool completed = false;
432 while (nSweeps < edgeFaceEscapeLimit_)
436 if (addedFaces.found(endFace))
441 // Add all face neighbours of face in the map
442 const labelList cf = addedFaces.toc();
447 const labelList& curNbrs = masterFaceFaces[cf[cfI]];
449 forAll (curNbrs, nbrI)
451 if (!curFaceMap.found(curNbrs[nbrI]))
453 curFaceMap.insert(curNbrs[nbrI]);
454 addedFaces.insert(curNbrs[nbrI]);
459 if (completed) break;
474 // It is impossible to reach the end from the start, probably
475 // due to disconnected domain. Do search in opposite direction
477 label nReverseSweeps = 0;
480 addedFaces.insert(endFace);
482 while (nReverseSweeps < edgeFaceEscapeLimit_)
486 if (addedFaces.found(startFace))
491 // Add all face neighbours of face in the map
492 const labelList cf = addedFaces.toc();
497 const labelList& curNbrs = masterFaceFaces[cf[cfI]];
499 forAll (curNbrs, nbrI)
501 if (!curFaceMap.found(curNbrs[nbrI]))
503 curFaceMap.insert(curNbrs[nbrI]);
504 addedFaces.insert(curNbrs[nbrI]);
509 if (completed) break;
535 // Create a map of edges the edge can interact with
536 labelHashSet curMasterEdgesMap
538 nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
541 const labelList curFaces = curFaceMap.toc();
542 // Pout << "curFaces: " << curFaces << endl;
543 forAll (curFaces, faceI)
545 // Pout<< "face: " << curFaces[faceI] << " "
546 // << masterPatch[curFaces[faceI]]
548 // << masterPatch.localFaces()[curFaces[faceI]]
550 const labelList& me = masterFaceEdges[curFaces[faceI]];
554 curMasterEdgesMap.insert(me[meI]);
558 const labelList curMasterEdges = curMasterEdgesMap.toc();
560 // For all master edges to intersect, skip the ones
561 // already used and cut the rest with a cutting plane. If
562 // the intersection point, falls inside of both edges, it
566 // The edge cutting code is repeated in
567 // slidingInterface::modifyMotionPoints. This is done for
568 // efficiency reasons and avoids multiple creation of cutting
569 // planes. Please update both simultaneously.
571 const point& a = projectedSlavePoints[curEdge.start()];
572 const point& b = projectedSlavePoints[curEdge.end()];
577 slaveLocalPoints[curEdge.start()]
578 + slavePointNormals[curEdge.start()] // Add start normal
579 + slaveLocalPoints[curEdge.end()]
580 + slavePointNormals[curEdge.end()] // Add end normal
584 plane cutPlane(a, b, c);
585 // Pout << "a: " << a << " b: " << b << " c: " << c << " plane: " << cutPlane << endl;
587 linePointRef curSlaveLine = curEdge.line(projectedSlavePoints);
588 const scalar curSlaveLineMag = curSlaveLine.mag();
589 // Pout << "curSlaveLine: " << curSlaveLine << endl;
590 forAll (curMasterEdges, masterEdgeI)
592 if (!curUme.found(curMasterEdges[masterEdgeI]))
600 const label cmeIndex = curMasterEdges[masterEdgeI];
601 const edge& cme = masterEdges[cmeIndex];
602 // Pout<< "Edge " << cmeIndex << " cme: " << cme << " line: " << cme.line(masterLocalPoints) << endl;
604 cutPlane.lineIntersect
606 cme.line(masterLocalPoints)
611 cutOnMaster > edgeEndCutoffTol_
612 && cutOnMaster < 1.0 - edgeEndCutoffTol_
615 // Master is cut, check the slave
616 point masterCutPoint =
617 masterLocalPoints[cme.start()]
618 + cutOnMaster*cme.vec(masterLocalPoints);
621 curSlaveLine.nearestDist(masterCutPoint);
625 // Strict checking of slave cut to avoid capturing
631 - curSlaveLine.start()
632 ) & curSlaveLine.vec()
633 )/sqr(curSlaveLineMag);
635 // Calculate merge tolerance from the
636 // target edge length
638 edgeCoPlanarTol_*mag(b - a);
639 // 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;
642 cutOnSlave > edgeEndCutoffTol_
643 && cutOnSlave < 1.0 - edgeEndCutoffTol_
644 && slaveCut.distance() < mergeTol
647 // Cut both master and slave. Add point
648 // to edge points The point is nominally
649 // added from the start of the master edge
650 // and added to the cut point zone
656 masterCutPoint, // point
657 masterMeshPoints[cme.start()],// m p
658 cutPointZoneID_.index(), // zone
662 // Pout << "Inserting point: " << newPoint << " as edge to edge intersection. Slave edge: " << edgeI << " " << curEdge << " master edge: " << cmeIndex << " " << cme << endl;
663 pointsIntoSlaveEdges[edgeI].append(newPoint);
664 pointsIntoMasterEdges[cmeIndex].append(newPoint);
666 // Add the point into the enriched patch map
673 // Record which two edges intersect to
677 newPoint, // Cut point index
682 masterMeshPoints[cme.start()],
683 masterMeshPoints[cme.end()]
687 slaveMeshPoints[curEdge.start()],
688 slaveMeshPoints[curEdge.end()]
695 Pout<< " " << newPoint << " = "
696 << masterCutPoint << " ";
703 // Intersection exists but it is too far
721 // Missed master edge
739 } // End if both ends missing
740 } // End for all slave edges
742 // Pout << "pointsIntoMasterEdges: " << pointsIntoMasterEdges << endl;
743 // Pout << "pointsIntoSlaveEdges: " << pointsIntoSlaveEdges << endl;
745 // Re-pack the points into edges lists
746 labelListList pime(pointsIntoMasterEdges.size());
748 forAll (pointsIntoMasterEdges, i)
750 pime[i].transfer(pointsIntoMasterEdges[i].shrink());
753 labelListList pise(pointsIntoSlaveEdges.size());
755 forAll (pointsIntoSlaveEdges, i)
757 pise[i].transfer(pointsIntoSlaveEdges[i].shrink());
760 // Prepare the enriched faces
761 cutPatch.calcEnrichedFaces
768 const faceList& cutFaces = cutPatch.cutFaces();
769 const labelList& cutFaceMaster = cutPatch.cutFaceMaster();
770 const labelList& cutFaceSlave = cutPatch.cutFaceSlave();
772 const labelList& masterFc = masterFaceCells();
773 const labelList& slaveFc = slaveFaceCells();
775 // Couple the interface
778 // 1) Go through the cut faces and check if the cut face is the same as the
779 // defining master or slave. If so, modify the appropriate
780 // face and mark the other for relegation into the face zone.
781 // 2) If different, mark both sides for relegation and insert the new face
784 boolList orphanedMaster(masterPatch.size(), false);
785 boolList orphanedSlave(slavePatch.size(), false);
787 forAll (cutFaces, faceI)
789 const face& curCutFace = cutFaces[faceI];
790 const label curMaster = cutFaceMaster[faceI];
791 const label curSlave = cutFaceSlave[faceI];
793 // Pout << "Doing insertion of face " << faceI << ": ";
795 // Check if the face has changed topologically
796 bool insertedFace = false;
800 // Face has got a master
801 if (curCutFace == masterPatch[curMaster])
803 // Face is equal to master. Modify master face.
804 // Pout << "Face is equal to master and is ";
806 // If the face has got both master and slave, it is an
807 // internal face; otherwise it is a patch face in the
808 // master patch. Keep it in the master face zone.
812 // Pout << "internal" << endl;
813 if (masterFc[curMaster] < slaveFc[curSlave])
815 // Cut face should point into slave.
816 // Be careful about flips in zone!
821 curCutFace, // new face
822 masterPatchAddr[curMaster], // master face id
823 masterFc[curMaster], // owner
824 slaveFc[curSlave], // neighbour
827 false, // remove from zone
828 masterFaceZoneID_.index(), // zone ID
829 masterPatchFlip[curMaster] // zone flip
832 // Pout << "modifying master face. Old master: " << masterPatch[curMaster] << " new face: " << curCutFace.reverseFace() << " own: " << masterFc[curMaster] << " nei: " << slaveFc[curSlave] << endl;
836 // Cut face should point into master. Flip required.
837 // Be careful about flips in zone!
842 curCutFace.reverseFace(), // new face
843 masterPatchAddr[curMaster], // master face id
844 slaveFc[curSlave], // owner
845 masterFc[curMaster], // neighbour
848 false, // remove from zone
849 masterFaceZoneID_.index(), // zone ID
850 !masterPatchFlip[curMaster] // zone flip
856 orphanedSlave[curSlave] = true;
860 // Pout << "master boundary" << endl;
865 curCutFace, // new face
866 masterPatchAddr[curMaster], // master face index
867 masterFc[curMaster], // owner
870 masterPatchID_.index(), // patch ID
871 false, // remove from zone
872 masterFaceZoneID_.index(), // zone ID
873 masterPatchFlip[curMaster] // zone flip
881 else if (curSlave >= 0)
883 // Face has got a slave
885 // Renumber the slave face into merged labels
886 face rsf(slavePatch[curSlave]);
890 Map<label>::const_iterator mpIter = pointMergeMap.find(rsf[i]);
892 if (mpIter != pointMergeMap.end())
898 if (curCutFace == rsf)
900 // Face is equal to slave. Modify slave face.
901 // Pout << "Face is equal to slave and is ";
905 // Pout << "regular internal" << endl;
906 if (masterFc[curMaster] < slaveFc[curSlave])
912 curCutFace, // new face
913 slavePatchAddr[curSlave], // master face id
914 masterFc[curMaster], // owner
915 slaveFc[curSlave], // neighbour
918 false, // remove from zone
919 slaveFaceZoneID_.index(), // zone ID
920 !slavePatchFlip[curMaster] // zone flip
926 // Cut face should point into master.
927 // Be careful about flips in zone!
928 // Pout << "flipped internal" << endl;
933 curCutFace.reverseFace(), // new face
934 slavePatchAddr[curSlave], // master face id
935 slaveFc[curSlave], // owner
936 masterFc[curMaster], // neighbour
939 false, // remove from zone
940 slaveFaceZoneID_.index(), // zone ID
941 slavePatchFlip[curSlave] // zone flip
947 orphanedMaster[curMaster] = true;
951 // Pout << "slave boundary" << endl;
956 curCutFace, // new face
957 slavePatchAddr[curSlave], // master face index
958 slaveFc[curSlave], // owner
961 slavePatchID_.index(), // patch ID
962 false, // remove from zone
963 slaveFaceZoneID_.index(), // zone ID
964 slavePatchFlip[curSlave] // zone flip
976 "void slidingInterface::coupleInterface("
977 "polyTopoChange& ref) const"
978 ) << "Face " << faceI << " in cut faces has neither a master "
979 << "nor a slave. Error in the cutting algorithm on modify."
980 << abort(FatalError);
985 // Face is different from both master and slave
986 // Pout << "Face different from both master and slave" << endl;
993 if (masterFc[curMaster] < slaveFc[curSlave])
995 // Pout << "Adding internal face " << curCutFace << " owner: " << masterFc[curMaster] << " slave: " << slaveFc[curSlave] << " master face: " << masterPatchAddr[curMaster] << endl;
996 // Cut face should point into slave.
1001 curCutFace, // new face
1002 masterFc[curMaster], // owner
1003 slaveFc[curSlave], // neighbour
1006 masterPatchAddr[curMaster], // master face id
1009 cutFaceZoneID_.index(), // zone ID
1016 // Cut face should point into master. Flip required.
1021 curCutFace.reverseFace(), // new face
1022 slaveFc[curSlave], // owner
1023 masterFc[curMaster], // neighbour
1026 masterPatchAddr[curMaster], // master face id
1029 cutFaceZoneID_.index(), // zone ID
1036 orphanedSlave[curSlave] = true;
1040 // Pout << "Adding solo master face " << curCutFace << " owner: " << masterFc[curMaster] << " master face: " << masterPatchAddr[curMaster] << endl;
1041 // Add master patch face
1046 curCutFace, // new face
1047 masterFc[curMaster], // owner
1051 masterPatchAddr[curMaster], // master face index
1053 masterPatchID_.index(), // patch ID
1054 cutFaceZoneID_.index(), // zone ID
1061 orphanedMaster[curMaster] = true;
1063 else if (curSlave >= 0)
1065 // Pout << "Adding solo slave face " << curCutFace << " owner: " << slaveFc[curSlave] << " master face: " << slavePatchAddr[curSlave] << endl;
1066 // Add slave patch face
1071 curCutFace, // new face
1072 slaveFc[curSlave], // owner
1076 slavePatchAddr[curSlave], // master face index
1078 slavePatchID_.index(), // patch ID
1079 cutFaceZoneID_.index(), // zone ID
1085 orphanedSlave[curSlave] = true;
1091 "void slidingInterface::coupleInterface("
1092 "polyTopoChange& ref) const"
1093 ) << "Face " << faceI << " in cut faces has neither a master "
1094 << "nor a slave. Error in the cutting algorithm on add."
1095 << abort(FatalError);
1100 // Move the orphaned faces into the face zone
1101 // Pout << "Orphaned master faces: " << orphanedMaster << endl;
1102 // Pout << "Orphaned slave faces: " << orphanedSlave << endl;
1104 label nOrphanedMasters = 0;
1106 forAll (orphanedMaster, faceI)
1108 if (orphanedMaster[faceI])
1112 // Recover original orientation
1117 masterPatch[faceI], // new face
1118 masterPatchAddr[faceI], // master face index
1123 false, // remove from zone
1124 masterFaceZoneID_.index(), // zone ID
1131 label nOrphanedSlaves = 0;
1133 forAll (orphanedSlave, faceI)
1135 if (orphanedSlave[faceI])
1139 // Recover original orientation
1144 slavePatch[faceI], // new face
1145 slavePatchAddr[faceI], // slave face index
1150 false, // remove from zone
1151 slaveFaceZoneID_.index(), // zone ID
1160 Pout<< "Number of orphaned faces: "
1161 << "master = " << nOrphanedMasters << " out of "
1162 << orphanedMaster.size()
1163 << " slave = " << nOrphanedSlaves << " out of "
1164 << orphanedSlave.size() << endl;
1167 // Finished coupling the plane of sliding interface.
1169 // Modify faces influenced by the sliding interface
1170 // These are the faces belonging to the master and slave cell
1171 // layer which have not already been modified.
1172 // Modification comes in three types:
1173 // 1) eliminate the points that have been removed when the old sliding
1174 // interface has been removed
1175 // 2) Merge the slave points that have hit points on the master patch
1176 // 3) Introduce new points resulting from the new sliding interface cut
1178 // Collect all affected faces
1182 // Grab the list of faces in the layer
1183 const labelList& masterStickOuts = masterStickOutFaces();
1185 // Pout << "masterStickOuts: " << masterStickOuts << endl;
1187 // Re-create the master stick-out faces
1188 forAll (masterStickOuts, faceI)
1190 // Renumber the face and remove additional points
1192 const label curFaceID = masterStickOuts[faceI];
1194 const face& oldRichFace = faces[curFaceID];
1196 bool changed = false;
1198 // Remove removed points from the face
1199 face oldFace(oldRichFace.size());
1202 forAll (oldRichFace, pointI)
1204 if (ref.pointRemoved(oldRichFace[pointI]))
1211 oldFace[nOldFace] = oldRichFace[pointI];
1216 oldFace.setSize(nOldFace);
1218 // Pout << "old rich master face: " << oldRichFace << " old face: " << oldFace << endl;
1219 DynamicList<label> newFaceLabels(2*oldFace.size());
1221 forAll (oldFace, pointI)
1223 if (masterMeshPointMap.found(oldFace[pointI]))
1225 // Point is in master patch. Add it
1227 // If the point is a direct hit, grab its label; otherwise
1228 // keep the original
1229 if (pointMergeMap.found(oldFace[pointI]))
1232 newFaceLabels.append
1234 pointMergeMap.find(oldFace[pointI])()
1239 newFaceLabels.append(oldFace[pointI]);
1242 // Find if there are additional points inserted onto the edge
1243 // between current point and the next point
1245 // 1) Find all the edges in the master patch coming
1246 // out of the current point.
1247 // 2) If the next point in the face to pick the right edge
1248 const label localFirstLabel =
1249 masterMeshPointMap.find(oldFace[pointI])();
1251 const labelList& curEdges = masterPointEdges[localFirstLabel];
1253 const label nextLabel = oldFace.nextLabel(pointI);
1255 Map<label>::const_iterator mmpmIter =
1256 masterMeshPointMap.find(nextLabel);
1258 if (mmpmIter != masterMeshPointMap.end())
1260 // Pout << "found label pair " << oldFace[pointI] << " and " << nextLabel;
1261 // Find the points on the edge between them
1262 const label localNextLabel = mmpmIter();
1264 forAll (curEdges, curEdgeI)
1268 masterEdges[curEdges[curEdgeI]].otherVertex
1275 // Pout << " found edge: " << curEdges[curEdgeI] << endl;
1277 // Get points on current edge
1278 const labelList& curPime = pime[curEdges[curEdgeI]];
1280 if (curPime.size() > 0)
1283 // Pout << "curPime: " << curPime << endl;
1284 // Insert the edge points into the face
1285 // in the correct order
1286 const point& startPoint =
1287 masterLocalPoints[localFirstLabel];
1290 masterLocalPoints[localNextLabel]
1295 scalarField edgePointWeights(curPime.size());
1297 forAll (curPime, curPimeI)
1299 edgePointWeights[curPimeI] =
1303 pointMap.find(curPime[curPimeI])()
1313 min(edgePointWeights) < 0
1314 || max(edgePointWeights) > 1
1319 "void slidingInterface::"
1321 "polyTopoChange& ref) const"
1322 ) << "Error in master stick-out edge "
1323 << "point collection."
1324 << abort(FatalError);
1328 // Go through the points and collect
1329 // them based on weights from lower to
1330 // higher. This gives the correct
1331 // order of points along the edge.
1335 passI < edgePointWeights.size();
1339 // Max weight can only be one, so
1340 // the sorting is done by
1342 label nextPoint = -1;
1345 forAll (edgePointWeights, wI)
1347 if (edgePointWeights[wI] < dist)
1349 dist = edgePointWeights[wI];
1354 // Insert the next point and reset
1355 // its weight to exclude it from
1357 newFaceLabels.append(curPime[nextPoint]);
1358 edgePointWeights[nextPoint] = GREAT;
1363 } // End of found edge
1364 } // End of looking through current edges
1369 newFaceLabels.append(oldFace[pointI]);
1375 if (newFaceLabels.size() < 3)
1379 "void slidingInterface::coupleInterface("
1380 "polyTopoChange& ref) const"
1381 ) << "Face " << curFaceID << " reduced to less than "
1382 << "3 points. Topological/cutting error A." << nl
1383 << "Old face: " << oldFace << " new face: " << newFaceLabels
1384 << abort(FatalError);
1387 // Get face zone and its flip
1388 label modifiedFaceZone = faceZones.whichZone(curFaceID);
1389 bool modifiedFaceZoneFlip = false;
1391 if (modifiedFaceZone >= 0)
1393 modifiedFaceZoneFlip =
1394 faceZones[modifiedFaceZone].flipMap()
1396 faceZones[modifiedFaceZone].whichFace(curFaceID)
1401 newFace.transfer(newFaceLabels.shrink());
1403 // Pout << "Modifying master stick-out face " << curFaceID << " old face: " << oldFace << " new face: " << newFace << endl;
1410 newFace, // modified face
1411 curFaceID, // label of face being modified
1412 own[curFaceID], // owner
1413 nei[curFaceID], // neighbour
1415 mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
1416 false, // remove from zone
1417 modifiedFaceZone, // zone for face
1418 modifiedFaceZoneFlip // face flip in zone
1424 // Pout << "Finished master side" << endl;
1428 // Grab the list of faces in the layer
1429 const labelList& slaveStickOuts = slaveStickOutFaces();
1431 // Pout << "slaveStickOuts: " << slaveStickOuts << endl;
1433 const Map<label>& rpm = retiredPointMap();
1435 // Re-create the slave stick-out faces
1437 forAll (slaveStickOuts, faceI)
1439 // Renumber the face and remove additional points
1440 const label curFaceID = slaveStickOuts[faceI];
1442 const face& oldRichFace = faces[curFaceID];
1444 bool changed = false;
1446 // Remove removed points from the face
1447 face oldFace(oldRichFace.size());
1450 forAll (oldRichFace, pointI)
1454 rpm.found(oldRichFace[pointI])
1455 || slaveMeshPointMap.found(oldRichFace[pointI])
1458 // Point definitely live. Add it
1459 oldFace[nOldFace] = oldRichFace[pointI];
1464 ref.pointRemoved(oldRichFace[pointI])
1465 || masterMeshPointMap.found(oldRichFace[pointI])
1468 // Point removed and not on slave patch
1469 // (first if takes care of that!) or
1470 // point belonging to master patch
1476 oldFace[nOldFace] = oldRichFace[pointI];
1481 oldFace.setSize(nOldFace);
1483 DynamicList<label> newFaceLabels(2*oldFace.size());
1485 // Pout << "old rich slave face: " << oldRichFace << " old face: " << oldFace << endl;
1486 forAll (oldFace, pointI)
1488 // Try to find the point in retired points
1489 label curP = oldFace[pointI];
1491 Map<label>::const_iterator rpmIter = rpm.find(oldFace[pointI]);
1493 if (rpmIter != rpm.end())
1499 if (slaveMeshPointMap.found(curP))
1501 // Point is in slave patch. Add it
1503 // If the point is a direct hit, grab its label; otherwise
1504 // keep the original
1505 if (pointMergeMap.found(curP))
1508 newFaceLabels.append
1510 pointMergeMap.find(curP)()
1515 newFaceLabels.append(curP);
1518 // Find if there are additional points inserted onto the edge
1519 // between current point and the next point
1521 // 1) Find all the edges in the slave patch coming
1522 // out of the current point.
1523 // 2) Use the next point in the face to pick the right edge
1525 const label localFirstLabel =
1526 slaveMeshPointMap.find(curP)();
1528 const labelList& curEdges = slavePointEdges[localFirstLabel];
1530 label nextLabel = oldFace.nextLabel(pointI);
1532 Map<label>::const_iterator rpmNextIter =
1533 rpm.find(nextLabel);
1535 if (rpmNextIter != rpm.end())
1537 nextLabel = rpmNextIter();
1540 Map<label>::const_iterator mmpmIter =
1541 slaveMeshPointMap.find(nextLabel);
1543 if (mmpmIter != slaveMeshPointMap.end())
1545 // Both points on the slave patch.
1546 // Find the points on the edge between them
1547 const label localNextLabel = mmpmIter();
1549 forAll (curEdges, curEdgeI)
1553 slaveEdges[curEdges[curEdgeI]].otherVertex
1560 // Pout << " found edge: " << curEdges[curEdgeI] << endl;
1562 // Get points on current edge
1563 const labelList& curPise = pise[curEdges[curEdgeI]];
1565 if (curPise.size() > 0)
1568 // Pout << "curPise: " << curPise << endl;
1569 // Insert the edge points into the face
1570 // in the correct order
1571 const point& startPoint =
1572 projectedSlavePoints[localFirstLabel];
1575 projectedSlavePoints[localNextLabel]
1580 scalarField edgePointWeights(curPise.size());
1582 forAll (curPise, curPiseI)
1584 edgePointWeights[curPiseI] =
1588 pointMap.find(curPise[curPiseI])()
1598 min(edgePointWeights) < 0
1599 || max(edgePointWeights) > 1
1604 "void slidingInterface::"
1606 "polyTopoChange& ref) const"
1607 ) << "Error in slave stick-out edge "
1608 << "point collection."
1609 << abort(FatalError);
1613 // Go through the points and collect
1614 // them based on weights from lower to
1615 // higher. This gives the correct
1616 // order of points along the edge.
1620 passI < edgePointWeights.size();
1624 // Max weight can only be one, so
1625 // the sorting is done by
1627 label nextPoint = -1;
1630 forAll (edgePointWeights, wI)
1632 if (edgePointWeights[wI] < dist)
1634 dist = edgePointWeights[wI];
1639 // Insert the next point and reset
1640 // its weight to exclude it from
1642 newFaceLabels.append(curPise[nextPoint]);
1643 edgePointWeights[nextPoint] = GREAT;
1649 } // End of found edge
1650 } // End of looking through current edges
1654 newFaceLabels.append(oldFace[pointI]);
1660 if (newFaceLabels.size() < 3)
1664 "void slidingInterface::coupleInterface("
1665 "polyTopoChange& ref) const"
1666 ) << "Face " << curFaceID << " reduced to less than "
1667 << "3 points. Topological/cutting error B." << nl
1668 << "Old face: " << oldFace << " new face: " << newFaceLabels
1669 << abort(FatalError);
1672 // Get face zone and its flip
1673 label modifiedFaceZone = faceZones.whichZone(curFaceID);
1674 bool modifiedFaceZoneFlip = false;
1676 if (modifiedFaceZone >= 0)
1678 modifiedFaceZoneFlip =
1679 faceZones[modifiedFaceZone].flipMap()
1681 faceZones[modifiedFaceZone].whichFace(curFaceID)
1686 newFace.transfer(newFaceLabels.shrink());
1688 // Pout << "Modifying slave stick-out face " << curFaceID << " old face: " << oldFace << " new face: " << newFace << endl;
1695 newFace, // modified face
1696 curFaceID, // label of face being modified
1697 own[curFaceID], // owner
1698 nei[curFaceID], // neighbour
1700 mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
1701 false, // remove from zone
1702 modifiedFaceZone, // zone for face
1703 modifiedFaceZoneFlip // face flip in zone
1709 // Activate and retire slave patch points
1710 // This needs to be done last, so that the map of removed points
1711 // does not get damaged by point modifications
1713 if (!retiredPointMapPtr_)
1717 "void slidingInterface::coupleInterface("
1718 "polyTopoChange& ref) const"
1719 ) << "Retired point map pointer not set."
1720 << abort(FatalError);
1723 Map<label>& addToRpm = *retiredPointMapPtr_;
1725 // Clear the old map
1728 label nRetiredPoints = 0;
1730 forAll (slaveMeshPoints, pointI)
1732 if (pointMergeMap.found(slaveMeshPoints[pointI]))
1734 // Retire the point - only used for supporting the detached
1742 slaveMeshPoints[pointI], // point ID
1743 points[slaveMeshPoints[pointI]], // point
1744 false, // remove from zone
1745 mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
1752 pointMergeMap.find(slaveMeshPoints[pointI])(),
1753 slaveMeshPoints[pointI]
1762 slaveMeshPoints[pointI], // point ID
1763 points[slaveMeshPoints[pointI]], // point
1764 false, // remove from zone
1765 mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
1774 Pout << "Retired " << nRetiredPoints << " out of "
1775 << slaveMeshPoints.size() << " points." << endl;
1778 // Grab cut face master and slave addressing
1779 if (cutFaceMasterPtr_) deleteDemandDrivenData(cutFaceMasterPtr_);
1780 cutFaceMasterPtr_ = new labelList(cutPatch.cutFaceMaster());
1782 if (cutFaceSlavePtr_) deleteDemandDrivenData(cutFaceSlavePtr_);
1783 cutFaceSlavePtr_ = new labelList(cutPatch.cutFaceSlave());
1785 // Finished coupling
1790 Pout<< "void slidingInterface::coupleInterface("
1791 << "polyTopoChange& ref) : "
1792 << "Finished coupling sliding interface " << name() << endl;
1797 // ************************************************************************* //