cutPoly: New polyhedral cutting routines and isoSurface algorithm

A set of routines for cutting polyhedra have been added. These can cut
polyhedral cells based on the adjacent point values and an iso-value
which defines the surface. The method operates directly on the
polyhedral cells; it does not decompose them into tetrahedra at any
point. The routines can compute the cut topology as well as integrals of
properties above and below the cut surface.

An iso-surface algorithm has been added based on these polyhedral
cutting routines. It is significantly more robust than the previous
algorithm, and produces compact surfaces equivalent to the previous
algorithm's maximum filtering level. It is also approximately 3 times
faster than the previous algorithm, and 10 times faster when run
repeatedly on the same set of cells (this is because some addressing is
cached and reused).

This algorithm is used by the 'isoSurface', 'distanceSurface' and
'cutPlane' sampled surfaces.

The 'cutPlane' sampled surface is a renaming of 'cuttingPlane' to make
it consistent with the corresponding packaged function. The name
'cuttingPlane' has been retained for backwards compatibility and can
still be used to select a 'cutPlane' surface. The legacy 'plane' surface
has been removed.

The 'average' keyword has been removed from specification of these
sampled surfaces as cell-centred values are no longer used in the
generation of or interpolation to an iso-surface. The 'filtering'
keyword has also been removed as it relates to options within the
previous algorithm. Zone support has been reinstated into the
'isoSurface' sampled surface. Interpolation to all these sampled
surfaces has been corrected to exactly match the user-selected
interpolation scheme, and the interpolation procedure no longer
unnecessarily re-generates data that is already available.

applications/test/tetTetOverlap/Test-tetTetOverlap.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,7 +32,7 @@ Description
 #include "tetPointRef.H"
 #include "OFstream.H"
 #include "meshTools.H"
-#include "cut.H"
+#include "cutTriTet.H"
 
 using namespace Foam;
 
@@ -87,9 +87,9 @@ int main(int argc, char *argv[])
     // Do intersection
     typedef DynamicList<FixedList<point, 4>> tetList;
     tetList tetsIn1, tetsIn2, tetsOut;
-    cut::appendOp<tetList> tetOpIn1(tetsIn1);
-    cut::appendOp<tetList> tetOpIn2(tetsIn2);
-    cut::appendOp<tetList> tetOpOut(tetsOut);
+    cutTriTet::appendOp<tetList> tetOpIn1(tetsIn1);
+    cutTriTet::appendOp<tetList> tetOpIn2(tetsIn2);
+    cutTriTet::appendOp<tetList> tetOpOut(tetsOut);
 
     const plane p0(tetB[1], tetB[3], tetB[2]);
     tetsIn1.clear();

applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C

@@ -41,7 +41,7 @@ Description
 #include "cellShape.H"
 #include "cellModeller.H"
 #include "DynamicField.H"
-#include "isoSurface.H"
+#include "cutPolyIsoSurface.H"
 #include "vtkSurfaceWriter.H"
 #include "syncTools.H"
 
@@ -707,13 +707,7 @@ int main(int argc, char *argv[])
             }
         }
 
-        isoSurface iso
-        (
-            mesh,
-            cellDistance,
-            pointDistance,
-            0
-        );
+        cutPolyIsoSurface iso(mesh, pointDistance, 0);
 
         isoFaces.setSize(iso.size());
         forAll(isoFaces, i)

etc/caseDicts/postProcessing/surface/cutPlaneSurface.cfg

@@ -12,7 +12,7 @@ surfaces
 (
     cutPlane
     {
-        type        cuttingPlane;
+        type        cutPlane;
         interpolate $interpolate;
         planeType   pointAndNormal;
         point       $point;

src/OpenFOAM/containers/HashTables/HashList/HashList.C

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "HashList.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    template<class Type, class Key, class Hash>
+    const Key HashList<Type, Key, Hash>::nullKey = Key::null;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class Key, class Hash>
+Foam::HashList<Type, Key, Hash>::HashList(const label size)
+:
+    List<Tuple2<Key, Type>>(size, Tuple2<Key, Type>(nullKey, Type()))
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Type, class Key, class Hash>
+bool Foam::HashList<Type, Key, Hash>::insert(const Key& k, const Type& t)
+{
+    List<Tuple2<Key, Type>>& map = *this;
+
+    const label n = map.size();
+
+    const unsigned h = Hash()(k);
+
+    for (label i = 0; i < n; i ++)
+    {
+        const label hi = (h + i) % n;
+
+        if (map[hi].first() == nullKey)
+        {
+            map[hi] = Tuple2<Key, Type>(k, t);
+            return true;
+        }
+
+        if (map[hi].first() == k)
+        {
+            return false;
+        }
+    }
+
+    FatalErrorInFunction
+        << "Hash list is full"
+        << exit(FatalError);
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+template<class Type, class Key, class Hash>
+const Type& Foam::HashList<Type, Key, Hash>::operator[](const Key& k) const
+{
+    const List<Tuple2<Key, Type>>& map = *this;
+
+    const label n = map.size();
+
+    const unsigned h = Hash()(k);
+
+    for (label i = 0; i < n; i ++)
+    {
+        const label hi = (h + i) % n;
+
+        if (map[hi].first() == k)
+        {
+            return map[hi].second();
+        }
+    }
+
+    FatalErrorInFunction
+        << "Hash list does not contain key \"" << k << "\""
+        << exit(FatalError);
+
+    return NullObjectRef<Type>();
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/containers/HashTables/HashList/HashList.H

@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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::HashList
+
+Description
+    HashList class. Like HashTable, but much less dynamic memory-y. Should be
+    faster for small sets of non-dynamic primitive types (labels, edges,
+    points, etc...). It is also much less functional at present. There is no
+    re-sizing, so you have to make sure it is constructed sufficiently large to
+    hold all the data that will ever be inserted into it.
+
+SourceFiles
+    HashList.C
+
+\*---------------------------------------------------------------------------*/
+
+#include "List.H"
+#include "Tuple2.H"
+#include "word.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifndef HashList_H
+#define HashList_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class HashList Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class Key=word, class Hash=string::hash>
+class HashList
+:
+    private List<Tuple2<Key, Type>>
+{
+public:
+
+    // Public Static Member Data
+
+        //- Null key value for unset elements in the list
+        static const Key nullKey;
+
+
+    // Constructors
+
+        //- Construct given a size
+        HashList(const label size);
+
+
+    // Member Functions
+
+        //- Insert into the hash list. Return true if the value was newly
+        //  inserted, or false if it was already there.
+        bool insert(const Key& k, const Type& t);
+
+
+    // Member Operators
+
+        //- Retrieve from the hash list
+        const Type& operator[](const Key& k) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "HashList.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/plane/plane.C

@@ -309,13 +309,19 @@ Foam::point Foam::plane::aPoint() const
 
 Foam::point Foam::plane::nearestPoint(const point& p) const
 {
-    return p - normal_*((p - point_) & normal_);
+    return p - normal_*signedDistance(p);
 }
 
 
 Foam::scalar Foam::plane::distance(const point& p) const
 {
-    return mag((p - point_) & normal_);
+    return mag(signedDistance(p));
+}
+
+
+Foam::scalar Foam::plane::signedDistance(const point& p) const
+{
+    return (p - point_) & normal_;
 }
 
 

src/OpenFOAM/meshes/primitiveShapes/plane/plane.H

@@ -194,6 +194,9 @@ public:
         //- Return distance from the given point to the plane
         scalar distance(const point& p) const;
 
+        //- Return signed distance from the given point to the plane
+        scalar signedDistance(const point& p) const;
+
         //- Return cut coefficient for plane and line defined by
         //  origin and direction
         scalar normalIntersect(const point& pnt0, const vector& dir) const;

src/finiteVolume/Make/files

@@ -290,6 +290,7 @@ $(interpolation)/interpolation/interpolations.C
 
 $(interpolation)/interpolationCell/makeInterpolationCell.C
 $(interpolation)/interpolationCellPatchConstrained/makeInterpolationCellPatchConstrained.C
+$(interpolation)/interpolationVolPointInterpolation/makeInterpolationVolPointInterpolation.C
 $(interpolation)/interpolationCellPoint/cellPointWeight/cellPointWeight.C
 $(interpolation)/interpolationCellPoint/makeInterpolationCellPoint.C
 $(interpolation)/interpolationCellPointFace/makeInterpolationCellPointFace.C

src/finiteVolume/cfdTools/general/levelSet/levelSet.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "levelSet.H"
-#include "cut.H"
+#include "cutTriTet.H"
 #include "polyMeshTetDecomposition.H"
 #include "tetIndices.H"
 
@@ -38,6 +38,9 @@ Foam::tmp<Foam::scalarField> Foam::levelSetFraction
     const bool above
 )
 {
+    typedef cutTriTet::noOp noOp;
+    typedef cutTriTet::volumeOp volumeOp;
+
     tmp<scalarField> tResult(new scalarField(mesh.nCells(), Zero));
     scalarField& result = tResult.ref();
 
@@ -69,15 +72,15 @@ Foam::tmp<Foam::scalarField> Foam::levelSetFraction
                     levelP[triIs[2]]
                 };
 
-            v += cut::volumeOp()(tet);
+            v += volumeOp()(tet);
 
             if (above)
             {
-                r += tetCut(tet, level, cut::volumeOp(), cut::noOp());
+                r += tetCut(tet, level, volumeOp(), noOp());
             }
             else
             {
-                r += tetCut(tet, level, cut::noOp(), cut::volumeOp());
+                r += tetCut(tet, level, noOp(), volumeOp());
             }
         }
 
@@ -96,6 +99,9 @@ Foam::tmp<Foam::scalarField> Foam::levelSetFraction
     const bool above
 )
 {
+    typedef cutTriTet::noOp noOp;
+    typedef cutTriTet::areaMagOp areaMagOp;
+
     tmp<scalarField> tResult(new scalarField(patch.size(), 0));
     scalarField& result = tResult.ref();
 
@@ -124,15 +130,15 @@ Foam::tmp<Foam::scalarField> Foam::levelSetFraction
                     levelP[e[1]]
                 };
 
-            a += cut::areaMagOp()(tri);
+            a += areaMagOp()(tri);
 
             if (above)
             {
-                r += triCut(tri, level, cut::areaMagOp(), cut::noOp());
+                r += triCut(tri, level, areaMagOp(), noOp());
             }
             else
             {
-                r += triCut(tri, level, cut::noOp(), cut::areaMagOp());
+                r += triCut(tri, level, noOp(), areaMagOp());
             }
         }