ENH: checkMesh: output vol fields of mesh quality. Fixes #466.

applications/utilities/mesh/manipulation/checkMesh/Make/files

@@ -1,3 +1,4 @@
+writeFields.C
 checkTools.C
 checkTopology.C
 checkGeometry.C

applications/utilities/mesh/manipulation/checkMesh/checkMesh.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,12 +57,18 @@ Usage
     Reconstruct all cellSets and faceSets geometry and write to postProcessing/
     directory according to surfaceFormat (e.g. vtk or ensight)
 
+    \param -writeAllFields \n
+    Writes all mesh quality measures as fields.
+
+    \param -writeFields '(\<fieldName\>)' \n
+    Writes selected mesh quality measures as fields.
+
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
 #include "timeSelector.H"
 #include "Time.H"
-#include "polyMesh.H"
+#include "fvMesh.H"
 #include "globalMeshData.H"
 #include "surfaceWriter.H"
 #include "vtkSetWriter.H"
@@ -71,6 +77,7 @@ Usage
 #include "checkTopology.H"
 #include "checkGeometry.H"
 #include "checkMeshQuality.H"
+#include "writeFields.H"
 
 using namespace Foam;
 
@@ -96,6 +103,17 @@ int main(int argc, char *argv[])
         "include extra topology checks"
     );
     argList::addBoolOption
+    (
+        "writeAllFields",
+        "write volFields with mesh quality parameters"
+    );
+    argList::addOption
+    (
+        "writeFields",
+        "wordList"
+        "write volFields with selected mesh quality parameters"
+    );
+    argList::addBoolOption
     (
         "meshQuality",
         "read user-defined mesh quality criterions from system/meshQualityDict"
@@ -110,7 +128,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     instantList timeDirs = timeSelector::select0(runTime, args);
-    #include "createNamedPolyMesh.H"
+    #include "createNamedMesh.H"
 
     const bool noTopology  = args.optionFound("noTopology");
     const bool allGeometry = args.optionFound("allGeometry");
@@ -119,6 +137,24 @@ int main(int argc, char *argv[])
 
     word surfaceFormat;
     const bool writeSets = args.optionReadIfPresent("writeSets", surfaceFormat);
+    HashSet<word> selectedFields;
+    bool writeFields = args.optionReadIfPresent
+    (
+        "writeFields",
+        selectedFields
+    );
+    if (!writeFields && args.optionFound("writeAllFields"))
+    {
+        selectedFields.insert("nonOrthoAngle");
+        selectedFields.insert("faceWeight");
+        selectedFields.insert("skewness");
+        selectedFields.insert("cellDeterminant");
+        selectedFields.insert("aspectRatio");
+        selectedFields.insert("cellShapes");
+        selectedFields.insert("cellVolume");
+        selectedFields.insert("cellVolumeRatio");
+    }
+
 
     if (noTopology)
     {
@@ -143,6 +179,11 @@ int main(int argc, char *argv[])
             << " representation"
             << " of all faceSets and cellSets." << nl << endl;
     }
+    if (selectedFields.size())
+    {
+        Info<< "Writing mesh quality as fields " << selectedFields << nl
+            << endl;
+    }
 
 
     autoPtr<IOdictionary> qualDict;
@@ -234,6 +275,10 @@ int main(int argc, char *argv[])
                 Info<< "\nFailed " << nFailedChecks << " mesh checks.\n"
                     << endl;
             }
+
+
+            // Write selected fields
+            Foam::writeFields(mesh, selectedFields);
         }
         else if (state == polyMesh::POINTS_MOVED)
         {
@@ -262,6 +307,10 @@ int main(int argc, char *argv[])
             {
                 Info<< "\nMesh OK.\n" << endl;
             }
+
+
+            // Write selected fields
+            Foam::writeFields(mesh, selectedFields);
         }
     }
 

applications/utilities/mesh/manipulation/checkMesh/writeFields.C

@@ -0,0 +1,360 @@
+#include "writeFields.H"
+#include "volFields.H"
+#include "polyMeshTools.H"
+#include "zeroGradientFvPatchFields.H"
+#include "syncTools.H"
+
+using namespace Foam;
+
+void maxFaceToCell
+(
+    const scalarField& faceData,
+    volScalarField& cellData
+)
+{
+    const cellList& cells = cellData.mesh().cells();
+
+    scalarField& cellFld = cellData.ref();
+
+    cellFld = -GREAT;
+    forAll(cells, cellI)
+    {
+        const cell& cFaces = cells[cellI];
+        forAll(cFaces, i)
+        {
+            cellFld[cellI] = max(cellFld[cellI], faceData[cFaces[i]]);
+        }
+    }
+
+    forAll(cellData.boundaryField(), patchI)
+    {
+        fvPatchScalarField& fvp = cellData.boundaryFieldRef()[patchI];
+
+        fvp = fvp.patch().patchSlice(faceData);
+    }
+    cellData.correctBoundaryConditions();
+}
+
+
+void minFaceToCell
+(
+    const scalarField& faceData,
+    volScalarField& cellData
+)
+{
+    const cellList& cells = cellData.mesh().cells();
+
+    scalarField& cellFld = cellData.ref();
+
+    cellFld = GREAT;
+    forAll(cells, cellI)
+    {
+        const cell& cFaces = cells[cellI];
+        forAll(cFaces, i)
+        {
+            cellFld[cellI] = min(cellFld[cellI], faceData[cFaces[i]]);
+        }
+    }
+
+    forAll(cellData.boundaryField(), patchI)
+    {
+        fvPatchScalarField& fvp = cellData.boundaryFieldRef()[patchI];
+
+        fvp = fvp.patch().patchSlice(faceData);
+    }
+    cellData.correctBoundaryConditions();
+}
+
+
+void Foam::writeFields
+(
+    const fvMesh& mesh,
+    const HashSet<word>& selectedFields
+)
+{
+    if (selectedFields.empty())
+    {
+        return;
+    }
+
+    Info<< "Writing fields with mesh quality parameters" << endl;
+
+    if (selectedFields.found("nonOrthoAngle"))
+    {
+        //- Face based orthogonality
+        const scalarField faceOrthogonality
+        (
+            polyMeshTools::faceOrthogonality
+            (
+                mesh,
+                mesh.faceAreas(),
+                mesh.cellCentres()
+            )
+        );
+
+        //- Face based angle
+        const scalarField nonOrthoAngle
+        (
+            radToDeg
+            (
+                Foam::acos(min(1.0, faceOrthogonality))
+            )
+        );
+
+        //- Cell field - max of either face
+        volScalarField cellNonOrthoAngle
+        (
+            IOobject
+            (
+                "nonOrthoAngle",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("angle", dimless, 0),
+            calculatedFvPatchScalarField::typeName
+        );
+        //- Take max
+        maxFaceToCell(nonOrthoAngle, cellNonOrthoAngle);
+        Info<< "    Writing non-orthogonality (angle) to "
+            << cellNonOrthoAngle.name() << endl;
+        cellNonOrthoAngle.write();
+    }
+
+    if (selectedFields.found("faceWeight"))
+    {
+        const scalarField faceWeights
+        (
+            polyMeshTools::faceWeights
+            (
+                mesh,
+                mesh.faceCentres(),
+                mesh.faceAreas(),
+                mesh.cellCentres()
+            )
+        );
+
+        volScalarField cellWeights
+        (
+            IOobject
+            (
+                "faceWeight",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("weight", dimless, 0),
+            calculatedFvPatchScalarField::typeName
+        );
+        //- Take min
+        minFaceToCell(faceWeights, cellWeights);
+        Info<< "    Writing face interpolation weights (0..0.5) to "
+            << cellWeights.name() << endl;
+        cellWeights.write();
+    }
+
+
+    // Skewness
+    // ~~~~~~~~
+
+    if (selectedFields.found("skewness"))
+    {
+        //- Face based skewness
+        const scalarField faceSkewness
+        (
+            polyMeshTools::faceSkewness
+            (
+                mesh,
+                mesh.points(),
+                mesh.faceCentres(),
+                mesh.faceAreas(),
+                mesh.cellCentres()
+            )
+        );
+
+        //- Cell field - max of either face
+        volScalarField cellSkewness
+        (
+            IOobject
+            (
+                "skewness",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("skewness", dimless, 0),
+            calculatedFvPatchScalarField::typeName
+        );
+        //- Take max
+        maxFaceToCell(faceSkewness, cellSkewness);
+        Info<< "    Writing face skewness to " << cellSkewness.name() << endl;
+        cellSkewness.write();
+    }
+
+
+    // cellDeterminant
+    // ~~~~~~~~~~~~~~~
+
+    if (selectedFields.found("cellDeterminant"))
+    {
+        volScalarField cellDeterminant
+        (
+            IOobject
+            (
+                "cellDeterminant",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("cellDeterminant", dimless, 0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        cellDeterminant.primitiveFieldRef() =
+            primitiveMeshTools::cellDeterminant
+            (
+                mesh,
+                mesh.geometricD(),
+                mesh.faceAreas(),
+                syncTools::getInternalOrCoupledFaces(mesh)
+            );
+        cellDeterminant.correctBoundaryConditions();
+        Info<< "    Writing cell determinant to "
+            << cellDeterminant.name() << endl;
+        cellDeterminant.write();
+    }
+
+
+    // Aspect ratio
+    // ~~~~~~~~~~~~
+    if (selectedFields.found("aspectRatio"))
+    {
+        volScalarField aspectRatio
+        (
+            IOobject
+            (
+                "aspectRatio",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("aspectRatio", dimless, 0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+
+
+        scalarField cellOpenness;
+        polyMeshTools::cellClosedness
+        (
+            mesh,
+            mesh.geometricD(),
+            mesh.faceAreas(),
+            mesh.cellVolumes(),
+            cellOpenness,
+            aspectRatio.ref()
+        );
+
+        aspectRatio.correctBoundaryConditions();
+        Info<< "    Writing aspect ratio to " << aspectRatio.name() << endl;
+        aspectRatio.write();
+    }
+
+
+    // cell type
+    // ~~~~~~~~~
+    if (selectedFields.found("cellShapes"))
+    {
+        volScalarField shape
+        (
+            IOobject
+            (
+                "cellShapes",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("cellShapes", dimless, 0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        const cellShapeList& cellShapes = mesh.cellShapes();
+        forAll(cellShapes, cellI)
+        {
+            const cellModel& model = cellShapes[cellI].model();
+            shape[cellI] = model.index();
+        }
+        shape.correctBoundaryConditions();
+        Info<< "    Writing cell shape (hex, tet etc.) to " << shape.name()
+            << endl;
+        shape.write();
+    }
+
+    if (selectedFields.found("cellVolume"))
+    {
+        volScalarField V
+        (
+            IOobject
+            (
+                "cellVolume",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("cellVolume", dimVolume, 0),
+            calculatedFvPatchScalarField::typeName
+        );
+        V.ref() = mesh.V();
+        Info<< "    Writing cell volume to " << V.name() << endl;
+        V.write();
+    }
+
+    if (selectedFields.found("cellVolumeRatio"))
+    {
+        const scalarField faceVolumeRatio
+        (
+            polyMeshTools::volRatio
+            (
+                mesh,
+                mesh.V()
+            )
+        );
+
+        volScalarField cellVolumeRatio
+        (
+            IOobject
+            (
+                "cellVolumeRatio",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("cellVolumeRatio", dimless, 0),
+            calculatedFvPatchScalarField::typeName
+        );
+        //- Take min
+        minFaceToCell(faceVolumeRatio, cellVolumeRatio);
+        Info<< "    Writing cell volume ratio to "
+            << cellVolumeRatio.name() << endl;
+        cellVolumeRatio.write();
+    }
+
+    Info<< endl;
+}

applications/utilities/mesh/manipulation/checkMesh/writeFields.H

@@ -0,0 +1,10 @@
+#include "fvMesh.H"
+
+namespace Foam
+{
+    void writeFields
+    (
+        const fvMesh&,
+        const HashSet<word>& selectedFields
+    );
+}