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

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C

@@ -122,7 +122,7 @@ Foam::DelaunayMesh<Triangulation>::DelaunayMesh
             if (indices.headerOk())
             {
                 vh->index() = indices[ptI];
-                vertexCount()++;
+                vertexCount_++;
             }
             else
             {

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H

@@ -201,7 +201,6 @@ public:
         inline void resetCellCount();
 
         inline label vertexCount() const;
-        inline label& vertexCount();
 
         inline void resetVertexCount();
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshI.H

@@ -124,12 +124,6 @@ Foam::label Foam::DelaunayMesh<Triangulation>::vertexCount() const
     return vertexCount_;
 }
 
-template<class Triangulation>
-Foam::label& Foam::DelaunayMesh<Triangulation>::vertexCount()
-{
-    return vertexCount_;
-}
-
 
 template<class Triangulation>
 void Foam::DelaunayMesh<Triangulation>::resetVertexCount()

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellAspectRatioControl/cellAspectRatioControl.C

@@ -40,11 +40,15 @@ Foam::cellAspectRatioControl::cellAspectRatioControl
         aspectRatioDict_.lookupOrDefault<vector>
         (
             "aspectRatioDirection",
-            vector(0, 0, 0)
+            vector::zero
         )
     )
 {
-    Info<< nl << "Cell Aspect Ratio Control" << nl
+    // Normalise the direction
+    aspectRatioDirection_ /= mag(aspectRatioDirection_) + SMALL;
+
+    Info<< nl
+        << "Cell Aspect Ratio Control" << nl
         << "    Ratio     : " << aspectRatio_ << nl
         << "    Direction : " << aspectRatioDirection_
         << endl;
@@ -66,20 +70,18 @@ void Foam::cellAspectRatioControl::updateCellSizeAndFaceArea
     scalar& targetCellSize
 ) const
 {
-    const scalar cosAngle = mag
+    const scalar cosAngle =
-    (
+        mag(vectorTools::cosPhi(alignmentDir, aspectRatioDirection_));
-        vectorTools::cosPhi(alignmentDir, aspectRatioDirection_)
-    );
 
     // Change target face area based on aspect ratio
+    targetFaceArea +=
         targetFaceArea
-        += targetFaceArea
        *(aspectRatio_ - 1.0)
        *(1.0 - cosAngle);
 
     // Change target cell size based on aspect ratio
+    targetCellSize +=
         targetCellSize
-        += targetCellSize
        *(aspectRatio_ - 1.0)
        *cosAngle;
 
@@ -95,12 +97,11 @@ void Foam::cellAspectRatioControl::updateDeltaVector
     vector& delta
 ) const
 {
-    const scalar cosAngle = mag
+    const scalar cosAngle =
-    (
+        mag(vectorTools::cosPhi(alignmentDir, aspectRatioDirection_));
-        vectorTools::cosPhi(alignmentDir, aspectRatioDirection_)
-    );
 
-    delta += 0.5
+    delta +=
+        0.5
        *delta
        *cosAngle
        *(targetCellSize/rABMag)

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellAspectRatioControl/cellAspectRatioControl.H

@@ -56,7 +56,7 @@ class cellAspectRatioControl
 
         const scalar aspectRatio_;
 
-        const vector aspectRatioDirection_;
+        vector aspectRatioDirection_;
 
 
     // Private Member Functions
@@ -73,10 +73,7 @@ public:
     // Constructors
 
         //- Construct from dictionary
-        cellAspectRatioControl
+        cellAspectRatioControl(const dictionary& motionDict);
-        (
-            const dictionary& motionDict
-        );
 
 
     //- Destructor

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

@@ -570,159 +570,8 @@ void Foam::conformalVoronoiMesh::buildCellSizeAndAlignmentMesh()
 }
 
 
-void Foam::conformalVoronoiMesh::storeSizesAndAlignments()
-{
-    DynamicList<Point> storePts(number_of_vertices());
-
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        vit++
-    )
-    {
-        if (vit->internalPoint())
-        {
-            storePts.append(vit->point());
-        }
-    }
-
-    storePts.shrink();
-
-    storeSizesAndAlignments(storePts);
-}
-
-
-void Foam::conformalVoronoiMesh::storeSizesAndAlignments
-(
-    const List<Point>& storePts
-)
-{
-//    timeCheck("Start of storeSizesAndAlignments");
-//
-//    Info << nl << "Store size and alignment" << endl;
-//
-//    sizeAndAlignmentLocations_.setSize(storePts.size());
-//
-//    storedSizes_.setSize(sizeAndAlignmentLocations_.size());
-//
-//    storedAlignments_.setSize(sizeAndAlignmentLocations_.size());
-//
-//    label i = 0;
-//
-//    //checkCellSizing();
-//
-//    for
-//    (
-//        List<Point>::const_iterator pit = storePts.begin();
-//        pit != storePts.end();
-//        ++pit
-//    )
-//    {
-//        pointFromPoint pt = topoint(*pit);
-//
-////        storedAlignments_[i] = requiredAlignment(pt);
-////
-////        storedSizes_[i] = cellShapeControls().cellSize(pt);
-//
-//        cellShapeControls().cellSizeAndAlignment
-//        (
-//            pt,
-//            storedSizes_[i],
-//            storedAlignments_[i]
-//        );
-//
-//        i++;
-//    }
-//
-//    timeCheck("Sizes and alignments calculated, build tree");
-//
-//    buildSizeAndAlignmentTree();
-//
-//    timeCheck("Size and alignment tree built");
-}
-
-
-void Foam::conformalVoronoiMesh::updateSizesAndAlignments
-(
-    const List<Point>& storePts
-)
-{
-    // This function is only used in serial, the background redistribution
-    // triggers this when unbalance is detected in parallel.
-
-    if
-    (
-        !Pstream::parRun()
-     && runTime_.run()
-     && runTime_.timeIndex()
-    )
-    {
-        storeSizesAndAlignments(storePts);
-
-        timeCheck("Updated sizes and alignments");
-    }
-}
-
-
-const Foam::indexedOctree<Foam::treeDataPoint>&
-Foam::conformalVoronoiMesh::sizeAndAlignmentTree() const
-{
-    if (sizeAndAlignmentTreePtr_.empty())
-    {
-        buildSizeAndAlignmentTree();
-    }
-
-    return sizeAndAlignmentTreePtr_();
-}
-
-
 void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
 {
-//    Info<< nl << "Looking up target cell alignment and size" << endl;
-//
-//    const indexedOctree<treeDataPoint>& tree = sizeAndAlignmentTree();
-//
-//    for
-//    (
-//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-//        vit != finite_vertices_end();
-//        vit++
-//    )
-//    {
-//        if
-//        (
-//            vit->internalOrBoundaryPoint()
-//         || vit->referredInternalOrBoundaryPoint()
-//        )
-//        {
-//            pointFromPoint pt = topoint(vit->point());
-//
-//            pointIndexHit info = tree.findNearest(pt, sqr(GREAT));
-//
-//            if (info.hit())
-//            {
-//                vit->alignment() = storedAlignments_[info.index()];
-//
-//                vit->targetCellSize() = storedSizes_[info.index()];
-//            }
-//            else
-//            {
-//                WarningIn
-//                (
-//                    "void "
-//                    "Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()"
-//                )
-//                    << "Point " << pt << " did not find a nearest point "
-//                    << " for alignment and size lookup." << endl;
-//
-//                vit->alignment() = cellShapeControls().cellAlignment(pt);
-//
-//                vit->targetCellSize() = cellShapeControls().cellSize(pt);
-//            }
-//        }
-//    }
-
     Info<< nl << "Calculating target cell alignment and size" << endl;
 
     for
@@ -744,44 +593,6 @@ void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
             );
         }
     }
-
-//    OFstream str(runTime_.path()/"alignments_internal.obj");
-//
-//    for
-//    (
-//        Finite_vertices_iterator vit = finite_vertices_begin();
-//        vit != finite_vertices_end();
-//        ++vit
-//    )
-//    {
-//        if (!vit->farPoint())
-//        {
-//            // Write alignments
-//            const tensor& alignment = vit->alignment();
-//            pointFromPoint pt = topoint(vit->point());
-//
-//            if
-//            (
-//                alignment.x() == triad::unset[0]
-//             || alignment.y() == triad::unset[0]
-//             || alignment.z() == triad::unset[0]
-//            )
-//            {
-//                Info<< "Bad alignment = " << vit->info();
-//
-//                vit->alignment() = tensor::I;
-//
-//                Info<< "New alignment = " << vit->info();
-//
-//                continue;
-//            }
-//
-//            meshTools::writeOBJ(str, pt, alignment.x() + pt);
-//            meshTools::writeOBJ(str, pt, alignment.y() + pt);
-//            meshTools::writeOBJ(str, pt, alignment.z() + pt);
-//        }
-//    }
-
 }
 
 
@@ -1039,10 +850,6 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
     featurePointLocations_(),
     edgeLocationTreePtr_(),
     surfacePtLocationTreePtr_(),
-    sizeAndAlignmentLocations_(),
-    storedSizes_(),
-    storedAlignments_(),
-    sizeAndAlignmentTreePtr_(),
     surfaceConformationVertices_(),
     initialPointsMethod_
     (
@@ -1133,10 +940,6 @@ void Foam::conformalVoronoiMesh::initialiseForMotion()
     // (potentially) redistributed.
     storeSurfaceConformation();
 
-    // Use storeSizesAndAlignments with no feed points because all background
-    // points may have been distributed.
-    storeSizesAndAlignments();
-
     // Report any Delaunay vertices that do not think that they are in the
     // domain the processor they are on.
     // reportProcessorOccupancy();
@@ -1836,8 +1639,6 @@ void Foam::conformalVoronoiMesh::move()
         writeMesh(time().timeName());
     }
 
-    updateSizesAndAlignments(pointsToInsert);
-
     Info<< nl
         << "Total displacement = " << totalDisp << nl
         << "Total distance = " << totalDist << nl

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

@@ -177,19 +177,6 @@ private:
 
         mutable DynamicList<Foam::point> existingSurfacePtLocations_;
 
-        //- Store locations where the cell size and alignments will be
-        //  pre-calculated and looked up
-        pointField sizeAndAlignmentLocations_;
-
-        //- Stored cell size at sizeAndAlignmentLocations_
-        scalarField storedSizes_;
-
-        //- Stored alignments at sizeAndAlignmentLocations_
-        tensorField storedAlignments_;
-
-        //- Search tree for size and alignment lookup points
-        mutable autoPtr<indexedOctree<treeDataPoint> > sizeAndAlignmentTreePtr_;
-
         //- Store the surface and feature edge conformation locations to be
         //  reinserted
         List<Vb> surfaceConformationVertices_;
@@ -271,40 +258,6 @@ private:
         //- Return the local maximum surface protrusion distance
         inline scalar maxSurfaceProtrusion(const Foam::point& pt) const;
 
-        //- Insert Point and return its auto-generated index
-        inline bool insertPoint
-        (
-            const Point& P,
-            const indexedVertexEnum::vertexType type
-        );
-
-        //- Insert Foam::point with specified index and type
-        inline bool insertPoint
-        (
-            const Foam::point& p,
-            const indexedVertexEnum::vertexType type
-        );
-
-        //- Insert Point with specified index, type and original processor
-        inline bool insertReferredPoint
-        (
-            const Point& P,
-            const label index,
-            const indexedVertexEnum::vertexType type,
-            const label processor
-        );
-
-        inline bool insertReferredPoint(const Vb& P);
-
-        //- Insert Foam::point with specified index, type and original processor
-        inline bool insertReferredPoint
-        (
-            const Foam::point& p,
-            const label index,
-            const indexedVertexEnum::vertexType type,
-            const label processor
-        );
-
         //- Insert Delaunay vertices using the CGAL range insertion method,
         //  optionally check processor occupancy and distribute to other
         //  processors
@@ -539,22 +492,6 @@ private:
 
         void buildCellSizeAndAlignmentMesh();
 
-        //- Store data for sizeAndAlignmentLocations_, storedSizes_ and
-        //  storedAlignments_ and initialise the sizeAndAlignmentTreePtr_,
-        //  determining the appropriate sizeAndAlignmentLocations_
-        //  automatically
-        void storeSizesAndAlignments();
-
-        //- Store data for sizeAndAlignmentLocations_, storedSizes_ and
-        //  storedAlignments_ and initialise the sizeAndAlignmentTreePtr_
-        void storeSizesAndAlignments(const List<Point>& storePts);
-
-        //- Restore the sizes and alignments if required
-        void updateSizesAndAlignments(const List<Point>& storePts);
-
-        //- Demand driven construction of octree for and alignment points
-        const indexedOctree<treeDataPoint>& sizeAndAlignmentTree() const;
-
         //- Set the size and alignment data for each vertex
         void setVertexSizeAndAlignment();
 
@@ -750,9 +687,6 @@ private:
             const DynamicList<Foam::point>& existingSurfacePtLocations
         ) const;
 
-        //- Build or rebuild the sizeAndAlignmentTree
-        void buildSizeAndAlignmentTree() const;
-
         //- Process the surface conformation locations to decide which surface
         //  and edge conformation locations to add
         void addSurfaceAndEdgeHits
@@ -1063,10 +997,6 @@ public:
 //            const List<scalar>& radiusSqrs
 //        ) const;
 
-        typedef K::Vector_3 CGALVector;
-
-        inline CGALVector toCGALVector(const Foam::vector& v) const;
-
 
         // Access
 

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

@@ -211,11 +211,6 @@ void Foam::conformalVoronoiMesh::checkDuals()
 
     List<Point> duals(number_of_finite_cells());
 
-    typedef CGAL::Exact_predicates_exact_constructions_kernel       EK2;
-    typedef CGAL::Regular_triangulation_euclidean_traits_3<EK2>     EK;
-    typedef CGAL::Cartesian_converter<baseK::Kernel, EK2>  To_exact;
-    typedef CGAL::Cartesian_converter<EK2, baseK::Kernel>  Back_from_exact;
-
 //    PackedBoolList bPoints(number_of_finite_cells());
 
 //    indexDualVertices(duals, bPoints);
@@ -1787,7 +1782,6 @@ void Foam::conformalVoronoiMesh::indexDualVertices
     OBJstream snapping1("snapToSurface1.obj");
     OBJstream snapping2("snapToSurface2.obj");
     OFstream tetToSnapTo("tetsToSnapTo.obj");
-    label offset = 0;
 
     for
     (

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

@@ -76,10 +76,6 @@ void Foam::conformalVoronoiMesh::conformToSurface()
             {
                 sync(decomposition_().procBounds());
             }
-
-            // Use storeSizesAndAlignments with no feed points because all
-            // background points may have been distributed.
-            storeSizesAndAlignments();
         }
 
         // Do not store the surface conformation until after it has been
@@ -2055,38 +2051,6 @@ void Foam::conformalVoronoiMesh::buildSurfacePtLocationTree
 }
 
 
-void Foam::conformalVoronoiMesh::buildSizeAndAlignmentTree() const
-{
-    if (sizeAndAlignmentLocations_.empty())
-    {
-        FatalErrorIn("buildSizeAndAlignmentTree()")
-            << "sizeAndAlignmentLocations empty, must be populated before "
-            << "sizeAndAlignmentTree can be built."
-            << exit(FatalError);
-    }
-
-    treeBoundBox overallBb
-    (
-        geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
-    );
-
-    overallBb.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-    overallBb.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-    sizeAndAlignmentTreePtr_.reset
-    (
-        new indexedOctree<treeDataPoint>
-        (
-            treeDataPoint(sizeAndAlignmentLocations_),
-            overallBb,  // overall search domain
-            10,         // max levels
-            20.0,       // maximum ratio of cubes v.s. cells
-            100.0       // max. duplicity; n/a since no bounding boxes.
-        )
-    );
-}
-
-
 void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
 (
     const Foam::point& vit,

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

@@ -236,94 +236,6 @@ inline Foam::scalar Foam::conformalVoronoiMesh::maxSurfaceProtrusion
 }
 
 
-inline bool Foam::conformalVoronoiMesh::insertPoint
-(
-    const Foam::point& p,
-    const indexedVertexEnum::vertexType type
-)
-{
-    return insertPoint(toPoint<Point>(p), type);
-}
-
-
-inline bool Foam::conformalVoronoiMesh::insertPoint
-(
-    const Point& P,
-    const indexedVertexEnum::vertexType type
-)
-{
-    uint nVert = number_of_vertices();
-
-    Vertex_handle vh = insert(P);
-
-    bool pointInserted = true;
-
-    if (nVert == number_of_vertices())
-    {
-        Pout<< "Failed to insert point : " << topoint(P)
-            << " of type " << type << endl;
-        pointInserted = false;
-    }
-    else
-    {
-        vh->index() = getNewVertexIndex();
-        vh->type() = type;
-    }
-
-    return pointInserted;
-}
-
-
-inline bool Foam::conformalVoronoiMesh::insertReferredPoint(const Vb& P)
-{
-    return insertReferredPoint(P.point(), P.index(), P.type(), P.procIndex());
-}
-
-
-inline bool Foam::conformalVoronoiMesh::insertReferredPoint
-(
-    const Foam::point& p,
-    const label index,
-    const indexedVertexEnum::vertexType type,
-    const label processor
-)
-{
-    return insertReferredPoint(toPoint<Point>(p), index, type, processor);
-}
-
-
-inline bool Foam::conformalVoronoiMesh::insertReferredPoint
-(
-    const Point& P,
-    const label index,
-    const indexedVertexEnum::vertexType type,
-    const label processor
-)
-{
-    uint nVert = number_of_vertices();
-
-    Vertex_handle vh = insert(P);
-
-    bool pointInserted = true;
-
-    if (nVert == number_of_vertices())
-    {
-        Pout<< "Failed to insert point " << topoint(P)
-            << " type: " << type << " index: " << index
-            << " proc: " << processor << endl;
-        pointInserted = false;
-    }
-    else
-    {
-        vh->index() = index;
-        vh->type() = type;
-        vh->procIndex() = processor;
-    }
-
-    return pointInserted;
-}
-
-
 inline void Foam::conformalVoronoiMesh::createPointPair
 (
     const scalar ppDist,
@@ -663,13 +575,6 @@ inline bool Foam::conformalVoronoiMesh::isProcBoundaryEdge
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::conformalVoronoiMesh::CGALVector
-Foam::conformalVoronoiMesh::toCGALVector(const Foam::vector& v) const
-{
-    return CGALVector(v.x(), v.y(), v.z());
-}
-
-
 inline const Foam::Time& Foam::conformalVoronoiMesh::time() const
 {
     return runTime_;

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H

@@ -108,17 +108,15 @@ int CGAL::indexedCell<Gt, Cb>::cellIndex() const
 
 
 #ifdef CGAL_INEXACT
+
     template<class Gt, class Cb>
     const Foam::point& CGAL::indexedCell<Gt, Cb>::dual()
     {
-    //    if (Foam::foamyHexMeshChecks::coplanarTet(*this, 1e-20) == 0)
-    //    {
-            // Do exact calc
-    //    }
-
         return reinterpret_cast<const Foam::point&>(this->circumcenter());
     }
+
 #else
+
     template<class Gt, class Cb>
     const Foam::point CGAL::indexedCell<Gt, Cb>::dual()
     {
@@ -131,6 +129,7 @@ int CGAL::indexedCell<Gt, Cb>::cellIndex() const
             CGAL::to_double(P.z())
         );
     }
+
 #endif
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H

@@ -25,8 +25,7 @@ Class
     pointFeatureEdgesTypes
 
 Description
-    struct for holding information on the types of feature edges attached to
+    Holds information on the types of feature edges attached to feature points.
-    feature points
 
 \*---------------------------------------------------------------------------*/
 

etc/caseDicts/foamyHexMeshDict

@@ -5,12 +5,19 @@
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      foamyHexMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #inputMode merge;
 
 surfaceConformation
 {
-    locationInMesh                      (0 0 0);
     pointPairDistanceCoeff              0.1;
     mixedFeaturePointPPDistanceCoeff    5.0;
     featurePointExclusionDistanceCoeff  0.65;
@@ -115,20 +122,18 @@ polyMeshFiltering
 }
 
 
+backgroundMeshDecomposition
+{
+    minLevels           1;
+    sampleResolution    4;
+    spanScale           20;
+    maxCellWeightCoeff  20;
+}
+
+
 meshQualityControls
 {
-    maxNonOrtho         65;
+    #include "meshQualityDict"
-    maxBoundarySkewness 50;
-    maxInternalSkewness 10;
-    maxConcave          80;
-    minVol              -1E30;
-    minTetQuality       1e-30;
-    minArea             -1;
-    minTwist            0.02;
-    minDeterminant      0.001;
-    minFaceWeight       0.02;
-    minVolRatio         0.01;
-    minTriangleTwist    -1;
 }
 
 

src/OpenFOAM/meshes/meshShapes/face/face.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -829,4 +829,26 @@ Foam::label Foam::face::trianglesQuads
 }
 
 
+Foam::label Foam::longestEdge(const face& f, const pointField& pts)
+{
+    const edgeList& eds = f.edges();
+
+    label longestEdgeI = -1;
+    scalar longestEdgeLength = -SMALL;
+
+    forAll(eds, edI)
+    {
+        scalar edgeLength = eds[edI].mag(pts);
+
+        if (edgeLength > longestEdgeLength)
+        {
+            longestEdgeI = edI;
+            longestEdgeLength = edgeLength;
+        }
+    }
+
+    return longestEdgeI;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/meshShapes/face/face.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -408,6 +408,12 @@ public:
 };
 
 
+// Global functions
+
+//- Find the longest edge on a face. Face point labels index into pts.
+label longestEdge(const face& f, const pointField& pts);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C

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

src/OpenFOAM/meshes/polyMesh/zones/pointZone/pointZone.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -154,8 +154,15 @@ bool Foam::pointZone::checkParallelSync(const bool report) const
 
     forAll(maxZone, pointI)
     {
-        // Check point in zone on both sides
+        // Check point in same (or no) zone on all processors
-        if (maxZone[pointI] != minZone[pointI])
+        if
+        (
+            (
+                maxZone[pointI] != -1
+             || minZone[pointI] != labelMax
+            )
+         && (maxZone[pointI] != minZone[pointI])
+        )
         {
             if (report && !error)
             {
@@ -167,7 +174,8 @@ bool Foam::pointZone::checkParallelSync(const bool report) const
                     << (minZone[pointI] == labelMax ? -1 : minZone[pointI])
                     << " on some processors and in zone "
                     << maxZone[pointI]
-                    << " on some other processors."
+                    << " on some other processors." << nl
+                    << "(suppressing further warnings)"
                     << endl;
             }
             error = true;

src/OpenFOAM/primitives/Pair/Pair.H

@@ -108,12 +108,6 @@ public:
             return this->operator[](1);
         }
 
-        //- Return reverse pair
-        inline Pair<Type> reversePair() const
-        {
-            return Pair<Type>(second(), first());
-        }
-
         //- Return other
         inline const Type& other(const Type& a) const
         {
@@ -147,11 +141,11 @@ public:
         //  - -1: same pair, but reversed order
         static inline int compare(const Pair<Type>& a, const Pair<Type>& b)
         {
-            if (a[0] == b[0] && a[1] == b[1])
+            if (a == b)
             {
                 return 1;
             }
-            else if (a[0] == b[1] && a[1] == b[0])
+            else if (a == reverse(b))
             {
                 return -1;
             }
@@ -160,21 +154,32 @@ public:
                 return 0;
             }
         }
+};
 
 
-    // Friend Operators
+template<class Type>
+Pair<Type> reverse(const Pair<Type>& p)
+{
+    return Pair<Type>(p.second(), p.first());
+}
 
-        friend bool operator==(const Pair<Type>& a, const Pair<Type>& b)
+
+template<class Type>
+bool operator==(const Pair<Type>& a, const Pair<Type>& b)
 {
     return (a.first() == b.first() && a.second() == b.second());
 }
 
-        friend bool operator!=(const Pair<Type>& a, const Pair<Type>& b)
+
+template<class Type>
+bool operator!=(const Pair<Type>& a, const Pair<Type>& b)
 {
     return !(a == b);
 }
 
-        friend bool operator<(const Pair<Type>& a, const Pair<Type>& b)
+
+template<class Type>
+bool operator<(const Pair<Type>& a, const Pair<Type>& b)
 {
     return
     (
@@ -186,7 +191,27 @@ public:
         )
     );
 }
-};
+
+
+template<class Type>
+bool operator<=(const Pair<Type>& a, const Pair<Type>& b)
+{
+    return !(b < a);
+}
+
+
+template<class Type>
+bool operator>(const Pair<Type>& a, const Pair<Type>& b)
+{
+    return (b < a);
+}
+
+
+template<class Type>
+bool operator>=(const Pair<Type>& a, const Pair<Type>& b)
+{
+    return !(a < b);
+}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/primitives/Tuple2/Tuple2.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,20 +47,6 @@ namespace Foam
 template<class Type1, class Type2>
 class Tuple2;
 
-template<class Type1, class Type2>
-inline bool operator==
-(
-    const Tuple2<Type1, Type2>&,
-    const Tuple2<Type1, Type2>&
-);
-
-template<class Type1, class Type2>
-inline bool operator!=
-(
-    const Tuple2<Type1, Type2>&,
-    const Tuple2<Type1, Type2>&
-);
-
 template<class Type1, class Type2>
 inline Istream& operator>>(Istream&, Tuple2<Type1, Type2>&);
 
@@ -129,27 +115,6 @@ public:
             return s_;
         }
 
-        //- Return reverse pair
-        inline Tuple2<Type2, Type1> reverseTuple2() const
-        {
-            return Tuple2<Type2, Type1>(second(), first());
-        }
-
-
-    // Friend Operators
-
-        friend bool operator== <Type1, Type2>
-        (
-            const Tuple2<Type1, Type2>& a,
-            const Tuple2<Type1, Type2>& b
-        );
-
-        friend bool operator!= <Type1, Type2>
-        (
-            const Tuple2<Type1, Type2>& a,
-            const Tuple2<Type1, Type2>& b
-        );
-
 
     // IOstream operators
 
@@ -169,6 +134,14 @@ public:
 };
 
 
+//- Return reverse of a tuple2
+template<class Type1, class Type2>
+inline Tuple2<Type2, Type1> reverse(const Tuple2<Type1, Type2>& t)
+{
+    return Tuple2<Type2, Type1>(t.second(), t.first());
+}
+
+
 template<class Type1, class Type2>
 inline bool operator==
 (

src/OpenFOAM/primitives/triad/triad.C

@@ -282,6 +282,11 @@ void Foam::triad::align(const vector& v)
 
 Foam::triad Foam::triad::sortxyz() const
 {
+    if (!this->set())
+    {
+        return *this;
+    }
+
     triad t;
 
     if

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C

@@ -42,32 +42,6 @@ defineTypeNameAndDebug(edgeCollapser, 0);
 }
 
 
-Foam::label Foam::edgeCollapser::longestEdge
-(
-    const face& f,
-    const pointField& pts
-)
-{
-    const edgeList& eds = f.edges();
-
-    label longestEdgeI = -1;
-    scalar longestEdgeLength = -SMALL;
-
-    forAll(eds, edI)
-    {
-        scalar edgeLength = eds[edI].mag(pts);
-
-        if (edgeLength > longestEdgeLength)
-        {
-            longestEdgeI = edI;
-            longestEdgeLength = edgeLength;
-        }
-    }
-
-    return longestEdgeI;
-}
-
-
 Foam::HashSet<Foam::label> Foam::edgeCollapser::checkBadFaces
 (
     const polyMesh& mesh,

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H

@@ -262,9 +262,6 @@ public:
 
         // Check
 
-            //- Find the longest edge in a face
-            static label longestEdge(const face& f, const pointField& pts);
-
             //- Calls motionSmoother::checkMesh and returns a set of bad faces
             static HashSet<label> checkBadFaces
             (

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -450,7 +450,7 @@ Foam::Map<Foam::labelPair>  Foam::meshRefinement::getZoneBafflePatches
                     labelPair patches = zPatches;
                     if (fZone.flipMap()[i])
                     {
-                       patches = patches.reversePair();
+                       patches = reverse(patches);
                     }
 
                     if (!bafflePatch.insert(faceI, patches))

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -605,8 +605,10 @@ void Foam::refinementSurfaces::setMinLevelFields
         // searchableBox (size=6)
         if (geom.regions().size() > 1 && geom.globalSize() > 10)
         {
-            // Representative local coordinates
+            // Representative local coordinates and bounding sphere
-            const pointField ctrs(geom.coordinates());
+            pointField ctrs;
+            scalarField radiusSqr;
+            geom.boundingSpheres(ctrs, radiusSqr);
 
             labelList minLevelField(ctrs.size(), -1);
             {
@@ -614,12 +616,7 @@ void Foam::refinementSurfaces::setMinLevelFields
                 // distributed surface where local indices differ from global
                 // ones (needed for getRegion call)
                 List<pointIndexHit> info;
-                geom.findNearest
+                geom.findNearest(ctrs, radiusSqr, info);
-                (
-                    ctrs,
-                    scalarField(ctrs.size(), sqr(GREAT)),
-                    info
-                );
 
                 // Get per element the region
                 labelList region;
@@ -628,7 +625,7 @@ void Foam::refinementSurfaces::setMinLevelFields
                 // From the region get the surface-wise refinement level
                 forAll(minLevelField, i)
                 {
-                    if (info[i].hit())
+                    if (info[i].hit()) //Note: should not be necessary
                     {
                         minLevelField[i] = minLevel(surfI, region[i]);
                     }

src/meshTools/searchableSurface/searchableBox.C

@@ -249,6 +249,41 @@ Foam::tmp<Foam::pointField> Foam::searchableBox::coordinates() const
 }
 
 
+void Foam::searchableBox::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    centres.setSize(size());
+    radiusSqr.setSize(size());
+    radiusSqr = 0.0;
+
+    const pointField pts(treeBoundBox::points());
+    const faceList& fcs = treeBoundBox::faces;
+
+    forAll(fcs, i)
+    {
+        const face& f = fcs[i];
+
+        centres[i] = f.centre(pts);
+        forAll(f, fp)
+        {
+            const point& pt = pts[f[fp]];
+
+            radiusSqr[i] = Foam::max
+            (
+                radiusSqr[i],
+                Foam::magSqr(pt-centres[i])
+            );
+        }
+    }
+
+    // Add a bit to make sure all points are tested inside
+    radiusSqr += Foam::sqr(SMALL);
+}
+
+
 Foam::tmp<Foam::pointField> Foam::searchableBox::points() const
 {
     return treeBoundBox::points();

src/meshTools/searchableSurface/searchableBox.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -129,6 +129,14 @@ public:
         //  Usually the element centres (should be of length size()).
         virtual tmp<pointField> coordinates() const;
 
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
+
         //- Get the points that define the surface.
         virtual tmp<pointField> points() const;
 

src/meshTools/searchableSurface/searchableCylinder.C

@@ -47,6 +47,23 @@ Foam::tmp<Foam::pointField> Foam::searchableCylinder::coordinates() const
 }
 
 
+void Foam::searchableCylinder::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    centres.setSize(1);
+    centres[0] = 0.5*(point1_ + point2_);
+
+    radiusSqr.setSize(1);
+    radiusSqr[0] = Foam::magSqr(point1_-centres[0]) + Foam::sqr(radius_);
+
+    // Add a bit to make sure all points are tested inside
+    radiusSqr += Foam::sqr(SMALL);
+}
+
+
 Foam::tmp<Foam::pointField> Foam::searchableCylinder::points() const
 {
     tmp<pointField> tPts(new pointField(2));

src/meshTools/searchableSurface/searchableCylinder.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -154,6 +154,14 @@ public:
         //  Usually the element centres (should be of length size()).
         virtual tmp<pointField> coordinates() const;
 
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
+
         //- Get the points that define the surface.
         virtual tmp<pointField> points() const;
 

src/meshTools/searchableSurface/searchablePlane.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -134,6 +134,20 @@ const Foam::wordList& Foam::searchablePlane::regions() const
 }
 
 
+void Foam::searchablePlane::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    centres.setSize(1);
+    centres[0] = refPoint();
+
+    radiusSqr.setSize(1);
+    radiusSqr[0] = Foam::sqr(GREAT);
+}
+
+
 void Foam::searchablePlane::findNearest
 (
     const pointField& samples,

src/meshTools/searchableSurface/searchablePlane.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -130,6 +130,15 @@ public:
             return tCtrs;
         }
 
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        //  Note: radius limited to sqr(GREAT)
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
+
         //- Get the points that define the surface.
         virtual tmp<pointField> points() const
         {

src/meshTools/searchableSurface/searchablePlate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -274,6 +274,71 @@ const Foam::wordList& Foam::searchablePlate::regions() const
 }
 
 
+Foam::tmp<Foam::pointField> Foam::searchablePlate::coordinates() const
+{
+    return tmp<pointField>(new pointField(1, origin_ + 0.5*span_));
+}
+
+
+void Foam::searchablePlate::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    centres.setSize(1);
+    centres[0] = origin_ + 0.5*span_;
+
+    radiusSqr.setSize(1);
+    radiusSqr[0] = Foam::magSqr(0.5*span_);
+
+    // Add a bit to make sure all points are tested inside
+    radiusSqr += Foam::sqr(SMALL);
+}
+
+
+Foam::tmp<Foam::pointField> Foam::searchablePlate::points() const
+{
+    tmp<pointField> tPts(new pointField(4));
+    pointField& pts = tPts();
+
+    pts[0] = origin_;
+    pts[2] = origin_ + span_;
+
+    if (span_.x() < span_.y() && span_.x() < span_.z())
+    {
+        pts[1] = origin_ + point(0, span_.y(), 0);
+        pts[3] = origin_ + point(0, 0, span_.z());
+    }
+    else if (span_.y() < span_.z())
+    {
+        pts[1] = origin_ + point(span_.x(), 0, 0);
+        pts[3] = origin_ + point(0, 0, span_.z());
+    }
+    else
+    {
+        pts[1] = origin_ + point(span_.x(), 0, 0);
+        pts[3] = origin_ + point(0, span_.y(), 0);
+    }
+
+    return tPts;
+}
+
+
+bool Foam::searchablePlate::overlaps(const boundBox& bb) const
+{
+    return
+    (
+        (origin_.x() + span_.x()) >= bb.min().x()
+     && origin_.x() <= bb.max().x()
+     && (origin_.y() + span_.y()) >= bb.min().y()
+     && origin_.y() <= bb.max().y()
+     && (origin_.z() + span_.z()) >= bb.min().z()
+     && origin_.z() <= bb.max().z()
+    );
+}
+
+
 void Foam::searchablePlate::findNearest
 (
     const pointField& samples,

src/meshTools/searchableSurface/searchablePlate.H

@@ -145,53 +145,21 @@ public:
 
         //- Get representative set of element coordinates
         //  Usually the element centres (should be of length size()).
-        virtual tmp<pointField> coordinates() const
+        virtual tmp<pointField> coordinates() const;
-        {
+
-            tmp<pointField> tCtrs(new pointField(1, origin_ + 0.5*span_));
+        //- Get bounding spheres (centre and radius squared), one per element.
-            return tCtrs;
+        //  Any point on element is guaranteed to be inside.
-        }
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
 
         //- Get the points that define the surface.
-        virtual tmp<pointField> points() const
+        virtual tmp<pointField> points() const;
-        {
-            tmp<pointField> tPts(new pointField(4));
-            pointField& pts = tPts();
-
-            pts[0] = origin_;
-            pts[2] = origin_ + span_;
-
-            if (span_.x() < span_.y() && span_.x() < span_.z())
-            {
-                pts[1] = origin_ + point(0, span_.y(), 0);
-                pts[3] = origin_ + point(0, 0, span_.z());
-            }
-            else if (span_.y() < span_.z())
-            {
-                pts[1] = origin_ + point(span_.x(), 0, 0);
-                pts[3] = origin_ + point(0, 0, span_.z());
-            }
-            else
-            {
-                pts[1] = origin_ + point(span_.x(), 0, 0);
-                pts[3] = origin_ + point(0, span_.y(), 0);
-            }
-
-            return tPts;
-        }
 
         //- Does any part of the surface overlap the supplied bound box?
-        virtual bool overlaps(const boundBox& bb) const
+        virtual bool overlaps(const boundBox& bb) const;
-        {
-            return
-            (
-                (origin_.x() + span_.x()) >= bb.min().x()
-             && origin_.x() <= bb.max().x()
-             && (origin_.y() + span_.y()) >= bb.min().y()
-             && origin_.y() <= bb.max().y()
-             && (origin_.z() + span_.z()) >= bb.min().z()
-             && origin_.z() <= bb.max().z()
-            );
-        }
 
 
         // Multiple point queries.

src/meshTools/searchableSurface/searchableSphere.C

@@ -185,6 +185,23 @@ const Foam::wordList& Foam::searchableSphere::regions() const
 
 
 
+void Foam::searchableSphere::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    centres.setSize(1);
+    centres[0] = centre_;
+
+    radiusSqr.setSize(1);
+    radiusSqr[0] = Foam::sqr(radius_);
+
+    // Add a bit to make sure all points are tested inside
+    radiusSqr += Foam::sqr(SMALL);
+}
+
+
 void Foam::searchableSphere::findNearest
 (
     const pointField& samples,

src/meshTools/searchableSurface/searchableSphere.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,7 +60,7 @@ private:
         //- Centre point
         const point centre_;
 
-        //- Radius squared
+        //- Radius
         const scalar radius_;
 
         //- Names of regions
@@ -139,6 +139,14 @@ public:
             return tCtrs;
         }
 
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
+
         //- Get the points that define the surface.
         virtual tmp<pointField> points() const
         {

src/meshTools/searchableSurface/searchableSurface.H

@@ -190,6 +190,14 @@ public:
         //  Usually the element centres (should be of length size()).
         virtual tmp<pointField> coordinates() const = 0;
 
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const = 0;
+
         //- Get the points that define the surface.
         virtual tmp<pointField> points() const = 0;
 

src/meshTools/searchableSurface/searchableSurfaceCollection.C

@@ -354,6 +354,42 @@ Foam::searchableSurfaceCollection::coordinates() const
 }
 
 
+void Foam::searchableSurfaceCollection::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    centres.setSize(size());
+    radiusSqr.setSize(centres.size());
+
+    // Append individual coordinates
+    label coordI = 0;
+
+    forAll(subGeom_, surfI)
+    {
+        scalar maxScale = cmptMax(scale_[surfI]);
+
+        pointField subCentres;
+        scalarField subRadiusSqr;
+        subGeom_[surfI].boundingSpheres(subCentres, subRadiusSqr);
+
+        forAll(subCentres, i)
+        {
+            centres[coordI++] = transform_[surfI].globalPosition
+            (
+                cmptMultiply
+                (
+                    subCentres[i],
+                    scale_[surfI]
+                )
+            );
+            radiusSqr[coordI++] = maxScale*subRadiusSqr[i];
+        }
+    }
+}
+
+
 Foam::tmp<Foam::pointField>
 Foam::searchableSurfaceCollection::points() const
 {

src/meshTools/searchableSurface/searchableSurfaceCollection.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -176,6 +176,14 @@ public:
         //  Usually the element centres (should be of length size()).
         virtual tmp<pointField> coordinates() const;
 
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
+
         //- Get the points that define the surface.
         virtual tmp<pointField> points() const;
 

src/meshTools/searchableSurface/searchableSurfaceWithGaps.H

@@ -156,6 +156,17 @@ public:
             return surface().coordinates();
         }
 
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const
+        {
+            surface().boundingSpheres(centres, radiusSqr);
+        }
+
         //- Get the points that define the surface.
         virtual tmp<pointField> points() const
         {

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -145,10 +145,12 @@ bool Foam::triSurfaceMesh::addFaceToEdge
 
 bool Foam::triSurfaceMesh::isSurfaceClosed() const
 {
+    const pointField& pts = triSurface::points();
+
     // Construct pointFaces. Let's hope surface has compact point
     // numbering ...
     labelListList pointFaces;
-    invertManyToMany(points()().size(), *this, pointFaces);
+    invertManyToMany(pts.size(), *this, pointFaces);
 
     // Loop over all faces surrounding point. Count edges emanating from point.
     // Every edge should be used by two faces exactly.
@@ -241,7 +243,9 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
     minQuality_(-1),
     surfaceClosed_(-1)
 {
-    bounds() = boundBox(points());
+    const pointField& pts = triSurface::points();
+
+    bounds() = boundBox(pts);
 }
 
 
@@ -287,7 +291,9 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
     minQuality_(-1),
     surfaceClosed_(-1)
 {
-    bounds() = boundBox(points());
+    const pointField& pts = triSurface::points();
+
+    bounds() = boundBox(pts);
 }
 
 
@@ -347,7 +353,9 @@ Foam::triSurfaceMesh::triSurfaceMesh
         triSurface::scalePoints(scaleFactor);
     }
 
-    bounds() = boundBox(points());
+    const pointField& pts = triSurface::points();
+
+    bounds() = boundBox(pts);
 
     // Have optional minimum quality for normal calculation
     if (dict.readIfPresent("minQuality", minQuality_) && minQuality_ > 0)
@@ -393,6 +401,34 @@ Foam::tmp<Foam::pointField> Foam::triSurfaceMesh::coordinates() const
 }
 
 
+void Foam::triSurfaceMesh::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    centres = coordinates();
+    radiusSqr.setSize(size());
+    radiusSqr = 0.0;
+
+    const pointField& pts = triSurface::points();
+
+    forAll(*this, faceI)
+    {
+        const labelledTri& f = triSurface::operator[](faceI);
+        const point& fc = centres[faceI];
+        forAll(f, fp)
+        {
+            const point& pt = pts[f[fp]];
+            radiusSqr[faceI] = max(radiusSqr[faceI], Foam::magSqr(fc-pt));
+        }
+    }
+
+    // Add a bit to make sure all points are tested inside
+    radiusSqr += Foam::sqr(SMALL);
+}
+
+
 Foam::tmp<Foam::pointField> Foam::triSurfaceMesh::points() const
 {
     return triSurface::points();
@@ -606,6 +642,7 @@ void Foam::triSurfaceMesh::getNormal
 ) const
 {
     const triSurface& s = static_cast<const triSurface&>(*this);
+    const pointField& pts = s.points();
 
     normal.setSize(info.size());
 
@@ -621,9 +658,9 @@ void Foam::triSurfaceMesh::getNormal
             if (info[i].hit())
             {
                 label faceI = info[i].index();
-                normal[i] = s[faceI].normal(points());
+                normal[i] = s[faceI].normal(pts);
 
-                scalar qual = s[faceI].tri(points()).quality();
+                scalar qual = s[faceI].tri(pts).quality();
 
                 if (qual < minQuality_)
                 {
@@ -633,11 +670,11 @@ void Foam::triSurfaceMesh::getNormal
                     forAll(fFaces, j)
                     {
                         label nbrI = fFaces[j];
-                        scalar nbrQual = s[nbrI].tri(points()).quality();
+                        scalar nbrQual = s[nbrI].tri(pts).quality();
                         if (nbrQual > qual)
                         {
                             qual = nbrQual;
-                            normal[i] = s[nbrI].normal(points());
+                            normal[i] = s[nbrI].normal(pts);
                         }
                     }
                 }
@@ -662,7 +699,7 @@ void Foam::triSurfaceMesh::getNormal
                 //normal[i] = faceNormals()[faceI];
 
                 //- Uncached
-                normal[i] = s[faceI].normal(points());
+                normal[i] = s[faceI].normal(pts);
                 normal[i] /= mag(normal[i]) + VSMALL;
             }
             else

src/meshTools/searchableSurface/triSurfaceMesh.H

@@ -187,6 +187,14 @@ public:
             //  Usually the element centres (should be of length size()).
             virtual tmp<pointField> coordinates() const;
 
+            //- Get bounding spheres (centre and radius squared). Any point
+            //  on surface is guaranteed to be inside.
+            virtual void boundingSpheres
+            (
+                pointField& centres,
+                scalarField& radiusSqr
+            ) const;
+
             //- Get the points that define the surface.
             virtual tmp<pointField> points() const;