mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
631 lines
15 KiB
C
631 lines
15 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2010-2010 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 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/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "polyMeshTetDecomposition.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
// Note: the use of this tolerance is ad-hoc, there may be extreme
|
|
// cases where the resulting tetrahedra still have particle tracking
|
|
// problems.
|
|
const Foam::scalar Foam::polyMeshTetDecomposition::minTetQuality = SMALL;
|
|
|
|
|
|
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
|
|
|
|
Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
|
|
(
|
|
const polyMesh& mesh,
|
|
label fI,
|
|
const point& nCc,
|
|
scalar tol,
|
|
bool report
|
|
)
|
|
{
|
|
const faceList& pFaces = mesh.faces();
|
|
const pointField& pPts = mesh.points();
|
|
const vectorField& pC = mesh.cellCentres();
|
|
const labelList& pOwner = mesh.faceOwner();
|
|
|
|
const face& f = pFaces[fI];
|
|
|
|
label oCI = pOwner[fI];
|
|
|
|
const point& oCc = pC[oCI];
|
|
|
|
forAll(f, faceBasePtI)
|
|
{
|
|
scalar thisBaseMinTetQuality = VGREAT;
|
|
|
|
const point& tetBasePt = pPts[f[faceBasePtI]];
|
|
|
|
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
|
|
{
|
|
label facePtI = (tetPtI + faceBasePtI) % f.size();
|
|
label otherFacePtI = f.fcIndex(facePtI);
|
|
|
|
List<scalar> tetQualities(2, 0.0);
|
|
|
|
{
|
|
// owner cell tet
|
|
label ptAI = f[facePtI];
|
|
label ptBI = f[otherFacePtI];
|
|
|
|
const point& pA = pPts[ptAI];
|
|
const point& pB = pPts[ptBI];
|
|
|
|
tetPointRef tet(oCc, tetBasePt, pA, pB);
|
|
|
|
tetQualities[0] = tet.quality();
|
|
}
|
|
|
|
{
|
|
// neighbour cell tet
|
|
label ptAI = f[otherFacePtI];
|
|
label ptBI = f[facePtI];
|
|
|
|
const point& pA = pPts[ptAI];
|
|
const point& pB = pPts[ptBI];
|
|
|
|
tetPointRef tet(nCc, tetBasePt, pA, pB);
|
|
|
|
tetQualities[1] = tet.quality();
|
|
}
|
|
|
|
if (min(tetQualities) < thisBaseMinTetQuality)
|
|
{
|
|
thisBaseMinTetQuality = min(tetQualities);
|
|
}
|
|
}
|
|
|
|
if (thisBaseMinTetQuality > tol)
|
|
{
|
|
return faceBasePtI;
|
|
}
|
|
|
|
}
|
|
|
|
// If a base point hasn't triggered a return by now, then there is
|
|
// non that can produce a good decomposition
|
|
return -1;
|
|
}
|
|
|
|
|
|
Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
|
|
(
|
|
const polyMesh& mesh,
|
|
label fI,
|
|
scalar tol,
|
|
bool report
|
|
)
|
|
{
|
|
return findSharedBasePoint
|
|
(
|
|
mesh,
|
|
fI,
|
|
mesh.cellCentres()[mesh.faceNeighbour()[fI]],
|
|
tol,
|
|
report
|
|
);
|
|
}
|
|
|
|
|
|
Foam::label Foam::polyMeshTetDecomposition::findBasePoint
|
|
(
|
|
const polyMesh& mesh,
|
|
label fI,
|
|
scalar tol,
|
|
bool report
|
|
)
|
|
{
|
|
const faceList& pFaces = mesh.faces();
|
|
const pointField& pPts = mesh.points();
|
|
const vectorField& pC = mesh.cellCentres();
|
|
const labelList& pOwner = mesh.faceOwner();
|
|
|
|
const face& f = pFaces[fI];
|
|
|
|
label cI = pOwner[fI];
|
|
|
|
bool own = (pOwner[fI] == cI);
|
|
|
|
const point& cC = pC[cI];
|
|
|
|
forAll(f, faceBasePtI)
|
|
{
|
|
scalar thisBaseMinTetQuality = VGREAT;
|
|
|
|
const point& tetBasePt = pPts[f[faceBasePtI]];
|
|
|
|
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
|
|
{
|
|
label facePtI = (tetPtI + faceBasePtI) % f.size();
|
|
label otherFacePtI = f.fcIndex(facePtI);
|
|
|
|
label ptAI = -1;
|
|
label ptBI = -1;
|
|
|
|
if (own)
|
|
{
|
|
ptAI = f[facePtI];
|
|
ptBI = f[otherFacePtI];
|
|
}
|
|
else
|
|
{
|
|
ptAI = f[otherFacePtI];
|
|
ptBI = f[facePtI];
|
|
}
|
|
|
|
const point& pA = pPts[ptAI];
|
|
const point& pB = pPts[ptBI];
|
|
|
|
tetPointRef tet(cC, tetBasePt, pA, pB);
|
|
|
|
scalar tetQuality = tet.quality();
|
|
|
|
if (tetQuality < thisBaseMinTetQuality)
|
|
{
|
|
thisBaseMinTetQuality = tetQuality;
|
|
}
|
|
}
|
|
|
|
if (thisBaseMinTetQuality > tol)
|
|
{
|
|
return faceBasePtI;
|
|
}
|
|
}
|
|
|
|
// If a base point hasn't triggered a return by now, then there is
|
|
// non that can produce a good decomposition
|
|
return -1;
|
|
}
|
|
|
|
|
|
Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts
|
|
(
|
|
const polyMesh& mesh,
|
|
scalar tol,
|
|
bool report
|
|
)
|
|
{
|
|
const labelList& pOwner = mesh.faceOwner();
|
|
|
|
// Find a suitable base point for each face, considering both
|
|
// cells for interface faces or those on coupled patches
|
|
|
|
labelList tetBasePtIs(mesh.nFaces(), -1);
|
|
|
|
label nInternalFaces = mesh.nInternalFaces();
|
|
|
|
for (label fI = 0; fI < nInternalFaces; fI++)
|
|
{
|
|
tetBasePtIs[fI] = findSharedBasePoint(mesh, fI, tol, report);
|
|
}
|
|
|
|
pointField neighbourCellCentres(mesh.nFaces() - nInternalFaces);
|
|
|
|
for(label faceI = nInternalFaces; faceI < mesh.nFaces(); faceI++)
|
|
{
|
|
neighbourCellCentres[faceI - nInternalFaces] =
|
|
mesh.cellCentres()[pOwner[faceI]];
|
|
}
|
|
|
|
syncTools::swapBoundaryFacePositions(mesh, neighbourCellCentres);
|
|
|
|
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
|
|
|
SubList<label> boundaryFaceTetBasePtIs
|
|
(
|
|
tetBasePtIs,
|
|
mesh.nFaces() - nInternalFaces,
|
|
nInternalFaces
|
|
);
|
|
|
|
for
|
|
(
|
|
label fI = nInternalFaces, bFI = 0;
|
|
fI < mesh.nFaces();
|
|
fI++, bFI++
|
|
)
|
|
{
|
|
label patchI =
|
|
mesh.boundaryMesh().patchID()[bFI];
|
|
|
|
if (patches[patchI].coupled())
|
|
{
|
|
const coupledPolyPatch& pp =
|
|
refCast<const coupledPolyPatch>(patches[patchI]);
|
|
|
|
if (pp.owner())
|
|
{
|
|
boundaryFaceTetBasePtIs[bFI] = findSharedBasePoint
|
|
(
|
|
mesh,
|
|
fI,
|
|
neighbourCellCentres[bFI],
|
|
tol,
|
|
report
|
|
);
|
|
}
|
|
else
|
|
{
|
|
// Assign -2, to distinguish from a failed base point
|
|
// find, which returns -1.
|
|
boundaryFaceTetBasePtIs[bFI] = -2;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
boundaryFaceTetBasePtIs[bFI] = findBasePoint
|
|
(
|
|
mesh,
|
|
fI,
|
|
tol,
|
|
report
|
|
);
|
|
}
|
|
}
|
|
|
|
// maxEqOp will replace the -2 values on the neighbour patches
|
|
// with the result from the owner base point find.
|
|
|
|
syncTools::syncBoundaryFaceList
|
|
(
|
|
mesh,
|
|
boundaryFaceTetBasePtIs,
|
|
maxEqOp<label>()
|
|
);
|
|
|
|
for
|
|
(
|
|
label fI = nInternalFaces, bFI = 0;
|
|
fI < mesh.nFaces();
|
|
fI++, bFI++
|
|
)
|
|
{
|
|
label& bFTetBasePtI = boundaryFaceTetBasePtIs[bFI];
|
|
|
|
if (bFTetBasePtI == -2)
|
|
{
|
|
FatalErrorIn
|
|
(
|
|
"Foam::labelList"
|
|
"Foam::polyMeshTetDecomposition::findFaceBasePts"
|
|
"("
|
|
"const polyMesh& mesh, "
|
|
"scalar tol, "
|
|
"bool report"
|
|
")"
|
|
)
|
|
<< "Coupled face base point exchange failure for face "
|
|
<< fI
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
if (bFTetBasePtI < 1)
|
|
{
|
|
// If the base point is -1, it should be left as such to
|
|
// indicate a problem, if it is 0, then no action is required.
|
|
|
|
continue;
|
|
}
|
|
|
|
label patchI = mesh.boundaryMesh().patchID()[bFI];
|
|
|
|
if (patches[patchI].coupled())
|
|
{
|
|
const coupledPolyPatch& pp =
|
|
refCast<const coupledPolyPatch>(patches[patchI]);
|
|
|
|
// Calculated base points on coupled faces are those of
|
|
// the owner patch face. They need to be reindexed to for
|
|
// the non-owner face, which has the opposite order.
|
|
|
|
// So, for fPtI_o != 0, fPtI_n = f.size() - fPtI_o
|
|
|
|
// i.e.:
|
|
|
|
// owner coupledPolyPatch face
|
|
// face (a b c d e f)
|
|
// fPtI 0 1 2 3 4 5
|
|
// +
|
|
// #
|
|
|
|
// neighbour coupledPolyPatch face
|
|
// face (a f e d c b)
|
|
// fPtI 0 1 2 3 4 5
|
|
// +
|
|
// #
|
|
// +: 6 - 1 = 5
|
|
// #: 6 - 2 = 4
|
|
|
|
if (!pp.owner())
|
|
{
|
|
bFTetBasePtI = mesh.faces()[fI].size() - bFTetBasePtI;
|
|
}
|
|
}
|
|
}
|
|
|
|
return tetBasePtIs;
|
|
}
|
|
|
|
|
|
bool Foam::polyMeshTetDecomposition::checkFaceTets
|
|
(
|
|
const polyMesh& mesh,
|
|
scalar tol,
|
|
const bool report,
|
|
labelHashSet* setPtr
|
|
)
|
|
{
|
|
const labelList& own = mesh.faceOwner();
|
|
const labelList& nei = mesh.faceNeighbour();
|
|
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
|
|
|
const vectorField& cc = mesh.cellCentres();
|
|
const vectorField& fc = mesh.faceCentres();
|
|
|
|
// Calculate coupled cell centre
|
|
pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
|
|
|
|
for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
|
|
{
|
|
neiCc[faceI - mesh.nInternalFaces()] = cc[own[faceI]];
|
|
}
|
|
|
|
syncTools::swapBoundaryFacePositions(mesh, neiCc);
|
|
|
|
const faceList& fcs = mesh.faces();
|
|
|
|
const pointField& p = mesh.points();
|
|
|
|
label nErrorTets = 0;
|
|
|
|
forAll(fcs, faceI)
|
|
{
|
|
const face& f = fcs[faceI];
|
|
|
|
forAll(f, fPtI)
|
|
{
|
|
scalar tetQual = tetPointRef
|
|
(
|
|
p[f[fPtI]],
|
|
p[f.nextLabel(fPtI)],
|
|
fc[faceI],
|
|
cc[own[faceI]]
|
|
).quality();
|
|
|
|
if (tetQual > -tol)
|
|
{
|
|
if (setPtr)
|
|
{
|
|
setPtr->insert(faceI);
|
|
}
|
|
|
|
nErrorTets++;
|
|
break; // no need to check other tets
|
|
}
|
|
}
|
|
|
|
if (mesh.isInternalFace(faceI))
|
|
{
|
|
// Create the neighbour tet - it will have positive volume
|
|
const face& f = fcs[faceI];
|
|
|
|
forAll(f, fPtI)
|
|
{
|
|
scalar tetQual = tetPointRef
|
|
(
|
|
p[f[fPtI]],
|
|
p[f.nextLabel(fPtI)],
|
|
fc[faceI],
|
|
cc[nei[faceI]]
|
|
).quality();
|
|
|
|
if (tetQual < tol)
|
|
{
|
|
if (setPtr)
|
|
{
|
|
setPtr->insert(faceI);
|
|
}
|
|
|
|
nErrorTets++;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (findSharedBasePoint(mesh, faceI, tol, report) == -1)
|
|
{
|
|
if (setPtr)
|
|
{
|
|
setPtr->insert(faceI);
|
|
}
|
|
|
|
nErrorTets++;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
label patchI = patches.patchID()[faceI - mesh.nInternalFaces()];
|
|
|
|
if (patches[patchI].coupled())
|
|
{
|
|
if
|
|
(
|
|
findSharedBasePoint
|
|
(
|
|
mesh,
|
|
faceI,
|
|
neiCc[faceI - mesh.nInternalFaces()],
|
|
tol,
|
|
report
|
|
) == -1
|
|
)
|
|
{
|
|
if (setPtr)
|
|
{
|
|
setPtr->insert(faceI);
|
|
}
|
|
|
|
nErrorTets++;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (findBasePoint(mesh, faceI, tol, report) == -1)
|
|
{
|
|
if (setPtr)
|
|
{
|
|
setPtr->insert(faceI);
|
|
}
|
|
|
|
nErrorTets++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
reduce(nErrorTets, sumOp<label>());
|
|
|
|
if (nErrorTets > 0)
|
|
{
|
|
if (report)
|
|
{
|
|
Info<< " ***Error in face tets: "
|
|
<< nErrorTets << " faces with low quality or negative volume "
|
|
<< "decomposition tets." << endl;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
if (report)
|
|
{
|
|
Info<< " Face tets OK." << endl;
|
|
}
|
|
|
|
return false;
|
|
}
|
|
}
|
|
|
|
|
|
Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
|
|
(
|
|
const polyMesh& mesh,
|
|
label fI,
|
|
label cI
|
|
)
|
|
{
|
|
const faceList& pFaces = mesh.faces();
|
|
const labelList& pOwner = mesh.faceOwner();
|
|
|
|
const face& f = pFaces[fI];
|
|
|
|
label nTets = f.size() - 2;
|
|
|
|
List<tetIndices> faceTets(nTets);
|
|
|
|
bool own = (pOwner[fI] == cI);
|
|
|
|
label tetBasePtI = mesh.tetBasePtIs()[fI];
|
|
|
|
if (tetBasePtI == -1)
|
|
{
|
|
FatalErrorIn
|
|
(
|
|
"Foam::List<Foam::FixedList<Foam::label, 4> >"
|
|
"Foam::Cloud<ParticleType>::"
|
|
"faceTetIndices(label fI, label cI) const"
|
|
)
|
|
<< "No base point for face " << fI << ", " << f
|
|
<< ", produces a valid tet decomposition."
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
|
|
{
|
|
tetIndices& faceTetIs = faceTets[tetPtI - 1];
|
|
|
|
label facePtI = (tetPtI + tetBasePtI) % f.size();
|
|
label otherFacePtI = f.fcIndex(facePtI);
|
|
|
|
faceTetIs.cell() = cI;
|
|
|
|
faceTetIs.face() = fI;
|
|
|
|
faceTetIs.faceBasePt() = tetBasePtI;
|
|
|
|
if (own)
|
|
{
|
|
faceTetIs.facePtA() = facePtI;
|
|
faceTetIs.facePtB() = otherFacePtI;
|
|
}
|
|
else
|
|
{
|
|
faceTetIs.facePtA() = otherFacePtI;
|
|
faceTetIs.facePtB() = facePtI;
|
|
}
|
|
|
|
faceTetIs.tetPt() = tetPtI;
|
|
}
|
|
|
|
return faceTets;
|
|
}
|
|
|
|
|
|
Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
|
|
(
|
|
const polyMesh& mesh,
|
|
label cI
|
|
)
|
|
{
|
|
const faceList& pFaces = mesh.faces();
|
|
const cellList& pCells = mesh.cells();
|
|
|
|
const cell& thisCell = pCells[cI];
|
|
|
|
label nTets = 0;
|
|
|
|
forAll(thisCell, cFI)
|
|
{
|
|
nTets += pFaces[thisCell[cFI]].size() - 2;
|
|
}
|
|
|
|
DynamicList<tetIndices> cellTets(nTets);
|
|
|
|
forAll(thisCell, cFI)
|
|
{
|
|
label fI = thisCell[cFI];
|
|
|
|
cellTets.append(faceTetIndices(mesh, fI, cI));
|
|
}
|
|
|
|
return cellTets;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|