ENH: region-wise self intersection for surfaceFeatureExtract (issue #450)

applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C

@@ -83,7 +83,7 @@ Description
 #include "edgeIntersections.H"
 #include "meshTools.H"
 #include "DynamicField.H"
-
+#include "Enum.H"
 
 #ifndef NO_CGAL
 
@@ -1514,8 +1514,8 @@ int main(int argc, char *argv[])
 {
     argList::noParallel();
     argList::validArgs.append("action");
-    argList::validArgs.append("surface file");
+    argList::validArgs.append("surfaceFile1");
-    argList::validArgs.append("surface file");
+    argList::validArgs.append("surfaceFile2");
 
     argList::addBoolOption
     (
@@ -1553,24 +1553,30 @@ int main(int argc, char *argv[])
         " 'mixed' (keep all)"
     );
 
+    argList::addNote
+    (
+        "Valid actions: \"intersection\", \"union\", \"difference\""
+    );
+
 
     #include "setRootCase.H"
     #include "createTime.H"
 
     const word action(args[1]);
 
-    const HashTable<booleanSurface::booleanOpType> validActions
+    const Enum<booleanSurface::booleanOpType> validActions
     {
-        {"intersection", booleanSurface::INTERSECTION},
+        { booleanSurface::INTERSECTION, "intersection" },
-        {"union", booleanSurface::UNION},
+        { booleanSurface::UNION, "union" },
-        {"difference", booleanSurface::DIFFERENCE}
+        { booleanSurface::DIFFERENCE, "difference" }
     };
 
-    if (!validActions.found(action))
+    if (!validActions.hasEnum(action))
     {
         FatalErrorInFunction
             << "Unsupported action " << action << endl
-            << "Supported actions:" << validActions.toc() << abort(FatalError);
+            << "Supported actions:" << validActions << nl
+            << abort(FatalError);
     }
 
 

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -149,8 +149,21 @@ int main(int argc, char *argv[])
         }
         Info<< "Output             : " << outputName << nl;
 
+        triSurfaceLoader::loadingOption loadingOption =
+            triSurfaceLoader::loadingOptionNames.lookupOrDefault
+            (
+                "loadingOption",
+                surfaceDict,
+                triSurfaceLoader::loadingOption::OFFSET_REGION
+            );
+
+        Info<<"loading with "
+            << triSurfaceLoader::loadingOptionNames[loadingOption]
+            << endl;
+
+
         // Load a single file, or load and combine multiple selected files
-        autoPtr<triSurface> surfPtr = loader.load();
+        autoPtr<triSurface> surfPtr = loader.load(loadingOption);
         if (!surfPtr.valid() || surfPtr().empty())
         {
             FatalErrorInFunction
@@ -390,14 +403,21 @@ int main(int argc, char *argv[])
             feMesh.add(addFeMesh);
         }
 
-        if (surfaceDict.lookupOrDefault<bool>("selfIntersection", false))
+        const surfaceIntersection::intersectionType selfIntersect =
-        {
+            surfaceIntersection::selfIntersectionNames.lookupOrDefault
-            // TODO: perturbance tolerance?
+            (
+                "intersectionMethod",
+                surfaceDict,
+                surfaceIntersection::NONE
+            );
 
+        if (selfIntersect != surfaceIntersection::NONE)
+        {
             triSurfaceSearch query(surf);
             surfaceIntersection intersect(query, surfaceDict);
 
-            intersect.mergePoints(5*SMALL);
+            // Remove rounding noise - could make adjustable
+            intersect.mergePoints(10*SMALL);
 
             labelPair sizeInfo
             (
@@ -413,14 +433,15 @@ int main(int argc, char *argv[])
                     intersect.cutEdges()
                 );
 
-                addMesh.mergePoints(5*SMALL);
                 feMesh.add(addMesh);
 
                 sizeInfo[0] = addMesh.points().size();
                 sizeInfo[1] = addMesh.edges().size();
             }
             Info<< nl
-                << "Self intersection:" << nl
+                << "intersection: "
+                << surfaceIntersection::selfIntersectionNames[selfIntersect]
+                << nl
                 << "    points : " << sizeInfo[0] << nl
                 << "    edges  : " << sizeInfo[1] << nl;
         }

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict

@@ -16,7 +16,7 @@ FoamFile
 
 surface1.stl
 {
-    // How to obtain raw features (extractFromFile | extractFromSurface | none)
+    // How to obtain raw features (none | extractFromFile | extractFromSurface)
     extractionMethod    extractFromSurface;
 
     // Mark edges whose adjacent surface normals are at an angle less
@@ -36,8 +36,8 @@ surface1.stl
         geometricTestOnly  yes;
     } */
 
-    // Generate additional features from self-intersect
+    // Generate additional intersection features (none | self | region)
-    selfIntersection    false;
+    intersectionMethod  none;
 
     // Tolerance for surface intersections
     tolerance           1e-3;
@@ -51,7 +51,7 @@ surface1.stl
 
 surface2.nas
 {
-    // How to obtain raw features (extractFromFile | extractFromSurface | none)
+    // How to obtain raw features (none | extractFromFile | extractFromSurface)
     extractionMethod    extractFromFile;
 
     extractFromFileCoeffs
@@ -114,8 +114,8 @@ surface2.nas
     // Out put the closeness of surface elements to other surface elements.
     closeness           no;
 
-    // Generate additional features from self-intersect
+    // Generate additional intersection features (none | self | region)
-    selfIntersection    false;
+    intersectionMethod  none;
 
     // Tolerance for surface intersections
     tolerance           1e-3;
@@ -148,8 +148,8 @@ dummyName
     // Base output name (optional)
     // output              surfaces;
 
-    // Generate additional features from self-intersect
+    // Generate additional intersection features (none | self | region)
-    selfIntersection    true;
+    intersectionMethod  self;
 
     // Tolerance for surface intersections
     tolerance           1e-3;
@@ -183,8 +183,8 @@ surfaces
     // Base output name (optional)
     // output              surfaces;
 
-    // Generate additional features from self-intersect
+    // Generate additional intersection features (none | self | region)
-    selfIntersection    true;
+    intersectionMethod  self;
 
     // Tolerance for surface intersections
     tolerance           1e-3;
@@ -193,8 +193,7 @@ surfaces
     noneCoeffs
     {
         includedAngle   0;
-    }
+    } */
-    */
 
     // Write options
 

src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C

@@ -27,6 +27,7 @@ License
 #include "triSurfaceSearch.H"
 #include "OBJstream.H"
 #include "labelPairHashes.H"
+#include "PackedBoolList.H"
 #include "triSurface.H"
 #include "pointIndexHit.H"
 #include "mergePoints.H"
@@ -40,6 +41,14 @@ namespace Foam
 defineTypeNameAndDebug(surfaceIntersection, 0);
 }
 
+const Foam::Enum<Foam::surfaceIntersection::intersectionType>
+Foam::surfaceIntersection::selfIntersectionNames
+{
+    { intersectionType::SELF, "self" },
+    { intersectionType::SELF_REGION, "region" },
+    { intersectionType::NONE, "none" },
+};
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -54,7 +63,7 @@ void Foam::surfaceIntersection::setOptions(const dictionary& dict)
 
 void Foam::surfaceIntersection::storeIntersection
 (
-    const enum originatingType cutFrom,
+    const enum intersectionType cutFrom,
     const labelList& facesA,
     const label faceB,
     const UList<point>& allCutPoints,
@@ -85,6 +94,7 @@ void Foam::surfaceIntersection::storeIntersection
                 break;
             }
             case surfaceIntersection::SELF:
+            case surfaceIntersection::SELF_REGION:
             {
                 // Lookup should be commutativity - use sorted order
                 if (faceA < faceB)
@@ -99,6 +109,10 @@ void Foam::surfaceIntersection::storeIntersection
                 }
                 break;
             }
+
+            case surfaceIntersection::NONE:
+                return;
+                break;
         }
 
 
@@ -245,7 +259,7 @@ void Foam::surfaceIntersection::classifyHit
     const triSurface& surf1,
     const scalarField& surf1PointTol,
     const triSurface& surf2,
-    const enum originatingType cutFrom,
+    const enum intersectionType cutFrom,
     const label edgeI,
     const pointIndexHit& pHit,
 
@@ -306,7 +320,11 @@ void Foam::surfaceIntersection::classifyHit
             // For self-intersection, we have tolerances for each point
             // (surf2 is actually surf1) so we shift the hit to coincide
             // identically.
-            if (cutFrom == surfaceIntersection::SELF)
+            if
+            (
+                cutFrom == surfaceIntersection::SELF
+             || cutFrom == surfaceIntersection::SELF_REGION
+            )
             {
                 const point& nearPt = surf1Pts[nearVert];
 
@@ -392,7 +410,15 @@ void Foam::surfaceIntersection::classifyHit
             // >0: store point/edge-cut. Attempt to create new edge.
             // <0: store point/edge-cut only
             int handling = (allowEdgeHits_ ? 1 : 0);
-            if (allowEdgeHits_ && cutFrom == surfaceIntersection::SELF)
+            if
+            (
+                allowEdgeHits_
+             &&
+                (
+                    cutFrom == surfaceIntersection::SELF
+                 || cutFrom == surfaceIntersection::SELF_REGION
+                )
+            )
             {
                 // The edge-edge intersection is hashed as an 'edge' to
                 // exploit the commutative lookup.
@@ -513,12 +539,18 @@ void Foam::surfaceIntersection::classifyHit
             switch (cutFrom)
             {
                 case surfaceIntersection::FIRST:
+                {
                     handling = 1;
                     break;
+                }
                 case surfaceIntersection::SECOND:
+                {
                     handling = -1;
                     break;
+                }
                 case surfaceIntersection::SELF:
+                case surfaceIntersection::SELF_REGION:
+                {
                     // The edge-edge intersection is hashed as an 'edge' to
                     // exploit the commutative lookup.
                     // Ie, only do the cut once
@@ -563,6 +595,12 @@ void Foam::surfaceIntersection::classifyHit
                             }
                         }
                     }
+
+                    break;
+                }
+
+                case surfaceIntersection::NONE:
+                    return;
                     break;
             }
 
@@ -744,7 +782,7 @@ void Foam::surfaceIntersection::doCutEdges
 (
     const triSurface& surf1,
     const triSurfaceSearch& querySurf2,
-    const enum originatingType cutFrom,
+    const enum intersectionType cutFrom,
 
     DynamicList<point>& allCutPoints,
     DynamicList<edge>& allCutEdges,
@@ -766,13 +804,21 @@ void Foam::surfaceIntersection::doCutEdges
     const indexedOctree<treeDataPrimitivePatch<triSurface>>& searchTree
         = querySurf2.tree();
 
-    if (cutFrom == surfaceIntersection::SELF)
+    if
+    (
+        cutFrom == surfaceIntersection::SELF
+     || cutFrom == surfaceIntersection::SELF_REGION
+    )
     {
         // An edge may intersect multiple faces
         // - mask out faces that have already been hit before trying again
         // - never intersect with faces attached to the edge itself
         DynamicList<label> maskFaces(32);
 
+        // Optionally prevent intersection within a single region.
+        // Like self-intersect, but only if regions are different
+        PackedBoolList maskRegions(32);
+
         treeDataTriSurface::findAllIntersectOp
             allIntersectOp(searchTree, maskFaces);
 
@@ -787,6 +833,15 @@ void Foam::surfaceIntersection::doCutEdges
             const point ptEnd =
                 surf1Pts[e.end()]   + 0.5*surf1PointTol[e.end()]*edgeVec;
 
+            maskRegions.clear();
+            if (cutFrom == surfaceIntersection::SELF_REGION)
+            {
+                for (auto& facei : surf1.edgeFaces()[edgeI])
+                {
+                    maskRegions.set(surf1[facei].region());
+                }
+            }
+
             // Never intersect with faces attached directly to the edge itself,
             // nor with faces attached to its end points. This mask contains
             // some duplicates, but filtering them out is less efficient.
@@ -809,6 +864,11 @@ void Foam::surfaceIntersection::doCutEdges
 
                 maskFaces.append(pHit.index());
 
+                if (maskRegions[surf1[pHit.index()].region()])
+                {
+                    continue;
+                }
+
                 classifyHit
                 (
                     surf1,
@@ -1133,6 +1193,27 @@ Foam::surfaceIntersection::surfaceIntersection
 {
     setOptions(dict);
 
+    const intersectionType cutFrom = selfIntersectionNames.lookupOrDefault
+    (
+        "intersectionMethod",
+        dict,
+        intersectionType::SELF
+    );
+
+    if (cutFrom == intersectionType::NONE)
+    {
+        if (debug)
+        {
+            Pout<< "Skipping self-intersection (selected: none)" << endl;
+        }
+
+        // Temporaries
+        facePairToEdge_.clear();
+        facePairToEdgeId_.clear();
+
+        return;
+    }
+
     const triSurface& surf1 = query1.surface();
 
     //
@@ -1154,7 +1235,7 @@ Foam::surfaceIntersection::surfaceIntersection
     (
         surf1,
         query1,
-        surfaceIntersection::SELF,
+        cutFrom,
         allCutPoints,
         allCutEdges,
         edgeCuts1
@@ -1402,17 +1483,16 @@ void Foam::surfaceIntersection::mergePoints(const scalar mergeDist)
     pointField newPoints;
     labelList pointMap;
 
-    const bool hasMerged = Foam::mergePoints
+    const label nMerged = Foam::mergePoints
     (
         cutPoints_,
         mergeDist,
         false,
         pointMap,
-        newPoints,
+        newPoints
-        vector::zero
     );
 
-    if (hasMerged)
+    if (nMerged)
     {
         cutPoints_.transfer(newPoints);
 

src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H

@@ -52,6 +52,7 @@ Description
         allowEdgeHits  | Edge-end cuts another edge     | bool   | true
         snap           | Snap near end-points           | bool   | true
         warnDegenerate | Number of warnings about degenerate edges | label | 0
+        intersectionMethod | Control for self intersection | (self,region,none)
     \endtable
 
 SourceFiles
@@ -71,6 +72,7 @@ SourceFiles
 #include "labelPairHashes.H"
 #include "pointIndexHit.H"
 #include "typeInfo.H"
+#include "Enum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -96,15 +98,25 @@ typedef EdgeMap<labelPairHashSet> edgelabelPairHashLookup;
 
 class surfaceIntersection
 {
-    // Private data
+public:
 
-        //- Surface intersection types for classify, doCutEdges
+    //- Surface intersection types for classify, doCutEdges
-        enum originatingType
+    enum intersectionType
-        {
+    {
-            FIRST,      //!< First surface
+        FIRST,          //!< First surface
-            SECOND,     //!< Second surface
+        SECOND,         //!< Second surface
-            SELF        //!< Self-intersection
+        SELF,           //!< Self-intersection
-        };
+        SELF_REGION,    //!< Self-intersection, region-wise only
+        NONE            //!< None = invalid (for input only)
+    };
+
+    //- The user-selectable self-intersection enumeration names
+    static const Enum<intersectionType> selfIntersectionNames;
+
+
+private:
+
+    // Private data
 
         //- Tolerance for intersections
         scalar tolerance_;
@@ -204,7 +216,7 @@ class surfaceIntersection
         //  Updates facePairToEdge_ and facePairToEdgeId_ (on the second hit)
         void storeIntersection
         (
-            const enum originatingType cutFrom,
+            const enum intersectionType cutFrom,
             const labelList& facesA,
             const label faceB,
             const UList<point>& allCutPoints,
@@ -219,7 +231,7 @@ class surfaceIntersection
             const triSurface& surf1,
             const scalarField& surf1PointTol,
             const triSurface& surf2,
-            const enum originatingType cutFrom,
+            const enum intersectionType cutFrom,
             const label edgeI,
             const pointIndexHit& pHit,
 
@@ -233,7 +245,7 @@ class surfaceIntersection
         (
             const triSurface& surf1,
             const triSurfaceSearch& querySurf2,
-            const enum originatingType cutFrom,
+            const enum intersectionType cutFrom,
 
             DynamicList<point>& allCutPoints,
             DynamicList<edge>& allCutEdges,