mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
00b94c4a40
- addresses problems if faces are already partly attached (eg, on a single vertex or along an edge).
691 lines
26 KiB
C
691 lines
26 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
|
|
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
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 3 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, see <http://www.gnu.org/licenses/>.
|
|
|
|
Description
|
|
Calculating cut faces of the enriched patch, together with the addressing
|
|
into master and slave patch.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "enrichedPatch.H"
|
|
#include "boolList.H"
|
|
#include "DynamicList.H"
|
|
#include "labelPair.H"
|
|
#include "primitiveMesh.H"
|
|
#include "edgeHashes.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
// If the cut face gets more than this number of points, it will be checked
|
|
const Foam::label Foam::enrichedPatch::maxFaceSizeDebug_ = 100;
|
|
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
// Index of debug signs:
|
|
// x - skip a point
|
|
// l - left turn
|
|
// r - right turn
|
|
|
|
void Foam::enrichedPatch::calcCutFaces() const
|
|
{
|
|
if (cutFacesPtr_ || cutFaceMasterPtr_ || cutFaceSlavePtr_)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Cut faces addressing already calculated."
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
const faceList& enFaces = enrichedFaces();
|
|
const labelList& mp = meshPoints();
|
|
const faceList& lf = localFaces();
|
|
const pointField& lp = localPoints();
|
|
const labelListList& pp = pointPoints();
|
|
// Pout<< "enFaces: " << enFaces << endl;
|
|
// Pout<< "lf: " << lf << endl;
|
|
// Pout<< "lp: " << lp << endl;
|
|
// Pout<< "pp: " << pp << endl;
|
|
const Map<labelList>& masterPointFaceAddr = masterPointFaces();
|
|
|
|
// Prepare the storage
|
|
DynamicList<face> cf(2*lf.size());
|
|
DynamicList<label> cfMaster(2*lf.size());
|
|
DynamicList<label> cfSlave(2*lf.size());
|
|
|
|
// Algorithm
|
|
// Go through all the faces
|
|
// 1) For each face, start at point zero and grab the edge zero.
|
|
// Both points are added into the cut face. Mark the first edge
|
|
// as used and position the current point as the end of the first
|
|
// edge and previous point as point zero.
|
|
// 2) Grab all possible points from the current point. Excluding
|
|
// the previous point find the point giving the biggest right
|
|
// turn. Add the point to the face and mark the edges as used.
|
|
// Continue doing this until the face is complete, i.e. the next point
|
|
// to add is the first point of the face.
|
|
// 3) Once the facet is completed, register its owner from the face
|
|
// it has been created from (remember that the first lot of faces
|
|
// inserted into the enriched faces list are the slaves, to allow
|
|
// matching of the other side).
|
|
// 4) If the facet is created from an original slave face, go
|
|
// through the master patch and try to identify the master face
|
|
// this facet belongs to. This is recognised by the fact that the
|
|
// faces consists exclusively of the points of the master face and
|
|
// the points projected onto the face.
|
|
|
|
// Create a set of edge usage parameters
|
|
edgeHashSet edgesUsedOnce(pp.size());
|
|
edgeHashSet edgesUsedTwice(pp.size()*primitiveMesh::edgesPerPoint_);
|
|
|
|
|
|
forAll(lf, facei)
|
|
{
|
|
const face& curLocalFace = lf[facei];
|
|
const face& curGlobalFace = enFaces[facei];
|
|
|
|
// Pout<< "Doing face " << facei
|
|
// << " local: " << curLocalFace
|
|
// << " or " << curGlobalFace
|
|
// << endl;
|
|
|
|
// if (facei < slavePatch_.size())
|
|
// {
|
|
// Pout<< "original slave: " << slavePatch_[facei]
|
|
// << " local: " << slavePatch_.localFaces()[facei] << endl;
|
|
// }
|
|
// else
|
|
// {
|
|
// Pout<< "original master: "
|
|
// << masterPatch_[facei - slavePatch_.size()] << " "
|
|
// << masterPatch_.localFaces()[facei - slavePatch_.size()]
|
|
// << endl;
|
|
// }
|
|
// {
|
|
// pointField facePoints = curLocalFace.points(lp);
|
|
// forAll(curLocalFace, pointi)
|
|
// {
|
|
// Pout<< "v " << facePoints[pointi].x() << " "
|
|
// << facePoints[pointi].y() << " "
|
|
// << facePoints[pointi].z() << endl;
|
|
// }
|
|
// }
|
|
|
|
// Track the usage of face edges. When all edges are used, the
|
|
// face decomposition is complete.
|
|
// The seed edges include all the edges of the original face + all edges
|
|
// of other faces that have been used in the construction of the
|
|
// facet. Edges from other faces can be considered as
|
|
// internal to the current face if used only once.
|
|
|
|
// Track the edge usage to avoid duplicate faces and reset it to unused
|
|
boolList usedFaceEdges(curLocalFace.size(), false);
|
|
|
|
SLList<edge> edgeSeeds;
|
|
|
|
// Insert the edges of current face into the seed list.
|
|
edgeList cfe = curLocalFace.edges();
|
|
for (const edge& e : cfe)
|
|
{
|
|
edgeSeeds.append(e);
|
|
}
|
|
|
|
// Grab face normal
|
|
vector normal = curLocalFace.normal(lp);
|
|
normal /= mag(normal);
|
|
|
|
while (edgeSeeds.size())
|
|
{
|
|
// Pout<< "edgeSeeds.size(): " << edgeSeeds.size() << endl;
|
|
|
|
const edge curEdge = edgeSeeds.removeHead();
|
|
|
|
// Locate the edge in current face
|
|
const label curEdgeWhich = curLocalFace.which(curEdge.start());
|
|
|
|
// Check if the edge is in current face and if it has been
|
|
// used already. If so, skip it
|
|
if
|
|
(
|
|
curEdgeWhich > -1
|
|
&& curLocalFace.nextLabel(curEdgeWhich) == curEdge.end()
|
|
)
|
|
{
|
|
// Edge is in the starting face.
|
|
// If unused, mark it as used; if used, skip it
|
|
if (usedFaceEdges[curEdgeWhich]) continue;
|
|
|
|
usedFaceEdges[curEdgeWhich] = true;
|
|
}
|
|
|
|
// If the edge has already been used twice, skip it
|
|
if (edgesUsedTwice.found(curEdge)) continue;
|
|
|
|
// Pout<< "Trying new edge (" << mp[curEdge.start()]
|
|
// << ", " << mp[curEdge.end()]
|
|
// << ") seed: " << curEdge
|
|
// << " used: " << edgesUsedTwice.found(curEdge)
|
|
// << endl;
|
|
|
|
// Estimate the size of cut face as twice the size of original face
|
|
DynamicList<label> cutFaceGlobalPoints(2*curLocalFace.size());
|
|
DynamicList<label> cutFaceLocalPoints(2*curLocalFace.size());
|
|
|
|
// Found unused edge.
|
|
label prevPointLabel = curEdge.start();
|
|
cutFaceGlobalPoints.append(mp[prevPointLabel]);
|
|
cutFaceLocalPoints.append(prevPointLabel);
|
|
// Pout<< "prevPointLabel: " << mp[prevPointLabel] << endl;
|
|
|
|
// Grab current point and append it to the list
|
|
label curPointLabel = curEdge.end();
|
|
point curPoint = lp[curPointLabel];
|
|
cutFaceGlobalPoints.append(mp[curPointLabel]);
|
|
cutFaceLocalPoints.append(curPointLabel);
|
|
|
|
bool completedCutFace = false;
|
|
|
|
label faceSizeDebug = maxFaceSizeDebug_;
|
|
|
|
do
|
|
{
|
|
// Grab the next point options
|
|
|
|
// Pout<< "curPointLabel: " << mp[curPointLabel] << endl;
|
|
|
|
const labelList& nextPoints = pp[curPointLabel];
|
|
|
|
// Pout<< "nextPoints: "
|
|
// << UIndirectList<label>(mp, nextPoints)
|
|
// << endl;
|
|
|
|
// Get the vector along the edge and the right vector
|
|
vector ahead = curPoint - lp[prevPointLabel];
|
|
ahead -= normal*(normal & ahead);
|
|
ahead /= mag(ahead);
|
|
|
|
vector right = normal ^ ahead;
|
|
right /= mag(right);
|
|
|
|
// Pout<< "normal: " << normal
|
|
// << " ahead: " << ahead
|
|
// << " right: " << right
|
|
// << endl;
|
|
|
|
scalar atanTurn = -GREAT;
|
|
label bestAtanPoint = -1;
|
|
|
|
forAll(nextPoints, nextI)
|
|
{
|
|
// Exclude the point we are coming from; there will always
|
|
// be more than one edge, so this is safe
|
|
if (nextPoints[nextI] != prevPointLabel)
|
|
{
|
|
// Pout<< "cur point: " << curPoint
|
|
// << " trying for point: "
|
|
// << mp[nextPoints[nextI]]
|
|
// << " " << lp[nextPoints[nextI]];
|
|
vector newDir = lp[nextPoints[nextI]] - curPoint;
|
|
// Pout<< " newDir: " << newDir
|
|
// << " mag: " << mag(newDir) << flush;
|
|
newDir -= normal*(normal & newDir);
|
|
scalar magNewDir = mag(newDir);
|
|
// Pout<< " corrected: " << newDir
|
|
// << " mag: " << mag(newDir) << flush;
|
|
|
|
if (magNewDir < SMALL)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "projection error: slave patch probably "
|
|
<< "does not project onto master. "
|
|
<< "Please switch on "
|
|
<< "enriched patch debug for more info"
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
newDir /= magNewDir;
|
|
|
|
const scalar curAtanTurn =
|
|
atan2(newDir & right, newDir & ahead);
|
|
|
|
// Pout<< " atan: " << curAtanTurn << endl;
|
|
|
|
if (curAtanTurn > atanTurn)
|
|
{
|
|
bestAtanPoint = nextPoints[nextI];
|
|
atanTurn = curAtanTurn;
|
|
}
|
|
} // end of prev point skip
|
|
} // end of next point selection
|
|
|
|
// Pout<< " bestAtanPoint: " << bestAtanPoint << " or "
|
|
// << mp[bestAtanPoint]
|
|
// << endl;
|
|
|
|
// Selected next best point.
|
|
|
|
// Pout<< "cutFaceGlobalPoints: "
|
|
// << cutFaceGlobalPoints
|
|
// << endl;
|
|
|
|
// Check if the edge about to be added has been used
|
|
// in the current face or twice in other faces. If
|
|
// so, the face is bad.
|
|
const label whichNextPoint = curLocalFace.which(curPointLabel);
|
|
|
|
if
|
|
(
|
|
whichNextPoint > -1
|
|
&& curLocalFace.nextLabel(whichNextPoint) == bestAtanPoint
|
|
&& usedFaceEdges[whichNextPoint]
|
|
)
|
|
{
|
|
// This edge is already used in current face
|
|
// face cannot be good; start on a new one
|
|
|
|
// Pout<< "Double usage in current face, cannot be good"
|
|
// << endl;
|
|
|
|
completedCutFace = true;
|
|
}
|
|
|
|
|
|
if (edgesUsedTwice.found(edge(curPointLabel, bestAtanPoint)))
|
|
{
|
|
// This edge is already used -
|
|
// face cannot be good; start on a new one
|
|
completedCutFace = true;
|
|
|
|
// Pout<< "Double usage elsewhere, cannot be good" << endl;
|
|
}
|
|
|
|
if (completedCutFace) continue;
|
|
|
|
// Insert the next best point into the list
|
|
if (mp[bestAtanPoint] == cutFaceGlobalPoints[0])
|
|
{
|
|
// Face is completed. Save it and move on to the next
|
|
// available edge
|
|
completedCutFace = true;
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< " local: " << cutFaceLocalPoints
|
|
<< " one side: " << facei;
|
|
}
|
|
|
|
// Append the face
|
|
face cutFaceGlobal;
|
|
cutFaceGlobal.transfer(cutFaceGlobalPoints);
|
|
|
|
cf.append(cutFaceGlobal);
|
|
|
|
face cutFaceLocal;
|
|
cutFaceLocal.transfer(cutFaceLocalPoints);
|
|
|
|
// Pout<< "\ncutFaceLocal: "
|
|
// << cutFaceLocal.points(lp)
|
|
// << endl;
|
|
|
|
// Go through all edges of the cut faces.
|
|
// If the edge corresponds to a starting face edge,
|
|
// mark the starting face edge as true
|
|
|
|
forAll(cutFaceLocal, cuti)
|
|
{
|
|
const edge curCutFaceEdge
|
|
(
|
|
cutFaceLocal[cuti],
|
|
cutFaceLocal.nextLabel(cuti)
|
|
);
|
|
|
|
// Increment the usage count using two hash sets
|
|
auto euoIter = edgesUsedOnce.find(curCutFaceEdge);
|
|
|
|
if (!euoIter.found())
|
|
{
|
|
// Pout<< "Found edge not used before: "
|
|
// << curCutFaceEdge
|
|
// << endl;
|
|
edgesUsedOnce.insert(curCutFaceEdge);
|
|
}
|
|
else
|
|
{
|
|
// Pout<< "Found edge used once: "
|
|
// << curCutFaceEdge
|
|
// << endl;
|
|
edgesUsedOnce.erase(euoIter);
|
|
edgesUsedTwice.insert(curCutFaceEdge);
|
|
}
|
|
|
|
const label curCutFaceEdgeWhich = curLocalFace.which
|
|
(
|
|
curCutFaceEdge.start()
|
|
);
|
|
|
|
if
|
|
(
|
|
curCutFaceEdgeWhich > -1
|
|
&& curLocalFace.nextLabel(curCutFaceEdgeWhich)
|
|
== curCutFaceEdge.end()
|
|
)
|
|
{
|
|
// Found edge in original face
|
|
|
|
// Pout<< "Found edge in orig face: "
|
|
// << curCutFaceEdge << ": "
|
|
// << curCutFaceEdgeWhich
|
|
// << endl;
|
|
|
|
usedFaceEdges[curCutFaceEdgeWhich] = true;
|
|
}
|
|
else
|
|
{
|
|
// Edge not in original face. Add it to seeds
|
|
|
|
// Pout<< "Found new edge seed: "
|
|
// << curCutFaceEdge
|
|
// << endl;
|
|
|
|
edgeSeeds.append(curCutFaceEdge.reverseEdge());
|
|
}
|
|
}
|
|
|
|
// Find out what the other side is
|
|
|
|
// Algorithm
|
|
|
|
// If the face is in the slave half of the
|
|
// enrichedFaces lists, it may be matched against
|
|
// the master face. It can be recognised by the
|
|
// fact that all its points belong to one master
|
|
// face. For this purpose:
|
|
// 1) Grab the master faces around the point zero
|
|
// of the cut face and collect all master faces
|
|
// around it.
|
|
// 2) Loop through the rest of cut face points and
|
|
// if the point refers to one of the faces marked
|
|
// by point zero, increment its count.
|
|
// 3) When completed, the face whose count is
|
|
// equal to the number of points in the cut face
|
|
// is the other side. If this is not the case, there is no
|
|
// face on the other side.
|
|
|
|
if (facei < slavePatch_.size())
|
|
{
|
|
const auto mpfAddrIter =
|
|
masterPointFaceAddr.cfind(cutFaceGlobal[0]);
|
|
|
|
bool otherSideFound = false;
|
|
|
|
if (mpfAddrIter.found())
|
|
{
|
|
bool miss = false;
|
|
|
|
// Create the label count using point zero
|
|
const labelList& masterFacesOfPZero = mpfAddrIter();
|
|
|
|
labelList hits(masterFacesOfPZero.size(), 1);
|
|
|
|
for
|
|
(
|
|
label pointi = 1;
|
|
pointi < cutFaceGlobal.size();
|
|
pointi++
|
|
)
|
|
{
|
|
const auto mpfAddrPointIter =
|
|
masterPointFaceAddr.cfind
|
|
(
|
|
cutFaceGlobal[pointi]
|
|
);
|
|
|
|
if (!mpfAddrPointIter.found())
|
|
{
|
|
// Point is off the master patch. Skip
|
|
miss = true;
|
|
break;
|
|
}
|
|
|
|
const labelList& curMasterFaces =
|
|
mpfAddrPointIter();
|
|
|
|
// For every current face, try to find it in the
|
|
// zero-list
|
|
for (const label facei : curMasterFaces)
|
|
{
|
|
forAll(masterFacesOfPZero, j)
|
|
{
|
|
if (facei == masterFacesOfPZero[j])
|
|
{
|
|
hits[j]++;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// If all point are found attempt matching
|
|
if (!miss)
|
|
{
|
|
forAll(hits, pointi)
|
|
{
|
|
if (hits[pointi] == cutFaceGlobal.size())
|
|
{
|
|
// Found other side.
|
|
otherSideFound = true;
|
|
|
|
cfMaster.append
|
|
(
|
|
masterFacesOfPZero[pointi]
|
|
);
|
|
|
|
cfSlave.append(facei);
|
|
|
|
// Reverse the face such that it
|
|
// points out of the master patch
|
|
cf.last().flip();
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< " other side: "
|
|
<< masterFacesOfPZero[pointi]
|
|
<< endl;
|
|
}
|
|
} // end of hits
|
|
} // end of for all hits
|
|
|
|
} // end of not miss
|
|
|
|
if (!otherSideFound || miss)
|
|
{
|
|
if (debug)
|
|
{
|
|
Pout<< " solo slave A" << endl;
|
|
}
|
|
|
|
cfMaster.append(-1);
|
|
cfSlave.append(facei);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// First point not in master patch
|
|
if (debug)
|
|
{
|
|
Pout<< " solo slave B" << endl;
|
|
}
|
|
|
|
cfMaster.append(-1);
|
|
cfSlave.append(facei);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (debug)
|
|
{
|
|
Pout<< " master side" << endl;
|
|
}
|
|
|
|
cfMaster.append(facei - slavePatch_.size());
|
|
cfSlave.append(-1);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Face not completed. Prepare for the next point search
|
|
prevPointLabel = curPointLabel;
|
|
curPointLabel = bestAtanPoint;
|
|
curPoint = lp[curPointLabel];
|
|
cutFaceGlobalPoints.append(mp[curPointLabel]);
|
|
cutFaceLocalPoints.append(curPointLabel);
|
|
|
|
if (debug || cutFaceGlobalPoints.size() > faceSizeDebug)
|
|
{
|
|
faceSizeDebug *= 2;
|
|
|
|
// Check for duplicate points in the face
|
|
forAll(cutFaceGlobalPoints, checkI)
|
|
{
|
|
for
|
|
(
|
|
label checkJ = checkI + 1;
|
|
checkJ < cutFaceGlobalPoints.size();
|
|
checkJ++
|
|
)
|
|
{
|
|
if
|
|
(
|
|
cutFaceGlobalPoints[checkI]
|
|
== cutFaceGlobalPoints[checkJ]
|
|
)
|
|
{
|
|
// Shrink local points for debugging
|
|
cutFaceLocalPoints.shrink();
|
|
|
|
face origFace;
|
|
face origFaceLocal;
|
|
if (facei < slavePatch_.size())
|
|
{
|
|
origFace = slavePatch_[facei];
|
|
origFaceLocal =
|
|
slavePatch_.localFaces()[facei];
|
|
}
|
|
else
|
|
{
|
|
origFace =
|
|
masterPatch_
|
|
[facei - slavePatch_.size()];
|
|
|
|
origFaceLocal =
|
|
masterPatch_.localFaces()
|
|
[facei - slavePatch_.size()];
|
|
}
|
|
|
|
FatalErrorInFunction
|
|
<< "Duplicate point found in cut face. "
|
|
<< "Error in the face cutting "
|
|
<< "algorithm for global face "
|
|
<< origFace << " local face "
|
|
<< origFaceLocal << nl
|
|
<< "Slave size: " << slavePatch_.size()
|
|
<< " Master size: "
|
|
<< masterPatch_.size()
|
|
<< " index: " << facei << ".\n"
|
|
<< "Face: " << curGlobalFace << nl
|
|
<< "Cut face: " << cutFaceGlobalPoints
|
|
<< " local: " << cutFaceLocalPoints
|
|
<< nl << "Points: "
|
|
<< face(cutFaceLocalPoints).points(lp)
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
}
|
|
} // end of debug
|
|
}
|
|
} while (!completedCutFace);
|
|
} // end of used edges
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< " Finished face " << facei << endl;
|
|
}
|
|
|
|
} // end of local faces
|
|
|
|
// Re-pack lists into compact storage
|
|
cutFacesPtr_ = new faceList();
|
|
cutFacesPtr_->transfer(cf);
|
|
|
|
cutFaceMasterPtr_ = new labelList();
|
|
cutFaceMasterPtr_->transfer(cfMaster);
|
|
|
|
cutFaceSlavePtr_ = new labelList();
|
|
cutFaceSlavePtr_->transfer(cfSlave);
|
|
}
|
|
|
|
|
|
void Foam::enrichedPatch::clearCutFaces()
|
|
{
|
|
deleteDemandDrivenData(cutFacesPtr_);
|
|
deleteDemandDrivenData(cutFaceMasterPtr_);
|
|
deleteDemandDrivenData(cutFaceSlavePtr_);
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
const Foam::faceList& Foam::enrichedPatch::cutFaces() const
|
|
{
|
|
if (!cutFacesPtr_)
|
|
{
|
|
calcCutFaces();
|
|
}
|
|
|
|
return *cutFacesPtr_;
|
|
}
|
|
|
|
|
|
const Foam::labelList& Foam::enrichedPatch::cutFaceMaster() const
|
|
{
|
|
if (!cutFaceMasterPtr_)
|
|
{
|
|
calcCutFaces();
|
|
}
|
|
|
|
return *cutFaceMasterPtr_;
|
|
}
|
|
|
|
|
|
const Foam::labelList& Foam::enrichedPatch::cutFaceSlave() const
|
|
{
|
|
if (!cutFaceSlavePtr_)
|
|
{
|
|
calcCutFaces();
|
|
}
|
|
|
|
return *cutFaceSlavePtr_;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|