View Issue Details

IDProjectCategoryView StatusLast Update
0001147OpenFOAMBugpublic2018-10-03 09:19
Reporteruser672Assigned Towill  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
PlatformLinuxOSScientific LinuxOS Version6.4
Fixed in Versiondev 
Summary0001147: Sample utility: line passing through internal cyclic boundary using midPoint set sampling definition isn't possible
DescriptionNo error happens if "face" set sampling definition is chosen. Same error for either a sampling functionObject, or using the postProcessing sample utility.
Steps To ReproduceIn the attached minimal case I've set up a square baffle in the centre of a channel of square cross-section, and a sampling functionObject to output the velocity along two axial lines; one going through the baffle and one running past it. The line running through the baffle fails if midPoint is used.
Additional InformationSame behaviour observed for both GCC-4.6 and Clang-4.3 builds.
Tagsbaffle, sample

Activities

user672

2014-02-01 19:47

 

sampleCyclic.tgz (3,373 bytes)

wyldckat

2015-03-29 21:29

updater  

bug1147.patch (12,741 bytes)   
diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index c71559d..861b7ec 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -61,7 +61,8 @@ Foam::particle::particle
     tetFaceI_(tetFaceI),
     tetPtI_(tetPtI),
     origProc_(Pstream::myProcNo()),
-    origId_(getNewParticleID())
+    origId_(getNewParticleID()),
+    markedToTeleport_(false)
 {}
 
 
@@ -81,7 +82,8 @@ Foam::particle::particle
     tetFaceI_(-1),
     tetPtI_(-1),
     origProc_(Pstream::myProcNo()),
-    origId_(getNewParticleID())
+    origId_(getNewParticleID()),
+    markedToTeleport_(false)
 {
     if (doCellFacePt)
     {
@@ -100,7 +102,8 @@ Foam::particle::particle(const particle& p)
     tetFaceI_(p.tetFaceI_),
     tetPtI_(p.tetPtI_),
     origProc_(p.origProc_),
-    origId_(p.origId_)
+    origId_(p.origId_),
+    markedToTeleport_(false)
 {}
 
 
@@ -114,7 +117,8 @@ Foam::particle::particle(const particle& p, const polyMesh& mesh)
     tetFaceI_(p.tetFaceI_),
     tetPtI_(p.tetPtI_),
     origProc_(p.origProc_),
-    origId_(p.origId_)
+    origId_(p.origId_),
+    markedToTeleport_(false)
 {}
 
 
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index 6cb1838..9d9b099 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -154,6 +154,8 @@ protected:
 
         //- Local particle id on originating processor
         label origId_;
+        
+        bool markedToTeleport_;
 
 
     // Private Member Functions
@@ -457,6 +459,9 @@ public:
             //- Is the particle on the boundary/(or outside the domain)?
             inline bool onBoundary() const;
 
+            inline bool needsTeleportHandled() const;
+            inline void resetTeleportHandler();
+            
             //- Is this global face an internal face?
             inline bool internalFace(const label faceI) const;
 
@@ -511,7 +516,7 @@ public:
             //  on entry 'stepFraction()' should be set to the fraction of the
             //  time-step at which the tracking starts.
             template<class TrackData>
-            scalar trackToFace(const vector& endPosition, TrackData& td);
+            scalar trackToFace(const vector& endPosition, TrackData& td, bool handleTeleport=false);
 
             //- Return the index of the face to be used in the interpolation
             //  routine
diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
index 803a200..38a0d0b 100644
--- a/src/lagrangian/basic/particle/particleI.H
+++ b/src/lagrangian/basic/particle/particleI.H
@@ -850,7 +850,19 @@ inline void Foam::particle::initCellFacePt()
 
 inline bool Foam::particle::onBoundary() const
 {
-    return faceI_ != -1 && faceI_ >= mesh_.nInternalFaces();
+    return faceI_ != -1 && faceI_ >= mesh_.nInternalFaces() && !markedToTeleport_;
+}
+
+
+inline bool Foam::particle::needsTeleportHandled() const
+{
+  return markedToTeleport_;
+}
+
+
+inline void Foam::particle::resetTeleportHandler()
+{
+  markedToTeleport_ = false;
 }
 
 
diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C
index 280165d..7e44e58 100644
--- a/src/lagrangian/basic/particle/particleTemplates.C
+++ b/src/lagrangian/basic/particle/particleTemplates.C
@@ -203,7 +203,8 @@ template<class TrackData>
 Foam::scalar Foam::particle::trackToFace
 (
     const vector& endPosition,
-    TrackData& td
+    TrackData& td,
+    bool handleTeleport
 )
 {
     typedef typename TrackData::cloudType cloudType;
@@ -448,6 +449,9 @@ Foam::scalar Foam::particle::trackToFace
             // This must be a cell face crossing
             faceI_ = tetFaceI_;
 
+            Info << "tetFaceI_: " << tetFaceI_ << endl;
+            Info << "tri, lambdaMin: " << triI << " " << lambdaMin << endl;
+
             // Set the faceHitTetIs to those for the current tet in case a
             // wall interaction is required with the cell face
             faceHitTetIs = tetIndices
@@ -576,14 +580,24 @@ Foam::scalar Foam::particle::trackToFace
                 patchI = patch(faceI_);
             }
 
+            Info << "origFaceI: " << origFaceI << endl;
+            Info << "faceI_: " << faceI_ << endl;
+            
             const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 
             if (isA<wedgePolyPatch>(patch))
             {
-                p.hitWedgePatch
-                (
-                    static_cast<const wedgePolyPatch&>(patch), td
-                );
+                if(!handleTeleport || markedToTeleport_)
+                {
+                    p.hitWedgePatch
+                    (
+                        static_cast<const wedgePolyPatch&>(patch), td
+                    );
+                }
+                else
+                {
+                    markedToTeleport_ = true;
+                }
             }
             else if (isA<symmetryPlanePolyPatch>(patch))
             {
@@ -601,26 +615,47 @@ Foam::scalar Foam::particle::trackToFace
             }
             else if (isA<cyclicPolyPatch>(patch))
             {
-                p.hitCyclicPatch
-                (
-                    static_cast<const cyclicPolyPatch&>(patch), td
-                );
+                if(!handleTeleport || markedToTeleport_)
+                {
+                  p.hitCyclicPatch
+                  (
+                      static_cast<const cyclicPolyPatch&>(patch), td
+                  );
+                }
+                else
+                {
+                    markedToTeleport_ = true;
+                }
             }
             else if (isA<cyclicAMIPolyPatch>(patch))
             {
-                p.hitCyclicAMIPatch
-                (
-                    static_cast<const cyclicAMIPolyPatch&>(patch),
-                    td,
-                    endPosition - position_
-                );
+                if(!handleTeleport || markedToTeleport_)
+                {
+                    p.hitCyclicAMIPatch
+                    (
+                        static_cast<const cyclicAMIPolyPatch&>(patch),
+                        td,
+                        endPosition - position_
+                    );
+                }
+                else
+                {
+                    markedToTeleport_ = true;
+                }
             }
             else if (isA<processorPolyPatch>(patch))
             {
-                p.hitProcessorPatch
-                (
-                    static_cast<const processorPolyPatch&>(patch), td
-                );
+                if(!handleTeleport || markedToTeleport_)
+                {
+                    p.hitProcessorPatch
+                    (
+                        static_cast<const processorPolyPatch&>(patch), td
+                    );
+                }
+                else
+                {
+                    markedToTeleport_ = true;
+                }
             }
             else if (isA<wallPolyPatch>(patch))
             {
@@ -633,9 +668,14 @@ Foam::scalar Foam::particle::trackToFace
             {
                 p.hitPatch(patch, td);
             }
+
+            Info << "faceI_ After: " << faceI_ << endl;
+          
         }
     }
 
+    Info << "A faceI_: " << faceI_ << endl;
+    
     if (lambdaMin < SMALL)
     {
         // Apply tracking correction towards tet centre.
@@ -708,6 +748,8 @@ Foam::scalar Foam::particle::trackToFace
         cloud.trackingRescue();
     }
 
+    Info << "V faceI_: " << faceI_ << endl;
+
     return trackFraction;
 }
 
diff --git a/src/sampling/sampledSet/face/faceOnlySet.C b/src/sampling/sampledSet/face/faceOnlySet.C
index 85fb7bc..8d0d19b 100644
--- a/src/sampling/sampledSet/face/faceOnlySet.C
+++ b/src/sampling/sampledSet/face/faceOnlySet.C
@@ -60,14 +60,19 @@ bool Foam::faceOnlySet::trackToBoundary
 
     // Alias
     const point& trackPt = singleParticle.position();
+    bool mustResetTeleport = false;
 
     while(true)
     {
         point oldPoint = trackPt;
 
-        singleParticle.trackToFace(end_, trackData);
-
-        if (singleParticle.face() != -1 && mag(oldPoint - trackPt) > smallDist)
+        singleParticle.trackToFace(end_, trackData, true);
+        
+        if (
+            singleParticle.face() != -1 && 
+            (mag(oldPoint - trackPt) > smallDist
+            || singleParticle.needsTeleportHandled())
+           )
         {
             // Reached face. Sample.
             samplingPts.append(trackPt);
@@ -86,6 +91,18 @@ bool Foam::faceOnlySet::trackToBoundary
             // Boundary reached.
             return true;
         }
+        if(singleParticle.needsTeleportHandled())
+        {
+            if(mustResetTeleport)
+            {
+                singleParticle.resetTeleportHandler();
+                mustResetTeleport = false;
+            }
+            else
+            {
+                mustResetTeleport = true;
+            }
+        }
     }
 }
 
diff --git a/src/sampling/sampledSet/midPoint/midPointSet.C b/src/sampling/sampledSet/midPoint/midPointSet.C
index 168d0bb..7e76bfe 100644
--- a/src/sampling/sampledSet/midPoint/midPointSet.C
+++ b/src/sampling/sampledSet/midPoint/midPointSet.C
@@ -62,11 +62,18 @@ void Foam::midPointSet::genSamples()
         {
             midPoints[midI] =
                 0.5*(operator[](sampleI) + operator[](sampleI+1));
+                
+            Info << "faces_[sampleI]: " << faces_[sampleI] << endl;
+            Info << "faces_[sampleI+1]: " << faces_[sampleI+1] << endl;
 
             label cell1 = getCell(faces_[sampleI], midPoints[midI]);
             label cell2 = getCell(faces_[sampleI+1], midPoints[midI]);
 
-            if (cell1 != cell2)
+            if (cell1 == -1 || cell2 == -1)
+            {
+               Info << "Ignoring intermediate jump" << endl;
+            }
+            else if (cell1 != cell2)
             {
                 FatalErrorIn("midPointSet::genSamples()")
                     << "  sampleI:" << sampleI
@@ -80,12 +87,14 @@ void Foam::midPointSet::genSamples()
                     << "  cell2:" << cell2
                     << abort(FatalError);
             }
+            else
+            {
+                midCells[midI] = cell1;
+                midSegments[midI] = segments_[sampleI];
+                midCurveDist[midI] = mag(midPoints[midI] - start());
+                midI++;
+            }
 
-            midCells[midI] = cell1;
-            midSegments[midI] = segments_[sampleI];
-            midCurveDist[midI] = mag(midPoints[midI] - start());
-
-            midI++;
             sampleI++;
         }
 
diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.C b/src/sampling/sampledSet/sampledSet/sampledSet.C
index 1442ca3..2283610 100644
--- a/src/sampling/sampledSet/sampledSet/sampledSet.C
+++ b/src/sampling/sampledSet/sampledSet/sampledSet.C
@@ -30,6 +30,11 @@ License
 #include "writer.H"
 #include "particle.H"
 
+#include "cyclicPolyPatch.H"
+#include "cyclicAMIPolyPatch.H"
+#include "processorPolyPatch.H"
+#include "wedgePolyPatch.H"
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -64,18 +69,33 @@ Foam::label Foam::sampledSet::getCell
             << abort(FatalError);
     }
 
-    if (faceI >= mesh().nInternalFaces())
+    if (!mesh().isInternalFace(faceI))
     {
         label cellI = getBoundaryCell(faceI);
 
         if (!mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
         {
-            FatalErrorIn
-            (
-                "sampledSet::getCell(const label, const point&)"
-            )   << "Found cell " << cellI << " using face " << faceI
-                << ". But cell does not contain point " << sample
-                << abort(FatalError);
+            label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+            const polyPatch& patch = mesh_.boundaryMesh()[patchI];
+            if(
+               isA<wedgePolyPatch>(patch) || 
+               isA<cyclicPolyPatch>(patch) || 
+               isA<cyclicAMIPolyPatch>(patch) ||
+               isA<processorPolyPatch>(patch)
+              )
+            {
+                return -1;
+            }
+            else
+            {
+
+                FatalErrorIn
+                (
+                    "sampledSet::getCell(const label, const point&)"
+                )   << "Found cell " << cellI << " using face " << faceI
+                    << ". But cell does not contain point " << sample
+                    << abort(FatalError);
+            }
         }
         return cellI;
     }
bug1147.patch (12,741 bytes)   

wyldckat

2015-03-29 21:53

updater   ~0004549

@wallace: Many thanks for the simple and quick test case! It's really just a matter of following your instructions!


@Henry: Attached is the file "bug1147.patch", which one crazy hack of a patch. Seriously consider it only a proof of concept, because I couldn't code anything barely more elegant than that in a short time.
And sorry, but I left in a few stray debug lines, to assess if things were working as intended.

The explanation for the reported bug is simple: the method "particle::trackToFace()" always teleports the tracking point to the other side of the cyclic baffle, which after going through the all of the correct steps, it results in "midPointSet::genSamples()" crashing midway of "sampledSet::getCell()", because the second face is always on the other side of the cyclic patch pair, due to the teleport.


Hopefully the following step-by-step makes it a bit clearer, with the details on how to the provided patch tries to overcome this problem:

  1- The story begins in "faceOnlySet::trackToBoundary()", which is where the collection of faces crossed by the tracking line is gathered. The hack here was to use the new option to properly handle teleporting, where:

    a- In the first hit, it stores the first face before teleport, along with a flag indicating that a teleport needs to be performed.

    b- In the second hit, it executes the teleport, stores the face from the other side and resets the flags related to the need for handling a teleport.

    c- The tricky part is how "faceOnlySet::trackToBoundary()" and "particle::trackToFace()" handle the teleportation. This is because in "particle::trackToFace()" it has to postpone the jump to the second pass, which is why it has to rely so heavily on the flag "markedToTeleport_".

  2- Later on, in "midPointSet::genSamples()" we have to handle the new scenario of "cell1" or "cell2" being -1, which occurs when it's a situation of having the two faces be related to the same special pair of patches.

  3- In "sampledSet::getCell()", when it's not an internal face and if it fails to find the cell the point belongs to, then double-check if the "faceI" is part of one of the special patches that do particle teleportation... which if true, then return "-1".

henry

2015-03-29 22:03

manager   ~0004550

This is a lot of change for this particular purpose. I am also not happy about the increase in storage for ALL particle/parcel types even though it is only needed for this particular post-processing use. Actually this is a general problem with the Lagrangian library that all data which might be possibly used for any application/model is stored even if it is not needed for the particular case. My feeling is that we need more levels in the class hierachy to handle all the options, i.e. in this case we need a derived variant of particle which stores and handles the "teleportation" data appropriately. I think this will need to be delayed until we start writing the new robust tracking algorithm which handles bad cell which decompose with negative volume pyramids.

MattijsJ

2018-09-28 11:29

reporter   ~0010076

Attached fixes the hanging in getting the boundary intersections (meshSearch.[CH]). It removes any geometric offset but instead makes sure it only checks a face once.

Thanks Will for the testcase and analysis.
meshSearch.C (20,510 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2018 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 "meshSearch.H"
#include "polyMesh.H"
#include "indexedOctree.H"
#include "DynamicList.H"
#include "demandDrivenData.H"
#include "treeDataCell.H"
#include "treeDataFace.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(meshSearch, 0);

    scalar meshSearch::tol_ = 1e-3;


    class findUniqueIntersectOp
    {
    public:
        const indexedOctree<treeDataFace>& tree_;
        DynamicList<label>& checkedIndices_;

    public:

        findUniqueIntersectOp
        (
            const indexedOctree<treeDataFace>& tree,
            DynamicList<label>& checkedIndices
        )
        :
            tree_(tree),
            checkedIndices_(checkedIndices)
        {}


        //- Calculate intersection of triangle with ray. Sets result
        //  accordingly
        bool operator()
        (
            const label index,
            const point& start,
            const point& end,
            point& intersectionPoint
        ) const
        {
            const treeDataFace& shape = tree_.shapes();
            const primitiveMesh& mesh = shape.mesh();

            if (findIndex(checkedIndices_, index) != -1)
            {
                return false;
            }
            checkedIndices_.append(index);

            const label facei = shape.faceLabels()[index];

            const vector dir(end - start);

            pointHit inter = mesh.faces()[facei].intersection
            (
                start,
                dir,
                mesh.faceCentres()[facei],
                mesh.points(),
                intersection::HALF_RAY
            );

            if (inter.hit() && inter.distance() <= 1)
            {
                // Note: no extra test on whether intersection is in front of us
                // since using half_ray
                intersectionPoint = inter.hitPoint();
                return true;
            }
            else
            {
                return false;
            }
        }
    };
}


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

bool Foam::meshSearch::findNearer
(
    const point& sample,
    const pointField& points,
    label& nearestI,
    scalar& nearestDistSqr
)
{
    bool nearer = false;

    forAll(points, pointi)
    {
        scalar distSqr = magSqr(points[pointi] - sample);

        if (distSqr < nearestDistSqr)
        {
            nearestDistSqr = distSqr;
            nearestI = pointi;
            nearer = true;
        }
    }

    return nearer;
}


bool Foam::meshSearch::findNearer
(
    const point& sample,
    const pointField& points,
    const labelList& indices,
    label& nearestI,
    scalar& nearestDistSqr
)
{
    bool nearer = false;

    forAll(indices, i)
    {
        label pointi = indices[i];

        scalar distSqr = magSqr(points[pointi] - sample);

        if (distSqr < nearestDistSqr)
        {
            nearestDistSqr = distSqr;
            nearestI = pointi;
            nearer = true;
        }
    }

    return nearer;
}


// tree based searching
Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
{
    const indexedOctree<treeDataCell>& tree = cellTree();

    pointIndexHit info = tree.findNearest
    (
        location,
        magSqr(tree.bb().max()-tree.bb().min())
    );

    if (!info.hit())
    {
        info = tree.findNearest(location, Foam::sqr(great));
    }
    return info.index();
}


Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
{
    const vectorField& centres = mesh_.cellCentres();

    label nearestIndex = 0;
    scalar minProximity = magSqr(centres[nearestIndex] - location);

    findNearer
    (
        location,
        centres,
        nearestIndex,
        minProximity
    );

    return nearestIndex;
}


Foam::label Foam::meshSearch::findNearestCellWalk
(
    const point& location,
    const label seedCelli
) const
{
    if (seedCelli < 0)
    {
        FatalErrorInFunction
            << "illegal seedCell:" << seedCelli << exit(FatalError);
    }

    // Walk in direction of face that decreases distance

    label curCelli = seedCelli;
    scalar distanceSqr = magSqr(mesh_.cellCentres()[curCelli] - location);

    bool closer;

    do
    {
        // Try neighbours of curCelli
        closer = findNearer
        (
            location,
            mesh_.cellCentres(),
            mesh_.cellCells()[curCelli],
            curCelli,
            distanceSqr
        );
    } while (closer);

    return curCelli;
}


Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
{
    // Search nearest cell centre.
    const indexedOctree<treeDataCell>& tree = cellTree();

    // Search with decent span
    pointIndexHit info = tree.findNearest
    (
        location,
        magSqr(tree.bb().max()-tree.bb().min())
    );

    if (!info.hit())
    {
        // Search with desperate span
        info = tree.findNearest(location, Foam::sqr(great));
    }


    // Now check any of the faces of the nearest cell
    const vectorField& centres = mesh_.faceCentres();
    const cell& ownFaces = mesh_.cells()[info.index()];

    label nearestFacei = ownFaces[0];
    scalar minProximity = magSqr(centres[nearestFacei] - location);

    findNearer
    (
        location,
        centres,
        ownFaces,
        nearestFacei,
        minProximity
    );

    return nearestFacei;
}


Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const
{
    const vectorField& centres = mesh_.faceCentres();

    label nearestFacei = 0;
    scalar minProximity = magSqr(centres[nearestFacei] - location);

    findNearer
    (
        location,
        centres,
        nearestFacei,
        minProximity
    );

    return nearestFacei;
}


Foam::label Foam::meshSearch::findNearestFaceWalk
(
    const point& location,
    const label seedFacei
) const
{
    if (seedFacei < 0)
    {
        FatalErrorInFunction
            << "illegal seedFace:" << seedFacei << exit(FatalError);
    }

    const vectorField& centres = mesh_.faceCentres();


    // Walk in direction of face that decreases distance

    label curFacei = seedFacei;
    scalar distanceSqr = magSqr(centres[curFacei] - location);

    while (true)
    {
        label betterFacei = curFacei;

        findNearer
        (
            location,
            centres,
            mesh_.cells()[mesh_.faceOwner()[curFacei]],
            betterFacei,
            distanceSqr
        );

        if (mesh_.isInternalFace(curFacei))
        {
            findNearer
            (
                location,
                centres,
                mesh_.cells()[mesh_.faceNeighbour()[curFacei]],
                betterFacei,
                distanceSqr
            );
        }

        if (betterFacei == curFacei)
        {
            break;
        }

        curFacei = betterFacei;
    }

    return curFacei;
}


Foam::label Foam::meshSearch::findCellLinear(const point& location) const
{
    bool cellFound = false;
    label n = 0;

    label celli = -1;

    while ((!cellFound) && (n < mesh_.nCells()))
    {
        if (mesh_.pointInCell(location, n, cellDecompMode_))
        {
            cellFound = true;
            celli = n;
        }
        else
        {
            n++;
        }
    }
    if (cellFound)
    {
        return celli;
    }
    else
    {
        return -1;
    }
}


Foam::label Foam::meshSearch::findCellWalk
(
    const point& location,
    const label seedCelli
) const
{
    if (seedCelli < 0)
    {
        FatalErrorInFunction
            << "illegal seedCell:" << seedCelli << exit(FatalError);
    }

    if (mesh_.pointInCell(location, seedCelli, cellDecompMode_))
    {
        return seedCelli;
    }

    // Walk in direction of face that decreases distance
    label curCelli = seedCelli;
    scalar nearestDistSqr = magSqr(mesh_.cellCentres()[curCelli] - location);

    while(true)
    {
        // Try neighbours of curCelli

        const cell& cFaces = mesh_.cells()[curCelli];

        label nearestCelli = -1;

        forAll(cFaces, i)
        {
            label facei = cFaces[i];

            if (mesh_.isInternalFace(facei))
            {
                label celli = mesh_.faceOwner()[facei];
                if (celli == curCelli)
                {
                    celli = mesh_.faceNeighbour()[facei];
                }

                // Check if this is the correct cell
                if (mesh_.pointInCell(location, celli, cellDecompMode_))
                {
                    return celli;
                }

                // Also calculate the nearest cell
                scalar distSqr = magSqr(mesh_.cellCentres()[celli] - location);

                if (distSqr < nearestDistSqr)
                {
                    nearestDistSqr = distSqr;
                    nearestCelli = celli;
                }
            }
        }

        if (nearestCelli == -1)
        {
            return -1;
        }

        // Continue with the nearest cell
        curCelli = nearestCelli;
    }

    return -1;
}


Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
(
    const point& location,
    const label seedFacei
) const
{
    if (seedFacei < 0)
    {
        FatalErrorInFunction
            << "illegal seedFace:" << seedFacei << exit(FatalError);
    }

    // Start off from seedFacei

    label curFacei = seedFacei;

    const face& f =  mesh_.faces()[curFacei];

    scalar minDist = f.nearestPoint
    (
        location,
        mesh_.points()
    ).distance();

    bool closer;

    do
    {
        closer = false;

        // Search through all neighbouring boundary faces by going
        // across edges

        label lastFacei = curFacei;

        const labelList& myEdges = mesh_.faceEdges()[curFacei];

        forAll(myEdges, myEdgeI)
        {
            const labelList& neighbours = mesh_.edgeFaces()[myEdges[myEdgeI]];

            // Check any face which uses edge, is boundary face and
            // is not curFacei itself.

            forAll(neighbours, nI)
            {
                label facei = neighbours[nI];

                if
                (
                    (facei >= mesh_.nInternalFaces())
                 && (facei != lastFacei)
                )
                {
                    const face& f =  mesh_.faces()[facei];

                    pointHit curHit = f.nearestPoint
                    (
                        location,
                        mesh_.points()
                    );

                    // If the face is closer, reset current face and distance
                    if (curHit.distance() < minDist)
                    {
                        minDist = curHit.distance();
                        curFacei = facei;
                        closer = true;  // a closer neighbour has been found
                    }
                }
            }
        }
    } while (closer);

    return curFacei;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::meshSearch::meshSearch
(
    const polyMesh& mesh,
    const polyMesh::cellDecomposition cellDecompMode
)
:
    mesh_(mesh),
    cellDecompMode_(cellDecompMode)
{
    if
    (
        cellDecompMode_ == polyMesh::FACE_DIAG_TRIS
     || cellDecompMode_ == polyMesh::CELL_TETS)
    {
        // Force construction of face diagonals
        (void)mesh.tetBasePtIs();
    }
}


Foam::meshSearch::meshSearch
(
    const polyMesh& mesh,
    const treeBoundBox& bb,
    const polyMesh::cellDecomposition cellDecompMode
)
:
    mesh_(mesh),
    cellDecompMode_(cellDecompMode)
{
    overallBbPtr_.reset(new treeBoundBox(bb));

    if
    (
        cellDecompMode_ == polyMesh::FACE_DIAG_TRIS
     || cellDecompMode_ == polyMesh::CELL_TETS
    )
    {
        // Force construction of face diagonals
        (void)mesh.tetBasePtIs();
    }
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::meshSearch::~meshSearch()
{
    clearOut();
}


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

const Foam::indexedOctree<Foam::treeDataFace>&
Foam::meshSearch::boundaryTree() const
{
    if (!boundaryTreePtr_.valid())
    {
        //
        // Construct tree
        //

        if (!overallBbPtr_.valid())
        {
            overallBbPtr_.reset
            (
                new treeBoundBox(mesh_.points())
            );

            treeBoundBox& overallBb = overallBbPtr_();

            // Extend slightly and make 3D
            overallBb = overallBb.extend(1e-4);
        }

        // all boundary faces (not just walls)
        labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
        forAll(bndFaces, i)
        {
            bndFaces[i] = mesh_.nInternalFaces() + i;
        }

        boundaryTreePtr_.reset
        (
            new indexedOctree<treeDataFace>
            (
                treeDataFace    // all information needed to search faces
                (
                    false,                      // do not cache bb
                    mesh_,
                    bndFaces                    // boundary faces only
                ),
                overallBbPtr_(),                // overall search domain
                8,                              // maxLevel
                10,                             // leafsize
                3.0                             // duplicity
            )
        );
    }

    return boundaryTreePtr_();
}


const Foam::indexedOctree<Foam::treeDataCell>&
Foam::meshSearch::cellTree() const
{
    if (!cellTreePtr_.valid())
    {
        //
        // Construct tree
        //

        if (!overallBbPtr_.valid())
        {
            overallBbPtr_.reset
            (
                new treeBoundBox(mesh_.points())
            );

            treeBoundBox& overallBb = overallBbPtr_();

            // Extend slightly and make 3D
            overallBb = overallBb.extend(1e-4);
        }

        cellTreePtr_.reset
        (
            new indexedOctree<treeDataCell>
            (
                treeDataCell
                (
                    false,          // not cache bb
                    mesh_,
                    cellDecompMode_ // cell decomposition mode for inside tests
                ),
                overallBbPtr_(),
                8,              // maxLevel
                10,             // leafsize
                6.0             // duplicity
            )
        );
    }

    return cellTreePtr_();
}


Foam::label Foam::meshSearch::findNearestCell
(
    const point& location,
    const label seedCelli,
    const bool useTreeSearch
) const
{
    if (seedCelli == -1)
    {
        if (useTreeSearch)
        {
            return findNearestCellTree(location);
        }
        else
        {
            return findNearestCellLinear(location);
        }
    }
    else
    {
        return findNearestCellWalk(location, seedCelli);
    }
}


Foam::label Foam::meshSearch::findNearestFace
(
    const point& location,
    const label seedFacei,
    const bool useTreeSearch
) const
{
    if (seedFacei == -1)
    {
        if (useTreeSearch)
        {
            return findNearestFaceTree(location);
        }
        else
        {
            return findNearestFaceLinear(location);
        }
    }
    else
    {
        return findNearestFaceWalk(location, seedFacei);
    }
}


Foam::label Foam::meshSearch::findCell
(
    const point& location,
    const label seedCelli,
    const bool useTreeSearch
) const
{
    // Find the nearest cell centre to this location
    if (seedCelli == -1)
    {
        if (useTreeSearch)
        {
            return cellTree().findInside(location);
        }
        else
        {
            return findCellLinear(location);
        }
    }
    else
    {
        return findCellWalk(location, seedCelli);
    }
}


Foam::label Foam::meshSearch::findNearestBoundaryFace
(
    const point& location,
    const label seedFacei,
    const bool useTreeSearch
) const
{
    if (seedFacei == -1)
    {
        if (useTreeSearch)
        {
            const indexedOctree<treeDataFace>& tree = boundaryTree();

            pointIndexHit info = boundaryTree().findNearest
            (
                location,
                magSqr(tree.bb().max()-tree.bb().min())
            );

            if (!info.hit())
            {
                info = boundaryTree().findNearest
                (
                    location,
                    Foam::sqr(great)
                );
            }

            return tree.shapes().faceLabels()[info.index()];
        }
        else
        {
            scalar minDist = great;

            label minFacei = -1;

            for
            (
                label facei = mesh_.nInternalFaces();
                facei < mesh_.nFaces();
                facei++
            )
            {
                const face& f =  mesh_.faces()[facei];

                pointHit curHit =
                    f.nearestPoint
                    (
                        location,
                        mesh_.points()
                    );

                if (curHit.distance() < minDist)
                {
                    minDist = curHit.distance();
                    minFacei = facei;
                }
            }
            return minFacei;
        }
    }
    else
    {
        return findNearestBoundaryFaceWalk(location, seedFacei);
    }
}


Foam::pointIndexHit Foam::meshSearch::intersection
(
    const point& pStart,
    const point& pEnd
) const
{
    pointIndexHit curHit = boundaryTree().findLine(pStart, pEnd);

    if (curHit.hit())
    {
        // Change index into octreeData into face label
        curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
    }
    return curHit;
}


Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
(
    const point& pStart,
    const point& pEnd
) const
{
    const indexedOctree<Foam::treeDataFace>& tree = boundaryTree();

    DynamicList<label> checkedFaces;
    findUniqueIntersectOp iop(tree, checkedFaces);

    DynamicList<pointIndexHit> hits;

    vector edgeVec = pEnd - pStart;
    edgeVec /= mag(edgeVec);

    point pt = pStart;

    pointIndexHit bHit;
    do
    {
        bHit = tree.findLine(pt, pEnd, iop);

        if (bHit.hit())
        {
            // Change index into octreeData into face label
            bHit.setIndex(tree.shapes().faceLabels()[bHit.index()]);

            hits.append(bHit);

            const vector& area = mesh_.faceAreas()[bHit.index()];

            scalar typDim = Foam::sqrt(mag(area));

            if ((mag(bHit.hitPoint() - pEnd)/typDim) < small)
            {
                break;
            }

            pt = bHit.hitPoint();
        }

    } while (bHit.hit());


    hits.shrink();

    return hits;
}


bool Foam::meshSearch::isInside(const point& p) const
{
    return (boundaryTree().getVolumeType(p) == volumeType::INSIDE);
}


void Foam::meshSearch::clearOut()
{
    boundaryTreePtr_.clear();
    cellTreePtr_.clear();
    overallBbPtr_.clear();
}


void Foam::meshSearch::correct()
{
    clearOut();
}


// ************************************************************************* //
meshSearch.C (20,510 bytes)   
meshSearch.H (8,416 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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/>.

Class
    Foam::meshSearch

Description
    Various (local, not parallel) searches on polyMesh;
    uses (demand driven) octree to search.

SourceFiles
    meshSearch.C

\*---------------------------------------------------------------------------*/

#ifndef meshSearch_H
#define meshSearch_H

#include "pointIndexHit.H"
#include "pointField.H"
#include "polyMesh.H"

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

namespace Foam
{

// Forward declaration of classes
class treeDataCell;
class treeDataFace;
template<class Type> class indexedOctree;
class treeBoundBox;

/*---------------------------------------------------------------------------*\
                           Class meshSearch Declaration
\*---------------------------------------------------------------------------*/

class meshSearch
{
    // Private data

        //- Reference to mesh
        const polyMesh& mesh_;

        //- Whether to use cell decomposition for all geometric tests
        const polyMesh::cellDecomposition cellDecompMode_;

        //- Data bounding box
        mutable autoPtr<treeBoundBox> overallBbPtr_;

        //- Demand driven octrees
        mutable autoPtr<indexedOctree<treeDataFace>> boundaryTreePtr_;
        mutable autoPtr<indexedOctree<treeDataCell>> cellTreePtr_;


    // Private Member Functions

        //- Updates nearestI, nearestDistSqr from any closer ones.
        static bool findNearer
        (
            const point& sample,
            const pointField& points,
            label& nearestI,
            scalar& nearestDistSqr
        );

        //- Updates nearestI, nearestDistSqr from any selected closer ones.
        static bool findNearer
        (
            const point& sample,
            const pointField& points,
            const labelList& indices,
            label& nearestI,
            scalar& nearestDistSqr
        );


        // Cells

            //- Nearest cell centre using octree
            label findNearestCellTree(const point&) const;

            //- Nearest cell centre going through all cells
            label findNearestCellLinear(const point&) const;

            //- Walk from seed. Does not 'go around' boundary, just returns
            //  last cell before boundary.
            label findNearestCellWalk(const point&, const label) const;

            //- Cell containing location. Linear search.
            label findCellLinear(const point&) const;

            //- Walk from seed. Does not 'go around' boundary, just returns
            //  last cell before boundary.
            label findCellWalk(const point&, const label) const;


        // Faces

            label findNearestFaceTree(const point&) const;

            label findNearestFaceLinear(const point&) const;

            label findNearestFaceWalk(const point&, const label) const;



        // Boundary faces

            //- Walk from seed to find nearest boundary face. Gets stuck in
            //  local minimum.
            label findNearestBoundaryFaceWalk
            (
                const point& location,
                const label seedFacei
            ) const;


        //- Disallow default bitwise copy construct
        meshSearch(const meshSearch&);

        //- Disallow default bitwise assignment
        void operator=(const meshSearch&);


public:

    // Declare name of the class and its debug switch
    ClassName("meshSearch");


    // Static data members

        //- Tolerance on linear dimensions
        static scalar tol_;


    // Constructors

        //- Construct from components. Constructs bb slightly bigger than
        //  mesh points bb.
        meshSearch
        (
            const polyMesh& mesh,
            const polyMesh::cellDecomposition = polyMesh::CELL_TETS
        );

        //- Construct with a custom bounding box. Any mesh element outside
        //  bb will not be found. Up to user to make sure bb
        //  extends slightly beyond wanted elements.
        meshSearch
        (
            const polyMesh& mesh,
            const treeBoundBox& bb,
            const polyMesh::cellDecomposition = polyMesh::CELL_TETS
        );

    //- Destructor
    ~meshSearch();


    // Member Functions

        // Access

            const polyMesh& mesh() const
            {
                return mesh_;
            }

            polyMesh::cellDecomposition decompMode() const
            {
                return cellDecompMode_;
            }

            //- Get (demand driven) reference to octree holding all
            //  boundary faces
            const indexedOctree<treeDataFace>& boundaryTree() const;

            //- Get (demand driven) reference to octree holding all cells
            const indexedOctree<treeDataCell>& cellTree() const;


        // Queries

            //- Find nearest cell in terms of cell centre.
            //  Options:
            //  - use octree
            //  - use linear search
            //  - if seed is provided walk. (uses findNearestCellWalk;
            //    does not handle holes in domain)
            label findNearestCell
            (
                const point& location,
                const label seedCelli = -1,
                const bool useTreeSearch = true
            ) const;

            label findNearestFace
            (
                const point& location,
                const label seedFacei = -1,
                const bool useTreeSearch = true
            ) const;

            //- Find cell containing location.
            //  If seed provided walks and falls back to linear/tree search.
            //  (so handles holes correctly)s
            //  Returns -1 if not in domain.
            label findCell
            (
                const point& location,
                const label seedCelli = -1,
                const bool useTreeSearch = true
            ) const;

            //- Find nearest boundary face
            //  If seed provided walks but then does not pass local minima
            //  in distance. Also does not jump from one connected region to
            //  the next.
            label findNearestBoundaryFace
            (
                const point& location,
                const label seedFacei = -1,
                const bool useTreeSearch = true
            ) const;

            //- Find first intersection of boundary in segment [pStart, pEnd]
            //  (so inclusive of endpoints). Always octree for now
            pointIndexHit intersection(const point& pStart, const point& pEnd)
            const;

            //- Find all intersections of boundary within segment pStart .. pEnd
            //  Always octree for now
            List<pointIndexHit> intersections
            (
                const point& pStart,
                const point& pEnd
            ) const;

            //- Determine inside/outside status
            bool isInside(const point&) const;


        //- Delete all storage
        void clearOut();

        //- Correct for mesh geom/topo changes
        void correct();
};


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

} // End namespace Foam

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

#endif

// ************************************************************************* //
meshSearch.H (8,416 bytes)   

will

2018-10-03 09:19

manager   ~0010081

Thanks for the fix, @MattijsJ. It needed a bit of modification, as marking all visited faces meant some similarly located hits (like the ones on the cyclics in the test case) were being missed. I changed it to consider only successfully hit faces. There's also an extra step in there to discount point-connected faces that the line projects through in the same direction. This gets rid of the duplicates when the line hits an edge or point.

This fix is now in dev as of 39ad49bc. This change, along with various modifications to Lagrangian since this bug was reported, I think this is now fixed.

Issue History

Date Modified Username Field Change
2014-02-01 19:47 user672 New Issue
2014-02-01 19:47 user672 File Added: sampleCyclic.tgz
2014-02-05 00:03 user672 Tag Attached: Added the used Allrun script
2014-02-05 00:03 user672 Tag Attached: baffle
2014-02-05 00:03 user672 Tag Attached: functionObject
2014-02-05 00:03 user672 Tag Attached: sample
2014-02-05 00:03 user672 Tag Detached: functionObject
2015-03-29 21:29 wyldckat File Added: bug1147.patch
2015-03-29 21:53 wyldckat Note Added: 0004549
2015-03-29 22:03 henry Note Added: 0004550
2018-07-10 11:16 administrator Tag Detached: Added the used Allrun script
2018-09-28 11:29 MattijsJ File Added: meshSearch.C
2018-09-28 11:29 MattijsJ File Added: meshSearch.H
2018-09-28 11:29 MattijsJ Note Added: 0010076
2018-10-03 09:19 will Assigned To => will
2018-10-03 09:19 will Status new => resolved
2018-10-03 09:19 will Resolution open => fixed
2018-10-03 09:19 will Fixed in Version => dev
2018-10-03 09:19 will Note Added: 0010081