Commit 68f65534 authored by Andrew Heather's avatar Andrew Heather
Browse files
parents 8293b56e b4263cca
......@@ -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];
......
......@@ -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));
......
......@@ -462,6 +462,7 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createPolyMesh.H"
scalar minLen(readScalar(IStringStream(args.additionalArgs()[0])()));
......
......@@ -439,6 +439,7 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createPolyMesh.H"
scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
......
......@@ -332,6 +332,7 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createPolyMesh.H"
bool overwrite = args.options().found("overwrite");
......
......@@ -56,6 +56,7 @@ int main(int argc, char *argv[])
argList::validArgs.append("cellSet");
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createMesh.H"
pointMesh pMesh(mesh);
......
......@@ -54,6 +54,7 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createPolyMesh.H"
word patchName(args.additionalArgs()[0]);
......
......@@ -53,6 +53,7 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createMesh.H"
bool overwrite = args.options().found("overwrite");
......
......@@ -532,6 +532,7 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createPolyMesh.H"
scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
......
......@@ -47,6 +47,7 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createPolyMesh.H"
scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
......
......@@ -346,6 +346,7 @@ int main(int argc, char *argv[])
# include "addTimeOptions.H"
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
bool patchFaces = args.options().found("patchFaces");
bool doCell = args.options().found("cell");
......
......@@ -61,6 +61,7 @@ int main(int argc, char *argv[])
argList::validOptions.insert("overwrite", "");
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createPolyMesh.H"
scalar thickness(readScalar(IStringStream(args.additionalArgs()[0])()));
......
......@@ -46,6 +46,7 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createPolyMesh.H"
bool overwrite = args.options().found("overwrite");
......
......@@ -75,6 +75,7 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
runTime.functionObjects().off();
# include "createPolyMesh.H"
Info<< "Mesh read in = "
......
printMeshStats.C
checkTopology.C
checkGeometry.C
checkMesh.C
EXE = $(FOAM_APPBIN)/checkMesh
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lmeshTools
#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;
}
#include "label.H"
namespace Foam
{
class polyMesh;
label checkGeometry
(
const polyMesh& mesh,
bool checkPointNearness,
bool checkCellDeterminant
);
}
#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)