Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev
This commit is contained in:
Henry
commit
a07ac2dc89
@ -24,7 +24,7 @@ EXE_INC = \
|
||||
-IPrintTable \
|
||||
-I../vectorTools
|
||||
|
||||
EXE_LIBS = \
|
||||
LIB_LIBS = \
|
||||
-lmeshTools \
|
||||
-ledgeMesh \
|
||||
-lfileFormats \
|
||||
|
||||
@ -30,6 +30,7 @@ License
|
||||
#include "indexedCellChecks.H"
|
||||
|
||||
#include "CGAL/Exact_predicates_exact_constructions_kernel.h"
|
||||
#include "CGAL/Gmpq.h"
|
||||
|
||||
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -803,7 +803,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
|
||||
if (allGeometry)
|
||||
{
|
||||
cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
|
||||
if (mesh.checkCellDeterminant(true, &cells, mesh.geometricD()))
|
||||
if (mesh.checkCellDeterminant(true, &cells))
|
||||
{
|
||||
noFailedChecks++;
|
||||
|
||||
|
||||
@ -994,7 +994,7 @@ DimensionSets
|
||||
acceleration acceleration [ m s^-2 ] 1.0;
|
||||
kinematicPressure kinematicPressure [ Pa density^-1 ] 1.0;
|
||||
|
||||
// Scaled units. Only allowed in dymensionedType (dimensionedScalar,
|
||||
// Scaled units. Only allowed in dimensionedType (dimensionedScalar,
|
||||
// dimensionedVector etc.) and UniformDimensionedField, not
|
||||
// in DimensionedField or GeometricField
|
||||
cm cm [ m ] 1e-2;
|
||||
|
||||
@ -427,6 +427,10 @@ $(polyMesh)/polyMeshInitMesh.C
|
||||
$(polyMesh)/polyMeshClear.C
|
||||
$(polyMesh)/polyMeshUpdate.C
|
||||
|
||||
polyMeshCheck = $(polyMesh)/polyMeshCheck
|
||||
$(polyMeshCheck)/polyMeshCheck.C
|
||||
$(polyMeshCheck)/polyMeshTools.C
|
||||
|
||||
primitiveMesh = meshes/primitiveMesh
|
||||
$(primitiveMesh)/primitiveMesh.C
|
||||
$(primitiveMesh)/primitiveMeshCellCells.C
|
||||
@ -447,9 +451,9 @@ $(primitiveMesh)/primitiveMeshCalcCellShapes.C
|
||||
|
||||
primitiveMeshCheck = $(primitiveMesh)/primitiveMeshCheck
|
||||
$(primitiveMeshCheck)/primitiveMeshCheck.C
|
||||
$(primitiveMeshCheck)/primitiveMeshCheckMotion.C
|
||||
$(primitiveMeshCheck)/primitiveMeshCheckPointNearness.C
|
||||
$(primitiveMeshCheck)/primitiveMeshCheckEdgeLength.C
|
||||
$(primitiveMeshCheck)/primitiveMeshTools.C
|
||||
|
||||
primitivePatch = $(primitiveMesh)/primitivePatch
|
||||
$(primitivePatch)/patchZones.C
|
||||
|
||||
@ -1129,7 +1129,7 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
|
||||
if (debug)
|
||||
{
|
||||
// Check mesh motion
|
||||
if (primitiveMesh::checkMeshMotion(points_, true))
|
||||
if (checkMeshMotion(points_, true))
|
||||
{
|
||||
Info<< "tmp<scalarField> polyMesh::movePoints"
|
||||
<< "(const pointField&) : "
|
||||
|
||||
@ -34,6 +34,7 @@ SourceFiles
|
||||
polyMeshFromShapeMesh.C
|
||||
polyMeshIO.C
|
||||
polyMeshUpdate.C
|
||||
polyMeshCheck.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
@ -224,7 +225,45 @@ private:
|
||||
cellList& cells
|
||||
);
|
||||
|
||||
// Geometry checks
|
||||
|
||||
//- Check non-orthogonality
|
||||
bool checkFaceOrthogonality
|
||||
(
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs,
|
||||
const bool report,
|
||||
const bool detailedReport,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Check face skewness
|
||||
bool checkFaceSkewness
|
||||
(
|
||||
const pointField& points,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs,
|
||||
const bool report,
|
||||
const bool detailedReport,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
bool checkEdgeAlignment
|
||||
(
|
||||
const pointField& p,
|
||||
const bool report,
|
||||
const Vector<label>& directions,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
bool checkCellDeterminant
|
||||
(
|
||||
const vectorField& faceAreas,
|
||||
const bool report,
|
||||
labelHashSet* setPtr,
|
||||
const Vector<label>& meshD
|
||||
) const;
|
||||
|
||||
public:
|
||||
|
||||
@ -535,6 +574,48 @@ public:
|
||||
void removeFiles() const;
|
||||
|
||||
|
||||
// Geometric checks. Selectively override primitiveMesh functionality.
|
||||
|
||||
//- Check boundary for closedness
|
||||
virtual bool checkClosedBoundary(const bool report = false) const;
|
||||
|
||||
//- Check non-orthogonality
|
||||
virtual bool checkFaceOrthogonality
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check face skewness
|
||||
virtual bool checkFaceSkewness
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check edge alignment for 1D/2D cases
|
||||
virtual bool checkEdgeAlignment
|
||||
(
|
||||
const bool report,
|
||||
const Vector<label>& directions,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
virtual bool checkCellDeterminant
|
||||
(
|
||||
const bool report,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Check mesh motion for correctness given motion points
|
||||
virtual bool checkMeshMotion
|
||||
(
|
||||
const pointField& newPoints,
|
||||
const bool report = false,
|
||||
const bool detailedReport = false
|
||||
) const;
|
||||
|
||||
|
||||
// Helper functions
|
||||
|
||||
//- Find the cell, tetFaceI and tetPtI for the given position
|
||||
|
||||
679
src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
Normal file
679
src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
Normal file
@ -0,0 +1,679 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
|
||||
\\/ 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 "polyMesh.H"
|
||||
#include "polyMeshTools.H"
|
||||
#include "unitConversion.H"
|
||||
#include "syncTools.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
bool Foam::polyMesh::checkFaceOrthogonality
|
||||
(
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs,
|
||||
const bool report,
|
||||
const bool detailedReport,
|
||||
labelHashSet* setPtr
|
||||
) const
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Info<< "bool polyMesh::checkFaceOrthogonality("
|
||||
<< "const bool, labelHashSet*) const: "
|
||||
<< "checking mesh non-orthogonality" << endl;
|
||||
}
|
||||
|
||||
const labelList& own = faceOwner();
|
||||
const labelList& nei = faceNeighbour();
|
||||
|
||||
// Calculate orthogonality for all internal and coupled boundary faces
|
||||
// (1 for uncoupled boundary faces)
|
||||
tmp<scalarField> tortho = polyMeshTools::faceOrthogonality
|
||||
(
|
||||
*this,
|
||||
fAreas,
|
||||
cellCtrs
|
||||
);
|
||||
const scalarField& ortho = tortho();
|
||||
|
||||
// Severe nonorthogonality threshold
|
||||
const scalar severeNonorthogonalityThreshold =
|
||||
::cos(degToRad(primitiveMesh::nonOrthThreshold_));
|
||||
|
||||
|
||||
scalar minDDotS = GREAT;
|
||||
scalar sumDDotS = 0.0;
|
||||
label nSummed = 0;
|
||||
label severeNonOrth = 0;
|
||||
label errorNonOrth = 0;
|
||||
|
||||
|
||||
// Statistics only for internal and masters of coupled faces
|
||||
PackedBoolList isMasterFace(syncTools::getInternalOrMasterFaces(*this));
|
||||
|
||||
forAll(ortho, faceI)
|
||||
{
|
||||
if (ortho[faceI] < severeNonorthogonalityThreshold)
|
||||
{
|
||||
if (ortho[faceI] > SMALL)
|
||||
{
|
||||
if (setPtr)
|
||||
{
|
||||
setPtr->insert(faceI);
|
||||
}
|
||||
|
||||
severeNonOrth++;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Error : non-ortho too large
|
||||
if (setPtr)
|
||||
{
|
||||
setPtr->insert(faceI);
|
||||
}
|
||||
if (detailedReport && errorNonOrth == 0)
|
||||
{
|
||||
// Non-orthogonality greater than 90 deg
|
||||
WarningIn
|
||||
(
|
||||
"polyMesh::checkFaceOrthogonality"
|
||||
"(const pointField&, const bool) const"
|
||||
) << "Severe non-orthogonality for face "
|
||||
<< faceI
|
||||
<< " between cells " << own[faceI]
|
||||
<< " and " << nei[faceI]
|
||||
<< ": Angle = " << radToDeg(::acos(ortho[faceI]))
|
||||
<< " deg." << endl;
|
||||
}
|
||||
|
||||
errorNonOrth++;
|
||||
}
|
||||
}
|
||||
|
||||
if (isMasterFace[faceI])
|
||||
{
|
||||
minDDotS = min(minDDotS, ortho[faceI]);
|
||||
sumDDotS += ortho[faceI];
|
||||
nSummed++;
|
||||
}
|
||||
}
|
||||
|
||||
reduce(minDDotS, minOp<scalar>());
|
||||
reduce(sumDDotS, sumOp<scalar>());
|
||||
reduce(nSummed, sumOp<label>());
|
||||
reduce(severeNonOrth, sumOp<label>());
|
||||
reduce(errorNonOrth, sumOp<label>());
|
||||
|
||||
if (debug || report)
|
||||
{
|
||||
if (nSummed > 0)
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Info<< " Mesh non-orthogonality Max: "
|
||||
<< radToDeg(::acos(minDDotS))
|
||||
<< " average: " << radToDeg(::acos(sumDDotS/nSummed))
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
|
||||
if (severeNonOrth > 0)
|
||||
{
|
||||
Info<< " *Number of severely non-orthogonal faces: "
|
||||
<< severeNonOrth << "." << endl;
|
||||
}
|
||||
}
|
||||
|
||||
if (errorNonOrth > 0)
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Info<< " ***Number of non-orthogonality errors: "
|
||||
<< errorNonOrth << "." << endl;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Info<< " Non-orthogonality check OK." << endl;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
bool Foam::polyMesh::checkFaceSkewness
|
||||
(
|
||||
const pointField& points,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs,
|
||||
const bool report,
|
||||
const bool detailedReport,
|
||||
labelHashSet* setPtr
|
||||
) const
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Info<< "bool polyMesh::checkFaceSkewnesss("
|
||||
<< "const bool, labelHashSet*) const: "
|
||||
<< "checking face skewness" << endl;
|
||||
}
|
||||
|
||||
const labelList& own = faceOwner();
|
||||
const labelList& nei = faceNeighbour();
|
||||
|
||||
// Warn if the skew correction vector is more than skewWarning times
|
||||
// larger than the face area vector
|
||||
|
||||
tmp<scalarField> tskew = polyMeshTools::faceSkewness
|
||||
(
|
||||
*this,
|
||||
points,
|
||||
fCtrs,
|
||||
fAreas,
|
||||
cellCtrs
|
||||
);
|
||||
const scalarField& skew = tskew();
|
||||
|
||||
scalar maxSkew = max(skew);
|
||||
label nWarnSkew = 0;
|
||||
|
||||
// Statistics only for all faces except slave coupled faces
|
||||
PackedBoolList isMasterFace(syncTools::getMasterFaces(*this));
|
||||
|
||||
forAll(skew, faceI)
|
||||
{
|
||||
// Check if the skewness vector is greater than the PN vector.
|
||||
// This does not cause trouble but is a good indication of a poor mesh.
|
||||
if (skew[faceI] > skewThreshold_)
|
||||
{
|
||||
if (setPtr)
|
||||
{
|
||||
setPtr->insert(faceI);
|
||||
}
|
||||
if (detailedReport && nWarnSkew == 0)
|
||||
{
|
||||
// Non-orthogonality greater than 90 deg
|
||||
if (isInternalFace(faceI))
|
||||
{
|
||||
WarningIn
|
||||
(
|
||||
"polyMesh::checkFaceSkewnesss"
|
||||
"(const pointField&, const bool) const"
|
||||
) << "Severe skewness " << skew[faceI]
|
||||
<< " for face " << faceI
|
||||
<< " between cells " << own[faceI]
|
||||
<< " and " << nei[faceI];
|
||||
}
|
||||
else
|
||||
{
|
||||
WarningIn
|
||||
(
|
||||
"polyMesh::checkFaceSkewnesss"
|
||||
"(const pointField&, const bool) const"
|
||||
) << "Severe skewness " << skew[faceI]
|
||||
<< " for boundary face " << faceI
|
||||
<< " on cell " << own[faceI];
|
||||
}
|
||||
}
|
||||
|
||||
if (isMasterFace[faceI])
|
||||
{
|
||||
nWarnSkew++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
reduce(maxSkew, maxOp<scalar>());
|
||||
reduce(nWarnSkew, sumOp<label>());
|
||||
|
||||
if (nWarnSkew > 0)
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Info<< " ***Max skewness = " << maxSkew
|
||||
<< ", " << nWarnSkew << " highly skew faces detected"
|
||||
" which may impair the quality of the results"
|
||||
<< endl;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Info<< " Max skewness = " << maxSkew << " OK." << endl;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Check 1D/2Dness of edges. Gets passed the non-empty directions and
|
||||
// checks all edges in the mesh whether they:
|
||||
// - have no component in a non-empty direction or
|
||||
// - are only in a singe non-empty direction.
|
||||
// Empty direction info is passed in as a vector of labels (synchronised)
|
||||
// which are 1 if the direction is non-empty, 0 if it is.
|
||||
bool Foam::polyMesh::checkEdgeAlignment
|
||||
(
|
||||
const pointField& p,
|
||||
const bool report,
|
||||
const Vector<label>& directions,
|
||||
labelHashSet* setPtr
|
||||
) const
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Info<< "bool polyMesh::checkEdgeAlignment("
|
||||
<< "const bool, const Vector<label>&, labelHashSet*) const: "
|
||||
<< "checking edge alignment" << endl;
|
||||
}
|
||||
|
||||
label nDirs = 0;
|
||||
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
|
||||
{
|
||||
if (directions[cmpt] == 1)
|
||||
{
|
||||
nDirs++;
|
||||
}
|
||||
else if (directions[cmpt] != 0)
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"polyMesh::checkEdgeAlignment"
|
||||
"(const bool, const Vector<label>&, labelHashSet*)"
|
||||
) << "directions should contain 0 or 1 but is now " << directions
|
||||
<< exit(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
if (nDirs == vector::nComponents)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
const faceList& fcs = faces();
|
||||
|
||||
EdgeMap<label> edgesInError;
|
||||
|
||||
forAll(fcs, faceI)
|
||||
{
|
||||
const face& f = fcs[faceI];
|
||||
|
||||
forAll(f, fp)
|
||||
{
|
||||
label p0 = f[fp];
|
||||
label p1 = f.nextLabel(fp);
|
||||
if (p0 < p1)
|
||||
{
|
||||
vector d(p[p1]-p[p0]);
|
||||
scalar magD = mag(d);
|
||||
|
||||
if (magD > ROOTVSMALL)
|
||||
{
|
||||
d /= magD;
|
||||
|
||||
// Check how many empty directions are used by the edge.
|
||||
label nEmptyDirs = 0;
|
||||
label nNonEmptyDirs = 0;
|
||||
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
|
||||
{
|
||||
if (mag(d[cmpt]) > 1e-6)
|
||||
{
|
||||
if (directions[cmpt] == 0)
|
||||
{
|
||||
nEmptyDirs++;
|
||||
}
|
||||
else
|
||||
{
|
||||
nNonEmptyDirs++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (nEmptyDirs == 0)
|
||||
{
|
||||
// Purely in ok directions.
|
||||
}
|
||||
else if (nEmptyDirs == 1)
|
||||
{
|
||||
// Ok if purely in empty directions.
|
||||
if (nNonEmptyDirs > 0)
|
||||
{
|
||||
edgesInError.insert(edge(p0, p1), faceI);
|
||||
}
|
||||
}
|
||||
else if (nEmptyDirs > 1)
|
||||
{
|
||||
// Always an error
|
||||
edgesInError.insert(edge(p0, p1), faceI);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
label nErrorEdges = returnReduce(edgesInError.size(), sumOp<label>());
|
||||
|
||||
if (nErrorEdges > 0)
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Info<< " ***Number of edges not aligned with or perpendicular to "
|
||||
<< "non-empty directions: " << nErrorEdges << endl;
|
||||
}
|
||||
|
||||
if (setPtr)
|
||||
{
|
||||
setPtr->resize(2*edgesInError.size());
|
||||
forAllConstIter(EdgeMap<label>, edgesInError, iter)
|
||||
{
|
||||
setPtr->insert(iter.key()[0]);
|
||||
setPtr->insert(iter.key()[1]);
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Info<< " All edges aligned with or perpendicular to "
|
||||
<< "non-empty directions." << endl;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
bool Foam::polyMesh::checkCellDeterminant
|
||||
(
|
||||
const vectorField& faceAreas,
|
||||
const bool report,
|
||||
labelHashSet* setPtr,
|
||||
const Vector<label>& meshD
|
||||
) const
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Info<< "bool polyMesh::checkCellDeterminant(const bool"
|
||||
<< ", labelHashSet*) const: "
|
||||
<< "checking for under-determined cells" << endl;
|
||||
}
|
||||
|
||||
// Determine number of dimensions and (for 2D) missing dimension
|
||||
label nDims = 0;
|
||||
label twoD = -1;
|
||||
for (direction dir = 0; dir < vector::nComponents; dir++)
|
||||
{
|
||||
if (meshD[dir] == 1)
|
||||
{
|
||||
nDims++;
|
||||
}
|
||||
else
|
||||
{
|
||||
twoD = dir;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
tmp<scalarField> tcellDeterminant = primitiveMeshTools::cellDeterminant
|
||||
(
|
||||
*this,
|
||||
meshD,
|
||||
faceAreas,
|
||||
syncTools::getInternalOrCoupledFaces(*this)
|
||||
);
|
||||
scalarField& cellDeterminant = tcellDeterminant();
|
||||
|
||||
|
||||
label nErrorCells = 0;
|
||||
scalar minDet = min(cellDeterminant);
|
||||
scalar sumDet = sum(cellDeterminant);
|
||||
|
||||
forAll (cellDeterminant, cellI)
|
||||
{
|
||||
if (cellDeterminant[cellI] < 1e-3)
|
||||
{
|
||||
if (setPtr)
|
||||
{
|
||||
setPtr->insert(cellI);
|
||||
}
|
||||
|
||||
nErrorCells++;
|
||||
}
|
||||
}
|
||||
|
||||
reduce(nErrorCells, sumOp<label>());
|
||||
reduce(minDet, minOp<scalar>());
|
||||
reduce(sumDet, sumOp<scalar>());
|
||||
label nSummed = returnReduce(cellDeterminant.size(), sumOp<label>());
|
||||
|
||||
if (debug || report)
|
||||
{
|
||||
if (nSummed > 0)
|
||||
{
|
||||
Info<< " Cell determinant (wellposedness) : minimum: " << minDet
|
||||
<< " average: " << sumDet/nSummed
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
|
||||
if (nErrorCells > 0)
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Info<< " ***Cells with small determinant found, number of cells: "
|
||||
<< nErrorCells << endl;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Info<< " Cell determinant check OK." << endl;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
bool Foam::polyMesh::checkClosedBoundary(const bool report) const
|
||||
{
|
||||
return primitiveMesh::checkClosedBoundary
|
||||
(
|
||||
faceAreas(),
|
||||
report,
|
||||
syncTools::getInternalOrCoupledFaces(*this)
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
bool Foam::polyMesh::checkFaceOrthogonality
|
||||
(
|
||||
const bool report,
|
||||
labelHashSet* setPtr
|
||||
) const
|
||||
{
|
||||
return checkFaceOrthogonality
|
||||
(
|
||||
faceAreas(),
|
||||
cellCentres(),
|
||||
report,
|
||||
false, // detailedReport
|
||||
setPtr
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
bool Foam::polyMesh::checkFaceSkewness
|
||||
(
|
||||
const bool report,
|
||||
labelHashSet* setPtr
|
||||
) const
|
||||
{
|
||||
return checkFaceSkewness
|
||||
(
|
||||
points(),
|
||||
faceCentres(),
|
||||
faceAreas(),
|
||||
cellCentres(),
|
||||
report,
|
||||
false, // detailedReport
|
||||
setPtr
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
bool Foam::polyMesh::checkEdgeAlignment
|
||||
(
|
||||
const bool report,
|
||||
const Vector<label>& directions,
|
||||
labelHashSet* setPtr
|
||||
) const
|
||||
{
|
||||
return checkEdgeAlignment
|
||||
(
|
||||
points(),
|
||||
report,
|
||||
directions,
|
||||
setPtr
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
bool Foam::polyMesh::checkCellDeterminant
|
||||
(
|
||||
const bool report,
|
||||
labelHashSet* setPtr
|
||||
) const
|
||||
{
|
||||
return checkCellDeterminant
|
||||
(
|
||||
faceAreas(),
|
||||
report,
|
||||
setPtr,
|
||||
geometricD()
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
bool Foam::polyMesh::checkMeshMotion
|
||||
(
|
||||
const pointField& newPoints,
|
||||
const bool report,
|
||||
const bool detailedReport
|
||||
) const
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Pout<< "bool polyMesh::checkMeshMotion("
|
||||
<< "const pointField&, const bool, const bool) const: "
|
||||
<< "checking mesh motion" << endl;
|
||||
}
|
||||
|
||||
vectorField fCtrs(nFaces());
|
||||
vectorField fAreas(nFaces());
|
||||
|
||||
makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
|
||||
|
||||
// Check cell volumes and calculate new cell centres
|
||||
vectorField cellCtrs(nCells());
|
||||
scalarField cellVols(nCells());
|
||||
|
||||
makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
|
||||
|
||||
// Check cell volumes
|
||||
bool error = checkCellVolumes
|
||||
(
|
||||
cellVols, // vols
|
||||
report, // report
|
||||
detailedReport, // detailedReport
|
||||
NULL // setPtr
|
||||
);
|
||||
|
||||
|
||||
// Check face areas
|
||||
bool areaError = checkFaceAreas
|
||||
(
|
||||
faceAreas(),
|
||||
report, // report
|
||||
detailedReport, // detailedReport,
|
||||
NULL // setPtr
|
||||
);
|
||||
error = error || areaError;
|
||||
|
||||
|
||||
// Check pyramid volumes
|
||||
bool pyrVolError = checkFacePyramids
|
||||
(
|
||||
newPoints,
|
||||
cellCtrs,
|
||||
report, // report,
|
||||
detailedReport, // detailedReport,
|
||||
-SMALL, // minPyrVol
|
||||
NULL // setPtr
|
||||
);
|
||||
error = error || pyrVolError;
|
||||
|
||||
|
||||
// Check face non-orthogonality
|
||||
bool nonOrthoError = checkFaceOrthogonality
|
||||
(
|
||||
fAreas,
|
||||
cellCtrs,
|
||||
report, // report
|
||||
detailedReport, // detailedReport
|
||||
NULL // setPtr
|
||||
);
|
||||
error = error || nonOrthoError;
|
||||
|
||||
|
||||
if (!error && (debug || report))
|
||||
{
|
||||
Pout<< "Mesh motion check OK." << endl;
|
||||
}
|
||||
|
||||
return error;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
190
src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
Normal file
190
src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
Normal file
@ -0,0 +1,190 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
|
||||
\\/ 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 "polyMeshTools.H"
|
||||
#include "syncTools.H"
|
||||
#include "pyramidPointFaceRef.H"
|
||||
#include "primitiveMeshTools.H"
|
||||
#include "polyMeshTools.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceOrthogonality
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const vectorField& areas,
|
||||
const vectorField& cc
|
||||
)
|
||||
{
|
||||
const labelList& own = mesh.faceOwner();
|
||||
const labelList& nei = mesh.faceNeighbour();
|
||||
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
|
||||
|
||||
tmp<scalarField> tortho(new scalarField(mesh.nFaces(), 1.0));
|
||||
scalarField& ortho = tortho();
|
||||
|
||||
// Internal faces
|
||||
forAll(nei, faceI)
|
||||
{
|
||||
ortho[faceI] = primitiveMeshTools::faceOrthogonality
|
||||
(
|
||||
cc[own[faceI]],
|
||||
cc[nei[faceI]],
|
||||
areas[faceI]
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// Coupled faces
|
||||
|
||||
pointField neighbourCc;
|
||||
syncTools::swapBoundaryCellList(mesh, cc, neighbourCc);
|
||||
|
||||
forAll(pbm, patchI)
|
||||
{
|
||||
const polyPatch& pp = pbm[patchI];
|
||||
if (pp.coupled())
|
||||
{
|
||||
forAll(pp, i)
|
||||
{
|
||||
label faceI = pp.start() + i;
|
||||
label bFaceI = faceI - mesh.nInternalFaces();
|
||||
|
||||
ortho[faceI] = primitiveMeshTools::faceOrthogonality
|
||||
(
|
||||
cc[own[faceI]],
|
||||
neighbourCc[bFaceI],
|
||||
areas[faceI]
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return tortho;
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const pointField& p,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs
|
||||
)
|
||||
{
|
||||
const labelList& own = mesh.faceOwner();
|
||||
const labelList& nei = mesh.faceNeighbour();
|
||||
const faceList& fcs = mesh.faces();
|
||||
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
|
||||
|
||||
tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
|
||||
scalarField& skew = tskew();
|
||||
|
||||
forAll(nei, faceI)
|
||||
{
|
||||
skew[faceI] = primitiveMeshTools::faceSkewness
|
||||
(
|
||||
mesh,
|
||||
p,
|
||||
fCtrs,
|
||||
fAreas,
|
||||
|
||||
faceI,
|
||||
cellCtrs[own[faceI]],
|
||||
cellCtrs[nei[faceI]]
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// Boundary faces: consider them to have only skewness error.
|
||||
// (i.e. treat as if mirror cell on other side)
|
||||
|
||||
pointField neighbourCc;
|
||||
syncTools::swapBoundaryCellList(mesh, cellCtrs, neighbourCc);
|
||||
|
||||
forAll(pbm, patchI)
|
||||
{
|
||||
const polyPatch& pp = pbm[patchI];
|
||||
if (pp.coupled())
|
||||
{
|
||||
forAll(pp, i)
|
||||
{
|
||||
label faceI = pp.start() + i;
|
||||
label bFaceI = faceI - mesh.nInternalFaces();
|
||||
|
||||
skew[faceI] = primitiveMeshTools::faceSkewness
|
||||
(
|
||||
mesh,
|
||||
p,
|
||||
fCtrs,
|
||||
fAreas,
|
||||
|
||||
faceI,
|
||||
cellCtrs[own[faceI]],
|
||||
neighbourCc[bFaceI]
|
||||
);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
forAll(pp, i)
|
||||
{
|
||||
label faceI = pp.start() + i;
|
||||
|
||||
vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
|
||||
|
||||
vector normal = fAreas[faceI];
|
||||
normal /= mag(normal) + VSMALL;
|
||||
vector d = normal*(normal & Cpf);
|
||||
|
||||
|
||||
// Skewness vector
|
||||
vector sv =
|
||||
Cpf
|
||||
- ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
|
||||
vector svHat = sv/(mag(sv) + VSMALL);
|
||||
|
||||
// Normalisation distance calculated as the approximate distance
|
||||
// from the face centre to the edge of the face in the direction
|
||||
// of the skewness
|
||||
scalar fd = 0.4*mag(d) + VSMALL;
|
||||
const face& f = fcs[faceI];
|
||||
forAll(f, pi)
|
||||
{
|
||||
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
|
||||
}
|
||||
|
||||
// Normalised skewness
|
||||
skew[faceI] = mag(sv)/fd;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return tskew;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
109
src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
Normal file
109
src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
Normal file
@ -0,0 +1,109 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
|
||||
\\/ 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/>.
|
||||
|
||||
Namespace
|
||||
Foam::polyMeshTools
|
||||
|
||||
Description
|
||||
Collection of static functions operating on polyMesh (mainly checks) so
|
||||
that need access to patch information.
|
||||
|
||||
SourceFiles
|
||||
polyMeshTools.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
#ifndef polyMeshTools_H
|
||||
#define polyMeshTools_H
|
||||
|
||||
#include "polyMesh.H"
|
||||
#include "primitiveMeshTools.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Namespace polyMeshTools Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class polyMeshTools
|
||||
:
|
||||
public primitiveMeshTools
|
||||
{
|
||||
|
||||
public:
|
||||
|
||||
//- Generate orthogonality field. (1 for fully orthogonal, < 1 for
|
||||
// non-orthogonal)
|
||||
static tmp<scalarField> faceOrthogonality
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs
|
||||
);
|
||||
|
||||
// static tmp<scalarField> faceOrthogonality(const polyMesh& mesh)
|
||||
// {
|
||||
// return faceOrthogonality
|
||||
// (
|
||||
// mesh,
|
||||
// mesh.faceAreas(),
|
||||
// mesh.cellCentres()
|
||||
// );
|
||||
// }
|
||||
|
||||
//- Generate skewness field
|
||||
static tmp<scalarField> faceSkewness
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const pointField& points,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs
|
||||
);
|
||||
|
||||
// static tmp<scalarField> faceSkewness(const polyMesh& mesh)
|
||||
// {
|
||||
// return faceSkewness
|
||||
// (
|
||||
// mesh,
|
||||
// mesh.points(),
|
||||
// mesh.faceCentres(),
|
||||
// mesh.faceAreas(),
|
||||
// mesh.cellCentres()
|
||||
// );
|
||||
// }
|
||||
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -27,7 +27,6 @@ License
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
// Determines for every point whether it is coupled and if so sets only one.
|
||||
Foam::PackedBoolList Foam::syncTools::getMasterPoints(const polyMesh& mesh)
|
||||
{
|
||||
PackedBoolList isMasterPoint(mesh.nPoints());
|
||||
@ -72,7 +71,6 @@ Foam::PackedBoolList Foam::syncTools::getMasterPoints(const polyMesh& mesh)
|
||||
}
|
||||
|
||||
|
||||
// Determines for every edge whether it is coupled and if so sets only one.
|
||||
Foam::PackedBoolList Foam::syncTools::getMasterEdges(const polyMesh& mesh)
|
||||
{
|
||||
PackedBoolList isMasterEdge(mesh.nEdges());
|
||||
@ -117,7 +115,6 @@ Foam::PackedBoolList Foam::syncTools::getMasterEdges(const polyMesh& mesh)
|
||||
}
|
||||
|
||||
|
||||
// Determines for every face whether it is coupled and if so sets only one.
|
||||
Foam::PackedBoolList Foam::syncTools::getMasterFaces(const polyMesh& mesh)
|
||||
{
|
||||
PackedBoolList isMasterFace(mesh.nFaces(), 1);
|
||||
@ -145,4 +142,66 @@ Foam::PackedBoolList Foam::syncTools::getMasterFaces(const polyMesh& mesh)
|
||||
}
|
||||
|
||||
|
||||
Foam::PackedBoolList Foam::syncTools::getInternalOrMasterFaces
|
||||
(
|
||||
const polyMesh& mesh
|
||||
)
|
||||
{
|
||||
PackedBoolList isMasterFace(mesh.nFaces(), 1);
|
||||
|
||||
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
||||
|
||||
forAll(patches, patchI)
|
||||
{
|
||||
const polyPatch& pp = patches[patchI];
|
||||
|
||||
if (pp.coupled())
|
||||
{
|
||||
if (!refCast<const coupledPolyPatch>(pp).owner())
|
||||
{
|
||||
forAll(pp, i)
|
||||
{
|
||||
isMasterFace.unset(pp.start()+i);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
forAll(pp, i)
|
||||
{
|
||||
isMasterFace.unset(pp.start()+i);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return isMasterFace;
|
||||
}
|
||||
|
||||
|
||||
Foam::PackedBoolList Foam::syncTools::getInternalOrCoupledFaces
|
||||
(
|
||||
const polyMesh& mesh
|
||||
)
|
||||
{
|
||||
PackedBoolList isMasterFace(mesh.nFaces(), 1);
|
||||
|
||||
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
||||
|
||||
forAll(patches, patchI)
|
||||
{
|
||||
const polyPatch& pp = patches[patchI];
|
||||
|
||||
if (!pp.coupled())
|
||||
{
|
||||
forAll(pp, i)
|
||||
{
|
||||
isMasterFace.unset(pp.start()+i);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return isMasterFace;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -579,15 +579,25 @@ public:
|
||||
|
||||
// Other
|
||||
|
||||
//- Get per point whether is it master (of a coupled set of points)
|
||||
//- Get per point whether it is uncoupled or a master of a
|
||||
// coupled set of points
|
||||
static PackedBoolList getMasterPoints(const polyMesh&);
|
||||
|
||||
//- Get per edge whether is it master (of a coupled set of edges)
|
||||
//- Get per edge whether it is uncoupled or a master of a
|
||||
// coupled set of edges
|
||||
static PackedBoolList getMasterEdges(const polyMesh&);
|
||||
|
||||
//- Get per face whether is it master (of a coupled set of faces)
|
||||
//- Get per face whether it is uncoupled or a master of a
|
||||
// coupled set of faces
|
||||
static PackedBoolList getMasterFaces(const polyMesh&);
|
||||
|
||||
//- Get per face whether it is internal or a master of a
|
||||
// coupled set of faces
|
||||
static PackedBoolList getInternalOrMasterFaces(const polyMesh&);
|
||||
|
||||
//- Get per face whether it is internal or coupled
|
||||
static PackedBoolList getInternalOrCoupledFaces(const polyMesh&);
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -43,9 +43,6 @@ SourceFiles
|
||||
primitiveMeshEdges.C
|
||||
primitiveMeshCellCentresAndVols.C
|
||||
primitiveMeshFaceCentresAndAreas.C
|
||||
primitiveMeshEdgeVectors.C
|
||||
primitiveMeshCheck.C
|
||||
primitiveMeshCheckMotion.C
|
||||
primitiveMeshFindCell.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
@ -69,6 +66,8 @@ SourceFiles
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
class PackedBoolList;
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class primitiveMesh Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
@ -225,6 +224,28 @@ class primitiveMesh
|
||||
const labelList&
|
||||
);
|
||||
|
||||
protected:
|
||||
|
||||
// Static data members
|
||||
|
||||
//- Static data to control mesh checking
|
||||
|
||||
//- Cell closedness warning threshold
|
||||
// set as the fraction of un-closed area to closed area
|
||||
static scalar closedThreshold_;
|
||||
|
||||
//- Aspect ratio warning threshold
|
||||
static scalar aspectThreshold_;
|
||||
|
||||
//- Non-orthogonality warning threshold in deg
|
||||
static scalar nonOrthThreshold_;
|
||||
|
||||
//- Skewness warning threshold
|
||||
static scalar skewThreshold_;
|
||||
|
||||
//- Threshold where faces are considered coplanar
|
||||
static scalar planarCosAngle_;
|
||||
|
||||
|
||||
// Geometrical calculations
|
||||
|
||||
@ -251,7 +272,7 @@ class primitiveMesh
|
||||
void calcEdgeVectors() const;
|
||||
|
||||
|
||||
// Helper functions for mesh checking
|
||||
// Mesh checking
|
||||
|
||||
//- Check if all points on face are shared with another face.
|
||||
bool checkDuplicateFaces
|
||||
@ -270,29 +291,104 @@ class primitiveMesh
|
||||
labelHashSet*
|
||||
) const;
|
||||
|
||||
//- Check boundary for closedness
|
||||
bool checkClosedBoundary
|
||||
(
|
||||
const vectorField&,
|
||||
const bool,
|
||||
const PackedBoolList&
|
||||
) const;
|
||||
|
||||
// Static data members
|
||||
//- Check cells for closedness
|
||||
bool checkClosedCells
|
||||
(
|
||||
const vectorField& faceAreas,
|
||||
const scalarField& cellVolumes,
|
||||
const bool report,
|
||||
labelHashSet* setPtr,
|
||||
labelHashSet* aspectSetPtr,
|
||||
const Vector<label>& meshD
|
||||
) const;
|
||||
|
||||
//- Static data to control mesh checking
|
||||
//- Check for negative face areas
|
||||
bool checkFaceAreas
|
||||
(
|
||||
const vectorField& faceAreas,
|
||||
const bool report,
|
||||
const bool detailedReport,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Cell closedness warning threshold
|
||||
// set as the fraction of un-closed area to closed area
|
||||
static scalar closedThreshold_;
|
||||
//- Check for negative cell volumes
|
||||
bool checkCellVolumes
|
||||
(
|
||||
const scalarField& vols,
|
||||
const bool report,
|
||||
const bool detailedReport,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Aspect ratio warning threshold
|
||||
static scalar aspectThreshold_;
|
||||
//- Check for non-orthogonality
|
||||
bool checkFaceOrthogonality
|
||||
(
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs,
|
||||
const bool report,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Non-orthogonality warning threshold in deg
|
||||
static scalar nonOrthThreshold_;
|
||||
//- Check face pyramid volume
|
||||
bool checkFacePyramids
|
||||
(
|
||||
const pointField& points,
|
||||
const vectorField& ctrs,
|
||||
const bool report,
|
||||
const bool detailedReport,
|
||||
const scalar minPyrVol,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Skewness warning threshold
|
||||
static scalar skewThreshold_;
|
||||
//- Check face skewness
|
||||
bool checkFaceSkewness
|
||||
(
|
||||
const pointField& points,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs,
|
||||
const bool report,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Threshold where faces are considered coplanar
|
||||
static scalar planarCosAngle_;
|
||||
//- Check face angles
|
||||
bool checkFaceAngles
|
||||
(
|
||||
const pointField& points,
|
||||
const vectorField& faceAreas,
|
||||
const bool report,
|
||||
const scalar maxDeg,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Check face warpage
|
||||
bool checkFaceFlatness
|
||||
(
|
||||
const pointField& points,
|
||||
const vectorField& faceCentres,
|
||||
const vectorField& faceAreas,
|
||||
const bool report,
|
||||
const scalar warnFlatness,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Check for concave cells by the planes of faces
|
||||
bool checkConcaveCells
|
||||
(
|
||||
const vectorField& fAreas,
|
||||
const pointField& fCentres,
|
||||
const bool report,
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
protected:
|
||||
|
||||
//- Construct null
|
||||
primitiveMesh();
|
||||
@ -502,29 +598,36 @@ public:
|
||||
|
||||
// Topological checks
|
||||
|
||||
//- Check face ordering
|
||||
virtual bool checkUpperTriangular
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check cell zip-up
|
||||
bool checkCellsZipUp
|
||||
virtual bool checkCellsZipUp
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check uniqueness of face vertices
|
||||
bool checkFaceVertices
|
||||
virtual bool checkFaceVertices
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check for unused points
|
||||
virtual bool checkPoints
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check face-face connectivity
|
||||
bool checkFaceFaces
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check face ordering
|
||||
bool checkUpperTriangular
|
||||
virtual bool checkFaceFaces
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
@ -534,10 +637,11 @@ public:
|
||||
// Geometric checks
|
||||
|
||||
//- Check boundary for closedness
|
||||
bool checkClosedBoundary(const bool report = false) const;
|
||||
virtual bool checkClosedBoundary(const bool report = false)
|
||||
const;
|
||||
|
||||
//- Check cells for closedness
|
||||
bool checkClosedCells
|
||||
virtual bool checkClosedCells
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL,
|
||||
@ -546,28 +650,28 @@ public:
|
||||
) const;
|
||||
|
||||
//- Check for negative face areas
|
||||
bool checkFaceAreas
|
||||
virtual bool checkFaceAreas
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check for negative cell volumes
|
||||
bool checkCellVolumes
|
||||
virtual bool checkCellVolumes
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check for non-orthogonality
|
||||
bool checkFaceOrthogonality
|
||||
virtual bool checkFaceOrthogonality
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check face pyramid volume
|
||||
bool checkFacePyramids
|
||||
virtual bool checkFacePyramids
|
||||
(
|
||||
const bool report = false,
|
||||
const scalar minPyrVol = -SMALL,
|
||||
@ -575,14 +679,14 @@ public:
|
||||
) const;
|
||||
|
||||
//- Check face skewness
|
||||
bool checkFaceSkewness
|
||||
virtual bool checkFaceSkewness
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check face angles
|
||||
bool checkFaceAngles
|
||||
virtual bool checkFaceAngles
|
||||
(
|
||||
const bool report = false,
|
||||
const scalar maxSin = 10, // In degrees
|
||||
@ -592,31 +696,16 @@ public:
|
||||
//- Check face warpage: decompose face and check ratio between
|
||||
// magnitude of sum of triangle areas and sum of magnitude of
|
||||
// triangle areas.
|
||||
bool checkFaceFlatness
|
||||
virtual bool checkFaceFlatness
|
||||
(
|
||||
const bool report,
|
||||
const scalar warnFlatness, // When to include in set.
|
||||
labelHashSet* setPtr
|
||||
) const;
|
||||
|
||||
//- Check edge alignment for 1D/2D cases
|
||||
bool checkEdgeAlignment
|
||||
(
|
||||
const bool report,
|
||||
const Vector<label>& directions,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check for unused points
|
||||
bool checkPoints
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check for point-point-nearness,
|
||||
// e.g. colocated points which may be part of baffles.
|
||||
bool checkPointNearness
|
||||
virtual bool checkPointNearness
|
||||
(
|
||||
const bool report,
|
||||
const scalar reportDistSqr,
|
||||
@ -624,23 +713,15 @@ public:
|
||||
) const;
|
||||
|
||||
//- Check edge length
|
||||
bool checkEdgeLength
|
||||
virtual bool checkEdgeLength
|
||||
(
|
||||
const bool report,
|
||||
const scalar minLenSqr,
|
||||
labelHashSet* setPtr = NULL
|
||||
) const;
|
||||
|
||||
//- Check cell determinant
|
||||
bool checkCellDeterminant
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL,
|
||||
const Vector<label>& solutionD = Vector<label>::one
|
||||
) const;
|
||||
|
||||
//- Check for concave cells by the planes of faces
|
||||
bool checkConcaveCells
|
||||
virtual bool checkConcaveCells
|
||||
(
|
||||
const bool report = false,
|
||||
labelHashSet* setPtr = NULL
|
||||
@ -649,22 +730,14 @@ public:
|
||||
|
||||
//- Check mesh topology for correctness.
|
||||
// Returns false for no error.
|
||||
bool checkTopology(const bool report = false) const;
|
||||
virtual bool checkTopology(const bool report = false) const;
|
||||
|
||||
//- Check mesh geometry (& implicitly topology) for correctness.
|
||||
// Returns false for no error.
|
||||
bool checkGeometry(const bool report = false) const;
|
||||
virtual bool checkGeometry(const bool report = false) const;
|
||||
|
||||
//- Check mesh for correctness. Returns false for no error.
|
||||
bool checkMesh(const bool report = false) const;
|
||||
|
||||
//- Check mesh motion for correctness given motion points
|
||||
bool checkMeshMotion
|
||||
(
|
||||
const pointField& newPoints,
|
||||
const bool report = false
|
||||
) const;
|
||||
|
||||
virtual bool checkMesh(const bool report = false) const;
|
||||
|
||||
//- Set the closedness ratio warning threshold
|
||||
static scalar setClosedThreshold(const scalar);
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,266 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
|
||||
\\/ 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/>.
|
||||
|
||||
Description
|
||||
Given a set of points, find out if the mesh resulting from point motion will
|
||||
be valid without actually changing the mesh.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "primitiveMesh.H"
|
||||
#include "pyramidPointFaceRef.H"
|
||||
#include "cell.H"
|
||||
#include "unitConversion.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
bool Foam::primitiveMesh::checkMeshMotion
|
||||
(
|
||||
const pointField& newPoints,
|
||||
const bool report
|
||||
) const
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Pout<< "bool primitiveMesh::checkMeshMotion("
|
||||
<< "const pointField& newPoints, const bool report) const: "
|
||||
<< "checking mesh motion" << endl;
|
||||
}
|
||||
|
||||
bool error = false;
|
||||
|
||||
const faceList& f = faces();
|
||||
|
||||
const labelList& own = faceOwner();
|
||||
const labelList& nei = faceNeighbour();
|
||||
|
||||
vectorField fCtrs(nFaces());
|
||||
vectorField fAreas(nFaces());
|
||||
|
||||
makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
|
||||
|
||||
// Check cell volumes and calculate new cell centres
|
||||
vectorField cellCtrs(nCells());
|
||||
scalarField cellVols(nCells());
|
||||
|
||||
makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
|
||||
|
||||
scalar minVolume = GREAT;
|
||||
label nNegVols = 0;
|
||||
|
||||
forAll(cellVols, cellI)
|
||||
{
|
||||
if (cellVols[cellI] < VSMALL)
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Pout<< "Zero or negative cell volume detected for cell "
|
||||
<< cellI << ". Volume = " << cellVols[cellI] << endl;
|
||||
}
|
||||
|
||||
nNegVols++;
|
||||
}
|
||||
|
||||
minVolume = min(minVolume, cellVols[cellI]);
|
||||
}
|
||||
|
||||
if (nNegVols > 0)
|
||||
{
|
||||
error = true;
|
||||
|
||||
Pout<< "Zero or negative cell volume in mesh motion in " << nNegVols
|
||||
<< " cells. Min volume: " << minVolume << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Pout<< "Min volume = " << minVolume
|
||||
<< ". Total volume = " << sum(cellVols)
|
||||
<< ". Cell volumes OK." << endl;
|
||||
}
|
||||
}
|
||||
|
||||
// Check face areas
|
||||
|
||||
scalar minArea = GREAT;
|
||||
label nNegAreas = 0;
|
||||
label nPyrErrors = 0;
|
||||
label nDotProductErrors = 0;
|
||||
|
||||
forAll(f, faceI)
|
||||
{
|
||||
const scalar a = Foam::mag(fAreas[faceI]);
|
||||
|
||||
if (a < VSMALL)
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
if (isInternalFace(faceI))
|
||||
{
|
||||
Pout<< "Zero or negative face area detected for "
|
||||
<< "internal face "<< faceI << " between cells "
|
||||
<< own[faceI] << " and " << nei[faceI]
|
||||
<< ". Face area magnitude = " << a << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
Pout<< "Zero or negative face area detected for "
|
||||
<< "boundary face " << faceI << " next to cell "
|
||||
<< own[faceI] << ". Face area magnitude = "
|
||||
<< a << endl;
|
||||
}
|
||||
}
|
||||
|
||||
nNegAreas++;
|
||||
}
|
||||
|
||||
minArea = min(minArea, a);
|
||||
|
||||
// Create the owner pyramid - it will have negative volume
|
||||
scalar pyrVol =
|
||||
pyramidPointFaceRef(f[faceI], cellCtrs[own[faceI]]).mag(newPoints);
|
||||
|
||||
if (pyrVol > SMALL)
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Pout<< "Negative pyramid volume: " << -pyrVol
|
||||
<< " for face " << faceI << " " << f[faceI]
|
||||
<< " and owner cell: " << own[faceI] << endl
|
||||
<< "Owner cell vertex labels: "
|
||||
<< cells()[own[faceI]].labels(f)
|
||||
<< endl;
|
||||
}
|
||||
|
||||
nPyrErrors++;
|
||||
}
|
||||
|
||||
if (isInternalFace(faceI))
|
||||
{
|
||||
// Create the neighbour pyramid - it will have positive volume
|
||||
scalar pyrVol =
|
||||
pyramidPointFaceRef
|
||||
(
|
||||
f[faceI],
|
||||
cellCtrs[nei[faceI]]
|
||||
).mag(newPoints);
|
||||
|
||||
if (pyrVol < -SMALL)
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Pout<< "Negative pyramid volume: " << pyrVol
|
||||
<< " for face " << faceI << " " << f[faceI]
|
||||
<< " and neighbour cell: " << nei[faceI] << nl
|
||||
<< "Neighbour cell vertex labels: "
|
||||
<< cells()[nei[faceI]].labels(f)
|
||||
<< endl;
|
||||
}
|
||||
|
||||
nPyrErrors++;
|
||||
}
|
||||
|
||||
const vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
|
||||
const vector& s = fAreas[faceI];
|
||||
scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
|
||||
|
||||
// Only write full message the first time
|
||||
if (dDotS < SMALL && nDotProductErrors == 0)
|
||||
{
|
||||
// Non-orthogonality greater than 90 deg
|
||||
WarningIn
|
||||
(
|
||||
"primitiveMesh::checkMeshMotion"
|
||||
"(const pointField& newPoints, const bool report) const"
|
||||
) << "Severe non-orthogonality in mesh motion for face "
|
||||
<< faceI
|
||||
<< " between cells " << own[faceI] << " and " << nei[faceI]
|
||||
<< ": Angle = " << radToDeg(::acos(dDotS))
|
||||
<< " deg." << endl;
|
||||
|
||||
nDotProductErrors++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (nNegAreas > 0)
|
||||
{
|
||||
error = true;
|
||||
|
||||
WarningIn
|
||||
(
|
||||
"primitiveMesh::checkMeshMotion"
|
||||
"(const pointField& newPoints, const bool report) const"
|
||||
) << "Zero or negative face area in mesh motion in " << nNegAreas
|
||||
<< " faces. Min area: " << minArea << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Pout<< "Min area = " << minArea
|
||||
<< ". Face areas OK." << endl;
|
||||
}
|
||||
}
|
||||
|
||||
if (nPyrErrors > 0)
|
||||
{
|
||||
Pout<< "Detected " << nPyrErrors
|
||||
<< " negative pyramid volume in mesh motion" << endl;
|
||||
|
||||
error = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Pout<< "Pyramid volumes OK." << endl;
|
||||
}
|
||||
}
|
||||
|
||||
if (nDotProductErrors > 0)
|
||||
{
|
||||
Pout<< "Detected " << nDotProductErrors
|
||||
<< " in non-orthogonality in mesh motion." << endl;
|
||||
|
||||
error = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (debug || report)
|
||||
{
|
||||
Pout<< "Non-orthogonality check OK." << endl;
|
||||
}
|
||||
}
|
||||
|
||||
if (!error && (debug || report))
|
||||
{
|
||||
Pout<< "Mesh motion check OK." << endl;
|
||||
}
|
||||
|
||||
return error;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,638 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
|
||||
\\/ 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 "primitiveMeshTools.H"
|
||||
#include "syncTools.H"
|
||||
#include "pyramidPointFaceRef.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
Foam::scalar Foam::primitiveMeshTools::faceSkewness
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& p,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& fAreas,
|
||||
|
||||
const label faceI,
|
||||
const point& ownCc,
|
||||
const point& neiCc
|
||||
)
|
||||
{
|
||||
vector Cpf = fCtrs[faceI] - ownCc;
|
||||
vector d = neiCc - ownCc;
|
||||
|
||||
// Skewness vector
|
||||
vector sv =
|
||||
Cpf
|
||||
- ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + SMALL))*d;
|
||||
vector svHat = sv/(mag(sv) + VSMALL);
|
||||
|
||||
// Normalisation distance calculated as the approximate distance
|
||||
// from the face centre to the edge of the face in the direction
|
||||
// of the skewness
|
||||
scalar fd = 0.2*mag(d) + VSMALL;
|
||||
const face& f = mesh.faces()[faceI];
|
||||
forAll(f, pi)
|
||||
{
|
||||
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
|
||||
}
|
||||
|
||||
// Normalised skewness
|
||||
return mag(sv)/fd;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::primitiveMeshTools::faceOrthogonality
|
||||
(
|
||||
const point& ownCc,
|
||||
const point& neiCc,
|
||||
const vector& s
|
||||
)
|
||||
{
|
||||
vector d = neiCc - ownCc;
|
||||
|
||||
return (d & s)/(mag(d)*mag(s) + VSMALL);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceOrthogonality
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const vectorField& areas,
|
||||
const vectorField& cc
|
||||
)
|
||||
{
|
||||
const labelList& own = mesh.faceOwner();
|
||||
const labelList& nei = mesh.faceNeighbour();
|
||||
|
||||
tmp<scalarField> tortho(new scalarField(mesh.nInternalFaces()));
|
||||
scalarField& ortho = tortho();
|
||||
|
||||
// Internal faces
|
||||
forAll(nei, faceI)
|
||||
{
|
||||
ortho[faceI] = faceOrthogonality
|
||||
(
|
||||
cc[own[faceI]],
|
||||
cc[nei[faceI]],
|
||||
areas[faceI]
|
||||
);
|
||||
}
|
||||
|
||||
return tortho;
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& p,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs
|
||||
)
|
||||
{
|
||||
const labelList& own = mesh.faceOwner();
|
||||
const labelList& nei = mesh.faceNeighbour();
|
||||
const faceList& fcs = mesh.faces();
|
||||
|
||||
tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
|
||||
scalarField& skew = tskew();
|
||||
|
||||
forAll(nei, faceI)
|
||||
{
|
||||
skew[faceI] = faceSkewness
|
||||
(
|
||||
mesh,
|
||||
p,
|
||||
fCtrs,
|
||||
fAreas,
|
||||
|
||||
faceI,
|
||||
cellCtrs[own[faceI]],
|
||||
cellCtrs[nei[faceI]]
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// Boundary faces: consider them to have only skewness error.
|
||||
// (i.e. treat as if mirror cell on other side)
|
||||
|
||||
for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
|
||||
{
|
||||
vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
|
||||
|
||||
vector normal = fAreas[faceI];
|
||||
normal /= mag(normal) + VSMALL;
|
||||
vector d = normal*(normal & Cpf);
|
||||
|
||||
|
||||
// Skewness vector
|
||||
vector sv =
|
||||
Cpf
|
||||
- ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
|
||||
vector svHat = sv/(mag(sv) + VSMALL);
|
||||
|
||||
// Normalisation distance calculated as the approximate distance
|
||||
// from the face centre to the edge of the face in the direction
|
||||
// of the skewness
|
||||
scalar fd = 0.4*mag(d) + VSMALL;
|
||||
const face& f = fcs[faceI];
|
||||
forAll(f, pi)
|
||||
{
|
||||
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
|
||||
}
|
||||
|
||||
// Normalised skewness
|
||||
skew[faceI] = mag(sv)/fd;
|
||||
}
|
||||
|
||||
return tskew;
|
||||
}
|
||||
|
||||
void Foam::primitiveMeshTools::facePyramidVolume
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& points,
|
||||
const vectorField& ctrs,
|
||||
|
||||
scalarField& ownPyrVol,
|
||||
scalarField& neiPyrVol
|
||||
)
|
||||
{
|
||||
const labelList& own = mesh.faceOwner();
|
||||
const labelList& nei = mesh.faceNeighbour();
|
||||
const faceList& f = mesh.faces();
|
||||
|
||||
ownPyrVol.setSize(mesh.nFaces());
|
||||
neiPyrVol.setSize(mesh.nInternalFaces());
|
||||
|
||||
forAll(f, faceI)
|
||||
{
|
||||
// Create the owner pyramid
|
||||
ownPyrVol[faceI] = -pyramidPointFaceRef
|
||||
(
|
||||
f[faceI],
|
||||
ctrs[own[faceI]]
|
||||
).mag(points);
|
||||
|
||||
if (mesh.isInternalFace(faceI))
|
||||
{
|
||||
// Create the neighbour pyramid - it will have positive volume
|
||||
neiPyrVol[faceI] = pyramidPointFaceRef
|
||||
(
|
||||
f[faceI],
|
||||
ctrs[nei[faceI]]
|
||||
).mag(points);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::primitiveMeshTools::cellClosedness
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const Vector<label>& meshD,
|
||||
const vectorField& areas,
|
||||
const scalarField& vols,
|
||||
|
||||
scalarField& openness,
|
||||
scalarField& aratio
|
||||
)
|
||||
{
|
||||
const labelList& own = mesh.faceOwner();
|
||||
const labelList& nei = mesh.faceNeighbour();
|
||||
|
||||
// Loop through cell faces and sum up the face area vectors for each cell.
|
||||
// This should be zero in all vector components
|
||||
|
||||
vectorField sumClosed(mesh.nCells(), vector::zero);
|
||||
vectorField sumMagClosed(mesh.nCells(), vector::zero);
|
||||
|
||||
forAll(own, faceI)
|
||||
{
|
||||
// Add to owner
|
||||
sumClosed[own[faceI]] += areas[faceI];
|
||||
sumMagClosed[own[faceI]] += cmptMag(areas[faceI]);
|
||||
}
|
||||
|
||||
forAll(nei, faceI)
|
||||
{
|
||||
// Subtract from neighbour
|
||||
sumClosed[nei[faceI]] -= areas[faceI];
|
||||
sumMagClosed[nei[faceI]] += cmptMag(areas[faceI]);
|
||||
}
|
||||
|
||||
|
||||
label nDims = 0;
|
||||
for (direction dir = 0; dir < vector::nComponents; dir++)
|
||||
{
|
||||
if (meshD[dir] == 1)
|
||||
{
|
||||
nDims++;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Check the sums
|
||||
openness.setSize(mesh.nCells());
|
||||
aratio.setSize(mesh.nCells());
|
||||
|
||||
forAll(sumClosed, cellI)
|
||||
{
|
||||
scalar maxOpenness = 0;
|
||||
|
||||
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
|
||||
{
|
||||
maxOpenness = max
|
||||
(
|
||||
maxOpenness,
|
||||
mag(sumClosed[cellI][cmpt])
|
||||
/(sumMagClosed[cellI][cmpt] + VSMALL)
|
||||
);
|
||||
}
|
||||
openness[cellI] = maxOpenness;
|
||||
|
||||
// Calculate the aspect ration as the maximum of Cartesian component
|
||||
// aspect ratio to the total area hydraulic area aspect ratio
|
||||
scalar minCmpt = VGREAT;
|
||||
scalar maxCmpt = -VGREAT;
|
||||
for (direction dir = 0; dir < vector::nComponents; dir++)
|
||||
{
|
||||
if (meshD[dir] == 1)
|
||||
{
|
||||
minCmpt = min(minCmpt, sumMagClosed[cellI][dir]);
|
||||
maxCmpt = max(maxCmpt, sumMagClosed[cellI][dir]);
|
||||
}
|
||||
}
|
||||
|
||||
scalar aspectRatio = maxCmpt/(minCmpt + VSMALL);
|
||||
if (nDims == 3)
|
||||
{
|
||||
aspectRatio = max
|
||||
(
|
||||
aspectRatio,
|
||||
1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
|
||||
);
|
||||
}
|
||||
|
||||
aratio[cellI] = aspectRatio;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceConcavity
|
||||
(
|
||||
const scalar maxSin,
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& p,
|
||||
const vectorField& faceAreas
|
||||
)
|
||||
{
|
||||
const faceList& fcs = mesh.faces();
|
||||
|
||||
vectorField faceNormals(faceAreas);
|
||||
faceNormals /= mag(faceNormals) + VSMALL;
|
||||
|
||||
tmp<scalarField> tfaceAngles(new scalarField(mesh.nFaces()));
|
||||
scalarField& faceAngles = tfaceAngles();
|
||||
|
||||
|
||||
forAll(fcs, faceI)
|
||||
{
|
||||
const face& f = fcs[faceI];
|
||||
|
||||
// Get edge from f[0] to f[size-1];
|
||||
vector ePrev(p[f.first()] - p[f.last()]);
|
||||
scalar magEPrev = mag(ePrev);
|
||||
ePrev /= magEPrev + VSMALL;
|
||||
|
||||
scalar maxEdgeSin = 0.0;
|
||||
|
||||
forAll(f, fp0)
|
||||
{
|
||||
// Get vertex after fp
|
||||
label fp1 = f.fcIndex(fp0);
|
||||
|
||||
// Normalized vector between two consecutive points
|
||||
vector e10(p[f[fp1]] - p[f[fp0]]);
|
||||
scalar magE10 = mag(e10);
|
||||
e10 /= magE10 + VSMALL;
|
||||
|
||||
if (magEPrev > SMALL && magE10 > SMALL)
|
||||
{
|
||||
vector edgeNormal = ePrev ^ e10;
|
||||
scalar magEdgeNormal = mag(edgeNormal);
|
||||
|
||||
if (magEdgeNormal < maxSin)
|
||||
{
|
||||
// Edges (almost) aligned -> face is ok.
|
||||
}
|
||||
else
|
||||
{
|
||||
// Check normal
|
||||
edgeNormal /= magEdgeNormal;
|
||||
|
||||
if ((edgeNormal & faceNormals[faceI]) < SMALL)
|
||||
{
|
||||
maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
ePrev = e10;
|
||||
magEPrev = magE10;
|
||||
}
|
||||
|
||||
faceAngles[faceI] = maxEdgeSin;
|
||||
}
|
||||
|
||||
return tfaceAngles;
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceFlatness
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& p,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& faceAreas
|
||||
)
|
||||
{
|
||||
const faceList& fcs = mesh.faces();
|
||||
|
||||
// Areas are calculated as the sum of areas. (see
|
||||
// primitiveMeshFaceCentresAndAreas.C)
|
||||
scalarField magAreas(mag(faceAreas));
|
||||
|
||||
tmp<scalarField> tfaceFlatness(new scalarField(mesh.nFaces(), 1.0));
|
||||
scalarField& faceFlatness = tfaceFlatness();
|
||||
|
||||
|
||||
forAll(fcs, faceI)
|
||||
{
|
||||
const face& f = fcs[faceI];
|
||||
|
||||
if (f.size() > 3 && magAreas[faceI] > VSMALL)
|
||||
{
|
||||
const point& fc = fCtrs[faceI];
|
||||
|
||||
// Calculate the sum of magnitude of areas and compare to magnitude
|
||||
// of sum of areas.
|
||||
|
||||
scalar sumA = 0.0;
|
||||
|
||||
forAll(f, fp)
|
||||
{
|
||||
const point& thisPoint = p[f[fp]];
|
||||
const point& nextPoint = p[f.nextLabel(fp)];
|
||||
|
||||
// Triangle around fc.
|
||||
vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
|
||||
sumA += mag(n);
|
||||
}
|
||||
|
||||
faceFlatness[faceI] = magAreas[faceI] / (sumA+VSMALL);
|
||||
}
|
||||
}
|
||||
|
||||
return tfaceFlatness;
|
||||
}
|
||||
|
||||
|
||||
//Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::edgeAlignment
|
||||
//(
|
||||
// const primitiveMesh& mesh,
|
||||
// const Vector<label>& meshD,
|
||||
// const pointField& p
|
||||
//)
|
||||
//{
|
||||
// label nDirs = 0;
|
||||
// for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
|
||||
// {
|
||||
// if (meshD[cmpt] == 1)
|
||||
// {
|
||||
// nDirs++;
|
||||
// }
|
||||
// else if (meshD[cmpt] != 0)
|
||||
// {
|
||||
// FatalErrorIn
|
||||
// (
|
||||
// "primitiveMeshTools::edgeAlignment"
|
||||
// "(const primitiveMesh&, const Vector<label>&"
|
||||
// ", const pointField&)"
|
||||
// ) << "directions should contain 0 or 1 but is now " << meshD
|
||||
// << exit(FatalError);
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// tmp<scalarField> tedgeAlignment(new scalarField(mesh.nFaces(), 0.0));
|
||||
// scalarField& edgeAlignment = tedgeAlignment();
|
||||
//
|
||||
// if (nDirs != vector::nComponents)
|
||||
// {
|
||||
// const faceList& fcs = mesh.faces();
|
||||
//
|
||||
// forAll(fcs, faceI)
|
||||
// {
|
||||
// const face& f = fcs[faceI];
|
||||
//
|
||||
// forAll(f, fp)
|
||||
// {
|
||||
// label p0 = f[fp];
|
||||
// label p1 = f.nextLabel(fp);
|
||||
// if (p0 < p1)
|
||||
// {
|
||||
// vector d(p[p1]-p[p0]);
|
||||
// scalar magD = mag(d);
|
||||
//
|
||||
// if (magD > ROOTVSMALL)
|
||||
// {
|
||||
// d /= magD;
|
||||
//
|
||||
// // Check how many empty directions are used by the
|
||||
// // edge.
|
||||
// label nEmptyDirs = 0;
|
||||
// label nNonEmptyDirs = 0;
|
||||
// for
|
||||
// (
|
||||
// direction cmpt=0;
|
||||
// cmpt<vector::nComponents;
|
||||
// cmpt++
|
||||
// )
|
||||
// {
|
||||
// if (mag(d[cmpt]) > 1e-6)
|
||||
// {
|
||||
// if (meshD[cmpt] == 0)
|
||||
// {
|
||||
// nEmptyDirs++;
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// nNonEmptyDirs++;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// if (nEmptyDirs == 0)
|
||||
// {
|
||||
// // Purely in ok directions.
|
||||
// }
|
||||
// else if (nEmptyDirs == 1)
|
||||
// {
|
||||
// // Ok if purely in empty directions.
|
||||
// if (nNonEmptyDirs > 0)
|
||||
// {
|
||||
// edgeAlignment[faceI] = max
|
||||
// (
|
||||
// edgeAlignment[faceI],
|
||||
// 1
|
||||
// );
|
||||
// }
|
||||
// }
|
||||
// else if (nEmptyDirs > 1)
|
||||
// {
|
||||
// // Always an error
|
||||
// edgeAlignment[faceI] = max
|
||||
// (
|
||||
// edgeAlignment[faceI],
|
||||
// 2
|
||||
// );
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// return tedgeAlignment;
|
||||
//}
|
||||
|
||||
|
||||
// Checks cells with 1 or less internal faces. Give numerical problems.
|
||||
Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::cellDeterminant
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const Vector<label>& meshD,
|
||||
const vectorField& faceAreas,
|
||||
const PackedBoolList& internalOrCoupledFace
|
||||
)
|
||||
{
|
||||
// Determine number of dimensions and (for 2D) missing dimension
|
||||
label nDims = 0;
|
||||
label twoD = -1;
|
||||
for (direction dir = 0; dir < vector::nComponents; dir++)
|
||||
{
|
||||
if (meshD[dir] == 1)
|
||||
{
|
||||
nDims++;
|
||||
}
|
||||
else
|
||||
{
|
||||
twoD = dir;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
tmp<scalarField> tcellDeterminant(new scalarField(mesh.nCells()));
|
||||
scalarField& cellDeterminant = tcellDeterminant();
|
||||
|
||||
const cellList& c = mesh.cells();
|
||||
|
||||
if (nDims == 1)
|
||||
{
|
||||
cellDeterminant = 1.0;
|
||||
}
|
||||
else
|
||||
{
|
||||
forAll (c, cellI)
|
||||
{
|
||||
const labelList& curFaces = c[cellI];
|
||||
|
||||
// Calculate local normalization factor
|
||||
scalar avgArea = 0;
|
||||
|
||||
label nInternalFaces = 0;
|
||||
|
||||
forAll(curFaces, i)
|
||||
{
|
||||
if (internalOrCoupledFace[curFaces[i]])
|
||||
{
|
||||
avgArea += mag(faceAreas[curFaces[i]]);
|
||||
|
||||
nInternalFaces++;
|
||||
}
|
||||
}
|
||||
|
||||
if (nInternalFaces == 0)
|
||||
{
|
||||
cellDeterminant[cellI] = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
avgArea /= nInternalFaces;
|
||||
|
||||
symmTensor areaTensor(symmTensor::zero);
|
||||
|
||||
forAll(curFaces, i)
|
||||
{
|
||||
if (internalOrCoupledFace[curFaces[i]])
|
||||
{
|
||||
areaTensor += sqr(faceAreas[curFaces[i]]/avgArea);
|
||||
}
|
||||
}
|
||||
|
||||
if (nDims == 2)
|
||||
{
|
||||
// Add the missing eigenvector (such that it does not
|
||||
// affect the determinant)
|
||||
if (twoD == 0)
|
||||
{
|
||||
areaTensor.xx() = 1;
|
||||
}
|
||||
else if (twoD == 1)
|
||||
{
|
||||
areaTensor.yy() = 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
areaTensor.zz() = 1;
|
||||
}
|
||||
}
|
||||
|
||||
cellDeterminant[cellI] = mag(det(areaTensor));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return tcellDeterminant;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,168 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
|
||||
\\/ 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/>.
|
||||
|
||||
Namespace
|
||||
Foam::primitiveMeshTools
|
||||
|
||||
Description
|
||||
Collection of static functions operating on primitiveMesh (mainly checks).
|
||||
|
||||
SourceFiles
|
||||
primitiveMeshTools.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
#ifndef primitiveMeshTools_H
|
||||
#define primitiveMeshTools_H
|
||||
|
||||
#include "primitiveMesh.H"
|
||||
#include "PackedBoolList.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Namespace primitiveMeshTools Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class primitiveMeshTools
|
||||
{
|
||||
|
||||
protected:
|
||||
|
||||
// Protected Member Functions
|
||||
|
||||
//- Skewness of single face
|
||||
static scalar faceSkewness
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& p,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& fAreas,
|
||||
|
||||
const label faceI,
|
||||
const point& ownCc,
|
||||
const point& neiCc
|
||||
);
|
||||
|
||||
//- Orthogonality of single face
|
||||
static scalar faceOrthogonality
|
||||
(
|
||||
const point& ownCc,
|
||||
const point& neiCc,
|
||||
const vector& s
|
||||
);
|
||||
|
||||
public:
|
||||
|
||||
//- Generate non-orthogonality field (internal faces only)
|
||||
static tmp<scalarField> faceOrthogonality
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs
|
||||
);
|
||||
|
||||
//- Generate face pyramid volume fields
|
||||
static void facePyramidVolume
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& points,
|
||||
const vectorField& cellCtrs,
|
||||
scalarField& ownPyrVol,
|
||||
scalarField& neiPyrVol
|
||||
);
|
||||
|
||||
//- Generate skewness field
|
||||
static tmp<scalarField> faceSkewness
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& points,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& fAreas,
|
||||
const vectorField& cellCtrs
|
||||
);
|
||||
|
||||
//- Generate cell openness and cell ascpect ratio field
|
||||
static void cellClosedness
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const Vector<label>& meshD,
|
||||
const vectorField& areas,
|
||||
const scalarField& vols,
|
||||
scalarField& openness,
|
||||
scalarField& aratio
|
||||
);
|
||||
|
||||
//- Generate face concavity field. Returns per face the (sin of the)
|
||||
// most concave angle between two consecutive edges
|
||||
static tmp<scalarField> faceConcavity
|
||||
(
|
||||
const scalar maxSin,
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& p,
|
||||
const vectorField& faceAreas
|
||||
);
|
||||
|
||||
//- Generate face flatness field. Compares the individual triangles'
|
||||
// normals against the face average normal. Between 0 (fully warped)
|
||||
// and 1 (fully flat)
|
||||
static tmp<scalarField> faceFlatness
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const pointField& p,
|
||||
const vectorField& fCtrs,
|
||||
const vectorField& faceAreas
|
||||
);
|
||||
|
||||
//- Generate edge alignment field. Is per face the minimum aligned edge
|
||||
// (does not use edge addressing)
|
||||
static tmp<scalarField> edgeAlignment
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const Vector<label>& directions,
|
||||
const pointField& p
|
||||
);
|
||||
|
||||
//- Generate cell determinant field
|
||||
static tmp<scalarField> cellDeterminant
|
||||
(
|
||||
const primitiveMesh& mesh,
|
||||
const Vector<label>& directions,
|
||||
const vectorField& faceAreas,
|
||||
const PackedBoolList& internalOrCoupledFace
|
||||
);
|
||||
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
Loading…
Reference in New Issue
Block a user