automatic.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2012-2017 OpenFOAM Foundation
    Copyright (C) 2016-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "automatic.H"
#include "addToRunTimeSelectionTable.H"
#include "triSurfaceMesh.H"
#include "foamVtkSurfaceWriter.H"
#include "primitivePatchInterpolation.H"
#include "Time.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    defineTypeNameAndDebug(automatic, 0);
    addToRunTimeSelectionTable(cellSizeCalculationType, automatic, dictionary);
}
 
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
void Foam::automatic::smoothField(triSurfaceScalarField& surf)
{
    label nSmoothingIterations = 10;
 
    for (label iter = 0; iter < nSmoothingIterations; ++iter)
    {
        const pointField& faceCentres = surface_.faceCentres();
 
        forAll(surf, sI)
        {
            const labelList& faceFaces = surface_.faceFaces()[sI];
 
            const point& fC = faceCentres[sI];
            const scalar value = surf[sI];
 
            scalar newValue = 0;
            scalar totalDist = 0;
 
            label nFaces = 0;
 
            forAll(faceFaces, fI)
            {
                const label faceLabel = faceFaces[fI];
                const point& faceCentre = faceCentres[faceLabel];
 
                const scalar faceValue = surf[faceLabel];
                const scalar distance = mag(faceCentre - fC);
 
                newValue += faceValue/(distance + SMALL);
 
                totalDist += 1.0/(distance + SMALL);
 
                if (value < faceValue)
                {
                    nFaces++;
                }
            }
 
            // Do not smooth out the peak values
            if (nFaces == faceFaces.size())
            {
                continue;
            }
 
            surf[sI] = newValue/totalDist;
        }
    }
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::automatic::automatic
(
    const dictionary& cellSizeCalcTypeDict,
    const triSurfaceMesh& surface,
    const scalar defaultCellSize
)
:
    cellSizeCalculationType
    (
        typeName,
        cellSizeCalcTypeDict,
        surface,
        defaultCellSize
    ),
    coeffsDict_(cellSizeCalcTypeDict.optionalSubDict(typeName + "Coeffs")),
    surfaceName_(surface.searchableSurface::name()),
 
    readCurvature_(coeffsDict_.get<bool>("curvature")),
    readFeatureProximity_(coeffsDict_.get<bool>("featureProximity")),
    readInternalCloseness_(coeffsDict_.get<bool>("internalCloseness")),
 
    curvatureFile_(coeffsDict_.get<word>("curvatureFile")),
    featureProximityFile_(coeffsDict_.get<word>("featureProximityFile")),
    internalClosenessFile_(coeffsDict_.get<word>("internalClosenessFile")),
 
    curvatureCellSizeCoeff_
    (
        coeffsDict_.get<scalar>("curvatureCellSizeCoeff")
    ),
    maximumCellSize_
    (
        coeffsDict_.get<scalar>("maximumCellSizeCoeff") * defaultCellSize
    )
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
Foam::tmp<Foam::triSurfacePointScalarField> Foam::automatic::load()
{
    Info<< indent
        << "Calculating cell size on surface: " << surfaceName_ << endl;
 
    tmp<triSurfacePointScalarField> tPointCellSize
    (
        new triSurfacePointScalarField
        (
            IOobject
            (
                surfaceName_ + ".cellSize",
                surface_.searchableSurface::time().constant(),
                "triSurface",
                surface_.searchableSurface::time(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            surface_,
            dimLength,
            scalarField(surface_.nPoints(), maximumCellSize_)
        )
    );
 
    triSurfacePointScalarField& pointCellSize = tPointCellSize.ref();
 
    if (readCurvature_)
    {
        Info<< indent
            << "Reading curvature         : " << curvatureFile_ << endl;
 
        triSurfacePointScalarField curvature
        (
            IOobject
            (
                curvatureFile_,
                surface_.searchableSurface::time().constant(),
                "triSurface",
                surface_.searchableSurface::time(),
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            ),
            surface_,
            dimLength,
            true
        );
 
        forAll(pointCellSize, pI)
        {
            pointCellSize[pI] =
                min
                (
                    1.0
                   /max
                    (
                        (1.0/curvatureCellSizeCoeff_)*mag(curvature[pI]),
                        1.0/maximumCellSize_
                    ),
                    pointCellSize[pI]
                );
        }
    }
 
    PrimitivePatchInterpolation
    <
        PrimitivePatch<::Foam::List<labelledTri>, pointField>
    > patchInterpolate(surface_);
 
    const Map<label>& meshPointMap = surface_.meshPointMap();
 
    if (readInternalCloseness_)
    {
        Info<< indent
            << "Reading internal closeness: " << internalClosenessFile_ << endl;
 
        triSurfaceScalarField internalCloseness
        (
            IOobject
            (
                internalClosenessFile_,
                surface_.searchableSurface::time().constant(),
                "triSurface",
                surface_.searchableSurface::time(),
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            ),
            surface_,
            dimLength,
            true
        );
 
        scalarField internalClosenessPointField
        (
            patchInterpolate.faceToPointInterpolate(internalCloseness)
        );
 
        forAll(pointCellSize, pI)
        {
            pointCellSize[pI] =
                min
                (
                    internalClosenessPointField[meshPointMap[pI]],
                    pointCellSize[pI]
                );
        }
    }
 
    if (readFeatureProximity_)
    {
        Info<< indent
            << "Reading feature proximity : " << featureProximityFile_ << endl;
 
        triSurfaceScalarField featureProximity
        (
            IOobject
            (
                featureProximityFile_,
                surface_.searchableSurface::time().constant(),
                "triSurface",
                surface_.searchableSurface::time(),
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            ),
            surface_,
            dimLength,
            true
        );
 
        scalarField featureProximityPointField
        (
            patchInterpolate.faceToPointInterpolate(featureProximity)
        );
 
        forAll(pointCellSize, pI)
        {
            pointCellSize[pI] =
                min
                (
                    featureProximityPointField[meshPointMap[pI]],
                    pointCellSize[pI]
                );
        }
    }
 
    //smoothField(surfaceCellSize);
 
    pointCellSize.write();
 
    if (debug)
    {
        faceList faces(surface_.size());
 
        forAll(surface_, fI)
        {
            faces[fI] = surface_.triSurface::operator[](fI);
        }
 
        vtk::surfaceWriter vtkWriter
        (
            surface_.points(),
            faces,
            (
                surface_.searchableSurface::time().constant()
              / "triSurface"
              / surfaceName_.nameLessExt() + "_cellSize"
            )
        );
 
 
        vtkWriter.writeGeometry();
 
        vtkWriter.beginPointData(1);
 
        vtkWriter.write("cellSize", pointCellSize);
    }
 
    return tPointCellSize;
}
 
 
// ************************************************************************* //