View Issue Details

IDProjectCategoryView StatusLast Update
0001487OpenFOAM[All Projects] Bugpublic2016-03-12 11:55
ReportermatteoLAssigned Tohenry 
PrioritynormalSeveritymajorReproducibilitysometimes
Status resolvedResolutionfixed 
PlatformunixOScentosOS Version(please specify)
Product Version 
Fixed in Version 
Summary0001487: Slice Sampling (interpolation)on ANSAgenerated mesh fails (glibc free) for large parallel meshes (>100 mln cells)
DescriptionHello,
I am using OF2.2 and I am having trouble running a sampling slice when working with a volume mesh generated with ANSA.

Here is my functionIncludeDict:

slices
{
type surfaces;
storeFilter false;
functionObjectLibs ( "libsampling.so" );
outputControl timeStep;
outputInterval 100;
setFormat raw;
surfaceFormat vtk;
interpolationScheme cellPoint;
fields U;
surfaces
(
X_Movie_0
{
type cuttingPlane;
planeType pointAndNormal;
pointAndNormalDict
{
normalVector (1.000000 0.000000 0.000000);
basePoint (-1.000000 -0.000000 -0.000000);
}
interpolate true;
mergeTol 1e-8;
}
);
}

When I run it (with execFlowFunctionObjects -latestTime ) I get the following error:
Surfaces slices: Reading surface description:
X_Movie_0

*** glibc detected *** execFlowFunctionObjects: double free or corruption (!prev): 0x00000000112848d0 ***
*** glibc detected *** execFlowFunctionObjects: free(): invalid next size (normal): 0x0000000010fe17e0 ***
*** glibc detected *** execFlowFunctionObjects: double free or corruption (!prev): 0x000000001087acd0 ***
*** glibc detected *** execFlowFunctionObjects: double free or corruption (!prev): 0x0000000013098d00 ***
*** glibc detected *** execFlowFunctionObjects: free(): invalid next size (normal): 0x00000000136d0880 ***
*** glibc detected *** execFlowFunctionObjects: double free or corruption (!prev): 0x0000000012675d60 ***
*** glibc detected *** execFlowFunctionObjects: double free or corruption (!prev): 0x0000000011f10520 ***
*** glibc detected *** execFlowFunctionObjects: malloc(): memory corruption: 0x0000000010bb2400 ***
*** glibc detected *** execFlowFunctionObjects: free(): invalid next size (normal): 0x00000000120b4550 ***
*** glibc detected *** execFlowFunctionObjects: free(): invalid next size (normal): 0x0000000013f0ea50 ***
======= Backtrace: =========
/lib64/libc.so.6(+0x750c6)[0x2b08c8f590c6]

Notice that:
1)The same dictionary works perfectly on a SnappyHexMesh generated mesh.
2)This dictionary works if I change the flag interpolate to false.
3)On a small test case, it works fine
4)On all my big mesh it fails

I have tried the other interpolation schemes options but those fail as well..
Steps To ReproduceIt Always happens on my big meshes. Unfortunately I cannot share a mesh with the problem because it is too big to upload and for confidentiality reasons.
Additional InformationI would be glad to give you any further information required or perform any test you need.
Thanks
TagsNo tags attached.

Activities

user4

2015-01-14 08:55

  ~0003537

- is your mesh ok?

    checkMesh -allGeometry

- does it work if you slightly perturb the basePoint - maybe you're cutting along edges or faces.

- set the mergeTol smaller

- switch on the debug flag for isoSurface

matteoL

2015-01-14 10:33

reporter  

StlFilesDebug.PNG (24,183 bytes)
StlFilesDebug.PNG (24,183 bytes)

matteoL

2015-01-14 10:41

reporter   ~0003543

Hello mattijs,
Thanks for your help.

here is the log of chekMesh:


Checking geometry...
    Overall domain bounding box (-28.796305 -42.500068 -1.027792) (54.959113 3.7873575e-05 20.510296)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (1.0127749e-15 -1.7754966e-16 -2.3397573e-15) OK.
    Max cell openness = 8.9573582e-16 OK.
    Max aspect ratio = 70.509793 OK.
    Minimum face area = 2.4086684e-10. Maximum face area = 1.7096958. Face area magnitudes OK.
    Min volume = 1.9372266e-14. Max volume = 2.057266. Total volume = 59463.08. Cell volumes OK.
    Mesh non-orthogonality Max: 85.118003 average: 16.411001
   *Number of severely non-orthogonal (> 70 degrees) faces: 278.
    Non-orthogonality check OK.
  <<Writing 278 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
 ***Max skewness = 5.105728, 24 highly skew faces detected which may impair the quality of the results
  <<Writing 24 skew faces to set skewFaces
    Coupled point location match (average 0) OK.
 ***Error in face tets: 1 faces with low quality or negative volume decomposition tets.
  <<Writing 1 faces with low quality or negative volume decomposition tets to set lowQualityTetFaces
   *Edges too small, min/max edge length = 3.0707654e-06 2.3819928, number too small: 2011612
  <<Writing 2011612 points on short edges to set shortEdges
  <<Writing 2600899 near (closer than 9.6359319e-05 apart) points to set nearPoints
   *There are 51431 faces with concave angles between consecutive edges. Max concave angle = 56.166339 degrees.
  <<Writing 51431 faces with concave angles to set concaveFaces
    Face flatness (1 = flat, 0 = butterfly) : average = 0.99977716 min = 0.59449512
   *There are 2 faces with ratio between projected and actual area < 0.8
    Minimum ratio (minimum flatness, maximum warpage) = 0.59449512
  <<Writing 2 warped faces to set warpedFaces
    Cell determinant (wellposedness) : minimum: 0 average: 3.9231018
 ***Cells with small determinant (< 0.001) found, number of cells: 50913
  <<Writing 50913 under-determined cells to set underdeterminedCells
 ***Concave cells (using face planes) found, number of cells: 2231003
  <<Writing 2231003 concave cells to set concaveCells

Failed 4 mesh checks.

I have reduced the mergetol to 1e-11, perturbed the basepoint and activated the debug flag. It is still crashing. I am uploading the log. With the debug it has created lots of stl files, most of them empty (I have uploaded a snapshot of my disk).

Any idea?
Thanks, again,
matteo

matteoL

2015-01-14 10:45

reporter  

WithDebug.txt.gz (37,034 bytes)

user4

2015-01-19 20:50

  ~0003567

From your traceback the problem is in the interpolation from volField to pointField (volPointInterpolation) and not related to the cutting plane. Do you see the same problem if e.g. you run foamToVTK?

- do you see the same problem in 2.3.x?
- what are your patch types?
- try running a small problem under valgrind (using e.g. the mpirunDebug script)

matteoL

2015-01-20 10:36

reporter   ~0003570

Hello,
foamToVTK works fine.

- yes.
- my patch types are:
wall, patch,cyclycAMI and processor
- unfortunately I can't run it on valgrind because I have those issues only on big cases. I haven't been able to generate a small case with such problem.

Indeed the problem is in the interpolation step, in fact if i run the sampling procedure without the interpolate option to TRUE everything works fine.

Thanks.

wyldckat

2015-02-21 20:52

updater   ~0003851

Greetings to all!

@Matteo: Any chance you can test this case with the latest bug fix version of OpenFOAM 2.3.x?
I ask this because with luck, this issue has been solved with the fix implemented for the bug report #1506: http://www.openfoam.org/mantisbt/view.php?id=1506

Best regards,
Bruno

user4

2015-03-11 12:21

  ~0004082

We can replicate your problem. It is due to a dangling pointer in the iso surface algorithm which is used by the cuttingPlane sampling method. We expect to have a fix for this soon.

As a workaround sometimes it helps to change the order of your functionObjects

matteoL

2015-08-20 15:09

reporter   ~0005277

Hello Mattijs,
Any news on this?

Thanks,
Matteo

wyldckat

2016-02-21 20:21

updater   ~0005976

This seems to have been fixed in the following commit on the OpenFOAM-history repository by Mattijs: https://github.com/OpenCFD/OpenFOAM-history/commit/a11a6cd044f4d5ce4ac34926be6f777bbc2bf5ec

@Matteo: Since you could not share the case with which this issue could be reproduced, can you please let us know if the commit I pointed out solves the problem?

If not, any chance you can provide a dummy test case that reproduces the same error?

MattijsJ

2016-03-11 17:14

reporter  

0001-BUG-isosurface-stores-reference-to-volPointInterpola.patch (67,596 bytes)
From 13d3ace2604b97361771b32b8d3e82f99cf0957f Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 11 Mar 2016 17:08:52 +0000
Subject: [PATCH] BUG: isosurface: stores reference to
 volPointInterpolation-cached field. If volPointInterpolation decides to
 remove field it invalidates the reference Resolves bug-report
 http://www.openfoam.org/mantisbt/view.php?id=1487

---
 .../sampledSurface/isoSurface/isoSurface.C         | 680 ++-------------------
 .../sampledSurface/isoSurface/isoSurface.H         |  80 +--
 .../sampledSurface/isoSurface/isoSurfaceCell.C     | 218 ++-----
 .../sampledSurface/isoSurface/isoSurfaceCell.H     |  45 --
 .../isoSurface/isoSurfaceCellTemplates.C           | 361 +++++------
 .../isoSurface/isoSurfaceTemplates.C               | 289 +++++----
 .../sampledSurface/isoSurface/sampledIsoSurface.C  |  72 ++-
 .../sampledSurface/isoSurface/sampledIsoSurface.H  |   1 -
 8 files changed, 470 insertions(+), 1276 deletions(-)

diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index a9d57ad..1ca392f 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -471,43 +471,6 @@ void Foam::isoSurface::calcCutTypes
 }
 
 
-// Return the two common points between two triangles
-Foam::labelPair Foam::isoSurface::findCommonPoints
-(
-    const labelledTri& tri0,
-    const labelledTri& tri1
-)
-{
-    labelPair common(-1, -1);
-
-    label fp0 = 0;
-    label fp1 = findIndex(tri1, tri0[fp0]);
-
-    if (fp1 == -1)
-    {
-        fp0 = 1;
-        fp1 = findIndex(tri1, tri0[fp0]);
-    }
-
-    if (fp1 != -1)
-    {
-        // So tri0[fp0] is tri1[fp1]
-
-        // Find next common point
-        label fp0p1 = tri0.fcIndex(fp0);
-        label fp1p1 = tri1.fcIndex(fp1);
-        label fp1m1 = tri1.rcIndex(fp1);
-
-        if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
-        {
-            common[0] = tri0[fp0];
-            common[1] = tri0[fp0p1];
-        }
-    }
-    return common;
-}
-
-
 // Caculate centre of surface.
 Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
 {
@@ -521,73 +484,6 @@ Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
 }
 
 
-// Replace surface (localPoints, localTris) with single point. Returns
-// point. Destructs arguments.
-Foam::pointIndexHit Foam::isoSurface::collapseSurface
-(
-    pointField& localPoints,
-    DynamicList<labelledTri, 64>& localTris
-)
-{
-    pointIndexHit info(false, vector::zero, localTris.size());
-
-    if (localTris.size() == 1)
-    {
-        const labelledTri& tri = localTris[0];
-        info.setPoint(tri.centre(localPoints));
-        info.setHit();
-    }
-    else if (localTris.size() == 2)
-    {
-        // Check if the two triangles share an edge.
-        const labelledTri& tri0 = localTris[0];
-        const labelledTri& tri1 = localTris[0];
-
-        labelPair shared = findCommonPoints(tri0, tri1);
-
-        if (shared[0] != -1)
-        {
-            info.setPoint
-            (
-                0.5
-              * (
-                    tri0.centre(localPoints)
-                  + tri1.centre(localPoints)
-                )
-            );
-            info.setHit();
-        }
-    }
-    else if (localTris.size())
-    {
-        // Check if single region. Rare situation.
-        triSurface surf
-        (
-            localTris,
-            geometricSurfacePatchList(0),
-            localPoints,
-            true
-        );
-        localTris.clearStorage();
-
-        labelList faceZone;
-        label nZones = surf.markZones
-        (
-            boolList(surf.nEdges(), false),
-            faceZone
-        );
-
-        if (nZones == 1)
-        {
-            info.setPoint(calcCentre(surf));
-            info.setHit();
-        }
-    }
-
-    return info;
-}
-
-
 // Determine per cell centre whether all the intersections get collapsed
 // to a single point
 void Foam::isoSurface::calcSnappedCc
@@ -1202,435 +1098,6 @@ bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
 }
 
 
-void Foam::isoSurface::calcAddressing
-(
-    const triSurface& surf,
-    List<FixedList<label, 3>>& faceEdges,
-    labelList& edgeFace0,
-    labelList& edgeFace1,
-    Map<labelList>& edgeFacesRest
-) const
-{
-    const pointField& points = surf.points();
-
-    pointField edgeCentres(3*surf.size());
-    label edgeI = 0;
-    forAll(surf, triI)
-    {
-        const labelledTri& tri = surf[triI];
-        edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
-        edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
-        edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
-    }
-
-    pointField mergedCentres;
-    labelList oldToMerged;
-    bool hasMerged = mergePoints
-    (
-        edgeCentres,
-        mergeDistance_,
-        false,
-        oldToMerged,
-        mergedCentres
-    );
-
-    if (debug)
-    {
-        Pout<< "isoSurface : detected "
-            << mergedCentres.size()
-            << " geometric edges on " << surf.size() << " triangles." << endl;
-    }
-
-    if (!hasMerged)
-    {
-        return;
-    }
-
-
-    // Determine faceEdges
-    faceEdges.setSize(surf.size());
-    edgeI = 0;
-    forAll(surf, triI)
-    {
-        faceEdges[triI][0] = oldToMerged[edgeI++];
-        faceEdges[triI][1] = oldToMerged[edgeI++];
-        faceEdges[triI][2] = oldToMerged[edgeI++];
-    }
-
-
-    // Determine edgeFaces
-    edgeFace0.setSize(mergedCentres.size());
-    edgeFace0 = -1;
-    edgeFace1.setSize(mergedCentres.size());
-    edgeFace1 = -1;
-    edgeFacesRest.clear();
-
-    // Overflow edge faces for geometric shared edges that turned
-    // out to be different anyway.
-    EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
-
-    forAll(oldToMerged, oldEdgeI)
-    {
-        label triI = oldEdgeI / 3;
-        label edgeI = oldToMerged[oldEdgeI];
-
-        if (edgeFace0[edgeI] == -1)
-        {
-            // First triangle for edge
-            edgeFace0[edgeI] = triI;
-        }
-        else
-        {
-            //- Check that the two triangles actually topologically
-            //  share an edge
-            const labelledTri& prevTri = surf[edgeFace0[edgeI]];
-            const labelledTri& tri = surf[triI];
-
-            label fp = oldEdgeI % 3;
-
-            edge e(tri[fp], tri[tri.fcIndex(fp)]);
-
-            label prevTriIndex = -1;
-
-            forAll(prevTri, i)
-            {
-                if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
-                {
-                    prevTriIndex = i;
-                    break;
-                }
-            }
-
-            if (prevTriIndex == -1)
-            {
-                // Different edge. Store for later.
-                EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
-
-                if (iter != extraEdgeFaces.end())
-                {
-                    labelList& eFaces = iter();
-                    label sz = eFaces.size();
-                    eFaces.setSize(sz+1);
-                    eFaces[sz] = triI;
-                }
-                else
-                {
-                    extraEdgeFaces.insert(e, labelList(1, triI));
-                }
-            }
-            else if (edgeFace1[edgeI] == -1)
-            {
-                edgeFace1[edgeI] = triI;
-            }
-            else
-            {
-                //WarningInFunction
-                //    << "Edge " << edgeI << " with centre "
-                //    << mergedCentres[edgeI]
-                //    << " used by more than two triangles: "
-                //    << edgeFace0[edgeI] << ", "
-                //    << edgeFace1[edgeI] << " and " << triI << endl;
-                Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
-
-                if (iter != edgeFacesRest.end())
-                {
-                    labelList& eFaces = iter();
-                    label sz = eFaces.size();
-                    eFaces.setSize(sz+1);
-                    eFaces[sz] = triI;
-                }
-                else
-                {
-                    edgeFacesRest.insert(edgeI, labelList(1, triI));
-                }
-            }
-        }
-    }
-
-    // Add extraEdgeFaces
-    edgeI = edgeFace0.size();
-
-    edgeFace0.setSize(edgeI + extraEdgeFaces.size());
-    edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
-
-    forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
-    {
-        const labelList& eFaces = iter();
-
-        // The current edge will become edgeI. Replace all occurrences in
-        // faceEdges
-        forAll(eFaces, i)
-        {
-            label triI = eFaces[i];
-            const labelledTri& tri = surf[triI];
-
-            FixedList<label, 3>& fEdges = faceEdges[triI];
-            forAll(tri, fp)
-            {
-                edge e(tri[fp], tri[tri.fcIndex(fp)]);
-
-                if (e == iter.key())
-                {
-                    fEdges[fp] = edgeI;
-                    break;
-                }
-            }
-        }
-
-
-        // Add face to edgeFaces
-
-        edgeFace0[edgeI] = eFaces[0];
-
-        if (eFaces.size() >= 2)
-        {
-            edgeFace1[edgeI] = eFaces[1];
-
-            if (eFaces.size() > 2)
-            {
-                edgeFacesRest.insert
-                (
-                    edgeI,
-                    SubList<label>(eFaces, eFaces.size()-2, 2)
-                );
-            }
-        }
-
-        edgeI++;
-    }
-}
-
-
-void Foam::isoSurface::walkOrientation
-(
-    const triSurface& surf,
-    const List<FixedList<label, 3>>& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const label seedTriI,
-    labelList& flipState
-)
-{
-    // Do walk for consistent orientation.
-    DynamicList<label> changedFaces(surf.size());
-
-    changedFaces.append(seedTriI);
-
-    while (changedFaces.size())
-    {
-        DynamicList<label> newChangedFaces(changedFaces.size());
-
-        forAll(changedFaces, i)
-        {
-            label triI = changedFaces[i];
-            const labelledTri& tri = surf[triI];
-            const FixedList<label, 3>& fEdges = faceEdges[triI];
-
-            forAll(fEdges, fp)
-            {
-                label edgeI = fEdges[fp];
-
-                // my points:
-                label p0 = tri[fp];
-                label p1 = tri[tri.fcIndex(fp)];
-
-                label nbrI =
-                (
-                    edgeFace0[edgeI] != triI
-                  ? edgeFace0[edgeI]
-                  : edgeFace1[edgeI]
-                );
-
-                if (nbrI != -1 && flipState[nbrI] == -1)
-                {
-                    const labelledTri& nbrTri = surf[nbrI];
-
-                    // nbr points
-                    label nbrFp = findIndex(nbrTri, p0);
-
-                    if (nbrFp == -1)
-                    {
-                        FatalErrorInFunction
-                            << "triI:" << triI
-                            << " tri:" << tri
-                            << " p0:" << p0
-                            << " p1:" << p1
-                            << " fEdges:" << fEdges
-                            << " edgeI:" << edgeI
-                            << " edgeFace0:" << edgeFace0[edgeI]
-                            << " edgeFace1:" << edgeFace1[edgeI]
-                            << " nbrI:" << nbrI
-                            << " nbrTri:" << nbrTri
-                            << abort(FatalError);
-                    }
-
-
-                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-
-                    bool sameOrientation = (p1 == nbrP1);
-
-                    if (flipState[triI] == 0)
-                    {
-                        flipState[nbrI] = (sameOrientation ? 0 : 1);
-                    }
-                    else
-                    {
-                        flipState[nbrI] = (sameOrientation ? 1 : 0);
-                    }
-                    newChangedFaces.append(nbrI);
-                }
-            }
-        }
-
-        changedFaces.transfer(newChangedFaces);
-    }
-}
-
-
-void Foam::isoSurface::orientSurface
-(
-    triSurface& surf,
-    const List<FixedList<label, 3>>& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const Map<labelList>& edgeFacesRest
-)
-{
-    // -1 : unvisited
-    //  0 : leave as is
-    //  1 : flip
-    labelList flipState(surf.size(), -1);
-
-    label seedTriI = 0;
-
-    while (true)
-    {
-        // Find first unvisited triangle
-        for
-        (
-            ;
-            seedTriI < surf.size() && flipState[seedTriI] != -1;
-            seedTriI++
-        )
-        {}
-
-        if (seedTriI == surf.size())
-        {
-            break;
-        }
-
-        // Note: Determine orientation of seedTriI?
-        // for now assume it is ok
-        flipState[seedTriI] = 0;
-
-        walkOrientation
-        (
-            surf,
-            faceEdges,
-            edgeFace0,
-            edgeFace1,
-            seedTriI,
-            flipState
-        );
-    }
-
-    // Do actual flipping
-    surf.clearOut();
-    forAll(surf, triI)
-    {
-        if (flipState[triI] == 1)
-        {
-            labelledTri tri(surf[triI]);
-
-            surf[triI][0] = tri[0];
-            surf[triI][1] = tri[2];
-            surf[triI][2] = tri[1];
-        }
-        else if (flipState[triI] == -1)
-        {
-            FatalErrorInFunction
-                << "problem" << abort(FatalError);
-        }
-    }
-}
-
-
-// Checks if triangle is connected through edgeI only.
-bool Foam::isoSurface::danglingTriangle
-(
-    const FixedList<label, 3>& fEdges,
-    const labelList& edgeFace1
-)
-{
-    label nOpen = 0;
-    forAll(fEdges, i)
-    {
-        if (edgeFace1[fEdges[i]] == -1)
-        {
-            nOpen++;
-        }
-    }
-
-    if (nOpen == 1 || nOpen == 2 || nOpen == 3)
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// Mark triangles to keep. Returns number of dangling triangles.
-Foam::label Foam::isoSurface::markDanglingTriangles
-(
-    const List<FixedList<label, 3>>& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const Map<labelList>& edgeFacesRest,
-    boolList& keepTriangles
-)
-{
-    keepTriangles.setSize(faceEdges.size());
-    keepTriangles = true;
-
-    label nDangling = 0;
-
-    // Remove any dangling triangles
-    forAllConstIter(Map<labelList>,  edgeFacesRest, iter)
-    {
-        // These are all the non-manifold edges. Filter out all triangles
-        // with only one connected edge (= this edge)
-
-        label edgeI = iter.key();
-        const labelList& otherEdgeFaces = iter();
-
-        // Remove all dangling triangles
-        if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
-        {
-            keepTriangles[edgeFace0[edgeI]] = false;
-            nDangling++;
-        }
-        if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
-        {
-            keepTriangles[edgeFace1[edgeI]] = false;
-            nDangling++;
-        }
-        forAll(otherEdgeFaces, i)
-        {
-            label triI = otherEdgeFaces[i];
-            if (danglingTriangle(faceEdges[triI], edgeFace1))
-            {
-                keepTriangles[triI] = false;
-                nDangling++;
-            }
-        }
-    }
-    return nDangling;
-}
-
-
 Foam::triSurface Foam::isoSurface::subsetMesh
 (
     const triSurface& s,
@@ -1953,57 +1420,58 @@ Foam::isoSurface::isoSurface
     }
 
 
+    {
+        DynamicList<point> triPoints(3*nCutCells_);
+        DynamicList<label> triMeshCells(nCutCells_);
 
-    DynamicList<point> triPoints(nCutCells_);
-    DynamicList<label> triMeshCells(nCutCells_);
-
-    generateTriPoints
-    (
-        cValsPtr_(),
-        pVals_,
+        generateTriPoints
+        (
+            cValsPtr_(),
+            pVals_,
 
-        meshC,
-        mesh_.points(),
+            meshC,
+            mesh_.points(),
 
-        snappedPoints,
-        snappedCc,
-        snappedPoint,
+            snappedPoints,
+            snappedCc,
+            snappedPoint,
 
-        triPoints,
-        triMeshCells
-    );
+            triPoints,      // 3 points of the triangle
+            triMeshCells    // per triangle the originating cell
+        );
 
-    if (debug)
-    {
-        Pout<< "isoSurface : generated " << triMeshCells.size()
-            << " unmerged triangles from " << triPoints.size()
-            << " unmerged points." << endl;
-    }
+        if (debug)
+        {
+            Pout<< "isoSurface : generated " << triMeshCells.size()
+                << " unmerged triangles from " << triPoints.size()
+                << " unmerged points." << endl;
+        }
 
 
-    // Merge points and compact out non-valid triangles
-    labelList triMap;           // merged to unmerged triangle
-    triSurface::operator=
-    (
-        stitchTriPoints
+        // Merge points and compact out non-valid triangles
+        labelList triMap;           // merged to unmerged triangle
+        triSurface::operator=
         (
-            true,               // check for duplicate tris
-            triPoints,
-            triPointMergeMap_,  // unmerged to merged point
-            triMap
-        )
-    );
+            stitchTriPoints
+            (
+                true,               // check for duplicate tris
+                triPoints,
+                triPointMergeMap_,  // unmerged to merged point
+                triMap
+            )
+        );
 
-    if (debug)
-    {
-        Pout<< "isoSurface : generated " << triMap.size()
-            << " merged triangles." << endl;
-    }
+        if (debug)
+        {
+            Pout<< "isoSurface : generated " << triMap.size()
+                << " merged triangles." << endl;
+        }
 
-    meshCells_.setSize(triMap.size());
-    forAll(triMap, i)
-    {
-        meshCells_[i] = triMeshCells[triMap[i]];
+        meshCells_.setSize(triMap.size());
+        forAll(triMap, i)
+        {
+            meshCells_[i] = triMeshCells[triMap[i]];
+        }
     }
 
     if (debug)
@@ -2019,70 +1487,6 @@ Foam::isoSurface::isoSurface
     }
 
 
-    if (false)
-    {
-        List<FixedList<label, 3>> faceEdges;
-        labelList edgeFace0, edgeFace1;
-        Map<labelList> edgeFacesRest;
-
-
-        while (true)
-        {
-            // Calculate addressing
-            calcAddressing
-            (
-                *this,
-                faceEdges,
-                edgeFace0,
-                edgeFace1,
-                edgeFacesRest
-            );
-
-            // See if any dangling triangles
-            boolList keepTriangles;
-            label nDangling = markDanglingTriangles
-            (
-                faceEdges,
-                edgeFace0,
-                edgeFace1,
-                edgeFacesRest,
-                keepTriangles
-            );
-
-            if (debug)
-            {
-                Pout<< "isoSurface : detected " << nDangling
-                    << " dangling triangles." << endl;
-            }
-
-            if (nDangling == 0)
-            {
-                break;
-            }
-
-            // Create face map (new to old)
-            labelList subsetTriMap(findIndices(keepTriangles, true));
-
-            labelList subsetPointMap;
-            labelList reversePointMap;
-            triSurface::operator=
-            (
-                subsetMesh
-                (
-                    *this,
-                    subsetTriMap,
-                    reversePointMap,
-                    subsetPointMap
-                )
-            );
-            meshCells_ = labelField(meshCells_, subsetTriMap);
-            inplaceRenumber(reversePointMap, triPointMergeMap_);
-        }
-
-        orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-    }
-
-
     if (debug)
     {
         fileName stlFile = mesh_.time().path() + ".stl";
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 967536b..6937bfb 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -193,20 +193,8 @@ class isoSurface
             const scalarField& pVals
         );
 
-        static labelPair findCommonPoints
-        (
-            const labelledTri&,
-            const labelledTri&
-        );
-
         static point calcCentre(const triSurface&);
 
-        static pointIndexHit collapseSurface
-        (
-            pointField& localPoints,
-            DynamicList<labelledTri, 64>& localTris
-        );
-
         //- Determine per cc whether all near cuts can be snapped to single
         //  point.
         void calcSnappedCc
@@ -323,6 +311,14 @@ class isoSurface
             DynamicList<label>& triMeshCells
         ) const;
 
+        template<class Type>
+        static tmp<Field<Type>> interpolate
+        (
+            const label nPoints,
+            const labelList& triPointMergeMap,
+            const DynamicList<Type>& unmergedValues
+        );
+
         triSurface stitchTriPoints
         (
             const bool checkDuplicates,
@@ -334,54 +330,6 @@ class isoSurface
         //- Check single triangle for (topological) validity
         static bool validTri(const triSurface&, const label);
 
-        //- Determine edge-face addressing
-        void calcAddressing
-        (
-            const triSurface& surf,
-            List<FixedList<label, 3>>& faceEdges,
-            labelList& edgeFace0,
-            labelList& edgeFace1,
-            Map<labelList>& edgeFacesRest
-        ) const;
-
-        //- Determine orientation
-        static void walkOrientation
-        (
-            const triSurface& surf,
-            const List<FixedList<label, 3>>& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const label seedTriI,
-            labelList& flipState
-        );
-
-        //- Orient surface
-        static void orientSurface
-        (
-            triSurface&,
-            const List<FixedList<label, 3>>& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const Map<labelList>& edgeFacesRest
-        );
-
-        //- Is triangle (given by 3 edges) not fully connected?
-        static bool danglingTriangle
-        (
-            const FixedList<label, 3>& fEdges,
-            const labelList& edgeFace1
-        );
-
-        //- Mark all non-fully connected triangles
-        static label markDanglingTriangles
-        (
-            const List<FixedList<label, 3>>& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const Map<labelList>& edgeFacesRest,
-            boolList& keepTriangles
-        );
-
         static triSurface subsetMesh
         (
             const triSurface& s,
@@ -392,6 +340,10 @@ class isoSurface
 
 public:
 
+    //- Declare friendship with isoSurfaceCell to share some functionality
+    friend class isoSurfaceCell;
+
+
     //- Runtime type information
     TypeName("isoSurface");
 
@@ -413,18 +365,12 @@ public:
 
     // Member Functions
 
-        //- For every face original cell in mesh
+        //- For every triangle the original cell in mesh
         const labelList& meshCells() const
         {
             return meshCells_;
         }
 
-        //- For every unmerged triangle point the point in the triSurface
-        const labelList& triPointMergeMap() const
-        {
-            return triPointMergeMap_;
-        }
-
         //- Interpolates cCoords,pCoords. Uses the references to the original
         //  fields used to create the iso surface.
         template<class Type>
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
index c2fec19..21359a8 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
@@ -1222,142 +1222,6 @@ void Foam::isoSurfaceCell::calcAddressing
 }
 
 
-//void Foam::isoSurfaceCell::walkOrientation
-//(
-//    const triSurface& surf,
-//    const List<FixedList<label, 3>>& faceEdges,
-//    const labelList& edgeFace0,
-//    const labelList& edgeFace1,
-//    const label seedTriI,
-//    labelList& flipState
-//)
-//{
-//    // Do walk for consistent orientation.
-//    DynamicList<label> changedFaces(surf.size());
-//
-//    changedFaces.append(seedTriI);
-//
-//    while (changedFaces.size())
-//    {
-//        DynamicList<label> newChangedFaces(changedFaces.size());
-//
-//        forAll(changedFaces, i)
-//        {
-//            label triI = changedFaces[i];
-//            const labelledTri& tri = surf[triI];
-//            const FixedList<label, 3>& fEdges = faceEdges[triI];
-//
-//            forAll(fEdges, fp)
-//            {
-//                label edgeI = fEdges[fp];
-//
-//                // my points:
-//                label p0 = tri[fp];
-//                label p1 = tri[tri.fcIndex(fp)];
-//
-//                label nbrI =
-//                (
-//                    edgeFace0[edgeI] != triI
-//                  ? edgeFace0[edgeI]
-//                  : edgeFace1[edgeI]
-//                );
-//
-//                if (nbrI != -1 && flipState[nbrI] == -1)
-//                {
-//                    const labelledTri& nbrTri = surf[nbrI];
-//
-//                    // nbr points
-//                    label nbrFp = findIndex(nbrTri, p0);
-//                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-//
-//                    bool sameOrientation = (p1 == nbrP1);
-//
-//                    if (flipState[triI] == 0)
-//                    {
-//                        flipState[nbrI] = (sameOrientation ? 0 : 1);
-//                    }
-//                    else
-//                    {
-//                        flipState[nbrI] = (sameOrientation ? 1 : 0);
-//                    }
-//                    newChangedFaces.append(nbrI);
-//                }
-//            }
-//        }
-//
-//        changedFaces.transfer(newChangedFaces);
-//    }
-//}
-//
-//
-//void Foam::isoSurfaceCell::orientSurface
-//(
-//    triSurface& surf,
-//    const List<FixedList<label, 3>>& faceEdges,
-//    const labelList& edgeFace0,
-//    const labelList& edgeFace1,
-//    const Map<labelList>& edgeFacesRest
-//)
-//{
-//    // -1 : unvisited
-//    //  0 : leave as is
-//    //  1 : flip
-//    labelList flipState(surf.size(), -1);
-//
-//    label seedTriI = 0;
-//
-//    while (true)
-//    {
-//        // Find first unvisited triangle
-//        for
-//        (
-//            ;
-//            seedTriI < surf.size() && flipState[seedTriI] != -1;
-//            seedTriI++
-//        )
-//        {}
-//
-//        if (seedTriI == surf.size())
-//        {
-//            break;
-//        }
-//
-//        // Note: Determine orientation of seedTriI?
-//        // for now assume it is ok
-//        flipState[seedTriI] = 0;
-//
-//        walkOrientation
-//        (
-//            surf,
-//            faceEdges,
-//            edgeFace0,
-//            edgeFace1,
-//            seedTriI,
-//            flipState
-//        );
-//    }
-//
-//    // Do actual flipping
-//    surf.clearOut();
-//    forAll(surf, triI)
-//    {
-//        if (flipState[triI] == 1)
-//        {
-//            labelledTri tri(surf[triI]);
-//
-//            surf[triI][0] = tri[0];
-//            surf[triI][1] = tri[2];
-//            surf[triI][2] = tri[1];
-//        }
-//        else if (flipState[triI] == -1)
-//        {
-//            FatalErrorInFunction
-//                << "problem" << abort(FatalError);
-//        }
-//    }
-//}
-
-
 // Checks if triangle is connected through edgeI only.
 bool Foam::isoSurfaceCell::danglingTriangle
 (
@@ -1605,57 +1469,59 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 
 
+    {
+        DynamicList<point> triPoints(nCutCells_);
+        DynamicList<label> triMeshCells(nCutCells_);
 
-    DynamicList<point> triPoints(nCutCells_);
-    DynamicList<label> triMeshCells(nCutCells_);
-
-    generateTriPoints
-    (
-        cVals,
-        pVals,
+        generateTriPoints
+        (
+            cVals,
+            pVals,
 
-        mesh_.cellCentres(),
-        mesh_.points(),
+            mesh_.cellCentres(),
+            mesh_.points(),
 
-        snappedPoints,
-        snappedCc,
-        snappedPoint,
+            snappedPoints,
+            snappedCc,
+            snappedPoint,
 
-        triPoints,
-        triMeshCells
-    );
+            triPoints,
+            triMeshCells
+        );
 
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
-            << " unmerged triangles." << endl;
-    }
+        if (debug)
+        {
+            Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
+                << " unmerged triangles." << endl;
+        }
 
-    // Merge points and compact out non-valid triangles
-    labelList triMap;           // merged to unmerged triangle
-    triSurface::operator=
-    (
-        stitchTriPoints
+        // Merge points and compact out non-valid triangles
+        labelList triMap;
+        triSurface::operator=
         (
-            regularise,         // check for duplicate tris
-            triPoints,
-            triPointMergeMap_,  // unmerged to merged point
-            triMap
-        )
-    );
+            stitchTriPoints
+            (
+                regularise,         // check for duplicate tris
+                triPoints,
+                triPointMergeMap_,  // unmerged to merged point
+                triMap              // merged to unmerged triangle
+            )
+        );
 
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : generated " << triMap.size()
-            << " merged triangles." << endl;
-    }
+        if (debug)
+        {
+            Pout<< "isoSurfaceCell : generated " << triMap.size()
+                << " merged triangles." << endl;
+        }
 
-    meshCells_.setSize(triMap.size());
-    forAll(triMap, i)
-    {
-        meshCells_[i] = triMeshCells[triMap[i]];
+        meshCells_.setSize(triMap.size());
+        forAll(triMap, i)
+        {
+            meshCells_[i] = triMeshCells[triMap[i]];
+        }
     }
 
+
     if (debug)
     {
         Pout<< "isoSurfaceCell : checking " << size()
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
index 034fbd3..21b54f2 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
@@ -214,19 +214,15 @@ class isoSurfaceCell
         void generateTriPoints
         (
             const DynamicList<Type>& snapped,
-            const scalar isoVal0,
             const scalar s0,
             const Type& p0,
             const label p0Index,
-            const scalar isoVal1,
             const scalar s1,
             const Type& p1,
             const label p1Index,
-            const scalar isoVal2,
             const scalar s2,
             const Type& p2,
             const label p2Index,
-            const scalar isoVal3,
             const scalar s3,
             const Type& p3,
             const label p3Index,
@@ -271,27 +267,6 @@ class isoSurfaceCell
             Map<labelList>& edgeFacesRest
         ) const;
 
-        ////- Determine orientation
-        //static void walkOrientation
-        //(
-        //    const triSurface& surf,
-        //    const List<FixedList<label, 3>>& faceEdges,
-        //    const labelList& edgeFace0,
-        //    const labelList& edgeFace1,
-        //    const label seedTriI,
-        //    labelList& flipState
-        //);
-
-        ////- Orient surface
-        //static void orientSurface
-        //(
-        //    triSurface&,
-        //    const List<FixedList<label, 3>>& faceEdges,
-        //    const labelList& edgeFace0,
-        //    const labelList& edgeFace1,
-        //    const Map<labelList>& edgeFacesRest
-        //);
-
         //- Is triangle (given by 3 edges) not fully connected?
         static bool danglingTriangle
         (
@@ -317,8 +292,6 @@ class isoSurfaceCell
             labelList& newToOldPoints
         );
 
-        //- Combine all triangles inside a cell into a minimal triangulation
-        void combineCellTriangles();
 
 public:
 
@@ -348,24 +321,6 @@ public:
             return meshCells_;
         }
 
-        //- For every unmerged triangle point the point in the triSurface
-        const labelList triPointMergeMap() const
-        {
-            return triPointMergeMap_;
-        }
-
-
-        //- Interpolates cCoords,pCoords. Takes the original fields
-        //  used to create the iso surface.
-        template<class Type>
-        tmp<Field<Type>> interpolate
-        (
-            const scalarField& cVals,
-            const scalarField& pVals,
-            const Field<Type>& cCoords,
-            const Field<Type>& pCoords
-        ) const;
-
         //- Interpolates cCoords,pCoords.
         template<class Type>
         tmp<Field<Type>> interpolate
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index 1eba206..ea26461 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -26,6 +26,7 @@ License
 #include "isoSurfaceCell.H"
 #include "polyMesh.H"
 #include "tetMatcher.H"
+#include "isoSurface.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -76,22 +77,18 @@ void Foam::isoSurfaceCell::generateTriPoints
 (
     const DynamicList<Type>& snapped,
 
-    const scalar isoVal0,
     const scalar s0,
     const Type& p0,
     const label p0Index,
 
-    const scalar isoVal1,
     const scalar s1,
     const Type& p1,
     const label p1Index,
 
-    const scalar isoVal2,
     const scalar s2,
     const Type& p2,
     const label p2Index,
 
-    const scalar isoVal3,
     const scalar s3,
     const Type& p3,
     const label p3Index,
@@ -117,167 +114,203 @@ void Foam::isoSurfaceCell::generateTriPoints
         triIndex |= 8;
     }
 
-    // Form the vertices of the triangles for each case
+    /* Form the vertices of the triangles for each case */
     switch (triIndex)
     {
         case 0x00:
         case 0x0F:
         break;
 
-        case 0x0E:
         case 0x01:
+        case 0x0E:
         {
-            // 0 is common point. Orient such that normal points in positive
-            // gradient direction
-            if (isoVal0 >= isoVal1)
-            {
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            }
-            else
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+            );
+            if (triIndex == 0x0E)
             {
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x0D:
         case 0x02:
+        case 0x0D:
         {
-            // 1 is common point
-            if (isoVal1 >= isoVal0)
-            {
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            }
-            else
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+            );
+
+            if (triIndex == 0x0D)
             {
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x0C:
         case 0x03:
+        case 0x0C:
         {
-            Type s02 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
-            Type s13 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
-
-            if (isoVal0 >= isoVal3)
-            {
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(s02);
-                pts.append(s13);
-                pts.append(s13);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s02);
-            }
-            else
+            Type p0p2 =
+                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type p1p3 =
+                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+            );
+            pts.append(p1p3);
+            pts.append(p0p2);
+
+            pts.append(p1p3);
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+            );
+            pts.append(p0p2);
+
+            if (triIndex == 0x0C)
             {
-                pts.append(s02);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(s13);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s13);
-                pts.append(s02);
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-5], pts[sz-4]);
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x0B:
         case 0x04:
+        case 0x0B:
         {
-            // 2 is common point
-            if (isoVal2 >= isoVal0)
+            pts.append
+            (
+                generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index)
+            );
+
+            if (triIndex == 0x0B)
             {
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x0A:
         case 0x05:
+        case 0x0A:
         {
-            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
-            if (isoVal3 >= isoVal0)
-            {
-                pts.append(s01);
-                pts.append(s23);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(s01);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s23);
-            }
-            else
+            Type p0p1 =
+                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type p2p3 =
+                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            pts.append(p0p1);
+            pts.append(p2p3);
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+            );
+
+            pts.append(p0p1);
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+            );
+            pts.append(p2p3);
+
+            if (triIndex == 0x0A)
             {
-                pts.append(s23);
-                pts.append(s01);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s01);
-                pts.append(s23);
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-5], pts[sz-4]);
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x09:
         case 0x06:
+        case 0x09:
         {
-            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
-            if (isoVal3 >= isoVal1)
-            {
-                pts.append(s01);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(s23);
-                pts.append(s01);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(s23);
-            }
-            else
+            Type p0p1 =
+                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type p2p3 =
+                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            pts.append(p0p1);
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+            );
+            pts.append(p2p3);
+
+            pts.append(p0p1);
+            pts.append(p2p3);
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+            );
+
+            if (triIndex == 0x09)
             {
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(s01);
-                pts.append(s23);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(s01);
-                pts.append(s23);
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-5], pts[sz-4]);
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
 
-        case 0x07:
         case 0x08:
+        case 0x07:
         {
-            // 3 is common point
-            if (isoVal3 >= isoVal0)
+            pts.append
+            (
+                generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index)
+            );
+            if (triIndex == 0x07)
             {
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+                // Flip normals
+                label sz = pts.size();
+                Swap(pts[sz-2], pts[sz-1]);
             }
         }
         break;
@@ -341,22 +374,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                     (
                         snappedPoints,
 
-                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
-                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
-                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
-                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -370,22 +399,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                     (
                         snappedPoints,
 
-                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
-                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
-                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
-                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -411,7 +436,6 @@ void Foam::isoSurfaceCell::generateTriPoints
                     }
 
                     label fp = f.fcIndex(fp0);
-
                     for (label i = 2; i < f.size(); i++)
                     {
                         label nextFp = f.fcIndex(fp);
@@ -425,22 +449,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
-                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
-                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
-                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
-                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],
@@ -454,22 +474,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
-                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
-                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
-                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
-                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],
@@ -495,7 +511,7 @@ void Foam::isoSurfaceCell::generateTriPoints
 
     if (countNotFoundTets > 0)
     {
-        WarningInFunction
+        WarningIn("Foam::isoSurfaceCell::generateTriPoints(..)")
             << "Could not find " << countNotFoundTets
             << " tet base points, which may lead to inverted triangles."
             << endl;
@@ -507,68 +523,14 @@ void Foam::isoSurfaceCell::generateTriPoints
 
 
 template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::isoSurfaceCell::interpolate
-(
-    const scalarField& cVals,
-    const scalarField& pVals,
-    const Field<Type>& cCoords,
-    const Field<Type>& pCoords
-) const
-{
-    DynamicList<Type> triPoints(nCutCells_);
-    DynamicList<label> triMeshCells(nCutCells_);
-
-    // Dummy snap data
-    DynamicList<Type> snappedPoints;
-    labelList snappedCc(mesh_.nCells(), -1);
-    labelList snappedPoint(mesh_.nPoints(), -1);
-
-
-    generateTriPoints
-    (
-        cVals,
-        pVals,
-
-        cCoords,
-        pCoords,
-
-        snappedPoints,
-        snappedCc,
-        snappedPoint,
-
-        triPoints,
-        triMeshCells
-    );
-
-
-    // One value per point
-    tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
-
-    forAll(triPoints, i)
-    {
-        label mergedPointI = triPointMergeMap_[i];
-
-        if (mergedPointI >= 0)
-        {
-            values[mergedPointI] = triPoints[i];
-        }
-    }
-
-    return tvalues;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
+Foam::tmp<Foam::Field<Type> >
 Foam::isoSurfaceCell::interpolate
 (
     const Field<Type>& cCoords,
     const Field<Type>& pCoords
 ) const
 {
-    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<Type> triPoints(3*nCutCells_);
     DynamicList<label> triMeshCells(nCutCells_);
 
     // Dummy snap data
@@ -576,7 +538,6 @@ Foam::isoSurfaceCell::interpolate
     labelList snappedCc(mesh_.nCells(), -1);
     labelList snappedPoint(mesh_.nPoints(), -1);
 
-
     generateTriPoints
     (
         cVals_,
@@ -593,22 +554,12 @@ Foam::isoSurfaceCell::interpolate
         triMeshCells
     );
 
-
-    // One value per point
-    tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues.ref();
-
-    forAll(triPoints, i)
-    {
-        label mergedPointI = triPointMergeMap_[i];
-
-        if (mergedPointI >= 0)
-        {
-            values[mergedPointI] = triPoints[i];
-        }
-    }
-
-    return tvalues;
+    return isoSurface::interpolate
+    (
+        points().size(),
+        triPointMergeMap_,
+        triPoints
+    );
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index ea92b65..295733c 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -120,7 +120,7 @@ Foam::isoSurface::adaptPatchFields
         {
             fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
             (
-                fld.boundaryField()[patchI]
+                sliceFld.boundaryField()[patchI]
             );
 
             const scalarField& w = mesh.weights().boundaryField()[patchI];
@@ -189,6 +189,8 @@ Type Foam::isoSurface::generatePoint
 }
 
 
+// Note: cannot use simpler isoSurfaceCell::generateTriPoints since
+//       the need here to sometimes pass in remote 'snappedPoints'
 template<class Type>
 void Foam::isoSurface::generateTriPoints
 (
@@ -240,8 +242,8 @@ void Foam::isoSurface::generateTriPoints
         case 0x0F:
         break;
 
-        case 0x0E:
         case 0x01:
+        case 0x0E:
             points.append
             (
                 generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
@@ -254,10 +256,16 @@ void Foam::isoSurface::generateTriPoints
             (
                 generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
             );
+            if (triIndex == 0x0E)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0D:
         case 0x02:
+        case 0x0D:
             points.append
             (
                 generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
@@ -270,97 +278,133 @@ void Foam::isoSurface::generateTriPoints
             (
                 generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
             );
+            if (triIndex == 0x0D)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0C:
         case 0x03:
-        {
-            Type tp1 =
-                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
-            Type tp2 =
-                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
+        case 0x0C:
+            {
+                Type p0p2 =
+                    generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
+                Type p1p3 =
+                    generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
 
-            points.append
-            (
-                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
-            );
-            points.append(tp1);
-            points.append(tp2);
-            points.append(tp2);
-            points.append
-            (
-                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
-            );
-            points.append(tp1);
-        }
+                points.append
+                (
+                    generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+                );
+                points.append(p1p3);
+                points.append(p0p2);
+
+                points.append(p1p3);
+                points.append
+                (
+                    generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+                );
+                points.append(p0p2);
+            }
+            if (triIndex == 0x0C)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-5], points[sz-4]);
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0B:
         case 0x04:
-        {
-            points.append
-            (
-                generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
-            );
-            points.append
-            (
-                generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
-            );
-            points.append
-            (
-                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
-            );
-        }
+        case 0x0B:
+            {
+                points.append
+                (
+                    generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
+                );
+                points.append
+                (
+                    generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
+                );
+                points.append
+                (
+                    generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
+                );
+            }
+            if (triIndex == 0x0B)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0A:
         case 0x05:
-        {
-            Type tp0 =
-                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
-            Type tp1 =
-                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+        case 0x0A:
+            {
+                Type p0p1 =
+                    generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+                Type p2p3 =
+                    generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+
+                points.append(p0p1);
+                points.append(p2p3);
+                points.append
+                (
+                    generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+                );
 
-            points.append(tp0);
-            points.append(tp1);
-            points.append
-            (
-                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
-            );
-            points.append(tp0);
-            points.append
-            (
-                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
-            );
-            points.append(tp1);
-        }
+                points.append(p0p1);
+                points.append
+                (
+                    generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+                );
+                points.append(p2p3);
+            }
+            if (triIndex == 0x0A)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-5], points[sz-4]);
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x09:
         case 0x06:
-        {
-            Type tp0 =
-                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
-            Type tp1 =
-                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+        case 0x09:
+            {
+                Type p0p1 =
+                    generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+                Type p2p3 =
+                    generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 
-            points.append(tp0);
-            points.append
-            (
-                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
-            );
-            points.append(tp1);
-            points.append(tp0);
-            points.append
-            (
-                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
-            );
-            points.append(tp1);
-        }
+                points.append(p0p1);
+                points.append
+                (
+                    generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
+                );
+                points.append(p2p3);
+
+                points.append(p0p1);
+                points.append(p2p3);
+                points.append
+                (
+                    generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
+                );
+            }
+            if (triIndex == 0x09)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-5], points[sz-4]);
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x07:
         case 0x08:
+        case 0x07:
             points.append
             (
                 generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0)
@@ -373,6 +417,12 @@ void Foam::isoSurface::generateTriPoints
             (
                 generatePoint(s3,p3,hasSnap3,snapP3,s1,p1,hasSnap1,snapP1)
             );
+            if (triIndex == 0x07)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
     }
 }
@@ -693,6 +743,43 @@ template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::isoSurface::interpolate
 (
+    const label nPoints,
+    const labelList& triPointMergeMap,
+    const DynamicList<Type>& unmergedValues
+)
+{
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(nPoints, pTraits<Type>::zero));
+    Field<Type>& values = tvalues.ref();
+    labelList nValues(values.size(), 0);
+
+    forAll(unmergedValues, i)
+    {
+        label mergedPointI = triPointMergeMap[i];
+
+        if (mergedPointI >= 0)
+        {
+            values[mergedPointI] += unmergedValues[i];
+            nValues[mergedPointI]++;
+        }
+    }
+
+    forAll(values, i)
+    {
+        if (nValues[i] > 0)
+        {
+            values[i] /= scalar(nValues[i]);
+        }
+    }
+
+    return tvalues;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::isoSurface::interpolate
+(
     const GeometricField<Type, fvPatchField, volMesh>& cCoords,
     const Field<Type>& pCoords
 ) const
@@ -707,7 +794,7 @@ Foam::isoSurface::interpolate
     >> c2(adaptPatchFields(cCoords));
 
 
-    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<Type> triPoints(3*nCutCells_);
     DynamicList<label> triMeshCells(nCutCells_);
 
     // Dummy snap data
@@ -731,52 +818,12 @@ Foam::isoSurface::interpolate
         triMeshCells
     );
 
-
-    // One value per point
-    tmp<Field<Type>> tvalues
+    return interpolate
     (
-        new Field<Type>(points().size(), pTraits<Type>::zero)
+        points().size(),
+        triPointMergeMap_,
+        triPoints
     );
-    Field<Type>& values = tvalues.ref();
-    labelList nValues(values.size(), 0);
-
-    forAll(triPoints, i)
-    {
-        label mergedPointI = triPointMergeMap_[i];
-
-        if (mergedPointI >= 0)
-        {
-            values[mergedPointI] += triPoints[i];
-            nValues[mergedPointI]++;
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "nValues:" << values.size() << endl;
-        label nMult = 0;
-        forAll(nValues, i)
-        {
-            if (nValues[i] == 0)
-            {
-                FatalErrorInFunction
-                    << "point:" << i << " nValues:" << nValues[i]
-                    << abort(FatalError);
-            }
-            else if (nValues[i] > 1)
-            {
-                nMult++;
-            }
-        }
-        Pout<< "Of which mult:" << nMult << endl;
-    }
-
-    forAll(values, i)
-    {
-        values[i] /= scalar(nValues[i]);
-    }
-
-    return tvalues;
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 3f8486e..d789abe 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -124,18 +124,52 @@ void Foam::sampledIsoSurface::getIsoFields() const
     // Get pointField
     // ~~~~~~~~~~~~~~
 
+    // In case of multiple iso values we don't want to calculate multiple e.g.
+    // "volPointInterpolate(p)" so register it and re-use it. This is the
+    // same as the 'cache' functionality from volPointInterpolate but
+    // unfortunately that one does not guarantee that the field pointer
+    // remain: e.g. some other functionObject might delete the cached version.
+    // (volPointInterpolation::interpolate with cache=false deletes any
+    //  registered one or if mesh.changing())
+
     if (!subMeshPtr_.valid())
     {
-        word pointFldName = "volPointInterpolate(" + isoField_ + ')';
+        const word pointFldName =
+            "volPointInterpolate_"
+          + type()
+          + "("
+          + isoField_
+          + ')';
 
         if (fvm.foundObject<pointScalarField>(pointFldName))
         {
             if (debug)
             {
                 InfoInFunction
-                    << "Lookup pointField " << pointFldName << endl;
+                    << "lookup pointField " << pointFldName << endl;
             }
-            pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
+            const pointScalarField& pfld = fvm.lookupObject<pointScalarField>
+            (
+                pointFldName
+            );
+
+            if (!pfld.upToDate(*volFieldPtr_))
+            {
+                if (debug)
+                {
+                    InfoInFunction
+                        << "updating pointField "
+                        << pointFldName << endl;
+                }
+                // Update the interpolated value
+                volPointInterpolation::New(fvm).interpolate
+                (
+                    *volFieldPtr_,
+                    const_cast<pointScalarField&>(pfld)
+                );
+            }
+
+            pointFieldPtr_ = &pfld;
         }
         else
         {
@@ -149,27 +183,20 @@ void Foam::sampledIsoSurface::getIsoFields() const
                     << endl;
             }
 
-            if
+            // Interpolate without cache. Note that we're registering it
+            // below so next time round it goes into the condition
+            // above.
+            tmp<pointScalarField> tpfld
             (
-                storedPointFieldPtr_.empty()
-             || (fvm.time().timeName() != storedPointFieldPtr_().instance())
-            )
-            {
-                if (debug)
-                {
-                    InfoInFunction
-                        << "Interpolating volField " << volFieldPtr_->name()
-                        << " to get pointField " << pointFldName << endl;
-                }
-
-                storedPointFieldPtr_.reset
+                volPointInterpolation::New(fvm).interpolate
                 (
-                    volPointInterpolation::New(fvm)
-                    .interpolate(*volFieldPtr_).ptr()
-                );
-                storedPointFieldPtr_->checkOut();
-                pointFieldPtr_ = storedPointFieldPtr_.operator->();
-            }
+                    *volFieldPtr_,
+                    pointFldName,
+                    false
+                )
+            );
+            pointFieldPtr_ = tpfld.ptr();
+            const_cast<pointScalarField*>(pointFieldPtr_)->store();
         }
 
 
@@ -427,7 +454,6 @@ Foam::sampledIsoSurface::sampledIsoSurface
     prevTimeIndex_(-1),
     storedVolFieldPtr_(NULL),
     volFieldPtr_(NULL),
-    storedPointFieldPtr_(NULL),
     pointFieldPtr_(NULL)
 {
     if (!sampledSurface::interpolate())
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index d0a0b0b..e128560 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -94,7 +94,6 @@ class sampledIsoSurface
             mutable const volScalarField* volFieldPtr_;
 
             //- Cached pointfield
-            mutable autoPtr<pointScalarField> storedPointFieldPtr_;
             mutable const pointScalarField* pointFieldPtr_;
 
             // And on subsetted mesh
-- 
2.7.2

MattijsJ

2016-03-11 17:18

reporter   ~0006031

Attached basically that OpenFOAM-history patch. Should also fix #905 (orientation of triangles)

What was happening is that
- isoSurface stores a reference to the pointField
- this pointField was cached by volPointInterpolation
- so if for any reason volPointInterpolation decided to create a new pointField the reference would become invalid.

I tried to replicate the problem but failed.

henry

2016-03-12 11:54

manager   ~0006032

Resolved in OpenFOAM-dev by commit 36f19c88c4362a28cfc7b13c70ab4f0431416997

Issue History

Date Modified Username Field Change
2015-01-13 15:41 matteoL New Issue
2015-01-14 08:55 user4 Note Added: 0003537
2015-01-14 10:33 matteoL File Added: StlFilesDebug.PNG
2015-01-14 10:41 matteoL Note Added: 0003543
2015-01-14 10:45 matteoL File Added: WithDebug.txt.gz
2015-01-19 20:50 user4 Note Added: 0003567
2015-01-20 10:36 matteoL Note Added: 0003570
2015-02-21 20:52 wyldckat Note Added: 0003851
2015-03-11 12:21 user4 Note Added: 0004082
2015-03-24 00:17 liuhuafei Issue cloned: 0001621
2015-08-20 15:09 matteoL Note Added: 0005277
2016-02-21 20:21 wyldckat Note Added: 0005976
2016-03-11 17:14 MattijsJ File Added: 0001-BUG-isosurface-stores-reference-to-volPointInterpola.patch
2016-03-11 17:18 MattijsJ Note Added: 0006031
2016-03-12 11:54 henry Note Added: 0006032
2016-03-12 11:54 henry Status new => resolved
2016-03-12 11:54 henry Resolution open => fixed
2016-03-12 11:54 henry Assigned To => henry