View Issue Details

IDProjectCategoryView StatusLast Update
0001872OpenFOAMBugpublic2015-10-17 13:59
Reporterwyldckat Assigned Tohenry  
PrioritylowSeverityfeatureReproducibilitysometimes
Status resolvedResolutionfixed 
Product Versiondev 
Summary0001872: Automatic direction skipping in fvMatrixSolve can lead to confusion in case of 3D empty face
DescriptionAttached are two files for updating the ones at "src/OpenFOAM/include":

 - createMesh.H
 - createNamedMesh.H

The addition is a warning in case "mesh.nSolutionD()<2". This is as an attempt to avoid the user getting disoriented when the solver is not reporting any residuals for any one of the momentum direction equations, even if the "momentumPredictor" is turned on.

I don't suggest this to be checked on all "solve" calls because that would slow down things a few nano seconds more than necessary.
Steps To ReproduceOnce in a blue moon I get a weird looking behaviour from OpenFOAM's solvers: none of the momentum direction equations are solved.

Today this happened again. It feels something like getting "2 plus 2 equals 3".
After checking that "momentumPredictor" isn't used in "simpleFoam", that the "solverPerformance" debug options are turned on, that the boundary conditions are all correct, I went looking into the code... because clearly something in the matrix solver was making it quit for all directions.

The reason: https://github.com/OpenFOAM/OpenFOAM-2.4.x/blob/master/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C#L147

   if (validComponents[cmpt] == -1) continue;

All directions were skipped, without any warning. I went straight for "checkMesh" and got the answer: all directions were compromised. The reason: a single unassigned face was set to "empty", even after using "subsetMesh" + "createPatch" + "changeDictionary"; the face didn't have its normal aligned with any direction.

To reproduce, use in the tutorial "incompressible/simpleFoam/pitzDaily" the following commands:

  blockMesh
  transformPoints -rotate '( (1 0 0) (1 1 1) )'
  simpleFoam

and no momentum direction equations shall be solved! Without any warning at all...
Additional InformationStill need to check if the multi-region solvers also need this check added to them...
TagsNo tags attached.

Activities

wyldckat

2015-10-16 23:34

updater  

createMesh.H (516 bytes)   
//
// createMesh.H
// ~~~~~~~~~~~~

    Foam::Info
        << "Create mesh for time = "
        << runTime.timeName() << Foam::nl << Foam::endl;

    Foam::fvMesh mesh
    (
        Foam::IOobject
        (
            Foam::fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            Foam::IOobject::MUST_READ
        )
    );

    if(mesh.nSolutionD()<2)
    {
        Foam::Warning
            << "The loaded mesh has less than 2 valid solution directions."
            << endl;
    }
createMesh.H (516 bytes)   

wyldckat

2015-10-16 23:34

updater  

createNamedMesh.H (872 bytes)   
//
// createNamedMesh.H
// ~~~~~~~~~~~~~~~~~

    Foam::word regionName;

    if (args.optionReadIfPresent("region", regionName))
    {
        Foam::Info
            << "Create mesh " << regionName << " for time = "
            << runTime.timeName() << Foam::nl << Foam::endl;
    }
    else
    {
        regionName = Foam::fvMesh::defaultRegion;
        Foam::Info
            << "Create mesh for time = "
            << runTime.timeName() << Foam::nl << Foam::endl;
    }

    Foam::fvMesh mesh
    (
        Foam::IOobject
        (
            regionName,
            runTime.timeName(),
            runTime,
            Foam::IOobject::MUST_READ
        )
    );

    if(mesh.nSolutionD()<2)
    {
        Foam::Warning
            << "The loaded mesh region '" << regionName
            << "' has less than 2 valid solution directions."
            << endl;
    }
createNamedMesh.H (872 bytes)   

wyldckat

2015-10-16 23:41

updater  

createFluidMeshes.H (830 bytes)   
    const wordList fluidNames(rp["fluid"]);

    PtrList<fvMesh> fluidRegions(fluidNames.size());

    forAll(fluidNames, i)
    {
        Info<< "Create fluid mesh for region " << fluidNames[i]
            << " for time = " << runTime.timeName() << nl << endl;

        fluidRegions.set
        (
            i,
            new fvMesh
            (
                IOobject
                (
                    fluidNames[i],
                    runTime.timeName(),
                    runTime,
                    IOobject::MUST_READ
                )
            )
        );

        if(fluidRegions[i].nSolutionD()<2)
        {
            Foam::Warning
                << "The loaded mesh region '" << fluidNames[i]
                << "' has less than 2 valid solution directions."
                << endl;
        }
    }
createFluidMeshes.H (830 bytes)   

wyldckat

2015-10-16 23:42

updater  

createSolidMeshes.H (1,055 bytes)   
    const wordList solidsNames(rp["solid"]);

    PtrList<fvMesh> solidRegions(solidsNames.size());

    forAll(solidsNames, i)
    {
        Info<< "Create solid mesh for region " << solidsNames[i]
            << " for time = " << runTime.timeName() << nl << endl;

        solidRegions.set
        (
            i,
            new fvMesh
            (
                IOobject
                (
                    solidsNames[i],
                    runTime.timeName(),
                    runTime,
                    IOobject::MUST_READ
                )
            )
        );

        if(solidRegions[i].nSolutionD()<2)
        {
            Foam::Warning
                << "The loaded mesh region '" << solidsNames[i]
                << "' has less than 2 valid solution directions."
                << endl;
        }

        // Force calculation of geometric properties to prevent it being done
        // later in e.g. some boundary evaluation
        //(void)solidRegions[i].weights();
        //(void)solidRegions[i].deltaCoeffs();
    }
createSolidMeshes.H (1,055 bytes)   

wyldckat

2015-10-16 23:46

updater   ~0005411

The attached file "createFluidMeshes.H" is for the folders "applications/solvers/heatTransfer/chtMultiRegionFoam/fluid" and "applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid".

The attached file "createSolidMeshes.H" is for the folder "applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid".

henry

2015-10-17 11:33

manager   ~0005412

Solving 1D cases is valid in OpenFOAM so I don't think that it should be a warning. We could provide the number of solution directions as information but given that the user would generally already know this it might be considered clutter.

In the circumstances you describe I would use checkMesh to check that the BCs are appropriate and the number of solution directions are as expected. What does checkMesh print for your case? If it is not sufficient to diagnose the problem then we should add more checks and output to checkMesh.

wyldckat

2015-10-17 12:05

updater   ~0005413

On the rotated "incompressible/simpleFoam/pitzDaily" case I mentioned about, checkMesh with the default options only reports this:

  ***Number of edges not aligned with or perpendicular to non-empty directions: 61966
   <<Writing 25012 points on non-aligned edges to set nonAlignedEdges

This does hint at the problem.

But it is a bit daunting when the solver doesn't solve any of the momentum directions, when at first sight everything seems fine, given that it can solve everything else.
The "best part" is that all other equations are solved either way, so this can also give a false sense of security, given that the simulation still runs and can in some cases "converge". It has continuity errors in this test with the rotated pitzDaily, but since I didn't change the orientation of U at the inlet, it's natural that it doesn't behave all that well.

Either way, if a warning can be given, it does make a lot more sense that it only warns about when no directions will be solved, instead of detecting if it's lesser than 2.

Placing the warning in the "createMesh" sections that are mostly not related to mesh manipulators seemed to me to be the less intrusive way, although it does lead to a bit more repeated code than I would like. But placing it in the constructor didn't look like it would be the right place.

henry

2015-10-17 12:20

manager   ~0005414

The momentum equation is not solved because the case is appaently 0D which is also a valid option in OpenFOAM. All the other equations are solved for oD cases. I don't see a need to warn that a case is 0D except in cases for which it should not be 0D.

So the only issue here is that your mesh has errors in it which need to fixed to run. The purpose of checkMesh is to check that the mesh is suitable to run on. The error message you get from checkMesh is serious and should be dealt with before proceeding to run on the case. If you would like more diagnostics from checkMesh concening the number of valid geometric directions and solution directions I would be happy for this to be added to checkMesh.

wyldckat

2015-10-17 12:40

updater  

checkGeometry.C (24,791 bytes)   
#include "checkGeometry.H"
#include "polyMesh.H"
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
#include "EdgeMap.H"
#include "wedgePolyPatch.H"
#include "unitConversion.H"
#include "polyMeshTetDecomposition.H"


// Find wedge with opposite orientation. Note: does not actually check that
// it is opposite, only that it has opposite normal and same axis
Foam::label Foam::findOppositeWedge
(
    const polyMesh& mesh,
    const wedgePolyPatch& wpp
)
{
    const polyBoundaryMesh& patches = mesh.boundaryMesh();

    scalar wppCosAngle = wpp.cosAngle();

    forAll(patches, patchI)
    {
        if
        (
            patchI != wpp.index()
         && patches[patchI].size()
         && isA<wedgePolyPatch>(patches[patchI])
        )
        {
            const wedgePolyPatch& pp =
                refCast<const wedgePolyPatch>(patches[patchI]);

            // Calculate (cos of) angle to wpp (not pp!) centre normal
            scalar ppCosAngle = wpp.centreNormal() & pp.n();

            if
            (
                pp.size() == wpp.size()
             && mag(pp.axis() & wpp.axis()) >= (1-1e-3)
             && mag(ppCosAngle - wppCosAngle) >= 1e-3
            )
            {
                return patchI;
            }
        }
    }
    return -1;
}


bool Foam::checkWedges
(
    const polyMesh& mesh,
    const bool report,
    const Vector<label>& directions,
    labelHashSet* setPtr
)
{
    // To mark edges without calculating edge addressing
    EdgeMap<label> edgesInError;

    const pointField& p = mesh.points();
    const faceList& fcs = mesh.faces();


    const polyBoundaryMesh& patches = mesh.boundaryMesh();
    forAll(patches, patchI)
    {
        if (patches[patchI].size() && isA<wedgePolyPatch>(patches[patchI]))
        {
            const wedgePolyPatch& pp =
                refCast<const wedgePolyPatch>(patches[patchI]);

            scalar wedgeAngle = acos(pp.cosAngle());

            if (report)
            {
                Info<< "    Wedge " << pp.name() << " with angle "
                    << radToDeg(wedgeAngle) << " degrees"
                    << endl;
            }

            // Find opposite
            label oppositePatchI = findOppositeWedge(mesh, pp);

            if (oppositePatchI == -1)
            {
                if (report)
                {
                    Info<< " ***Cannot find opposite wedge for wedge "
                        << pp.name() << endl;
                }
                return true;
            }

            const wedgePolyPatch& opp =
                refCast<const wedgePolyPatch>(patches[oppositePatchI]);


            if (mag(opp.axis() & pp.axis()) < (1-1e-3))
            {
                if (report)
                {
                    Info<< " ***Wedges do not have the same axis."
                        << " Encountered " << pp.axis()
                        << " on patch " << pp.name()
                        << " which differs from " << opp.axis()
                        << " on opposite wedge patch" << opp.axis()
                        << endl;
                }
                return true;
            }



            // Mark edges on wedgePatches
            forAll(pp, i)
            {
                const face& f = pp[i];
                forAll(f, fp)
                {
                    label p0 = f[fp];
                    label p1 = f.nextLabel(fp);
                    edgesInError.insert(edge(p0, p1), -1);  // non-error value
                }
            }


            // Check that wedge patch is flat
            const point& p0 = p[pp.meshPoints()[0]];
            forAll(pp.meshPoints(), i)
            {
                const point& pt = p[pp.meshPoints()[i]];
                scalar d = mag((pt - p0) & pp.n());

                if (d > ROOTSMALL)
                {
                    if (report)
                    {
                        Info<< " ***Wedge patch " << pp.name() << " not planar."
                            << " Point " << pt << " is not in patch plane by "
                            << d << " metre."
                            << endl;
                    }
                    return true;
                }
            }
        }
    }



    // Check all non-wedge faces
    label nEdgesInError = 0;

    forAll(fcs, faceI)
    {
        const face& f = fcs[faceI];

        forAll(f, fp)
        {
            label p0 = f[fp];
            label p1 = f.nextLabel(fp);
            if (p0 < p1)
            {
                vector d(p[p1]-p[p0]);
                scalar magD = mag(d);

                if (magD > ROOTVSMALL)
                {
                    d /= magD;

                    // Check how many empty directions are used by the edge.
                    label nEmptyDirs = 0;
                    label nNonEmptyDirs = 0;
                    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
                    {
                        if (mag(d[cmpt]) > 1e-6)
                        {
                            if (directions[cmpt] == 0)
                            {
                                nEmptyDirs++;
                            }
                            else
                            {
                                nNonEmptyDirs++;
                            }
                        }
                    }

                    if (nEmptyDirs == 0)
                    {
                        // Purely in ok directions.
                    }
                    else if (nEmptyDirs == 1)
                    {
                        // Ok if purely in empty directions.
                        if (nNonEmptyDirs > 0)
                        {
                            if (edgesInError.insert(edge(p0, p1), faceI))
                            {
                                nEdgesInError++;
                            }
                        }
                    }
                    else if (nEmptyDirs > 1)
                    {
                        // Always an error
                        if (edgesInError.insert(edge(p0, p1), faceI))
                        {
                            nEdgesInError++;
                        }
                    }
                }
            }
        }
    }

    label nErrorEdges = returnReduce(nEdgesInError, sumOp<label>());

    if (nErrorEdges > 0)
    {
        if (report)
        {
            Info<< " ***Number of edges not aligned with or perpendicular to "
                << "non-empty directions: " << nErrorEdges << endl;
        }

        if (setPtr)
        {
            setPtr->resize(2*nEdgesInError);
            forAllConstIter(EdgeMap<label>, edgesInError, iter)
            {
                if (iter() >= 0)
                {
                    setPtr->insert(iter.key()[0]);
                    setPtr->insert(iter.key()[1]);
                }
            }
        }

        return true;
    }
    else
    {
        if (report)
        {
            Info<< "    All edges aligned with or perpendicular to "
                << "non-empty directions." << endl;
        }
        return false;
    }
}


namespace Foam
{
    //- Default transformation behaviour for position
    class transformPositionList
    {
    public:

        //- Transform patch-based field
        void operator()
        (
            const coupledPolyPatch& cpp,
            List<pointField>& pts
        ) const
        {
            // Each element of pts is all the points in the face. Convert into
            // lists of size cpp to transform.

            List<pointField> newPts(pts.size());
            forAll(pts, faceI)
            {
                newPts[faceI].setSize(pts[faceI].size());
            }

            label index = 0;
            while (true)
            {
                label n = 0;

                // Extract for every face the i'th position
                pointField ptsAtIndex(pts.size(), vector::zero);
                forAll(cpp, faceI)
                {
                    const pointField& facePts = pts[faceI];
                    if (facePts.size() > index)
                    {
                        ptsAtIndex[faceI] = facePts[index];
                        n++;
                    }
                }

                if (n == 0)
                {
                    break;
                }

                // Now ptsAtIndex will have for every face either zero or
                // the position of the i'th vertex. Transform.
                cpp.transformPosition(ptsAtIndex);

                // Extract back from ptsAtIndex into newPts
                forAll(cpp, faceI)
                {
                    pointField& facePts = newPts[faceI];
                    if (facePts.size() > index)
                    {
                        facePts[index] = ptsAtIndex[faceI];
                    }
                }

                index++;
            }

            pts.transfer(newPts);
        }
    };
}


bool Foam::checkCoupledPoints
(
    const polyMesh& mesh,
    const bool report,
    labelHashSet* setPtr
)
{
    const pointField& p = mesh.points();
    const faceList& fcs = mesh.faces();
    const polyBoundaryMesh& patches = mesh.boundaryMesh();

    // Zero'th point on coupled faces
    //pointField nbrZeroPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
    List<pointField> nbrPoints(fcs.size() - mesh.nInternalFaces());

    // Exchange zero point
    forAll(patches, patchI)
    {
        if (patches[patchI].coupled())
        {
            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
            (
                patches[patchI]
            );

            forAll(cpp, i)
            {
                label bFaceI = cpp.start() + i - mesh.nInternalFaces();
                const face& f = cpp[i];
                nbrPoints[bFaceI].setSize(f.size());
                forAll(f, fp)
                {
                    const point& p0 = p[f[fp]];
                    nbrPoints[bFaceI][fp] = p0;
                }
            }
        }
    }
    syncTools::syncBoundaryFaceList
    (
        mesh,
        nbrPoints,
        eqOp<pointField>(),
        transformPositionList()
    );

    // Compare to local ones. Use same tolerance as for matching
    label nErrorFaces = 0;
    scalar avgMismatch = 0;
    label nCoupledPoints = 0;

    forAll(patches, patchI)
    {
        if (patches[patchI].coupled())
        {
            const coupledPolyPatch& cpp =
                refCast<const coupledPolyPatch>(patches[patchI]);

            if (cpp.owner())
            {
                scalarField smallDist
                (
                    cpp.calcFaceTol
                    (
                        //cpp.matchTolerance(),
                        cpp,
                        cpp.points(),
                        cpp.faceCentres()
                    )
                );

                forAll(cpp, i)
                {
                    label bFaceI = cpp.start() + i - mesh.nInternalFaces();
                    const face& f = cpp[i];

                    if (f.size() != nbrPoints[bFaceI].size())
                    {
                        FatalErrorIn
                        (
                            "Foam::checkCoupledPoints\n"
                            "(\n"
                            "   const polyMesh&, const bool, labelHashSet*\n"
                            ")\n"
                        )   << "Local face size : " << f.size()
                            << " does not equal neighbour face size : "
                            << nbrPoints[bFaceI].size()
                            << abort(FatalError);
                    }

                    label fp = 0;
                    forAll(f, j)
                    {
                        const point& p0 = p[f[fp]];
                        scalar d = mag(p0 - nbrPoints[bFaceI][j]);

                        if (d > smallDist[i])
                        {
                            if (setPtr)
                            {
                                setPtr->insert(cpp.start()+i);
                            }
                            nErrorFaces++;

                            break;
                        }

                        avgMismatch += d;
                        nCoupledPoints++;

                        fp = f.rcIndex(fp);
                    }
                }
            }
        }
    }

    reduce(nErrorFaces, sumOp<label>());
    reduce(avgMismatch, maxOp<scalar>());
    reduce(nCoupledPoints, sumOp<label>());

    if (nCoupledPoints > 0)
    {
        avgMismatch /= nCoupledPoints;
    }

    if (nErrorFaces > 0)
    {
        if (report)
        {
            Info<< "  **Error in coupled point location: "
                << nErrorFaces
                << " faces have their 0th or consecutive vertex not opposite"
                << " their coupled equivalent. Average mismatch "
                << avgMismatch << "."
                << endl;
        }

        return true;
    }
    else
    {
        if (report)
        {
            Info<< "    Coupled point location match (average "
                << avgMismatch << ") OK." << endl;
        }

        return false;
    }
}


Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
{
    label noFailedChecks = 0;

    Info<< "\nChecking geometry..." << endl;

    // Get a small relative length from the bounding box
    const boundBox& globalBb = mesh.bounds();

    Info<< "    Overall domain bounding box "
        << globalBb.min() << " " << globalBb.max() << endl;


    // Min length
    scalar minDistSqr = magSqr(1e-6 * globalBb.span());

    // Non-empty directions
    const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one)/2;
    Info<< "    Mesh (non-empty, non-wedge) directions " << validDirs << endl;

    const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one)/2;
    Info<< "    Mesh (non-empty) directions " << solDirs
        << ", " << mesh.nSolutionD() << "D solution space." << endl;

    if (mesh.nGeometricD() < 3)
    {
        pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);

        if
        (
            (
                validDirs != solDirs
             && checkWedges(mesh, true, validDirs, &nonAlignedPoints)
            )
         || (
                validDirs == solDirs
             && mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints)
            )
        )
        {
            noFailedChecks++;
            label nNonAligned = returnReduce
            (
                nonAlignedPoints.size(),
                sumOp<label>()
            );

            if (nNonAligned > 0)
            {
                Info<< "  <<Writing " << nNonAligned
                    << " points on non-aligned edges to set "
                    << nonAlignedPoints.name() << endl;
                nonAlignedPoints.instance() = mesh.pointsInstance();
                nonAlignedPoints.write();
            }
        }
    }

    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,
                mesh.geometricD()
            )
        )
        {
            noFailedChecks++;

            label nNonClosed = returnReduce(cells.size(), sumOp<label>());

            if (nNonClosed > 0)
            {
                Info<< "  <<Writing " << nNonClosed
                    << " non closed cells to set " << cells.name() << endl;
                cells.instance() = mesh.pointsInstance();
                cells.write();
            }
        }

        label nHighAspect = returnReduce(aspectCells.size(), sumOp<label>());

        if (nHighAspect > 0)
        {
            Info<< "  <<Writing " << nHighAspect
                << " cells with high aspect ratio to set "
                << aspectCells.name() << endl;
            aspectCells.instance() = mesh.pointsInstance();
            aspectCells.write();
        }
    }

    {
        faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces()/100+1);
        if (mesh.checkFaceAreas(true, &faces))
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info<< "  <<Writing " << nFaces
                    << " zero area faces to set " << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
            }
        }
    }

    {
        cellSet cells(mesh, "zeroVolumeCells", mesh.nCells()/100+1);
        if (mesh.checkCellVolumes(true, &cells))
        {
            noFailedChecks++;

            label nCells = returnReduce(cells.size(), sumOp<label>());

            if (nCells > 0)
            {
                Info<< "  <<Writing " << nCells
                    << " zero volume cells to set " << cells.name() << endl;
                cells.instance() = mesh.pointsInstance();
                cells.write();
            }
        }
    }

    {
        faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100+1);
        if (mesh.checkFaceOrthogonality(true, &faces))
        {
            noFailedChecks++;
        }

        label nFaces = returnReduce(faces.size(), sumOp<label>());

        if (nFaces > 0)
        {
            Info<< "  <<Writing " << nFaces
                << " non-orthogonal faces to set " << faces.name() << endl;
            faces.instance() = mesh.pointsInstance();
            faces.write();
        }
    }

    {
        faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces()/100 + 1);
        if (mesh.checkFacePyramids(true, -SMALL, &faces))
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info<< "  <<Writing " << nFaces
                    << " faces with incorrect orientation to set "
                    << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
            }
        }
    }

    {
        faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
        if (mesh.checkFaceSkewness(true, &faces))
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info<< "  <<Writing " << nFaces
                    << " skew faces to set " << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
            }
        }
    }

    {
        faceSet faces(mesh, "coupledFaces", mesh.nFaces()/100 + 1);
        if (checkCoupledPoints(mesh, true, &faces))
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info<< "  <<Writing " << nFaces
                    << " faces with incorrectly matched 0th (or consecutive)"
                    << " vertex to set "
                    << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
            }
        }
    }

    if (allGeometry)
    {
        faceSet faces(mesh, "lowQualityTetFaces", mesh.nFaces()/100+1);
        if
        (
            polyMeshTetDecomposition::checkFaceTets
            (
                mesh,
                polyMeshTetDecomposition::minTetQuality,
                true,
                &faces
            )
        )
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info<< "  <<Writing " << nFaces
                    << " faces with low quality or negative volume "
                    << "decomposition tets to set " << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
            }
        }
    }

    if (allGeometry)
    {
        // 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++;

            label nPoints = returnReduce(points.size(), sumOp<label>());

            if (nPoints > 0)
            {
                Info<< "  <<Writing " << nPoints
                    << " points on short edges to set " << points.name()
                    << endl;
                points.instance() = mesh.pointsInstance();
                points.write();
            }
        }

        label nEdgeClose = returnReduce(points.size(), sumOp<label>());

        if (mesh.checkPointNearness(false, minDistSqr, &points))
        {
            //noFailedChecks++;

            label nPoints = returnReduce(points.size(), sumOp<label>());

            if (nPoints > nEdgeClose)
            {
                pointSet nearPoints(mesh, "nearPoints", points);
                Info<< "  <<Writing " << nPoints
                    << " near (closer than " << Foam::sqrt(minDistSqr)
                    << " apart) points to set " << nearPoints.name() << endl;
                nearPoints.instance() = mesh.pointsInstance();
                nearPoints.write();
            }
        }
    }

    if (allGeometry)
    {
        faceSet faces(mesh, "concaveFaces", mesh.nFaces()/100 + 1);
        if (mesh.checkFaceAngles(true, 10, &faces))
        {
            //noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info<< "  <<Writing " << nFaces
                    << " faces with concave angles to set " << faces.name()
                    << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
            }
        }
    }

    if (allGeometry)
    {
        faceSet faces(mesh, "warpedFaces", mesh.nFaces()/100 + 1);
        if (mesh.checkFaceFlatness(true, 0.8, &faces))
        {
            //noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info<< "  <<Writing " << nFaces
                    << " warped faces to set " << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
            }
        }
    }

    if (allGeometry)
    {
        cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
        if (mesh.checkCellDeterminant(true, &cells))
        {
            noFailedChecks++;

            label nCells = returnReduce(cells.size(), sumOp<label>());

            Info<< "  <<Writing " << nCells
                << " under-determined cells to set " << cells.name() << endl;
            cells.instance() = mesh.pointsInstance();
            cells.write();
        }
    }

    if (allGeometry)
    {
        cellSet cells(mesh, "concaveCells", mesh.nCells()/100);
        if (mesh.checkConcaveCells(true, &cells))
        {
            noFailedChecks++;

            label nCells = returnReduce(cells.size(), sumOp<label>());

            Info<< "  <<Writing " << nCells
                << " concave cells to set " << cells.name() << endl;
            cells.instance() = mesh.pointsInstance();
            cells.write();
        }
    }

    if (allGeometry)
    {
        faceSet faces(mesh, "lowWeightFaces", mesh.nFaces()/100);
        if (mesh.checkFaceWeight(true, 0.05, &faces))
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            Info<< "  <<Writing " << nFaces
                << " faces with low interpolation weights to set "
                << faces.name() << endl;
            faces.instance() = mesh.pointsInstance();
            faces.write();
        }
    }

    if (allGeometry)
    {
        faceSet faces(mesh, "lowVolRatioFaces", mesh.nFaces()/100);
        if (mesh.checkVolRatio(true, 0.01, &faces))
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            Info<< "  <<Writing " << nFaces
                << " faces with low volume ratio cells to set "
                << faces.name() << endl;
            faces.instance() = mesh.pointsInstance();
            faces.write();
        }
    }

    return noFailedChecks;
}
checkGeometry.C (24,791 bytes)   

wyldckat

2015-10-17 12:44

updater   ~0005416

Attached the file "checkGeometry.C" for the folder "applications/utilities/mesh/manipulation/checkMesh/".

The proposed change will affect this line in the output of checkMesh:

    Mesh (non-empty) directions (1 1 0)

Resulting in these kinds of output:

    Mesh (non-empty) directions (1 1 0), 2D solution space.

    Mesh (non-empty) directions (0 0 0), 0D solution space.


Using the expression "solution space" sounds weird, but it seems to be the one that conveys a good indication of what the directions/dimensions stand for, given that the 4th dimension usually is "time"...

henry

2015-10-17 13:16

manager   ~0005417

The number of geometric directions and solutions directions may be different, e.g. in the case of a wedge geometry. Does "solution space" convey geometric directions or solutions directions?

henry

2015-10-17 13:59

manager   ~0005418

I have added the number of solution and geometric directions to the output from checkGeometry.C in a manner consistent with the naming of the corresponding functions in OpenFOAM. I hope this sufficient.

OpenFOAM-dev: commit a10397d10fd68411fa008e616b0fb95b2fba5ecf

Issue History

Date Modified Username Field Change
2015-10-16 23:16 wyldckat New Issue
2015-10-16 23:17 wyldckat Steps to Reproduce Updated
2015-10-16 23:18 wyldckat Steps to Reproduce Updated
2015-10-16 23:34 wyldckat File Added: createMesh.H
2015-10-16 23:34 wyldckat File Added: createNamedMesh.H
2015-10-16 23:41 wyldckat File Added: createFluidMeshes.H
2015-10-16 23:42 wyldckat File Added: createSolidMeshes.H
2015-10-16 23:46 wyldckat Note Added: 0005411
2015-10-17 11:33 henry Note Added: 0005412
2015-10-17 12:05 wyldckat Note Added: 0005413
2015-10-17 12:20 henry Note Added: 0005414
2015-10-17 12:40 wyldckat File Added: checkGeometry.C
2015-10-17 12:44 wyldckat Note Added: 0005416
2015-10-17 13:16 henry Note Added: 0005417
2015-10-17 13:59 henry Note Added: 0005418
2015-10-17 13:59 henry Status new => resolved
2015-10-17 13:59 henry Resolution open => fixed
2015-10-17 13:59 henry Assigned To => henry