interpolationCellPointWallModified.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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 "interpolationCellPointWallModified.H"
#include "syncTools.H"
#include "volPointInterpolation.H"
#include "wallPolyPatch.H"
#include "zeroGradientFvPatchFields.H"
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
template<class Type>
template<class TYPE>
Foam::tmp<Foam::PointField<TYPE>>
Foam::interpolationCellPointWallModified<Type>::calcPointField
(
    const VolField<TYPE>& psi
) const
{
    FatalErrorInFunction
        << typeName << " interpolation is only defined for vector fields"
        << exit(FatalError);
 
    return tmp<PointField<TYPE>>(nullptr);
}
 
 
template<class Type>
Foam::tmp<Foam::pointVectorField>
Foam::interpolationCellPointWallModified<Type>::calcPointField
(
    const volVectorField& psi
) const
{
    using namespace constant::mathematical;
 
    const fvMesh& mesh = psi.mesh();
 
 
    // Create the point field
    tmp<pointVectorField> tPsip
    (
        pointVectorField::New
        (
            "volPointInterpolateWallModified(" + psi.name() + ')',
            pointMesh::New(mesh),
            dimensioned<Type>("zero", psi.dimensions(), Zero)
        )
    );
    pointVectorField& psip = tPsip.ref();
 
 
    // Interpolate to the points with wall patches extrapolated
    {
        wordList patchTypes(psi.boundaryField().size());
        forAll(patchTypes, patchi)
        {
            if (isA<wallPolyPatch>(mesh.boundaryMesh()[patchi]))
            {
                patchTypes[patchi] = zeroGradientFvPatchVectorField::typeName;
            }
            else
            {
                patchTypes[patchi] = calculatedFvPatchVectorField::typeName;
            }
        }
        volVectorField psiExtrapolated
        (
            IOobject
            (
                psi.name() + "Extrapolated",
                mesh.time().name(),
                mesh
            ),
            psi,
            patchTypes
        );
        psiExtrapolated.correctBoundaryConditions();
        volPointInterpolation::New(mesh).interpolate(psiExtrapolated, psip);
    }
 
 
    // Generate point normals across the entire boundary
    pointField pointNormals(mesh.nPoints(), vector::zero);
    {
        scalarField pointCount(mesh.nPoints(), 0);
 
        forAll(mesh.boundaryMesh(), patchi)
        {
            const polyPatch& patch = mesh.boundaryMesh()[patchi];
 
            forAll(patch, patchFacei)
            {
                const face& f = patch[patchFacei];
                const vector& n = patch.faceNormals()[patchFacei];
 
                forAll(f, i)
                {
                    pointNormals[f[i]] += n;
                    pointCount[f[i]] += 1;
                }
            }
        }
 
        syncTools::syncPointList
        (
            mesh,
            pointNormals,
            plusEqOp<vector>(),
            vector::zero
        );
 
        syncTools::syncPointList
        (
            mesh,
            pointCount,
            plusEqOp<scalar>(),
            scalar(0)
        );
 
        pointNormals /= max(pointCount, small);
    }
 
 
    // Calculate the rotation necessary from the vector to the (negative) point
    // normal needed to make the interpolated field point into the mesh
    scalarField theta0(mesh.nPoints(), -vGreat), theta1(mesh.nPoints(), vGreat);
    scalar maxVHatDotN = - vGreat;
    forAll(mesh.boundaryMesh(), patchi)
    {
        const polyPatch& patch = mesh.boundaryMesh()[patchi];
        if (isA<wallPolyPatch>(patch))
        {
            forAll(patch, patchFacei)
            {
                const label facei = patch.start() + patchFacei;
                const face& f = patch[patchFacei];
 
                for (label i = 1; i < f.size() - 1; ++ i)
                {
                    const triFace tri
                    (
                        tetIndices(mesh.faceOwner()[facei], facei, i)
                       .faceTriIs(mesh)
                    );
 
                    const vector n = tri.normal(mesh.points());
 
                    forAll(tri, triPointI)
                    {
                        const label pointi = tri[triPointI];
 
                        const vector& v = psip[pointi];
 
                        const scalar vHatDotN = normalised(v) & n;
                        maxVHatDotN = max(maxVHatDotN, vHatDotN);
 
                        const vector a =
                            normalised(v ^ pointNormals[pointi]);
 
                        const scalar C = v & n, S = (v ^ a) & n;
 
                        const scalar theta = atan2(C, - S);
 
                        theta0[pointi] = max(theta0[pointi], theta);
                        theta1[pointi] = min(theta1[pointi], theta + pi);
                    }
                }
            }
        }
    }
 
    if (debug)
    {
        reduce(maxVHatDotN, maxOp<scalar>());
        Info<< typeName << ": Maximum in-to-wall dot product before = "
            << maxVHatDotN << endl;
    }
 
    syncTools::syncPointList
    (
        mesh,
        theta0,
        maxEqOp<scalar>(),
        scalar(0)
    );
 
    syncTools::syncPointList
    (
        mesh,
        theta1,
        minEqOp<scalar>(),
        scalar(0)
    );
 
 
    if (debug > 1)
    {
        pointVectorField(psip.name() + "Before", psip).write();
    }
 
 
    // Apply the rotations so that the interpolated field points into the mesh
    forAll(mesh.boundaryMesh(), patchi)
    {
        const polyPatch& patch = mesh.boundaryMesh()[patchi];
        if (isA<wallPolyPatch>(patch))
        {
            forAll(patch.meshPoints(), patchPointi)
            {
                const label pointi = patch.meshPoints()[patchPointi];
 
                vector& v = psip[pointi];
 
                if (theta0[pointi] <= 0 && theta1[pointi] >= 0)
                {
                    continue;
                }
 
                if (theta0[pointi] >= theta1[pointi])
                {
                    v = Zero;
                    continue;
                }
 
                const scalar theta =
                    theta0[pointi] > 0 ? theta0[pointi] : theta1[pointi];
 
                const scalar c = cos(theta), s = sin(theta);
 
                const scalar scale = max(c, 0); // or mag(theta)/(pi/2) ...
 
                const vector a = normalised(v ^ pointNormals[pointi]);
 
                v = scale*(tensor::I*c - (*a)*s + sqr(a)*(1 - c)) & v;
 
                theta0[pointi] = theta1[pointi] = 0;
            }
        }
    }
 
 
    // Report the field-normal dot products
    if (debug)
    {
        maxVHatDotN = - vGreat;
 
        forAll(mesh.boundaryMesh(), patchi)
        {
            const polyPatch& patch = mesh.boundaryMesh()[patchi];
            if (isA<wallPolyPatch>(patch))
            {
                forAll(patch, patchFacei)
                {
                    const label facei = patch.start() + patchFacei;
                    const face& f = patch[patchFacei];
 
                    for (label i = 1; i < f.size() - 1; ++ i)
                    {
                        const triFace tri
                        (
                            tetIndices(mesh.faceOwner()[facei], facei, i)
                           .faceTriIs(mesh)
                        );
 
                        const vector n = tri.normal(mesh.points());
 
                        forAll(tri, triPointI)
                        {
                            const label pointi = tri[triPointI];
 
                            const vector& v = psip[pointi];
 
                            const scalar vHatDotN = normalised(v) & n;
                            maxVHatDotN = max(maxVHatDotN, vHatDotN);
                        }
                    }
                }
            }
        }
 
        reduce(maxVHatDotN, maxOp<scalar>());
        Info<< typeName << ": Maximum in-to-wall dot product after = "
            << maxVHatDotN << endl;
    }
 
 
    // Write the interpolated field
    if (debug > 1)
    {
        psip.write();
    }
 
 
    return tPsip;
}
 
 
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
 
template<class Type>
Foam::interpolationCellPointWallModified<Type>::
interpolationCellPointWallModified
(
    const VolField<Type>& psi
)
:
    interpolationCellPoint<Type>(psi, calcPointField(psi))
{}
 
 
template<class Type>
Foam::interpolationCellPointWallModified<Type>::
interpolationCellPointWallModified
(
    const interpolationCellPointWallModified<Type>& i
)
:
    interpolationCellPoint<Type>(i)
{}
 
 
// ************************************************************************* //