Merge branch 'cvm'

applications/utilities/mesh/generation/cvMesh/Make/files

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

applications/utilities/mesh/generation/cvMesh/Make/options

@@ -0,0 +1,26 @@
+EXE_DEBUG = -DFULLDEBUG -g -O0
+EXE_FROUNDING_MATH = -frounding-math
+EXE_NDEBUG = -DNDEBUG
+
+include $(GENERAL_RULES)/CGAL
+
+EXE_INC = \
+    ${EXE_FROUNDING_MATH} \
+    ${EXE_NDEBUG} \
+    ${CGAL_INC} \
+    -I$(LIB_SRC)/mesh/conformalVoronoiMesh/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude
+
+EXE_LIBS = \
+    $(CGAL_LIBS) \
+    -lconformalVoronoiMesh \
+    -lmeshTools \
+    -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
+    -ledgeMesh \
+    -ltriSurface \
+    -ldynamicMesh

applications/utilities/mesh/generation/cvMesh/cvMesh.C

@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
+  \\    /    O peration      |
+   \\  /     A nd            | Copyright (C) 2008-2011 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
+    cvMesh
+
+Description
+    Conformal Voronoi automatic mesh generator
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "conformalVoronoiMesh.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+    runTime.functionObjects().off();
+
+    IOdictionary cvMeshDict
+    (
+        IOobject
+        (
+            "cvMeshDict",
+            runTime.system(),
+            runTime,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    conformalVoronoiMesh mesh(runTime, cvMeshDict);
+
+    while (runTime.loop())
+    {
+        Info<< nl << "Time = " << runTime.timeName() << endl;
+
+        mesh.move();
+
+        Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    mesh.writeMesh(runTime.constant());
+
+    Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+    Info<< nl << "End" << nl << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/cvMesh/cvMeshDict

@@ -0,0 +1,496 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+    class           dictionary;
+    object          cvMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+/*
+
+Control dictionary for cvMesh - polyhedral mesh generator.
+
+cvMesh phases:
+1. fill volume with initial points (initialPoints subdictionary). An option
+   is to reread from previous set of points.
+
+2. internal point motion (motionControl subdictionary)
+
+3. every once in a while add point duplets/triplets to conform to
+   surfaces and features (surfaceConformation subdictionary)
+
+4. back to 2
+
+5. construct polyMesh.
+    - filter (polyMeshFiltering subdictionary)
+    - check (meshQualityControls subdictionary) and undo filtering
+
+
+See also cvControls.H in the conformalVoronoiMesh library
+
+*/
+
+
+// Important:
+// ----------
+// Any scalar with a name <name>Coeff specifies a value that will be implemented
+// as a faction of the local target cell size
+// Any scalar with a name <name>Size specifies an absolute size.
+
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces need to be (almost) closed - use closedTriSurfaceMesh
+// if they are not topologically closed. Surfaces need to be oriented so
+// the space to be meshed is always on the inside of all surfaces. Use e.g.
+// surfaceOrient.
+geometry
+{
+    // Surface to conform to
+    flange.obj
+    {
+        type triSurfaceMesh;
+    }
+
+    // Surface used for cell size control
+    ref7.stl
+    {
+        name ref7;
+        type triSurfaceMesh;
+    }
+
+    // Surface used for cell size control
+    tunnel_APPROACH_INLET.obj
+    {
+        type triSurfaceMesh;
+    }
+}
+
+
+// Controls for conforming to the surfaces.
+surfaceConformation
+{
+    // A point inside surfaces that is inside mesh.
+    locationInMesh (0 0 0);
+
+    // How far apart are point-duplets generated. Balance this between
+    // - very low distance: little chance of interference from other
+    //   surfaces
+    // - largish distance: less non-orthogonality in final cell
+    //   (circumcentre far away from centroid)
+    pointPairDistanceCoeff      0.1;
+
+    // Mixed feature points - connected to both inside and outside edges.
+    // Recreated by inserting triplets of points to recreate a single edge.
+    // Done for all edges emanating from point. triplets of points get inserted
+    // mixedFeaturePointPPDistanceCoeff distance away from feature point.
+    // Set to a negative number to disable.
+    mixedFeaturePointPPDistanceCoeff    5.0; //-1;
+
+    // Distance to a feature point within which surface and edge
+    // conformation points are excluded - fraction of the local target
+    // cell size
+    featurePointExclusionDistanceCoeff  0.4;
+
+    // Distance to an existing feature edge conformation location
+    // within which other edge conformation location are excluded -
+    // fraction of the local target cell size
+    featureEdgeExclusionDistanceCoeff   0.2;
+
+    // Optimisation: do not check for surface intersection (of dual edges)
+    // for points near to surface.
+    surfaceSearchDistanceCoeff  2.5;
+
+    // Maximum allowable protrusion through the surface before
+    // conformation points are added - fraction of the local target
+    // cell size. These small protusions are (hopefully) done by mesh filtering
+    // instead.
+    maxSurfaceProtrusionCoeff   0.1;
+
+    // If feature edge with large angle (so more than 125 degrees) introduce
+    // additional points to create two half angled cells (= mitering).
+    maxQuadAngle 125;
+
+    // Frequency to redo surface conformation (expensive).
+    surfaceConformationRebuildFrequency 10;
+
+    // Initial and intermediate controls
+    coarseConformationControls
+    {
+        // Initial conformation
+        initial
+        {
+            // We've got a point poking through the surface. Don't do any
+            // surface conformation if near feature edge (since feature edge
+            // conformation should have priority)
+
+            // distance to search for near feature edges
+            edgeSearchDistCoeff         1.1;
+
+            // Proximity to a feature edge where a surface hit is
+            //  not created, only the edge conformation is created
+            // - fraction of the local target cell size. Coarse
+            // conformation, initial protrusion tests.
+            surfacePtReplaceDistCoeff   0.5;
+        }
+
+        // Same for iterations
+        iteration
+        {
+            edgeSearchDistCoeff         1.25;
+            surfacePtReplaceDistCoeff   0.7;
+        }
+
+        // Stop either at maxIterations or if the number of surface pokes
+        // is very small (iterationToInitialHitRatioLimit * initial number)
+        // Note: perhaps iterationToInitialHitRatioLimit should be absolute
+        // count?
+        maxIterations 15;
+
+        iterationToInitialHitRatioLimit 0.001;
+    }
+
+    // Final (at endTime) controls
+    fineConformationControls
+    {
+        initial
+        {
+            edgeSearchDistCoeff         1.1;
+            surfacePtReplaceDistCoeff   0.5;
+        }
+
+        iteration
+        {
+            edgeSearchDistCoeff         1.25;
+            surfacePtReplaceDistCoeff   0.7;
+        }
+
+        maxIterations 15;
+
+        iterationToInitialHitRatioLimit 0.001;
+    }
+
+    // Geometry to mesh to
+    geometryToConformTo
+    {
+        flange.obj
+        {
+            featureMethod extendedFeatureEdgeMesh; // or none;
+            extendedFeatureEdgeMesh "flange.extendedFeatureEdgeMesh";
+        }
+    }
+
+    // Additional features.
+    additionalFeatures
+    {
+        //line
+        //{
+        //    featureMethod extendedFeatureEdgeMesh; // or none;
+        //    extendedFeatureEdgeMesh "line.extendedFeatureEdgeMesh";
+        //}
+    }
+}
+
+
+// Controls for seeding initial points and general control of the target
+// cell size (used everywhere)
+initialPoints
+{
+    // Do not place point closer than minimumSurfaceDistanceCoeff
+    // to the surface. Is fraction of local target cell size (see below)
+    minimumSurfaceDistanceCoeff 0.55;
+
+    initialPointsMethod             autoDensity;
+    // initialPointsMethod          uniformGrid;
+    // initialPointsMethod          bodyCentredCubic;
+    // initialPointsMethod          pointFile;
+
+
+    // Take boundbox of all geometry. Sample with this box. If too much
+    // samples in box (due to target cell size) split box.
+    autoDensityDetails
+    {
+        // Initial number of refinement levels. Needs to be enough to pick
+        // up features due to size ratio. If not enough it will take longer
+        // to determine point seeding.
+        minLevels 4;
+        // Split box if ratio of min to max cell size larger than maxSizeRatio
+        maxSizeRatio 5.0;
+        // Per box sample 5x5x5 internally
+        sampleResolution 5;
+        // Additionally per face of the box sample 5x5
+        surfaceSampleResolution 5;
+    }
+
+    uniformGridDetails
+    {
+        // Absolute cell size.
+        initialCellSize         0.0015;
+        randomiseInitialGrid    yes;
+        randomPerturbationCoeff 0.02;
+    }
+
+    bodyCentredCubicDetails
+    {
+        initialCellSize         0.0015;
+        randomiseInitialGrid    no;
+        randomPerturbationCoeff 0.1;
+    }
+
+    pointFileDetails
+    {
+        // Reads points from file. Still rejects points that are too
+        // close to the surface (minimumSurfaceDistanceCoeff) or on the
+        // wrong side of the surfaces.
+        pointFile              "constant/internalDelaunayVertices";
+    }
+}
+
+
+// Control size of voronoi cells i.e. distance between points. This
+// determines the target cell size which is used everywhere.
+// It determines the cell size given a location. It then uses all
+// the rules
+// - defaultCellSize
+// - cellSizeControlGeometry
+// to determine target cell size. Rule with highest priority wins. If same
+// priority smallest cell size wins.
+motionControl
+{
+    // Absolute cell size of back ground mesh. This is the maximum cell size.
+    defaultCellSize     0.00075;
+
+    // Assign a priority to all requests for cell sizes, the highest overrules.
+    defaultPriority     0;
+
+    cellSizeControlGeometry
+    {
+        ref7_outside
+        {
+            // optional name of geometry
+            surface             ref7;
+            priority            1;
+            mode                outside;//inside/bothSides
+            cellSizeFunction    linearDistance;
+
+            // cellSizeFunctions:
+            //  uniform : uniform size
+            //  uniformDistance : fixed size for all within distance
+            //  linearSpatial : grading in specified direction only
+            //  linearDistance : vary linearly as distance to surface
+            //  surfaceOffsetLinearDistance : constant close to surface then
+            //  fade like linearDistance
+
+            // Vary from surfaceCellSize (close to the surface) to
+            // distanceCellSize (further than 'distance')
+            linearDistanceCoeffs
+            {
+                surfaceCellSize     1e-5;    // absolute size
+                distanceCellSize    $defaultCellSize;
+                distance            1.0;
+            }
+        }
+
+        tunnel_APPROACH_INLET.obj
+        {
+            priority            1;
+            mode                bothSides;
+            cellSizeFunction    surfaceOffsetLinearDistance;
+            // Constant within a certain distance then linear fade away.
+            // Good for layers.
+            surfaceOffsetLinearDistanceCoeffs
+            {
+                surfaceCellSize     1e-5;
+                distanceCellSize    $defaultCellSize;
+                surfaceOffset       0.1;
+                totalDistance       1.0;
+            }
+        }
+    }
+
+    // Underrelaxation for point motion. Simulated annealing: starts off at 1
+    // and lowers to 0 (at simulation endTime) to converge points.
+    // adaptiveLinear is preferred choice.
+    // Points move by e.g. 10% of tet size.
+    relaxationModel     adaptiveLinear; //rampHoldFall
+
+    adaptiveLinearCoeffs
+    {
+        relaxationStart 1.0;
+        relaxationEnd   0.0;
+    }
+
+    // Output lots and lots of .obj files
+    objOutput                   no;
+
+    // Timing and memory usage.
+    timeChecks                  yes;
+
+    // Number of rays in plane parallel to nearest surface. Used to detect
+    // next closest surfaces. Used to work out alignment (three vectors)
+    // to surface.
+    // Note that only the initial points (from the seeding) calculate this
+    // information so if these are not fine enough the alignment will
+    // not be correct. (any points added during the running will lookup
+    // this information from the nearest initial point since it is
+    // expensive)
+    alignmentSearchSpokes       36;
+
+    // For each delaunay edge (between two vertices, becomes
+    // the Voronoi face normal) snap to the alignment direction if within
+    // alignmentAcceptanceAngle. Slightly > 45 is a good choice - prevents
+    // flipping.
+    alignmentAcceptanceAngle    48;
+
+    // How often to rebuild the alignment info (expensive)
+    sizeAndAlignmentRebuildFrequency 20;
+
+    // When to insert points. Not advisable change to
+    // these settings.
+    pointInsertionCriteria
+    {
+        // If edge larger than 1.75 target cell size
+        // (so tets too large/stretched) insert point
+        cellCentreDistCoeff     1.75;
+        // Do not insert point if voronoi face (on edge) very small.
+        faceAreaRatioCoeff      0.0025;
+        // Insert point only if edge closely aligned to local alignment
+        // direction.
+        acceptanceAngle         21.5;
+    }
+
+    // Opposite: remove point if mesh too compressed. Do not change these
+    // settings.
+    pointRemovalCriteria
+    {
+        cellCentreDistCoeff  0.65;
+    }
+
+    // How to determine the point motion. All edges got some direction.
+    // Sum all edge contributions to determine point motion. Weigh by
+    // face area so motion is preferentially determined by large faces
+    // (or more importantly ignore contribution from small faces).
+    // Do not change these settings.
+    faceAreaWeightModel piecewiseLinearRamp;
+
+    piecewiseLinearRampCoeffs
+    {
+        lowerAreaFraction       0.5;
+        upperAreaFraction       1.0;
+    }
+}
+
+
+// After simulation, when converting to polyMesh, filter out small faces/edges.
+// Do not change. See cvControls.H
+polyMeshFiltering
+{
+    // Write the underlying Delaunay tet mesh at output time
+    writeTetDualMesh                        false;  //true;
+
+    // Upper limit on the size of faces to be filtered.
+    // fraction of the local target cell size
+    filterSizeCoeff                         0.2;
+
+    // Upper limit on how close two dual vertices can be before
+    // being merged, fraction of the local target cell size
+    mergeClosenessCoeff                     1e-9;
+
+    // To not filter: set maxNonOrtho to 1 (so check fails) and then
+    // set continueFilteringOnBadInitialPolyMesh to false.
+    continueFilteringOnBadInitialPolyMesh   true;
+
+    // When a face is "bad", what fraction should the filterSizeCoeff be
+    // reduced by.  Recursive, so for a filterCount value of fC, the
+    // filterSizeCoeff is reduced by pow(filterErrorReductionCoeff, fC)
+    filterErrorReductionCoeff               0.5;
+
+    // Maximum number of filterCount applications before a face
+    // is not attempted to be filtered
+    filterCountSkipThreshold                4;
+
+    // Maximum number of permissible iterations of the face collapse
+    // algorithm.  The value to choose will be related the maximum number
+    // of points on a face that is to be collapsed and how many faces
+    // around it need to be collapsed.
+    maxCollapseIterations                   25;
+
+    // Maximum number of times an to allow an equal faceSet to be
+    // returned from the face quality assessment before stopping iterations
+    // to break an infinitie loop.
+    maxConsecutiveEqualFaceSets             5;
+    // Remove little steps (almost perp to surface) by collapsing face.
+    surfaceStepFaceAngle                    80;
+    // Do not collapse face to edge if should become edges
+    edgeCollapseGuardFraction               0.3;
+    // Only collapse face to point if high aspect ratio
+    maxCollapseFaceToPointSideLengthCoeff   0.35;
+}
+
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo. Same as in snappyHexMeshDict
+meshQualityControls
+{
+    //- Maximum non-orthogonality allowed. Set to 180 to disable.
+    maxNonOrtho 65;
+
+    //- Max skewness allowed. Set to <0 to disable.
+    maxBoundarySkewness 50;
+    maxInternalSkewness 10;
+
+    //- Max concaveness allowed. Is angle (in degrees) below which concavity
+    //  is allowed. 0 is straight face, <0 would be convex face.
+    //  Set to 180 to disable.
+    maxConcave 80;
+
+    //- Minimum pyramid volume. Is absolute volume of cell pyramid.
+    //  Set to a sensible fraction of the smallest cell volume expected.
+    //  Set to very negative number (e.g. -1E30) to disable.
+    minVol -1E30;
+
+    //- Minimum quality of the tet formed by the
+    //  variable base point minimum decomposition triangles and
+    //  the cell centre (so not face-centre decomposition).
+    //  This has to be a positive number for tracking
+    //  to work. Set to very negative number (e.g. -1E30) to
+    //  disable.
+    //     <0 = inside out tet,
+    //      0 = flat tet
+    //      1 = regular tet
+    minTetQuality 1e-30;
+
+    //- Minimum absolute face area. Set to <0 to disable.
+    minArea -1;
+
+    //- Minimum face twist. Set to <-1 to disable. dot product of face normal
+    //- and face centre triangles normal
+    minTwist 0.02;
+
+    //- minimum normalised cell determinant
+    //- 1 = hex, <= 0 = folded or flattened illegal cell
+    minDeterminant 0.001;
+
+    //- minFaceWeight (0 -> 0.5)
+    minFaceWeight 0.02;
+
+    //- minVolRatio (0 -> 1)
+    minVolRatio 0.01;
+
+    //must be >0 for Fluent compatibility
+    minTriangleTwist -1;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/manipulation/checkMesh/checkMesh.C

@@ -91,6 +91,26 @@ int main(int argc, char *argv[])
             // Clear mesh before checking
             mesh.clearOut();
 
+            pointIOField overrideCCs
+            (
+                IOobject
+                (
+                    "cellCentres",
+                    mesh.pointsInstance(),
+                    polyMesh::meshSubDir,
+                    runTime,
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::NO_WRITE
+                )
+            );
+
+            if (overrideCCs.headerOk())
+            {
+                Info<< "Read " << overrideCCs.size() << " cell centres" << endl;
+
+                mesh.overrideCellCentres(overrideCCs);
+            }
+
             // Reconstruct globalMeshData
             mesh.globalData();
 

applications/utilities/mesh/manipulation/topoSet/topoSet.C

@@ -208,7 +208,8 @@ int main(int argc, char *argv[])
                     IOobject::MUST_READ
                 );
                 Info<< "Read set " << setName << " with size "
-                    << currentSet().size() << endl;
+                    << returnReduce(currentSet().size(), sumOp<label>())
+                    << endl;
             }
 
 
@@ -292,7 +293,8 @@ int main(int argc, char *argv[])
             if (currentSet.valid())
             {
                 Info<< "    Set " << currentSet().name()
-                    << " now size " << currentSet().size()
+                    << " now size "
+                    << returnReduce(currentSet().size(), sumOp<label>())
                     << endl;
             }
         }

applications/utilities/postProcessing/dataConversion/foamToTetDualMesh/Make/files

@@ -0,0 +1,4 @@
+
+foamToTetDualMesh.C
+
+EXE = $(FOAM_USER_APPBIN)/foamToTetDualMesh

applications/utilities/postProcessing/dataConversion/foamToTetDualMesh/Make/options

@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume

applications/utilities/postProcessing/dataConversion/foamToTetDualMesh/foamToTetDualMesh.C

@@ -0,0 +1,313 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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
+    foamToTetDualMesh
+
+Description
+    Converts polyMesh results to tetDualMesh.
+
+\*---------------------------------------------------------------------------*/
+
+
+#include "argList.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "pointFields.H"
+#include "Time.H"
+#include "IOobjectList.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class ReadGeoField, class MappedGeoField>
+void ReadAndMapFields
+(
+    const fvMesh& mesh,
+    const IOobjectList& objects,
+    const fvMesh& tetDualMesh,
+    const labelList& map,
+    const typename MappedGeoField::value_type& nullValue,
+    PtrList<MappedGeoField>& tetFields
+)
+{
+    typedef typename MappedGeoField::value_type Type;
+
+    // Search list of objects for wanted type
+    IOobjectList fieldObjects(objects.lookupClass(ReadGeoField::typeName));
+
+    tetFields.setSize(fieldObjects.size());
+
+    label i = 0;
+    forAllConstIter(IOobjectList, fieldObjects, iter)
+    {
+        Info<< "Converting " << ReadGeoField::typeName << ' ' << iter.key()
+            << endl;
+
+        ReadGeoField readField(*iter(), mesh);
+
+        tetFields.set
+        (
+            i,
+            new MappedGeoField
+            (
+                IOobject
+                (
+                    readField.name(),
+                    readField.instance(),
+                    readField.local(),
+                    tetDualMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE,
+                    readField.registerObject()
+                ),
+                pointMesh::New(tetDualMesh),
+                dimensioned<Type>
+                (
+                    "zero",
+                    readField.dimensions(),
+                    pTraits<Type>::zero
+                )
+            )
+        );
+
+        Field<Type>& fld = tetFields[i].internalField();
+
+        // Map from read field. Set unmapped entries to nullValue.
+        fld.setSize(map.size(), nullValue);
+        forAll(map, pointI)
+        {
+            label index = map[pointI];
+
+            if (index > 0)
+            {
+                label cellI = index-1;
+                fld[pointI] = readField[cellI];
+            }
+            else if (index < 0)
+            {
+                label faceI = -index-1;
+                label bFaceI = faceI - mesh.nInternalFaces();
+                if (bFaceI >= 0)
+                {
+                    label patchI = mesh.boundaryMesh().patchID()[bFaceI];
+                    label localFaceI = mesh.boundaryMesh()[patchI].whichFace
+                    (
+                        faceI
+                    );
+                    fld[pointI] = readField.boundaryField()[patchI][localFaceI];
+                }
+                //else
+                //{
+                //    FatalErrorIn("ReadAndMapFields(..)")
+                //        << "Face " << faceI << " from index " << index
+                //        << " is not a boundary face." << abort(FatalError);
+                //}
+
+            }
+            //else
+            //{
+            //    WarningIn("ReadAndMapFields(..)")
+            //        << "Point " << pointI << " at "
+            //        << tetDualMesh.points()[pointI]
+            //        << " has no dual correspondence." << endl;
+            //}
+        }
+        tetFields[i].correctBoundaryConditions();
+
+        i++;
+    }
+}
+
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "addOverwriteOption.H"
+#   include "addTimeOptions.H"
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+    // Get times list
+    instantList Times = runTime.times();
+#   include "checkTimeOptions.H"
+    runTime.setTime(Times[startTime], startTime);
+
+
+    // Read the mesh
+#   include "createMesh.H"
+
+    // Read the tetDualMesh
+    Info<< "Create tetDualMesh for time = "
+        << runTime.timeName() << nl << endl;
+
+    fvMesh tetDualMesh
+    (
+        IOobject
+        (
+            "tetDualMesh",
+            runTime.timeName(),
+            runTime,
+            IOobject::MUST_READ
+        )
+    );
+    // From tet vertices to poly cells/faces
+    const labelIOList pointDualAddressing
+    (
+        IOobject
+        (
+            "pointDualAddressing",
+            tetDualMesh.facesInstance(),
+            tetDualMesh.meshSubDir,
+            tetDualMesh,
+            IOobject::MUST_READ
+        )
+    );
+
+    if (pointDualAddressing.size() != tetDualMesh.nPoints())
+    {
+            FatalErrorIn(args.executable())
+                << "Size " << pointDualAddressing.size()
+                << " of addressing map " << pointDualAddressing.objectPath()
+                << " differs from number of points in mesh "
+                << tetDualMesh.nPoints()
+                << exit(FatalError);
+    }
+
+
+    // Some stats on addressing
+    label nCells = 0;
+    label nPatchFaces = 0;
+    label nUnmapped = 0;
+    forAll(pointDualAddressing, pointI)
+    {
+        label index = pointDualAddressing[pointI];
+
+        if (index > 0)
+        {
+            nCells++;
+        }
+        else if (index == 0)
+        {
+            nUnmapped++;
+        }
+        else
+        {
+            label faceI = -index-1;
+            if (faceI < mesh.nInternalFaces())
+            {
+                FatalErrorIn(args.executable())
+                    << "Face " << faceI << " from index " << index
+                    << " is not a boundary face."
+                    << " nInternalFaces:" << mesh.nInternalFaces()
+                    << exit(FatalError);
+            }
+            else
+            {
+                nPatchFaces++;
+            }
+        }
+    }
+
+    reduce(nCells, sumOp<label>());
+    reduce(nPatchFaces, sumOp<label>());
+    reduce(nUnmapped, sumOp<label>());
+    Info<< "tetDualMesh points : " << tetDualMesh.nPoints()
+        << " of which mapped to" << nl
+        << "    cells       : " << nCells << nl
+        << "    patch faces : " << nPatchFaces << nl
+        << "    not mapped  : " << nUnmapped << nl
+        << endl;
+
+
+    // Read objects in time directory
+    IOobjectList objects(mesh, runTime.timeName());
+
+    // Read vol fields, interpolate onto tet points
+    PtrList<pointScalarField> psFlds;
+    ReadAndMapFields<volScalarField, pointScalarField>
+    (
+        mesh,
+        objects,
+        tetDualMesh,
+        pointDualAddressing,
+        pTraits<scalar>::zero,  // nullValue
+        psFlds
+    );
+
+    PtrList<pointVectorField> pvFlds;
+    ReadAndMapFields<volVectorField, pointVectorField>
+    (
+        mesh,
+        objects,
+        tetDualMesh,
+        pointDualAddressing,
+        pTraits<vector>::zero,  // nullValue
+        pvFlds
+    );
+
+    PtrList<pointSphericalTensorField> pstFlds;
+    ReadAndMapFields<volSphericalTensorField, pointSphericalTensorField>
+    (
+        mesh,
+        objects,
+        tetDualMesh,
+        pointDualAddressing,
+        pTraits<sphericalTensor>::zero,  // nullValue
+        pstFlds
+    );
+
+    PtrList<pointSymmTensorField> psymmtFlds;
+    ReadAndMapFields<volSymmTensorField, pointSymmTensorField>
+    (
+        mesh,
+        objects,
+        tetDualMesh,
+        pointDualAddressing,
+        pTraits<symmTensor>::zero,  // nullValue
+        psymmtFlds
+    );
+
+    PtrList<pointTensorField> ptFlds;
+    ReadAndMapFields<volTensorField, pointTensorField>
+    (
+        mesh,
+        objects,
+        tetDualMesh,
+        pointDualAddressing,
+        pTraits<tensor>::zero,  // nullValue
+        ptFlds
+    );
+
+    tetDualMesh.objectRegistry::write();
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/mdInitialise/mdInitialise.C

@@ -2,7 +2,7 @@
  =========                   |
  \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
   \\    /    O peration      |
-   \\  /     A nd            | Copyright (C) 2008-2010 OpenCFD Ltd.
+   \\  /     A nd            | Copyright (C) 2008-2011 OpenCFD Ltd.
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License

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,462 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2011 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 extendedFeatureEdgeMesh 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 "extendedFeatureEdgeMesh.H"
+#include "triSurfaceSearch.H"
+#include "OFstream.H"
+#include "booleanSurface.H"
+#include "edgeIntersections.H"
+#include "meshTools.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::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);
+    }
+
+    fileName surf1Name(args[2]);
+    Info<< "Reading surface " << surf1Name << endl;
+    triSurface surf1(surf1Name);
+
+    Info<< surf1Name << " statistics:" << endl;
+    surf1.writeStats(Info);
+    Info<< endl;
+
+    fileName surf2Name(args[3]);
+    Info<< "Reading surface " << surf2Name << endl;
+    triSurface surf2(surf2Name);
+
+    Info<< surf2Name << " statistics:" << endl;
+    surf2.writeStats(Info);
+    Info<< 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);
+
+    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");
+
+    scalar scale = 0.05*min
+    (
+        querySurf1.tree().bb().mag(),
+        querySurf2.tree().bb().mag()
+    );
+
+    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());
+
+        {
+            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
+
+    extendedFeatureEdgeMesh feMesh
+    (
+        IOobject
+        (
+            sFeatFileName + ".extendedFeatureEdgeMesh",
+            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;
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceCheck/surfaceCheck.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -315,32 +315,12 @@ int main(int argc, char *argv[])
             }
             else
             {
-                triPointRef tri
+                triQ[faceI] = triPointRef
                 (
                     surf.points()[f[0]],
                     surf.points()[f[1]],
                     surf.points()[f[2]]
-                );
-
-                vector ba(tri.b() - tri.a());
-                ba /= mag(ba) + VSMALL;
-
-                vector ca(tri.c() - tri.a());
-                ca /= mag(ca) + VSMALL;
-
-                if (mag(ba&ca) > 1-1E-3)
-                {
-                    triQ[faceI] = SMALL;
-                }
-                else
-                {
-                    triQ[faceI] = triPointRef
-                    (
-                        surf.points()[f[0]],
-                        surf.points()[f[1]],
-                        surf.points()[f[2]]
-                    ).quality();
-                }
+                ).quality();
             }
         }
 

applications/utilities/surface/surfaceClean/collapseBase.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,8 +66,8 @@ static void writeRegionOBJ
 
     triSurface regionSurf(surf.subsetMesh(include, pointMap, faceMap));
 
-    //Pout<< "Region " << regionI << " surface:" << nl;
-    //regionSurf.writeStats(Pout);
+    Pout<< "Region " << regionI << " surface:" << nl;
+    regionSurf.writeStats(Pout);
 
     regionSurf.write(regionName);
 
@@ -97,7 +97,7 @@ static void splitTri
     DynamicList<labelledTri>& tris
 )
 {
-    label oldNTris = tris.size();
+    //label oldNTris = tris.size();
 
     label fp = findIndex(f, e[0]);
     label fp1 = f.fcIndex(fp);
@@ -175,7 +175,7 @@ static void splitTri
     }
     else
     {
-        FatalErrorIn("splitTri")
+        FatalErrorIn("splitTri(..)")
             << "Edge " << e << " not part of triangle " << f
             << " fp:" << fp
             << " fp1:" << fp1
@@ -183,13 +183,13 @@ static void splitTri
             << abort(FatalError);
     }
 
-    Pout<< "Split face " << f << " along edge " << e
-        << " into triangles:" << endl;
-
-    for (label i = oldNTris; i < tris.size(); i++)
-    {
-        Pout<< "   " << tris[i] << nl;
-    }
+    //Pout<< "Split face " << f << " along edge " << e
+    //    << " into triangles:" << endl;
+    //
+    //for (label i = oldNTris; i < tris.size(); i++)
+    //{
+    //    Pout<< "   " << tris[i] << nl;
+    //}
 }
 
 
@@ -206,14 +206,14 @@ static bool insertSorted
 {
     if (findIndex(sortedVerts, vertI) != -1)
     {
-        FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
+        FatalErrorIn("insertSorted(..)") << "Trying to insert vertex " << vertI
             << " which is already in list of sorted vertices "
             << sortedVerts << abort(FatalError);
     }
 
     if (weight <= 0 || weight >= 1)
     {
-        FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
+        FatalErrorIn("insertSorted(..)") << "Trying to insert vertex " << vertI
             << " with illegal weight " << weight
             << " into list of sorted vertices "
             << sortedVerts << abort(FatalError);
@@ -228,7 +228,7 @@ static bool insertSorted
 
         if (mag(w - weight) < SMALL)
         {
-            WarningIn("insertSorted")
+            WarningIn("insertSorted(..)")
                 << "Trying to insert weight " << weight << " which is close to"
                 << " existing weight " << w << " in " << sortedWeights
                 << endl;
@@ -263,64 +263,103 @@ static bool insertSorted
 }
 
 
+// Is triangle candidate for collapse? Small height or small quality
+bool isSliver
+(
+    const triSurface& surf,
+    const scalar minLen,
+    const scalar minQuality,
+    const label faceI,
+    const label edgeI
+)
+{
+    const pointField& localPoints = surf.localPoints();
+
+    // Check
+    // - opposite vertex projects onto base edge
+    // - normal distance is small
+    // - or triangle quality is small
+
+    label opposite0 =
+        triSurfaceTools::oppositeVertex
+        (
+            surf,
+            faceI,
+            edgeI
+        );
+
+    const edge& e = surf.edges()[edgeI];
+    const labelledTri& f = surf[faceI];
+
+    pointHit pHit =
+        e.line(localPoints).nearestDist
+        (
+            localPoints[opposite0]
+        );
+
+    if
+    (
+        pHit.hit()
+     && (
+            pHit.distance() < minLen
+         || f.tri(surf.points()).quality() < minQuality
+        )
+    )
+    {
+        // Remove faceI and split all other faces using this
+        // edge. This is done by 'replacing' the edgeI with the
+        // opposite0 vertex
+        //Pout<< "Splitting face " << faceI << " since distance "
+        //    << pHit.distance()
+        //    << " from vertex " << opposite0
+        //    << " to edge " << edgeI
+        //    << "  points "
+        //    << localPoints[e[0]]
+        //    << localPoints[e[1]]
+        //    << " is too small or triangle quality "
+        //    << f.tri(surf.points()).quality()
+        //    << " too small." << endl;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 // Mark all faces that are going to be collapsed.
 // faceToEdge: per face -1 or the base edge of the face.
 static void markCollapsedFaces
 (
     const triSurface& surf,
     const scalar minLen,
+    const scalar minQuality,
     labelList& faceToEdge
 )
 {
     faceToEdge.setSize(surf.size());
     faceToEdge = -1;
 
-    const pointField& localPoints = surf.localPoints();
     const labelListList& edgeFaces = surf.edgeFaces();
 
     forAll(edgeFaces, edgeI)
     {
-        const edge& e = surf.edges()[edgeI];
-
         const labelList& eFaces = surf.edgeFaces()[edgeI];
 
         forAll(eFaces, i)
         {
             label faceI = eFaces[i];
 
-            // Check distance of vertex to edge.
-            label opposite0 =
-                triSurfaceTools::oppositeVertex
-                (
-                    surf,
-                    faceI,
-                    edgeI
-                );
+            bool isCandidate = isSliver(surf, minLen, minQuality, faceI, edgeI);
 
-            pointHit pHit =
-                e.line(localPoints).nearestDist
-                (
-                    localPoints[opposite0]
-                );
-
-            if (pHit.hit() && pHit.distance() < minLen)
+            if (isCandidate)
             {
-                // Remove faceI and split all other faces using this
-                // edge. This is done by 'replacing' the edgeI with the
-                // opposite0 vertex
-                Pout<< "Splitting face " << faceI << " since distance "
-                    << pHit.distance()
-                    << " from vertex " << opposite0
-                    << " to edge " << edgeI
-                    << "  points "
-                    << localPoints[e[0]]
-                    << localPoints[e[1]]
-                    << " is too small" << endl;
-
                 // Mark face as being collapsed
                 if (faceToEdge[faceI] != -1)
                 {
-                    FatalErrorIn("markCollapsedFaces")
+                    FatalErrorIn("markCollapsedFaces(..)")
                         << "Cannot collapse face " << faceI << " since "
                         << " is marked to be collapsed both to edge "
                         << faceToEdge[faceI] << " and " << edgeI
@@ -347,7 +386,7 @@ static void markRegion
 {
     if (faceToEdge[faceI] == -1 || collapseRegion[faceI] != -1)
     {
-        FatalErrorIn("markRegion")
+        FatalErrorIn("markRegion(..)")
             << "Problem : crossed into uncollapsed/regionized face"
             << abort(FatalError);
     }
@@ -383,7 +422,7 @@ static void markRegion
                 }
                 else if (collapseRegion[nbrFaceI] != regionI)
                 {
-                    FatalErrorIn("markRegion")
+                    FatalErrorIn("markRegion(..)")
                         << "Edge:" << edgeI << " between face " << faceI
                         << " with region " << regionI
                         << " and face " << nbrFaceI
@@ -411,8 +450,8 @@ static label markRegions
     {
         if (collapseRegion[faceI] == -1 && faceToEdge[faceI] != -1)
         {
-            Pout<< "markRegions : Marking region:" << regionI
-                << " starting from face " << faceI << endl;
+            //Pout<< "markRegions : Marking region:" << regionI
+            //    << " starting from face " << faceI << endl;
 
             // Collapsed face. Mark connected region with current region number
             markRegion(surf, faceToEdge, regionI++, faceI, collapseRegion);
@@ -728,7 +767,12 @@ static void getSplitVerts
 }
 
 
-label collapseBase(triSurface& surf, const scalar minLen)
+label collapseBase
+(
+    triSurface& surf,
+    const scalar minLen,
+    const scalar minQuality
+)
 {
     label nTotalSplit = 0;
 
@@ -743,7 +787,7 @@ label collapseBase(triSurface& surf, const scalar minLen)
         labelList faceToEdge(surf.size(), -1);
 
         // Calculate faceToEdge (face collapses)
-        markCollapsedFaces(surf, minLen, faceToEdge);
+        markCollapsedFaces(surf, minLen, minQuality, faceToEdge);
 
 
         // Find regions of connected collapsed faces
@@ -754,8 +798,8 @@ label collapseBase(triSurface& surf, const scalar minLen)
 
         label nRegions = markRegions(surf, faceToEdge, collapseRegion);
 
-        Pout<< "Detected " << nRegions << " regions of faces to be collapsed"
-            << nl << endl;
+        //Pout<< "Detected " << nRegions << " regions of faces to be collapsed"
+        //    << nl << endl;
 
         // Pick up all vertices on outside of region
         labelListList outsideVerts
@@ -772,10 +816,10 @@ label collapseBase(triSurface& surf, const scalar minLen)
         {
             spanPoints[regionI] = getSpanPoints(surf, outsideVerts[regionI]);
 
-            Pout<< "For region " << regionI << " found extrema at points "
-                << surf.localPoints()[spanPoints[regionI][0]]
-                << surf.localPoints()[spanPoints[regionI][1]]
-                << endl;
+            //Pout<< "For region " << regionI << " found extrema at points "
+            //    << surf.localPoints()[spanPoints[regionI][0]]
+            //    << surf.localPoints()[spanPoints[regionI][1]]
+            //    << endl;
 
             // Project all non-span points onto the span edge.
             projectNonSpanPoints
@@ -787,21 +831,21 @@ label collapseBase(triSurface& surf, const scalar minLen)
                 orderedWeights[regionI]
             );
 
-            Pout<< "For region:" << regionI
-                << " span:" << spanPoints[regionI]
-                << " orderedVerts:" << orderedVertices[regionI]
-                << " orderedWeights:" << orderedWeights[regionI]
-                << endl;
+            //Pout<< "For region:" << regionI
+            //    << " span:" << spanPoints[regionI]
+            //    << " orderedVerts:" << orderedVertices[regionI]
+            //    << " orderedWeights:" << orderedWeights[regionI]
+            //    << endl;
 
-            writeRegionOBJ
-            (
-                surf,
-                regionI,
-                collapseRegion,
-                outsideVerts[regionI]
-            );
+            //writeRegionOBJ
+            //(
+            //    surf,
+            //    regionI,
+            //    collapseRegion,
+            //    outsideVerts[regionI]
+            //);
 
-            Pout<< endl;
+            //Pout<< endl;
         }
 
 
@@ -864,20 +908,19 @@ label collapseBase(triSurface& surf, const scalar minLen)
                     // Split edge using splitVerts. All non-collapsed triangles
                     // using edge will get split.
 
-
-                    {
-                        const pointField& localPoints = surf.localPoints();
-                        Pout<< "edge " << edgeI << ' ' << e
-                            << "  points "
-                            << localPoints[e[0]] << ' ' << localPoints[e[1]]
-                            << " split into edges with extra points:"
-                            << endl;
-                        forAll(splitVerts, i)
-                        {
-                            Pout<< "    " << splitVerts[i] << " weight "
-                                << splitWeights[i] << nl;
-                        }
-                    }
+                    //{
+                    //    const pointField& localPoints = surf.localPoints();
+                    //    Pout<< "edge " << edgeI << ' ' << e
+                    //        << "  points "
+                    //        << localPoints[e[0]] << ' ' << localPoints[e[1]]
+                    //        << " split into edges with extra points:"
+                    //        << endl;
+                    //    forAll(splitVerts, i)
+                    //    {
+                    //        Pout<< "    " << splitVerts[i] << " weight "
+                    //            << splitWeights[i] << nl;
+                    //    }
+                    //}
 
                     const labelList& eFaces = surf.edgeFaces()[edgeI];
 
@@ -914,7 +957,8 @@ label collapseBase(triSurface& surf, const scalar minLen)
             }
         }
 
-        Pout<< "collapseBase : splitting " << nSplit << " triangles"
+        Info<< "collapseBase : collapsing " << nSplit
+            << " triangles by splitting their base edge."
             << endl;
 
         nTotalSplit += nSplit;
@@ -927,15 +971,15 @@ label collapseBase(triSurface& surf, const scalar minLen)
         // Pack the triangles
         newTris.shrink();
 
-        Pout<< "Resetting surface from " << surf.size() << " to "
-            << newTris.size() << " triangles" << endl;
+        //Pout<< "Resetting surface from " << surf.size() << " to "
+        //    << newTris.size() << " triangles" << endl;
         surf = triSurface(newTris, surf.patches(), surf.localPoints());
 
-        {
-            fileName fName("bla" + name(iter) + ".obj");
-            Pout<< "Writing surf to " << fName << endl;
-            surf.write(fName);
-        }
+        //{
+        //    fileName fName("bla" + name(iter) + ".obj");
+        //    Pout<< "Writing surf to " << fName << endl;
+        //    surf.write(fName);
+        //}
 
         iter++;
     }

applications/utilities/surface/surfaceClean/collapseBase.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,9 +41,14 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//- Keep collapsing all triangles whose height is < minLen.
+//- Keep collapsing all triangles whose height is < minLen or quality < minQ.
 //  Returns number of triangles collapsed.
-label collapseBase(triSurface& surf, const scalar minLen);
+label collapseBase
+(
+    triSurface& surf,
+    const scalar minLen,
+    const scalar minQuality
+);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/surface/surfaceClean/collapseEdge.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -121,8 +121,8 @@ label collapseEdge(triSurface& surf, const scalar minLen)
                         pointMap[v1] = v;
                         newPoints[v] = 0.5*(localPoints[v1] + localPoints[v]);
 
-                        Pout<< "Collapsing triange " << faceI << " to edge mid "
-                            << newPoints[v] << endl;
+                        //Pout<< "Collapsing triange " << faceI
+                        //    << " to edge mid " << newPoints[v] << endl;
 
                         nCollapsed++;
                         okToCollapse[faceI] = false;
@@ -136,7 +136,8 @@ label collapseEdge(triSurface& surf, const scalar minLen)
             }
         }
 
-        Pout<< "collapseEdge : collapsing " << nCollapsed << " triangles"
+        Info<< "collapseEdge : collapsing " << nCollapsed
+            << " triangles to a single edge."
             << endl;
 
         nTotalCollapsed += nCollapsed;

applications/utilities/surface/surfaceClean/surfaceClean.C

@@ -50,6 +50,7 @@ int main(int argc, char *argv[])
     argList::noParallel();
     argList::validArgs.append("surfaceFile");
     argList::validArgs.append("min length");
+    argList::validArgs.append("min quality");
     argList::validArgs.append("output surfaceFile");
     argList::addBoolOption
     (
@@ -60,10 +61,13 @@ int main(int argc, char *argv[])
 
     const fileName inFileName = args[1];
     const scalar minLen = args.argRead<scalar>(2);
-    const fileName outFileName = args[3];
+    const scalar minQuality = args.argRead<scalar>(3);
+    const fileName outFileName = args[4];
 
     Info<< "Reading surface " << inFileName << nl
-        << "Collapsing all triangles with edges or heights < " << minLen << nl
+        << "Collapsing all triangles with" << nl
+        << "    edges or heights < " << minLen << nl
+        << "    quality          < " << minQuality << nl
         << "Writing result to " << outFileName << nl << endl;
 
 
@@ -90,7 +94,7 @@ int main(int argc, char *argv[])
     }
     while (true)
     {
-        label nSplitEdge = collapseBase(surf, minLen);
+        label nSplitEdge = collapseBase(surf, minLen, minQuality);
 
         if (nSplitEdge == 0)
         {

applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H

@@ -0,0 +1,203 @@
+#ifndef CGAL_PSURF_RINGS_H_
+#define CGAL_PSURF_RINGS_H_
+
+// This file adapted from
+//     CGAL-3.7/examples/Jet_fitting_3//PolyhedralSurf_rings.h
+// Licensed under CGAL-3.7/LICENSE.FREE_USE
+
+// Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
+// Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
+// Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
+// Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
+// Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
+// (Israel).  All rights reserved.
+
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+
+#include <cassert>
+
+using namespace std;
+
+template
+<
+    class TPoly,
+    class VertexPropertyMap
+>
+
+class T_PolyhedralSurf_rings
+{
+
+protected:
+
+  //Polyhedron
+  typedef typename TPoly::Vertex Vertex;
+  typedef typename TPoly::Halfedge Halfedge;
+  typedef typename TPoly::Facet Facet;
+  typedef typename TPoly::Halfedge_around_vertex_circulator
+      Halfedge_around_vertex_circulator;
+  typedef typename TPoly::Vertex_iterator Vertex_iterator;
+
+  //vertex indices are initialised to -1
+  static void reset_ring_indices(std::vector < Vertex * >&vces,
+                                 VertexPropertyMap& vpm);
+
+  // i >= 1; from a start vertex on the current i-1 ring, push non-visited
+  // neighbors of start in the nextRing and set indices to i. Also add these
+  // vertices in all.
+  static void push_neighbours_of(Vertex * start, int ith,
+                          std::vector < Vertex * >&nextRing,
+                          std::vector < Vertex * >&all,
+                          VertexPropertyMap& vpm);
+
+  // i >= 1, from a currentRing i-1, collect all neighbors, set indices
+  // to i and store them in nextRing and all.
+  static void collect_ith_ring(int ith,
+                        std::vector < Vertex * >&currentRing,
+                        std::vector < Vertex * >&nextRing,
+                        std::vector < Vertex * >&all,
+                        VertexPropertyMap& vpm);
+
+ public:
+
+  // collect i>=1 rings : all neighbours up to the ith ring,
+  static void
+    collect_i_rings(Vertex* v,
+                    int ring_i,
+                    std::vector < Vertex * >& all,
+                    VertexPropertyMap& vpm);
+
+  //collect enough rings (at least 1), to get at least min_nb of neighbors
+  static void
+    collect_enough_rings(Vertex* v,
+                         unsigned int min_nb,
+                         std::vector < Vertex * >& all,
+                         VertexPropertyMap& vpm);
+};
+
+
+////IMPLEMENTATION////////////////////////////////////////////////////
+
+template < class TPoly , class VertexPropertyMap>
+void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+push_neighbours_of(Vertex * start, int ith,
+                   std::vector < Vertex * >&nextRing,
+                   std::vector < Vertex * >&all,
+                   VertexPropertyMap& vpm)
+{
+  Vertex *v;
+  Halfedge_around_vertex_circulator
+    hedgeb = start->vertex_begin(), hedgee = hedgeb;
+
+ CGAL_For_all(hedgeb, hedgee)
+  {
+    v = &*(hedgeb->opposite()->vertex());
+    if(get(vpm, v) != -1) continue;//if visited: next
+
+    put(vpm, v, ith);
+    nextRing.push_back(v);
+    all.push_back(v);
+  }
+}
+
+
+template <class TPoly, class VertexPropertyMap>
+void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+collect_ith_ring(int ith, std::vector < Vertex * >&currentRing,
+                 std::vector < Vertex * >&nextRing,
+                 std::vector < Vertex * >&all,
+                 VertexPropertyMap& vpm)
+{
+  typename std::vector < Vertex * >::iterator
+    itb =    currentRing.begin(), ite = currentRing.end();
+
+  CGAL_For_all(itb, ite) push_neighbours_of(*itb, ith, nextRing, all, vpm);
+}
+
+
+template <class TPoly, class VertexPropertyMap>
+  void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+reset_ring_indices(std::vector < Vertex * >&vces,
+                   VertexPropertyMap& vpm)
+{
+  typename std::vector < Vertex * >::iterator
+    itb = vces.begin(), ite = vces.end();
+  CGAL_For_all(itb, ite)  put(vpm, *itb, -1);
+}
+
+
+template <class TPoly, class VertexPropertyMap>
+  void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+collect_i_rings(Vertex* v,
+                int ring_i,
+                std::vector < Vertex * >& all,
+                VertexPropertyMap& vpm)
+{
+  std::vector<Vertex*> current_ring, next_ring;
+  std::vector<Vertex*> *p_current_ring, *p_next_ring;
+  assert(ring_i >= 1);
+  //initialize
+  p_current_ring = &current_ring;
+  p_next_ring = &next_ring;
+  put(vpm, v, 0);
+  current_ring.push_back(v);
+  all.push_back(v);
+
+  for (int i=1; i<=ring_i; i++)
+    {
+      collect_ith_ring(i, *p_current_ring, *p_next_ring, all, vpm);
+      //next round must be launched from p_nextRing...
+      p_current_ring->clear();
+      std::swap(p_current_ring, p_next_ring);
+    }
+  //clean up
+  reset_ring_indices(all, vpm);
+}
+
+
+template <class TPoly, class VertexPropertyMap>
+  void T_PolyhedralSurf_rings <TPoly, VertexPropertyMap>::
+collect_enough_rings(Vertex* v,
+                     unsigned int min_nb,
+                     std::vector < Vertex * >& all,
+                     VertexPropertyMap& vpm)
+{
+  std::vector<Vertex*> current_ring, next_ring;
+  std::vector<Vertex*> *p_current_ring, *p_next_ring;
+
+  //initialize
+  p_current_ring = &current_ring;
+  p_next_ring = &next_ring;
+  put(vpm, v, 0);
+  current_ring.push_back(v);
+  all.push_back(v);
+
+  int i = 1;
+
+  while ( (all.size() < min_nb) &&  (p_current_ring->size() != 0) )
+    {
+      collect_ith_ring(i, *p_current_ring, *p_next_ring, all, vpm);
+      //next round must be launched from p_nextRing...
+      p_current_ring->clear();
+      std::swap(p_current_ring, p_next_ring);
+      i++;
+    }
+  //clean up
+  reset_ring_indices(all, vpm);
+}
+
+
+#endif

applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "buildCGALPolyhedron.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::buildCGALPolyhedron::buildCGALPolyhedron
+(
+    const Foam::triSurface& surf
+)
+:
+    CGAL::Modifier_base<HalfedgeDS>(),
+    surf_(surf)
+{}
+
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::buildCGALPolyhedron::~buildCGALPolyhedron()
+{}
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+void Foam::buildCGALPolyhedron::operator()
+(
+    HalfedgeDS& hds
+)
+{
+    typedef HalfedgeDS::Traits     Traits;
+    typedef Traits::Point_3        Point;
+
+    // Postcondition: `hds' is a valid polyhedral surface.
+    CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> B(hds, false);
+
+    B.begin_surface
+    (
+        surf_.points().size(),      // n points
+        surf_.size(),               // n facets
+        2*surf_.edges().size()      // n halfedges
+    );
+
+    forAll(surf_.points(), pI)
+    {
+        const Foam::point& p = surf_.points()[pI];
+
+        B.add_vertex(Point(p.x(), p.y(), p.z()));
+    }
+
+    forAll(surf_, fI)
+    {
+        B.begin_facet();
+
+        B.add_vertex_to_facet(surf_[fI][0]);
+        B.add_vertex_to_facet(surf_[fI][1]);
+        B.add_vertex_to_facet(surf_[fI][2]);
+
+        B.end_facet();
+    }
+
+    B.end_surface();
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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/>.
+
+Class
+    Foam::buildCGALPolyhedron
+
+Description
+    Convert a triSurface into a CGAL Polyhedron
+
+SourceFiles
+    buildCGALPolyhedron.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef buildCGALPolyhedron_H
+#define buildCGALPolyhedron_H
+
+#include "triSurface.H"
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Polyhedron_incremental_builder_3.h>
+#include <CGAL/Polyhedron_3.h>
+
+typedef CGAL::Simple_cartesian<double>               Kernel;
+typedef CGAL::Polyhedron_3<Kernel>                   Polyhedron;
+typedef Polyhedron::HalfedgeDS                       HalfedgeDS;
+typedef Polyhedron::Vertex                           Vertex;
+typedef Polyhedron::Vertex_iterator                  Vertex_iterator;
+typedef Kernel::Point_3                              Point_3;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class buildCGALPolyhedron Declaration
+\*---------------------------------------------------------------------------*/
+
+class buildCGALPolyhedron
+:
+    public CGAL::Modifier_base<HalfedgeDS>
+{
+    // Private data
+
+        //- Reference to triSurface to convert
+        const Foam::triSurface& surf_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        buildCGALPolyhedron(const buildCGALPolyhedron&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const buildCGALPolyhedron&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct with reference to triSurface
+        buildCGALPolyhedron(const triSurface& surf);
+
+
+    //- Destructor
+    ~buildCGALPolyhedron();
+
+
+    // Member Operators
+
+        //- operator() of this `modifier' called by delegate function of
+        //  Polyhedron
+        void operator()(HalfedgeDS& hds);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/Make/files

@@ -1,3 +1,4 @@
 surfaceFeatureExtract.C
+CGALPolyhedron/buildCGALPolyhedron.C
 
 EXE = $(FOAM_APPBIN)/surfaceFeatureExtract

applications/utilities/surface/surfaceFeatureExtract/Make/options

@@ -1,9 +1,29 @@
+EXE_FROUNDING_MATH = -frounding-math
+EXE_NDEBUG = -DNDEBUG
+USE_F2C = -DCGAL_USE_F2C
+
+include $(GENERAL_RULES)/CGAL
+
 EXE_INC = \
+    ${EXE_FROUNDING_MATH} \
+    ${EXE_NDEBUG} \
+    ${USE_F2C} \
+    ${CGAL_INC} \
+    -ICGALPolyhedron \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
+    $(CGAL_LIBS) \
+    -L$(CGAL_ARCH_PATH)/lib \
+    -llapack \
+    -lblas \
+    -lCGAL \
     -lmeshTools \
     -ledgeMesh \
-    -ltriSurface
+    -ltriSurface \
+    -lsampling

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -40,7 +40,19 @@ Description
 #include "treeBoundBox.H"
 #include "meshTools.H"
 #include "OFstream.H"
+#include "triSurfaceMesh.H"
+#include "vtkSurfaceWriter.H"
+#include "triSurfaceFields.H"
 #include "unitConversion.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "unitConversion.H"
+
+#include "buildCGALPolyhedron.H"
+#include "CGALPolyhedronRings.H"
+#include <CGAL/Monge_via_jet_fitting.h>
+#include <CGAL/Lapack/Linear_algebra_lapack.h>
+#include <CGAL/property_map.h>
 
 using namespace Foam;
 
@@ -91,6 +103,201 @@ void deleteBox
 }
 
 
+void drawHitProblem
+(
+    label fI,
+    const triSurface& surf,
+    const pointField& start,
+    const pointField& faceCentres,
+    const pointField& end,
+    const List<pointIndexHit>& hitInfo
+)
+{
+    Info<< nl << "# findLineAll did not hit its own face."
+        << nl << "# fI " << fI
+        << nl << "# start " << start[fI]
+        << nl << "# f centre " << faceCentres[fI]
+        << nl << "# end " << end[fI]
+        << nl << "# hitInfo " << hitInfo
+        << endl;
+
+    meshTools::writeOBJ(Info, start[fI]);
+    meshTools::writeOBJ(Info, faceCentres[fI]);
+    meshTools::writeOBJ(Info, end[fI]);
+
+    Info<< "l 1 2 3" << endl;
+
+    meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
+    meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
+    meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
+
+    Info<< "f 4 5 6" << endl;
+
+    forAll(hitInfo, hI)
+    {
+        label hFI = hitInfo[hI].index();
+
+        meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
+        meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
+        meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
+
+        Info<< "f "
+            << 3*hI + 7 << " "
+            << 3*hI + 8 << " "
+            << 3*hI + 9
+            << endl;
+    }
+}
+
+
+scalarField calcCurvature(const triSurface& surf)
+{
+    scalarField k(surf.points().size(), 0);
+
+    Polyhedron P;
+
+    buildCGALPolyhedron convert(surf);
+    P.delegate(convert);
+
+    // Info<< "Created CGAL Polyhedron with " << label(P.size_of_vertices())
+    //     << " vertices and " << label(P.size_of_facets())
+    //     << " facets. " << endl;
+
+    // The rest of this function adapted from
+    //     CGAL-3.7/examples/Jet_fitting_3/Mesh_estimation.cpp
+    // Licensed under CGAL-3.7/LICENSE.FREE_USE
+
+    // Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
+    // Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
+    // Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
+    // Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
+    // Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
+    // (Israel).  All rights reserved.
+
+    // Permission is hereby granted, free of charge, to any person obtaining a
+    // copy of this software and associated documentation files (the
+    // "Software"), to deal in the Software without restriction, including
+    // without limitation the rights to use, copy, modify, merge, publish,
+    // distribute, sublicense, and/or sell copies of the Software, and to permit
+    // persons to whom the Software is furnished to do so, subject to the
+    // following conditions:
+
+    // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+    // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+    // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+    // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT
+    // OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
+    // THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+     //Vertex property map, with std::map
+    typedef std::map<Vertex*, int> Vertex2int_map_type;
+    typedef boost::associative_property_map< Vertex2int_map_type >
+        Vertex_PM_type;
+    typedef T_PolyhedralSurf_rings<Polyhedron, Vertex_PM_type > Poly_rings;
+
+    typedef CGAL::Monge_via_jet_fitting<Kernel>         Monge_via_jet_fitting;
+    typedef Monge_via_jet_fitting::Monge_form           Monge_form;
+
+    std::vector<Point_3> in_points;  //container for data points
+
+    // default parameter values and global variables
+    unsigned int d_fitting = 2;
+    unsigned int d_monge = 2;
+    unsigned int min_nb_points = (d_fitting + 1)*(d_fitting + 2)/2;
+
+    //initialize the tag of all vertices to -1
+    Vertex_iterator vitb = P.vertices_begin();
+    Vertex_iterator vite = P.vertices_end();
+
+    Vertex2int_map_type vertex2props;
+    Vertex_PM_type vpm(vertex2props);
+
+    CGAL_For_all(vitb, vite)
+    {
+        put(vpm, &(*vitb), -1);
+    }
+
+    vite = P.vertices_end();
+
+    label vertI = 0;
+
+    for (vitb = P.vertices_begin(); vitb != vite; vitb++)
+    {
+        //initialize
+        Vertex* v = &(*vitb);
+
+        //gather points around the vertex using rings
+        // From: gather_fitting_points(v, in_points, vpm);
+        {
+            std::vector<Vertex*> gathered;
+            in_points.clear();
+
+            Poly_rings::collect_enough_rings(v, min_nb_points, gathered, vpm);
+
+            //store the gathered points
+            std::vector<Vertex*>::iterator itb = gathered.begin();
+            std::vector<Vertex*>::iterator ite = gathered.end();
+
+            CGAL_For_all(itb, ite)
+            {
+                in_points.push_back((*itb)->point());
+            }
+        }
+
+        //skip if the nb of points is to small
+        if ( in_points.size() < min_nb_points )
+        {
+            std::cerr
+                << "not enough pts for fitting this vertex"
+                << in_points.size()
+                << std::endl;
+
+            continue;
+        }
+
+        // perform the fitting
+        Monge_via_jet_fitting monge_fit;
+
+        Monge_form monge_form = monge_fit
+        (
+            in_points.begin(),
+            in_points.end(),
+            d_fitting,
+            d_monge
+        );
+
+        // std::cout<< monge_form << std::endl;
+
+        // std::cout<< "k1 " << monge_form.principal_curvatures(0) << std::endl;
+        // std::cout<< "k2 " << monge_form.principal_curvatures(1) << std::endl;
+
+        // std::vector<Point_3>::iterator itbp = in_points.begin();
+        // std::vector<Point_3>::iterator itep = in_points.end();
+
+        // std::cout << "in_points list : " << std::endl;
+
+        // for (; itbp != itep; itbp++)
+        // {
+        //     std::cout << *itbp << std::endl;
+        // }
+
+        // std::cout << "--- vertex " <<  vertI
+        //     << " : " << v->point() << std::endl
+        //     << "number of points used : " << in_points.size()
+        //     << std::endl;
+
+        k[vertI++] = Foam::sqrt
+        (
+            sqr(monge_form.principal_curvatures(0))
+          + sqr(monge_form.principal_curvatures(1))
+        );
+    }
+
+    return k;
+}
+
+
 // Unmark non-manifold edges if individual triangles are not features
 void unmarkBaffles
 (
@@ -111,7 +318,7 @@ void unmarkBaffles
         {
             label i0 = eFaces[0];
             //const labelledTri& f0 = surf[i0];
-            const vector& n0 = surf.faceNormals()[i0];
+            const Foam::vector& n0 = surf.faceNormals()[i0];
 
             //Pout<< "edge:" << edgeI << " n0:" << n0 << endl;
 
@@ -120,7 +327,7 @@ void unmarkBaffles
             for (label i = 1; i < eFaces.size(); i++)
             {
                 //const labelledTri& f = surf[i];
-                const vector& n = surf.faceNormals()[eFaces[i]];
+                const Foam::vector& n = surf.faceNormals()[eFaces[i]];
 
                 //Pout<< "    mag(n&n0): " << mag(n&n0) << endl;
 
@@ -193,15 +400,63 @@ int main(int argc, char *argv[])
         "writeObj",
         "write extendedFeatureEdgeMesh obj files"
     );
+    argList::addBoolOption
+    (
+        "writeVTK",
+        "write extendedFeatureEdgeMesh vtk files"
+    );
+    argList::addBoolOption
+    (
+        "calcCurvature",
+        "calculate curvature and closeness fields"
+    );
+    argList::addOption
+    (
+        "closeness",
+        "scalar",
+        "span to look for surface closeness"
+    );
+    argList::addOption
+    (
+        "featureProximity",
+        "scalar",
+        "distance to look for close features"
+    );
+    argList::addBoolOption
+    (
+        "writeVTK",
+        "write surface property VTK files"
+    );
+    argList::addBoolOption
+    (
+        "manifoldEdgesOnly",
+        "remove any non-manifold (open or more than two connected faces) edges"
+    );
+
 
 #   include "setRootCase.H"
 #   include "createTime.H"
 
-    Info<< "Feature line extraction is only valid on closed manifold surfaces."
-        << endl;
+    bool writeVTK = args.optionFound("writeVTK");
 
     bool writeObj = args.optionFound("writeObj");
 
+    bool curvature = args.optionFound("curvature");
+
+    if (curvature && env("FOAM_SIGFPE"))
+    {
+        WarningIn(args.executable())
+            << "Detected floating point exception trapping (FOAM_SIGFPE)."
+            << " This might give" << nl
+            << "    problems when calculating curvature on straight angles"
+            << " (infinite curvature)" << nl
+            << "    Switch it off in case of problems." << endl;
+    }
+
+
+    Info<< "Feature line extraction is only valid on closed manifold surfaces."
+        << endl;
+
     const fileName surfFileName = args[1];
     const fileName outFileName  = args[2];
 
@@ -221,6 +476,12 @@ int main(int argc, char *argv[])
     surf.writeStats(Info);
     Info<< endl;
 
+    faceList faces(surf.size());
+
+    forAll(surf, fI)
+    {
+        faces[fI] = surf[fI].triFaceFace();
+    }
 
     // Either construct features from surface&featureangle or read set.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -334,11 +595,25 @@ int main(int argc, char *argv[])
         );
     }
 
+    if (args.optionFound("manifoldEdgesOnly"))
+    {
+        Info<< "Removing all non-manifold edges" << endl;
+
+        forAll(edgeStat, edgeI)
+        {
+            if (surf.edgeFaces()[edgeI].size() != 2)
+            {
+                edgeStat[edgeI] = surfaceFeatures::NONE;
+            }
+        }
+    }
+
+
     surfaceFeatures newSet(surf);
     newSet.setFromStatus(edgeStat);
 
-    Info<< endl << "Writing trimmed features to " << outFileName << endl;
-    newSet.write(outFileName);
+    //Info<< endl << "Writing trimmed features to " << outFileName << endl;
+    //newSet.write(outFileName);
 
     // Info<< endl << "Writing edge objs." << endl;
     // newSet.writeObj("final");
@@ -397,6 +672,358 @@ int main(int argc, char *argv[])
         bfeMesh.regIOobject::write();
     }
 
+    triSurfaceMesh searchSurf
+    (
+        IOobject
+        (
+            sFeatFileName + ".closeness",
+            runTime.constant(),
+            "extendedFeatureEdgeMesh",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf
+    );
+
+    if (!curvature)
+    {
+        Info<< "End\n" << endl;
+
+        return 0;
+    }
+
+    // Find close features
+
+    // // Dummy trim operation to mark features
+    // labelList featureEdgeIndexing = newSet.trimFeatures(-GREAT, 0);
+
+    // scalarField surfacePtFeatureIndex(surf.points().size(), -1);
+
+    // forAll(newSet.featureEdges(), eI)
+    // {
+    //     const edge& e = surf.edges()[newSet.featureEdges()[eI]];
+
+    //     surfacePtFeatureIndex[surf.meshPoints()[e.start()]] =
+    //     featureEdgeIndexing[newSet.featureEdges()[eI]];
+
+    //     surfacePtFeatureIndex[surf.meshPoints()[e.end()]] =
+    //     featureEdgeIndexing[newSet.featureEdges()[eI]];
+    // }
+
+    // if (writeVTK)
+    // {
+    //     vtkSurfaceWriter().write
+    //     (
+    //         runTime.constant()/"triSurface",    // outputDir
+    //         sFeatFileName,                      // surfaceName
+    //         surf.points(),
+    //         faces,
+    //         "surfacePtFeatureIndex",            // fieldName
+    //         surfacePtFeatureIndex,
+    //         true,                               // isNodeValues
+    //         true                                // verbose
+    //     );
+    // }
+
+    // Random rndGen(343267);
+
+    // treeBoundBox surfBB
+    // (
+    //     treeBoundBox(searchSurf.bounds()).extend(rndGen, 1e-4)
+    // );
+
+    // surfBB.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    // surfBB.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+    // indexedOctree<treeDataEdge> ftEdTree
+    // (
+    //     treeDataEdge
+    //     (
+    //         false,
+    //         surf.edges(),
+    //         surf.localPoints(),
+    //         newSet.featureEdges()
+    //     ),
+    //     surfBB,
+    //     8,      // maxLevel
+    //     10,     // leafsize
+    //     3.0     // duplicity
+    // );
+
+    // labelList nearPoints = ftEdTree.findBox
+    // (
+    //     treeBoundBox
+    //     (
+    //         sPt - featureSearchSpan*Foam::vector::one,
+    //         sPt + featureSearchSpan*Foam::vector::one
+    //     )
+    // );
+
+    Info<< "Examine curvature, feature proximity and internal and "
+        << "external closeness." << endl;
+
+    // Internal and external closeness
+
+    // Prepare start and end points for intersection tests
+
+    const vectorField& normals = searchSurf.faceNormals();
+
+    scalar span = searchSurf.bounds().mag();
+
+    args.optionReadIfPresent("closeness", span);
+
+    scalar externalAngleTolerance = 10;
+    scalar externalToleranceCosAngle = Foam::cos
+    (
+        degToRad(180 - externalAngleTolerance)
+    );
+
+    scalar internalAngleTolerance = 45;
+    scalar internalToleranceCosAngle = Foam::cos
+    (
+        degToRad(180 - internalAngleTolerance)
+    );
+
+    Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle << nl
+        << "internalToleranceCosAngle: " << internalToleranceCosAngle
+        << endl;
+
+    // Info<< "span " << span << endl;
+
+    pointField start = searchSurf.faceCentres() - span*normals;
+    pointField end = searchSurf.faceCentres() + span*normals;
+    const pointField& faceCentres = searchSurf.faceCentres();
+
+    List<List<pointIndexHit> > allHitInfo;
+
+    // Find all intersections (in order)
+    searchSurf.findLineAll(start, end, allHitInfo);
+
+    scalarField internalCloseness(start.size(), GREAT);
+    scalarField externalCloseness(start.size(), GREAT);
+
+    forAll(allHitInfo, fI)
+    {
+        const List<pointIndexHit>& hitInfo = allHitInfo[fI];
+
+        if (hitInfo.size() < 1)
+        {
+            drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
+
+            // FatalErrorIn(args.executable())
+            //     << "findLineAll did not hit its own face."
+            //     << exit(FatalError);
+        }
+        else if (hitInfo.size() == 1)
+        {
+            if (!hitInfo[0].hit())
+            {
+                // FatalErrorIn(args.executable())
+                //     << "findLineAll did not hit any face."
+                //     << exit(FatalError);
+            }
+            else if (hitInfo[0].index() != fI)
+            {
+                drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
+
+                // FatalErrorIn(args.executable())
+                //     << "findLineAll did not hit its own face."
+                //     << exit(FatalError);
+            }
+        }
+        else
+        {
+            label ownHitI = -1;
+
+            forAll(hitInfo, hI)
+            {
+                // Find the hit on the triangle that launched the ray
+
+                if (hitInfo[hI].index() == fI)
+                {
+                    ownHitI = hI;
+
+                    break;
+                }
+            }
+
+            if (ownHitI < 0)
+            {
+                drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
+
+                // FatalErrorIn(args.executable())
+                //     << "findLineAll did not hit its own face."
+                //     << exit(FatalError);
+            }
+            else if (ownHitI == 0)
+            {
+                // There are no internal hits, the first hit is the closest
+                // external hit
+
+                if
+                (
+                    (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
+                  < externalToleranceCosAngle
+                )
+                {
+                    externalCloseness[fI] = mag
+                    (
+                        faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
+                    );
+                }
+            }
+            else if (ownHitI == hitInfo.size() - 1)
+            {
+                // There are no external hits, the last but one hit is the
+                // closest internal hit
+
+                if
+                (
+                    (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
+                  < internalToleranceCosAngle
+                )
+                {
+                    internalCloseness[fI] = mag
+                    (
+                        faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
+                    );
+                }
+            }
+            else
+            {
+                if
+                (
+                    (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
+                  < externalToleranceCosAngle
+                )
+                {
+                    externalCloseness[fI] = mag
+                    (
+                        faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
+                    );
+                }
+
+                if
+                (
+                    (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
+                  < internalToleranceCosAngle
+                )
+                {
+                    internalCloseness[fI] = mag
+                    (
+                        faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
+                    );
+                }
+            }
+        }
+    }
+
+    triSurfaceScalarField internalClosenessField
+    (
+        IOobject
+        (
+            sFeatFileName + ".internalCloseness",
+            runTime.constant(),
+            "extendedFeatureEdgeMesh",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        internalCloseness
+    );
+
+    internalClosenessField.write();
+
+    triSurfaceScalarField externalClosenessField
+    (
+        IOobject
+        (
+            sFeatFileName + ".externalCloseness",
+            runTime.constant(),
+            "extendedFeatureEdgeMesh",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        externalCloseness
+    );
+
+    externalClosenessField.write();
+
+    scalarField k = calcCurvature(surf);
+
+    // Modify the curvature values on feature edges and points to be zero.
+
+    forAll(newSet.featureEdges(), fEI)
+    {
+        const edge& e = surf.edges()[newSet.featureEdges()[fEI]];
+
+        k[surf.meshPoints()[e.start()]] = 0.0;
+        k[surf.meshPoints()[e.end()]] = 0.0;
+    }
+
+    triSurfacePointScalarField kField
+    (
+        IOobject
+        (
+            sFeatFileName + ".curvature",
+            runTime.constant(),
+            "extendedFeatureEdgeMesh",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        k
+    );
+
+    kField.write();
+
+    if (writeVTK)
+    {
+        vtkSurfaceWriter().write
+        (
+            runTime.constant()/"triSurface",    // outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "internalCloseness",                // fieldName
+            internalCloseness,
+            false,                              // isNodeValues
+            true                                // verbose
+        );
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constant()/"triSurface",    // outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "externalCloseness",                // fieldName
+            externalCloseness,
+            false,                              // isNodeValues
+            true                                // verbose
+        );
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constant()/"triSurface",    // outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "curvature",                        // fieldName
+            k,
+            true,                               // isNodeValues
+            true                                // verbose
+        );
+    }
+
     Info<< "End\n" << endl;
 
     return 0;

applications/utilities/surface/surfaceSplitByTopology/Make/files

@@ -0,0 +1,4 @@
+
+surfaceSplitByTopology.C
+
+EXE = $(FOAM_APPBIN)/surfaceSplitByTopology

applications/utilities/surface/surfaceSplitByTopology/Make/options

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

applications/utilities/surface/surfaceSplitByTopology/surfaceSplitByTopology.C

@@ -0,0 +1,220 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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/>.
+
+Description
+
+    Strips any baffle parts of a surface.  A baffle region is one which is
+    reached by walking from an open edge, and stopping when a multiply connected
+    edge is reached.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "triSurface.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    argList::validArgs.clear();
+    argList::validArgs.append("input surface file");
+    argList::validArgs.append("output surface file");
+    argList args(argc, argv);
+
+    fileName surfFileName(args.additionalArgs()[0]);
+    Info<< "Reading surface from " << surfFileName << endl;
+
+    fileName outFileName(args.additionalArgs()[1]);
+    fileName outFileBaseName = outFileName.lessExt();
+    word outExtension = outFileName.ext();
+
+    // Load surface
+    triSurface surf(surfFileName);
+
+    bool anyZoneRemoved = false;
+
+    label iterationNo = 0;
+    label iterationLimit = 10;
+
+    Info<< "Splitting off baffle parts " << endl;
+
+    do
+    {
+        anyZoneRemoved = false;
+
+        labelList faceZone;
+
+        const labelListList& edFaces = surf.edgeFaces();
+        const labelListList& faceEds = surf.faceEdges();
+
+        boolList multipleEdges(edFaces.size(), false);
+
+        forAll(multipleEdges, i)
+        {
+            if (edFaces[i].size() > 2)
+            {
+                multipleEdges[i] = true;
+            }
+        }
+
+        label nZones = surf.markZones(multipleEdges, faceZone);
+
+        if (nZones < 2)
+        {
+            break;
+        }
+
+        boolList nonBaffle(faceZone.size(), true);
+        boolList baffle(faceZone.size(), true);
+        labelList pointMap;
+        labelList faceMap;
+
+
+        for (label z = 0; z < nZones; z++)
+        {
+            bool keepZone = true;
+
+            forAll(faceZone, f)
+            {
+                if (faceZone[f] == z)
+                {
+                    forAll (faceEds[f], fe)
+                    {
+                        if (edFaces[faceEds[f][fe]].size() < 2)
+                        {
+                            keepZone = false;
+
+                            anyZoneRemoved = true;
+
+                            break;
+                        }
+                    }
+                }
+
+                if (!keepZone)
+                {
+                    break;
+                }
+            }
+
+            forAll(faceZone, f)
+            {
+                if (faceZone[f] == z)
+                {
+                    nonBaffle[f] = keepZone;
+                    baffle[f] = !keepZone;
+                }
+            }
+        }
+
+        Info<< "    Iteration " << iterationNo << endl;
+
+        triSurface baffleSurf = surf.subsetMesh(baffle, pointMap, faceMap);
+
+        if (baffleSurf.size())
+        {
+            fileName bafflePartFileName =
+            outFileBaseName
+          + "_bafflePart_"
+          + name(iterationNo)
+          + "." + outExtension;
+
+            Info<< "    Writing baffle part to " << bafflePartFileName << endl;
+
+            baffleSurf.write(bafflePartFileName);
+        }
+
+        surf = surf.subsetMesh(nonBaffle, pointMap, faceMap);
+
+        if (iterationNo == iterationLimit)
+        {
+            WarningIn("surfaceSplitByTopology")
+            << "Iteration limit of " << iterationLimit << "reached" << endl;
+        }
+
+        iterationNo++;
+
+    } while (anyZoneRemoved && iterationNo < iterationLimit);
+
+    Info<< "Writing new surface to " << outFileName << endl;
+
+    surf.write(outFileName);
+
+    labelList faceZone;
+
+    const labelListList& edFaces = surf.edgeFaces();
+
+    boolList multipleEdges(edFaces.size(), false);
+
+    forAll(multipleEdges, i)
+    {
+        if (edFaces[i].size() > 2)
+        {
+            multipleEdges[i] = true;
+        }
+    }
+
+    label nZones = surf.markZones(multipleEdges, faceZone);
+
+    Info<< "Splitting remaining multiply connected parts" << endl;
+
+    for (label z = 0; z < nZones; z++)
+    {
+
+        boolList include(faceZone.size(), false);
+        labelList pointMap;
+        labelList faceMap;
+
+        forAll(faceZone, f)
+        {
+            if (faceZone[f] == z)
+            {
+                include[f] = true;
+            }
+        }
+
+        triSurface zoneSurf = surf.subsetMesh(include, pointMap, faceMap);
+
+
+        fileName remainingPartFileName =
+            outFileBaseName
+          + "_multiplePart_"
+          + name(z)
+          + "." + outExtension;
+
+        Info<< "    Writing mulitple part "
+            << z << " to " << remainingPartFileName << endl;
+
+        zoneSurf.write(remainingPartFileName);
+    }
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //