signed distance surface error

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -45,6 +45,11 @@ namespace Foam
 
 void Foam::distanceSurface::createGeometry()
 {
+    if (debug)
+    {
+        Pout<< "distanceSurface::createGeometry :updating geometry." << endl;
+    }
+
     // Clear any stored topologies
     facesPtr_.clear();
 
@@ -67,7 +72,7 @@ void Foam::distanceSurface::createGeometry()
                 false
             ),
             fvm,
-            dimensionedScalar("zero", dimless/dimTime, 0)
+            dimensionedScalar("zero", dimLength, 0)
         )
     );
     volScalarField& cellDistance = cellDistancePtr_();
@@ -157,6 +162,7 @@ void Foam::distanceSurface::createGeometry()
         }
     }
 
+
     // On processor patches the mesh.C() will already be the cell centre
     // on the opposite side so no need to swap cellDistance.
 
@@ -164,23 +170,70 @@ void Foam::distanceSurface::createGeometry()
     // Distance to points
     pointDistance_.setSize(fvm.nPoints());
     {
+        const pointField& pts = fvm.points();
+
         List<pointIndexHit> nearest;
         surfPtr_().findNearest
         (
-            fvm.points(),
-            scalarField(fvm.nPoints(), GREAT),
+            pts,
+            scalarField(pts.size(), GREAT),
             nearest
         );
-        forAll(pointDistance_, pointI)
+
+        if (signed_)
         {
-            pointDistance_[pointI] = Foam::mag
-            (
-                nearest[pointI].hitPoint()
-              - fvm.points()[pointI]
-            );
+            vectorField normal;
+            surfPtr_().getNormal(nearest, normal);
+
+            forAll(nearest, i)
+            {
+                vector d(pts[i]-nearest[i].hitPoint());
+
+                if ((d&normal[i]) > 0)
+                {
+                    pointDistance_[i] = Foam::mag(d);
+                }
+                else
+                {
+                    pointDistance_[i] = -Foam::mag(d);
+                }
+            }
+        }
+        else
+        {
+            forAll(nearest, i)
+            {
+                pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint());
+            }
         }
     }
 
+
+    if (debug)
+    {
+        Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
+        cellDistance.write();
+        pointScalarField pDist
+        (
+            IOobject
+            (
+                "pointDistance",
+                fvm.time().timeName(),
+                fvm.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            pointMesh::New(fvm),
+            dimensionedScalar("zero", dimLength, 0)
+        );
+        pDist.internalField() = pointDistance_;
+
+        Pout<< "Writing point distance:" << pDist.objectPath() << endl;
+        pDist.write();
+    }
+
+
     //- Direct from cell field and point field.
     isoSurfPtr_.reset
     (
@@ -196,6 +249,7 @@ void Foam::distanceSurface::createGeometry()
     if (debug)
     {
         print(Pout);
+        Pout<< endl;
     }
 }
 
@@ -264,6 +318,13 @@ bool Foam::distanceSurface::needsUpdate() const
 
 bool Foam::distanceSurface::expire()
 {
+    if (debug)
+    {
+        Pout<< "distanceSurface::expire :"
+            << " have-facesPtr_:" << facesPtr_.valid()
+            << " needsUpdate_:" << needsUpdate_ << endl;
+    }
+
     // Clear any stored topologies
     facesPtr_.clear();
 
@@ -280,6 +341,13 @@ bool Foam::distanceSurface::expire()
 
 bool Foam::distanceSurface::update()
 {
+    if (debug)
+    {
+        Pout<< "distanceSurface::update :"
+            << " have-facesPtr_:" << facesPtr_.valid()
+            << " needsUpdate_:" << needsUpdate_ << endl;
+    }
+
     if (!needsUpdate_)
     {
         return false;

src/sampling/sampledSurface/isoSurface/isoSurface.C

@@ -58,6 +58,31 @@ Foam::scalar Foam::isoSurface::isoFraction
 }
 
 
+bool Foam::isoSurface::isEdgeOfFaceCut
+(
+    const scalarField& pVals,
+    const face& f,
+    const bool ownLower,
+    const bool neiLower
+) const
+{
+    forAll(f, fp)
+    {
+        bool fpLower = (pVals[f[fp]] < iso_);
+        if
+        (
+            (fpLower != ownLower)
+         || (fpLower != neiLower)
+         || (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
+        )
+        {
+            return true;
+        }
+    }
+    return false;
+}
+
+
 // Determine for every face/cell whether it (possibly) generates triangles.
 void Foam::isoSurface::calcCutTypes
 (
@@ -84,22 +109,13 @@ void Foam::isoSurface::calcCutTypes
         }
         else
         {
-            // Mesh edge.
+            // See if any mesh edge is cut by looping over all the edges of the
+            // face.
             const face f = mesh_.faces()[faceI];
 
-            forAll(f, fp)
+            if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
             {
-                bool fpLower = (pVals[f[fp]] < iso_);
-                if
-                (
-                    (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
-                 || (fpLower != ownLower)
-                 || (fpLower != neiLower)
-                )
-                {
-                    faceCutType_[faceI] = CUT;
-                    break;
-                }
+                faceCutType_[faceI] = CUT;
             }
         }
     }
@@ -117,22 +133,13 @@ void Foam::isoSurface::calcCutTypes
             {
                 bool ownLower = (cVals[own[faceI]] < iso_);
 
-                // Mesh edge.
                 const face f = mesh_.faces()[faceI];
 
-                forAll(f, fp)
+                if (isEdgeOfFaceCut(pVals, f, ownLower, ownLower))
                 {
-                    bool fpLower = (pVals[f[fp]] < iso_);
-                    if
-                    (
-                        (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
-                     || (fpLower != ownLower)
-                    )
-                    {
-                        faceCutType_[faceI] = CUT;
-                        break;
-                    }
+                    faceCutType_[faceI] = CUT;
                 }
+
                 faceI++;
             }
         }
@@ -152,19 +159,9 @@ void Foam::isoSurface::calcCutTypes
                     // Mesh edge.
                     const face f = mesh_.faces()[faceI];
 
-                    forAll(f, fp)
+                    if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
                     {
-                        bool fpLower = (pVals[f[fp]] < iso_);
-                        if
-                        (
-                            (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
-                         || (fpLower != ownLower)
-                         || (fpLower != neiLower)
-                        )
-                        {
-                            faceCutType_[faceI] = CUT;
-                            break;
-                        }
+                        faceCutType_[faceI] = CUT;
                     }
                 }
                 faceI++;
@@ -1355,6 +1352,30 @@ Foam::isoSurface::isoSurface
     iso_(iso),
     mergeDistance_(mergeTol*mesh_.bounds().mag())
 {
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const labelList& own = mesh_.faceOwner();
+
+    // Check
+    forAll(patches, patchI)
+    {
+        if (isA<emptyPolyPatch>(patches[patchI]))
+        {
+            FatalErrorIn
+            (
+                "isoSurface::isoSurface\n"
+                "(\n"
+                "    const volScalarField& cVals,\n"
+                "    const scalarField& pVals,\n"
+                "    const scalar iso,\n"
+                "    const bool regularise,\n"
+                "    const scalar mergeTol\n"
+                ")\n"
+            )   << "Iso surfaces not supported on case with empty patches."
+                << exit(FatalError);
+        }
+    }
+
+
     // Determine if any cut through face/cell
     calcCutTypes(cVals, pVals);
 
@@ -1364,8 +1385,6 @@ Foam::isoSurface::isoSurface
     PackedList<1> isBoundaryPoint(mesh_.nPoints());
 
     labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    const labelList& own = mesh_.faceOwner();
 
     forAll(patches, patchI)
     {

src/sampling/sampledSurface/isoSurface/isoSurface.H

@@ -31,6 +31,9 @@ Description
     G.M. Treece, R.W. Prager and A.H. Gee.
 
     Note:
+    - not possible on patches of type 'empty'. There are no values on
+    'empty' patch fields so even the api would have to change 
+     (no volScalarField as argument). Too messy.
     - in parallel the regularisation (coarsening) always takes place
     and slightly different surfaces will be created compared to non-parallel.
     The surface will still be continuous though!
@@ -123,6 +126,15 @@ class isoSurface
             const scalar s1
         ) const;
 
+        //- Check if any edge of a face is cut
+        bool isEdgeOfFaceCut
+        (
+            const scalarField& pVals,
+            const face& f,
+            const bool ownLower,
+            const bool neiLower
+        ) const;
+
         //- Set faceCutType,cellCutType.
         void calcCutTypes
         (
@@ -217,7 +229,7 @@ class isoSurface
         ) const;
 
         template<class Type>
-        label generateTriPoints
+        label generateFaceTriPoints
         (
             const volScalarField& cVals,
             const scalarField& pVals,

src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C

@@ -200,7 +200,7 @@ void Foam::isoSurface::generateTriPoints
 
 
 template<class Type>
-Foam::label Foam::isoSurface::generateTriPoints
+Foam::label Foam::isoSurface::generateFaceTriPoints
 (
     const volScalarField& cVals,
     const scalarField& pVals,
@@ -288,6 +288,19 @@ void Foam::isoSurface::generateTriPoints
     const labelList& own = mesh_.faceOwner();
     const labelList& nei = mesh_.faceNeighbour();
 
+    if
+    (
+        (cVals.size() != mesh_.nCells())
+     || (pVals.size() != mesh_.nPoints())
+     || (cCoords.size() != mesh_.nCells())
+     || (pCoords.size() != mesh_.nPoints())
+     || (snappedCc.size() != mesh_.nCells())
+     || (snappedPoint.size() != mesh_.nPoints())
+    )
+    {
+        FatalErrorIn("isoSurface::generateTriPoints(..)")
+            << "Incorrect size." << abort(FatalError);
+    }
 
     // Determine neighbouring snap status
     labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
@@ -319,7 +332,7 @@ void Foam::isoSurface::generateTriPoints
     {
         if (faceCutType_[faceI] != NOTCUT)
         {
-            generateTriPoints
+            generateFaceTriPoints
             (
                 cVals,
                 pVals,
@@ -357,7 +370,7 @@ void Foam::isoSurface::generateTriPoints
                 {
                     if (faceCutType_[faceI] != NOTCUT)
                     {
-                        generateTriPoints
+                        generateFaceTriPoints
                         (
                             cVals,
                             pVals,
@@ -384,14 +397,14 @@ void Foam::isoSurface::generateTriPoints
         }
         else if (isA<emptyPolyPatch>(pp))
         {
-            // Assume zero-gradient.
+            // Assume zero-gradient. But what about coordinates?
             label faceI = pp.start();
 
             forAll(pp, i)
             {
                 if (faceCutType_[faceI] != NOTCUT)
                 {
-                    generateTriPoints
+                    generateFaceTriPoints
                     (
                         cVals,
                         pVals,
@@ -423,7 +436,7 @@ void Foam::isoSurface::generateTriPoints
             {
                 if (faceCutType_[faceI] != NOTCUT)
                 {
-                    generateTriPoints
+                    generateFaceTriPoints
                     (
                         cVals,
                         pVals,