isoSurface instead of isoSurfaceCell

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -30,7 +30,8 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "isoSurfaceCell.H"
+#include "isoSurface.H"
+//#include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -42,19 +43,46 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::distanceSurface::createGeometry() const
+void Foam::distanceSurface::createGeometry()
 {
     // Clear any stored topo
     facesPtr_.clear();
 
+
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
     // Distance to cell centres
-    scalarField cellDistance(mesh().nCells());
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    cellDistancePtr_.reset
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "cellDistance",
+                fvm.time().timeName(),
+                fvm.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fvm,
+            dimensionedScalar("zero", dimless/dimTime, 0)
+        )
+    );
+    volScalarField& cellDistance = cellDistancePtr_();
+
+    // Internal field
     {
+        const pointField& cc = fvm.C();
+        scalarField& fld = cellDistance.internalField();
+
         List<pointIndexHit> nearest;
         surfPtr_().findNearest
         (
-            mesh().cellCentres(),
+            cc,
-            scalarField(mesh().nCells(), GREAT),
+            scalarField(cc.size(), GREAT),
             nearest
         );
 
@@ -63,98 +91,109 @@ void Foam::distanceSurface::createGeometry() const
             vectorField normal;
             surfPtr_().getNormal(nearest, normal);
 
-            forAll(cellDistance, cellI)
+            forAll(nearest, i)
             {
-                vector d(mesh().cellCentres()[cellI]-nearest[cellI].hitPoint());
+                vector d(cc[i]-nearest[i].hitPoint());
 
-                if ((d&normal[cellI]) > 0)
+                if ((d&normal[i]) > 0)
                 {
-                    cellDistance[cellI] = Foam::mag(d);
+                    fld[i] = Foam::mag(d);
                 }
                 else
                 {
-                    cellDistance[cellI] = -Foam::mag(d);
+                    fld[i] = -Foam::mag(d);
                 }
             }
         }
         else
         {
-            forAll(cellDistance, cellI)
+            forAll(nearest, i)
             {
-                cellDistance[cellI] = Foam::mag
+                fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
-                (
-                    nearest[cellI].hitPoint()
-                  - mesh().cellCentres()[cellI]
-                );
             }
         }
     }
 
+    // Patch fields
+    {
+        forAll(fvm.C().boundaryField(), patchI)
+        {
+            const pointField& cc = fvm.C().boundaryField()[patchI];
+            fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
+
+            List<pointIndexHit> nearest;
+            surfPtr_().findNearest
+            (
+                cc,
+                scalarField(cc.size(), GREAT),
+                nearest
+            );
+
+            if (signed_)
+            {
+                vectorField normal;
+                surfPtr_().getNormal(nearest, normal);
+
+                forAll(nearest, i)
+                {
+                    vector d(cc[i]-nearest[i].hitPoint());
+
+                    if ((d&normal[i]) > 0)
+                    {
+                        fld[i] = Foam::mag(d);
+                    }
+                    else
+                    {
+                        fld[i] = -Foam::mag(d);
+                    }
+                }
+            }
+            else
+            {
+                forAll(nearest, i)
+                {
+                    fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
+                }
+            }
+        }
+    }
+
+    // On processor patches the mesh.C() will already be the cell centre
+    // on the opposite side so no need to swap cellDistance.
+
+
     // Distance to points
-    scalarField pointDistance(mesh().nPoints());
+    pointDistance_.setSize(fvm.nPoints());
     {
         List<pointIndexHit> nearest;
         surfPtr_().findNearest
         (
-            mesh().points(),
+            fvm.points(),
-            scalarField(mesh().nPoints(), GREAT),
+            scalarField(fvm.nPoints(), GREAT),
             nearest
         );
-        forAll(pointDistance, pointI)
+        forAll(pointDistance_, pointI)
         {
-            pointDistance[pointI] = Foam::mag
+            pointDistance_[pointI] = Foam::mag
             (
                 nearest[pointI].hitPoint()
-              - mesh().points()[pointI]
+              - fvm.points()[pointI]
             );
         }
     }
 
     //- Direct from cell field and point field.
-    const isoSurfaceCell iso
+    isoSurfPtr_.reset
     (
-        mesh(),
+        new isoSurface
-        cellDistance,
+        (
-        pointDistance,
+            cellDistance,
-        distance_,
+            pointDistance_,
-        regularise_
+            distance_,
+            regularise_
+        )
     );
 
-    ////- From point field and interpolated cell.
-    //scalarField cellAvg(mesh().nCells(), scalar(0.0));
-    //labelField nPointCells(mesh().nCells(), 0);
-    //{
-    //    for (label pointI = 0; pointI < mesh().nPoints(); pointI++)
-    //    {
-    //        const labelList& pCells = mesh().pointCells(pointI);
-    //
-    //        forAll(pCells, i)
-    //        {
-    //            label cellI = pCells[i];
-    //
-    //            cellAvg[cellI] += pointDistance[pointI];
-    //            nPointCells[cellI]++;
-    //        }
-    //    }
-    //}
-    //forAll(cellAvg, cellI)
-    //{
-    //    cellAvg[cellI] /= nPointCells[cellI];
-    //}
-    //
-    //const isoSurface iso
-    //(
-    //    mesh(),
-    //    cellAvg,
-    //    pointDistance,
-    //    distance_,
-    //    regularise_
-    //);
-
-
-    const_cast<distanceSurface&>(*this).triSurface::operator=(iso);
-    meshCells_ = iso.meshCells();
-
     if (debug)
     {
         print(Pout);
@@ -193,9 +232,8 @@ Foam::distanceSurface::distanceSurface
     signed_(readBool(dict.lookup("signed"))),
     regularise_(dict.lookupOrDefault("regularise", true)),
     zoneName_(word::null),
-    facesPtr_(NULL),
+    isoSurfPtr_(NULL),
-    storedTimeIndex_(-1),
+    facesPtr_(NULL)
-    meshCells_(0)
 {
 //    label zoneId = -1;
 //    if (dict.readIfPresent("zone", zoneName_))

src/sampling/sampledSurface/distanceSurface/distanceSurface.H

@@ -37,8 +37,9 @@ SourceFiles
 #define distanceSurface_H
 
 #include "sampledSurface.H"
-#include "triSurface.H"
 #include "searchableSurface.H"
+//#include "isoSurfaceCell.H"
+#include "isoSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,8 +52,7 @@ namespace Foam
 
 class distanceSurface
 :
-    public sampledSurface,
+    public sampledSurface
-    public triSurface
 {
     // Private data
 
@@ -71,23 +71,24 @@ class distanceSurface
         //- zone name (if restricted to zones)
         word zoneName_;
 
+
+        //- Distance to cell centres
+        autoPtr<volScalarField> cellDistancePtr_;
+
+        //- Distance to points
+        scalarField pointDistance_;
+
+        //- Constructed iso surface
+        autoPtr<isoSurface> isoSurfPtr_;
+
         //- triangles converted to faceList
         mutable autoPtr<faceList> facesPtr_;
 
 
-        // Recreated for every isoSurface
-
-            //- Time at last call
-            mutable label storedTimeIndex_;
-
-            //- For every triangle the original cell in mesh
-            mutable labelList meshCells_;
-
-
     // Private Member Functions
 
-        //- Create iso surface (if time has changed)
+        //- Create iso surface
-        void createGeometry() const;
+        void createGeometry();
 
         //- sample field on faces
         template <class Type>
@@ -129,7 +130,7 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return triSurface::points();
+            return surface().points();
         }
 
         //- Faces of surface
@@ -137,7 +138,7 @@ public:
         {
             if (!facesPtr_.valid())
             {
-                const triSurface& s = *this;
+                const triSurface& s = surface();
 
                 facesPtr_.reset(new faceList(s.size()));
 
@@ -153,6 +154,10 @@ public:
         //- Correct for mesh movement and/or field changes
         virtual void correct(const bool meshChanged);
 
+        const isoSurface& surface() const
+        {
+            return isoSurfPtr_();
+        }
 
         //- sample field on surface
         virtual tmp<scalarField> sample

src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C

@@ -39,7 +39,7 @@ Foam::distanceSurface::sampleField
     const GeometricField<Type, fvPatchField, volMesh>& vField
 ) const
 {
-    return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
+    return tmp<Field<Type> >(new Field<Type>(vField, surface().meshCells()));
 }
 
 
@@ -50,33 +50,25 @@ Foam::distanceSurface::interpolateField
     const interpolation<Type>& interpolator
 ) const
 {
-    // One value per point
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
 
-    boolList pointDone(points().size(), false);
+    // Get fields to sample. Assume volPointInterpolation!
+    const GeometricField<Type, fvPatchField, volMesh>& volFld =
+        interpolator.psi();
 
-    forAll(faces(), cutFaceI)
+    tmp<GeometricField<Type, pointPatchField, pointMesh> > pointFld
-    {
+    (
-        const face& f = faces()[cutFaceI];
+        volPointInterpolation::New(fvm).interpolate(volFld)
+    );
 
-        forAll(f, faceVertI)
+    // Sample.
-        {
+    return surface().interpolate
-            label pointI = f[faceVertI];
+    (
-
+        cellDistancePtr_(),
-            if (!pointDone[pointI])
+        pointDistance_,
-            {
+        volFld,
-                values[pointI] = interpolator.interpolate
+        pointFld()
-                (
+    );
-                    points()[pointI],
-                    meshCells_[cutFaceI]
-                );
-                pointDone[pointI] = true;
-            }
-        }
-    }
-
-    return tvalues;
 }