ENH: isoSurfaceTopo: replacement for isoSurfaceCell.

src/sampling/Make/files

@@ -29,6 +29,7 @@ surface/cutting/cuttingSurfaceBaseSelection.C
 surface/distanceSurface/distanceSurface.C
 surface/isoSurface/isoSurface.C
 surface/isoSurface/isoSurfaceCell.C
+surface/isoSurface/isoSurfaceTopo.C
 surface/thresholdCellFaces/thresholdCellFaces.C
 surface/triSurfaceMesh/discreteSurface.C
 
@@ -44,6 +45,7 @@ sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C
 sampledSurface/sampledPlane/sampledPlane.C
 sampledSurface/isoSurface/sampledIsoSurface.C
 sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
 sampledSurface/distanceSurface/sampledDistanceSurface.C
 sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
 sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C

@@ -0,0 +1,328 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledIsoSurfaceTopo.H"
+#include "dictionary.H"
+#include "volFields.H"
+#include "volPointInterpolation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "isoSurfaceTopo.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sampledIsoSurfaceTopo, 0);
+    addNamedToRunTimeSelectionTable
+    (
+        sampledSurface,
+        sampledIsoSurfaceTopo,
+        word,
+        isoSurfaceTopo
+    );
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
+{
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+    // No update needed
+    if (fvm.time().timeIndex() == prevTimeIndex_)
+    {
+        return false;
+    }
+
+    prevTimeIndex_ = fvm.time().timeIndex();
+
+    // Clear derived data
+    sampledSurface::clearGeom();
+
+    // Use field from database, or try to read it in
+
+    const auto* cellFldPtr = fvm.findObject<volScalarField>(isoField_);
+
+    if (debug)
+    {
+        if (cellFldPtr)
+        {
+            InfoInFunction << "Lookup " << isoField_ << endl;
+        }
+        else
+        {
+            InfoInFunction
+                << "Reading " << isoField_
+                << " from time " << fvm.time().timeName()
+                << endl;
+        }
+    }
+
+    // For holding the volScalarField read in.
+    autoPtr<volScalarField> fieldReadPtr;
+
+    if (!cellFldPtr)
+    {
+        // Bit of a hack. Read field and store.
+
+        fieldReadPtr = autoPtr<volScalarField>::New
+        (
+            IOobject
+            (
+                isoField_,
+                fvm.time().timeName(),
+                fvm,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fvm
+        );
+    }
+
+    const volScalarField& cellFld =
+        (fieldReadPtr.valid() ? *fieldReadPtr : *cellFldPtr);
+
+    auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
+
+    //- Direct from cell field and point field. Gives bad continuity.
+    isoSurfaceTopo surf
+    (
+        fvm,
+        cellFld.primitiveField(),
+        tpointFld().primitiveField(),
+        isoVal_,
+        (regularise_ ? isoSurfaceTopo::DIAGCELL : isoSurfaceTopo::NONE)
+    );
+
+    MeshedSurface<face>& mySurface = const_cast<sampledIsoSurfaceTopo&>(*this);
+
+    mySurface.transfer(static_cast<meshedSurface&>(surf));
+    meshCells_ = std::move(surf.meshCells());
+
+    // triangulate uses remapFaces()
+    // - this is somewhat less efficient since it recopies the faces
+    // that we just created, but we probably don't want to do this
+    // too often anyhow.
+    if (triangulate_)
+    {
+        labelList faceMap;
+        mySurface.triangulate(faceMap);
+        meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+    }
+
+    if (debug)
+    {
+        Pout<< "sampledIsoSurfaceTopo::updateGeometry() : constructed iso:"
+            << nl
+            << "    regularise     : " << regularise_ << nl
+            << "    triangulate    : " << triangulate_ << nl
+            << "    isoField       : " << isoField_ << nl
+            << "    isoValue       : " << isoVal_ << nl
+            << "    points         : " << points().size() << nl
+            << "    faces          : " << MeshStorage::size() << nl
+            << "    cut cells      : " << meshCells_.size() << endl;
+    }
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    MeshStorage(),
+    isoField_(dict.get<word>("isoField")),
+    isoVal_(dict.get<scalar>("isoValue")),
+    regularise_(dict.lookupOrDefault("regularise", true)),
+    triangulate_(dict.lookupOrDefault("triangulate", false)),
+    prevTimeIndex_(-1),
+    meshCells_()
+{
+    if (triangulate_ && !regularise_)
+    {
+        FatalIOErrorInFunction(dict) << "Cannot both use regularise"
+            << " and triangulate" << exit(FatalIOError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfaceTopo::~sampledIsoSurfaceTopo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::sampledIsoSurfaceTopo::needsUpdate() const
+{
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+    return fvm.time().timeIndex() != prevTimeIndex_;
+}
+
+
+bool Foam::sampledIsoSurfaceTopo::expire()
+{
+    // Clear derived data
+    sampledSurface::clearGeom();
+
+    // Already marked as expired
+    if (prevTimeIndex_ == -1)
+    {
+        return false;
+    }
+
+    // Force update
+    prevTimeIndex_ = -1;
+    return true;
+}
+
+
+bool Foam::sampledIsoSurfaceTopo::update()
+{
+    return updateGeometry();
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<scalar>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<vector>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<sphericalTensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<symmTensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurfaceTopo::sample
+(
+    const interpolation<tensor>& sampler
+) const
+{
+    return sampleOnFaces(sampler);
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurfaceTopo::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return sampleOnPoints(interpolator);
+}
+
+
+void Foam::sampledIsoSurfaceTopo::print(Ostream& os) const
+{
+    os  << "sampledIsoSurfaceTopo: " << name() << " :"
+        << "  field:" << isoField_
+        << "  value:" << isoVal_;
+        //<< "  faces:" << faces().size()   // possibly no geom yet
+        //<< "  points:" << points().size();
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H

@@ -0,0 +1,291 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sampledIsoSurfaceTopo
+
+Description
+    A sampledSurface defined by a surface of iso value.
+    To be used in sampleSurfaces / functionObjects. Recalculates iso surface
+    only if time changes.
+
+    This is often embedded as part of a sampled surfaces function object.
+
+Usage
+    Example of function object partial specification:
+    \verbatim
+    surfaces
+    (
+        surface1
+        {
+            type    isoSurfaceTopo;
+            isoField        p;
+            isoValue        0.0;
+        }
+    );
+    \endverbatim
+
+    Where the sub-entries comprise:
+    \table
+        Property | Description                             | Required | Default
+        type     | isoSurfaceTopo                          | yes      |
+        isoField | field name for obtaining iso-surface    | yes      |
+        isoValue | value of iso-surface                    | yes      |
+        regularise | filter faces                          | no       | true
+        triangulate | triangulate faces (if regularise)    | no       | false
+    \endtable
+
+Note
+    Does not currently support cell zones.
+
+SourceFiles
+    sampledIsoSurfaceTopo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledIsoSurfaceTopo_H
+#define sampledIsoSurfaceTopo_H
+
+#include "sampledSurface.H"
+#include "MeshedSurface.H"
+#include "MeshedSurfacesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class sampledIsoSurfaceTopo Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledIsoSurfaceTopo
+:
+    public sampledSurface,
+    public MeshedSurface<face>
+{
+    // Private typedefs for convenience
+    typedef MeshedSurface<face> MeshStorage;
+
+    // Private data
+
+        //- Field to get isoSurface of
+        const word isoField_;
+
+        //- Iso value
+        const scalar isoVal_;
+
+        //- Whether to coarse
+        const bool regularise_;
+
+        //- Whether to triangulate
+        const bool triangulate_;
+
+    // Recreated for every isoSurface
+
+        //- Time at last call, also track it surface needs an update
+        mutable label prevTimeIndex_;
+
+        //- For every triangle/face the original cell in mesh
+        mutable labelList meshCells_;
+
+
+    // Private Member Functions
+
+        //- Create iso surface (if time has changed)
+        //  Do nothing (and return false) if no update was needed
+        bool updateGeometry() const;
+
+        //- Sample volume field onto surface faces
+        template<class Type>
+        tmp<Field<Type>> sampleOnFaces
+        (
+            const interpolation<Type>& sampler
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        template<class Type>
+        tmp<Field<Type>> sampleOnPoints
+        (
+            const interpolation<Type>& interpolator
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sampledIsoSurfaceTopo");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        sampledIsoSurfaceTopo
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~sampledIsoSurfaceTopo();
+
+
+    // Member Functions
+
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surface as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if surface was already marked as expired.
+        virtual bool expire();
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was needed
+        virtual bool update();
+
+
+        //- Points of surface
+        virtual const pointField& points() const
+        {
+            return MeshStorage::points();
+        }
+
+        //- Faces of surface
+        virtual const faceList& faces() const
+        {
+            return *this;
+        }
+
+        //- Const access to per-face zone/region information
+        virtual const labelList& zoneIds() const
+        {
+            return Foam::emptyLabelList;
+        }
+
+        //- Face area magnitudes
+        virtual const vectorField& Sf() const
+        {
+            return MeshStorage::Sf();
+        }
+
+        //- Face area magnitudes
+        virtual const scalarField& magSf() const
+        {
+            return MeshStorage::magSf();
+        }
+
+        //- Face centres
+        virtual const vectorField& Cf() const
+        {
+            return MeshStorage::Cf();
+        }
+
+
+    // Sample
+
+        //- Sample volume field onto surface faces
+        virtual tmp<scalarField> sample
+        (
+            const interpolation<scalar>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<vectorField> sample
+        (
+            const interpolation<vector>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<sphericalTensorField> sample
+        (
+            const interpolation<sphericalTensor>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<symmTensorField> sample
+        (
+            const interpolation<symmTensor>& sampler
+        ) const;
+
+        //- Sample volume field onto surface faces
+        virtual tmp<tensorField> sample
+        (
+            const interpolation<tensor>& sampler
+        ) const;
+
+
+    // Interpolate
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<scalarField> interpolate
+        (
+            const interpolation<scalar>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<vectorField> interpolate
+        (
+            const interpolation<vector>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<sphericalTensorField> interpolate
+        (
+            const interpolation<sphericalTensor>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<symmTensorField> interpolate
+        (
+            const interpolation<symmTensor>& interpolator
+        ) const;
+
+        //- Interpolate volume field onto surface points
+        virtual tmp<tensorField> interpolate
+        (
+            const interpolation<tensor>& interpolator
+        ) const;
+
+        //- Write
+        virtual void print(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "sampledIsoSurfaceTopoTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledIsoSurfaceTopo.H"
+#include "isoSurface.H"
+#include "volFieldsFwd.H"
+#include "pointFields.H"
+#include "volPointInterpolation.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledIsoSurfaceTopo::sampleOnFaces
+(
+    const interpolation<Type>& sampler
+) const
+{
+    updateGeometry();  // Recreate geometry if time has changed
+
+    return sampledSurface::sampleOnFaces
+    (
+        sampler,
+        meshCells_,
+        faces(),
+        points()
+    );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledIsoSurfaceTopo::sampleOnPoints
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    updateGeometry();  // Recreate geometry if time has changed
+
+    const labelList& elements = meshCells_;
+
+    // One value per point
+    auto tvalues = tmp<Field<Type>>::New(points().size());
+    auto& values = tvalues.ref();
+
+    const faceList& fcs = faces();
+    const pointField& pts = points();
+
+    bitSet pointDone(points().size());
+
+    forAll(faces(), cutFacei)
+    {
+        const face& f = fcs[cutFacei];
+        const label celli = elements[cutFacei];
+
+        for (const label pointi : f)
+        {
+            if (pointDone.set(pointi))
+            {
+                values[pointi] = interpolator.interpolate
+                (
+                    pts[pointi],
+                    celli
+                );
+            }
+        }
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //

src/sampling/surface/isoSurface/isoSurfaceTopo.H

@@ -0,0 +1,266 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::isoSurfaceTopo
+
+Description
+    Marching tet iso surface algorithm with optional filtering to keep only
+    points originating from mesh edges.
+
+SourceFiles
+    isoSurfaceTopo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoSurfaceTopo_H
+#define isoSurfaceTopo_H
+
+#include "labelPair.H"
+#include "pointIndexHit.H"
+#include "PackedBoolList.H"
+#include "MeshedSurface.H"
+#include "edgeList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+class tetMatcher;
+
+/*---------------------------------------------------------------------------*\
+                       Class isoSurfaceTopo Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoSurfaceTopo
+:
+    public MeshedSurface<face>
+{
+    // Private typedefs for convenience
+    typedef MeshedSurface<face> MeshStorage;
+
+public:
+
+        enum filterType
+        {
+            NONE,       // No filtering
+            DIAGCELL,   // Remove points from face-diagonal and pyramid
+                        // (vertex to cell-centre) edges
+            CELL        // Only remove points from pyramid edges
+        };
+
+
+private:
+
+    // Private data
+
+        enum cellCutType
+        {
+            NOTCUT,     // Not cut
+            SPHERE,     // All edges to cell centre cut
+            CUT         // Normal cut
+        };
+
+
+        //- Reference to mesh
+        const polyMesh& mesh_;
+
+        const scalarField& cVals_;
+
+        const scalarField& pVals_;
+
+        //- Iso value
+        const scalar iso_;
+
+        //- Per point: originating mesh vertex/cc. See encoding above
+        edgeList pointToVerts_;
+
+        //- For every face the original cell in mesh
+        labelList meshCells_;
+
+        //- For every point the originating face in mesh
+        labelList pointToFace_;
+
+
+    // Private Member Functions
+
+        //- Does any edge of triangle cross iso value?
+        bool isTriCut
+        (
+            const triFace& tri,
+            const scalarField& pointValues
+        ) const;
+
+        //- Determine whether cell is cut
+        cellCutType calcCutType
+        (
+            const bool isTet,
+            const label
+        ) const;
+
+        //- Determine for all mesh whether cell is cut
+        label calcCutTypes
+        (
+            tetMatcher& tet,
+            List<cellCutType>& cellCutTypes
+        );
+
+        //- Generate single point on edge
+        label generatePoint
+        (
+            const label facei,
+            const bool edgeIsDiag,
+            const edge& vertices,
+
+            DynamicList<edge>& pointToVerts,
+            DynamicList<label>& pointToFace,
+            DynamicList<bool>& pointFromDiag,
+            EdgeMap<label>& vertsToPoint
+        ) const;
+
+        //- Generate triangles from tet
+        void generateTriPoints
+        (
+            const label facei,
+            const FixedList<scalar, 4>& s,
+            const FixedList<point, 4>& p,
+            const FixedList<label, 4>& pIndex,
+            const FixedList<bool, 6>& edgeIsDiag,
+
+            DynamicList<edge>& pointToVerts,
+            DynamicList<label>& pointToFace,
+            DynamicList<bool>& pointFromDiag,
+
+            EdgeMap<label>& vertsToPoint,
+            DynamicList<label>& verts
+        ) const;
+
+        //- Generate triangles from cell
+        void generateTriPoints
+        (
+            const polyMesh& mesh,
+            const label celli,
+            const bool isTet,
+
+            DynamicList<edge>& pointToVerts,
+            DynamicList<label>& pointToFace,
+            DynamicList<bool>& pointFromDiag,
+
+            EdgeMap<label>& vertsToPoint,
+            DynamicList<label>& verts,
+            DynamicList<label>& faceLabels
+        ) const;
+
+
+        // Simplification
+
+            void triangulateOutside
+            (
+                const bool filterDiag,
+                const PrimitivePatch<face, SubList, const pointField&>& pp,
+                const boolList& pointFromDiag,
+                const labelList& pointToFace,
+                const label cellID,
+
+                DynamicList<face>& compactFaces,
+                DynamicList<label>& compactCellIDs
+            ) const;
+
+            MeshStorage removeInsidePoints
+            (
+                const bool filterDiag,
+                const MeshStorage& s,
+                const boolList& pointFromDiag,
+                const labelList& pointToFace,
+                const labelList& start,              // Per cell:starting tri
+                DynamicList<label>& pointCompactMap, // Per point the original
+                DynamicList<label>& compactCellIDs   // Per face the cellID
+            ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("isoSurfaceTopo");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        isoSurfaceTopo
+        (
+            const polyMesh& mesh,
+            const scalarField& cellValues,
+            const scalarField& pointValues,
+            const scalar iso,
+            const filterType filter = DIAGCELL
+        );
+
+
+    // Member Functions
+
+        //- For every face original cell in mesh
+        const labelList& meshCells() const
+        {
+            return meshCells_;
+        }
+
+        //- For every point originating face (pyramid) in mesh
+        const labelList& pointToFace() const
+        {
+            return pointToFace_;
+        }
+
+        //- Per point: originating mesh vertex/cc. See encoding above<76>
+        const edgeList& pointToVerts() const
+        {
+            return pointToVerts_;
+        }
+
+        //- Interpolates cCoords,pCoords.
+        template<class Type>
+        tmp<Field<Type>> interpolate
+        (
+            const Field<Type>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "isoSurfaceTopoTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C

@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::isoSurfaceTopo::interpolate
+(
+    const Field<Type>& cellCoords,
+    const Field<Type>& pointCoords
+) const
+{
+    tmp<Field<Type>> tfld(new Field<Type>(pointToVerts_.size()));
+    Field<Type>& fld = tfld.ref();
+
+    forAll(pointToVerts_, i)
+    {
+        scalar s0;
+        Type p0;
+        {
+            label v0 = pointToVerts_[i][0];
+            if (v0 < mesh_.nPoints())
+            {
+                s0 = pVals_[v0];
+                p0 = pointCoords[v0];
+            }
+            else
+            {
+                label celli = v0-mesh_.nPoints();
+                s0 = cVals_[celli];
+                p0 = cellCoords[celli];
+            }
+        }
+
+        scalar s1;
+        Type p1;
+        {
+            label v1 = pointToVerts_[i][1];
+            if (v1 < mesh_.nPoints())
+            {
+                s1 = pVals_[v1];
+                p1 = pointCoords[v1];
+            }
+            else
+            {
+                label celli = v1-mesh_.nPoints();
+                s1 = cVals_[celli];
+                p1 = cellCoords[celli];
+            }
+        }
+
+        scalar d = s1-s0;
+        if (mag(d) > VSMALL)
+        {
+            scalar s = (iso_-s0)/d;
+            fld[i] = s*p1+(1.0-s)*p0;
+        }
+        else
+        {
+            fld[i] = 0.5*(p0+p1);
+        }
+    }
+
+    return tfld;
+}
+
+
+// ************************************************************************* //