ENH: patchCloudSet: new sampledSet for nearest point on cloud

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -49,6 +49,9 @@ formatOptions
 //      cell          : use cell-centre value only; constant over cells (default)
 //      cellPoint     : use cell-centre and vertex values
 //      cellPointFace : use cell-centre, vertex and face values.
+//      pointMVC      : use point values only (Mean Value Coordinates)
+//      cellPatchConstrained : use cell-centre except on boundary faces where
+//        it uses the boundary value. For use with e.g. patchCloudSet.
 // 1] vertex values determined from neighbouring cell-centre values
 // 2] face values determined using the current face interpolation scheme
 //    for the field (linear, gamma, etc.)
@@ -83,6 +86,7 @@ fields
 //      uniform, face, midPoint, midPointAndFace : start and end coordinate
 //      uniform: extra number of sampling points
 //      polyLine, cloud: list of coordinates
+//      patchCloud: list of coordinates and set of patches to look for nearest
 sets
 (
     lineX1
@@ -113,8 +117,21 @@ sets
         points  ((0.049 0.049 0.00501)(0.051 0.049 0.00501));
     }
 
+    somePatchPoints
+    {
+        // Sample nearest points on selected patches. Use with
+        // interpolations:
+        // - cell (cell value)
+        // - cellPatchConstrained (boundary value)
+        // - cellPoint (interpolated boundary value)
+        type    patchCloud;
+        axis    xyz;
+        points  ((0.049 0.099 0.005)(0.051 0.054 0.005));
+        patches (".*Wall.*");
+    }
 );
 
+
 // Surface sampling definition
 //
 // 1] patches are not triangulated by default
@@ -241,4 +258,5 @@ surfaces
     }
 );
 
+
 // *********************************************************************** //

src/finiteVolume/Make/files

@@ -201,6 +201,7 @@ interpolation = interpolation/interpolation
 $(interpolation)/interpolation/interpolations.C
 
 $(interpolation)/interpolationCell/makeInterpolationCell.C
+$(interpolation)/interpolationCellPatchConstrained/makeInterpolationCellPatchConstrained.C
 $(interpolation)/interpolationCellPoint/cellPointWeight/cellPointWeight.C
 $(interpolation)/interpolationCellPoint/makeInterpolationCellPoint.C
 $(interpolation)/interpolationCellPointFace/makeInterpolationCellPointFace.C

src/finiteVolume/interpolation/interpolation/interpolationCellPatchConstrained/interpolationCellPatchConstrained.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+     \\/     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 "interpolationCellPatchConstrained.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
+
+template<class Type>
+interpolationCellPatchConstrained<Type>::interpolationCellPatchConstrained
+(
+    const GeometricField<Type, fvPatchField, volMesh>& psi
+)
+:
+    interpolation<Type>(psi)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Type interpolationCellPatchConstrained<Type>::interpolate
+(
+    const vector& pt,
+    const label cellI,
+    const label faceI
+) const
+{
+    if (faceI >= 0 && faceI >= this->psi_.mesh().nInternalFaces())
+    {
+        // Use boundary value
+        const polyBoundaryMesh& pbm = this->psi_.mesh().boundaryMesh();
+        label patchI = pbm.patchID()[faceI-this->psi_.mesh().nInternalFaces()];
+        label patchFaceI = pbm[patchI].whichFace(faceI);
+
+        return this->psi_.boundaryField()[patchI][patchFaceI];
+    }
+    else
+    {
+        return this->psi_[cellI];
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/interpolation/interpolation/interpolationCellPatchConstrained/interpolationCellPatchConstrained.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     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/>.
+
+Class
+    Foam::interpolationCellPatchConstrained
+
+Description
+    Uses the cell value for any point in the cell apart from a boundary face
+    where it uses the boundary value directly.
+    Note: will not work on an empty patch.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interpolationCellPatchConstrained_H
+#define interpolationCellPatchConstrained_H
+
+#include "interpolation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class interpolationCellPatchConstrained Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class interpolationCellPatchConstrained
+:
+    public interpolation<Type>
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("cellPatchConstrained");
+
+
+    // Constructors
+
+        //- Construct from components
+        interpolationCellPatchConstrained
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& psi
+        );
+
+
+    // Member Functions
+
+        //- Interpolate field to the given point in the given cell
+        Type interpolate
+        (
+            const vector& position,
+            const label cellI,
+            const label faceI = -1
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "interpolationCellPatchConstrained.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/interpolation/interpolation/interpolationCellPatchConstrained/makeInterpolationCellPatchConstrained.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     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 "interpolationCellPatchConstrained.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeInterpolation(interpolationCellPatchConstrained);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/sampling/Make/files

@@ -5,6 +5,7 @@ probes/probesFunctionObject/probesFunctionObject.C
 
 sampledSet/cloud/cloudSet.C
 sampledSet/coordSet/coordSet.C
+sampledSet/patchCloud/patchCloudSet.C
 sampledSet/polyLine/polyLineSet.C
 sampledSet/face/faceOnlySet.C
 sampledSet/midPoint/midPointSet.C

src/sampling/sampledSet/patchCloud/patchCloudSet.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     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 "patchCloudSet.H"
+#include "polyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointIndexHit.H"
+#include "Tuple2.H"
+#include "treeBoundBox.H"
+#include "indexedOctree.H"
+#include "treeDataFace.H"
+#include "Time.H"
+#include "meshTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(patchCloudSet, 0);
+    addToRunTimeSelectionTable(sampledSet, patchCloudSet, word);
+
+
+    //- Helper class for finding nearest
+    // Nearest:
+    //  - point+local index
+    //  - sqr(distance)
+    //  - processor
+    typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
+
+    class nearestEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first().hit())
+            {
+                if (!x.first().hit())
+                {
+                    x = y;
+                }
+                else if (y.second().first() < x.second().first())
+                {
+                    x = y;
+                }
+            }
+        }
+    };
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::patchCloudSet::calcSamples
+(
+    DynamicList<point>& samplingPts,
+    DynamicList<label>& samplingCells,
+    DynamicList<label>& samplingFaces,
+    DynamicList<label>& samplingSegments,
+    DynamicList<scalar>& samplingCurveDist
+) const
+{
+    if (debug)
+    {
+        Info<< "patchCloudSet : sampling on patches :" << endl;
+    }
+
+    // Construct search tree for all patch faces.
+    label sz = 0;
+    forAllConstIter(labelHashSet, patchSet_, iter)
+    {
+        const polyPatch& pp = mesh().boundaryMesh()[iter.key()];
+
+        sz += pp.size();
+
+        if (debug)
+        {
+            Info<< "    " << pp.name() << " size " << pp.size() << endl;
+        }
+    }
+
+    labelList patchFaces(sz);
+    sz = 0;
+    treeBoundBox bb(point::max, point::min);
+    forAllConstIter(labelHashSet, patchSet_, iter)
+    {
+        const polyPatch& pp = mesh().boundaryMesh()[iter.key()];
+
+        forAll(pp, i)
+        {
+            patchFaces[sz++] = pp.start()+i;
+        }
+
+        const boundBox patchBb(pp.points(), pp.meshPoints());
+
+        bb.min() = min(bb.min(), patchBb.min());
+        bb.max() = max(bb.max(), patchBb.max());
+    }
+
+    // Not very random
+    Random rndGen(123456);
+    // Make bb asymetric just to avoid problems on symmetric meshes
+    bb = bb.extend(rndGen, 1E-4);
+
+    bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+
+    indexedOctree<treeDataFace> patchTree
+    (
+        treeDataFace    // all information needed to search faces
+        (
+            false,      // do not cache bb
+            mesh(),
+            patchFaces  // boundary faces only
+        ),
+        bb,             // overall search domain
+        8,              // maxLevel
+        10,             // leafsize
+        3.0             // duplicity
+    );
+
+
+
+    // All the info for nearest. Construct to miss
+    List<nearInfo> nearest(sampleCoords_.size());
+
+    forAll(sampleCoords_, sampleI)
+    {
+        const point& sample = sampleCoords_[sampleI];
+
+        pointIndexHit& nearInfo = nearest[sampleI].first();
+
+        // Find the nearest locally
+        nearInfo = patchTree.findNearest(sample, magSqr(bb.span()));
+
+        // Fill in the distance field and the processor field
+        if (!nearInfo.hit())
+        {
+            nearest[sampleI].second().first() = Foam::sqr(GREAT);
+            nearest[sampleI].second().second() = Pstream::myProcNo();
+        }
+        else
+        {
+            // Set nearest to mesh face label
+            nearInfo.setIndex(patchFaces[nearInfo.index()]);
+
+            nearest[sampleI].second().first() = magSqr
+            (
+                nearInfo.hitPoint()
+              - sample
+            );
+            nearest[sampleI].second().second() = Pstream::myProcNo();
+        }
+    }
+
+
+    // Find nearest.
+    Pstream::listCombineGather(nearest, nearestEqOp());
+    Pstream::listCombineScatter(nearest);
+
+
+    if (debug && Pstream::master())
+    {
+        OFstream str
+        (
+            mesh().time().path()
+          / name()
+          + "_nearest.obj"
+        );
+        Info<< "Dumping mapping as lines from supplied points to"
+            << " nearest patch face to file " << str.name() << endl;
+
+        label vertI = 0;
+
+        forAll(nearest, i)
+        {
+            meshTools::writeOBJ(str, sampleCoords_[i]);
+            vertI++;
+            meshTools::writeOBJ(str, nearest[i].first().hitPoint());
+            vertI++;
+            str << "l " << vertI-1 << ' ' << vertI << nl;
+        }
+    }
+
+
+    // Store the sampling locations on the nearest processor
+    forAll(nearest, sampleI)
+    {
+        if (nearest[sampleI].second().second() == Pstream::myProcNo())
+        {
+            const pointIndexHit& nearInfo = nearest[sampleI].first();
+
+            label faceI = nearInfo.index();
+
+            samplingPts.append(nearInfo.hitPoint());
+            samplingCells.append(mesh().faceOwner()[faceI]);
+            samplingFaces.append(faceI);
+            samplingSegments.append(0);
+            samplingCurveDist.append(1.0 * sampleI);
+        }
+    }
+}
+
+
+void Foam::patchCloudSet::genSamples()
+{
+    // Storage for sample points
+    DynamicList<point> samplingPts;
+    DynamicList<label> samplingCells;
+    DynamicList<label> samplingFaces;
+    DynamicList<label> samplingSegments;
+    DynamicList<scalar> samplingCurveDist;
+
+    calcSamples
+    (
+        samplingPts,
+        samplingCells,
+        samplingFaces,
+        samplingSegments,
+        samplingCurveDist
+    );
+
+    samplingPts.shrink();
+    samplingCells.shrink();
+    samplingFaces.shrink();
+    samplingSegments.shrink();
+    samplingCurveDist.shrink();
+
+    setSamples
+    (
+        samplingPts,
+        samplingCells,
+        samplingFaces,
+        samplingSegments,
+        samplingCurveDist
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchCloudSet::patchCloudSet
+(
+    const word& name,
+    const polyMesh& mesh,
+    meshSearch& searchEngine,
+    const word& axis,
+    const List<point>& sampleCoords,
+    const labelHashSet& patchSet
+)
+:
+    sampledSet(name, mesh, searchEngine, axis),
+    sampleCoords_(sampleCoords),
+    patchSet_(patchSet)
+{
+    genSamples();
+
+    if (debug)
+    {
+        write(Info);
+    }
+}
+
+
+Foam::patchCloudSet::patchCloudSet
+(
+    const word& name,
+    const polyMesh& mesh,
+    meshSearch& searchEngine,
+    const dictionary& dict
+)
+:
+    sampledSet(name, mesh, searchEngine, dict),
+    sampleCoords_(dict.lookup("points")),
+    patchSet_
+    (
+        mesh.boundaryMesh().patchSet
+        (
+            wordReList(dict.lookup("patches"))
+        )
+    )
+{
+    genSamples();
+
+    if (debug)
+    {
+        write(Info);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::patchCloudSet::~patchCloudSet()
+{}