zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-12-28 03:37:59 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge commit 'OpenCFD/master' into olesenm

Browse Source
This commit is contained in:
Mark Olesen
2008-10-09 16:13:41 +02:00
parent a60a620b80 b4263cca93
commit 7a77e7ad3e
36 changed files with 1877 additions and 1279 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
applications/utilities/mesh/manipulation/checkMesh.save/Make/files
View File

@ -1,6 +0,0 @@
printMeshStats.C
checkTopology.C
checkGeometry.C
checkMesh.C
EXE = $(FOAM_APPBIN)/checkMesh

5
applications/utilities/mesh/manipulation/checkMesh.save/Make/options
View File

@ -1,5 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lmeshTools

220
applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.C
View File

@ -1,220 +0,0 @@
#include "checkGeometry.H"
#include "polyMesh.H"
#include "globalMeshData.H"
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
Foam::label Foam::checkGeometry
(
const polyMesh& mesh,
bool checkPointNearness,
bool checkCellDeterminant
)
{
label noFailedChecks = 0;
Info<< "\nChecking geometry..." << endl;
boundBox bb(mesh.points());
Pout<< " Domain bounding box: "
<< bb.min() << " " << bb.max() << endl;
// Get a small relative length from the bounding box
const boundBox& globalBb = mesh.globalData().bb();
if (Pstream::parRun())
{
Info<< " Overall domain bounding box: "
<< globalBb.min() << " " << globalBb.max() << endl;
}
// Min length
scalar minDistSqr = magSqr(1e-6*(globalBb.max() - globalBb.min()));
if (mesh.checkClosedBoundary(true)) noFailedChecks++;
{
cellSet cells(mesh, "nonClosedCells", mesh.nCells()/100+1);
cellSet aspectCells(mesh, "highAspectRatioCells", mesh.nCells()/100+1);
if (mesh.checkClosedCells(true, &cells, &aspectCells))
{
noFailedChecks++;
if (cells.size() > 0)
{
Pout<< " <<Writing " << cells.size()
<< " non closed cells to set " << cells.name() << endl;
cells.write();
}
}
if (aspectCells.size() > 0)
{
Pout<< " <<Writing " << aspectCells.size()
<< " cells with high aspect ratio to set "
<< aspectCells.name() << endl;
aspectCells.write();
}
}
{
faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFaceAreas(true, &faces))
{
noFailedChecks++;
if (faces.size() > 0)
{
Pout<< " <<Writing " << faces.size()
<< " zero area faces to set " << faces.name() << endl;
faces.write();
}
}
}
{
cellSet cells(mesh, "zeroVolumeCells", mesh.nCells()/100 + 1);
if (mesh.checkCellVolumes(true, &cells))
{
noFailedChecks++;
if (cells.size() > 0)
{
Pout<< " <<Writing " << cells.size()
<< " zero volume cells to set " << cells.name() << endl;
cells.write();
}
}
}
{
faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFaceOrthogonality(true, &faces))
{
noFailedChecks++;
}
if (faces.size() > 0)
{
Pout<< " <<Writing " << faces.size()
<< " non-orthogonal faces to set " << faces.name() << endl;
faces.write();
}
}
{
faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFacePyramids(true, -SMALL, &faces))
{
noFailedChecks++;
if (faces.size() > 0)
{
Pout<< " <<Writing " << faces.size()
<< " faces with incorrect orientation to set "
<< faces.name() << endl;
faces.write();
}
}
}
{
faceSet faces(mesh, "skewFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFaceSkewness(true, &faces))
{
noFailedChecks++;
if (faces.size() > 0)
{
Pout<< " <<Writing " << faces.size()
<< " skew faces to set " << faces.name() << endl;
faces.write();
}
}
}
if (checkPointNearness)
{
// Note use of nPoints since don't want edge construction.
pointSet points(mesh, "shortEdges", mesh.nPoints()/1000 + 1);
if (mesh.checkEdgeLength(true, minDistSqr, &points))
{
//noFailedChecks++;
if (points.size() > 0)
{
Pout<< " <<Writing " << points.size()
<< " points on short edges to set " << points.name()
<< endl;
points.write();
}
}
label nEdgeClose = points.size();
if (mesh.checkPointNearness(false, minDistSqr, &points))
{
//noFailedChecks++;
if (points.size() > nEdgeClose)
{
pointSet nearPoints(mesh, "nearPoints", points);
Pout<< " <<Writing " << nearPoints.size()
<< " near (closer than " << Foam::sqrt(minDistSqr)
<< " apart) points to set " << nearPoints.name() << endl;
nearPoints.write();
}
}
}
{
faceSet faces(mesh, "concaveFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFaceAngles(true, 10, &faces))
{
//noFailedChecks++;
if (faces.size() > 0)
{
Pout<< " <<Writing " << faces.size()
<< " faces with concave angles to set " << faces.name()
<< endl;
faces.write();
}
}
}
{
faceSet faces(mesh, "warpedFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFaceFlatness(true, 0.8, &faces))
{
//noFailedChecks++;
if (faces.size() > 0)
{
Pout<< " <<Writing " << faces.size()
<< " warped faces to set " << faces.name() << endl;
faces.write();
}
}
}
if (checkCellDeterminant)
{
cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
if (mesh.checkCellDeterminant(true, &cells))
{
noFailedChecks++;
Pout<< " <<Writing " << cells.size()
<< " under-determines cells to set " << cells.name() << endl;
cells.write();
}
}
return noFailedChecks;
}

13
applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.H
View File

@ -1,13 +0,0 @@
#include "label.H"
namespace Foam
{
class polyMesh;
label checkGeometry
(
const polyMesh& mesh,
bool checkPointNearness,
bool checkCellDeterminant
);
}

152
applications/utilities/mesh/manipulation/checkMesh.save/checkMesh.C
View File

@ -1,152 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
checkMesh
Description
Checks validity of a mesh
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "polyMesh.H"
#include "globalMeshData.H"
#include "printMeshStats.H"
#include "checkTopology.H"
#include "checkGeometry.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "addTimeOptionsNoConstant.H"
argList::validOptions.insert("fullTopology", "");
argList::validOptions.insert("pointNearness", "");
argList::validOptions.insert("cellDeterminant", "");
# include "setRootCase.H"
# include "createTime.H"
// Get times list
instantList Times = runTime.times();
# include "checkTimeOptionsNoConstant.H"
runTime.setTime(Times[startTime], startTime);
# include "createPolyMesh.H"
bool firstCheck = true;
for (label i=startTime; i<endTime; i++)
{
runTime.setTime(Times[i], i);
polyMesh::readUpdateState state = mesh.readUpdate();
if
(
firstCheck
|| state == polyMesh::TOPO_CHANGE
|| state == polyMesh::TOPO_PATCH_CHANGE
)
{
firstCheck = false;
Info<< "Time = " << runTime.timeName() << nl << endl;
// Clear mesh before checking
mesh.clearOut();
// Reconstruct globalMeshData
mesh.globalData();
printMeshStats(mesh);
label noFailedChecks = 0;
noFailedChecks += checkTopology
(
mesh,
args.options().found("fullTopology")
);
noFailedChecks += checkGeometry
(
mesh,
args.options().found("pointNearness"),
args.options().found("cellDeterminant")
);
reduce(noFailedChecks, sumOp<label>());
if (noFailedChecks == 0)
{
Info<< "\nMesh OK."
<< nl << endl;
}
else
{
Info<< "\nFailed " << noFailedChecks << " mesh checks."
<< nl << endl;
}
}
else if (state == polyMesh::POINTS_MOVED)
{
label noFailedChecks = checkGeometry
(
mesh,
args.options().found("pointNearness"),
args.options().found("cellDeterminant")
);
reduce(noFailedChecks, sumOp<label>());
if (noFailedChecks == 0)
{
Info << "\nMesh OK."
<< nl << endl;
}
else
{
Info<< "\nFailed " << noFailedChecks << " mesh checks."
<< nl << endl;
}
}
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

236
applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.C
View File

@ -1,236 +0,0 @@
#include "checkTopology.H"
#include "polyMesh.H"
#include "Time.H"
#include "regionSplit.H"
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
#include "IOmanip.H"
Foam::label Foam::checkTopology(const polyMesh& mesh, bool fullTopology)
{
label noFailedChecks = 0;
Pout<< "Checking topology..." << endl;
// Check if the boundary definition is unique
mesh.boundaryMesh().checkDefinition(true);
// Check if the boundary processor patches are correct
mesh.boundaryMesh().checkParallelSync(true);
{
pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
if (mesh.checkPoints(true, &points))
{
noFailedChecks++;
Pout<< " <<Writing " << points.size()
<< " unused points to set " << points.name() << endl;
points.write();
}
}
{
faceSet faces(mesh, "upperTriangularFace", mesh.nFaces()/100);
if (mesh.checkUpperTriangular(true, &faces))
{
noFailedChecks++;
Pout<< " <<Writing " << faces.size()
<< " unordered faces to set " << faces.name() << endl;
faces.write();
}
}
{
cellSet cells(mesh, "zipUpCells", mesh.nCells()/100);
if (mesh.checkCellsZipUp(true, &cells))
{
noFailedChecks++;
Pout<< " <<Writing " << cells.size()
<< " cells with over used edges to set " << cells.name()
<< endl;
cells.write();
}
}
{
faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
if (mesh.checkFaceVertices(true, &faces))
{
noFailedChecks++;
Pout<< " <<Writing " << faces.size()
<< " faces with out-of-range vertices to set " << faces.name()
<< endl;
faces.write();
}
}
{
faceSet faces(mesh, "edgeFaces", mesh.nFaces()/100);
if (mesh.checkFaceFaces(true, &faces))
{
noFailedChecks++;
Pout<< " <<Writing " << faces.size()
<< " faces with incorrect edges to set " << faces.name()
<< endl;
faces.write();
}
}
{
regionSplit rs(mesh);
if (rs.nRegions() == 1)
{
Info<< " Number of regions: " << rs.nRegions() << " (OK)."
<< endl;
}
else
{
Info<< " *Number of regions: " << rs.nRegions() << endl;
Info<< " The mesh has multiple regions which are not connected "
"by any face." << endl
<< " <<Writing region information to "
<< mesh.time().timeName()/"cellToRegion"
<< endl;
labelIOList ctr
(
IOobject
(
"cellToRegion",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
rs
);
ctr.write();
}
}
if (!Pstream::parRun())
{
Pout<< "\nChecking patch topology for multiply connected surfaces ..."
<< endl;
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Non-manifold points
pointSet points
(
mesh,
"nonManifoldPoints",
mesh.nPoints()/100
);
Pout.setf(ios_base::left);
Pout<< " "
<< setw(20) << "Patch"
<< setw(9) << "Faces"
<< setw(9) << "Points"
<< " Surface" << endl;
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
primitivePatch::surfaceTopo pTyp = pp.surfaceType();
if (pp.size() == 0)
{
Pout<< " "
<< setw(20) << pp.name()
<< setw(9) << pp.size()
<< setw(9) << pp.nPoints()
<< " ok (empty)" << endl;
}
else if (pTyp == primitivePatch::MANIFOLD)
{
if (pp.checkPointManifold(true, &points))
{
Pout<< " "
<< setw(20) << pp.name()
<< setw(9) << pp.size()
<< setw(9) << pp.nPoints()
<< " multiply connected (shared point)" << endl;
}
else
{
Pout<< " "
<< setw(20) << pp.name()
<< setw(9) << pp.size()
<< setw(9) << pp.nPoints()
<< " ok (closed singly connected surface)" << endl;
}
// Add points on non-manifold edges to make set complete
pp.checkTopology(false, &points);
}
else
{
pp.checkTopology(false, &points);
if (pTyp == primitivePatch::OPEN)
{
Pout<< " "
<< setw(20) << pp.name()
<< setw(9) << pp.size()
<< setw(9) << pp.nPoints()
<< " ok (not multiply connected)" << endl;
}
else
{
Pout<< " "
<< setw(20) << pp.name()
<< setw(9) << pp.size()
<< setw(9) << pp.nPoints()
<< " multiply connected surface (shared edge)"
<< endl;
}
}
}
if (points.size() > 0)
{
Pout<< " <<Writing " << points.size()
<< " conflicting points to set "
<< points.name() << endl;
points.write();
}
//Pout.setf(ios_base::right);
}
// Force creation of all addressing if requested.
// Errors will be reported as required
if (fullTopology)
{
mesh.cells();
mesh.faces();
mesh.edges();
mesh.points();
mesh.faceOwner();
mesh.faceNeighbour();
mesh.cellCells();
mesh.edgeCells();
mesh.pointCells();
mesh.edgeFaces();
mesh.pointFaces();
mesh.cellEdges();
mesh.faceEdges();
mesh.pointEdges();
}
return noFailedChecks;
}

8
applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.H
View File

@ -1,8 +0,0 @@
#include "label.H"
namespace Foam
{
class polyMesh;
label checkTopology(const polyMesh& mesh, bool fullTopology);
}

95
applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.C
View File

@ -1,95 +0,0 @@
#include "printMeshStats.H"
#include "polyMesh.H"
#include "globalMeshData.H"
#include "hexMatcher.H"
#include "wedgeMatcher.H"
#include "prismMatcher.H"
#include "pyrMatcher.H"
#include "tetWedgeMatcher.H"
#include "tetMatcher.H"
void Foam::printMeshStats(const polyMesh& mesh)
{
Pout<< "Mesh stats" << nl
<< " points: " << mesh.points().size() << nl
<< " faces: " << mesh.faces().size() << nl
<< " internal faces: " << mesh.faceNeighbour().size() << nl
<< " cells: " << mesh.cells().size() << nl
<< " boundary patches: " << mesh.boundaryMesh().size() << nl
<< " point zones: " << mesh.pointZones().size() << nl
<< " face zones: " << mesh.faceZones().size() << nl
<< " cell zones: " << mesh.cellZones().size() << nl
<< endl;
if (Pstream::parRun())
{
const globalMeshData& parData = mesh.globalData();
Info<< "Overall stats" << nl
<< " points: " << parData.nTotalPoints() << nl
<< " faces: " << parData.nTotalFaces() << nl
<< " cells: " << parData.nTotalCells() << nl
<< endl;
}
// Construct shape recognizers
hexMatcher hex;
prismMatcher prism;
wedgeMatcher wedge;
pyrMatcher pyr;
tetWedgeMatcher tetWedge;
tetMatcher tet;
// Counters for different cell types
label nHex = 0;
label nWedge = 0;
label nPrism = 0;
label nPyr = 0;
label nTet = 0;
label nTetWedge = 0;
label nUnknown = 0;
for(label cellI = 0; cellI < mesh.nCells(); cellI++)
{
if (hex.isA(mesh, cellI))
{
nHex++;
}
else if (tet.isA(mesh, cellI))
{
nTet++;
}
else if (pyr.isA(mesh, cellI))
{
nPyr++;
}
else if (prism.isA(mesh, cellI))
{
nPrism++;
}
else if (wedge.isA(mesh, cellI))
{
nWedge++;
}
else if (tetWedge.isA(mesh, cellI))
{
nTetWedge++;
}
else
{
nUnknown++;
}
}
Pout<< "Number of cells of each type: " << nl
<< " hexahedra: " << nHex << nl
<< " prisms: " << nPrism << nl
<< " wedges: " << nWedge << nl
<< " pyramids: " << nPyr << nl
<< " tet wedges: " << nTetWedge << nl
<< " tetrahedra: " << nTet << nl
<< " polyhedra: " << nUnknown
<< nl << endl;
}

6
applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.H
View File

@ -1,6 +0,0 @@
namespace Foam
{
class polyMesh;
void printMeshStats(const polyMesh& mesh);
}

82
applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
View File

@ -39,10 +39,15 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
if (mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints))
{
noFailedChecks++;
label nNonAligned = returnReduce
(
nonAlignedPoints.size(),
sumOp<label>()
);
if (nonAlignedPoints.size() > 0)
if (nNonAligned > 0)
{
Pout<< " <<Writing " << nonAlignedPoints.size()
Info<< " <<Writing " << nNonAligned
<< " points on non-aligned edges to set "
<< nonAlignedPoints.name() << endl;
nonAlignedPoints.write();
@ -59,16 +64,21 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
noFailedChecks++;
if (cells.size() > 0)
label nNonClosed = returnReduce(cells.size(), sumOp<label>());
if (nNonClosed > 0)
{
Pout<< " <<Writing " << cells.size()
Info<< " <<Writing " << nNonClosed
<< " non closed cells to set " << cells.name() << endl;
cells.write();
}
}
if (aspectCells.size() > 0)
label nHighAspect = returnReduce(aspectCells.size(), sumOp<label>());
if (nHighAspect > 0)
{
Pout<< " <<Writing " << aspectCells.size()
Info<< " <<Writing " << nHighAspect
<< " cells with high aspect ratio to set "
<< aspectCells.name() << endl;
aspectCells.write();
@ -81,9 +91,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
noFailedChecks++;
if (faces.size() > 0)
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Pout<< " <<Writing " << faces.size()
Info<< " <<Writing " << nFaces
<< " zero area faces to set " << faces.name() << endl;
faces.write();
}
@ -96,9 +108,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
noFailedChecks++;
if (cells.size() > 0)
label nCells = returnReduce(cells.size(), sumOp<label>());
if (nCells > 0)
{
Pout<< " <<Writing " << cells.size()
Info<< " <<Writing " << nCells
<< " zero volume cells to set " << cells.name() << endl;
cells.write();
}
@ -112,9 +126,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
noFailedChecks++;
}
if (faces.size() > 0)
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Pout<< " <<Writing " << faces.size()
Info<< " <<Writing " << nFaces
<< " non-orthogonal faces to set " << faces.name() << endl;
faces.write();
}
@ -127,9 +143,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
noFailedChecks++;
if (faces.size() > 0)
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Pout<< " <<Writing " << faces.size()
Info<< " <<Writing " << nFaces
<< " faces with incorrect orientation to set "
<< faces.name() << endl;
faces.write();
@ -143,9 +161,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
noFailedChecks++;
if (faces.size() > 0)
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Pout<< " <<Writing " << faces.size()
Info<< " <<Writing " << nFaces
<< " skew faces to set " << faces.name() << endl;
faces.write();
}
@ -160,25 +180,29 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
//noFailedChecks++;
if (points.size() > 0)
label nPoints = returnReduce(points.size(), sumOp<label>());
if (nPoints > 0)
{
Pout<< " <<Writing " << points.size()
Info<< " <<Writing " << nPoints
<< " points on short edges to set " << points.name()
<< endl;
points.write();
}
}
label nEdgeClose = points.size();
label nEdgeClose = returnReduce(points.size(), sumOp<label>());
if (mesh.checkPointNearness(false, minDistSqr, &points))
{
//noFailedChecks++;
if (points.size() > nEdgeClose)
label nPoints = returnReduce(points.size(), sumOp<label>());
if (nPoints > nEdgeClose)
{
pointSet nearPoints(mesh, "nearPoints", points);
Pout<< " <<Writing " << nearPoints.size()
Info<< " <<Writing " << nPoints
<< " near (closer than " << Foam::sqrt(minDistSqr)
<< " apart) points to set " << nearPoints.name() << endl;
nearPoints.write();
@ -193,9 +217,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
//noFailedChecks++;
if (faces.size() > 0)
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Pout<< " <<Writing " << faces.size()
Info<< " <<Writing " << nFaces
<< " faces with concave angles to set " << faces.name()
<< endl;
faces.write();
@ -210,9 +236,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
//noFailedChecks++;
if (faces.size() > 0)
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Pout<< " <<Writing " << faces.size()
Info<< " <<Writing " << nFaces
<< " warped faces to set " << faces.name() << endl;
faces.write();
}
@ -226,7 +254,9 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
noFailedChecks++;
Pout<< " <<Writing " << cells.size()
label nCells = returnReduce(cells.size(), sumOp<label>());
Info<< " <<Writing " << nCells
<< " under-determined cells to set " << cells.name() << endl;
cells.write();
}

25
applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
View File

@ -16,7 +16,7 @@ Foam::label Foam::checkTopology
{
label noFailedChecks = 0;
Pout<< "Checking topology..." << endl;
Info<< "Checking topology..." << endl;
// Check if the boundary definition is unique
mesh.boundaryMesh().checkDefinition(true);
@ -30,7 +30,9 @@ Foam::label Foam::checkTopology
{
noFailedChecks++;
Pout<< " <<Writing " << points.size()
label nPoints = returnReduce(points.size(), sumOp<label>());
Info<< " <<Writing " << nPoints
<< " unused points to set " << points.name() << endl;
points.write();
}
@ -42,9 +44,12 @@ Foam::label Foam::checkTopology
{
noFailedChecks++;
}
if (faces.size() > 0)
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Pout<< " <<Writing " << faces.size()
Info<< " <<Writing " << nFaces
<< " unordered faces to set " << faces.name() << endl;
faces.write();
}
@ -57,7 +62,9 @@ Foam::label Foam::checkTopology
{
noFailedChecks++;
Pout<< " <<Writing " << cells.size()
label nCells = returnReduce(cells.size(), sumOp<label>());
Info<< " <<Writing " << nCells
<< " cells with over used edges to set " << cells.name()
<< endl;
cells.write();
@ -70,7 +77,9 @@ Foam::label Foam::checkTopology
{
noFailedChecks++;
Pout<< " <<Writing " << faces.size()
label nFaces = returnReduce(faces.size(), sumOp<label>());
Info<< " <<Writing " << nFaces
<< " faces with out-of-range or duplicate vertices to set "
<< faces.name() << endl;
faces.write();
@ -84,7 +93,9 @@ Foam::label Foam::checkTopology
{
noFailedChecks++;
Pout<< " <<Writing " << faces.size()
label nFaces = returnReduce(faces.size(), sumOp<label>());
Info<< " <<Writing " << nFaces
<< " faces with incorrect edges to set " << faces.name()
<< endl;
faces.write();

91
applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
View File

@ -12,43 +12,72 @@
void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
{
Pout<< "Mesh stats" << nl
<< " points: " << mesh.points().size() << nl;
Info<< "Mesh stats" << nl
<< " points: "
<< returnReduce(mesh.points().size(), sumOp<label>()) << nl;
if (mesh.nInternalPoints() != -1)
label nInternalPoints = returnReduce
(
mesh.nInternalPoints(),
sumOp<label>()
);
if (nInternalPoints != -Pstream::nProcs())
{
Pout<< " internal points: " << mesh.nInternalPoints() << nl;
Info<< " internal points: " << nInternalPoints << nl;
if (returnReduce(mesh.nInternalPoints(), minOp<label>()) == -1)
{
WarningIn("Foam::printMeshStats(const polyMesh&, const bool)")
<< "Some processors have their points sorted into internal"
<< " and external and some do not." << endl
<< "This can cause problems later on." << endl;
}
}
if (allTopology && mesh.nInternalPoints() != -1)
if (allTopology && nInternalPoints != -Pstream::nProcs())
{
Pout<< " edges: " << mesh.nEdges() << nl
<< " internal edges: " << mesh.nInternalEdges() << nl
label nEdges = returnReduce(mesh.nEdges(), sumOp<label>());
label nInternalEdges = returnReduce
(
mesh.nInternalEdges(),
sumOp<label>()
);
label nInternal1Edges = returnReduce
(
mesh.nInternal1Edges(),
sumOp<label>()
);
label nInternal0Edges = returnReduce
(
mesh.nInternal0Edges(),
sumOp<label>()
);
Info<< " edges: " << nEdges << nl
<< " internal edges: " << nInternalEdges << nl
<< " internal edges using one boundary point: "
<< mesh.nInternal1Edges()-mesh.nInternal0Edges() << nl
<< nInternal1Edges-nInternal0Edges << nl
<< " internal edges using two boundary points: "
<< mesh.nInternalEdges()-mesh.nInternal1Edges() << nl;
<< nInternalEdges-nInternal1Edges << nl;
}
Pout<< " faces: " << mesh.faces().size() << nl
<< " internal faces: " << mesh.faceNeighbour().size() << nl
<< " cells: " << mesh.cells().size() << nl
<< " boundary patches: " << mesh.boundaryMesh().size() << nl
<< " point zones: " << mesh.pointZones().size() << nl
<< " face zones: " << mesh.faceZones().size() << nl
<< " cell zones: " << mesh.cellZones().size() << nl
Info<< " faces: "
<< returnReduce(mesh.faces().size(), sumOp<label>()) << nl
<< " internal faces: "
<< returnReduce(mesh.faceNeighbour().size(), sumOp<label>()) << nl
<< " cells: "
<< returnReduce(mesh.cells().size(), sumOp<label>()) << nl
<< " boundary patches: "
<< returnReduce(mesh.boundaryMesh().size(), sumOp<label>()) << nl
<< " point zones: "
<< returnReduce(mesh.pointZones().size(), sumOp<label>()) << nl
<< " face zones: "
<< returnReduce(mesh.faceZones().size(), sumOp<label>()) << nl
<< " cell zones: "
<< returnReduce(mesh.cellZones().size(), sumOp<label>()) << nl
<< endl;
if (Pstream::parRun())
{
const globalMeshData& parData = mesh.globalData();
Info<< "Overall stats" << nl
<< " points: " << parData.nTotalPoints() << nl
<< " faces: " << parData.nTotalFaces() << nl
<< " cells: " << parData.nTotalCells() << nl
<< endl;
}
// Construct shape recognizers
hexMatcher hex;
@ -99,7 +128,15 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
}
}
Pout<< "Number of cells of each type: " << nl
reduce(nHex,sumOp<label>());
reduce(nPrism,sumOp<label>());
reduce(nWedge,sumOp<label>());
reduce(nPyr,sumOp<label>());
reduce(nTetWedge,sumOp<label>());
reduce(nTet,sumOp<label>());
reduce(nUnknown,sumOp<label>());
Info<< "Overall number of cells of each type:" << nl
<< " hexahedra: " << nHex << nl
<< " prisms: " << nPrism << nl
<< " wedges: " << nWedge << nl

375
applications/utilities/mesh/manipulation/createPatch/createPatch.C
View File

@ -34,6 +34,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "cyclicPolyPatch.H"
#include "syncTools.H"
#include "argList.H"
#include "polyMesh.H"
@ -256,27 +257,6 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
);
}
// cycPatch.writeOBJ
// (
// prefix+cycPatch.name()+"_half0.obj",
// SubList<face>
// (
// cycPatch,
// halfSize
// ),
// cycPatch.points()
// );
// cycPatch.writeOBJ
// (
// prefix+cycPatch.name()+"_half1.obj",
// SubList<face>
// (
// cycPatch,
// halfSize,
// halfSize
// ),
// cycPatch.points()
// );
// Lines between corresponding face centres
OFstream str(prefix+cycPatch.name()+"_match.obj");
@ -289,7 +269,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
vertI++;
label nbrFaceI = halfSize + faceI;
const point& fc1 = mesh.faceCentres()[cycPatch.start()+nbrFaceI];
const point& fc1 =
mesh.faceCentres()[cycPatch.start()+nbrFaceI];
meshTools::writeOBJ(str, fc1);
vertI++;
@ -300,6 +281,247 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
}
void separateList
(
const vectorField& separation,
UList<vector>& field
)
{
if (separation.size() == 1)
{
// Single value for all.
forAll(field, i)
{
field[i] += separation[0];
}
}
else if (separation.size() == field.size())
{
forAll(field, i)
{
field[i] += separation[i];
}
}
else
{
FatalErrorIn
(
"separateList(const vectorField&, UList<vector>&)"
) << "Sizes of field and transformation not equal. field:"
<< field.size() << " transformation:" << separation.size()
<< abort(FatalError);
}
}
// Synchronise points on both sides of coupled boundaries.
template <class CombineOp>
void syncPoints
(
const polyMesh& mesh,
pointField& points,
const CombineOp& cop,
const point& nullValue
)
{
if (points.size() != mesh.nPoints())
{
FatalErrorIn
(
"syncPoints"
"(const polyMesh&, pointField&, const CombineOp&, const point&)"
) << "Number of values " << points.size()
<< " is not equal to the number of points in the mesh "
<< mesh.nPoints() << abort(FatalError);
}
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Is there any coupled patch with transformation?
bool hasTransformation = false;
if (Pstream::parRun())
{
// Send
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if
(
isA<processorPolyPatch>(pp)
&& pp.nPoints() > 0
&& refCast<const processorPolyPatch>(pp).owner()
)
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
// Get data per patchPoint in neighbouring point numbers.
pointField patchInfo(procPatch.nPoints(), nullValue);
const labelList& meshPts = procPatch.meshPoints();
const labelList& nbrPts = procPatch.neighbPoints();
forAll(nbrPts, pointI)
{
label nbrPointI = nbrPts[pointI];
if (nbrPointI >= 0 && nbrPointI < patchInfo.size())
{
patchInfo[nbrPointI] = points[meshPts[pointI]];
}
}
OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
toNbr << patchInfo;
}
}
// Receive and set.
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if
(
isA<processorPolyPatch>(pp)
&& pp.nPoints() > 0
&& !refCast<const processorPolyPatch>(pp).owner()
)
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
pointField nbrPatchInfo(procPatch.nPoints());
{
// We do not know the number of points on the other side
// so cannot use Pstream::read.
IPstream fromNbr
(
Pstream::blocking,
procPatch.neighbProcNo()
);
fromNbr >> nbrPatchInfo;
}
// Null any value which is not on neighbouring processor
nbrPatchInfo.setSize(procPatch.nPoints(), nullValue);
if (!procPatch.parallel())
{
hasTransformation = true;
transformList(procPatch.forwardT(), nbrPatchInfo);
}
else if (procPatch.separated())
{
hasTransformation = true;
separateList(-procPatch.separation(), nbrPatchInfo);
}
const labelList& meshPts = procPatch.meshPoints();
forAll(meshPts, pointI)
{
label meshPointI = meshPts[pointI];
points[meshPointI] = nbrPatchInfo[pointI];
}
}
}
}
// Do the cyclics.
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (isA<cyclicPolyPatch>(pp))
{
const cyclicPolyPatch& cycPatch =
refCast<const cyclicPolyPatch>(pp);
const edgeList& coupledPoints = cycPatch.coupledPoints();
const labelList& meshPts = cycPatch.meshPoints();
pointField half0Values(coupledPoints.size());
forAll(coupledPoints, i)
{
const edge& e = coupledPoints[i];
label point0 = meshPts[e[0]];
half0Values[i] = points[point0];
}
if (!cycPatch.parallel())
{
hasTransformation = true;
transformList(cycPatch.reverseT(), half0Values);
}
else if (cycPatch.separated())
{
hasTransformation = true;
const vectorField& v = cycPatch.coupledPolyPatch::separation();
separateList(v, half0Values);
}
forAll(coupledPoints, i)
{
const edge& e = coupledPoints[i];
label point1 = meshPts[e[1]];
points[point1] = half0Values[i];
}
}
}
//- Note: hasTransformation is only used for warning messages so
// reduction not strictly nessecary.
//reduce(hasTransformation, orOp<bool>());
// Synchronize multiple shared points.
const globalMeshData& pd = mesh.globalData();
if (pd.nGlobalPoints() > 0)
{
if (hasTransformation)
{
WarningIn
(
"syncPoints"
"(const polyMesh&, pointField&, const CombineOp&, const point&)"
) << "There are decomposed cyclics in this mesh with"
<< " transformations." << endl
<< "This is not supported. The result will be incorrect"
<< endl;
}
// Values on shared points.
pointField sharedPts(pd.nGlobalPoints(), nullValue);
forAll(pd.sharedPointLabels(), i)
{
label meshPointI = pd.sharedPointLabels()[i];
// Fill my entries in the shared points
sharedPts[pd.sharedPointAddr()[i]] = points[meshPointI];
}
// Combine on master.
Pstream::listCombineGather(sharedPts, cop);
Pstream::listCombineScatter(sharedPts);
// Now we will all have the same information. Merge it back with
// my local information.
forAll(pd.sharedPointLabels(), i)
{
label meshPointI = pd.sharedPointLabels()[i];
points[meshPointI] = sharedPts[pd.sharedPointAddr()[i]];
}
}
}
// Main program:
int main(int argc, char *argv[])
@ -393,25 +615,28 @@ int main(int argc, char *argv[])
label destPatchI = patches.findPatchID(patchName);
word patchType(dict.lookup("type"));
if (destPatchI == -1)
{
dictionary patchDict(dict.subDict("dictionary"));
destPatchI = allPatches.size();
Info<< "Adding new patch " << patchName
<< " of type " << patchType
<< " as patch " << destPatchI << endl;
<< " as patch " << destPatchI
<< " from " << patchDict << endl;
patchDict.remove("nFaces");
patchDict.add("nFaces", 0);
patchDict.remove("startFace");
patchDict.add("startFace", startFaceI);
// Add an empty patch.
allPatches.append
(
polyPatch::New
(
patchType,
patchName,
0, // size
startFaceI, // start
patchDict,
destPatchI,
patches
).ptr()
@ -557,16 +782,100 @@ int main(int argc, char *argv[])
}
else
{
Info<< "Synchronising points." << nl << endl;
// This is a bit tricky. Both normal and position might be out and
// current separation also includes the normal
// ( separation_ = (nf&(Cr - Cf))*nf ).
// For processor patches:
// - disallow multiple separation/transformation. This basically
// excludes decomposed cyclics. Use the (probably 0) separation
// to align the points.
// For cyclic patches:
// - for separated ones use our own recalculated offset vector
// - for rotational ones use current one.
forAll(mesh.boundaryMesh(), patchI)
{
const polyPatch& pp = mesh.boundaryMesh()[patchI];
if (pp.size() && isA<coupledPolyPatch>(pp))
{
const coupledPolyPatch& cpp =
refCast<const coupledPolyPatch>(pp);
if (cpp.separated())
{
Info<< "On coupled patch " << pp.name()
<< " separation[0] was "
<< cpp.separation()[0] << endl;
if (isA<cyclicPolyPatch>(pp))
{
const cyclicPolyPatch& cycpp =
refCast<const cyclicPolyPatch>(pp);
if (cycpp.transform() == cyclicPolyPatch::TRANSLATIONAL)
{
Info<< "On cyclic translation patch " << pp.name()
<< " forcing uniform separation of "
<< cycpp.separationVector() << endl;
const_cast<vectorField&>(cpp.separation()) =
pointField(1, cycpp.separationVector());
}
else
{
const_cast<vectorField&>(cpp.separation()) =
pointField
(
1,
pp[pp.size()/2].centre(mesh.points())
- pp[0].centre(mesh.points())
);
}
}
else
{
const_cast<vectorField&>(cpp.separation())
.setSize(1);
}
Info<< "On coupled patch " << pp.name()
<< " forcing uniform separation of "
<< cpp.separation() << endl;
}
else if (!cpp.parallel())
{
Info<< "On coupled patch " << pp.name()
<< " forcing uniform rotation of "
<< cpp.forwardT()[0] << endl;
const_cast<tensorField&>
(
cpp.forwardT()
).setSize(1);
const_cast<tensorField&>
(
cpp.reverseT()
).setSize(1);
Info<< "On coupled patch " << pp.name()
<< " forcing uniform rotation of "
<< cpp.forwardT() << endl;
}
}
}
Info<< "Synchronising points." << endl;
pointField newPoints(mesh.points());
syncTools::syncPointList
syncPoints
(
mesh,
newPoints,
nearestEqOp(), // cop
point(GREAT, GREAT, GREAT), // nullValue
true // applySeparation
nearestEqOp(),
point(GREAT, GREAT, GREAT)
);
scalarField diff(mag(newPoints-mesh.points()));

63
applications/utilities/mesh/manipulation/createPatch/createPatchDict
View File

@ -1,58 +1,77 @@
/*--------------------------------*- C++ -*----------------------------------*\
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\ / O peration | Version: 1.0 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object createPatchDict;
version 2.0;
format ascii;
root "";
case "";
instance "system";
local "";
class dictionary;
object createPatcheDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Tolerance used in matching faces. Absolute tolerance is span of
// face times this factor.
matchTolerance 1E-3;
matchTolerance 1E-6;
// Do a synchronisation of coupled points.
pointSync true;
// Patches to create.
// If no patches does a coupled point and face synchronisation anyway.
patches
(
{
// Name of new patch
name leftRight0;
name sidePatches;
// Type of new patch
type cyclic;
// Dictionary for new patch
dictionary
{
type cyclic;
// Optional: used when matching and synchronising points.
//transform translational;
//separationVector (-2289 0 0);
}
// How to construct: either 'patches' or 'set'
constructFrom patches;
// If constructFrom = patches : names of patches
patches (half0 half1);
//patches (periodic-1 periodic-2);
patches (outlet-side1 outlet-side2);
// If constructFrom = set : name of faceSet
set f0;
}
{
name bottom;
type patch;
constructFrom set;
patches (half0 half1);
set bottomFaces;
}
//{
// name bottom;
// // Dictionary for new patch
// dictionary
// {
// type patch;
// }
//
// constructFrom set;
//
// patches (half0 half1);
//
// set bottomFaces;
//}
);