Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

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

@@ -20,7 +20,8 @@ EXE_LIBS = \
     $(CGAL_LIBS) \
     -lconformalVoronoiMesh \
     -lmeshTools \
-    -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
+    -ldecompositionMethods \
+    -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
     -ledgeMesh \
     -ltriSurface \
     -ldynamicMesh

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

@@ -19,7 +19,6 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lmeshTools \
-    -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
     -ledgeMesh \
     -ltriSurface \
     -ldynamicMesh

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C

@@ -63,7 +63,7 @@ Foam::fieldFromFile::fieldFromFile
 
 Foam::triSurfaceScalarField Foam::fieldFromFile::load()
 {
-    Info<< "Loading: " << fileName_ << endl;
+    Info<< indent << "Loading: " << fileName_ << endl;
 
     triSurfaceScalarField surfaceCellSize
     (

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/surfaceCellSizeFunction/surfaceCellSizeFunction.C

@@ -46,7 +46,10 @@ Foam::surfaceCellSizeFunction::surfaceCellSizeFunction
     dictionary(surfaceCellSizeFunctionDict),
     surface_(surface),
     coeffsDict_(subDict(type + "Coeffs")),
-    refinementFactor_(readScalar(lookup("refinementFactor")))
+    refinementFactor_
+    (
+        lookupOrDefault<scalar>("refinementFactor", 1.0)
+    )
 {}
 
 

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

@@ -44,6 +44,142 @@ const Foam::scalar Foam::conformalVoronoiMesh::tolParallel = 1e-3;
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
+Foam::scalar Foam::conformalVoronoiMesh::requiredSize
+(
+    const Foam::point& pt
+) const
+{
+    pointIndexHit surfHit;
+    label hitSurface;
+
+    DynamicList<scalar> cellSizeHits;
+
+    geometryToConformTo_.findSurfaceNearest
+    (
+        pt,
+        sqr(GREAT),
+        surfHit,
+        hitSurface
+    );
+
+    if (!surfHit.hit())
+    {
+        FatalErrorIn
+        (
+            "Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment"
+        )
+            << "findSurfaceNearest did not find a hit across the surfaces."
+            << exit(FatalError) << endl;
+    }
+
+    cellSizeHits.append(cellSizeControl().cellSize(pt));
+
+    // Primary alignment
+
+    vectorField norm(1);
+
+    allGeometry_[hitSurface].getNormal
+    (
+        List<pointIndexHit>(1, surfHit),
+        norm
+    );
+
+    const vector np = norm[0];
+
+    // Generate equally spaced 'spokes' in a circle normal to the
+    // direction from the vertex to the closest point on the surface
+    // and look for a secondary intersection.
+
+    const vector d = surfHit.hitPoint() - pt;
+
+    const tensor Rp = rotationTensor(vector(0,0,1), np);
+
+    const label s = cvMeshControls().alignmentSearchSpokes();
+
+    const scalar spanMag = geometryToConformTo_.globalBounds().mag();
+
+    scalar totalDist = 0;
+
+    for (label i = 0; i < s; i++)
+    {
+        vector spoke
+        (
+            Foam::cos(i*constant::mathematical::twoPi/s),
+            Foam::sin(i*constant::mathematical::twoPi/s),
+            0
+        );
+
+        spoke *= spanMag;
+
+        spoke = Rp & spoke;
+
+        pointIndexHit spokeHit;
+
+        label spokeSurface = -1;
+
+        // internal spoke
+
+        geometryToConformTo_.findSurfaceNearestIntersection
+        (
+            pt,
+            pt + spoke,
+            spokeHit,
+            spokeSurface
+        );
+
+        if (spokeHit.hit())
+        {
+            const Foam::point& hitPt = spokeHit.hitPoint();
+
+            scalar spokeHitDistance = mag(hitPt - pt);
+
+            cellSizeHits.append
+            (
+                cellSizeControl().cellSize(hitPt)
+            );
+
+            totalDist += spokeHitDistance;
+        }
+
+        //external spoke
+
+        Foam::point mirrorPt = pt + 2*d;
+
+        geometryToConformTo_.findSurfaceNearestIntersection
+        (
+            mirrorPt,
+            mirrorPt + spoke,
+            spokeHit,
+            spokeSurface
+        );
+
+        if (spokeHit.hit())
+        {
+            const Foam::point& hitPt = spokeHit.hitPoint();
+
+            scalar spokeHitDistance = mag(hitPt - mirrorPt);
+
+            cellSizeHits.append
+            (
+                cellSizeControl().cellSize(hitPt)
+            );
+
+            totalDist += spokeHitDistance;
+        }
+    }
+
+    scalar cellSize = 0;
+
+    forAll(cellSizeHits, hitI)
+    {
+        cellSize += cellSizeHits[hitI];
+    }
+
+    return cellSize/cellSizeHits.size();
+    //return cellSizeControl().cellSize(pt);
+}
+
+
 Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
 (
     const Foam::point& pt
@@ -828,10 +964,7 @@ void Foam::conformalVoronoiMesh::storeSizesAndAlignments
     {
         sizeAndAlignmentLocations_[i] = topoint(*pit);
 
-        storedSizes_[i] = cellSizeControl().cellSize
+        storedSizes_[i] = requiredSize(sizeAndAlignmentLocations_[i]);
-        (
-            sizeAndAlignmentLocations_[i]
-        );
 
         storedAlignments_[i] = requiredAlignment(sizeAndAlignmentLocations_[i]);
 

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -268,6 +268,9 @@ private:
         //- Return the local maximum surface protrusion distance
         inline scalar maxSurfaceProtrusion(const Foam::point& pt) const;
 
+        //- Return the required cell size at the given location
+        scalar requiredSize(const Foam::point& pt) const;
+
         //- Return the required alignment directions at the given location
         tensor requiredAlignment(const Foam::point& pt) const;
 
@@ -533,6 +536,11 @@ private:
             const Delaunay::Finite_edges_iterator& eit
         ) const;
 
+        boolList dualFaceBoundaryPoints
+        (
+            const Delaunay::Finite_edges_iterator& eit
+        ) const;
+
         //- Finds the maximum filterCount of the dual vertices
         //  (Delaunay cells) that form the dual face produced by the
         //  supplied edge
@@ -839,6 +847,13 @@ private:
             const Delaunay::Finite_facets_iterator& fit
         ) const;
 
+        //- Merge adjacent edges that are not attached to other faces
+        label mergeNearlyParallelEdges
+        (
+            const pointField& pts,
+            const scalar maxCosAngle
+        );
+
         //- Merge vertices that are very close together
         void mergeCloseDualVertices
         (

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

@@ -212,6 +212,20 @@ void Foam::conformalVoronoiMesh::calcDualMesh
                     Info<< nl << "Collapsing unnecessary faces" << endl;
 
                     collapseFaces(points, boundaryPts, deferredCollapseFaces);
+
+                    const scalar maxCosAngle
+                        = cos(degToRad(cvMeshControls().edgeMergeAngle()));
+
+                    Info<< nl << "Merging adjacent edges which have an angle "
+                        << "of greater than "
+                        << cvMeshControls().edgeMergeAngle() << ": " << endl;
+
+                    label nRemovedEdges =
+                        mergeNearlyParallelEdges(points, maxCosAngle);
+
+                    reduce(nRemovedEdges, sumOp<label>());
+
+                    Info<< "    Merged " << nRemovedEdges << " edges" << endl;
                 }
 
                 labelHashSet wrongFaces = checkPolyMeshQuality(points);
@@ -652,6 +666,113 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
 }
 
 
+Foam::label Foam::conformalVoronoiMesh::mergeNearlyParallelEdges
+(
+    const pointField& pts,
+    const scalar maxCosAngle
+)
+{
+    List<HashSet<label> > pointFaceCount(number_of_cells());
+    labelList pointNeighbour(number_of_cells(), -1);
+
+    Map<label> dualPtIndexMap;
+
+    for
+    (
+        Delaunay::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_handle c = eit->first;
+        Vertex_handle vA = c->vertex(eit->second);
+        Vertex_handle vB = c->vertex(eit->third);
+
+        if (isBoundaryDualFace(eit))
+        {
+            const face f = buildDualFace(eit);
+
+            forAll(f, pI)
+            {
+                const label pIndex = f[pI];
+
+                const label prevPointI = f.prevLabel(pI);
+                const label nextPointI = f.nextLabel(pI);
+
+                pointFaceCount[pIndex].insert(prevPointI);
+                pointFaceCount[pIndex].insert(nextPointI);
+                pointNeighbour[pIndex] = nextPointI;
+            }
+        }
+        else if
+        (
+            vA->internalOrBoundaryPoint()
+         || vB->internalOrBoundaryPoint()
+        )
+        {
+            const face f = buildDualFace(eit);
+            const boolList faceBoundaryPoints = dualFaceBoundaryPoints(eit);
+
+            forAll(f, pI)
+            {
+                const label pIndex = f[pI];
+
+                const label prevPointI = f.prevLabel(pI);
+                const label nextPointI = f.nextLabel(pI);
+
+                pointFaceCount[pIndex].insert(prevPointI);
+                pointFaceCount[pIndex].insert(nextPointI);
+
+                if (faceBoundaryPoints[pI] == false)
+                {
+                    pointNeighbour[pIndex] = nextPointI;
+                }
+                else
+                {
+                    if (faceBoundaryPoints[prevPointI] == true)
+                    {
+                        pointNeighbour[pIndex] = prevPointI;
+                    }
+                    else if (faceBoundaryPoints[nextPointI] == true)
+                    {
+                        pointNeighbour[pIndex] = nextPointI;
+                    }
+                    else
+                    {
+                        pointNeighbour[pIndex] = pIndex;
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(pointFaceCount, pI)
+    {
+        if (pointFaceCount[pI].size() == 2)
+        {
+            List<vector> edges(2, vector(0, 0, 0));
+
+            label count = 0;
+            forAllConstIter(HashSet<label>, pointFaceCount[pI], iter)
+            {
+                edges[count] = pts[pI] - pts[iter.key()];
+                edges[count] /= mag(edges[count]) + VSMALL;
+                count++;
+            }
+
+            if (mag(edges[0] & edges[1]) > maxCosAngle)
+            {
+                dualPtIndexMap.insert(pI, pointNeighbour[pI]);
+            }
+        }
+    }
+
+    reindexDualVertices(dualPtIndexMap);
+
+    return dualPtIndexMap.size();
+}
+
+
 void Foam::conformalVoronoiMesh::smoothSurface
 (
     pointField& pts,
@@ -696,7 +817,6 @@ void Foam::conformalVoronoiMesh::smoothSurface
     } while (nCollapsedFaces > 0);
 
     // Force all points of boundary faces to be on the surface
-
     for
     (
         Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -917,6 +1037,7 @@ void Foam::conformalVoronoiMesh::collapseFaces
 
         mergeCloseDualVertices(pts, boundaryPts);
 
+
         if (nCollapsedFaces > 0)
         {
             Info<< "    Collapsed " << nCollapsedFaces << " faces" << endl;
@@ -932,6 +1053,7 @@ void Foam::conformalVoronoiMesh::collapseFaces
         }
 
     } while (nCollapsedFaces > 0);
+
 }
 
 
@@ -1064,6 +1186,25 @@ Foam::conformalVoronoiMesh::collapseFace
 
     bool allowEarlyCollapseToPoint = true;
 
+
+//    // Quick exit
+//    label smallEdges = 0;
+//    const edgeList& fEdges = f.edges();
+//    forAll(fEdges, eI)
+//    {
+//        const edge& e = fEdges[eI];
+//
+//        if (e.mag(pts) < 0.2*targetFaceSize)
+//        {
+//            smallEdges++;
+//        }
+//    }
+//    if (smallEdges == 0)
+//    {
+//        return fcmNone;
+//    }
+
+
     // if (maxFC > cvMeshControls().filterCountSkipThreshold() - 3)
     // {
     //     limitToQuadsOrTris = true;
@@ -1071,6 +1212,7 @@ Foam::conformalVoronoiMesh::collapseFace
     //     allowEarlyCollapseToPoint = false;
     // }
 
+
     collapseSizeLimitCoeff *= pow
     (
         cvMeshControls().filterErrorReductionCoeff(),
@@ -1158,6 +1300,91 @@ Foam::conformalVoronoiMesh::collapseFace
         }
     }
 
+
+    scalar maxDist = 0;
+    scalar minDist = GREAT;
+//
+//    if (f.size() <= 3)
+//    {
+//        const edgeList& fEdges = f.edges();
+//
+//        forAll(fEdges, eI)
+//        {
+//            const edge& e = fEdges[eI];
+//            const scalar d = e.mag(pts);
+//
+//            if (d > maxDist)
+//            {
+//                maxDist = d;
+//                collapseAxis = e.vec(pts);
+//            }
+//            else if (d < minDist && d != 0)
+//            {
+//                minDist = d;
+//            }
+//        }
+//    }
+//    else
+//    {
+//        forAll(f, pI)
+//        {
+//            for (label i = pI + 1; i < f.size(); ++i)
+//            {
+//                if
+//                (
+//                    f[i] != f.nextLabel(pI)
+//                 && f[i] != f.prevLabel(pI)
+//                )
+//                {
+//                    scalar d = mag(pts[f[pI]] - pts[f[i]]);
+//
+//                    if (d > maxDist)
+//                    {
+//                        maxDist = d;
+//                        collapseAxis = pts[f[pI]] - pts[f[i]];
+//                    }
+//                    else if (d < minDist && d != 0)
+//                    {
+//                        minDist = d;
+//                    }
+//                }
+//            }
+//        }
+//    }
+//
+//    const edgeList& fEdges = f.edges();
+//
+//    scalar perimeter = 0;
+//
+//    forAll(fEdges, eI)
+//    {
+//        const edge& e = fEdges[eI];
+//        const scalar d = e.mag(pts);
+//
+//        perimeter += d;
+//
+////        collapseAxis += e.vec(pts);
+//
+//        if (d > maxDist)
+//        {
+//            collapseAxis = e.vec(pts);
+//            maxDist = d;
+//        }
+//        else if (d < minDist && d != 0)
+//        {
+//            minDist = d;
+//        }
+//    }
+//
+//    collapseAxis /= mag(collapseAxis);
+//
+////    Info<< f.size() << " " << minDist << " " << maxDist << " | "
+////        << collapseAxis << endl;
+//
+//    aspectRatio = maxDist/minDist;
+//
+//    aspectRatio = min(aspectRatio, sqr(perimeter)/(16.0*fA));
+
     if (magSqr(collapseAxis) < VSMALL)
     {
         WarningIn
@@ -1170,9 +1397,9 @@ Foam::conformalVoronoiMesh::collapseFace
         // Output face and collapse axis for visualisation
 
         Pout<< "# Aspect ratio = " << aspectRatio << nl
-            << "# inertia = " << J << nl
+//            << "# inertia = " << J << nl
-            << "# determinant = " << detJ << nl
+//            << "# determinant = " << detJ << nl
-            << "# eigenvalues = " << eigenValues(J) << nl
+//            << "# eigenvalues = " << eigenValues(J) << nl
             << "# collapseAxis = " << collapseAxis << nl
             << "# facePts = " << facePts << nl
             << endl;
@@ -1280,8 +1507,6 @@ Foam::conformalVoronoiMesh::collapseFace
     {
         scalar guardFraction = cvMeshControls().edgeCollapseGuardFraction();
 
-        cvMeshControls().maxCollapseFaceToPointSideLengthCoeff();
-
         if
         (
             allowEarlyCollapseToPoint
@@ -1337,6 +1562,8 @@ Foam::conformalVoronoiMesh::collapseFace
             Foam::point collapseToPt =
                 collapseAxis*(sum(dNeg)/dNeg.size() - dShift) + fC;
 
+//            DynamicList<label> faceBoundaryPts(f.size());
+
             forAll(facePtsNeg, fPtI)
             {
                 if (boundaryPts[facePtsNeg[fPtI]] == true)
@@ -1351,9 +1578,36 @@ Foam::conformalVoronoiMesh::collapseFace
                     collapseToPt = pts[collapseToPtI];
 
                     break;
+
+//                    faceBoundaryPts.append(facePtsNeg[fPtI]);
                 }
             }
 
+//            if (!faceBoundaryPts.empty())
+//            {
+//                if (faceBoundaryPts.size() == 2)
+//                {
+//                    collapseToPtI = faceBoundaryPts[0];
+//
+//                    collapseToPt =
+//                        0.5
+//                       *(
+//                            pts[faceBoundaryPts[0]]
+//                          + pts[faceBoundaryPts[1]]
+//                        );
+//                }
+//                else if (faceBoundaryPts.size() < f.size())
+//                {
+//                    face bFace(faceBoundaryPts);
+//
+//                    collapseToPtI = faceBoundaryPts.first();
+//
+//                    collapseToPt = bFace.centre(pts);
+//                }
+//            }
+//
+//            faceBoundaryPts.clear();
+
             // ...otherwise arbitrarily choosing the most distant
             // point as the index to collapse to.
 
@@ -1384,9 +1638,34 @@ Foam::conformalVoronoiMesh::collapseFace
                     collapseToPt = pts[collapseToPtI];
 
                     break;
+
+//                    faceBoundaryPts.append(facePtsNeg[fPtI]);
                 }
             }
 
+//            if (!faceBoundaryPts.empty())
+//            {
+//                if (faceBoundaryPts.size() == 2)
+//                {
+//                    collapseToPtI = faceBoundaryPts[0];
+//
+//                    collapseToPt =
+//                        0.5
+//                       *(
+//                            pts[faceBoundaryPts[0]]
+//                          + pts[faceBoundaryPts[1]]
+//                        );
+//                }
+//                else if (faceBoundaryPts.size() < f.size())
+//                {
+//                    face bFace(faceBoundaryPts);
+//
+//                    collapseToPtI = faceBoundaryPts.first();
+//
+//                    collapseToPt = bFace.centre(pts);
+//                }
+//            }
+
             // ...otherwise arbitrarily choosing the most distant
             // point as the index to collapse to.
 
@@ -1406,6 +1685,8 @@ Foam::conformalVoronoiMesh::collapseFace
 
             Foam::point collapseToPt = fC;
 
+            DynamicList<label> faceBoundaryPts(f.size());
+
             forAll(facePts, fPtI)
             {
                 if (boundaryPts[facePts[fPtI]] == true)
@@ -1415,11 +1696,32 @@ Foam::conformalVoronoiMesh::collapseFace
                     // use the first boundary point encountered if
                     // there are multiple boundary points.
 
-                    collapseToPtI = facePts[fPtI];
+//                    collapseToPtI = facePts[fPtI];
+//
+//                    collapseToPt = pts[collapseToPtI];
+//
+//                    break;
 
-                    collapseToPt = pts[collapseToPtI];
+                    faceBoundaryPts.append(facePts[fPtI]);
+                }
+            }
 
-                    break;
+            if (!faceBoundaryPts.empty())
+            {
+                if (faceBoundaryPts.size() == 2)
+                {
+                    collapseToPtI = faceBoundaryPts[0];
+
+                    collapseToPt =
+                        0.5*(pts[faceBoundaryPts[0]] + pts[faceBoundaryPts[1]]);
+                }
+                else if (faceBoundaryPts.size() < f.size())
+                {
+                    face bFace(faceBoundaryPts);
+
+                    collapseToPtI = faceBoundaryPts.first();
+
+                    collapseToPt = bFace.centre(pts);
                 }
             }
 
@@ -1791,9 +2093,9 @@ void Foam::conformalVoronoiMesh::checkCellSizing()
         }
 
         Info<< "    Automatically re-sizing " << cellsToResize.size()
-            << " cells that are attached to the bad faces: DISABLED" << endl;
+            << " cells that are attached to the bad faces: " << endl;
 
-        //cellSizeControl_.setCellSizes(cellsToResize);
+        cellSizeControl_.setCellSizes(cellsToResize);
     }
 
     timeCheck("End of Cell Sizing");
@@ -1858,7 +2160,10 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findOffsetPatchFaces
         }
     }
 
+    if (cvMeshControls().objOutput())
+    {
         offsetBoundaryCells.write();
+    }
 
     return offsetBoundaryCells;
 }
@@ -2092,13 +2397,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
 
             pts[dualVertI] = cit->dual();
 
-            if
+            if (cit->boundaryDualVertex())
-            (
-                !cit->vertex(0)->internalOrBoundaryPoint()
-             || !cit->vertex(1)->internalOrBoundaryPoint()
-             || !cit->vertex(2)->internalOrBoundaryPoint()
-             || !cit->vertex(3)->internalOrBoundaryPoint()
-            )
             {
                 // This is a boundary dual vertex
                 boundaryPts[dualVertI] = true;

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

@@ -719,12 +719,15 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
 
                         const Foam::point& p = infoList[hitI].hitPoint();
 
-                        const scalar separationDistance
+                        const scalar separationDistance =
-                            = mag(p - info.hitPoint());
+                            mag(p - info.hitPoint());
 
-                        const scalar minSepDist
+                        const scalar minSepDist =
-                            = sqr(cvMeshControls().removalDistCoeff()
+                            sqr
-                             *targetCellSize(p));
+                            (
+                                cvMeshControls().removalDistCoeff()
+                               *targetCellSize(p)
+                            );
 
                         // Reject the point if it is too close to another
                         // surface point.

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H

@@ -366,6 +366,49 @@ inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
 }
 
 
+inline Foam::List<bool> Foam::conformalVoronoiMesh::dualFaceBoundaryPoints
+(
+    const Delaunay::Finite_edges_iterator& eit
+) const
+{
+    Cell_circulator ccStart = incident_cells(*eit);
+    Cell_circulator cc1 = ccStart;
+    Cell_circulator cc2 = cc1;
+
+    // Advance the second circulator so that it always stays on the next
+    // cell around the edge;
+    cc2++;
+
+    DynamicList<bool> tmpFaceBoundaryPoints;
+
+    do
+    {
+        label cc1I = cc1->cellIndex();
+
+        label cc2I = cc2->cellIndex();
+
+        if (cc1I != cc2I)
+        {
+            if (cc1->boundaryDualVertex())
+            {
+                tmpFaceBoundaryPoints.append(true);
+            }
+            else
+            {
+                tmpFaceBoundaryPoints.append(false);
+            }
+        }
+
+        cc1++;
+
+        cc2++;
+
+    } while (cc1 != ccStart);
+
+    return tmpFaceBoundaryPoints;
+}
+
+
 inline Foam::List<Foam::label> Foam::conformalVoronoiMesh::processorsAttached
 (
     const Delaunay::Finite_facets_iterator& fit

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

@@ -291,6 +291,11 @@ Foam::cvControls::cvControls
         filteringDict.lookup("mergeClosenessCoeff")
     );
 
+    edgeMergeAngle_ = readScalar
+    (
+        filteringDict.lookup("edgeMergeAngle")
+    );
+
     continueFilteringOnBadInitialPolyMesh_ = Switch
     (
         filteringDict.lookupOrDefault<Switch>

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H

@@ -234,6 +234,11 @@ class cvControls
         //  being merged, fraction of the local target cell size
         scalar mergeClosenessCoeff_;
 
+        //- If the angle between two dual edges that are connected by a single
+        //  point is less than this angle, then the edges will be merged into a
+        //  single edge.
+        scalar edgeMergeAngle_;
+
         //- If the mesh quality criteria cannot be satisfied, continue
         //  with filtering anyway?
         Switch continueFilteringOnBadInitialPolyMesh_;
@@ -405,6 +410,9 @@ public:
             //- Return the mergeClosenessCoeff
             inline scalar mergeClosenessCoeff() const;
 
+            //- Return the edgeMergeAngle
+            inline scalar edgeMergeAngle() const;
+
             //- Return the continueFilteringOnBadInitialPolyMesh Switch
             inline Switch continueFilteringOnBadInitialPolyMesh() const;
 

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H

@@ -166,6 +166,12 @@ inline Foam::scalar Foam::cvControls::mergeClosenessCoeff() const
 }
 
 
+inline Foam::scalar Foam::cvControls::edgeMergeAngle() const
+{
+    return edgeMergeAngle_;
+}
+
+
 inline Foam::Switch
 Foam::cvControls::continueFilteringOnBadInitialPolyMesh() const
 {

applications/utilities/surface/surfaceFeatureExtract/Allwmake

@@ -3,30 +3,25 @@ cd ${0%/*} || exit 1    # run from this directory
 
 set -x
 
-if [ ! -e "Make/files" ] || [ ! -e "Make/options" ]
+if [ -n "$CGAL_ARCH_PATH" ]
 then
-    mkdir -p Make
-
-    if [ -n "$CGAL_ARCH_PATH" ]
-    then
-        cp -r MakeWithCGAL/* Make
-
     echo
-        echo Compiling surfaceFeatureExtract with CGAL support for curvature
+    echo "Compiling surfaceFeatureExtract with CGAL curvature support"
     echo
 
-        wmake
+    wmake "ENABLE_CURVATURE=-DENABLE_CURVATURE \
-    else
+           EXE_FROUNDING_MATH=-frounding-math \
-        cp -r MakeWithoutCGAL/* Make
+           USE_F2C=-DCGAL_USE_F2C \
-
+           CGAL_LIBDIR=-L$CGAL_ARCH_PATH/lib \
-        echo
+           LAPACK_LIB=-llapack \
-        echo Compiling surfaceFeatureExtract without CGAL support for curvature
+           BLAS_LIB=-lblas \
-        echo
+           CGAL_LIB=-lCGAL"
-
-        wmake
-    fi
 else
-    echo surfaceFeatureExtract already has a Make folder
+    echo
+    echo "Compiling surfaceFeatureExtract without CGAL curvature support"
+    echo
+
+    wmake
 fi
 
 

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

@@ -80,18 +80,56 @@ public:
     // Constructors
 
         //- Construct with reference to triSurface
-        buildCGALPolyhedron(const triSurface& surf);
+        explicit buildCGALPolyhedron(const triSurface& surf)
+        :
+            CGAL::Modifier_base<HalfedgeDS>(),
+            surf_(surf)
+        {}
 
 
     //- Destructor
-    ~buildCGALPolyhedron();
+    ~buildCGALPolyhedron(){}
 
 
     // Member Operators
 
         //- operator() of this `modifier' called by delegate function of
         //  Polyhedron
-        void operator()(HalfedgeDS& hds);
+        void 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/Make/options

@@ -1,4 +1,11 @@
+include $(GENERAL_RULES)/CGAL
+
 EXE_INC = \
+    ${ENABLE_CURVATURE}\
+    ${EXE_FROUNDING_MATH} \
+    ${USE_F2C} \
+    ${CGAL_INC} \
+    -ICGALPolyhedron \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
@@ -7,6 +14,11 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
+    $(CGAL_LIBS) \
+    ${CGAL_LIBDIR} \
+    ${LAPACK_LIB} \
+    ${BLAS_LIB} \
+    ${CGAL_LIB} \
     -lmeshTools \
     -ledgeMesh \
     -ltriSurface \

applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files

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

applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options

@@ -1,29 +0,0 @@
-EXE_FROUNDING_MATH = -frounding-math
-EXE_NDEBUG = -DNDEBUG
-USE_F2C = -DCGAL_USE_F2C
-include $(GENERAL_RULES)/CGAL
-
-EXE_INC = \
-    -DENABLE_CURVATURE \
-    ${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)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
-    $(CGAL_LIBS) \
-    -L$(CGAL_ARCH_PATH)/lib \
-    -llapack \
-    -lblas \
-    -lCGAL \
-    -lmeshTools \
-    -ledgeMesh \
-    -ltriSurface \
-    -lsampling

applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files

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

applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options

@@ -1,13 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/edgeMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
-    -lmeshTools \
-    -ledgeMesh \
-    -ltriSurface \
-    -lsampling

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -74,30 +74,6 @@ scalarField calcCurvature(const triSurface& surf)
 
     // 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;
@@ -181,8 +157,8 @@ scalarField calcCurvature(const triSurface& surf)
 //                 << monge_fit.condition_number() << nl << std::endl;
 
         // Use the maximum curvature to give smaller cell sizes later.
-        k[vertI++]
+        k[vertI++] =
-            = max
+            max
             (
                 mag(monge_form.principal_curvatures(0)),
                 mag(monge_form.principal_curvatures(1))
@@ -198,8 +174,10 @@ bool edgesConnected(const edge& e1, const edge& e2)
 {
     if
     (
-        e1.start() == e2.start() || e1.start() == e2.end()
+        e1.start() == e2.start()
-     || e1.end() == e2.start() || e1.end() == e2.end()
+     || e1.start() == e2.end()
+     || e1.end() == e2.start()
+     || e1.end() == e2.end()
     )
     {
         return true;
@@ -289,8 +267,8 @@ scalar calcProximityOfFeatureEdges
                     // Don't refine if the edges are connected to each other
                     if (!edgesConnected(e1, e2))
                     {
-                        scalar curDist
+                        scalar curDist =
-                            = mag(pHit1.hitPoint() - pHit2.hitPoint());
+                            mag(pHit1.hitPoint() - pHit2.hitPoint());
 
                         minDist = min(curDist, minDist);
                     }
@@ -503,96 +481,77 @@ int main(int argc, char *argv[])
         "extract and write surface features to file"
     );
     argList::noParallel();
-    argList::validArgs.append("surface");
-    argList::validArgs.append("output set");
 
     argList::addOption
     (
-        "includedAngle",
+        "dict",
-        "degrees",
+        "word",
-        "construct feature set from included angle [0..180]"
+        "name of dictionary to provide feature extraction information"
     );
-    argList::addOption
-    (
-        "set",
-        "name",
-        "use existing feature set from file"
-    );
-    argList::addOption
-    (
-        "minLen",
-        "scalar",
-        "remove features shorter than the specified cumulative length"
-    );
-    argList::addOption
-    (
-        "minElem",
-        "int",
-        "remove features with fewer than the specified number of edges"
-    );
-    argList::addOption
-    (
-        "subsetBox",
-        "((x0 y0 z0)(x1 y1 z1))",
-        "remove edges outside specified bounding box"
-    );
-    argList::addOption
-    (
-        "deleteBox",
-        "((x0 y0 z0)(x1 y1 z1))",
-        "remove edges within specified bounding box"
-    );
-    argList::addBoolOption
-    (
-        "writeObj",
-        "write extendedFeatureEdgeMesh obj files"
-    );
-    argList::addBoolOption
-    (
-        "writeVTK",
-        "write extendedFeatureEdgeMesh vtk files"
-    );
-    argList::addBoolOption
-    (
-        "closeness",
-        "span to look for surface closeness"
-    );
-    argList::addOption
-    (
-        "featureProximity",
-        "scalar",
-        "distance to look for close features"
-    );
-    argList::addBoolOption
-    (
-        "manifoldEdgesOnly",
-        "remove any non-manifold (open or more than two connected faces) edges"
-    );
-    argList::addOption
-    (
-        "plane",
-        "(nx ny nz)(z0 y0 z0)",
-        "use a plane to create feature edges for 2D meshing"
-    );
-
-#   ifdef ENABLE_CURVATURE
-    argList::addBoolOption
-    (
-        "curvature",
-        "calculate curvature and closeness fields"
-    );
-#   endif
-
     
 #   include "setRootCase.H"
 #   include "createTime.H"
 
-    bool writeVTK = args.optionFound("writeVTK");
+    word dictName
+    (
+        args.optionLookupOrDefault<word>("dict", "surfaceFeatureExtractDict")
+    );
 
-    bool writeObj = args.optionFound("writeObj");
+    Info<< "Reading " << dictName << nl << endl;
 
-    bool curvature = args.optionFound("curvature");
+    IOdictionary dict
+    (
+        IOobject
+        (
+            dictName,
+            runTime.system(),
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
 
+    forAllConstIter(dictionary, dict, iter)
+    {
+        dictionary surfaceDict = dict.subDict(iter().keyword());
+
+        const fileName surfFileName = iter().keyword();
+        const fileName sFeatFileName = surfFileName.lessExt().name();
+
+        Info<< "Surface            : " << surfFileName << nl << endl;
+
+        const Switch writeVTK =
+            surfaceDict.lookupOrAddDefault<Switch>("writeVTK", "off");
+        const Switch writeObj =
+            surfaceDict.lookupOrAddDefault<Switch>("writeObj", "off");
+        const Switch writeFeatureEdgeMesh =
+            surfaceDict.lookupOrAddDefault<Switch>
+            (
+                "writeFeatureEdgeMesh",
+                "off"
+            );
+
+        const Switch curvature =
+            surfaceDict.lookupOrAddDefault<Switch>("curvature", "off");
+        const Switch featureProximity =
+            surfaceDict.lookupOrAddDefault<Switch>("featureProximity", "off");
+        const Switch closeness =
+            surfaceDict.lookupOrAddDefault<Switch>("closeness", "off");
+
+        const word extractionMethod = surfaceDict.lookup("extractionMethod");
+
+
+#ifndef ENABLE_CURVATURE
+        if (curvature)
+        {
+            WarningIn(args.executable())
+                << "Curvature calculation has been requested but "
+                << args.executable() << " has not " << nl
+                << "    been compiled with CGAL. "
+                << "Skipping the curvature calculation." << endl;
+        }
+#else
         if (curvature && env("FOAM_SIGFPE"))
         {
             WarningIn(args.executable())
@@ -602,25 +561,16 @@ int main(int argc, char *argv[])
                 << " (infinite curvature)" << nl
                 << "    Switch it off in case of problems." << endl;
         }
+#endif
 
 
-    Info<< "Feature line extraction is only valid on closed manifold surfaces."
+        Info<< nl << "Feature line extraction is only valid on closed manifold "
-        << endl;
+            << "surfaces." << endl;
-
-    const fileName surfFileName = args[1];
-    const fileName outFileName  = args[2];
-
-    Info<< "Surface            : " << surfFileName << nl
-        << "Output feature set : " << outFileName << nl
-        << endl;
-
-    fileName sFeatFileName = surfFileName.lessExt().name();
-
 
         // Read
         // ~~~~
 
-    triSurface surf(surfFileName);
+        triSurface surf("constant/triSurface/" + surfFileName);
 
         Info<< "Statistics:" << endl;
         surf.writeStats(Info);
@@ -633,38 +583,52 @@ int main(int argc, char *argv[])
             faces[fI] = surf[fI].triFaceFace();
         }
 
-    // Either construct features from surface&featureangle or read set.
+
+        // Either construct features from surface & featureAngle or read set.
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         surfaceFeatures set(surf);
 
-    if (args.optionFound("set"))
+        if (extractionMethod == "extractFromFile")
         {
-        const fileName setName = args["set"];
+            const fileName featureEdgeFile =
+                surfaceDict.subDict("extractFromFile").lookup
+                (
+                    "featureEdgeFile"
+                );
 
-        Info<< "Reading existing feature set from file " << setName << endl;
+            edgeMesh eMesh(featureEdgeFile);
 
-        set = surfaceFeatures(surf, setName);
+            // Sometimes duplicate edges are present. Remove them.
+            eMesh.mergeEdges();
+
+            Info<< nl << "Reading existing feature edges from file "
+                << featureEdgeFile << endl;
+
+            set = surfaceFeatures(surf, eMesh.points(), eMesh.edges());
         }
-    else if (args.optionFound("includedAngle"))
+        else if (extractionMethod == "extractFromSurface")
         {
-        const scalar includedAngle = args.optionRead<scalar>("includedAngle");
+            const scalar includedAngle =
+                readScalar
+                (
+                    surfaceDict.subDict("extractFromSurface").lookup
+                    (
+                        "includedAngle"
+                    )
+                );
 
-        Info<< "Constructing feature set from included angle " << includedAngle
+            Info<< nl << "Constructing feature set from included angle "
-            << endl;
+                << includedAngle << endl;
 
             set = surfaceFeatures(surf, includedAngle);
-
-        // Info<< nl << "Writing initial features" << endl;
-        // set.write("initial.fSet");
-        // set.writeObj("initial");
         }
         else
         {
             FatalErrorIn(args.executable())
                 << "No initial feature set. Provide either one"
-            << " of -set (to read existing set)" << nl
+                << " of extractFromFile (to read existing set)" << nl
-            << " or -includedAngle (to new set construct from angle)"
+                << " or extractFromSurface (to construct new set from angle)"
                 << exit(FatalError);
         }
 
@@ -682,23 +646,25 @@ int main(int argc, char *argv[])
         // Trim set
         // ~~~~~~~~
 
-    scalar minLen = -GREAT;
+        if (surfaceDict.isDict("trimFeatures"))
-    if (args.optionReadIfPresent("minLen", minLen))
         {
-        Info<< "Removing features of length < " << minLen << endl;
+            dictionary trimDict = surfaceDict.subDict("trimFeatures");
-    }
 
-    label minElem = 0;
+            scalar minLen =
-    if (args.optionReadIfPresent("minElem", minElem))
+                trimDict.lookupOrAddDefault<scalar>("minLen", -GREAT);
-    {
+
-        Info<< "Removing features with number of edges < " << minElem << endl;
+            label minElem = trimDict.lookupOrAddDefault<label>("minElem", 0);
-    }
 
             // Trim away small groups of features
             if (minElem > 0 || minLen > 0)
             {
+                Info<< "Removing features of length < "
+                    << minLen << endl;
+                Info<< "Removing features with number of edges < "
+                    << minElem << endl;
+
                 set.trimFeatures(minLen, minElem);
-        Info<< endl << "Removed small features" << endl;
+            }
         }
 
 
@@ -708,12 +674,13 @@ int main(int argc, char *argv[])
         // Convert to marked edges, points
         List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
 
-    if (args.optionFound("subsetBox"))
+        if (surfaceDict.isDict("subsetFeatures"))
         {
-        treeBoundBox bb
+            dictionary subsetDict = surfaceDict.subDict("subsetFeatures");
-        (
+
-            args.optionLookup("subsetBox")()
+            if (subsetDict.found("insideBox"))
-        );
+            {
+                treeBoundBox bb(subsetDict.lookup("insideBox")());
 
                 Info<< "Removing all edges outside bb " << bb << endl;
                 dumpBox(bb, "subsetBox.obj");
@@ -726,12 +693,9 @@ int main(int argc, char *argv[])
                     edgeStat
                 );
             }
-    else if (args.optionFound("deleteBox"))
+            else if (subsetDict.found("outsideBox"))
             {
-        treeBoundBox bb
+                treeBoundBox bb(subsetDict.lookup("outsideBox")());
-        (
-            args.optionLookup("deleteBox")()
-        );
 
                 Info<< "Removing all edges inside bb " << bb << endl;
                 dumpBox(bb, "deleteBox.obj");
@@ -745,7 +709,10 @@ int main(int argc, char *argv[])
                 );
             }
 
-    if (args.optionFound("manifoldEdgesOnly"))
+            const Switch manifoldEdges =
+                subsetDict.lookupOrAddDefault<Switch>("manifoldEdges", "no");
+
+            if (manifoldEdges)
             {
                 Info<< "Removing all non-manifold edges" << endl;
 
@@ -758,9 +725,9 @@ int main(int argc, char *argv[])
                 }
             }
 
-    if (args.optionFound("plane"))
+            if (subsetDict.found("plane"))
             {
-        plane cutPlane(args.optionLookup("plane")());
+                plane cutPlane(subsetDict.lookup("plane")());
 
                 deleteEdges
                 (
@@ -770,22 +737,23 @@ int main(int argc, char *argv[])
                 );
 
                 Info<< "Only edges that intersect the plane with normal "
-            << cutPlane.normal() << " and base point " << cutPlane.refPoint()
+                    << cutPlane.normal()
+                    << " and base point " << cutPlane.refPoint()
                     << " will be included as feature edges."<< endl;
             }
+        }
 
 
         surfaceFeatures newSet(surf);
         newSet.setFromStatus(edgeStat);
 
-    //Info<< endl << "Writing trimmed features to " << outFileName << endl;
+        if (writeObj)
-    //newSet.write(outFileName);
+        {
-
+            newSet.writeObj("final");
-    // Info<< endl << "Writing edge objs." << endl;
+        }
-    // newSet.writeObj("final");
 
         Info<< nl
-        << "Final feature set:" << nl
+            << "Final feature set after trimming and subsetting:" << nl
             << "    feature points : " << newSet.featurePoints().size() << nl
             << "    feature edges  : " << newSet.featureEdges().size() << nl
             << "    of which" << nl
@@ -802,8 +770,8 @@ int main(int argc, char *argv[])
             sFeatFileName + ".extendedFeatureEdgeMesh"
         );
 
-    Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
+        Info<< nl << "Writing extendedFeatureEdgeMesh to "
-        << endl;
+            << feMesh.objectPath() << endl;
 
         if (writeObj)
         {
@@ -813,6 +781,7 @@ int main(int argc, char *argv[])
         feMesh.write();
 
         // Write a featureEdgeMesh for backwards compatibility
+        if (writeFeatureEdgeMesh)
         {
             featureEdgeMesh bfeMesh
             (
@@ -917,10 +886,10 @@ int main(int argc, char *argv[])
     //     )
     // );
 
-    if (args.optionFound("closeness"))
+        if (closeness)
         {
-        Info<< nl << "Extracting internal and external closeness of surface."
+            Info<< nl << "Extracting internal and external closeness of "
-            << endl;
+                << "surface." << endl;
 
             // Internal and external closeness
 
@@ -930,21 +899,22 @@ int main(int argc, char *argv[])
 
             scalar span = searchSurf.bounds().mag();
 
-        //args.optionReadIfPresent("closeness", span);
-
             scalar externalAngleTolerance = 10;
-        scalar externalToleranceCosAngle = Foam::cos
+            scalar externalToleranceCosAngle =
+                Foam::cos
                 (
                     degToRad(180 - externalAngleTolerance)
                 );
 
             scalar internalAngleTolerance = 45;
-        scalar internalToleranceCosAngle = Foam::cos
+            scalar internalToleranceCosAngle =
+                Foam::cos
                 (
                     degToRad(180 - internalAngleTolerance)
                 );
 
-        Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle << nl
+            Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle
+                << nl
                 << "internalToleranceCosAngle: " << internalToleranceCosAngle
                 << endl;
 
@@ -984,7 +954,15 @@ int main(int argc, char *argv[])
                     }
                     else if (hitInfo[0].index() != fI)
                     {
-                    drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
+                        drawHitProblem
+                        (
+                            fI,
+                            surf,
+                            start,
+                            faceCentres,
+                            end,
+                            hitInfo
+                        );
 
                         // FatalErrorIn(args.executable())
                         //     << "findLineAll did not hit its own face."
@@ -1009,7 +987,15 @@ int main(int argc, char *argv[])
 
                     if (ownHitI < 0)
                     {
-                    drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
+                        drawHitProblem
+                        (
+                            fI,
+                            surf,
+                            start,
+                            faceCentres,
+                            end,
+                            hitInfo
+                        );
 
                         // FatalErrorIn(args.executable())
                         //     << "findLineAll did not hit its own face."
@@ -1017,35 +1003,45 @@ int main(int argc, char *argv[])
                     }
                     else if (ownHitI == 0)
                     {
-                    // There are no internal hits, the first hit is the closest
+                        // There are no internal hits, the first hit is the
-                    // external hit
+                        // closest external hit
 
                         if
                         (
-                        (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
+                            (
+                                normals[fI]
+                              & normals[hitInfo[ownHitI + 1].index()]
+                            )
                           < externalToleranceCosAngle
                         )
                         {
-                        externalCloseness[fI] = mag
+                            externalCloseness[fI] =
+                                mag
                                 (
-                            faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
+                                    faceCentres[fI]
+                                  - hitInfo[ownHitI + 1].hitPoint()
                                 );
                         }
                     }
                     else if (ownHitI == hitInfo.size() - 1)
                     {
-                    // There are no external hits, the last but one hit is the
+                        // There are no external hits, the last but one hit is
-                    // closest internal hit
+                        // the closest internal hit
 
                         if
                         (
-                        (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
+                            (
+                                normals[fI]
+                              & normals[hitInfo[ownHitI - 1].index()]
+                            )
                           < internalToleranceCosAngle
                         )
                         {
-                        internalCloseness[fI] = mag
+                            internalCloseness[fI] =
+                                mag
                                 (
-                            faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
+                                    faceCentres[fI]
+                                  - hitInfo[ownHitI - 1].hitPoint()
                                 );
                         }
                     }
@@ -1053,25 +1049,35 @@ int main(int argc, char *argv[])
                     {
                         if
                         (
-                        (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
+                            (
+                                normals[fI]
+                              & normals[hitInfo[ownHitI + 1].index()]
+                            )
                           < externalToleranceCosAngle
                         )
                         {
-                        externalCloseness[fI] = mag
+                            externalCloseness[fI] =
+                                mag
                                 (
-                            faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
+                                    faceCentres[fI]
+                                  - hitInfo[ownHitI + 1].hitPoint()
                                 );
                         }
 
                         if
                         (
-                        (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
+                            (
+                                normals[fI]
+                              & normals[hitInfo[ownHitI - 1].index()]
+                            )
                           < internalToleranceCosAngle
                         )
                         {
-                        internalCloseness[fI] = mag
+                            internalCloseness[fI] =
+                                mag
                                 (
-                            faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
+                                    faceCentres[fI]
+                                  - hitInfo[ownHitI - 1].hitPoint()
                                 );
                         }
                     }
@@ -1144,9 +1150,10 @@ int main(int argc, char *argv[])
 
 
 #ifdef ENABLE_CURVATURE
-    if (args.optionFound("curvature"))
+        if (curvature)
         {
-        Info<< nl << "Extracting curvature of surface at the points." << endl;
+            Info<< nl << "Extracting curvature of surface at the points."
+                << endl;
 
             scalarField k = calcCurvature(surf);
 
@@ -1196,13 +1203,13 @@ int main(int argc, char *argv[])
 #endif
 
 
-    if (args.optionFound("featureProximity"))
+        if (featureProximity)
         {
-        Info<< nl << "Extracting proximity of close feature points and edges "
+            Info<< nl << "Extracting proximity of close feature points and "
-            << "to the surface" << endl;
+                << "edges to the surface" << endl;
 
             const scalar searchDistance =
-            args.optionRead<scalar>("featureProximity");
+                readScalar(surfaceDict.lookup("maxFeatureProximity"));
 
             const scalar radiusSqr = sqr(searchDistance);
 
@@ -1269,6 +1276,9 @@ int main(int argc, char *argv[])
             }
         }
 
+        Info<< endl;
+    }
+
     Info<< "End\n" << endl;
 
     return 0;

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict

@@ -0,0 +1,86 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      surfaceFeatureExtractDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+surface1.stl
+{
+    // extractFromFile || extractFromSurface
+    extractionMethod    extractFromFile;
+
+    extractFromFile
+    {
+        // Load from an existing feature edge file
+        featureEdgeFile "constant/triSurface/featureEdges.nas";
+    }
+    
+    trimFeatures
+    {
+        // Remove features with fewer than the specified number of edges
+        minElem         0;
+        
+        // Remove features shorter than the specified cumulative length
+        minLen          0.0;
+    }  
+    
+    subsetFeatures
+    {        
+        // Use a plane to select feature edges
+        // (normal)(basePoint)
+        plane               (1 0 0)(0 0 0);
+        
+        // Select feature edges using a box
+        // (minPt)(maxPt)
+        insideBox           (0 0 0)(1 1 1);
+        outsideBox          (0 0 0)(1 1 1);
+        
+        // Remove any non-manifold (open or > 2 connected faces) edges
+        manifoldEdges       no;
+    }
+    
+    // Output the curvature of the surface
+    curvature               no;
+    
+    // Output the proximity of feature points and edges to each other
+    featureProximity        no;
+    // The maximum search distance to use when looking for other feature
+    // points and edges
+    maxFeatureProximity     1;
+    
+    // Out put the closeness of surface elements to other surface elements.
+    closeness               no;
+
+    // Write options
+    writeVTK                no;
+    writeObj                yes;
+    writeFeatureEdgeMesh    no;
+}
+
+
+surface2.nas
+{
+    extractionMethod    extractFromSurface;
+
+    extractFromSurface
+    {
+        // Mark edges whose adjacent surface normals are at an angle less
+        // than includedAngle as features
+        // - 0  : selects no edges
+        // - 180: selects all edges
+        includedAngle   120;
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C

@@ -2832,7 +2832,7 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
 ) const
 {
     label nearestShapeI = -1;
-    point nearestPoint;
+    point nearestPoint = vector::zero;
 
     if (nodes_.size())
     {
@@ -2847,10 +2847,6 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
             nearestPoint
         );
     }
-    else
-    {
-        nearestPoint = vector::zero;
-    }
 
     return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
 }

src/edgeMesh/edgeFormats/nas/NASedgeFormat.C

@@ -113,6 +113,17 @@ bool Foam::fileFormats::NASedgeFormat::read
             // discard groupID
             dynEdges.append(e);
         }
+        else if (cmd == "PLOTEL")
+        {
+            edge e;
+
+            // label groupId = readLabel(IStringStream(line.substr(16,8))());
+            e[0] = readLabel(IStringStream(line.substr(16,8))());
+            e[1] = readLabel(IStringStream(line.substr(24,8))());
+
+            // discard groupID
+            dynEdges.append(e);
+        }
         else if (cmd == "GRID")
         {
             label index = readLabel(IStringStream(line.substr(8,8))());

src/edgeMesh/edgeMesh.C

@@ -28,6 +28,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "addToMemberFunctionSelectionTable.H"
 #include "ListOps.H"
+#include "EdgeMap.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -327,10 +328,7 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
         }
 
         // Compact using a hashtable and commutative hash of edge.
-        HashTable<label, edge, Hash<edge> > edgeToLabel
+        EdgeMap<label> edgeToLabel(2*edges_.size());
-        (
-            2*edges_.size()
-        );
 
         label newEdgeI = 0;
 
@@ -349,13 +347,7 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
 
         edges_.setSize(newEdgeI);
 
-        for
+        forAllConstIter(EdgeMap<label>, edgeToLabel, iter)
-        (
-            HashTable<label, edge, Hash<edge> >::const_iterator iter =
-                edgeToLabel.begin();
-            iter != edgeToLabel.end();
-            ++iter
-        )
         {
             edges_[iter()] = iter.key();
         }
@@ -363,4 +355,35 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
 }
 
 
+void Foam::edgeMesh::mergeEdges()
+{
+    EdgeMap<label> existingEdges(2*edges_.size());
+
+    label curEdgeI = 0;
+    forAll(edges_, edgeI)
+    {
+        const edge& e = edges_[edgeI];
+
+        if (existingEdges.insert(e, curEdgeI))
+        {
+            curEdgeI++;
+        }
+    }
+
+    if (debug)
+    {
+        Info<< "Merging duplicate edges: "
+            << edges_.size() - existingEdges.size()
+            << " edges will be deleted." << endl;
+    }
+
+    edges_.setSize(existingEdges.size());
+
+    forAllConstIter(EdgeMap<label>, existingEdges, iter)
+    {
+        edges_[iter()] = iter.key();
+    }
+}
+
+
 // ************************************************************************* //

src/edgeMesh/edgeMesh.H

@@ -246,6 +246,9 @@ public:
         //- Merge common points (points within mergeDist)
         void mergePoints(const scalar mergeDist);
 
+        //- Merge similar edges
+        void mergeEdges();
+
 
     // Write
 

src/finiteVolume/Make/files

@@ -124,6 +124,7 @@ $(derivedFvPatchFields)/mappedFixedPushedInternalValue/mappedFixedPushedInternal
 $(derivedFvPatchFields)/mappedFixedValue/mappedFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C
 $(derivedFvPatchFields)/mappedFlowRate/mappedFlowRateFvPatchVectorField.C
+$(derivedFvPatchFields)/fixedMean/fixedMeanFvPatchFields.C
 $(derivedFvPatchFields)/fan/fanFvPatchFields.C
 $(derivedFvPatchFields)/fanPressure/fanPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/buoyantPressure/buoyantPressureFvPatchScalarField.C

src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedMeanFvPatchField.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<Type>(p, iF),
+    meanValue_(pTraits<Type>::zero)
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fixedMeanFvPatchField<Type>& ptf,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchField<Type>(ptf, p, iF, mapper),
+    meanValue_(ptf.meanValue_)
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<Type>(p, iF, dict),
+    meanValue_(pTraits<Type>(dict.lookup("meanValue")))
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fixedMeanFvPatchField<Type>& ptf
+)
+:
+    fixedValueFvPatchField<Type>(ptf),
+    meanValue_(ptf.meanValue_)
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fixedMeanFvPatchField<Type>& ptf,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<Type>(ptf, iF),
+    meanValue_(ptf.meanValue_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void fixedMeanFvPatchField<Type>::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    Field<Type> newValues(this->patchInternalField());
+
+    Type meanValuePsi =
+        gSum(this->patch().magSf()*newValues)
+       /gSum(this->patch().magSf());
+
+    if (mag(meanValue_) > SMALL && mag(meanValuePsi)/mag(meanValue_) > 0.5)
+    {
+        newValues *= mag(meanValue_)/mag(meanValuePsi);
+    }
+    else
+    {
+        newValues += (meanValue_ - meanValuePsi);
+    }
+
+    this->operator==(newValues);
+
+    fixedValueFvPatchField<Type>::updateCoeffs();
+}
+
+
+template<class Type>
+void fixedMeanFvPatchField<Type>::write(Ostream& os) const
+{
+    fvPatchField<Type>::write(os);
+    os.writeKeyword("meanValue") << meanValue_ << token::END_STATEMENT << nl;
+    this->writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fixedMeanFvPatchField
+
+Description
+    Extrapolates field to the patch using the near-cell values and adjusts
+    the distribution to match the specified meanValue.
+
+SourceFiles
+    fixedMeanFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedMeanFvPatchField_H
+#define fixedMeanFvPatchField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class fixedMeanFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class fixedMeanFvPatchField
+:
+    public fixedValueFvPatchField<Type>
+{
+
+protected:
+
+    // Protected data
+
+        //- MeanValue value the field is adjusted to maintain
+        Type meanValue_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fixedMean");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        fixedMeanFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        fixedMeanFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given fixedMeanFvPatchField
+        //  onto a new patch
+        fixedMeanFvPatchField
+        (
+            const fixedMeanFvPatchField<Type>&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        fixedMeanFvPatchField
+        (
+            const fixedMeanFvPatchField<Type>&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<Type> > clone() const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new fixedMeanFvPatchField<Type>(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        fixedMeanFvPatchField
+        (
+            const fixedMeanFvPatchField<Type>&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<Type> > clone
+        (
+            const DimensionedField<Type, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new fixedMeanFvPatchField<Type>(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "fixedMeanFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.C

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedMeanFvPatchFields.H"
+#include "volMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makePatchFields(fixedMean);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedMeanFvPatchFields_H
+#define fixedMeanFvPatchFields_H
+
+#include "fixedMeanFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(fixedMean);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFieldsFwd.H

@@ -23,67 +23,28 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "buildCGALPolyhedron.H"
+#ifndef fixedMeanFvPatchFieldsFwd_H
+#define fixedMeanFvPatchFieldsFwd_H
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+#include "fieldTypes.H"
 
-Foam::buildCGALPolyhedron::buildCGALPolyhedron
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-(
-    const Foam::triSurface& surf
-)
-:
-    CGAL::Modifier_base<HalfedgeDS>(),
-    surf_(surf)
-{}
 
-
+namespace Foam
-
-// * * * * * * * * * * * * * * * * 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
+template<class Type> class fixedMeanFvPatchField;
-    (
-        surf_.points().size(),      // n points
-        surf_.size(),               // n facets
-        2*surf_.edges().size()      // n halfedges
-    );
 
-    forAll(surf_.points(), pI)
+makePatchTypeFieldTypedefs(fixedMean);
-    {
-        const Foam::point& p = surf_.points()[pI];
 
-        B.add_vertex(Point(p.x(), p.y(), p.z()));
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-    }
 
-    forAll(surf_, fI)
+} // End namespace Foam
-    {
-        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();
-}
 
+#endif
 
 // ************************************************************************* //

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C

@@ -33,11 +33,15 @@ License
 #include "OFstream.H"
 #include "IFstream.H"
 #include "unitConversion.H"
+#include "EdgeMap.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::surfaceFeatures, 0);
 
+const Foam::scalar Foam::surfaceFeatures::parallelTolerance =
+    sin(degToRad(1.0));
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -170,7 +174,10 @@ void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
 
 
 //construct feature points where more than 2 feature edges meet
-void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
+void Foam::surfaceFeatures::calcFeatPoints
+(
+    const List<edgeStatus>& edgeStat
+)
 {
     DynamicList<label> featurePoints(surf_.nPoints()/1000);
 
@@ -200,6 +207,60 @@ void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
 }
 
 
+void Foam::surfaceFeatures::classifyFeatureAngles
+(
+    const labelListList& edgeFaces,
+    List<edgeStatus>& edgeStat,
+    const scalar minCos
+) const
+{
+    const vectorField& faceNormals = surf_.faceNormals();
+    const pointField& points = surf_.points();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() != 2)
+        {
+            // Non-manifold. What to do here? Is region edge? external edge?
+            edgeStat[edgeI] = REGION;
+        }
+        else
+        {
+            label face0 = eFaces[0];
+            label face1 = eFaces[1];
+
+            if (surf_[face0].region() != surf_[face1].region())
+            {
+                edgeStat[edgeI] = REGION;
+            }
+            else
+            {
+                if ((faceNormals[face0] & faceNormals[face1]) < minCos)
+                {
+
+                    // Check if convex or concave by looking at angle
+                    // between face centres and normal
+                    vector f0Tof1 =
+                        surf_[face1].centre(points)
+                      - surf_[face0].centre(points);
+
+                    if ((f0Tof1 & faceNormals[face0]) > 0.0)
+                    {
+                        edgeStat[edgeI] = INTERNAL;
+                    }
+                    else
+                    {
+                        edgeStat[edgeI] = EXTERNAL;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 // Returns next feature edge connected to pointI with correct value.
 Foam::label Foam::surfaceFeatures::nextFeatEdge
 (
@@ -436,6 +497,80 @@ Foam::surfaceFeatures::surfaceFeatures
 }
 
 
+Foam::surfaceFeatures::surfaceFeatures
+(
+    const triSurface& surf,
+    const pointField& points,
+    const edgeList& edges,
+    const scalar mergeTol
+)
+:
+    surf_(surf),
+    featurePoints_(0),
+    featureEdges_(0),
+    externalStart_(0),
+    internalStart_(0)
+{
+    // Match edge mesh edges with the triSurface edges
+
+    const labelListList& surfEdgeFaces = surf_.edgeFaces();
+    const edgeList& surfEdges = surf_.edges();
+
+    scalar mergeTolSqr = sqr(mergeTol);
+
+    EdgeMap<label> dynFeatEdges(2*edges.size());
+    DynamicList<labelList> dynFeatureEdgeFaces(edges.size());
+
+    labelList edgeLabel;
+
+    nearestFeatEdge
+    (
+        edges,
+        points,
+        mergeTolSqr,
+        edgeLabel     // label of surface edge or -1
+    );
+
+    label count = 0;
+    forAll(edgeLabel, sEdgeI)
+    {
+        const label sEdge = edgeLabel[sEdgeI];
+
+        if (sEdge == -1)
+        {
+            continue;
+        }
+
+        dynFeatEdges.insert(surfEdges[sEdge], count++);
+        dynFeatureEdgeFaces.append(surfEdgeFaces[sEdge]);
+    }
+
+    // Find whether an edge is external or internal
+    List<edgeStatus> edgeStat(dynFeatEdges.size(), NONE);
+
+    classifyFeatureAngles(dynFeatureEdgeFaces, edgeStat, GREAT);
+
+    // Transfer the edge status to a list encompassing all edges in the surface
+    // so that calcFeatPoints can be used.
+    List<edgeStatus> allEdgeStat(surf_.nEdges(), NONE);
+
+    forAll(allEdgeStat, eI)
+    {
+        EdgeMap<label>::const_iterator iter = dynFeatEdges.find(surfEdges[eI]);
+
+        if (iter != dynFeatEdges.end())
+        {
+            allEdgeStat[eI] = edgeStat[iter()];
+        }
+    }
+
+    edgeStat.clear();
+    dynFeatEdges.clear();
+
+    setFromStatus(allEdgeStat);
+}
+
+
 //- Construct as copy
 Foam::surfaceFeatures::surfaceFeatures(const surfaceFeatures& sf)
 :
@@ -496,54 +631,10 @@ void Foam::surfaceFeatures::findFeatures(const scalar includedAngle)
 {
     scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
 
-    const labelListList& edgeFaces = surf_.edgeFaces();
-    const vectorField& faceNormals = surf_.faceNormals();
-    const pointField& points = surf_.points();
-
     // Per edge whether is feature edge.
     List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
 
-    forAll(edgeFaces, edgeI)
+    classifyFeatureAngles(surf_.edgeFaces(), edgeStat, minCos);
-    {
-        const labelList& eFaces = edgeFaces[edgeI];
-
-        if (eFaces.size() != 2)
-        {
-            // Non-manifold. What to do here? Is region edge? external edge?
-            edgeStat[edgeI] = REGION;
-        }
-        else
-        {
-            label face0 = eFaces[0];
-            label face1 = eFaces[1];
-
-            if (surf_[face0].region() != surf_[face1].region())
-            {
-                edgeStat[edgeI] = REGION;
-            }
-            else
-            {
-                if ((faceNormals[face0] & faceNormals[face1]) < minCos)
-                {
-
-                    // Check if convex or concave by looking at angle
-                    // between face centres and normal
-                    vector f0Tof1 =
-                        surf_[face1].centre(points)
-                        - surf_[face0].centre(points);
-
-                    if ((f0Tof1 & faceNormals[face0]) > 0.0)
-                    {
-                        edgeStat[edgeI] = INTERNAL;
-                    }
-                    else
-                    {
-                        edgeStat[edgeI] = EXTERNAL;
-                    }
-                }
-            }
-        }
-    }
 
     setFromStatus(edgeStat);
 }
@@ -1147,6 +1238,9 @@ void Foam::surfaceFeatures::nearestSurfEdge
 
     const pointField& localPoints = surf_.localPoints();
 
+    treeBoundBox searchDomain(localPoints);
+    searchDomain.inflate(0.1);
+
     indexedOctree<treeDataEdge> ppTree
     (
         treeDataEdge
@@ -1156,7 +1250,7 @@ void Foam::surfaceFeatures::nearestSurfEdge
             localPoints,
             selectedEdges
         ),          // all information needed to do geometric checks
-        treeBoundBox(localPoints),  // overall search domain
+        searchDomain,  // overall search domain
         8,      // maxLevel
         10,     // leafsize
         3.0     // duplicity
@@ -1216,6 +1310,8 @@ void Foam::surfaceFeatures::nearestSurfEdge
     pointOnEdge.setSize(selectedSampleEdges.size());
     pointOnFeature.setSize(selectedSampleEdges.size());
 
+    treeBoundBox searchDomain(surf_.localPoints());
+
     indexedOctree<treeDataEdge> ppTree
     (
         treeDataEdge
@@ -1225,7 +1321,7 @@ void Foam::surfaceFeatures::nearestSurfEdge
             surf_.localPoints(),
             selectedEdges
         ),              // all information needed to do geometric checks
-        treeBoundBox(surf_.localPoints()),  // overall search domain
+        searchDomain,   // overall search domain
         8,              // maxLevel
         10,             // leafsize
         3.0             // duplicity
@@ -1262,6 +1358,67 @@ void Foam::surfaceFeatures::nearestSurfEdge
 }
 
 
+void Foam::surfaceFeatures::nearestFeatEdge
+(
+    const edgeList& edges,
+    const pointField& points,
+    scalar searchSpanSqr,   // Search span
+    labelList& edgeLabel
+) const
+{
+    edgeLabel = labelList(surf_.nEdges(), -1);
+
+    treeBoundBox searchDomain(points);
+    searchDomain.inflate(0.1);
+
+    indexedOctree<treeDataEdge> ppTree
+    (
+        treeDataEdge
+        (
+            false,
+            edges,
+            points,
+            identity(edges.size())
+        ),          // all information needed to do geometric checks
+        searchDomain,  // overall search domain
+        8,      // maxLevel
+        10,     // leafsize
+        3.0     // duplicity
+    );
+
+    const edgeList& surfEdges = surf_.edges();
+    const pointField& surfLocalPoints = surf_.localPoints();
+
+    forAll(surfEdges, edgeI)
+    {
+        const edge& sample = surfEdges[edgeI];
+
+        const point& startPoint = surfLocalPoints[sample.start()];
+        const point& midPoint = sample.centre(surfLocalPoints);
+
+        pointIndexHit infoMid = ppTree.findNearest
+        (
+            midPoint,
+            searchSpanSqr
+        );
+
+        if (infoMid.hit())
+        {
+            const vector surfEdgeDir = midPoint - startPoint;
+
+            const edge& featEdge = edges[infoMid.index()];
+            const vector featEdgeDir = featEdge.vec(points);
+
+            // Check that the edges are nearly parallel
+            if (mag(surfEdgeDir ^ featEdgeDir) < parallelTolerance)
+            {
+                edgeLabel[edgeI] = edgeI;
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 void Foam::surfaceFeatures::operator=(const surfaceFeatures& rhs)

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H

@@ -95,6 +95,11 @@ private:
         {}
     };
 
+    // Static data
+
+        //- Tolerance for determining whether two vectors are parallel
+        static const scalar parallelTolerance;
+
 
     // Private data
 
@@ -130,6 +135,14 @@ private:
         //- Construct feature points where more than 2 feature edges meet
         void calcFeatPoints(const List<edgeStatus>&);
 
+        //- Classify the angles of the feature edges
+        void classifyFeatureAngles
+        (
+            const labelListList& edgeFaces,
+            List<edgeStatus>& edgeStat,
+            const scalar minCos
+        ) const;
+
         //- Choose next unset feature edge.
         label nextFeatEdge
         (
@@ -186,6 +199,15 @@ public:
         //- Construct from file
         surfaceFeatures(const triSurface&, const fileName& fName);
 
+        //- Construct from pointField and edgeList (edgeMesh)
+        surfaceFeatures
+        (
+            const triSurface&,
+            const pointField& points,
+            const edgeList& edges,
+            const scalar mergeTol = 1e-6
+        );
+
         //- Construct as copy
         surfaceFeatures(const surfaceFeatures&);
 
@@ -327,7 +349,6 @@ public:
                 pointField& edgePoint
             ) const;
 
-
             //- Find nearest surface edge (out of selectedEdges) for each
             //  sample edge.
             //  Sets:
@@ -347,6 +368,16 @@ public:
                 pointField& pointOnFeature  // point on sample edge
             ) const;
 
+            //- Find nearest feature edge to each surface edge. Uses the
+            //  mid-point of the surface edges.
+            void nearestFeatEdge
+            (
+                const edgeList& edges,
+                const pointField& points,
+                scalar searchSpanSqr,
+                labelList& edgeLabel
+            ) const;
+
 
         // Write