Commit cec4b184 authored by Mark Olesen's avatar Mark Olesen
Browse files

Merge remote branch 'OpenCFD/master' into olesenm

parents f1a0179d 553b439f
......@@ -327,14 +327,17 @@ meshQualityControls
// Set to 180 to disable.
maxConcave 80;
//- Minimum projected area v.s. actual area. Set to -1 to disable.
minFlatness 0.5;
//- Minimum pyramid volume. Is absolute volume of cell pyramid.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minVol 1e-13;
//- Minimum tet volume. Is absolute volume of the tet formed by the
// face-centre decomposition triangle and the cell centre.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minTetVol 1e-20;
//- Minimum face area. Set to <0 to disable.
minArea -1;
......
......@@ -159,6 +159,25 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
}
}
{
faceSet faces(mesh, "wrongOrientedTriangleFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFaceTets(true, 0, &faces))
{
noFailedChecks++;
label nFaces = returnReduce(faces.size(), sumOp<label>());
if (nFaces > 0)
{
Info<< " <<Writing " << nFaces
<< " faces with incorrectly orientated triangles to set "
<< faces.name() << endl;
faces.instance() = mesh.pointsInstance();
faces.write();
}
}
}
{
faceSet faces(mesh, "skewFaces", mesh.nFaces()/100 + 1);
if (mesh.checkFaceSkewness(true, &faces))
......
......@@ -22,6 +22,9 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
splitMeshRegions
Description
Splits mesh into multiple regions.
......@@ -32,6 +35,7 @@ Description
- any face inbetween differing cellZones (-cellZones)
Output is:
- volScalarField with regions as different scalars (-detectOnly) or
- mesh with multiple regions or
- mesh with cells put into cellZones (-makeCellZones)
......@@ -47,7 +51,13 @@ Description
- useCellZonesOnly does not do a walk and uses the cellZones only. Use
this if you don't mind having disconnected domains in a single region.
This option requires all cells to be in one (and one only) cellZone.
- writes maps like decomposePar back to original mesh:
- pointRegionAddressing : for every point in this region the point in
the original mesh
- cellRegionAddressing : ,, cell ,, cell ,,
- faceRegionAddressing : ,, face ,, face in
the original mesh + 'turning index'. For a face in the same orientation
this is the original facelabel+1, for a turned face this is -facelabel-1
\*---------------------------------------------------------------------------*/
#include "SortableList.H"
......
......@@ -334,6 +334,14 @@ int main(int argc, char *argv[])
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
PtrList<surfaceScalarField> surfaceScalarFields;
readFields(mesh, objects, surfaceScalarFields);
PtrList<surfaceVectorField> surfaceVectorFields;
readFields(mesh, objects, surfaceVectorFields);
PtrList<surfaceSphericalTensorField> surfaceSphericalTensorFields;
readFields(mesh, objects, surfaceSphericalTensorFields);
PtrList<surfaceSymmTensorField> surfaceSymmTensorFields;
readFields(mesh, objects, surfaceSymmTensorFields);
PtrList<surfaceTensorField> surfaceTensorFields;
readFields(mesh, objects, surfaceTensorFields);
// Construct the point fields
......@@ -619,6 +627,10 @@ int main(int argc, char *argv[])
|| volSymmTensorFields.size()
|| volTensorFields.size()
|| surfaceScalarFields.size()
|| surfaceVectorFields.size()
|| surfaceSphericalTensorFields.size()
|| surfaceSymmTensorFields.size()
|| surfaceTensorFields.size()
)
{
labelIOList faceProcAddressing
......@@ -650,6 +662,10 @@ int main(int argc, char *argv[])
fieldDecomposer.decomposeFields(volTensorFields);
fieldDecomposer.decomposeFields(surfaceScalarFields);
fieldDecomposer.decomposeFields(surfaceVectorFields);
fieldDecomposer.decomposeFields(surfaceSphericalTensorFields);
fieldDecomposer.decomposeFields(surfaceSymmTensorFields);
fieldDecomposer.decomposeFields(surfaceTensorFields);
}
......
......@@ -191,6 +191,10 @@ int main(int argc, char *argv[])
|| objects.lookupClass(volSymmTensorField::typeName).size()
|| objects.lookupClass(volTensorField::typeName).size()
|| objects.lookupClass(surfaceScalarField::typeName).size()
|| objects.lookupClass(surfaceVectorField::typeName).size()
|| objects.lookupClass(surfaceSphericalTensorField::typeName).size()
|| objects.lookupClass(surfaceSymmTensorField::typeName).size()
|| objects.lookupClass(surfaceTensorField::typeName).size()
)
{
Info<< "Reconstructing FV fields" << nl << endl;
......@@ -235,6 +239,26 @@ int main(int argc, char *argv[])
objects,
selectedFields
);
fvReconstructor.reconstructFvSurfaceFields<vector>
(
objects,
selectedFields
);
fvReconstructor.reconstructFvSurfaceFields<sphericalTensor>
(
objects,
selectedFields
);
fvReconstructor.reconstructFvSurfaceFields<symmTensor>
(
objects,
selectedFields
);
fvReconstructor.reconstructFvSurfaceFields<tensor>
(
objects,
selectedFields
);
}
else
{
......
......@@ -61,8 +61,7 @@ Foam::polyBoundaryMesh::polyBoundaryMesh
:
polyPatchList(),
regIOobject(io),
mesh_(mesh),
neighbourEdgesPtr_(NULL)
mesh_(mesh)
{
if (readOpt() == IOobject::MUST_READ)
{
......@@ -110,17 +109,14 @@ Foam::polyBoundaryMesh::polyBoundaryMesh
:
polyPatchList(size),
regIOobject(io),
mesh_(pm),
neighbourEdgesPtr_(NULL)
mesh_(pm)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::polyBoundaryMesh::~polyBoundaryMesh()
{
deleteDemandDrivenData(neighbourEdgesPtr_);
}
{}
void Foam::polyBoundaryMesh::clearGeom()
......@@ -134,7 +130,8 @@ void Foam::polyBoundaryMesh::clearGeom()
void Foam::polyBoundaryMesh::clearAddressing()
{
deleteDemandDrivenData(neighbourEdgesPtr_);
neighbourEdgesPtr_.clear();
patchIDPtr_.clear();
forAll (*this, patchi)
{
......@@ -201,10 +198,10 @@ Foam::polyBoundaryMesh::neighbourEdges() const
<< " boundaries." << endl;
}
if (!neighbourEdgesPtr_)
if (!neighbourEdgesPtr_.valid())
{
neighbourEdgesPtr_ = new List<labelPairList>(size());
List<labelPairList>& neighbourEdges = *neighbourEdgesPtr_;
neighbourEdgesPtr_.reset(new List<labelPairList>(size()));
List<labelPairList>& neighbourEdges = neighbourEdgesPtr_();
// Initialize.
label nEdgePairs = 0;
......@@ -320,7 +317,36 @@ Foam::polyBoundaryMesh::neighbourEdges() const
}
}
return *neighbourEdgesPtr_;
return neighbourEdgesPtr_();
}
const Foam::labelList& Foam::polyBoundaryMesh::patchID() const
{
if (!patchIDPtr_.valid())
{
patchIDPtr_.reset
(
new labelList
(
mesh_.nFaces()
- mesh_.nInternalFaces()
)
);
labelList& patchID = patchIDPtr_();
const polyBoundaryMesh& bm = *this;
forAll(bm, patchI)
{
label bFaceI = bm[patchI].start() - mesh_.nInternalFaces();
forAll(bm[patchI], i)
{
patchID[bFaceI++] = patchI;
}
}
}
return patchIDPtr_();
}
......@@ -654,7 +680,8 @@ void Foam::polyBoundaryMesh::movePoints(const pointField& p)
void Foam::polyBoundaryMesh::updateMesh()
{
deleteDemandDrivenData(neighbourEdgesPtr_);
neighbourEdgesPtr_.clear();
patchIDPtr_.clear();
PstreamBuffers pBufs(Pstream::defaultCommsType);
......
......@@ -67,8 +67,10 @@ class polyBoundaryMesh
//- Reference to mesh
const polyMesh& mesh_;
mutable autoPtr<labelList> patchIDPtr_;
//- Edges of neighbouring patches
mutable List<labelPairList>* neighbourEdgesPtr_;
mutable autoPtr<List<labelPairList> > neighbourEdgesPtr_;
// Private Member Functions
......@@ -157,6 +159,9 @@ public:
//- Return patch index for a given face label
label whichPatch(const label faceIndex) const;
//- Per boundary face label the patch index
const labelList& patchID() const;
//- Return the set of patch IDs corresponding to the given list of names
// Wild cards are expanded.
labelHashSet patchSet(const wordList&) const;
......
......@@ -556,7 +556,6 @@ public:
//- Check for negative cell volumes
bool checkCellVolumes
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
......@@ -576,6 +575,14 @@ public:
labelHashSet* setPtr = NULL
) const;
//- Check face-decomposition tet volume
bool checkFaceTets
(
const bool report = false,
const scalar minTetVol = 0,
labelHashSet* setPtr = NULL
) const;
//- Check face skewness
bool checkFaceSkewness
(
......
......@@ -26,6 +26,7 @@ License
#include "primitiveMesh.H"
#include "pyramidPointFaceRef.H"
#include "tetPointRef.H"
#include "ListOps.H"
#include "unitConversion.H"
#include "SortableList.H"
......@@ -401,7 +402,7 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
<< "checking mesh non-orthogonality" << endl;
}
// for all internal faces check theat the d dot S product is positive
// for all internal faces check that the d dot S product is positive
const vectorField& centres = cellCentres();
const vectorField& areas = faceAreas();
......@@ -591,6 +592,114 @@ bool Foam::primitiveMesh::checkFacePyramids
}
bool Foam::primitiveMesh::checkFaceTets
(
const bool report,
const scalar minTetVol,
labelHashSet* setPtr
) const
{
if (debug)
{
Info<< "bool primitiveMesh::checkFaceTets("
<< "const bool, const scalar, labelHashSet*) const: "
<< "checking face orientation" << endl;
}
// check whether face area vector points to the cell with higher label
const vectorField& cc = cellCentres();
const vectorField& fc = faceCentres();
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const faceList& fcs = faces();
const pointField& p = points();
label nErrorPyrs = 0;
forAll (fcs, faceI)
{
// Create the owner tets - they will have negative volume
const face& f = fcs[faceI];
forAll(f, fp)
{
scalar tetVol = tetPointRef
(
p[f[fp]],
p[f.nextLabel(fp)],
fc[faceI],
cc[own[faceI]]
).mag();
if (tetVol > -minTetVol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorPyrs++;
break; // no need to check other tets
}
}
if (isInternalFace(faceI))
{
// Create the neighbour tet - it will have positive volume
const face& f = fcs[faceI];
forAll(f, fp)
{
scalar tetVol = tetPointRef
(
p[f[fp]],
p[f.nextLabel(fp)],
fc[faceI],
cc[nei[faceI]]
).mag();
if (tetVol < minTetVol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
nErrorPyrs++;
break;
}
}
}
}
reduce(nErrorPyrs, sumOp<label>());
if (nErrorPyrs > 0)
{
if (debug || report)
{
Info<< " ***Error in face tets: "
<< nErrorPyrs << " faces have incorrectly oriented face"
<< " decomposition triangles." << endl;
}
return true;
}
else
{
if (debug || report)
{
Info<< " Face tets OK." << endl;
}
return false;
}
}
bool Foam::primitiveMesh::checkFaceSkewness
(
const bool report,
......
......@@ -69,6 +69,10 @@ bool Foam::motionSmoother::checkMesh
(
readScalar(dict.lookup("minVol", true))
);
const scalar minTetVol
(
readScalar(dict.lookup("minTetVol", true))
);
const scalar maxConcave
(
readScalar(dict.lookup("maxConcave", true))
......@@ -105,7 +109,6 @@ bool Foam::motionSmoother::checkMesh
(
readScalar(dict.lookup("minDeterminant", true))
);
label nWrongFaces = 0;
Info<< "Checking faces in error :" << endl;
......@@ -158,6 +161,30 @@ bool Foam::motionSmoother::checkMesh
nWrongFaces = nNewWrongFaces;
}
if (minTetVol > -GREAT)
{
polyMeshGeometry::checkFaceTets
(
report,
minTetVol,
mesh,
mesh.cellCentres(),
mesh.faceCentres(),
mesh.points(),
checkFaces,
baffles,
&wrongFaces
);
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
Info<< " faces with face-decomposition tet volume < "
<< setw(5) << minTetVol << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
if (maxConcave < 180.0-SMALL)
{
polyMeshGeometry::checkFaceAngles
......@@ -417,6 +444,10 @@ bool Foam::motionSmoother::checkMesh
(
readScalar(dict.lookup("minVol", true))
);
const scalar minTetVol
(
readScalar(dict.lookup("minTetVol", true))
);
const scalar maxConcave
(
readScalar(dict.lookup("maxConcave", true))
......@@ -501,6 +532,27 @@ bool Foam::motionSmoother::checkMesh
nWrongFaces = nNewWrongFaces;
}
if (minTetVol > -GREAT)
{
meshGeom.checkFaceTets
(
report,
minTetVol,
meshGeom.mesh().points(),
checkFaces,
baffles,
&wrongFaces
);
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
Info<< " faces with face-decomposition tet volume < "
<< setw(5) << minTetVol << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
if (maxConcave < 180.0-SMALL)
{
meshGeom.checkFaceAngles
......
......@@ -26,6 +26,7 @@ License
#include "polyMeshGeometry.H"
#include "pyramidPointFaceRef.H"
#include "tetPointRef.H"
#include "syncTools.H"
#include "unitConversion.H"
......@@ -308,6 +309,59 @@ Foam::scalar Foam::polyMeshGeometry::calcSkewness
}
// Create the neighbour pyramid - it will have positive volume
bool Foam::polyMeshGeometry::checkFaceTet
(
const polyMesh& mesh,
const bool report,
const scalar minTetVol,
const pointField& p,
const label faceI,
const point& fc, // face centre
const point& cc, // cell centre
labelHashSet* setPtr
)
{
const face& f = mesh.faces()[faceI];
forAll(f, fp)
{
scalar tetVol = tetPointRef
(
p[f[fp]],
p[f.nextLabel(fp)],
fc,
cc
).mag();
if (tetVol < minTetVol)
{
if (report)
{
Pout<< "bool polyMeshGeometry::checkFaceTets("
<< "const bool, const scalar, const pointField&"
<< ", const pointField&, const labelList&,"
<< " labelHashSet*): "
<< "face " << faceI
<< " has a triangle that points the wrong way."
<< endl
<< "Tet volume: " << tetVol
<< " Face " << faceI
<< endl;
}
if (setPtr)
{
setPtr->insert(faceI);
}
return true;
}
}
return false;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
......@@ -719,6 +773,139 @@ bool Foam::polyMeshGeometry::checkFacePyramids
}
bool Foam::polyMeshGeometry::checkFaceTets
(
const bool report,
const scalar minTetVol,
const polyMesh& mesh,
const vectorField& cellCentres,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
const List&