MappedFileFilterField.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) 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 "MappedFileFilterField.H"
#include "oneField.H"
#include "MeshedSurfaces.H"
#include "MinMax.H"
#include "indexedOctree.H"
#include "treeDataPoint.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
namespace PatchFunction1Types
{
    defineTypeNameAndDebug(FilterField, 0);
}
}
 
 
const Foam::Enum
<
    Foam::PatchFunction1Types::FilterField::RBF_type
>
Foam::PatchFunction1Types::FilterField::RBF_typeNames_
{
    { RBF_type::RBF_linear, "linear" },
    { RBF_type::RBF_quadratic, "quadratic" },
    { RBF_type::RBF_linear, "default" },
};
 
 
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
namespace Foam
{
 
//- Linear basis function: 1 - (r/r0)
static inline scalar linearWeight
(
    const point& p,
    const point& p0,
    const scalar radiusSqr
)
{
    return (1 - ::sqrt(p.distSqr(p0) / radiusSqr));
}
 
 
//- Quadratic basis function: 1 - (r/r0)^2
static inline scalar quadraticWeight
(
    const point& p,
    const point& p0,
    const scalar radiusSqr
)
{
    return (1 - p.distSqr(p0) / radiusSqr);
}
 
 
//- Construct search tree for points
static autoPtr<indexedOctree<treeDataPoint>>
createTree(const pointField& points)
{
    // Avoid bounding box error in case of 2D cases
    treeBoundBox bb(points);
    bb.grow(1e-4);
 
    const int debug = PatchFunction1Types::FilterField::debug;
 
    const int oldDebug = indexedOctreeBase::debug;
 
    if (debug & 2)
    {
        indexedOctreeBase::debug = 1;
    }
 
    auto treePtr = autoPtr<indexedOctree<treeDataPoint>>::New
    (
        treeDataPoint(points),
        bb,             // overall search domain
        points.size(),  // maxLevel
        16,             // leafsize
        1.0             // duplicity
    );
 
    indexedOctreeBase::debug = oldDebug;
 
    if (debug & 2)
    {
        const auto& tree = treePtr();
 
        OFstream os("indexedOctree.obj");
        tree.writeOBJ(os);
        Info<< "=================" << endl;
        Info<< "have " << tree.nodes().size() << " nodes" << nl;
        Info<< "=================" << endl;
    }
 
    return treePtr;
}
 
} // End namespace Foam
 
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
template
<
    class TreeType,
    class RadiusSqrOp,
    class BasisFunctionOp,
    class PointWeightFieldType
>
void Foam::PatchFunction1Types::FilterField::buildWeightsImpl
(
    const TreeType& tree,
    const RadiusSqrOp& radiusSqrOp,
    const BasisFunctionOp& basisFuncOp,
    const PointWeightFieldType& pointWeightFld
)
{
    tmp<pointField> tpoints = tree.shapes().centres();
 
    const auto& points = tpoints();
 
    const label nPoints = points.size();
 
    weights_.clear();
    addressing_.clear();
 
    weights_.resize(nPoints);
    addressing_.resize(nPoints);
 
    labelHashSet neighbours(2*128);
 
    // Catch degenerate case where weight/addressing not actually needed
    bool usesNeighbours = false;
 
    for (label pointi = 0; pointi < nPoints; ++pointi)
    {
        const point& p0 = points[pointi];
        auto& currAddr = addressing_[pointi];
        auto& currWeight = weights_[pointi];
 
        // Search with local sphere
        const scalar radiusSqr = radiusSqrOp(pointi);
 
        tree.findSphere(p0, radiusSqr, neighbours);
 
        // Paranoid, enforce identity weighting as a minimum
        if (neighbours.size() < 2)
        {
            currAddr.resize(1, pointi);
            currWeight.resize(1, scalar(1));
            continue;
        }
 
        usesNeighbours = true;
 
        currAddr = neighbours.sortedToc();
        currWeight.resize_nocopy(currAddr.size());
 
        scalar sumWeight = 0;
 
        forAll(currAddr, i)
        {
            const point& p = points[currAddr[i]];
 
            currWeight[i] = basisFuncOp(p, p0, radiusSqr);
 
            // Eg, face area weighting.
            // - cast required for oneField()
            currWeight[i] *= static_cast<scalar>(pointWeightFld[i]);
            sumWeight += currWeight[i];
        }
 
        for (auto& w : currWeight)
        {
            w /= sumWeight;
        }
    }
 
    if (!usesNeighbours)
    {
        // No addressing/weighting required
        reset();
    }
 
    if (debug && !addressing_.empty())
    {
        labelMinMax limits(addressing_[0].size());
 
        label total = 0;
 
        for (const labelList& addr : addressing_)
        {
            const label n = addr.size();
 
            limits.add(n);
            total += n;
        }
 
        Pout<< "Weight neighbours: min=" << limits.min()
            << " avg=" << (total / addressing_.size())
            << " max=" << limits.max() << endl;
    }
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::PatchFunction1Types::FilterField::FilterField
(
    const pointField& points,
    const scalar radius,
    const RBF_type interp
)
{
    reset(points, radius, interp);
}
 
 
Foam::PatchFunction1Types::FilterField::FilterField
(
    const meshedSurface& geom,
    const scalar radius,
    const bool relative,
    const RBF_type interp
)
{
    reset(geom, radius, relative, interp);
}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
void Foam::PatchFunction1Types::FilterField::reset()
{
    addressing_.clear();
    weights_.clear();
}
 
 
void Foam::PatchFunction1Types::FilterField::reset
(
    const pointField& points,
    const scalar radius,
    const RBF_type interp
)
{
    scalarField dummyWeights;
 
    if (debug)
    {
        Pout<< "Apply " << RBF_typeNames_[interp] << " filter,"
            << " radius=" << radius << nl
            << "Create tree..." << endl;
    }
 
    autoPtr<indexedOctree<treeDataPoint>> treePtr
        = createTree(points);
 
    const scalar radiusSqr = sqr(radius);
 
    auto weightFunc = linearWeight;
 
    if (interp == RBF_type::RBF_quadratic)
    {
        weightFunc = quadraticWeight;
    }
 
 
    // Use (constant) absolute size
    buildWeightsImpl
    (
        treePtr(),
        [=](label index) -> scalar { return radiusSqr; },
        weightFunc,
        oneField()
    );
}
 
 
void Foam::PatchFunction1Types::FilterField::reset
(
    const meshedSurface& geom,
    const scalar radius,
    const bool relative,
    const RBF_type interp
)
{
    if (!geom.nFaces())
    {
        // Received geometry without any faces eg, boundaryData
        reset(geom.points(), radius, interp);
        return;
    }
 
    if (debug)
    {
        Pout<< "Apply " << RBF_typeNames_[interp] << " filter,"
            << (relative ? " relative" : "") << " radius=" << radius << nl
            << "Create tree..." << endl;
    }
 
    autoPtr<indexedOctree<treeDataPoint>> treePtr
        = createTree(geom.faceCentres());
 
    const scalar radiusSqr = sqr(radius);
 
    auto weightFunc = linearWeight;
 
    if (interp == RBF_type::RBF_quadratic)
    {
        weightFunc = quadraticWeight;
    }
 
 
    if (relative)
    {
        // Use relative size
        buildWeightsImpl
        (
            treePtr(),
            [&](label index) -> scalar
            {
                return (radiusSqr * geom.sphere(index));
            },
            weightFunc,
            geom.magSf()
        );
    }
    else
    {
        // Use (constant) absolute size
        buildWeightsImpl
        (
            treePtr(),
            [=](label index) -> scalar { return radiusSqr; },
            weightFunc,
            geom.magSf()
        );
    }
}
 
 
// ************************************************************************* //