View Issue Details

IDProjectCategoryView StatusLast Update
0001147OpenFOAM[All Projects] Bugpublic2018-10-03 09:19
Reporteruser672Assigned Towill 
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
PlatformLinuxOSScientific LinuxOS Version6.4
Product Version 
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)
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