Files
openfoam/src/dynamicMesh/slidingInterface/coupleSlidingInterface.C
2008-06-25 15:01:46 +02:00

1798 lines
68 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "slidingInterface.H"
#include "polyTopoChange.H"
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "enrichedPatch.H"
#include "DynamicList.H"
#include "pointHit.H"
#include "triPointRef.H"
#include "plane.H"
#include "polyTopoChanger.H"
#include "polyAddPoint.H"
#include "polyAddFace.H"
#include "polyModifyPoint.H"
#include "polyModifyFace.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::scalar Foam::slidingInterface::edgeCoPlanarTol_ = 0.8;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Index of debug signs:
// Index of debug signs:
// . - loop of the edge-to-face interaction detection
// x - reversal of direction in edge-to-face interaction detection
// + - complete edge-to-face interaction detection
// z - incomplete edge-to-face interaction detection. This may be OK for edges
// crossing from one to the other side of multiply connected master patch
// e - adding a point on edge
// f - adding a point on face
// . - collecting edges off another face for edge-to-edge cut
// + - finished collection of edges
// * - cut both master and slave
// n - cutting new edge
// t - intersection exists but it is outside of tolerance
// x - missed slave edge in cut
// - - missed master edge in cut
// u - edge already used in cutting
void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
{
if (debug)
{
Pout<< "void slidingInterface::coupleInterface"
<< "(polyTopoChange& ref) : "
<< "Coupling sliding interface " << name() << endl;
}
const polyMesh& mesh = topoChanger().mesh();
const pointField& points = mesh.points();
const faceList& faces = mesh.faces();
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const faceZoneMesh& faceZones = mesh.faceZones();
const primitiveFacePatch& masterPatch =
faceZones[masterFaceZoneID_.index()]();
const labelList& masterPatchAddr = faceZones[masterFaceZoneID_.index()];
const boolList& masterPatchFlip =
faceZones[masterFaceZoneID_.index()].flipMap();
const primitiveFacePatch& slavePatch =
faceZones[slaveFaceZoneID_.index()]();
const labelList& slavePatchAddr = faceZones[slaveFaceZoneID_.index()];
const boolList& slavePatchFlip =
faceZones[slaveFaceZoneID_.index()].flipMap();
const edgeList& masterEdges = masterPatch.edges();
const labelListList& masterPointEdges = masterPatch.pointEdges();
const labelList& masterMeshPoints = masterPatch.meshPoints();
const pointField& masterLocalPoints = masterPatch.localPoints();
const labelListList& masterFaceFaces = masterPatch.faceFaces();
const labelListList& masterFaceEdges = masterPatch.faceEdges();
const Map<label>& masterMeshPointMap = masterPatch.meshPointMap();
const edgeList& slaveEdges = slavePatch.edges();
const labelListList& slavePointEdges = slavePatch.pointEdges();
const labelList& slaveMeshPoints = slavePatch.meshPoints();
const pointField& slaveLocalPoints = slavePatch.localPoints();
const Map<label>& slaveMeshPointMap = slavePatch.meshPointMap();
const vectorField& slavePointNormals = slavePatch.pointNormals();
// Collect projection addressing
if
(
!(
slavePointPointHitsPtr_
&& slavePointEdgeHitsPtr_
&& slavePointFaceHitsPtr_
&& masterPointEdgeHitsPtr_
)
)
{
FatalErrorIn
(
"void slidingInterface::coupleInterface("
"polyTopoChange& ref) const"
) << "Point projection addressing not available."
<< abort(FatalError);
}
const labelList& slavePointPointHits = *slavePointPointHitsPtr_;
const labelList& slavePointEdgeHits = *slavePointEdgeHitsPtr_;
const List<objectHit>& slavePointFaceHits = *slavePointFaceHitsPtr_;
const labelList& masterPointEdgeHits = *masterPointEdgeHitsPtr_;
const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
// Create enriched patch
enrichedPatch cutPatch
(
masterPatch,
slavePatch,
slavePointPointHits,
slavePointEdgeHits,
slavePointFaceHits
);
// Get reference to list of added point for the enriched patch.
// This will be used for point addition
Map<point>& pointMap = cutPatch.pointMap();
// Get reference to the list of merged points
Map<label>& pointMergeMap = cutPatch.pointMergeMap();
// Create mapping for every merged point of the slave patch
forAll (slavePointPointHits, pointI)
{
if (slavePointPointHits[pointI] >= 0)
{
// Pout << "Inserting point merge pair: " << slaveMeshPoints[pointI] << " : " << masterMeshPoints[slavePointPointHits[pointI]] << endl;
pointMergeMap.insert
(
slaveMeshPoints[pointI],
masterMeshPoints[slavePointPointHits[pointI]]
);
}
}
// Collect the list of used edges for every slave edge
List<labelHashSet> usedMasterEdges(slaveEdges.size());
// Collect of slave point hits
forAll (slavePointPointHits, pointI)
{
// For point hits, add all point-edges from master side to all point
const labelList& curSlaveEdges = slavePointEdges[pointI];
if (slavePointPointHits[pointI] > -1)
{
const labelList& curMasterEdges =
masterPointEdges[slavePointPointHits[pointI]];
// Mark all current master edges as used for all the current slave
// edges
forAll (curSlaveEdges, slaveEdgeI)
{
labelHashSet& sm = usedMasterEdges[curSlaveEdges[slaveEdgeI]];
forAll (curMasterEdges, masterEdgeI)
{
sm.insert(curMasterEdges[masterEdgeI]);
}
}
}
else if (slavePointEdgeHits[pointI] > -1)
{
// For edge hits, add the master edge
forAll (curSlaveEdges, slaveEdgeI)
{
usedMasterEdges[curSlaveEdges[slaveEdgeI]].insert
(
slavePointEdgeHits[pointI]
);
}
}
}
// Collect off master point hits
// For every master point that has hit an edge, add all edges coming from
// that point to the slave edge that has been hit
forAll (masterPointEdgeHits, masterPointI)
{
if (masterPointEdgeHits[masterPointI] > -1)
{
const labelList& curMasterEdges = masterPointEdges[masterPointI];
labelHashSet& sm =
usedMasterEdges[masterPointEdgeHits[masterPointI]];
forAll (curMasterEdges, masterEdgeI)
{
sm.insert(curMasterEdges[masterEdgeI]);
}
}
}
// Pout << "used edges: " << endl;
// forAll (usedMasterEdges, edgeI)
// {
// Pout << "edge: " << edgeI << " used: " << usedMasterEdges[edgeI].toc() << endl;
// }
// For every master and slave edge make a list of points to be added into
// that edge.
List<DynamicList<label> > pointsIntoMasterEdges(masterEdges.size());
List<DynamicList<label> > pointsIntoSlaveEdges(slaveEdges.size());
// Add all points from the slave patch that have hit the edge
forAll (slavePointEdgeHits, pointI)
{
if (slavePointEdgeHits[pointI] > -1)
{
// Create a new point on the master edge
point edgeCutPoint =
masterEdges[slavePointEdgeHits[pointI]].line
(
masterLocalPoints
).nearestDist(projectedSlavePoints[pointI]).hitPoint();
label newPoint =
ref.setAction
(
polyAddPoint
(
edgeCutPoint, // point
slaveMeshPoints[pointI], // master point
cutPointZoneID_.index(), // zone for point
true // supports a cell
)
);
// Pout << "Inserting merge pair off edge: " << slaveMeshPoints[pointI] << " " << newPoint << " cut point: " << edgeCutPoint << " orig: " << slaveLocalPoints[pointI] << " proj: " << projectedSlavePoints[pointI] << endl;
// Add the new edge point into the merge map
pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);
pointsIntoMasterEdges[slavePointEdgeHits[pointI]].append
(
newPoint
);
// Add the point into the enriched patch map
pointMap.insert
(
newPoint,
edgeCutPoint
);
if (debug)
{
Pout<< "e";
// Pout<< newPoint << " = " << edgeCutPoint << endl;
}
}
}
if (debug)
{
Pout<< endl;
}
// Add all points from the slave patch that have hit a face
forAll (slavePointFaceHits, pointI)
{
if
(
slavePointPointHits[pointI] < 0
&& slavePointEdgeHits[pointI] < 0
&& slavePointFaceHits[pointI].hit()
)
{
label newPoint =
ref.setAction
(
polyAddPoint
(
projectedSlavePoints[pointI], // point
slaveMeshPoints[pointI], // master point
cutPointZoneID_.index(), // zone for point
true // supports a cell
)
);
// Pout << "Inserting merge pair off face: " << slaveMeshPoints[pointI] << " " << newPoint << endl;
// Add the new edge point into the merge map
pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);
// Add the point into the enriched patch map
pointMap.insert
(
newPoint,
projectedSlavePoints[pointI]
);
if (debug)
{
Pout<< "f: " << newPoint << " = "
<< projectedSlavePoints[pointI] << endl;
}
}
}
forAll (masterPointEdgeHits, pointI)
{
if (masterPointEdgeHits[pointI] > -1)
{
pointsIntoSlaveEdges[masterPointEdgeHits[pointI]].append
(
masterMeshPoints[pointI]
);
}
}
// Cut all slave edges.
// Collect all master edges the slave edge interacts with. Skip
// all the edges already marked as used. For every unused edge,
// calculate the cut and insert the new point into the master and
// slave edge.
// For the edge selection algorithm, see, comment in
// slidingInterfaceProjectPoints.C.
// Edge cutting algoritm:
// As the master patch defines the cutting surface, the newly
// inserted point needs to be on the master edge. Also, in 3-D
// the pair of edges generally misses each other rather than
// intersect. Therefore, the intersection is calculated using the
// plane the slave edge defines during projection. The plane is
// defined by the centrepoint of the edge in the original
// configuration and the projected end points. In case of flat
// geometries (when the three points are colinear), the plane is
// defined by the two projected end-points and the slave edge
// normal used as the in-plane vector. When the intersection
// point is created, it is added as a new point for both the
// master and the slave edge.
// In order to be able to re-create the points on edges in mesh
// motion without the topology change, the edge pair used to
// create the cut point needs to be recorded in
// cutPointEdgePairMap
// Note. "Processing slave edges" code is repeated twice in the
// sliding intergace coupling in order to allow the point
// projection to be done separately from the actual cutting.
// Please change consistently with slidingInterfaceProjectPoints.C
//
if (debug)
{
Pout << "Processing slave edges " << endl;
}
if (!cutPointEdgePairMapPtr_)
{
FatalErrorIn
(
"void slidingInterface::coupleInterface("
"polyTopoChange& ref) const"
) << "Cut point edge pair map pointer not set."
<< abort(FatalError);
}
Map<Pair<edge> >& addToCpepm = *cutPointEdgePairMapPtr_;
// Clear the old map
addToCpepm.clear();
// Create a map of faces the edge can interact with
labelHashSet curFaceMap
(
nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
);
labelHashSet addedFaces(2*primitiveMesh::edgesPerFace_);
forAll (slaveEdges, edgeI)
{
const edge& curEdge = slaveEdges[edgeI];
if
(
slavePointFaceHits[curEdge.start()].hit()
|| slavePointFaceHits[curEdge.end()].hit()
)
{
labelHashSet& curUme = usedMasterEdges[edgeI];
// Pout<< "Doing edge " << edgeI << " curEdge: " << curEdge << " curUme: " << curUme << endl;
// Clear the maps
curFaceMap.clear();
addedFaces.clear();
// Grab the faces for start and end points.
const label startFace =
slavePointFaceHits[curEdge.start()].hitObject();
const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
// Pout << "startFace: " << slavePointFaceHits[curEdge.start()] << " endFace: " << slavePointFaceHits[curEdge.end()] << endl;
// Insert the start face into the list
curFaceMap.insert(startFace);
addedFaces.insert(startFace);
// Pout << "curFaceMap: " << curFaceMap.toc() << endl;
label nSweeps = 0;
bool completed = false;
while (nSweeps < edgeFaceEscapeLimit_)
{
nSweeps++;
if (addedFaces.found(endFace))
{
completed = true;
}
// Add all face neighbours of face in the map
const labelList cf = addedFaces.toc();
addedFaces.clear();
forAll (cf, cfI)
{
const labelList& curNbrs = masterFaceFaces[cf[cfI]];
forAll (curNbrs, nbrI)
{
if (!curFaceMap.found(curNbrs[nbrI]))
{
curFaceMap.insert(curNbrs[nbrI]);
addedFaces.insert(curNbrs[nbrI]);
}
}
}
if (completed) break;
if (debug)
{
Pout << ".";
}
}
if (!completed)
{
if (debug)
{
Pout << "x";
}
// It is impossible to reach the end from the start, probably
// due to disconnected domain. Do search in opposite direction
label nReverseSweeps = 0;
addedFaces.clear();
addedFaces.insert(endFace);
while (nReverseSweeps < edgeFaceEscapeLimit_)
{
nReverseSweeps++;
if (addedFaces.found(startFace))
{
completed = true;
}
// Add all face neighbours of face in the map
const labelList cf = addedFaces.toc();
addedFaces.clear();
forAll (cf, cfI)
{
const labelList& curNbrs = masterFaceFaces[cf[cfI]];
forAll (curNbrs, nbrI)
{
if (!curFaceMap.found(curNbrs[nbrI]))
{
curFaceMap.insert(curNbrs[nbrI]);
addedFaces.insert(curNbrs[nbrI]);
}
}
}
if (completed) break;
if (debug)
{
Pout << ".";
}
}
}
if (completed)
{
if (debug)
{
Pout << "+ ";
}
}
else
{
if (debug)
{
Pout << "z ";
}
}
// Collect the edges
// Create a map of edges the edge can interact with
labelHashSet curMasterEdgesMap
(
nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
);
const labelList curFaces = curFaceMap.toc();
// Pout << "curFaces: " << curFaces << endl;
forAll (curFaces, faceI)
{
// Pout<< "face: " << curFaces[faceI] << " "
// << masterPatch[curFaces[faceI]]
// << " local: "
// << masterPatch.localFaces()[curFaces[faceI]]
// << endl;
const labelList& me = masterFaceEdges[curFaces[faceI]];
forAll (me, meI)
{
curMasterEdgesMap.insert(me[meI]);
}
}
const labelList curMasterEdges = curMasterEdgesMap.toc();
// For all master edges to intersect, skip the ones
// already used and cut the rest with a cutting plane. If
// the intersection point, falls inside of both edges, it
// is valid.
// Note.
// The edge cutting code is repeated in
// slidingInterface::modifyMotionPoints. This is done for
// efficiency reasons and avoids multiple creation of cutting
// planes. Please update both simultaneously.
const point& a = projectedSlavePoints[curEdge.start()];
const point& b = projectedSlavePoints[curEdge.end()];
point c =
0.5*
(
slaveLocalPoints[curEdge.start()]
+ slavePointNormals[curEdge.start()] // Add start normal
+ slaveLocalPoints[curEdge.end()]
+ slavePointNormals[curEdge.end()] // Add end normal
);
// Create the plane
plane cutPlane(a, b, c);
// Pout << "a: " << a << " b: " << b << " c: " << c << " plane: " << cutPlane << endl;
linePointRef curSlaveLine = curEdge.line(projectedSlavePoints);
const scalar curSlaveLineMag = curSlaveLine.mag();
// Pout << "curSlaveLine: " << curSlaveLine << endl;
forAll (curMasterEdges, masterEdgeI)
{
if (!curUme.found(curMasterEdges[masterEdgeI]))
{
// New edge
if (debug)
{
Pout << "n";
}
const label cmeIndex = curMasterEdges[masterEdgeI];
const edge& cme = masterEdges[cmeIndex];
// Pout<< "Edge " << cmeIndex << " cme: " << cme << " line: " << cme.line(masterLocalPoints) << endl;
scalar cutOnMaster =
cutPlane.lineIntersect
(
cme.line(masterLocalPoints)
);
if
(
cutOnMaster > edgeEndCutoffTol_
&& cutOnMaster < 1.0 - edgeEndCutoffTol_
)
{
// Master is cut, check the slave
point masterCutPoint =
masterLocalPoints[cme.start()]
+ cutOnMaster*cme.vec(masterLocalPoints);
pointHit slaveCut =
curSlaveLine.nearestDist(masterCutPoint);
if (slaveCut.hit())
{
// Strict checking of slave cut to avoid capturing
// end points.
scalar cutOnSlave =
(
(
slaveCut.hitPoint()
- curSlaveLine.start()
) & curSlaveLine.vec()
)/sqr(curSlaveLineMag);
// Calculate merge tolerance from the
// target edge length
scalar mergeTol =
edgeCoPlanarTol_*mag(b - a);
// 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;
if
(
cutOnSlave > edgeEndCutoffTol_
&& cutOnSlave < 1.0 - edgeEndCutoffTol_
&& slaveCut.distance() < mergeTol
)
{
// Cut both master and slave. Add point
// to edge points The point is nominally
// added from the start of the master edge
// and added to the cut point zone
label newPoint =
ref.setAction
(
polyAddPoint
(
masterCutPoint, // point
masterMeshPoints[cme.start()],// m p
cutPointZoneID_.index(), // zone
true // active
)
);
// Pout << "Inserting point: " << newPoint << " as edge to edge intersection. Slave edge: " << edgeI << " " << curEdge << " master edge: " << cmeIndex << " " << cme << endl;
pointsIntoSlaveEdges[edgeI].append(newPoint);
pointsIntoMasterEdges[cmeIndex].append(newPoint);
// Add the point into the enriched patch map
pointMap.insert
(
newPoint,
masterCutPoint
);
// Record which two edges intersect to
// create cut point
addToCpepm.insert
(
newPoint, // Cut point index
Pair<edge>
(
edge
(
masterMeshPoints[cme.start()],
masterMeshPoints[cme.end()]
), // Master edge
edge
(
slaveMeshPoints[curEdge.start()],
slaveMeshPoints[curEdge.end()]
)// Slave edge
)
);
if (debug)
{
Pout<< " " << newPoint << " = "
<< masterCutPoint << " ";
}
}
else
{
if (debug)
{
// Intersection exists but it is too far
Pout << "t";
}
}
}
else
{
if (debug)
{
// Missed slave edge
Pout << "x";
}
}
}
else
{
if (debug)
{
// Missed master edge
Pout << "-";
}
}
}
else
{
if (debug)
{
Pout << "u";
}
}
}
if (debug)
{
Pout << endl;
}
} // End if both ends missing
} // End for all slave edges
// Pout << "pointsIntoMasterEdges: " << pointsIntoMasterEdges << endl;
// Pout << "pointsIntoSlaveEdges: " << pointsIntoSlaveEdges << endl;
// Re-pack the points into edges lists
labelListList pime(pointsIntoMasterEdges.size());
forAll (pointsIntoMasterEdges, i)
{
pime[i].transfer(pointsIntoMasterEdges[i].shrink());
}
labelListList pise(pointsIntoSlaveEdges.size());
forAll (pointsIntoSlaveEdges, i)
{
pise[i].transfer(pointsIntoSlaveEdges[i].shrink());
}
// Prepare the enriched faces
cutPatch.calcEnrichedFaces
(
pime,
pise,
projectedSlavePoints
);
const faceList& cutFaces = cutPatch.cutFaces();
const labelList& cutFaceMaster = cutPatch.cutFaceMaster();
const labelList& cutFaceSlave = cutPatch.cutFaceSlave();
const labelList& masterFc = masterFaceCells();
const labelList& slaveFc = slaveFaceCells();
// Couple the interface
// Algorithm:
// 1) Go through the cut faces and check if the cut face is the same as the
// defining master or slave. If so, modify the appropriate
// face and mark the other for relegation into the face zone.
// 2) If different, mark both sides for relegation and insert the new face
boolList orphanedMaster(masterPatch.size(), false);
boolList orphanedSlave(slavePatch.size(), false);
forAll (cutFaces, faceI)
{
const face& curCutFace = cutFaces[faceI];
const label curMaster = cutFaceMaster[faceI];
const label curSlave = cutFaceSlave[faceI];
// Pout << "Doing insertion of face " << faceI << ": ";
// Check if the face has changed topologically
bool insertedFace = false;
if (curMaster >= 0)
{
// Face has got a master
if (curCutFace == masterPatch[curMaster])
{
// Face is equal to master. Modify master face.
// Pout << "Face is equal to master and is ";
// If the face has got both master and slave, it is an
// internal face; otherwise it is a patch face in the
// master patch. Keep it in the master face zone.
if (curSlave >= 0)
{
// Pout << "internal" << endl;
if (masterFc[curMaster] < slaveFc[curSlave])
{
// Cut face should point into slave.
// Be careful about flips in zone!
ref.setAction
(
polyModifyFace
(
curCutFace, // new face
masterPatchAddr[curMaster], // master face id
masterFc[curMaster], // owner
slaveFc[curSlave], // neighbour
false, // flux flip
-1, // patch ID
false, // remove from zone
masterFaceZoneID_.index(), // zone ID
masterPatchFlip[curMaster] // zone flip
)
);
// Pout << "modifying master face. Old master: " << masterPatch[curMaster] << " new face: " << curCutFace.reverseFace() << " own: " << masterFc[curMaster] << " nei: " << slaveFc[curSlave] << endl;
}
else
{
// Cut face should point into master. Flip required.
// Be careful about flips in zone!
ref.setAction
(
polyModifyFace
(
curCutFace.reverseFace(), // new face
masterPatchAddr[curMaster], // master face id
slaveFc[curSlave], // owner
masterFc[curMaster], // neighbour
true, // flux flip
-1, // patch ID
false, // remove from zone
masterFaceZoneID_.index(), // zone ID
!masterPatchFlip[curMaster] // zone flip
)
);
}
// Orphan the slave
orphanedSlave[curSlave] = true;
}
else
{
// Pout << "master boundary" << endl;
ref.setAction
(
polyModifyFace
(
curCutFace, // new face
masterPatchAddr[curMaster], // master face index
masterFc[curMaster], // owner
-1, // neighbour
false, // flux flip
masterPatchID_.index(), // patch ID
false, // remove from zone
masterFaceZoneID_.index(), // zone ID
masterPatchFlip[curMaster] // zone flip
)
);
}
insertedFace = true;
}
}
else if (curSlave >= 0)
{
// Face has got a slave
// Renumber the slave face into merged labels
face rsf(slavePatch[curSlave]);
forAll (rsf, i)
{
Map<label>::const_iterator mpIter = pointMergeMap.find(rsf[i]);
if (mpIter != pointMergeMap.end())
{
rsf[i] = mpIter();
}
}
if (curCutFace == rsf)
{
// Face is equal to slave. Modify slave face.
// Pout << "Face is equal to slave and is ";
if (curMaster >= 0)
{
// Pout << "regular internal" << endl;
if (masterFc[curMaster] < slaveFc[curSlave])
{
ref.setAction
(
polyModifyFace
(
curCutFace, // new face
slavePatchAddr[curSlave], // master face id
masterFc[curMaster], // owner
slaveFc[curSlave], // neighbour
true, // flux flip
-1, // patch ID
false, // remove from zone
slaveFaceZoneID_.index(), // zone ID
!slavePatchFlip[curMaster] // zone flip
)
);
}
else
{
// Cut face should point into master.
// Be careful about flips in zone!
// Pout << "flipped internal" << endl;
ref.setAction
(
polyModifyFace
(
curCutFace.reverseFace(), // new face
slavePatchAddr[curSlave], // master face id
slaveFc[curSlave], // owner
masterFc[curMaster], // neighbour
false, // flux flip
-1, // patch ID
false, // remove from zone
slaveFaceZoneID_.index(), // zone ID
slavePatchFlip[curSlave] // zone flip
)
);
}
// Orphan the master
orphanedMaster[curMaster] = true;
}
else
{
// Pout << "slave boundary" << endl;
ref.setAction
(
polyModifyFace
(
curCutFace, // new face
slavePatchAddr[curSlave], // master face index
slaveFc[curSlave], // owner
-1, // neighbour
false, // flux flip
slavePatchID_.index(), // patch ID
false, // remove from zone
slaveFaceZoneID_.index(), // zone ID
slavePatchFlip[curSlave] // zone flip
)
);
}
insertedFace = true;
}
}
else
{
FatalErrorIn
(
"void slidingInterface::coupleInterface("
"polyTopoChange& ref) const"
) << "Face " << faceI << " in cut faces has neither a master "
<< "nor a slave. Error in the cutting algorithm on modify."
<< abort(FatalError);
}
if (!insertedFace)
{
// Face is different from both master and slave
// Pout << "Face different from both master and slave" << endl;
if (curMaster >= 0)
{
if (curSlave >= 0)
{
// Add internal face
if (masterFc[curMaster] < slaveFc[curSlave])
{
// Pout << "Adding internal face " << curCutFace << " owner: " << masterFc[curMaster] << " slave: " << slaveFc[curSlave] << " master face: " << masterPatchAddr[curMaster] << endl;
// Cut face should point into slave.
ref.setAction
(
polyAddFace
(
curCutFace, // new face
masterFc[curMaster], // owner
slaveFc[curSlave], // neighbour
-1, // master point
-1, // master edge
masterPatchAddr[curMaster], // master face id
false, // flux flip
-1, // patch ID
cutFaceZoneID_.index(), // zone ID
false // zone flip
)
);
}
else
{
// Cut face should point into master. Flip required.
ref.setAction
(
polyAddFace
(
curCutFace.reverseFace(), // new face
slaveFc[curSlave], // owner
masterFc[curMaster], // neighbour
-1, // master point
-1, // master edge
masterPatchAddr[curMaster], // master face id
true, // flux flip
-1, // patch ID
cutFaceZoneID_.index(), // zone ID
true // zone flip
)
);
}
// Orphan slave
orphanedSlave[curSlave] = true;
}
else
{
// Pout << "Adding solo master face " << curCutFace << " owner: " << masterFc[curMaster] << " master face: " << masterPatchAddr[curMaster] << endl;
// Add master patch face
ref.setAction
(
polyAddFace
(
curCutFace, // new face
masterFc[curMaster], // owner
-1, // neighbour
-1, // master point
-1, // master edge
masterPatchAddr[curMaster], // master face index
false, // flux flip
masterPatchID_.index(), // patch ID
cutFaceZoneID_.index(), // zone ID
false // zone flip
)
);
}
// Orphan master
orphanedMaster[curMaster] = true;
}
else if (curSlave >= 0)
{
// Pout << "Adding solo slave face " << curCutFace << " owner: " << slaveFc[curSlave] << " master face: " << slavePatchAddr[curSlave] << endl;
// Add slave patch face
ref.setAction
(
polyAddFace
(
curCutFace, // new face
slaveFc[curSlave], // owner
-1, // neighbour
-1, // master point
-1, // master edge
slavePatchAddr[curSlave], // master face index
false, // flux flip
slavePatchID_.index(), // patch ID
cutFaceZoneID_.index(), // zone ID
false // zone flip
)
);
// Orphan slave
orphanedSlave[curSlave] = true;
}
else
{
FatalErrorIn
(
"void slidingInterface::coupleInterface("
"polyTopoChange& ref) const"
) << "Face " << faceI << " in cut faces has neither a master "
<< "nor a slave. Error in the cutting algorithm on add."
<< abort(FatalError);
}
}
}
// Move the orphaned faces into the face zone
// Pout << "Orphaned master faces: " << orphanedMaster << endl;
// Pout << "Orphaned slave faces: " << orphanedSlave << endl;
label nOrphanedMasters = 0;
forAll (orphanedMaster, faceI)
{
if (orphanedMaster[faceI])
{
nOrphanedMasters++;
// Recover original orientation
ref.setAction
(
polyModifyFace
(
masterPatch[faceI], // new face
masterPatchAddr[faceI], // master face index
-1, // owner
-1, // neighbour
false, // flux flip
-1, // patch ID
false, // remove from zone
masterFaceZoneID_.index(), // zone ID
false // zone flip
)
);
}
}
label nOrphanedSlaves = 0;
forAll (orphanedSlave, faceI)
{
if (orphanedSlave[faceI])
{
nOrphanedSlaves++;
// Recover original orientation
ref.setAction
(
polyModifyFace
(
slavePatch[faceI], // new face
slavePatchAddr[faceI], // slave face index
-1, // owner
-1, // neighbour
false, // flux flip
-1, // patch ID
false, // remove from zone
slaveFaceZoneID_.index(), // zone ID
false // zone flip
)
);
}
}
if (debug)
{
Pout<< "Number of orphaned faces: "
<< "master = " << nOrphanedMasters << " out of "
<< orphanedMaster.size()
<< " slave = " << nOrphanedSlaves << " out of "
<< orphanedSlave.size() << endl;
}
// Finished coupling the plane of sliding interface.
// Modify faces influenced by the sliding interface
// These are the faces belonging to the master and slave cell
// layer which have not already been modified.
// Modification comes in three types:
// 1) eliminate the points that have been removed when the old sliding
// interface has been removed
// 2) Merge the slave points that have hit points on the master patch
// 3) Introduce new points resulting from the new sliding interface cut
// Collect all affected faces
// Master side
// Grab the list of faces in the layer
const labelList& masterStickOuts = masterStickOutFaces();
// Pout << "masterStickOuts: " << masterStickOuts << endl;
// Re-create the master stick-out faces
forAll (masterStickOuts, faceI)
{
// Renumber the face and remove additional points
const label curFaceID = masterStickOuts[faceI];
const face& oldRichFace = faces[curFaceID];
bool changed = false;
// Remove removed points from the face
face oldFace(oldRichFace.size());
label nOldFace = 0;
forAll (oldRichFace, pointI)
{
if (ref.pointRemoved(oldRichFace[pointI]))
{
changed = true;
}
else
{
// Point off patch
oldFace[nOldFace] = oldRichFace[pointI];
nOldFace++;
}
}
oldFace.setSize(nOldFace);
// Pout << "old rich master face: " << oldRichFace << " old face: " << oldFace << endl;
DynamicList<label> newFaceLabels(2*oldFace.size());
forAll (oldFace, pointI)
{
if (masterMeshPointMap.found(oldFace[pointI]))
{
// Point is in master patch. Add it
// If the point is a direct hit, grab its label; otherwise
// keep the original
if (pointMergeMap.found(oldFace[pointI]))
{
changed = true;
newFaceLabels.append
(
pointMergeMap.find(oldFace[pointI])()
);
}
else
{
newFaceLabels.append(oldFace[pointI]);
}
// Find if there are additional points inserted onto the edge
// between current point and the next point
// Algorithm:
// 1) Find all the edges in the master patch coming
// out of the current point.
// 2) If the next point in the face to pick the right edge
const label localFirstLabel =
masterMeshPointMap.find(oldFace[pointI])();
const labelList& curEdges = masterPointEdges[localFirstLabel];
const label nextLabel = oldFace.nextLabel(pointI);
Map<label>::const_iterator mmpmIter =
masterMeshPointMap.find(nextLabel);
if (mmpmIter != masterMeshPointMap.end())
{
// Pout << "found label pair " << oldFace[pointI] << " and " << nextLabel;
// Find the points on the edge between them
const label localNextLabel = mmpmIter();
forAll (curEdges, curEdgeI)
{
if
(
masterEdges[curEdges[curEdgeI]].otherVertex
(
localFirstLabel
)
== localNextLabel
)
{
// Pout << " found edge: " << curEdges[curEdgeI] << endl;
// Get points on current edge
const labelList& curPime = pime[curEdges[curEdgeI]];
if (curPime.size() > 0)
{
changed = true;
// Pout << "curPime: " << curPime << endl;
// Insert the edge points into the face
// in the correct order
const point& startPoint =
masterLocalPoints[localFirstLabel];
vector e =
masterLocalPoints[localNextLabel]
- startPoint;
e /= magSqr(e);
scalarField edgePointWeights(curPime.size());
forAll (curPime, curPimeI)
{
edgePointWeights[curPimeI] =
(
e
& (
pointMap.find(curPime[curPimeI])()
- startPoint
)
);
}
if (debug)
{
if
(
min(edgePointWeights) < 0
|| max(edgePointWeights) > 1
)
{
FatalErrorIn
(
"void slidingInterface::"
"coupleInterface("
"polyTopoChange& ref) const"
) << "Error in master stick-out edge "
<< "point collection."
<< abort(FatalError);
}
}
// Go through the points and collect
// them based on weights from lower to
// higher. This gives the correct
// order of points along the edge.
for
(
label passI = 0;
passI < edgePointWeights.size();
passI++
)
{
// Max weight can only be one, so
// the sorting is done by
// elimination.
label nextPoint = -1;
scalar dist = 2;
forAll (edgePointWeights, wI)
{
if (edgePointWeights[wI] < dist)
{
dist = edgePointWeights[wI];
nextPoint = wI;
}
}
// Insert the next point and reset
// its weight to exclude it from
// future picks
newFaceLabels.append(curPime[nextPoint]);
edgePointWeights[nextPoint] = GREAT;
}
}
break;
} // End of found edge
} // End of looking through current edges
}
}
else
{
newFaceLabels.append(oldFace[pointI]);
}
}
if (changed)
{
if (newFaceLabels.size() < 3)
{
FatalErrorIn
(
"void slidingInterface::coupleInterface("
"polyTopoChange& ref) const"
) << "Face " << curFaceID << " reduced to less than "
<< "3 points. Topological/cutting error A." << nl
<< "Old face: " << oldFace << " new face: " << newFaceLabels
<< abort(FatalError);
}
// Get face zone and its flip
label modifiedFaceZone = faceZones.whichZone(curFaceID);
bool modifiedFaceZoneFlip = false;
if (modifiedFaceZone >= 0)
{
modifiedFaceZoneFlip =
faceZones[modifiedFaceZone].flipMap()
[
faceZones[modifiedFaceZone].whichFace(curFaceID)
];
}
face newFace;
newFace.transfer(newFaceLabels.shrink());
// Pout << "Modifying master stick-out face " << curFaceID << " old face: " << oldFace << " new face: " << newFace << endl;
// Modify the face
ref.setAction
(
polyModifyFace
(
newFace, // modified face
curFaceID, // label of face being modified
own[curFaceID], // owner
nei[curFaceID], // neighbour
false, // face flip
mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
false, // remove from zone
modifiedFaceZone, // zone for face
modifiedFaceZoneFlip // face flip in zone
)
);
}
}
// Pout << "Finished master side" << endl;
// Slave side
// Grab the list of faces in the layer
const labelList& slaveStickOuts = slaveStickOutFaces();
// Pout << "slaveStickOuts: " << slaveStickOuts << endl;
const Map<label>& rpm = retiredPointMap();
// Re-create the slave stick-out faces
forAll (slaveStickOuts, faceI)
{
// Renumber the face and remove additional points
const label curFaceID = slaveStickOuts[faceI];
const face& oldRichFace = faces[curFaceID];
bool changed = false;
// Remove removed points from the face
face oldFace(oldRichFace.size());
label nOldFace = 0;
forAll (oldRichFace, pointI)
{
if
(
rpm.found(oldRichFace[pointI])
|| slaveMeshPointMap.found(oldRichFace[pointI])
)
{
// Point definitely live. Add it
oldFace[nOldFace] = oldRichFace[pointI];
nOldFace++;
}
else if
(
ref.pointRemoved(oldRichFace[pointI])
|| masterMeshPointMap.found(oldRichFace[pointI])
)
{
// Point removed and not on slave patch
// (first if takes care of that!) or
// point belonging to master patch
changed = true;
}
else
{
// Point off patch
oldFace[nOldFace] = oldRichFace[pointI];
nOldFace++;
}
}
oldFace.setSize(nOldFace);
DynamicList<label> newFaceLabels(2*oldFace.size());
// Pout << "old rich slave face: " << oldRichFace << " old face: " << oldFace << endl;
forAll (oldFace, pointI)
{
// Try to find the point in retired points
label curP = oldFace[pointI];
Map<label>::const_iterator rpmIter = rpm.find(oldFace[pointI]);
if (rpmIter != rpm.end())
{
changed = true;
curP = rpmIter();
}
if (slaveMeshPointMap.found(curP))
{
// Point is in slave patch. Add it
// If the point is a direct hit, grab its label; otherwise
// keep the original
if (pointMergeMap.found(curP))
{
changed = true;
newFaceLabels.append
(
pointMergeMap.find(curP)()
);
}
else
{
newFaceLabels.append(curP);
}
// Find if there are additional points inserted onto the edge
// between current point and the next point
// Algorithm:
// 1) Find all the edges in the slave patch coming
// out of the current point.
// 2) Use the next point in the face to pick the right edge
const label localFirstLabel =
slaveMeshPointMap.find(curP)();
const labelList& curEdges = slavePointEdges[localFirstLabel];
label nextLabel = oldFace.nextLabel(pointI);
Map<label>::const_iterator rpmNextIter =
rpm.find(nextLabel);
if (rpmNextIter != rpm.end())
{
nextLabel = rpmNextIter();
}
Map<label>::const_iterator mmpmIter =
slaveMeshPointMap.find(nextLabel);
if (mmpmIter != slaveMeshPointMap.end())
{
// Both points on the slave patch.
// Find the points on the edge between them
const label localNextLabel = mmpmIter();
forAll (curEdges, curEdgeI)
{
if
(
slaveEdges[curEdges[curEdgeI]].otherVertex
(
localFirstLabel
)
== localNextLabel
)
{
// Pout << " found edge: " << curEdges[curEdgeI] << endl;
// Get points on current edge
const labelList& curPise = pise[curEdges[curEdgeI]];
if (curPise.size() > 0)
{
changed = true;
// Pout << "curPise: " << curPise << endl;
// Insert the edge points into the face
// in the correct order
const point& startPoint =
projectedSlavePoints[localFirstLabel];
vector e =
projectedSlavePoints[localNextLabel]
- startPoint;
e /= magSqr(e);
scalarField edgePointWeights(curPise.size());
forAll (curPise, curPiseI)
{
edgePointWeights[curPiseI] =
(
e
& (
pointMap.find(curPise[curPiseI])()
- startPoint
)
);
}
if (debug)
{
if
(
min(edgePointWeights) < 0
|| max(edgePointWeights) > 1
)
{
FatalErrorIn
(
"void slidingInterface::"
"coupleInterface("
"polyTopoChange& ref) const"
) << "Error in slave stick-out edge "
<< "point collection."
<< abort(FatalError);
}
}
// Go through the points and collect
// them based on weights from lower to
// higher. This gives the correct
// order of points along the edge.
for
(
label passI = 0;
passI < edgePointWeights.size();
passI++
)
{
// Max weight can only be one, so
// the sorting is done by
// elimination.
label nextPoint = -1;
scalar dist = 2;
forAll (edgePointWeights, wI)
{
if (edgePointWeights[wI] < dist)
{
dist = edgePointWeights[wI];
nextPoint = wI;
}
}
// Insert the next point and reset
// its weight to exclude it from
// future picks
newFaceLabels.append(curPise[nextPoint]);
edgePointWeights[nextPoint] = GREAT;
}
}
break;
}
} // End of found edge
} // End of looking through current edges
}
else
{
newFaceLabels.append(oldFace[pointI]);
}
}
if (changed)
{
if (newFaceLabels.size() < 3)
{
FatalErrorIn
(
"void slidingInterface::coupleInterface("
"polyTopoChange& ref) const"
) << "Face " << curFaceID << " reduced to less than "
<< "3 points. Topological/cutting error B." << nl
<< "Old face: " << oldFace << " new face: " << newFaceLabels
<< abort(FatalError);
}
// Get face zone and its flip
label modifiedFaceZone = faceZones.whichZone(curFaceID);
bool modifiedFaceZoneFlip = false;
if (modifiedFaceZone >= 0)
{
modifiedFaceZoneFlip =
faceZones[modifiedFaceZone].flipMap()
[
faceZones[modifiedFaceZone].whichFace(curFaceID)
];
}
face newFace;
newFace.transfer(newFaceLabels.shrink());
// Pout << "Modifying slave stick-out face " << curFaceID << " old face: " << oldFace << " new face: " << newFace << endl;
// Modify the face
ref.setAction
(
polyModifyFace
(
newFace, // modified face
curFaceID, // label of face being modified
own[curFaceID], // owner
nei[curFaceID], // neighbour
false, // face flip
mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
false, // remove from zone
modifiedFaceZone, // zone for face
modifiedFaceZoneFlip // face flip in zone
)
);
}
}
// Activate and retire slave patch points
// This needs to be done last, so that the map of removed points
// does not get damaged by point modifications
if (!retiredPointMapPtr_)
{
FatalErrorIn
(
"void slidingInterface::coupleInterface("
"polyTopoChange& ref) const"
) << "Retired point map pointer not set."
<< abort(FatalError);
}
Map<label>& addToRpm = *retiredPointMapPtr_;
// Clear the old map
addToRpm.clear();
label nRetiredPoints = 0;
forAll (slaveMeshPoints, pointI)
{
if (pointMergeMap.found(slaveMeshPoints[pointI]))
{
// Retire the point - only used for supporting the detached
// slave patch
nRetiredPoints++;
ref.setAction
(
polyModifyPoint
(
slaveMeshPoints[pointI], // point ID
points[slaveMeshPoints[pointI]], // point
false, // remove from zone
mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
false // in a cell
)
);
addToRpm.insert
(
pointMergeMap.find(slaveMeshPoints[pointI])(),
slaveMeshPoints[pointI]
);
}
else
{
ref.setAction
(
polyModifyPoint
(
slaveMeshPoints[pointI], // point ID
points[slaveMeshPoints[pointI]], // point
false, // remove from zone
mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
true // in a cell
)
);
}
}
if (debug)
{
Pout << "Retired " << nRetiredPoints << " out of "
<< slaveMeshPoints.size() << " points." << endl;
}
// Grab cut face master and slave addressing
if (cutFaceMasterPtr_) deleteDemandDrivenData(cutFaceMasterPtr_);
cutFaceMasterPtr_ = new labelList(cutPatch.cutFaceMaster());
if (cutFaceSlavePtr_) deleteDemandDrivenData(cutFaceSlavePtr_);
cutFaceSlavePtr_ = new labelList(cutPatch.cutFaceSlave());
// Finished coupling
attached_ = true;
if (debug)
{
Pout<< "void slidingInterface::coupleInterface("
<< "polyTopoChange& ref) : "
<< "Finished coupling sliding interface " << name() << endl;
}
}
// ************************************************************************* //