View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0001506OpenFOAMBugpublic2015-02-01 19:282015-02-05 18:14
Reporterwyldckat Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilitysometimes
Status resolvedResolutionfixed 
Summary0001506: Problems with sampling unbased triangles
DescriptionThere are 2 issues that will be reported in this report, all related one way or another:

Issue 1: There is a line that is redefined without any particular reason in "src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C", namely this line:

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

This line is in this particular method:

    template<class Type>
    void Foam::isoSurfaceCell::generateTriPoints
    (
        const scalarField& cVals,
        const scalarField& pVals,

        const Field<Type>& cCoords,
        const Field<Type>& pCoords,

        const DynamicList<Type>& snappedPoints,
        const labelList& snappedCc,
        const labelList& snappedPoint,

        DynamicList<Type>& triPoints,
        DynamicList<label>& triMeshCells
    ) const


Issue 2: This same variable "cFaces" is used unchecked by this same method "generateTriPoints". For example, it's straight out used like this:

    const face& f0 = mesh_.faces()[cFaces[0]];

    //...
    
    label faceI = cFaces[cFaceI];
    const face& f = mesh_.faces()[faceI];

    const label fp0 = mesh_.tetBasePtIs()[faceI];
    
    
The code seems innocent enough, but there are a few situations where "cFaces[*]" is actually "-1". This can be proved by turning on this flag in "controlDict" (in the global/user one, not in the case's one):

    DebugSwitches
    {
        polyMesh 1;
    }

The following warning message will be issued when the aforementioned "cFaces[*]" have the value "-1":

    void polyMesh::initMesh() : initialising primitiveMesh
    --> FOAM Warning :
        From function const labelList& polyMesh::tetBasePtIs() const
        in file meshes/polyMesh/polyMesh.C at line 868
        Tet base point indices not available. Forcing storage of base points.

This error message will be issued by the utility "sample", because it has to request for the tetrahedral decomposition of the mesh (if I remember correctly).

For example, this can occur when this check is reported by "checkMesh -allTopology -allGeometry":

    ***Error in face tets: 2 faces with low quality or negative volume decomposition tets.
      <<Writing 2 faces with low quality or negative volume decomposition tets to set lowQualityTetFaces

This means that these 2 faces will have a base point index of "-1", because they were not decomposed with success.
Steps To ReproduceAttached is the tarball "isolated_problem.tar.gz", which has a subset mesh of a bigger case, including some results, for which the check mesh will give a lot of warnings and which has the 2 aforementioned damaged faces.

First, activate the global debug switch "polyMesh", so that you can see the warning.

After unpacking the case and going into the folder "isolated_problem", run:

  checkMesh -allTopology -allGeometry

  ./runSamples.sh

  
Assuming neither one of the 2 "sample" runs crash on Linux (although they will on Windows), then open the resulting VTK files that are in the resulting "postProcessing/surfaces/30/" folder and prepared to be amazed ;)
Attached is the image "isosurfaces_bug.png" that shows an example of this problem. The section cut is really near the problem faces and theoretically should also access the troubled cells.
Additional InformationThis issue was detected back on the 26th of October 2012, with OpenFOAM 2.1.x, but at the time it was diagnosed that this particular problem was something that would be deemed as "cannot reproduce on Linux", because we were only able to reproduce this on Windows. So we ended up postponing reporting it until we had time and/or a good reason to do so.

In the past few months, a couple of bug reports were spotted that seemed to hint at related errors to this one that we detected. Which is why only today the planets aligned and I managed to write up this bug report. Said bug reports are:
 - http://www.openfoam.org/mantisbt/view.php?id=1487
 - http://www.openfoam.org/mantisbt/view.php?id=1235 (Note: not entirely related, but sample-triangle referenced)
TagsNo tags attached.

Activities

wyldckat

2015-02-01 19:28

updater  

isolated_problem.tar.gz (24,999 bytes)

wyldckat

2015-02-01 19:28

updater  

isosurfaces_bug.png (29,005 bytes)   
isosurfaces_bug.png (29,005 bytes)   

wyldckat

2015-02-01 19:38

updater   ~0003643

Last edited: 2015-02-01 20:42

 Forgot to mention that in the image "isosurfaces_bug.png", that thing sticking out across diagonally is not meant to exist.


As for a solution for this problem, I'm not certain which one of the following is the best solution for OpenFOAM:

  a. That warning message that shows up with the debug switch "polyMesh", will actually cause an error and stopping when "sample" or any related functionality is used.

  b. Always check if every single "faceId" for this sampling mechanism will check if the value is sane (within [0, maxIndex]), or at least not smaller than 0. And if not sane, simply don't sample it and issue out a warning accordingly.

  c. Implemented a mechanism that can break apart damaged cells and still be able to create a barely working tetrahedral cell group or triangle face(s).


Either way, in a few minutes I'll post a proposition patch for "b", to be tested on issue 1487. It won't be efficient, but it can help diagnose if it's the problem for 1487...


*edit*: apparently issue #1487 is not directly related, since this case only demonstrates problems with sampling iso-surfaces.

wyldckat

2015-02-01 20:39

updater  

example_patch_for_isoSurfaceCellTemplates.patch (1,499 bytes)    
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index 4a88c3d..b592b51 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -303,6 +303,7 @@ void Foam::isoSurfaceCell::generateTriPoints
 ) const
 {
     tetMatcher tet;
+    label countNotFoundTets = 0;
 
     forAll(mesh_.cells(), cellI)
     {
@@ -404,7 +405,15 @@ void Foam::isoSurfaceCell::generateTriPoints
 
                     const label fp0 = mesh_.tetBasePtIs()[faceI];
 
+                    //skip undefined tetrahedral cells
+                    if(fp0 == -1)
+                    {
+                        countNotFoundTets++;
+                        continue;
+                    }
+
                     label fp = f.fcIndex(fp0);
+
                     for (label i = 2; i < f.size(); i++)
                     {
                         label nextFp = f.fcIndex(fp);
@@ -485,6 +494,15 @@ void Foam::isoSurfaceCell::generateTriPoints
             }
         }
     }
+    
+    if(countNotFoundTets>0)
+    {
+        WarningIn("Foam::isoSurfaceCell::generateTriPoints(...) const")
+            << "Did not find " << countNotFoundTets
+            << " tet base point indices, which has likelly lead to incomplete"
+            << " sample data. "
+            << endl;
+    }
 
     triPoints.shrink();
     triMeshCells.shrink();
example_patch_for_isoSurfaceCellTemplates.patch (1,499 bytes)   

wyldckat

2015-02-01 20:41

updater   ~0003644

Actually, I double checked just now and the problem is usually this line:

  const label fp0 = mesh_.tetBasePtIs()[faceI];

since it's "tetBasePtIs()" that not always has a clean list.


The attached patch "example_patch_for_isoSurfaceCellTemplates.patch" pretty much covered the issue for the attached example case, but there are several other occurrences of "tetBasePtIs()" which are not checked.

user4

2015-02-02 14:14

  ~0003651

The alternative is to fall back to using e.g. f[0] and generate inverted triangles.

henry

2015-02-02 14:36

manager   ~0003652

Which approach should be implemented to avoid problem?

wyldckat

2015-02-02 14:46

updater   ~0003653

If I understood Mattijs correctly, it could be a new "d" possible solution, something like this:

  const label fp0 = mesh_.tetBasePtIs()[faceI];

  if(fp0<0)
  {
    fp0 = 0; //better than nothing
    countNotFoundTets++;
  }


where "countNotFoundTets" is the variable introduced in the proposed patch.

Although I don't know if this isn't too risky, since picking the 0th vertex of the cell "when in doubt" might be too bad for getting an triangle at all.

If nobody else beats me to it, I'll try to test this possibility at the end of the day with the attached test case.

user4

2015-02-02 16:50

  ~0003654

Correct. I'd rather have an iso surface with some inverted triangles than a holey surface.

wyldckat

2015-02-02 20:02

updater  

example_patch_for_isoSurfaceCellTemplates_v2.patch (1,928 bytes)    
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index 4a88c3d..f234e5f 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -303,6 +303,7 @@ void Foam::isoSurfaceCell::generateTriPoints
 ) const
 {
     tetMatcher tet;
+    label countNotFoundTets = 0;
 
     forAll(mesh_.cells(), cellI)
     {
@@ -395,16 +396,25 @@ void Foam::isoSurfaceCell::generateTriPoints
             }
             else
             {
-                const cell& cFaces = mesh_.cells()[cellI];
+                //Repeated line
+                //const cell& cFaces = mesh_.cells()[cellI];
 
                 forAll(cFaces, cFaceI)
                 {
                     label faceI = cFaces[cFaceI];
                     const face& f = mesh_.faces()[faceI];
 
-                    const label fp0 = mesh_.tetBasePtIs()[faceI];
+                    label fp0 = mesh_.tetBasePtIs()[faceI];
+
+                    //skip undefined tetrahedral cells
+                    if(fp0 < 0)
+                    {
+                        fp0 = 0;
+                        countNotFoundTets++;
+                    }
 
                     label fp = f.fcIndex(fp0);
+
                     for (label i = 2; i < f.size(); i++)
                     {
                         label nextFp = f.fcIndex(fp);
@@ -485,6 +495,15 @@ void Foam::isoSurfaceCell::generateTriPoints
             }
         }
     }
+    
+    if(countNotFoundTets>0)
+    {
+        WarningIn("Foam::isoSurfaceCell::generateTriPoints(...) const")
+            << "Did not find " << countNotFoundTets
+            << " tet base point indices, which has likelly lead to inverted"
+            << " triangles. "
+            << endl;
+    }
 
     triPoints.shrink();
     triMeshCells.shrink();
example_patch_for_isoSurfaceCellTemplates_v2.patch (1,928 bytes)   

wyldckat

2015-02-02 20:02

updater  

comparison_v2_v1.png (109,198 bytes)   
comparison_v2_v1.png (109,198 bytes)   

wyldckat

2015-02-02 20:12

updater   ~0003655

Attached is the patch "example_patch_for_isoSurfaceCellTemplates_v2.patch", based on what Mattijs indicated.

In addition, the attached image "comparison_v2_v1.png" that compares the resulting 4 iso-surfaces resulting from the code patched with the latest proposed code patch versus the same iso-surfaces that were generated with the previous code patch.

And to complete the analysis, looks like this time it was a somewhat lucky triangle order, since the image "comparison_v2_v1_normals.png" shows the normal vectors coming out of the triangles.


Although, in this way of "fixing", a pre-filtered list such as "mesh_.tetBasePtIsBounded()" would likely be the most efficient way of fixing this... although not memory wise, since it would mean that yet another list would exist on RAM...

wyldckat

2015-02-02 20:13

updater  

comparison_v2_v1_normals.png (89,669 bytes)   
comparison_v2_v1_normals.png (89,669 bytes)   

henry

2015-02-05 18:14

manager   ~0003698

Resolved by commit 31b5ea0e3a530f865f1d00c9e81267bdf9a5c2d6

Issue History

Date Modified Username Field Change
2015-02-01 19:28 wyldckat New Issue
2015-02-01 19:28 wyldckat File Added: isolated_problem.tar.gz
2015-02-01 19:28 wyldckat File Added: isosurfaces_bug.png
2015-02-01 19:38 wyldckat Note Added: 0003643
2015-02-01 20:39 wyldckat File Added: example_patch_for_isoSurfaceCellTemplates.patch
2015-02-01 20:41 wyldckat Note Added: 0003644
2015-02-01 20:42 wyldckat Note Edited: 0003643
2015-02-02 14:14 user4 Note Added: 0003651
2015-02-02 14:36 henry Note Added: 0003652
2015-02-02 14:46 wyldckat Note Added: 0003653
2015-02-02 16:50 user4 Note Added: 0003654
2015-02-02 20:02 wyldckat File Added: example_patch_for_isoSurfaceCellTemplates_v2.patch
2015-02-02 20:02 wyldckat File Added: comparison_v2_v1.png
2015-02-02 20:12 wyldckat Note Added: 0003655
2015-02-02 20:13 wyldckat File Added: comparison_v2_v1_normals.png
2015-02-05 18:14 henry Note Added: 0003698
2015-02-05 18:14 henry Status new => resolved
2015-02-05 18:14 henry Resolution open => fixed
2015-02-05 18:14 henry Assigned To => henry