diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C b/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C
index 4db4faac505f40e5b0b8a40d152a48843bf9d119..5a17820680fc602afaad19b9c931fb21d1c8c171 100644
--- a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C
+++ b/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C
@@ -157,9 +157,8 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
if(thermalCreep_)
{
- const GeometricField<scalar, fvPatchField, volMesh>& vsfT =
- this->db().objectRegistry::
- lookupObject<GeometricField<scalar, fvPatchField, volMesh> >("T");
+ const volScalarField& vsfT =
+ this->db().objectRegistry::lookupObject<volScalarField>("T");
label patchi = this->patch().index();
const fvPatchScalarField& pT = vsfT.boundaryField()[patchi];
Field<vector> gradpT = fvc::grad(vsfT)().boundaryField()[patchi];
diff --git a/applications/solvers/incompressible/oodles/oodles.C b/applications/solvers/incompressible/oodles/oodles.C
index eb5fa65dcf6bd67668c26dd3d2026e658f711f7e..8c51065db1d06490ed8fba17107553c52f0b2be0 100644
--- a/applications/solvers/incompressible/oodles/oodles.C
+++ b/applications/solvers/incompressible/oodles/oodles.C
@@ -66,6 +66,9 @@ int main(int argc, char *argv[])
+ sgsModel->divDevBeff(U)
);
+ // Optionally ensure diagonal-dominance of the momentum matrix
+ UEqn.relax();
+
if (momentumPredictor)
{
solve(UEqn == -fvc::grad(p));
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/Make/files b/applications/utilities/mesh/manipulation/checkMesh.save/Make/files
deleted file mode 100644
index f0b7c166946761f10023900d29f2bb92e60b67b6..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/Make/files
+++ /dev/null
@@ -1,6 +0,0 @@
-printMeshStats.C
-checkTopology.C
-checkGeometry.C
-checkMesh.C
-
-EXE = $(FOAM_APPBIN)/checkMesh
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/Make/options b/applications/utilities/mesh/manipulation/checkMesh.save/Make/options
deleted file mode 100644
index 54c035b8f55d183c1ad02bc372398feceaf31718..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/Make/options
+++ /dev/null
@@ -1,5 +0,0 @@
-EXE_INC = \
- -I$(LIB_SRC)/meshTools/lnInclude
-
-EXE_LIBS = \
- -lmeshTools
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.C
deleted file mode 100644
index 9eb4ef63ad7d362ebdb66931c288ea67d81d76f2..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.C
+++ /dev/null
@@ -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;
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.H b/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.H
deleted file mode 100644
index 791b007820d76fa8b4dff2626c8fcdccfe02a8ee..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.H
+++ /dev/null
@@ -1,13 +0,0 @@
-#include "label.H"
-
-namespace Foam
-{
- class polyMesh;
-
- label checkGeometry
- (
- const polyMesh& mesh,
- bool checkPointNearness,
- bool checkCellDeterminant
- );
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkMesh.C b/applications/utilities/mesh/manipulation/checkMesh.save/checkMesh.C
deleted file mode 100644
index ec598b4c215b45114ed9e9dc7268582c6d688dc5..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkMesh.C
+++ /dev/null
@@ -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);
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.C
deleted file mode 100644
index 885c5f36aa9fa94021fc4f5667d60e29443dd63c..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.C
+++ /dev/null
@@ -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;
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.H b/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.H
deleted file mode 100644
index 6db7e0d3c502623329afc29badb2692a15104f22..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.H
+++ /dev/null
@@ -1,8 +0,0 @@
-#include "label.H"
-
-namespace Foam
-{
- class polyMesh;
-
- label checkTopology(const polyMesh& mesh, bool fullTopology);
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.C b/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.C
deleted file mode 100644
index 30da859db629c8b227571343cf919395b9b6eaf5..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.C
+++ /dev/null
@@ -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;
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.H b/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.H
deleted file mode 100644
index 2da12011883452de3d88f4b1df938906b0cf081b..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.H
+++ /dev/null
@@ -1,6 +0,0 @@
-namespace Foam
-{
- class polyMesh;
-
- void printMeshStats(const polyMesh& mesh);
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 6b1a10d90e6f7e1b771a7672e0f5bafdd28bc800..8eb1f3f05fbb1a92ef5c889418e0741796e7dd72 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -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();
}
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index 85dc0bae62720f8bf7f8fd1219624cc176c29289..ac924136e1cf104e13e9a1fdc861e9ad6ba9c911 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -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();
diff --git a/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C b/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
index 9c608b689d7958ba94344a5208124cdafcc2d3c2..9a2c9396afba0a645483f86eec3d65239372cb43 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
@@ -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
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
index a44966198217a88a5169c38d1fbd0b3e149fa1c5..01df2c2f284329f542596fd141e33c2318660b94 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
@@ -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()));
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatchDict b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
index 1153d0204d7db267d52ed8830556e52fe1010032..5f3597f21af7ebb740e5dac56e5076c26fee4953 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatchDict
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
@@ -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;
+ //}
);
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options b/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
index d38cd8b1801b429e5f42c7b0ba2ab58ba2cd1d11..b90b6fdd4e1c3b114664c9123a92ff043fca7244 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
@@ -1,7 +1,9 @@
EXE_INC = \
+ -I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
+ -lmeshTools \
-lfiniteVolume \
-lsampling
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
index 74315de31a44e057b2c64017cf9d2c343081048d..cdba8c0201dd0ce865a5750ea7b584990f584d71 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
@@ -25,156 +25,269 @@ License
\*---------------------------------------------------------------------------*/
#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)
-:
- 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])
+template<>
+const char* Foam::NamedEnum<Foam::vector::components, 3>::names[] =
{
- // initialise the layers
- cumNy_[0] = ny_[0];
-
- for (label j=1; j<ny_.size(); j++)
- {
- nLayers_ += ny_[j];
- cumNy_[j] = ny_[j]+cumNy_[j-1];
- }
-}
-
+ "x",
+ "y",
+ "z"
+};
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+const Foam::NamedEnum<Foam::vector::components, 3>
+ Foam::channelIndex::vectorComponentsNames_;
-channelIndex::~channelIndex()
-{}
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-scalarField channelIndex::collapse
+// Determines face blocking
+void Foam::channelIndex::walkOppositeFaces
(
- const volScalarField& vsf,
- const bool asymmetric
-) const
+ const polyMesh& mesh,
+ const labelList& startFaces,
+ boolList& blockedFace
+)
{
- scalarField cs(nLayers(), 0.0);
+ const cellList& cells = mesh.cells();
+ const faceList& faces = mesh.faces();
+ label nBnd = mesh.nFaces() - mesh.nInternalFaces();
+
+ DynamicList<label> frontFaces(startFaces);
+ forAll(frontFaces, i)
+ {
+ label faceI = frontFaces[i];
+ blockedFace[faceI] = true;
+ }
- forAll(cs, j)
+ while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
{
- // sweep over all cells in this layer
- for (label i=0; i<nx(); i++)
+ // Transfer across.
+ boolList isFrontBndFace(nBnd, false);
+ forAll(frontFaces, i)
{
- for (label k=0; k<nz(); k++)
+ label faceI = frontFaces[i];
+
+ if (!mesh.isInternalFace(faceI))
{
- cs[j] += vsf[operator()(i,j,k)];
+ isFrontBndFace[faceI-mesh.nInternalFaces()] = true;
}
}
+ syncTools::swapBoundaryFaceList(mesh, isFrontBndFace, false);
- // and divide by the number of cells in the layer
- cs[j] /= scalar(nx()*nz());
- }
+ // Add
+ forAll(isFrontBndFace, i)
+ {
+ label faceI = mesh.nInternalFaces()+i;
+ if (isFrontBndFace[i] && !blockedFace[faceI])
+ {
+ blockedFace[faceI] = true;
+ frontFaces.append(faceI);
+ }
+ }
- if (symmetric_)
- {
- label nlb2 = nLayers()/2;
+ // Transfer across cells
+ DynamicList<label> newFrontFaces(frontFaces.size());
- if (asymmetric)
+ forAll(frontFaces, i)
{
- for (label j=0; j<nlb2; j++)
+ label faceI = frontFaces[i];
+
{
- cs[j] = 0.5*(cs[j] - cs[nLayers() - j - 1]);
+ 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);
+ }
+ }
}
}
- else
+
+ frontFaces.transfer(newFrontFaces);
+ }
+}
+
+
+// Calculate regions.
+void Foam::channelIndex::calcLayeredRegions
+(
+ const polyMesh& mesh,
+ const labelList& startFaces
+)
+{
+ boolList blockedFace(mesh.nFaces(), false);
+ walkOppositeFaces
+ (
+ mesh,
+ startFaces,
+ blockedFace
+ );
+
+
+ if (false)
+ {
+ OFstream str(mesh.time().path()/"blockedFaces.obj");
+ label vertI = 0;
+ forAll(blockedFace, faceI)
{
- for (label j=0; j<nlb2; j++)
+ if (blockedFace[faceI])
{
- cs[j] = 0.5*(cs[j] + cs[nLayers() - j - 1]);
+ 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();
}
}
-
- cs.setSize(nlb2);
}
- return cs;
+
+ // 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);
+ }
}
-scalarField channelIndex::y
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::channelIndex::channelIndex
(
- const volVectorField& cellCentres
-) const
+ const polyMesh& mesh,
+ const dictionary& dict
+)
+:
+ symmetric_(readBool(dict.lookup("symmetric"))),
+ dir_(vectorComponentsNames_.read(dict.lookup("component")))
{
- if (symmetric_)
+ const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+ const wordList patchNames(dict.lookup("patches"));
+
+ label nFaces = 0;
+
+ forAll(patchNames, i)
{
- scalarField Y(nLayers()/2);
+ label patchI = patches.findPatchID(patchNames[i]);
- for (label j=0; j<nLayers()/2; j++)
+ if (patchI == -1)
{
- Y[j] = cellCentres[operator()(0, j, 0)].y();
+ FatalErrorIn("channelIndex::channelIndex(const polyMesh&)")
+ << "Illegal patch " << patchNames[i]
+ << ". Valid patches are " << patches.name()
+ << exit(FatalError);
}
- return Y;
+ nFaces += patches[patchI].size();
}
- else
+
+ labelList startFaces(nFaces);
+ nFaces = 0;
+
+ forAll(patchNames, i)
{
- scalarField Y(nLayers());
+ const polyPatch& pp = patches[patches.findPatchID(patchNames[i])];
- for (label j=0; j<nLayers(); j++)
+ forAll(pp, j)
{
- Y[j] = cellCentres[operator()(0, j, 0)].y();
+ startFaces[nFaces++] = pp.start()+j;
}
-
- return Y;
}
+
+ // Calculate regions.
+ calcLayeredRegions(mesh, startFaces);
}
-// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
-
-label channelIndex::operator()
+Foam::channelIndex::channelIndex
(
- const label Jx,
- const label Jy,
- const label Jz
-) const
+ const polyMesh& mesh,
+ const labelList& startFaces,
+ const bool symmetric,
+ const direction dir
+)
+:
+ symmetric_(symmetric),
+ dir_(dir)
{
- label index(0);
+ // Calculate regions.
+ calcLayeredRegions(mesh, startFaces);
+}
- // 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++;
- }
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
- index += Jx + nx_*tmpJy + nx_*ny_[j]*Jz;
- return index;
-}
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
index 500d13f4eb00469ee0e1599db2535efa025116e5..796e20160c4bb6a53c5352d403059035d0edc9c0 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
@@ -26,19 +26,25 @@ Class
Foam::channelIndex
Description
- does indexing for a standard (Christer-style) channel.
-Assumes that the blocks are aranged in the y direction.
+ Does averaging of fields over layers of cells. Assumes layered mesh.
SourceFiles
channelIndex.C
- channelIndexIO.C
\*---------------------------------------------------------------------------*/
#ifndef 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
{
+
// Private data
- IOdictionary indexingDict_;
+ static const NamedEnum<vector::components, 3> vectorComponentsNames_;
- const label nx_;
- const labelList ny_;
- const label nz_;
+ //- Is mesh symmetric
const bool symmetric_;
- labelList cumNy_;
- label nLayers_;
+ //- direction to sort
+ 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
+ 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
channelIndex(const channelIndex&);
void operator=(const channelIndex&);
@@ -71,55 +103,53 @@ public:
// Constructors
- channelIndex(const fvMesh& m);
-
+ //- Construct from dictionary
+ channelIndex(const polyMesh&, const dictionary&);
- // Destructor
-
- ~channelIndex();
+ //- Construct from supplied starting faces
+ channelIndex
+ (
+ const polyMesh& mesh,
+ const labelList& startFaces,
+ const bool symmetric,
+ const direction dir
+ );
// Member Functions
// Access
- //- return number of layers
- label nLayers() 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_;
- }
+ //- Sum field per region
+ template<class T>
+ Field<T> regionSum(const Field<T>& cellField) const;
//- 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;
//- return the field of Y locations from the cell centres
- scalarField y
- (
- const volVectorField& cellCentres
- ) const;
+ const scalarField& y() const
+ {
+ return y_;
+ }
+};
- // Member Operators
- label operator()(const label, const label, const label) const;
-};
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+# include "channelIndexTemplates.C"
+#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..9fc8ff03609b0681b67d664a54b70293b493091d
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
@@ -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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H b/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
index ec047e58627ce52b4edd1b33cf2e21ffec7d596f..b3bf5594110835698535275ed2f58b647a7e238d 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
@@ -1,16 +1,16 @@
scalarField UMeanXvalues = channelIndexing.collapse
(
- UMean.component(vector::X)
+ UMean.component(vector::X)()
);
scalarField UMeanYvalues = channelIndexing.collapse
(
- UMean.component(vector::Y)
+ UMean.component(vector::Y)()
);
scalarField UMeanZvalues = channelIndexing.collapse
(
- UMean.component(vector::Z)
+ UMean.component(vector::Z)()
);
scalarField RxxValues = channelIndexing.collapse(Rxx);
@@ -38,7 +38,7 @@
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, urmsValues, "u", UMean.path(), gFormat);
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
index 2de4c534a8f1751e09c29b8abdac0d61b0c4ee85..08051cb1da9bf18a47c54d1dfe845949d212b4d9 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
@@ -67,7 +67,20 @@ int main(int argc, char *argv[])
const word& gFormat = runTime.graphFormat();
// 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 (label i=startTime; i<endTime; i++)
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict
index e69c0386ef95781462a1b01424c2d8eee53f1e2c..d63cd656675c6ebe4aecf1608dca9b561a5a3fe8 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict
@@ -15,14 +15,14 @@ FoamFile
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-Nx 40;
-Ny
-(
- 25
- 25
-);
-Nz 30;
+// Seed patches to start layering from
+patches (bottomWall);
+// 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;
// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.C b/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.C
new file mode 100644
index 0000000000000000000000000000000000000000..8662ec578b8ece2f4bdcff0a80c2abdf8b9dc0c4
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.C
@@ -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_;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.H b/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.H
new file mode 100644
index 0000000000000000000000000000000000000000..15a26a3aaa7142714929938c9f0dad6f2b3d84c7
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.H
@@ -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
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/sumDataI.H b/applications/utilities/postProcessing/miscellaneous/postChannel/sumDataI.H
new file mode 100644
index 0000000000000000000000000000000000000000..107cba063dccdc7006cb43da53913ce11a29e6dc
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/sumDataI.H
@@ -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);
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index a74458fe9fce68936ec738ba775192f3f883fc5d..e7d7567a138da742b2eb21630b1b0866aad718ad 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -552,6 +552,8 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
nextPatchStart += bm[patchI].size();
}
+ reduce(boundaryError, orOp<bool>());
+
if (boundaryError)
{
if (debug || report)
@@ -565,7 +567,7 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
{
if (debug || report)
{
- Pout << " Boundary definition OK." << endl;
+ Info << " Boundary definition OK." << endl;
}
return false;
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
index f4e073d4e14b1d46e12938234ac6616b9161e7db..73aff8a1fd2c711402e5ddd3a201583166df211b 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
@@ -263,6 +263,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl
<< " (half)size:" << Cf.size() << nl
<< " absTol:" << absTol << nl
+ //<< " smallDist:" << smallDist << nl
<< " sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl;
}
@@ -316,6 +317,13 @@ void Foam::coupledPolyPatch::calcTransformTensors
{
forwardT_.setSize(1);
reverseT_.setSize(1);
+
+ if (debug)
+ {
+ Pout<< " rotation " << sum(mag(forwardT_ - forwardT_[0]))
+ << " more than local tolerance " << error
+ << ". Assuming uniform rotation." << endl;
+ }
}
}
else
@@ -384,7 +392,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
if (debug)
{
- Pout<< " separation_:" << separation_ << nl
+ Pout<< " separation_:" << separation_.size() << nl
<< " forwardT size:" << forwardT_.size() << endl;
}
}
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
index 31faea5669f6cb899fb41f8aaaa81180771ff215..171fb099f1a08ce774a19e2421c94ee83249cdb4 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
@@ -127,6 +127,22 @@ void Foam::cyclicPolyPatch::calcTransforms()
Pout<< "cyclicPolyPatch::calcTransforms : Writing half1"
<< " faces to OBJ file " << nm1 << endl;
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());
@@ -397,8 +413,6 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
anchors0 = getAnchorPoints(half0Faces, pp.points());
half1Ctrs = calcFaceCentres(half1Faces, pp.points());
- vector n0 = vector::zero;
- vector n1 = vector::zero;
switch (transform_)
{
case ROTATIONAL:
@@ -406,12 +420,46 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
label face0 = getConsistentRotationFace(half0Ctrs);
label face1 = getConsistentRotationFace(half1Ctrs);
- n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
- n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
+ vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
+ vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
n0 /= mag(n0) + 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;
}
+ //- 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:
{
// Assumes that cyclic is planar. This is also the initial
@@ -420,56 +468,68 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
// Determine the face with max area on both halves. These
// two faces are used to determine the transformation tensors
label max0I = findMaxArea(pp.points(), half0Faces);
- n0 = half0Faces[max0I].normal(pp.points());
+ vector n0 = half0Faces[max0I].normal(pp.points());
n0 /= mag(n0) + VSMALL;
label max1I = findMaxArea(pp.points(), half1Faces);
- n1 = half1Faces[max1I].normal(pp.points());
+ vector n1 = half1Faces[max1I].normal(pp.points());
n1 /= mag(n1) + VSMALL;
- }
- }
- if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
- {
- if (debug)
- {
- Pout<< "cyclicPolyPatch::getCentresAndAnchors : Rotation :"
- << " n0:" << n0 << " n1:" << n1 << endl;
- }
+ 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 (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.
+ // 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 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());
+ 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;
- }
+ if (debug)
+ {
+ Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+ << " Detected translation :"
+ << " n0:" << n0 << " n1:" << n1
+ << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
+ }
- half0Ctrs += ctr1 - ctr0;
- anchors0 += ctr1 - ctr0;
+ half0Ctrs += ctr1 - ctr0;
+ anchors0 += ctr1 - ctr0;
+ }
+ break;
+ }
}
+
// Calculate typical distance per face
tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
}
@@ -615,7 +675,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(0.9),
transform_(UNKNOWN),
rotationAxis_(vector::zero),
- rotationCentre_(point::zero)
+ rotationCentre_(point::zero),
+ separationVector_(vector::zero)
{
calcTransforms();
}
@@ -635,7 +696,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(0.9),
transform_(UNKNOWN),
rotationAxis_(vector::zero),
- rotationCentre_(point::zero)
+ rotationCentre_(point::zero),
+ separationVector_(vector::zero)
{
dict.readIfPresent("featureCos", featureCos_);
@@ -650,9 +712,14 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
dict.lookup("rotationCentre") >> rotationCentre_;
break;
}
+ case TRANSLATIONAL:
+ {
+ dict.lookup("separationVector") >> separationVector_;
+ break;
+ }
default:
{
- // no additioanl info required
+ // no additional info required
}
}
}
@@ -673,7 +740,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(pp.featureCos_),
transform_(pp.transform_),
rotationAxis_(pp.rotationAxis_),
- rotationCentre_(pp.rotationCentre_)
+ rotationCentre_(pp.rotationCentre_),
+ separationVector_(pp.separationVector_)
{
calcTransforms();
}
@@ -694,7 +762,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
featureCos_(pp.featureCos_),
transform_(pp.transform_),
rotationAxis_(pp.rotationAxis_),
- rotationCentre_(pp.rotationCentre_)
+ rotationCentre_(pp.rotationCentre_),
+ separationVector_(pp.separationVector_)
{
calcTransforms();
}
@@ -1322,6 +1391,8 @@ void Foam::cyclicPolyPatch::write(Ostream& os) const
{
os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL]
<< token::END_STATEMENT << nl;
+ os.writeKeyword("separationVector") << separationVector_
+ << token::END_STATEMENT << nl;
break;
}
default:
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
index f9c4e66dbb2c69cb6e2cdae369ae1c80f35b7bb4..7db828985fb8e37dcb78820e0d1dec0989484893 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
@@ -101,11 +101,18 @@ private:
//- Type of transformation - rotational or translational
transformType transform_;
- //- Axis of rotation for rotational cyclics
- vector rotationAxis_;
+ // For rotation
- //- point on axis of rotation for rotational cyclics
- point rotationCentre_;
+ //- Axis of rotation for rotational cyclics
+ vector rotationAxis_;
+
+ //- point on axis of rotation for rotational cyclics
+ point rotationCentre_;
+
+ // For translation
+
+ //- Translation vector
+ vector separationVector_;
// Private member functions
@@ -267,66 +274,94 @@ public:
const edgeList& coupledEdges() const;
- vector separation(const label facei) const
- {
- if (facei < size()/2)
- {
- return coupledPolyPatch::separation()[0];
- }
- else
+
+ // Transformation
+
+ vector separation(const label facei) const
{
- return -coupledPolyPatch::separation()[0];
+ if (facei < size()/2)
+ {
+ return coupledPolyPatch::separation()[0];
+ }
+ else
+ {
+ return -coupledPolyPatch::separation()[0];
+ }
}
- }
- const tensor& transformT(const label facei) const
- {
- if (facei < size()/2)
+ const tensor& transformT(const label facei) const
{
- return reverseT()[0];
+ if (facei < size()/2)
+ {
+ return reverseT()[0];
+ }
+ else
+ {
+ return forwardT()[0];
+ }
}
- else
+
+ template<class T>
+ T transform(const T& t, const label facei) const
{
- return forwardT()[0];
+ if (parallel())
+ {
+ return t;
+ }
+ else
+ {
+ return Foam::transform(transformT(facei), t);
+ }
}
- }
- template<class T>
- T transform(const T& t, const label facei) const
- {
- if (parallel())
+ label transformLocalFace(const label facei) const
{
- return t;
+ if (facei < size()/2)
+ {
+ return facei + size()/2;
+ }
+ else
+ {
+ return facei - size()/2;
+ }
}
- else
+
+ label transformGlobalFace(const label facei) const
{
- return Foam::transform(transformT(facei), t);
+ if (facei - start() < size()/2)
+ {
+ return facei + size()/2;
+ }
+ else
+ {
+ return facei - size()/2;
+ }
}
- }
- label transformLocalFace(const label facei) const
- {
- if (facei < size()/2)
+ //- Type of transform
+ transformType transform() const
{
- return facei + size()/2;
+ return transform_;
}
- else
+
+ //- Axis of rotation for rotational cyclics
+ const vector& rotationAxis() const
{
- return facei - size()/2;
+ return rotationAxis_;
}
- }
- label transformGlobalFace(const label facei) const
- {
- if (facei - start() < size()/2)
+ //- point on axis of rotation for rotational cyclics
+ const point& rotationCentre() const
{
- return facei + size()/2;
+ return rotationCentre_;
}
- else
+
+ //- Translation vector for translational cyclics
+ const vector& separationVector() const
{
- return facei - size()/2;
+ return separationVector_;
}
- }
+
//- Initialize ordering for primitivePatch. Does not
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
index 0caebcd4218570ca9ecf90cacceb33022f597317..79f9107e97f61b4b9f77a364256dca9c021f3620 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
@@ -113,7 +113,9 @@ Foam::quadraticFitData::quadraticFitData
{
interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
}
+
interpPolySize.write();
+
if (debug)
{
Info<< "quadraticFitData::quadraticFitData() :"
@@ -228,9 +230,9 @@ Foam::label Foam::quadraticFitData::calcFit
pz /= scale;
label is = 0;
- B[ip][is++] = wts[ip];
- B[ip][is++] = wts[ip]*px;
+ B[ip][is++] = wts[0]*wts[ip];
+ B[ip][is++] = wts[0]*wts[ip]*px;
B[ip][is++] = wts[ip]*sqr(px);
if (dim_ >= 2)
@@ -254,46 +256,96 @@ Foam::label Foam::quadraticFitData::calcFit
scalarList singVals(minSize_);
label nSVDzeros = 0;
+ const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
+ mesh().surfaceInterpolation::weights();
+
bool goodFit = false;
- for(scalar tol = SMALL; tol < 0.1 && !goodFit; tol *= 10)
+ for(int iIt = 0; iIt < 10 && !goodFit; iIt++)
{
- SVD svd(B, tol);
+ SVD svd(B, SMALL);
+
+ scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0];
+ scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1];
+
+ //goodFit = (fit0 > 0 && fit1 > 0);
- scalar fit0 = wts[0]*svd.VSinvUt()[0][0];
- scalar fit1 = wts[1]*svd.VSinvUt()[0][1];
+ goodFit =
+ (mag(fit0 - w[faci])/w[faci] < 0.5)
+ && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < 0.5);
- goodFit = sign(fit0) == sign(fit1);
+ //scalar w0Err = fit0/w[faci];
+ //scalar w1Err = fit1/(1 - w[faci]);
+
+ //goodFit =
+ // (w0Err > 0.5 && w0Err < 1.5)
+ // && (w1Err > 0.5 && w1Err < 1.5);
if (goodFit)
{
fit_[faci][0] = fit0;
fit_[faci][1] = fit1;
- for(label i = 2; i < stencilSize; i++)
+
+ for(label i=2; i<stencilSize; i++)
{
- fit_[faci][i] = wts[i]*svd.VSinvUt()[0][i];
+ fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
}
singVals = svd.S();
nSVDzeros = svd.nZeros();
}
- }
- if (!goodFit)
- {
- FatalErrorIn
- (
- "quadraticFitData::calcFit(const pointField&, const label)"
- ) << "For face " << faci << endl
- << "Fit not good even once tol >= 0.1"
- << exit(FatalError);
- }
+ else // (not good fit so increase weight in the centre and for linear)
+ {
+ wts[0] *= 10;
+ wts[1] *= 10;
- const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
- mesh().surfaceInterpolation::weights();
+ for(label i = 0; i < B.n(); i++)
+ {
+ B[i][0] *= 10;
+ B[i][1] *= 10;
+ }
+
+ for(label j = 0; j < B.m(); j++)
+ {
+ B[0][j] *= 10;
+ B[1][j] *= 10;
+ }
+ }
+ }
- // remove the uncorrected linear coefficients
+ //static const scalar alpha = 1.5;
+ //static const scalar beta = alpha/0.5;
- fit_[faci][0] -= w[faci];
- fit_[faci][1] -= 1 - w[faci];
+ if (goodFit)
+ {
+ // scalar limiter =
+ // max
+ // (
+ // min
+ // (
+ // min(alpha - beta*mag(fit_[faci][0] - w[faci])/w[faci], 1),
+ // min(alpha - beta*mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]), 1)
+ // ), 0
+ // );
+
+ // Remove the uncorrected linear coefficients
+ fit_[faci][0] -= w[faci];
+ fit_[faci][1] -= 1 - w[faci];
+
+ // if (limiter < 0.99)
+ // {
+ // for(label i = 0; i < stencilSize; i++)
+ // {
+ // fit_[faci][i] *= limiter;
+ // }
+ // }
+ }
+ else
+ {
+ Pout<< "Could not fit face " << faci
+ << " " << fit_[faci][0] << " " << w[faci]
+ << " " << fit_[faci][1] << " " << 1 - w[faci]<< endl;
+ fit_[faci] = 0;
+ }
return minSize_ - nSVDzeros;
}
diff --git a/src/postProcessing/fieldAverage/fieldAverageItem/fieldAverageItem.H b/src/postProcessing/fieldAverage/fieldAverageItem/fieldAverageItem.H
index 1bcb5689bbb45027b259e28d1a4274c08b5fa867..a8b8830016fc84d82b6d13b907582b80d96e3cac 100644
--- a/src/postProcessing/fieldAverage/fieldAverageItem/fieldAverageItem.H
+++ b/src/postProcessing/fieldAverage/fieldAverageItem/fieldAverageItem.H
@@ -102,15 +102,6 @@ private:
baseType base_;
- // Private Member Functions
-
- //- Disallow default bitwise copy construct
-// fieldAverageItem(const fieldAverageItem&);
-
- //- Disallow default bitwise assignment
-// void operator=(const fieldAverageItem&);
-
-
public:
// Constructors
@@ -175,6 +166,29 @@ public:
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
diff --git a/src/thermophysicalModels/basic/hThermo/hThermo.C b/src/thermophysicalModels/basic/hThermo/hThermo.C
index 19e43d2410876416c43b7e93c1f3a546799b5cb6..7e32f7a0636288169e7bc36f02a3b9140c20efba 100644
--- a/src/thermophysicalModels/basic/hThermo/hThermo.C
+++ b/src/thermophysicalModels/basic/hThermo/hThermo.C
@@ -28,15 +28,10 @@ License
#include "fvMesh.H"
#include "fixedValueFvPatchFields.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class MixtureType>
-hThermo<MixtureType>::hThermo(const fvMesh& mesh)
+Foam::hThermo<MixtureType>::hThermo(const fvMesh& mesh)
:
basicThermo(mesh),
MixtureType(*this, mesh),
@@ -56,9 +51,12 @@ hThermo<MixtureType>::hThermo(const fvMesh& mesh)
hBoundaryTypes()
)
{
- forAll(h_, celli)
+ scalarField& hCells = h_.internalField();
+ const scalarField& TCells = T_.internalField();
+
+ forAll(hCells, celli)
{
- h_[celli] = this->cellMixture(celli).H(T_[celli]);
+ hCells[celli] = this->cellMixture(celli).H(TCells[celli]);
}
forAll(h_.boundaryField(), patchi)
@@ -76,25 +74,33 @@ hThermo<MixtureType>::hThermo(const fvMesh& mesh)
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class MixtureType>
-hThermo<MixtureType>::~hThermo()
+Foam::hThermo<MixtureType>::~hThermo()
{}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class MixtureType>
-void hThermo<MixtureType>::calculate()
+void Foam::hThermo<MixtureType>::calculate()
{
- forAll(T_, celli)
+ const scalarField& hCells = h_.internalField();
+ const scalarField& pCells = p_.internalField();
+
+ scalarField& TCells = T_.internalField();
+ scalarField& psiCells = psi_.internalField();
+ scalarField& muCells = mu_.internalField();
+ scalarField& alphaCells = alpha_.internalField();
+
+ forAll(TCells, celli)
{
const typename MixtureType::thermoType& mixture_ =
this->cellMixture(celli);
- T_[celli] = mixture_.TH(h_[celli], T_[celli]);
- psi_[celli] = mixture_.psi(p_[celli], T_[celli]);
+ TCells[celli] = mixture_.TH(hCells[celli], TCells[celli]);
+ psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
- mu_[celli] = mixture_.mu(T_[celli]);
- alpha_[celli] = mixture_.alpha(T_[celli]);
+ muCells[celli] = mixture_.mu(TCells[celli]);
+ alphaCells[celli] = mixture_.alpha(TCells[celli]);
}
forAll(T_.boundaryField(), patchi)
@@ -143,7 +149,7 @@ void hThermo<MixtureType>::calculate()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class MixtureType>
-void hThermo<MixtureType>::correct()
+void Foam::hThermo<MixtureType>::correct()
{
if (debug)
{
@@ -163,7 +169,7 @@ void hThermo<MixtureType>::correct()
template<class MixtureType>
-tmp<scalarField> hThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::h
(
const scalarField& T,
const labelList& cells
@@ -182,7 +188,7 @@ tmp<scalarField> hThermo<MixtureType>::h
template<class MixtureType>
-tmp<scalarField> hThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::h
(
const scalarField& T,
const label patchi
@@ -201,7 +207,7 @@ tmp<scalarField> hThermo<MixtureType>::h
template<class MixtureType>
-tmp<scalarField> hThermo<MixtureType>::Cp
+Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::Cp
(
const scalarField& T,
const label patchi
@@ -220,7 +226,7 @@ tmp<scalarField> hThermo<MixtureType>::Cp
template<class MixtureType>
-tmp<volScalarField> hThermo<MixtureType>::Cp() const
+Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cp() const
{
const fvMesh& mesh = T_.mesh();
@@ -258,7 +264,7 @@ tmp<volScalarField> hThermo<MixtureType>::Cp() const
template<class MixtureType>
-tmp<volScalarField> hThermo<MixtureType>::Cv() const
+Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cv() const
{
const fvMesh& mesh = T_.mesh();
@@ -303,7 +309,7 @@ tmp<volScalarField> hThermo<MixtureType>::Cv() const
template<class MixtureType>
-bool hThermo<MixtureType>::read()
+bool Foam::hThermo<MixtureType>::read()
{
if (basicThermo::read())
{
@@ -317,8 +323,4 @@ bool hThermo<MixtureType>::read()
}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
// ************************************************************************* //
diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C
index 2d72a4f6887e5b03f7df977ed0ceb1bde5e9a8b9..064e7c3ebafeb0f67be03f7d547df9e1d2c31f55 100644
--- a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C
+++ b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C
@@ -28,22 +28,20 @@ License
#include "fvMesh.H"
#include "fixedValueFvPatchFields.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class MixtureType>
-hMixtureThermo<MixtureType>::hMixtureThermo(const fvMesh& mesh)
+Foam::hMixtureThermo<MixtureType>::hMixtureThermo(const fvMesh& mesh)
:
hCombustionThermo(mesh),
MixtureType(*this, mesh)
{
- forAll(h_, celli)
+ scalarField& hCells = h_.internalField();
+ const scalarField& TCells = T_.internalField();
+
+ forAll(hCells, celli)
{
- h_[celli] = this->cellMixture(celli).H(T_[celli]);
+ hCells[celli] = this->cellMixture(celli).H(TCells[celli]);
}
forAll(h_.boundaryField(), patchi)
@@ -61,25 +59,33 @@ hMixtureThermo<MixtureType>::hMixtureThermo(const fvMesh& mesh)
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class MixtureType>
-hMixtureThermo<MixtureType>::~hMixtureThermo()
+Foam::hMixtureThermo<MixtureType>::~hMixtureThermo()
{}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class MixtureType>
-void hMixtureThermo<MixtureType>::calculate()
+void Foam::hMixtureThermo<MixtureType>::calculate()
{
- forAll(T_, celli)
+ const scalarField& hCells = h_.internalField();
+ const scalarField& pCells = p_.internalField();
+
+ scalarField& TCells = T_.internalField();
+ scalarField& psiCells = psi_.internalField();
+ scalarField& muCells = mu_.internalField();
+ scalarField& alphaCells = alpha_.internalField();
+
+ forAll(TCells, celli)
{
const typename MixtureType::thermoType& mixture_ =
this->cellMixture(celli);
- T_[celli] = mixture_.TH(h_[celli], T_[celli]);
- psi_[celli] = mixture_.psi(p_[celli], T_[celli]);
+ TCells[celli] = mixture_.TH(hCells[celli], TCells[celli]);
+ psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
- mu_[celli] = mixture_.mu(T_[celli]);
- alpha_[celli] = mixture_.alpha(T_[celli]);
+ muCells[celli] = mixture_.mu(TCells[celli]);
+ alphaCells[celli] = mixture_.alpha(TCells[celli]);
}
forAll(T_.boundaryField(), patchi)
@@ -128,7 +134,7 @@ void hMixtureThermo<MixtureType>::calculate()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class MixtureType>
-tmp<scalarField> hMixtureThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::h
(
const scalarField& T,
const labelList& cells
@@ -147,7 +153,7 @@ tmp<scalarField> hMixtureThermo<MixtureType>::h
template<class MixtureType>
-tmp<scalarField> hMixtureThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::h
(
const scalarField& T,
const label patchi
@@ -166,7 +172,7 @@ tmp<scalarField> hMixtureThermo<MixtureType>::h
template<class MixtureType>
-tmp<scalarField> hMixtureThermo<MixtureType>::Cp
+Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::Cp
(
const scalarField& T,
const label patchi
@@ -186,7 +192,7 @@ tmp<scalarField> hMixtureThermo<MixtureType>::Cp
template<class MixtureType>
-tmp<volScalarField> hMixtureThermo<MixtureType>::Cp() const
+Foam::tmp<Foam::volScalarField> Foam::hMixtureThermo<MixtureType>::Cp() const
{
const fvMesh& mesh = T_.mesh();
@@ -224,7 +230,7 @@ tmp<volScalarField> hMixtureThermo<MixtureType>::Cp() const
template<class MixtureType>
-void hMixtureThermo<MixtureType>::correct()
+void Foam::hMixtureThermo<MixtureType>::correct()
{
if (debug)
{
@@ -244,7 +250,7 @@ void hMixtureThermo<MixtureType>::correct()
template<class MixtureType>
-bool hMixtureThermo<MixtureType>::read()
+bool Foam::hMixtureThermo<MixtureType>::read()
{
if (hCombustionThermo::read())
{
@@ -258,8 +264,4 @@ bool hMixtureThermo<MixtureType>::read()
}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
// ************************************************************************* //
diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
index 04bf08e33027859a4b216fc6638cdf8465e40f0a..0086e8cd147c2824097599d087c7c6537e0787b3 100644
--- a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
+++ b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
@@ -28,23 +28,23 @@ License
#include "fvMesh.H"
#include "fixedValueFvPatchFields.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class MixtureType>
-hhuMixtureThermo<MixtureType>::hhuMixtureThermo(const fvMesh& mesh)
+Foam::hhuMixtureThermo<MixtureType>::hhuMixtureThermo(const fvMesh& mesh)
:
hhuCombustionThermo(mesh),
MixtureType(*this, mesh)
{
- forAll(h_, celli)
+ scalarField& hCells = h_.internalField();
+ scalarField& huCells = hu_.internalField();
+ const scalarField& TCells = T_.internalField();
+ const scalarField& TuCells = Tu_.internalField();
+
+ forAll(hCells, celli)
{
- h_[celli] = this->cellMixture(celli).H(T_[celli]);
- hu_[celli] = this->cellReactants(celli).H(Tu_[celli]);
+ hCells[celli] = this->cellMixture(celli).H(TCells[celli]);
+ huCells[celli] = this->cellReactants(celli).H(TuCells[celli]);
}
forAll(h_.boundaryField(), patchi)
@@ -71,27 +71,38 @@ hhuMixtureThermo<MixtureType>::hhuMixtureThermo(const fvMesh& mesh)
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class MixtureType>
-hhuMixtureThermo<MixtureType>::~hhuMixtureThermo()
+Foam::hhuMixtureThermo<MixtureType>::~hhuMixtureThermo()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class MixtureType>
-void hhuMixtureThermo<MixtureType>::calculate()
+void Foam::hhuMixtureThermo<MixtureType>::calculate()
{
- forAll(T_, celli)
+ const scalarField& hCells = h_.internalField();
+ const scalarField& huCells = hu_.internalField();
+ const scalarField& pCells = p_.internalField();
+
+ scalarField& TCells = T_.internalField();
+ scalarField& TuCells = Tu_.internalField();
+ scalarField& psiCells = psi_.internalField();
+ scalarField& muCells = mu_.internalField();
+ scalarField& alphaCells = alpha_.internalField();
+
+ forAll(TCells, celli)
{
- const typename MixtureType::thermoType& mixture_ =
+ const typename MixtureType::thermoType& mixture_ =
this->cellMixture(celli);
- T_[celli] = mixture_.TH(h_[celli], T_[celli]);
- psi_[celli] = mixture_.psi(p_[celli], T_[celli]);
+ TCells[celli] = mixture_.TH(hCells[celli], TCells[celli]);
+ psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
- mu_[celli] = mixture_.mu(T_[celli]);
- alpha_[celli] = mixture_.alpha(T_[celli]);
+ muCells[celli] = mixture_.mu(TCells[celli]);
+ alphaCells[celli] = mixture_.alpha(TCells[celli]);
- Tu_[celli] = this->cellReactants(celli).TH(hu_[celli], Tu_[celli]);
+ TuCells[celli] =
+ this->cellReactants(celli).TH(huCells[celli], TuCells[celli]);
}
forAll(T_.boundaryField(), patchi)
@@ -144,7 +155,7 @@ void hhuMixtureThermo<MixtureType>::calculate()
template<class MixtureType>
-void hhuMixtureThermo<MixtureType>::correct()
+void Foam::hhuMixtureThermo<MixtureType>::correct()
{
if (debug)
{
@@ -163,7 +174,7 @@ void hhuMixtureThermo<MixtureType>::correct()
}
template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::h
(
const scalarField& T,
const labelList& cells
@@ -182,7 +193,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::h
template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::h
(
const scalarField& T,
const label patchi
@@ -201,7 +212,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::h
template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::Cp
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::Cp
(
const scalarField& T,
const label patchi
@@ -221,7 +232,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::Cp
template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::Cp() const
+Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Cp() const
{
const fvMesh& mesh = T_.mesh();
@@ -260,7 +271,7 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::Cp() const
template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::hu
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::hu
(
const scalarField& Tu,
const labelList& cells
@@ -279,7 +290,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::hu
template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::hu
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::hu
(
const scalarField& Tu,
const label patchi
@@ -298,7 +309,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::hu
template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::Tb() const
+Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Tb() const
{
tmp<volScalarField> tTb
(
@@ -342,7 +353,8 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::Tb() const
template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::psiu() const
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::psiu() const
{
tmp<volScalarField> tpsiu
(
@@ -388,7 +400,8 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::psiu() const
template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::psib() const
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::psib() const
{
tmp<volScalarField> tpsib
(
@@ -426,7 +439,8 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::psib() const
forAll(ppsib, facei)
{
ppsib[facei] =
- this->patchFaceReactants(patchi, facei).psi(pp[facei], pTb[facei]);
+ this->patchFaceReactants
+ (patchi, facei).psi(pp[facei], pTb[facei]);
}
}
@@ -435,7 +449,7 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::psib() const
template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::muu() const
+Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::muu() const
{
tmp<volScalarField> tmuu
(
@@ -478,7 +492,7 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::muu() const
template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::mub() const
+Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::mub() const
{
tmp<volScalarField> tmub
(
@@ -521,7 +535,7 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::mub() const
template<class MixtureType>
-bool hhuMixtureThermo<MixtureType>::read()
+bool Foam::hhuMixtureThermo<MixtureType>::read()
{
if (hhuCombustionThermo::read())
{
@@ -535,8 +549,4 @@ bool hhuMixtureThermo<MixtureType>::read()
}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
// ************************************************************************* //
diff --git a/tutorials/channelOodles/channel395/constant/postChannelDict b/tutorials/channelOodles/channel395/constant/postChannelDict
index e69c0386ef95781462a1b01424c2d8eee53f1e2c..d63cd656675c6ebe4aecf1608dca9b561a5a3fe8 100644
--- a/tutorials/channelOodles/channel395/constant/postChannelDict
+++ b/tutorials/channelOodles/channel395/constant/postChannelDict
@@ -15,14 +15,14 @@ FoamFile
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-Nx 40;
-Ny
-(
- 25
- 25
-);
-Nz 30;
+// Seed patches to start layering from
+patches (bottomWall);
+// 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;
// ************************************************************************* //