Skip to content
Snippets Groups Projects
writeFields.C 17.4 KiB
Newer Older
  • Learn to ignore specific revisions
  • #include "writeFields.H"
    #include "volFields.H"
    
    #include "surfaceFields.H"
    
    #include "polyMeshTools.H"
    #include "zeroGradientFvPatchFields.H"
    #include "syncTools.H"
    
    #include "tetPointRef.H"
    
    #include "regionSplit.H"
    
    #include "wallDist.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 minFaceToCell
    (
        const surfaceScalarField& faceData,
        volScalarField& cellData,
        const bool correctBoundaryConditions
    )
    {
        scalarField& cellFld = cellData.ref();
    
        cellFld = GREAT;
    
        const labelUList& own = cellData.mesh().owner();
        const labelUList& nei = cellData.mesh().neighbour();
    
        // Internal faces
        forAll(own, facei)
        {
            cellFld[own[facei]] = min(cellFld[own[facei]], faceData[facei]);
            cellFld[nei[facei]] = min(cellFld[nei[facei]], faceData[facei]);
        }
    
        // Patch faces
        forAll(faceData.boundaryField(), patchi)
        {
            const fvsPatchScalarField& fvp = faceData.boundaryField()[patchi];
            const labelUList& fc = fvp.patch().faceCells();
    
            forAll(fc, i)
            {
                cellFld[fc[i]] = min(cellFld[fc[i]], fvp[i]);
            }
        }
    
        volScalarField::Boundary& bfld = cellData.boundaryFieldRef();
    
        forAll(bfld, patchi)
        {
            bfld[patchi] = faceData.boundaryField()[patchi];
        }
        if (correctBoundaryConditions)
        {
            cellData.correctBoundaryConditions();
        }
    }
    
    
    
    void Foam::writeFields
    (
        const fvMesh& mesh,
    
        const wordHashSet& 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(scalar(1), faceOrthogonality))
    
                )
            );
    
            //- Cell field - max of either face
            volScalarField cellNonOrthoAngle
            (
                IOobject
                (
                    "nonOrthoAngle",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
    
                    IOobject::AUTO_WRITE,
                    false
    
                dimensionedScalar(dimless, Zero),
    
                calculatedFvPatchScalarField::typeName
            );
            //- Take max
            maxFaceToCell(nonOrthoAngle, cellNonOrthoAngle);
            Info<< "    Writing non-orthogonality (angle) to "
                << cellNonOrthoAngle.name() << endl;
            cellNonOrthoAngle.write();
        }
    
        if (selectedFields.found("faceWeight"))
        {
            volScalarField cellWeights
            (
                IOobject
                (
                    "faceWeight",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
    
                    IOobject::AUTO_WRITE,
                    false
    
                dimensionedScalar(dimless, Zero),
    
                wordList                                // wanted bc types
                (
                    mesh.boundary().size(),
                    calculatedFvPatchScalarField::typeName
                ),
                mesh.weights().boundaryField().types()  // current bc types
    
            minFaceToCell(mesh.weights(), cellWeights, false);
    
            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,
                    false
    
                dimensionedScalar(dimless, Zero),
    
                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(dimless, Zero),
    
                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(dimless, Zero),
    
                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(dimless, Zero),
    
                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(dimVolume, Zero),
    
                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,
                    false
    
                dimensionedScalar(dimless, Zero),
    
                calculatedFvPatchScalarField::typeName
            );
            //- Take min
            minFaceToCell(faceVolumeRatio, cellVolumeRatio);
            Info<< "    Writing cell volume ratio to "
                << cellVolumeRatio.name() << endl;
            cellVolumeRatio.write();
        }
    
    
        // minTetVolume
        if (selectedFields.found("minTetVolume"))
        {
            volScalarField minTetVolume
            (
                IOobject
                (
                    "minTetVolume",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE,
                    false
                ),
                mesh,
                dimensionedScalar("minTetVolume", dimless, GREAT),
                zeroGradientFvPatchScalarField::typeName
            );
    
    
            const labelList& own = mesh.faceOwner();
            const labelList& nei = mesh.faceNeighbour();
            const pointField& p = mesh.points();
            forAll(own, facei)
            {
                const face& f = mesh.faces()[facei];
                const point& fc = mesh.faceCentres()[facei];
    
                {
                    const point& ownCc = mesh.cellCentres()[own[facei]];
                    scalar& ownVol = minTetVolume[own[facei]];
                    forAll(f, fp)
                    {
                        scalar tetQual = tetPointRef
                        (
                            p[f[fp]],
                            p[f.nextLabel(fp)],
                            ownCc,
                            fc
                        ).quality();
                        ownVol = min(ownVol, tetQual);
                    }
                }
                if (mesh.isInternalFace(facei))
                {
                    const point& neiCc = mesh.cellCentres()[nei[facei]];
                    scalar& neiVol = minTetVolume[nei[facei]];
                    forAll(f, fp)
                    {
                        scalar tetQual = tetPointRef
                        (
                            p[f[fp]],
                            p[f.nextLabel(fp)],
                            fc,
                            neiCc
                        ).quality();
                        neiVol = min(neiVol, tetQual);
                    }
                }
            }
    
            minTetVolume.correctBoundaryConditions();
            Info<< "    Writing minTetVolume to " << minTetVolume.name() << endl;
            minTetVolume.write();
        }
    
    
        // minPyrVolume
        if (selectedFields.found("minPyrVolume"))
        {
            volScalarField minPyrVolume
            (
                IOobject
                (
                    "minPyrVolume",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE,
                    false
                ),
                mesh,
                dimensionedScalar("minPyrVolume", dimless, GREAT),
                zeroGradientFvPatchScalarField::typeName
            );
    
            // Get owner and neighbour pyr volumes
            scalarField ownPyrVol(mesh.nFaces());
            scalarField neiPyrVol(mesh.nInternalFaces());
            primitiveMeshTools::facePyramidVolume
            (
                mesh,
                mesh.points(),
                mesh.cellCentres(),
    
                ownPyrVol,
                neiPyrVol
            );
    
            // Get min pyr vol per cell
            scalarField& cellFld = minPyrVolume.ref();
            cellFld = GREAT;
    
            const labelUList& own = mesh.owner();
            const labelUList& nei = mesh.neighbour();
    
            // Internal faces
            forAll(own, facei)
            {
                cellFld[own[facei]] = min(cellFld[own[facei]], ownPyrVol[facei]);
                cellFld[nei[facei]] = min(cellFld[nei[facei]], neiPyrVol[facei]);
            }
    
            // Patch faces
            for (const auto& fvp : minPyrVolume.boundaryField())
            {
                const labelUList& fc = fvp.patch().faceCells();
    
                forAll(fc, i)
                {
                    const label meshFacei = fvp.patch().start();
                    cellFld[fc[i]] = min(cellFld[fc[i]], ownPyrVol[meshFacei]);
                }
            }
    
            minPyrVolume.correctBoundaryConditions();
            Info<< "    Writing minPyrVolume to " << minPyrVolume.name() << endl;
            minPyrVolume.write();
        }
    
    
        if (selectedFields.found("cellRegion"))
        {
            volScalarField cellRegion
            (
                IOobject
                (
                    "cellRegion",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
    
                    IOobject::AUTO_WRITE,
                    false
    
                dimensionedScalar(dimless, Zero),
    
                calculatedFvPatchScalarField::typeName
            );
    
            regionSplit rs(mesh);
            forAll(rs, celli)
            {
                cellRegion[celli] = rs[celli];
            }
            cellRegion.correctBoundaryConditions();
            Info<< "    Writing cell region to " << cellRegion.name() << endl;
            cellRegion.write();
        }
    
        if (selectedFields.found("wallDistance"))
        {
    
            // See if wallDist.method entry in fvSchemes before calling factory
            // method of wallDist. Have 'failing' version of wallDist::New instead?
            const dictionary& schemesDict =
                static_cast<const fvSchemes&>(mesh).schemesDict();
            if (schemesDict.found("wallDist"))
            {
                if (schemesDict.subDict("wallDist").found("method"))
                {
                    // Wall distance
                    volScalarField y("wallDistance", wallDist::New(mesh).y());
                    Info<< "    Writing wall distance to " << y.name() << endl;
                    y.write();
    
                    // Wall-reflection vectors
                    //const volVectorField& n = wallDist::New(mesh).n();
                    //Info<< "    Writing wall normal to " << n.name() << endl;
                    //n.write();
                }
            }
    
        if (selectedFields.found("cellZone"))
        {
            volScalarField cellZone
            (
                IOobject
                (
                    "cellZone",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE,
                    false
                ),
                mesh,
                dimensionedScalar(scalar(-1)),
                calculatedFvPatchScalarField::typeName
            );
    
            const cellZoneMesh& czs = mesh.cellZones();
            for (const auto& zone : czs)
            {
                UIndirectList<scalar>(cellZone, zone) = zone.index();
            }
    
            cellZone.correctBoundaryConditions();
            Info<< "    Writing cell zoning to " << cellZone.name() << endl;
            cellZone.write();
        }
        if (selectedFields.found("faceZone"))
        {
    
            // Determine for each face the zone index (scalar for ease of
            // manipulation)
            scalarField zoneID(mesh.nFaces(), -1);
            const faceZoneMesh& czs = mesh.faceZones();
            for (const auto& zone : czs)
            {
                UIndirectList<scalar>(zoneID, zone) = zone.index();
            }
    
    
            // Split into internal and boundary values
    
            surfaceScalarField faceZone
            (
                IOobject
                (
                    "faceZone",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE,
                    false
                ),
                mesh,
                dimensionedScalar(scalar(-1)),
                calculatedFvsPatchScalarField::typeName
            );
    
    
            faceZone.primitiveFieldRef() =
                SubField<scalar>(zoneID, mesh.nInternalFaces());
            surfaceScalarField::Boundary& bfld = faceZone.boundaryFieldRef();
            for (auto& pfld : bfld)
    
                const fvPatch& fvp = pfld.patch();
                pfld == SubField<scalar>(zoneID, fvp.size(), fvp.start());
    
            }
    
            //faceZone.correctBoundaryConditions();
            Info<< "    Writing face zoning to " << faceZone.name() << endl;
            faceZone.write();
        }