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

Merge branch 'master' of ssh://hunt/home/hunt2/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
henry
2008-10-09 14:18:23 +01:00
parent 45422cda59 12bec33d04
commit f86db703c3
31 changed files with 1697 additions and 1168 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C
View File

@ -157,9 +157,8 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
if(thermalCreep_) if(thermalCreep_)
{ {
const GeometricField<scalar, fvPatchField, volMesh>& vsfT = const volScalarField& vsfT =
this->db().objectRegistry:: this->db().objectRegistry::lookupObject<volScalarField>("T");
lookupObject<GeometricField<scalar, fvPatchField, volMesh> >("T");
label patchi = this->patch().index(); label patchi = this->patch().index();
const fvPatchScalarField& pT = vsfT.boundaryField()[patchi]; const fvPatchScalarField& pT = vsfT.boundaryField()[patchi];
Field<vector> gradpT = fvc::grad(vsfT)().boundaryField()[patchi]; Field<vector> gradpT = fvc::grad(vsfT)().boundaryField()[patchi];

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

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

@ -16,7 +16,7 @@ Foam::label Foam::checkTopology
{ {
label noFailedChecks = 0; label noFailedChecks = 0;
Pout<< "Checking topology..." << endl; Info<< "Checking topology..." << endl;
// Check if the boundary definition is unique // Check if the boundary definition is unique
mesh.boundaryMesh().checkDefinition(true); mesh.boundaryMesh().checkDefinition(true);
@ -30,7 +30,9 @@ Foam::label Foam::checkTopology
{ {
noFailedChecks++; noFailedChecks++;
Pout<< " <<Writing " << points.size() label nPoints = returnReduce(points.size(), sumOp<label>());
Info<< " <<Writing " << nPoints
<< " unused points to set " << points.name() << endl; << " unused points to set " << points.name() << endl;
points.write(); points.write();
} }
@ -42,9 +44,12 @@ Foam::label Foam::checkTopology
{ {
noFailedChecks++; 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; << " unordered faces to set " << faces.name() << endl;
faces.write(); faces.write();
} }
@ -57,7 +62,9 @@ Foam::label Foam::checkTopology
{ {
noFailedChecks++; noFailedChecks++;
Pout<< " <<Writing " << cells.size() label nCells = returnReduce(cells.size(), sumOp<label>());
Info<< " <<Writing " << nCells
<< " cells with over used edges to set " << cells.name() << " cells with over used edges to set " << cells.name()
<< endl; << endl;
cells.write(); cells.write();
@ -70,7 +77,9 @@ Foam::label Foam::checkTopology
{ {
noFailedChecks++; 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 with out-of-range or duplicate vertices to set "
<< faces.name() << endl; << faces.name() << endl;
faces.write(); faces.write();
@ -84,7 +93,9 @@ Foam::label Foam::checkTopology
{ {
noFailedChecks++; noFailedChecks++;
Pout<< " <<Writing " << faces.size() label nFaces = returnReduce(faces.size(), sumOp<label>());
Info<< " <<Writing " << nFaces
<< " faces with incorrect edges to set " << faces.name() << " faces with incorrect edges to set " << faces.name()
<< endl; << endl;
faces.write(); 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) void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
{ {
Pout<< "Mesh stats" << nl Info<< "Mesh stats" << nl
<< " points: " << mesh.points().size() << 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 label nEdges = returnReduce(mesh.nEdges(), sumOp<label>());
<< " internal edges: " << mesh.nInternalEdges() << nl 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: " << " internal edges using one boundary point: "
<< mesh.nInternal1Edges()-mesh.nInternal0Edges() << nl << nInternal1Edges-nInternal0Edges << nl
<< " internal edges using two boundary points: " << " internal edges using two boundary points: "
<< mesh.nInternalEdges()-mesh.nInternal1Edges() << nl; << nInternalEdges-nInternal1Edges << nl;
} }
Pout<< " faces: " << mesh.faces().size() << nl Info<< " faces: "
<< " internal faces: " << mesh.faceNeighbour().size() << nl << returnReduce(mesh.faces().size(), sumOp<label>()) << nl
<< " cells: " << mesh.cells().size() << nl << " internal faces: "
<< " boundary patches: " << mesh.boundaryMesh().size() << nl << returnReduce(mesh.faceNeighbour().size(), sumOp<label>()) << nl
<< " point zones: " << mesh.pointZones().size() << nl << " cells: "
<< " face zones: " << mesh.faceZones().size() << nl << returnReduce(mesh.cells().size(), sumOp<label>()) << nl
<< " cell zones: " << mesh.cellZones().size() << 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; << 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 // Construct shape recognizers
hexMatcher hex; 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 << " hexahedra: " << nHex << nl
<< " prisms: " << nPrism << nl << " prisms: " << nPrism << nl
<< " wedges: " << nWedge << 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 "syncTools.H"
#include "argList.H" #include "argList.H"
#include "polyMesh.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 // Lines between corresponding face centres
OFstream str(prefix+cycPatch.name()+"_match.obj"); OFstream str(prefix+cycPatch.name()+"_match.obj");
@ -289,7 +269,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
vertI++; vertI++;
label nbrFaceI = halfSize + faceI; label nbrFaceI = halfSize + faceI;
const point& fc1 = mesh.faceCentres()[cycPatch.start()+nbrFaceI]; const point& fc1 =
mesh.faceCentres()[cycPatch.start()+nbrFaceI];
meshTools::writeOBJ(str, fc1); meshTools::writeOBJ(str, fc1);
vertI++; 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: // Main program:
int main(int argc, char *argv[]) int main(int argc, char *argv[])
@ -393,25 +615,28 @@ int main(int argc, char *argv[])
label destPatchI = patches.findPatchID(patchName); label destPatchI = patches.findPatchID(patchName);
word patchType(dict.lookup("type"));
if (destPatchI == -1) if (destPatchI == -1)
{ {
dictionary patchDict(dict.subDict("dictionary"));
destPatchI = allPatches.size(); destPatchI = allPatches.size();
Info<< "Adding new patch " << patchName Info<< "Adding new patch " << patchName
<< " of type " << patchType << " as patch " << destPatchI
<< " as patch " << destPatchI << endl; << " from " << patchDict << endl;
patchDict.remove("nFaces");
patchDict.add("nFaces", 0);
patchDict.remove("startFace");
patchDict.add("startFace", startFaceI);
// Add an empty patch. // Add an empty patch.
allPatches.append allPatches.append
( (
polyPatch::New polyPatch::New
( (
patchType,
patchName, patchName,
0, // size patchDict,
startFaceI, // start
destPatchI, destPatchI,
patches patches
).ptr() ).ptr()
@ -557,16 +782,100 @@ int main(int argc, char *argv[])
} }
else 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; Info<< "Synchronising points." << endl;
pointField newPoints(mesh.points()); pointField newPoints(mesh.points());
syncTools::syncPointList
syncPoints
( (
mesh, mesh,
newPoints, newPoints,
nearestEqOp(), // cop nearestEqOp(),
point(GREAT, GREAT, GREAT), // nullValue point(GREAT, GREAT, GREAT)
true // applySeparation
); );
scalarField diff(mag(newPoints-mesh.points())); 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 | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 | | \\ / O peration | Version: 1.0 |
| \\ / A nd | Web: http://www.OpenFOAM.org | | \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | | | \\/ M anipulation | |
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
FoamFile FoamFile
{ {
version 2.0; version 2.0;
format ascii; format ascii;
class dictionary;
object createPatchDict; root "";
case "";
instance "system";
local "";
class dictionary;
object createPatcheDict;
} }
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Tolerance used in matching faces. Absolute tolerance is span of // Tolerance used in matching faces. Absolute tolerance is span of
// face times this factor. // face times this factor.
matchTolerance 1E-3; matchTolerance 1E-6;
// Do a synchronisation of coupled points. // Do a synchronisation of coupled points.
pointSync true; pointSync true;
// Patches to create. // Patches to create.
// If no patches does a coupled point and face synchronisation anyway. // If no patches does a coupled point and face synchronisation anyway.
patches patches
( (
{ {
// Name of new patch // Name of new patch
name leftRight0; name sidePatches;
// Type of new patch // Dictionary for new patch
type cyclic; dictionary
{
type cyclic;
// Optional: used when matching and synchronising points.
//transform translational;
//separationVector (-2289 0 0);
}
// How to construct: either 'patches' or 'set' // How to construct: either 'patches' or 'set'
constructFrom patches; constructFrom patches;
// If constructFrom = patches : names of 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 // If constructFrom = set : name of faceSet
set f0; set f0;
} }
{ //{
name bottom; // name bottom;
type patch; // // Dictionary for new patch
// dictionary
constructFrom set; // {
// type patch;
patches (half0 half1); // }
//
set bottomFaces; // constructFrom set;
} //
// patches (half0 half1);
//
// set bottomFaces;
//}
); );

2
applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
View File

@ -1,7 +1,9 @@
EXE_INC = \ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude -I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \ EXE_LIBS = \
-lmeshTools \
-lfiniteVolume \ -lfiniteVolume \
-lsampling -lsampling

381
applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
View File

@ -25,156 +25,269 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "channelIndex.H" #include "channelIndex.H"
#include "boolList.H"
#include "syncTools.H"
#include "OFstream.H"
#include "meshTools.H"
#include "Time.H"
#include "SortableList.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * * //
channelIndex::channelIndex(const fvMesh& m) template<>
: const char* Foam::NamedEnum<Foam::vector::components, 3>::names[] =
indexingDict_
(
IOobject
(
"postChannelDict",
m.time().constant(),
m,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
nx_(readLabel(indexingDict_.lookup("Nx"))),
ny_(indexingDict_.lookup("Ny")),
nz_(readLabel(indexingDict_.lookup("Nz"))),
symmetric_
(
readBool(indexingDict_.lookup("symmetric"))
),
cumNy_(ny_.size()),
nLayers_(ny_[0])
{ {
// initialise the layers "x",
cumNy_[0] = ny_[0]; "y",
"z"
};
for (label j=1; j<ny_.size(); j++) const Foam::NamedEnum<Foam::vector::components, 3>
Foam::channelIndex::vectorComponentsNames_;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Determines face blocking
void Foam::channelIndex::walkOppositeFaces
(
const polyMesh& mesh,
const labelList& startFaces,
boolList& blockedFace
)
{
const cellList& cells = mesh.cells();
const faceList& faces = mesh.faces();
label nBnd = mesh.nFaces() - mesh.nInternalFaces();
DynamicList<label> frontFaces(startFaces);
forAll(frontFaces, i)
{ {
nLayers_ += ny_[j]; label faceI = frontFaces[i];
cumNy_[j] = ny_[j]+cumNy_[j-1]; blockedFace[faceI] = true;
}
while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
{
// Transfer across.
boolList isFrontBndFace(nBnd, false);
forAll(frontFaces, i)
{
label faceI = frontFaces[i];
if (!mesh.isInternalFace(faceI))
{
isFrontBndFace[faceI-mesh.nInternalFaces()] = true;
}
}
syncTools::swapBoundaryFaceList(mesh, isFrontBndFace, false);
// Add
forAll(isFrontBndFace, i)
{
label faceI = mesh.nInternalFaces()+i;
if (isFrontBndFace[i] && !blockedFace[faceI])
{
blockedFace[faceI] = true;
frontFaces.append(faceI);
}
}
// Transfer across cells
DynamicList<label> newFrontFaces(frontFaces.size());
forAll(frontFaces, i)
{
label faceI = frontFaces[i];
{
const cell& ownCell = cells[mesh.faceOwner()[faceI]];
label oppositeFaceI = ownCell.opposingFaceLabel(faceI, faces);
if (oppositeFaceI == -1)
{
FatalErrorIn("channelIndex::walkOppositeFaces(..)")
<< "Face:" << faceI << " owner cell:" << ownCell
<< " is not a hex?" << abort(FatalError);
}
else
{
if (!blockedFace[oppositeFaceI])
{
blockedFace[oppositeFaceI] = true;
newFrontFaces.append(oppositeFaceI);
}
}
}
if (mesh.isInternalFace(faceI))
{
const cell& neiCell = mesh.cells()[mesh.faceNeighbour()[faceI]];
label oppositeFaceI = neiCell.opposingFaceLabel(faceI, faces);
if (oppositeFaceI == -1)
{
FatalErrorIn("channelIndex::walkOppositeFaces(..)")
<< "Face:" << faceI << " neighbour cell:" << neiCell
<< " is not a hex?" << abort(FatalError);
}
else
{
if (!blockedFace[oppositeFaceI])
{
blockedFace[oppositeFaceI] = true;
newFrontFaces.append(oppositeFaceI);
}
}
}
}
frontFaces.transfer(newFrontFaces);
} }
} }
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // Calculate regions.
void Foam::channelIndex::calcLayeredRegions
(
const polyMesh& mesh,
const labelList& startFaces
)
{
boolList blockedFace(mesh.nFaces(), false);
walkOppositeFaces
(
mesh,
startFaces,
blockedFace
);
channelIndex::~channelIndex()
{} if (false)
{
OFstream str(mesh.time().path()/"blockedFaces.obj");
label vertI = 0;
forAll(blockedFace, faceI)
{
if (blockedFace[faceI])
{
const face& f = mesh.faces()[faceI];
forAll(f, fp)
{
meshTools::writeOBJ(str, mesh.points()[f[fp]]);
}
str<< 'f';
forAll(f, fp)
{
str << ' ' << vertI+fp+1;
}
str << nl;
vertI += f.size();
}
}
}
// Do analysis for connected regions
cellRegion_.reset(new regionSplit(mesh, blockedFace));
Info<< "Detected " << cellRegion_().nRegions() << " layers." << nl << endl;
// Sum number of entries per region
regionCount_ = regionSum(scalarField(mesh.nCells(), 1.0));
// Average cell centres to determine ordering.
pointField regionCc
(
regionSum(mesh.cellCentres())
/ regionCount_
);
SortableList<scalar> sortComponent(regionCc.component(dir_));
sortMap_ = sortComponent.indices();
y_ = sortComponent;
if (symmetric_)
{
y_.setSize(cellRegion_().nRegions()/2);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::channelIndex::channelIndex
(
const polyMesh& mesh,
const dictionary& dict
)
:
symmetric_(readBool(dict.lookup("symmetric"))),
dir_(vectorComponentsNames_.read(dict.lookup("component")))
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
const wordList patchNames(dict.lookup("patches"));
label nFaces = 0;
forAll(patchNames, i)
{
label patchI = patches.findPatchID(patchNames[i]);
if (patchI == -1)
{
FatalErrorIn("channelIndex::channelIndex(const polyMesh&)")
<< "Illegal patch " << patchNames[i]
<< ". Valid patches are " << patches.name()
<< exit(FatalError);
}
nFaces += patches[patchI].size();
}
labelList startFaces(nFaces);
nFaces = 0;
forAll(patchNames, i)
{
const polyPatch& pp = patches[patches.findPatchID(patchNames[i])];
forAll(pp, j)
{
startFaces[nFaces++] = pp.start()+j;
}
}
// Calculate regions.
calcLayeredRegions(mesh, startFaces);
}
Foam::channelIndex::channelIndex
(
const polyMesh& mesh,
const labelList& startFaces,
const bool symmetric,
const direction dir
)
:
symmetric_(symmetric),
dir_(dir)
{
// Calculate regions.
calcLayeredRegions(mesh, startFaces);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalarField channelIndex::collapse
(
const volScalarField& vsf,
const bool asymmetric
) const
{
scalarField cs(nLayers(), 0.0);
forAll(cs, j)
{
// sweep over all cells in this layer
for (label i=0; i<nx(); i++)
{
for (label k=0; k<nz(); k++)
{
cs[j] += vsf[operator()(i,j,k)];
}
}
// and divide by the number of cells in the layer
cs[j] /= scalar(nx()*nz());
}
if (symmetric_)
{
label nlb2 = nLayers()/2;
if (asymmetric)
{
for (label j=0; j<nlb2; j++)
{
cs[j] = 0.5*(cs[j] - cs[nLayers() - j - 1]);
}
}
else
{
for (label j=0; j<nlb2; j++)
{
cs[j] = 0.5*(cs[j] + cs[nLayers() - j - 1]);
}
}
cs.setSize(nlb2);
}
return cs;
}
scalarField channelIndex::y
(
const volVectorField& cellCentres
) const
{
if (symmetric_)
{
scalarField Y(nLayers()/2);
for (label j=0; j<nLayers()/2; j++)
{
Y[j] = cellCentres[operator()(0, j, 0)].y();
}
return Y;
}
else
{
scalarField Y(nLayers());
for (label j=0; j<nLayers(); j++)
{
Y[j] = cellCentres[operator()(0, j, 0)].y();
}
return Y;
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
label channelIndex::operator()
(
const label Jx,
const label Jy,
const label Jz
) const
{
label index(0);
// count up `full' layers in the mesh
label j(0);
label tmpJy(Jy);
while(Jy >= cumNy_[j])
{
index += nx_*ny_[j]*nz_;
tmpJy -= ny_[j];
j++;
}
index += Jx + nx_*tmpJy + nx_*ny_[j]*Jz;
return index;
}
// ************************************************************************* // // ************************************************************************* //

114
applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
View File

@ -26,19 +26,25 @@ Class
Foam::channelIndex Foam::channelIndex
Description Description
does indexing for a standard (Christer-style) channel. Does averaging of fields over layers of cells. Assumes layered mesh.
Assumes that the blocks are aranged in the y direction.
SourceFiles SourceFiles
channelIndex.C channelIndex.C
channelIndexIO.C
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#ifndef channelIndex_H #ifndef channelIndex_H
#define channelIndex_H #define channelIndex_H
#include "fvCFD.H" #include "regionSplit.H"
#include "direction.H"
#include "scalarField.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
@ -47,21 +53,47 @@ SourceFiles
class channelIndex class channelIndex
{ {
// Private data // Private data
IOdictionary indexingDict_; static const NamedEnum<vector::components, 3> vectorComponentsNames_;
const label nx_; //- Is mesh symmetric
const labelList ny_;
const label nz_;
const bool symmetric_; const bool symmetric_;
labelList cumNy_; //- direction to sort
label nLayers_; const direction dir_;
//- Per cell the global region
autoPtr<regionSplit> cellRegion_;
//- Per global region the number of cells (scalarField so we can use
// field algebra)
scalarField regionCount_;
//- From sorted region back to unsorted global region
labelList sortMap_;
//- Sorted component of cell centres
scalarField y_;
// Private Member Functions // Private Member Functions
void walkOppositeFaces
(
const polyMesh& mesh,
const labelList& startFaces,
boolList& blockedFace
);
void calcLayeredRegions
(
const polyMesh& mesh,
const labelList& startFaces
);
//- Disallow default bitwise copy construct and assignment //- Disallow default bitwise copy construct and assignment
channelIndex(const channelIndex&); channelIndex(const channelIndex&);
void operator=(const channelIndex&); void operator=(const channelIndex&);
@ -71,56 +103,54 @@ public:
// Constructors // Constructors
channelIndex(const fvMesh& m); //- Construct from dictionary
channelIndex(const polyMesh&, const dictionary&);
//- Construct from supplied starting faces
// Destructor channelIndex
(
~channelIndex(); const polyMesh& mesh,
const labelList& startFaces,
const bool symmetric,
const direction dir
);
// Member Functions // Member Functions
// Access // Access
//- return number of layers //- Sum field per region
label nLayers() const template<class T>
{ Field<T> regionSum(const Field<T>& cellField) const;
return nLayers_;
}
//- number of cells in X direction
label nx() const
{
return nx_;
}
//- number of cells in Z direction
label nz() const
{
return nz_;
}
//- collapse a field to a line //- collapse a field to a line
scalarField collapse template<class T>
Field<T> collapse
( (
const volScalarField& vsf, const Field<T>& vsf,
const bool asymmetric=false const bool asymmetric=false
) const; ) const;
//- return the field of Y locations from the cell centres //- return the field of Y locations from the cell centres
scalarField y const scalarField& y() const
( {
const volVectorField& cellCentres return y_;
) const; }
// Member Operators
label operator()(const label, const label, const label) const;
}; };
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "channelIndexTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif #endif

101
applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C Normal file
View File

@ -0,0 +1,101 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "channelIndex.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class T>
Foam::Field<T> Foam::channelIndex::regionSum(const Field<T>& cellField) const
{
Field<T> regionField(cellRegion_().nRegions(), pTraits<T>::zero);
forAll(cellRegion_(), cellI)
{
regionField[cellRegion_()[cellI]] += cellField[cellI];
}
// Global sum
Pstream::listCombineGather(regionField, plusEqOp<T>());
Pstream::listCombineScatter(regionField);
return regionField;
}
template<class T>
Foam::Field<T> Foam::channelIndex::collapse
(
const Field<T>& cellField,
const bool asymmetric
) const
{
// Average and order
Field<T> regionField
(
regionSum(cellField)
/ regionCount_,
sortMap_
);
// Symmetry?
if (symmetric_)
{
label nlb2 = cellRegion_().nRegions()/2;
if (asymmetric)
{
for (label j=0; j<nlb2; j++)
{
regionField[j] =
0.5
* (
regionField[j]
- regionField[cellRegion_().nRegions() - j - 1]
);
}
}
else
{
for (label j=0; j<nlb2; j++)
{
regionField[j] =
0.5
* (
regionField[j]
+ regionField[cellRegion_().nRegions() - j - 1]
);
}
}
regionField.setSize(nlb2);
}
return regionField;
}
// ************************************************************************* //

8
applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
View File

@ -1,16 +1,16 @@
scalarField UMeanXvalues = channelIndexing.collapse scalarField UMeanXvalues = channelIndexing.collapse
( (
UMean.component(vector::X) UMean.component(vector::X)()
); );
scalarField UMeanYvalues = channelIndexing.collapse scalarField UMeanYvalues = channelIndexing.collapse
( (
UMean.component(vector::Y) UMean.component(vector::Y)()
); );
scalarField UMeanZvalues = channelIndexing.collapse scalarField UMeanZvalues = channelIndexing.collapse
( (
UMean.component(vector::Z) UMean.component(vector::Z)()
); );
scalarField RxxValues = channelIndexing.collapse(Rxx); scalarField RxxValues = channelIndexing.collapse(Rxx);
@ -38,7 +38,7 @@
0.5*(sqr(urmsValues) + sqr(vrmsValues) + sqr(wrmsValues)); 0.5*(sqr(urmsValues) + sqr(vrmsValues) + sqr(wrmsValues));
scalarField y = channelIndexing.y(mesh.C()); const scalarField& y = channelIndexing.y();
makeGraph(y, UMeanXvalues, "Uf", UMean.path(), gFormat); makeGraph(y, UMeanXvalues, "Uf", UMean.path(), gFormat);
makeGraph(y, urmsValues, "u", UMean.path(), gFormat); makeGraph(y, urmsValues, "u", UMean.path(), gFormat);

15
applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
View File

@ -67,7 +67,20 @@ int main(int argc, char *argv[])
const word& gFormat = runTime.graphFormat(); const word& gFormat = runTime.graphFormat();
// Setup channel indexing for averaging over channel down to a line // Setup channel indexing for averaging over channel down to a line
channelIndex channelIndexing(mesh);
IOdictionary channelDict
(
IOobject
(
"postChannelDict",
mesh.time().constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
channelIndex channelIndexing(mesh, channelDict);
// For each time step read all fields // For each time step read all fields
for (label i=startTime; i<endTime; i++) for (label i=startTime; i<endTime; i++)

14
applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict
View File

@ -15,14 +15,14 @@ FoamFile
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Nx 40; // Seed patches to start layering from
Ny patches (bottomWall);
(
25
25
);
Nz 30;
// Direction in which the layers are
component y;
// Is the mesh symmetric? If so average(symmetric fields) or
// subtract(asymmetric) contributions from both halves
symmetric true; symmetric true;
// ************************************************************************* // // ************************************************************************* //

49
applications/utilities/postProcessing/miscellaneous/postChannel/sumData.C Normal file
View File

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "sumData.H"
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<
(
Foam::Ostream& os,
const Foam::sumData& wDist
)
{
return os
<< wDist.oldFace_ << token::SPACE
<< wDist.sum_ << token::SPACE << wDist.count_;
}
Foam::Istream& Foam::operator>>(Foam::Istream& is, Foam::sumData& wDist)
{
return is >> wDist.oldFace_ >> wDist.sum_ >> wDist.count_;
}
// ************************************************************************* //

200
applications/utilities/postProcessing/miscellaneous/postChannel/sumData.H Normal file
View File

@ -0,0 +1,200 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
Class
Foam::sumData
Description
Sums data while walking across cells. Used in collapsing fields.
SourceFiles
sumDataI.H
sumData.C
\*---------------------------------------------------------------------------*/
#ifndef sumData_H
#define sumData_H
#include "point.H"
#include "tensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class polyPatch;
class polyMesh;
/*---------------------------------------------------------------------------*\
Class sumData Declaration
\*---------------------------------------------------------------------------*/
class sumData
{
// Private data
//- Previous face
label oldFace_;
//- summed data
scalar sum_;
//- number of items summed
label count_;
public:
// Constructors
//- Construct null
inline sumData();
//- Construct from count
inline sumData
(
const label oldFace,
const scalar sum,
const label count
);
// Member Functions
// Access
inline label oldFace() const
{
return oldFace_;
}
inline scalar sum() const
{
return sum_;
}
inline label count() const
{
return count_;
}
// Needed by FaceCellWave
//- Check whether origin has been changed at all or
// still contains original (invalid) value.
inline bool valid() const;
//- Check for identical geometrical data. Used for cyclics checking.
inline bool sameGeometry
(
const polyMesh&,
const sumData&,
const scalar
) const;
//- Convert any absolute coordinates into relative to (patch)face
// centre
inline void leaveDomain
(
const polyMesh&,
const polyPatch&,
const label patchFaceI,
const point& faceCentre
);
//- Reverse of leaveDomain
inline void enterDomain
(
const polyMesh&,
const polyPatch&,
const label patchFaceI,
const point& faceCentre
);
//- Apply rotation matrix to any coordinates
inline void transform
(
const polyMesh&,
const tensor&
);
//- Influence of neighbouring face.
inline bool updateCell
(
const polyMesh&,
const label thisCellI,
const label neighbourFaceI,
const sumData& neighbourInfo,
const scalar tol
);
//- Influence of neighbouring cell.
inline bool updateFace
(
const polyMesh&,
const label thisFaceI,
const label neighbourCellI,
const sumData& neighbourInfo,
const scalar tol
);
//- Influence of different value on same face.
inline bool updateFace
(
const polyMesh&,
const label thisFaceI,
const sumData& neighbourInfo,
const scalar tol
);
// Member Operators
// Needed for List IO
inline bool operator==(const sumData&) const;
inline bool operator!=(const sumData&) const;
// IOstream Operators
friend Ostream& operator<<(Ostream&, const sumData&);
friend Istream& operator>>(Istream&, sumData&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "sumDataI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

227
applications/utilities/postProcessing/miscellaneous/postChannel/sumDataI.H Normal file
View File

@ -0,0 +1,227 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "polyMesh.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Null constructor
inline Foam::sumData::sumData()
:
oldFace_(-1),
sum_(0.0),
count_(0)
{}
// Construct from components
inline Foam::sumData::sumData
(
const label oldFace,
const scalar sum,
const label count
)
:
oldFace_(oldFace),
sum_(sum),
count_(count)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline bool Foam::sumData::valid() const
{
return oldFace_ != -1;
}
// No geometric data so never any problem on cyclics
inline bool Foam::sumData::sameGeometry
(
const polyMesh&,
const sumData&,
const scalar
) const
{
return true;
}
// No geometric data.
inline void Foam::sumData::leaveDomain
(
const polyMesh&,
const polyPatch& patch,
const label patchFaceI,
const point& faceCentre
)
{}
// No geometric data.
inline void Foam::sumData::transform
(
const polyMesh&,
const tensor& rotTensor
)
{}
// No geometric data.
inline void Foam::sumData::enterDomain
(
const polyMesh&,
const polyPatch& patch,
const label patchFaceI,
const point& faceCentre
)
{
oldFace_ = -1;
}
// Update cell with neighbouring face information
inline bool Foam::sumData::updateCell
(
const polyMesh&,
const label thisCellI,
const label neighbourFaceI,
const sumData& neighbourInfo,
const scalar tol
)
{
if (!valid())
{
FatalErrorIn("sumData::updateCell") << "problem"
<< abort(FatalError);
return false;
}
if (count_ == 0)
{
sum_ += neighbourInfo.sum();
count_ = neighbourInfo.count() + 1;
oldFace_ = neighbourFaceI;
return true;
}
else
{
return false;
}
}
// Update face with neighbouring cell information
inline bool Foam::sumData::updateFace
(
const polyMesh& mesh,
const label thisFaceI,
const label neighbourCellI,
const sumData& neighbourInfo,
const scalar tol
)
{
// From cell to its faces.
// Check that face is opposite the previous one.
const cell& cFaces = mesh.cells()[neighbourCellI];
label wantedFaceI = cFaces.opposingFaceLabel
(
neighbourInfo.oldFace(),
mesh.faces()
);
if (thisFaceI == wantedFaceI)
{
if (count_ != 0)
{
FatalErrorIn("sumData::updateFace") << "problem"
<< abort(FatalError);
return false;
}
sum_ += neighbourInfo.sum();
count_ = neighbourInfo.count();
oldFace_ = thisFaceI;
return true;
}
else
{
return false;
}
}
// Update face with coupled face information
inline bool Foam::sumData::updateFace
(
const polyMesh&,
const label thisFaceI,
const sumData& neighbourInfo,
const scalar tol
)
{
// From face to face (e.g. coupled faces)
if (count_ == 0)
{
sum_ += neighbourInfo.sum();
count_ = neighbourInfo.count();
oldFace_ = thisFaceI;
return true;
}
else
{
return false;
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline bool Foam::sumData::operator==(const Foam::sumData& rhs)
const
{
return
oldFace() == rhs.oldFace()
&& sum() == rhs.sum()
&& count() == rhs.count();
}
inline bool Foam::sumData::operator!=(const Foam::sumData& rhs)
const
{
return !(*this == rhs);
}
// ************************************************************************* //

4
src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
View File

@ -552,6 +552,8 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
nextPatchStart += bm[patchI].size(); nextPatchStart += bm[patchI].size();
} }
reduce(boundaryError, orOp<bool>());
if (boundaryError) if (boundaryError)
{ {
if (debug || report) if (debug || report)
@ -565,7 +567,7 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
{ {
if (debug || report) if (debug || report)
{ {
Pout << " Boundary definition OK." << endl; Info << " Boundary definition OK." << endl;
} }
return false; return false;

10
src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
View File

@ -263,6 +263,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl
<< " (half)size:" << Cf.size() << nl << " (half)size:" << Cf.size() << nl
<< " absTol:" << absTol << nl << " absTol:" << absTol << nl
//<< " smallDist:" << smallDist << nl
<< " sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl; << " sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl;
} }
@ -316,6 +317,13 @@ void Foam::coupledPolyPatch::calcTransformTensors
{ {
forwardT_.setSize(1); forwardT_.setSize(1);
reverseT_.setSize(1); reverseT_.setSize(1);
if (debug)
{
Pout<< " rotation " << sum(mag(forwardT_ - forwardT_[0]))
<< " more than local tolerance " << error
<< ". Assuming uniform rotation." << endl;
}
} }
} }
else else
@ -384,7 +392,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
if (debug) if (debug)
{ {
Pout<< " separation_:" << separation_ << nl Pout<< " separation_:" << separation_.size() << nl
<< " forwardT size:" << forwardT_.size() << endl; << " forwardT size:" << forwardT_.size() << endl;
} }
} }

173
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
View File

@ -127,6 +127,22 @@ void Foam::cyclicPolyPatch::calcTransforms()
Pout<< "cyclicPolyPatch::calcTransforms : Writing half1" Pout<< "cyclicPolyPatch::calcTransforms : Writing half1"
<< " faces to OBJ file " << nm1 << endl; << " faces to OBJ file " << nm1 << endl;
writeOBJ(nm1, half1, half1.points()); writeOBJ(nm1, half1, half1.points());
OFstream str(casePath/name()+"_half0_to_half1.obj");
label vertI = 0;
Pout<< "cyclicPolyPatch::calcTransforms :"
<< " Writing coupled face centres as lines to " << str.name()
<< endl;
forAll(half0Ctrs, i)
{
const point& p0 = half0Ctrs[i];
str << "v " << p0.x() << ' ' << p0.y() << ' ' << p0.z() << nl;
vertI++;
const point& p1 = half1Ctrs[i];
str << "v " << p1.x() << ' ' << p1.y() << ' ' << p1.z() << nl;
vertI++;
str << "l " << vertI-1 << ' ' << vertI << nl;
}
} }
vectorField half0Normals(half0.size()); vectorField half0Normals(half0.size());
@ -397,8 +413,6 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
anchors0 = getAnchorPoints(half0Faces, pp.points()); anchors0 = getAnchorPoints(half0Faces, pp.points());
half1Ctrs = calcFaceCentres(half1Faces, pp.points()); half1Ctrs = calcFaceCentres(half1Faces, pp.points());
vector n0 = vector::zero;
vector n1 = vector::zero;
switch (transform_) switch (transform_)
{ {
case ROTATIONAL: case ROTATIONAL:
@ -406,12 +420,46 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
label face0 = getConsistentRotationFace(half0Ctrs); label face0 = getConsistentRotationFace(half0Ctrs);
label face1 = getConsistentRotationFace(half1Ctrs); label face1 = getConsistentRotationFace(half1Ctrs);
n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_); vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_); vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
n0 /= mag(n0) + VSMALL; n0 /= mag(n0) + VSMALL;
n1 /= mag(n1) + VSMALL; n1 /= mag(n1) + VSMALL;
if (debug)
{
Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " Specified rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl;
}
// Rotation (around origin)
const tensor reverseT(rotationTensor(n0, -n1));
// Rotation
forAll(half0Ctrs, faceI)
{
half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
}
break; break;
} }
//- Problem: usually specified translation is not accurate enough
//- to get proper match so keep automatic determination over here.
//case TRANSLATIONAL:
//{
// // Transform 0 points.
//
// if (debug)
// {
// Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
// << "Specified translation : " << separationVector_ << endl;
// }
//
// half0Ctrs += separationVector_;
// anchors0 += separationVector_;
// break;
//}
default: default:
{ {
// Assumes that cyclic is planar. This is also the initial // Assumes that cyclic is planar. This is also the initial
@ -420,55 +468,67 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
// Determine the face with max area on both halves. These // Determine the face with max area on both halves. These
// two faces are used to determine the transformation tensors // two faces are used to determine the transformation tensors
label max0I = findMaxArea(pp.points(), half0Faces); label max0I = findMaxArea(pp.points(), half0Faces);
n0 = half0Faces[max0I].normal(pp.points()); vector n0 = half0Faces[max0I].normal(pp.points());
n0 /= mag(n0) + VSMALL; n0 /= mag(n0) + VSMALL;
label max1I = findMaxArea(pp.points(), half1Faces); label max1I = findMaxArea(pp.points(), half1Faces);
n1 = half1Faces[max1I].normal(pp.points()); vector n1 = half1Faces[max1I].normal(pp.points());
n1 /= mag(n1) + VSMALL; n1 /= mag(n1) + VSMALL;
if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
{
if (debug)
{
Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " Detected rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl;
}
// Rotation (around origin)
const tensor reverseT(rotationTensor(n0, -n1));
// Rotation
forAll(half0Ctrs, faceI)
{
half0Ctrs[faceI] = Foam::transform
(
reverseT,
half0Ctrs[faceI]
);
anchors0[faceI] = Foam::transform
(
reverseT,
anchors0[faceI]
);
}
}
else
{
// Parallel translation. Get average of all used points.
primitiveFacePatch half0(half0Faces, pp.points());
const pointField& half0Pts = half0.localPoints();
const point ctr0(sum(half0Pts)/half0Pts.size());
primitiveFacePatch half1(half1Faces, pp.points());
const pointField& half1Pts = half1.localPoints();
const point ctr1(sum(half1Pts)/half1Pts.size());
if (debug)
{
Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
<< " Detected translation :"
<< " n0:" << n0 << " n1:" << n1
<< " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
}
half0Ctrs += ctr1 - ctr0;
anchors0 += ctr1 - ctr0;
}
break;
} }
} }
if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
{
if (debug)
{
Pout<< "cyclicPolyPatch::getCentresAndAnchors : Rotation :"
<< " n0:" << n0 << " n1:" << n1 << endl;
}
// Rotation (around origin)
const tensor reverseT(rotationTensor(n0, -n1));
// Rotation
forAll(half0Ctrs, faceI)
{
half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
}
}
else
{
// Parallel translation. Get average of all used points.
primitiveFacePatch half0(half0Faces, pp.points());
const pointField& half0Pts = half0.localPoints();
const point ctr0(sum(half0Pts)/half0Pts.size());
primitiveFacePatch half1(half1Faces, pp.points());
const pointField& half1Pts = half1.localPoints();
const point ctr1(sum(half1Pts)/half1Pts.size());
if (debug)
{
Pout<< "cyclicPolyPatch::getCentresAndAnchors : Translation :"
<< " n0:" << n0 << " n1:" << n1
<< " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
}
half0Ctrs += ctr1 - ctr0;
anchors0 += ctr1 - ctr0;
}
// Calculate typical distance per face // Calculate typical distance per face
tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs); tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
@ -615,7 +675,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(0.9), featureCos_(0.9),
transform_(UNKNOWN), transform_(UNKNOWN),
rotationAxis_(vector::zero), rotationAxis_(vector::zero),
rotationCentre_(point::zero) rotationCentre_(point::zero),
separationVector_(vector::zero)
{ {
calcTransforms(); calcTransforms();
} }
@ -635,7 +696,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(0.9), featureCos_(0.9),
transform_(UNKNOWN), transform_(UNKNOWN),
rotationAxis_(vector::zero), rotationAxis_(vector::zero),
rotationCentre_(point::zero) rotationCentre_(point::zero),
separationVector_(vector::zero)
{ {
dict.readIfPresent("featureCos", featureCos_); dict.readIfPresent("featureCos", featureCos_);
@ -650,9 +712,14 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
dict.lookup("rotationCentre") >> rotationCentre_; dict.lookup("rotationCentre") >> rotationCentre_;
break; break;
} }
case TRANSLATIONAL:
{
dict.lookup("separationVector") >> separationVector_;
break;
}
default: default:
{ {
// no additioanl info required // no additional info required
} }
} }
} }
@ -673,7 +740,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(pp.featureCos_), featureCos_(pp.featureCos_),
transform_(pp.transform_), transform_(pp.transform_),
rotationAxis_(pp.rotationAxis_), rotationAxis_(pp.rotationAxis_),
rotationCentre_(pp.rotationCentre_) rotationCentre_(pp.rotationCentre_),
separationVector_(pp.separationVector_)
{ {
calcTransforms(); calcTransforms();
} }
@ -694,7 +762,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(pp.featureCos_), featureCos_(pp.featureCos_),
transform_(pp.transform_), transform_(pp.transform_),
rotationAxis_(pp.rotationAxis_), rotationAxis_(pp.rotationAxis_),
rotationCentre_(pp.rotationCentre_) rotationCentre_(pp.rotationCentre_),
separationVector_(pp.separationVector_)
{ {
calcTransforms(); calcTransforms();
} }
@ -1322,6 +1391,8 @@ void Foam::cyclicPolyPatch::write(Ostream& os) const
{ {
os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL] os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL]
<< token::END_STATEMENT << nl; << token::END_STATEMENT << nl;
os.writeKeyword("separationVector") << separationVector_
<< token::END_STATEMENT << nl;
break; break;
} }
default: default:

139
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
View File

@ -101,11 +101,18 @@ private:
//- Type of transformation - rotational or translational //- Type of transformation - rotational or translational
transformType transform_; transformType transform_;
//- Axis of rotation for rotational cyclics // For rotation
vector rotationAxis_;
//- point on axis of rotation for rotational cyclics //- Axis of rotation for rotational cyclics
point rotationCentre_; vector rotationAxis_;
//- point on axis of rotation for rotational cyclics
point rotationCentre_;
// For translation
//- Translation vector
vector separationVector_;
// Private member functions // Private member functions
@ -267,66 +274,94 @@ public:
const edgeList& coupledEdges() const; const edgeList& coupledEdges() const;
vector separation(const label facei) const
{
if (facei < size()/2)
{
return coupledPolyPatch::separation()[0];
}
else
{
return -coupledPolyPatch::separation()[0];
}
}
const tensor& transformT(const label facei) const // Transformation
{
if (facei < size()/2)
{
return reverseT()[0];
}
else
{
return forwardT()[0];
}
}
template<class T> vector separation(const label facei) const
T transform(const T& t, const label facei) const
{
if (parallel())
{ {
return t; if (facei < size()/2)
{
return coupledPolyPatch::separation()[0];
}
else
{
return -coupledPolyPatch::separation()[0];
}
} }
else
{
return Foam::transform(transformT(facei), t);
}
}
label transformLocalFace(const label facei) const const tensor& transformT(const label facei) const
{
if (facei < size()/2)
{ {
return facei + size()/2; if (facei < size()/2)
{
return reverseT()[0];
}
else
{
return forwardT()[0];
}
} }
else
{
return facei - size()/2;
}
}
label transformGlobalFace(const label facei) const template<class T>
{ T transform(const T& t, const label facei) const
if (facei - start() < size()/2)
{ {
return facei + size()/2; if (parallel())
{
return t;
}
else
{
return Foam::transform(transformT(facei), t);
}
} }
else
label transformLocalFace(const label facei) const
{ {
return facei - size()/2; if (facei < size()/2)
{
return facei + size()/2;
}
else
{
return facei - size()/2;
}
} }
}
label transformGlobalFace(const label facei) const
{
if (facei - start() < size()/2)
{
return facei + size()/2;
}
else
{
return facei - size()/2;
}
}
//- Type of transform
transformType transform() const
{
return transform_;
}
//- Axis of rotation for rotational cyclics
const vector& rotationAxis() const
{
return rotationAxis_;
}
//- point on axis of rotation for rotational cyclics
const point& rotationCentre() const
{
return rotationCentre_;
}
//- Translation vector for translational cyclics
const vector& separationVector() const
{
return separationVector_;
}
//- Initialize ordering for primitivePatch. Does not //- Initialize ordering for primitivePatch. Does not

32
src/postProcessing/fieldAverage/fieldAverageItem/fieldAverageItem.H
View File

@ -102,15 +102,6 @@ private:
baseType base_; baseType base_;
// Private Member Functions
//- Disallow default bitwise copy construct
// fieldAverageItem(const fieldAverageItem&);
//- Disallow default bitwise assignment
// void operator=(const fieldAverageItem&);
public: public:
// Constructors // Constructors
@ -175,6 +166,29 @@ public:
void operator=(const fieldAverageItem&); void operator=(const fieldAverageItem&);
// Friend Operators
friend bool operator==
(
const fieldAverageItem& a,
const fieldAverageItem& b
)
{
return
a.fieldName_ == b.fieldName_
&& a.mean_ == b.mean_
&& a.prime2Mean_ == b.prime2Mean_
&& a.base_ == b.base_;
}
friend bool operator!=
(
const fieldAverageItem& a,
const fieldAverageItem& b
)
{
return !(a == b);
}
// IOstream Operators // IOstream Operators

14
tutorials/channelOodles/channel395/constant/postChannelDict
View File

@ -15,14 +15,14 @@ FoamFile
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Nx 40; // Seed patches to start layering from
Ny patches (bottomWall);
(
25
25
);
Nz 30;
// Direction in which the layers are
component y;
// Is the mesh symmetric? If so average(symmetric fields) or
// subtract(asymmetric) contributions from both halves
symmetric true; symmetric true;
// ************************************************************************* // // ************************************************************************* //