ENH: surfaceBoooleanFeatures utility.

applications/utilities/surface/surfaceBooleanFeatures/Make/files

@@ -0,0 +1,3 @@
+surfaceBooleanFeatures.C
+
+EXE = $(FOAM_APPBIN)/surfaceBooleanFeatures

applications/utilities/surface/surfaceBooleanFeatures/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -ltriSurface \
+    -ledgeMesh \
+    -lmeshTools

applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C

@@ -0,0 +1,464 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     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/>.
+
+Application
+    surfaceBooleanFeatures
+
+Description
+
+    Generates the featureEdgeMesh for the interface between a boolean operation
+    on two surfaces.  Assumes that the orientation of the surfaces is correct:
+
+    + if the operation is union or intersection, that both surface's normals
+      (n) have the same orientation with respect to a point, i.e. surfaces and b
+      are orientated the same with respect to point x:
+
+    @verbatim
+       _______
+      |       |--> n
+      |    ___|___             x
+      |a  |   |   |--> n
+      |___|___|  b|
+          |       |
+          |_______|
+
+    @endverbatim
+
+    + if the operation is a subtraction, the surfaces should be oppositely
+    oriented with respect to a point, i.e. for (a - b), then b's orientation
+    should be such that x is "inside", and a's orientation such that x is
+    "outside"
+
+    @verbatim
+       _______
+      |       |--> n
+      |    ___|___             x
+      |a  |   |   |
+      |___|___|  b|
+          |  n <--|
+          |_______|
+
+    @endverbatim
+
+    When the operation is peformed - for union, all of the edges generates where
+    one surfaces cuts another are all "internal" for union, and "external" for
+    intersection, b - a and a - b.  This has been assumed, formal (dis)proof is
+    invited.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triSurface.H"
+#include "argList.H"
+#include "Time.H"
+#include "featureEdgeMesh.H"
+#include "triSurfaceSearch.H"
+#include "OFstream.H"
+#include "booleanSurface.H"
+#include "edgeIntersections.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Keep on shuffling surface points until no more degenerate intersections.
+// Moves both surfaces and updates set of edge cuts.
+bool intersectSurfaces
+(
+    triSurface& surf1,
+    edgeIntersections& edgeCuts1,
+    triSurface& surf2,
+    edgeIntersections& edgeCuts2
+)
+{
+    bool hasMoved1 = false;
+    bool hasMoved2 = false;
+
+    for (label iter = 0; iter < 10; iter++)
+    {
+        Info<< "Determining intersections of surf1 edges with surf2"
+            << " faces" << endl;
+
+        // Determine surface1 edge intersections. Allow surface to be moved.
+
+        // Number of iterations needed to resolve degenerates
+        label nIters1 = 0;
+        {
+            triSurfaceSearch querySurf2(surf2);
+
+            scalarField surf1PointTol
+            (
+                1e-3*edgeIntersections::minEdgeLength(surf1)
+            );
+
+            // Determine raw intersections
+            edgeCuts1 = edgeIntersections
+            (
+                surf1,
+                querySurf2,
+                surf1PointTol
+            );
+
+            // Shuffle a bit to resolve degenerate edge-face hits
+            {
+                pointField points1(surf1.points());
+
+                nIters1 =
+                    edgeCuts1.removeDegenerates
+                    (
+                        5,              // max iterations
+                        surf1,
+                        querySurf2,
+                        surf1PointTol,
+                        points1         // work array
+                    );
+
+                if (nIters1 != 0)
+                {
+                    // Update geometric quantities
+                    surf1.movePoints(points1);
+                    hasMoved1 = true;
+                }
+            }
+        }
+
+        Info<< "Determining intersections of surf2 edges with surf1"
+            << " faces" << endl;
+
+        label nIters2 = 0;
+        {
+            triSurfaceSearch querySurf1(surf1);
+
+            scalarField surf2PointTol
+            (
+                1e-3*edgeIntersections::minEdgeLength(surf2)
+            );
+
+            // Determine raw intersections
+            edgeCuts2 = edgeIntersections
+            (
+                surf2,
+                querySurf1,
+                surf2PointTol
+            );
+
+            // Shuffle a bit to resolve degenerate edge-face hits
+            {
+                pointField points2(surf2.points());
+
+                nIters2 =
+                    edgeCuts2.removeDegenerates
+                    (
+                        5,              // max iterations
+                        surf2,
+                        querySurf1,
+                        surf2PointTol,
+                        points2         // work array
+                    );
+
+                if (nIters2 != 0)
+                {
+                    // Update geometric quantities
+                    surf2.movePoints(points2);
+                    hasMoved2 = true;
+                }
+            }
+        }
+
+        if (nIters1 == 0 && nIters2 == 0)
+        {
+            Info<< "** Resolved all intersections to be proper edge-face pierce"
+                << endl;
+            break;
+        }
+    }
+
+    if (hasMoved1)
+    {
+        fileName newFile("surf1.obj");
+        Info<< "Surface 1 has been moved. Writing to " << newFile
+            << endl;
+        surf1.write(newFile);
+    }
+
+    if (hasMoved2)
+    {
+        fileName newFile("surf2.obj");
+        Info<< "Surface 2 has been moved. Writing to " << newFile
+            << endl;
+        surf2.write(newFile);
+    }
+
+    return hasMoved1 || hasMoved2;
+}
+
+
+int main(int argc, char *argv[])
+{
+    argList::validArgs.clear();
+    argList::noParallel();
+    argList::validArgs.append("action");
+    argList::validArgs.append("surface file");
+    argList::validArgs.append("surface file");
+    argList::validArgs.append("output file");
+
+    argList::addBoolOption
+    (
+        "perturb",
+        "Perturb surface points to escape degenerate intersections"
+    );
+
+    argList::addBoolOption
+    (
+        "invertedSpace",
+        "do the surfaces have inverted space orientation, "
+        "i.e. a point at infinity is considered inside. "
+        "This is only sensible for union and intersection."
+    );
+
+    #   include "setRootCase.H"
+    #   include "createTime.H"
+
+    word action(args.args()[1]);
+
+    HashTable<booleanSurface::booleanOpType> validActions;
+    validActions.insert("intersection", booleanSurface::INTERSECTION);
+    validActions.insert("union", booleanSurface::UNION);
+    validActions.insert("difference", booleanSurface::DIFFERENCE);
+
+    if (!validActions.found(action))
+    {
+        FatalErrorIn(args.executable())
+            << "Unsupported action " << action << endl
+            << "Supported actions:" << validActions.toc() << exit(FatalError);
+    }
+
+    Info<< "Reading surface 1 .." << endl;
+    fileName surf1Name(args[2]);
+    triSurface surf1(surf1Name);
+
+    Info<< "Surface 1 statistics:" << endl;
+    surf1.writeStats(Info);
+    Info<< endl;
+
+    Info<< "Reading surface 2 .." << endl;
+    fileName surf2Name(args[3]);
+    triSurface surf2(surf2Name);
+
+    Info<< "Surface 2 statistics:" << endl;
+    surf2.writeStats(Info);
+    Info<< endl;
+
+    fileName outFileName(args[4]);
+    Info<< "Output file name " << outFileName << endl;
+
+    Info<< "Intersecting surface 1 and 2" << endl;
+
+    edgeIntersections edge1Cuts;
+    edgeIntersections edge2Cuts;
+
+    bool invertedSpace = args.optionFound("invertedSpace");
+
+    if (invertedSpace && validActions[action] == booleanSurface::DIFFERENCE)
+    {
+        FatalErrorIn(args.executable())
+            << "Inverted space only makes sense for union or intersection."
+            << exit(FatalError);
+    }
+
+    if (args.optionFound("perturb"))
+    {
+        intersectSurfaces
+        (
+            surf1,
+            edge1Cuts,
+            surf2,
+            edge2Cuts
+        );
+    }
+    else
+    {
+        triSurfaceSearch querySurf2(surf2);
+
+        Info<< "Determining intersections of surf1 edges with surf2 faces"
+            << endl;
+
+        edge1Cuts = edgeIntersections
+        (
+            surf1,
+            querySurf2,
+            1e-3*edgeIntersections::minEdgeLength(surf1)
+        );
+
+        triSurfaceSearch querySurf1(surf1);
+
+        Info<< "Determining intersections of surf2 edges with surf1 faces"
+            << endl;
+
+        edge2Cuts = edgeIntersections
+        (
+            surf2,
+            querySurf1,
+            1e-3*edgeIntersections::minEdgeLength(surf2)
+        );
+    }
+
+    // Determine intersection edges
+    surfaceIntersection inter(surf1, edge1Cuts, surf2, edge2Cuts);
+
+    OFstream intFile(outFileName);
+
+    fileName sFeatFileName =
+        surf1Name.lessExt().name()
+      + "_"
+      + surf2Name.lessExt().name()
+      + "_"
+      + action;
+
+    label nFeatEds = inter.cutEdges().size();
+
+    vectorField normals(2*nFeatEds, vector::zero);
+    vectorField edgeDirections(nFeatEds, vector::zero);
+    labelListList edgeNormals(nFeatEds, labelList(2, -1));
+
+    triSurfaceSearch querySurf1(surf1);
+    triSurfaceSearch querySurf2(surf2);
+
+    OFstream normalFile(sFeatFileName + "_normals.obj");
+
+    forAll(inter.cutEdges(), i)
+    {
+        const edge& fE(inter.cutEdges()[i]);
+
+        point fEC = fE.centre(inter.cutPoints());
+
+        pointIndexHit nearest1 = querySurf1.tree().findNearest(fEC, sqr(GREAT));
+        pointIndexHit nearest2 = querySurf2.tree().findNearest(fEC, sqr(GREAT));
+
+        normals[2*i] = surf1.faceNormals()[nearest1.index()];
+        normals[2*i + 1] = surf2.faceNormals()[nearest2.index()];
+
+        edgeNormals[i][0] = 2*i;
+        edgeNormals[i][1] = 2*i + 1;
+
+        edgeDirections[i] = fE.vec(inter.cutPoints());
+
+        {
+            scalar scale = 3*fE.mag(inter.cutPoints());
+
+            meshTools::writeOBJ(normalFile, inter.cutPoints()[fE.start()]);
+            meshTools::writeOBJ(normalFile, inter.cutPoints()[fE.end()]);
+
+            normalFile<< "l " << (5*i) + 1 << " " << (5*i) + 2<< endl;
+
+            meshTools::writeOBJ(normalFile, fEC);
+            meshTools::writeOBJ(normalFile, fEC + scale*normals[2*i]);
+            meshTools::writeOBJ(normalFile, fEC + scale*normals[2*i + 1]);
+
+            normalFile<< "l " << (5*i) + 3 << " " << (5*i) + 4 << endl;
+            normalFile<< "l " << (5*i) + 3 << " " << (5*i) + 5 << endl;
+        }
+    }
+
+    label internalStart = -1;
+
+    if (validActions[action] == booleanSurface::UNION)
+    {
+        if (!invertedSpace)
+        {
+            // All edges are internal
+            internalStart = 0;
+        }
+        else
+        {
+            // All edges are external
+            internalStart = nFeatEds;
+        }
+    }
+    else if (validActions[action] == booleanSurface::INTERSECTION)
+    {
+        if (!invertedSpace)
+        {
+            // All edges are external
+            internalStart = nFeatEds;
+        }
+        else
+        {
+            // All edges are internal
+            internalStart = 0;
+        }
+    }
+    else if (validActions[action] == booleanSurface::DIFFERENCE)
+    {
+        // All edges are external
+        internalStart = nFeatEds;
+    }
+    else
+    {
+        FatalErrorIn(args.executable())
+            << "Unsupported booleanSurface:booleanOpType and space "
+            << action << " " << invertedSpace
+            << abort(FatalError);
+    }
+
+    // There are no feature points supported by surfaceIntersection
+    // Flat, open or multiple edges are assumed to be impossible
+    // Region edges are not explicitly supported by surfaceIntersection
+
+    featureEdgeMesh feMesh
+    (
+        IOobject
+        (
+            sFeatFileName + ".featureEdgeMesh",
+            runTime.constant(),
+            "featureEdgeMesh",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        inter.cutPoints(),
+        inter.cutEdges(),
+        0,                  // concaveStart,
+        0,                  // mixedStart,
+        0,                  // nonFeatureStart,
+        internalStart,      // internalStart,
+        nFeatEds,           // flatStart,
+        nFeatEds,           // openStart,
+        nFeatEds,           // multipleStart,
+        normals,
+        edgeDirections,
+        edgeNormals,
+        labelListList(0),   // featurePointNormals,
+        labelList(0)        // regionEdges
+    );
+
+    feMesh.write();
+
+    feMesh.writeObj(sFeatFileName);
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/meshTools/triSurface/booleanOps/booleanSurface/booleanSurface.C

@@ -997,9 +997,9 @@ Foam::booleanSurface::booleanSurface
 
 
     // Depending on operation include certain faces.
-    //  AND: faces on inside of 1 and of 2
+    //  INTERSECTION: faces on inside of 1 and of 2
-    //  OR: faces on outside of 1 and of 2
+    //  UNION: faces on outside of 1 and of 2
-    //  XOR: faces on outside of 1 and inside of 2
+    //  DIFFERENCE: faces on outside of 1 and inside of 2
 
     boolList include(combinedSurf.size(), false);
 
@@ -1018,30 +1018,30 @@ Foam::booleanSurface::booleanSurface
         }
         else if (side[faceI] == OUTSIDE)
         {
-            if (booleanOp == booleanSurface::OR)
+            if (booleanOp == booleanSurface::UNION)
             {
                 include[faceI] = true;
             }
-            else if (booleanOp == booleanSurface::AND)
+            else if (booleanOp == booleanSurface::INTERSECTION)
             {
                 include[faceI] = false;
             }
-            else    // xor
+            else    // difference
             {
                 include[faceI] = (faceI < cutSurf1.size()); // face from surf1
             }
         }
         else    // inside
         {
-            if (booleanOp == booleanSurface::OR)
+            if (booleanOp == booleanSurface::UNION)
             {
                 include[faceI] = false;
             }
-            else if (booleanOp == booleanSurface::AND)
+            else if (booleanOp == booleanSurface::INTERSECTION)
             {
                 include[faceI] = true;
             }
-            else    // xor
+            else    // difference
             {
                 include[faceI] = (faceI >= cutSurf1.size()); // face from surf2
             }

src/meshTools/triSurface/booleanOps/booleanSurface/booleanSurface.H

@@ -30,7 +30,7 @@ Description
     Called 'boolean' since the volume of resulting surface will encompass
     the volumes of the original surface according to some boolean operation:
     - all which is in surface1 AND in surface2 (intersection)
-    - all which is in surface1 AND NOT in surface2 ('chop' out surface2)
+    - all which is in surface1 AND NOT in surface2 (surface1 minus surface2)
     - all which is in surface1 OR in surface2 (union)
 
     Algorithm:
@@ -163,11 +163,11 @@ public:
         //- Enumeration listing the possible volume operator types
         enum booleanOpType
         {
-            OR,     // Union of volumes
+            UNION,        // Union of volumes
-            AND,    // Intersection of volumes
+            INTERSECTION, // Intersection of volumes
-            XOR,    // Difference of volumes
+            DIFFERENCE,   // Difference of volumes
-            ALL     // Special: does not subset combined surface. (Produces
+            ALL           // Special: does not subset combined
-                    // multiply connected surface)
+                          // surface. (Produces multiply connected surface)
         };