CalcPatchToPatchWeights.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) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2020-2022 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 "PatchToPatchInterpolation.H"
#include "objectHit.H"
#include "pointHit.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

template<class FromPatch, class ToPatch>
void PatchToPatchInterpolation<FromPatch, ToPatch>::calcPointAddressing() const
{
    // Calculate pointWeights

    pointWeightsPtr_.reset(new FieldField<Field, scalar>(toPatch_.nPoints()));
    auto& pointWeights = *pointWeightsPtr_;

    pointDistancePtr_.reset(new scalarField(toPatch_.nPoints(), GREAT));
    auto& pointDistance = *pointDistancePtr_;

    const pointField& fromPatchPoints = fromPatch_.localPoints();
    const auto& fromPatchFaces = fromPatch_.localFaces();

    const pointField& toPatchPoints = toPatch_.localPoints();
    const vectorField& projectionDirection = toPatch_.pointNormals();
    const edgeList& toPatchEdges = toPatch_.edges();
    const labelListList& toPatchPointEdges = toPatch_.pointEdges();

    if (debug)
    {
        Info<< "projecting points" << endl;
    }

    List<objectHit> proj =
        toPatch_.projectPoints(fromPatch_, projectionDirection, alg_, dir_);

    pointAddressingPtr_.reset(new labelList(proj.size(), -1));
    auto& pointAddressing = *pointAddressingPtr_;

    bool doWeights = false;

    forAll(pointAddressing, pointi)
    {
        doWeights = false;

        const auto& hitFace = fromPatchFaces[proj[pointi].hitObject()];

        point hitPoint = Zero;

        if (proj[pointi].hit())
        {
            // A hit exists
            doWeights = true;

            pointAddressing[pointi] = proj[pointi].hitObject();

            pointHit curHit =
                hitFace.ray
                (
                    toPatchPoints[pointi],
                    projectionDirection[pointi],
                    fromPatchPoints,
                    alg_,
                    dir_
                );

            // Grab distance to target
            if (dir_ == intersection::CONTACT_SPHERE)
            {
                pointDistance[pointi] =
                    hitFace.contactSphereDiameter
                    (
                        toPatchPoints[pointi],
                        projectionDirection[pointi],
                        fromPatchPoints
                    );
            }
            else
            {
                pointDistance[pointi] = curHit.distance();
            }

            // Grab hit point
            hitPoint = curHit.hitPoint();
        }
        else if (projectionTol_ > SMALL)
        {
            // Check for a near miss
            pointHit ph =
                hitFace.ray
                (
                    toPatchPoints[pointi],
                    projectionDirection[pointi],
                    fromPatchPoints,
                    alg_,
                    dir_
                );

            scalar dist =
                Foam::mag
                (
                    toPatchPoints[pointi]
                  + projectionDirection[pointi]*ph.distance()
                  - ph.missPoint()
                );

            // Calculate the local tolerance
            scalar minEdgeLength = GREAT;

            // Do shortest edge of hit object
            edgeList hitFaceEdges =
                fromPatchFaces[proj[pointi].hitObject()].edges();

            forAll(hitFaceEdges, edgeI)
            {
                minEdgeLength =
                    min
                    (
                        minEdgeLength,
                        hitFaceEdges[edgeI].mag(fromPatchPoints)
                    );
            }

            const labelList& curEdges = toPatchPointEdges[pointi];

            forAll(curEdges, edgeI)
            {
                minEdgeLength =
                    min
                    (
                        minEdgeLength,
                        toPatchEdges[curEdges[edgeI]].mag(toPatchPoints)
                    );
            }

            if (dist < minEdgeLength*projectionTol_)
            {
                // This point is being corrected
                doWeights = true;

                pointAddressing[pointi] = proj[pointi].hitObject();

                // Grab nearest point on face as hit point
                hitPoint = ph.missPoint();

                // Grab distance to target
                if (dir_ == intersection::CONTACT_SPHERE)
                {
                    pointDistance[pointi] =
                        hitFace.contactSphereDiameter
                        (
                            toPatchPoints[pointi],
                            projectionDirection[pointi],
                            fromPatchPoints
                        );
                }
                else
                {
                    pointDistance[pointi] =
                        (
                            projectionDirection[pointi]
                            /mag(projectionDirection[pointi])
                        )
                      & (hitPoint - toPatchPoints[pointi]);
                }
            }
        }

        if (doWeights)
        {
            // Set interpolation pointWeights
            const pointField hitFacePoints
            (
                hitFace.points(fromPatchPoints)
            );

            auto& pointiWeights =
                pointWeights.emplace_set(pointi, hitFacePoints.size());

            scalar sumWeight = 0;

            forAll(hitFacePoints, masterPointi)
            {
                const point& p = hitFacePoints[masterPointi];

                const scalar w =
                (
                    1.0 / (hitPoint.dist(p) + VSMALL)
                );

                pointiWeights[masterPointi] = w;
                sumWeight += w;
            }

            pointiWeights /= sumWeight;
        }
        else
        {
            pointWeights.emplace_set(pointi);
        }
    }
}


template<class FromPatch, class ToPatch>
void PatchToPatchInterpolation<FromPatch, ToPatch>::calcFaceAddressing() const
{
    faceWeightsPtr_.reset(new FieldField<Field, scalar>(toPatch_.size()));
    auto& faceWeights = *faceWeightsPtr_;

    faceDistancePtr_.reset(new scalarField(toPatch_.size(), GREAT));
    auto& faceDistance = *faceDistancePtr_;

    if (debug)
    {
        Info<< "projecting face centres" << endl;
    }

    const pointField& fromPatchPoints = fromPatch_.points();
    const auto& fromPatchFaces = fromPatch_;
    const labelListList& fromPatchFaceFaces = fromPatch_.faceFaces();

    vectorField fromPatchFaceCentres(fromPatchFaces.size());

    forAll(fromPatchFaceCentres, facei)
    {
        fromPatchFaceCentres[facei] =
            fromPatchFaces[facei].centre(fromPatchPoints);
    }

    const pointField& toPatchPoints = toPatch_.points();
    const auto& toPatchFaces = toPatch_;

    const vectorField& projectionDirection = toPatch_.faceNormals();

    List<objectHit> proj =
        toPatch_.projectFaceCentres
        (
            fromPatch_,
            projectionDirection,
            alg_,
            dir_
        );

    faceAddressingPtr_.reset(new labelList(proj.size(), -1));
    auto& faceAddressing = *faceAddressingPtr_;

    forAll(faceAddressing, facei)
    {
        if (proj[facei].hit())
        {
            // A hit exists
            faceAddressing[facei] = proj[facei].hitObject();

            const auto& hitFace = fromPatchFaces[faceAddressing[facei]];

            pointHit curHit =
                hitFace.ray
                (
                    toPatchFaces[facei].centre(toPatchPoints),
                    projectionDirection[facei],
                    fromPatchPoints,
                    alg_,
                    dir_
                );

            // grab distance to target
            faceDistance[facei] = curHit.distance();

            // grab face centre of the hit face
            const point& hitFaceCentre =
                fromPatchFaceCentres[faceAddressing[facei]];

            // grab neighbours of hit face
            const labelList& neighbours =
                fromPatchFaceFaces[faceAddressing[facei]];

            const point& hitPoint = curHit.hitPoint();

            scalar m = hitPoint.dist(hitFaceCentre);

            if
            (
                m < directHitTol                            // Direct hit
             || neighbours.empty()
            )
            {
                auto& faceiWeights = faceWeights.emplace_set(facei, 1);
                faceiWeights[0] = 1.0;
            }
            else
            {
                // set interpolation faceWeights

                // The first coefficient corresponds to the centre face.
                // The rest is ordered in the same way as the faceFaces list.

                auto& faceiWeights =
                    faceWeights.emplace_set(facei, neighbours.size() + 1);

                faceiWeights[0] = 1.0/m;

                scalar sumWeight = faceiWeights[0];

                forAll(neighbours, nbri)
                {
                    const point& p = fromPatchFaceCentres[neighbours[nbri]];

                    const scalar w =
                    (
                        1.0 / (hitPoint.dist(p) + VSMALL)
                    );

                    faceWeights[nbri+1] = w;
                    sumWeight += w;
                }

                faceiWeights /= sumWeight;
            }
        }
        else
        {
            faceWeights.emplace_set(facei);
        }
    }
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //