Added densityWeightedStochastic initial points method which places random points

with a density corresponding to that of the desired mesh.  Added and
restructured queries to make this possible.

Also, for all initialPointsMethods, made the minimumSurfaceDistance into
minimumSurfaceDistanceCoeff so that the guard distance for the surface is a
function of the local cell size.

src/conformalVoronoiMesh/Make/files

@@ -20,6 +20,7 @@ $(cellSiseFunctions)/linearSpatial/linearSpatial.C
 initialPointsMethod/initialPointsMethod/initialPointsMethod.C
 initialPointsMethod/uniformGrid/uniformGrid.C
 initialPointsMethod/pointFile/pointFile.C
+initialPointsMethod/densityWeightedStochastic/densityWeightedStochastic.C
 
 relaxationModel/relaxationModel/relaxationModel.C
 relaxationModel/adaptiveLinear/adaptiveLinear.C

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C

@@ -213,6 +213,40 @@ Foam::scalar Foam::cellSizeControlSurfaces::cellSize
 }
 
 
+Foam::scalarField Foam::cellSizeControlSurfaces::cellSize
+(
+    const pointField& pts,
+    const List<bool>& isSurfacePoint
+) const
+{
+    if (pts.size() != isSurfacePoint.size())
+    {
+        FatalErrorIn
+        (
+            "Foam::cellSizeControlSurfaces::cellSizeControlSurfaces \
+             ( \
+                 const pointField& pts, \
+                 const List<bool>& isSurfacePoint \
+             ) \
+             const"
+        )   << "Size of pointField (" << pts.size()
+            << ") and List<bool> (" << isSurfacePoint.size()
+            << ") do not match." << nl
+            << exit(FatalError);
+    }
+
+    scalarField cellSizes(pts.size());
+
+    forAll(pts, i)
+    {
+        cellSizes[i] = cellSize(pts[i], isSurfacePoint[i]);
+    }
+
+    return cellSizes;
+}
+
+
+
 // ************************************************************************* //
 
 

src/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H

@@ -131,6 +131,13 @@ public:
                 const point& pt,
                 bool isSurfacePoint = false
             ) const;
+
+            //- Return the cell size at the given locations
+            scalarField cellSize
+            (
+                const pointField& pts,
+                const List<bool>& isSurfacePoint
+            ) const;
 };
 
 

src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -2135,6 +2135,8 @@ void Foam::conformalVoronoiMesh::move()
         }
     }
 
+    timeCheck();
+
     Info<< nl << "    Looping over all dual faces" << endl;
 
     vectorField cartesianDirections(3);

src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C

@@ -45,7 +45,8 @@ Foam::conformationSurfaces::conformationSurfaces
     baffleSurfaces_(),
     patchNames_(0),
     patchOffsets_(),
-    bounds_()
+    bounds_(),
+    referenceVolumeTypes_(0)
 {
     const dictionary& surfacesDict
     (
@@ -166,50 +167,11 @@ Foam::conformationSurfaces::conformationSurfaces
 
     bounds_ = searchableSurfacesQueries::bounds(allGeometry_, surfaces_);
 
-    if(cvMesh_.cvMeshControls().objOutput())
-    {
-        writeFeatureObj("cvMesh");
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::conformationSurfaces::~conformationSurfaces()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::Field<bool> Foam::conformationSurfaces::inside
-(
-    const pointField& samplePts
-) const
-{
-    return wellInside(samplePts, 0.0);
-}
-
-
-Foam::Field<bool> Foam::conformationSurfaces::outside
-(
-    const pointField& samplePts
-) const
-{
-    return wellOutside(samplePts, 0.0);
-}
-
-
-Foam::Field<bool> Foam::conformationSurfaces::wellInside
-(
-    const pointField& samplePts,
-    const scalar dist2
-) const
-{
     // Look at all surfaces at determine whether the locationInMesh point is
     // inside or outside each, to establish a signature for the domain to be
     // meshed.
 
-    List<searchableSurface::volumeType> referenceVolumeTypes
+    referenceVolumeTypes_.setSize
     (
         surfaces_.size(),
         searchableSurface::UNKNOWN
@@ -231,10 +193,49 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
 
             surface.getVolumeType(pts, vTypes);
 
-            referenceVolumeTypes[s] = vTypes[0];
+            referenceVolumeTypes_[s] = vTypes[0];
         }
     }
 
+    if(cvMesh_.cvMeshControls().objOutput())
+    {
+        writeFeatureObj("cvMesh");
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::conformationSurfaces::~conformationSurfaces()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::Field<bool> Foam::conformationSurfaces::inside
+(
+    const pointField& samplePts
+) const
+{
+    return wellInside(samplePts, scalarField(samplePts.size(), 0.0));
+}
+
+
+Foam::Field<bool> Foam::conformationSurfaces::outside
+(
+    const pointField& samplePts
+) const
+{
+    return wellOutside(samplePts, scalarField(samplePts.size(), 0.0));
+}
+
+
+Foam::Field<bool> Foam::conformationSurfaces::wellInside
+(
+    const pointField& samplePts,
+    const scalarField& testDistSqr
+) const
+{
     List<List<searchableSurface::volumeType> > surfaceVolumeTests
     (
         surfaces_.size(),
@@ -264,8 +265,6 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
 
     //Check if the points are inside the surface by the given distance squared
 
-    scalarField testDistSqr(insidePoints.size(), dist2);
-
     labelList hitSurfaces;
 
     List<pointIndexHit> hitInfo;
@@ -293,7 +292,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
 
         forAll(surfaces_, s)
         {
-            if (surfaceVolumeTests[s][i] != referenceVolumeTypes[s])
+            if (surfaceVolumeTests[s][i] != referenceVolumeTypes_[s])
             {
                 insidePoints[i] = false;
 
@@ -306,10 +305,20 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
 }
 
 
+bool Foam::conformationSurfaces::wellInside
+(
+    const point& samplePt,
+    scalar testDistSqr
+) const
+{
+    return wellInside(pointField(1, samplePt), scalarField(1, testDistSqr))[0];
+}
+
+
 Foam::Field<bool> Foam::conformationSurfaces::wellOutside
 (
     const pointField& samplePts,
-    const scalar dist2
+    const scalarField& testDistSqr
 ) const
 {
     notImplemented("Field<bool> Foam::conformationSurfaces::wellOutside");
@@ -318,6 +327,16 @@ Foam::Field<bool> Foam::conformationSurfaces::wellOutside
 }
 
 
+bool Foam::conformationSurfaces::wellOutside
+(
+    const point& samplePt,
+    scalar testDistSqr
+) const
+{
+    return wellOutside(pointField(1, samplePt), scalarField(1, testDistSqr))[0];
+}
+
+
 bool Foam::conformationSurfaces::findSurfaceAnyIntersection
 (
     const point& start,

src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H

@@ -94,6 +94,10 @@ class conformationSurfaces
         //- The overall boundBox of all of the surfaces to be conformed to
         boundBox bounds_;
 
+        //- The pattern/signature of volumeTypes representing a point in the
+        // domain to be meshed
+        List<searchableSurface::volumeType> referenceVolumeTypes_;
+
 
     // Private Member Functions
 
@@ -151,19 +155,31 @@ public:
             Field<bool> outside(const pointField& samplePts) const;
 
             //- Check if point is inside surfaces to conform to by at least
-            //  dist2
+            //  testDistSqr
             Field<bool> wellInside
             (
                 const pointField& samplePts,
-                const scalar dist2
+                const scalarField& testDistSqr
+            ) const;
+
+            bool wellInside
+            (
+                const point& samplePt,
+                scalar testDistSqr
             ) const;
 
             //- Check if point is outside surfaces to conform to by at least
-            // dist2
+            //  testDistSqr
             Field<bool> wellOutside
             (
                 const pointField& samplePts,
-                const scalar dist2
+                const scalarField& testDistSqr
+            ) const;
+
+            bool wellOutside
+            (
+                const point& samplePt,
+                scalar testDistSqr
             ) const;
 
             // Finding if the line joining start and end intersects the surface

src/conformalVoronoiMesh/initialPointsMethod/densityWeightedStochastic/densityWeightedStochastic.C

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "densityWeightedStochastic.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(densityWeightedStochastic, 0);
+addToRunTimeSelectionTable
+(
+    initialPointsMethod,
+    densityWeightedStochastic,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+densityWeightedStochastic::densityWeightedStochastic
+(
+    const dictionary& initialPointsDict,
+    const conformalVoronoiMesh& cvMesh
+)
+:
+    initialPointsMethod(typeName, initialPointsDict, cvMesh),
+    totalVolume_(readScalar(detailsDict().lookup("totalVolume"))),
+    maxDensity_
+    (
+        1.0/pow3(readScalar(detailsDict().lookup("minCellSize")))
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+std::vector<Vb::Point> densityWeightedStochastic::initialPoints() const
+{
+    const boundBox& bb = cvMesh_.geometryToConformTo().bounds();
+
+    Random rndGen(5234986);
+
+    std::vector<Vb::Point> initialPoints;
+
+    scalar volumeAdded = 0.0;
+
+    const point& min = bb.min();
+
+    vector span = bb.span();
+
+    while (volumeAdded < totalVolume_)
+    {
+        point p =
+            min
+          + vector
+            (
+                span.x()*rndGen.scalar01(),
+                span.y()*rndGen.scalar01(),
+                span.z()*rndGen.scalar01()
+            );
+
+        scalar localSize = cvMesh_.cellSizeControl().cellSize(p);
+
+        scalar localDensity = 1/pow3(max(localSize, VSMALL));
+
+        // Accept possible placements proportional to the relative local density
+        if (localDensity/maxDensity_ > rndGen.scalar01())
+        {
+            // Determine if the point is "wellInside" the domain
+            if
+            (
+                cvMesh_.geometryToConformTo().wellInside
+                (
+                    p,
+                    minimumSurfaceDistanceCoeffSqr_*sqr(localSize)
+                )
+            )
+            {
+                initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
+
+                volumeAdded += 1/localDensity;
+            }
+        }
+    }
+
+    return initialPoints;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/conformalVoronoiMesh/initialPointsMethod/densityWeightedStochastic/densityWeightedStochastic.H

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::densityWeightedStochastic
+
+Description
+
+SourceFiles
+    densityWeightedStochastic.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef densityWeightedStochastic_H
+#define densityWeightedStochastic_H
+
+#include "Switch.H"
+#include "Random.H"
+#include "initialPointsMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class densityWeightedStochastic Declaration
+\*---------------------------------------------------------------------------*/
+
+class densityWeightedStochastic
+:
+    public initialPointsMethod
+{
+
+private:
+
+    // Private data
+
+        //- The total volume to be filled
+        scalar totalVolume_;
+
+        //- 1/(minimum cell size)^3, gives the maximum permissible point
+        //  density
+        scalar maxDensity_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("densityWeightedStochastic");
+
+    // Constructors
+
+        //- Construct from components
+        densityWeightedStochastic
+        (
+            const dictionary& initialPointsDict,
+            const conformalVoronoiMesh& cvMesh
+        );
+
+
+    //- Destructor
+    virtual ~densityWeightedStochastic()
+    {}
+
+
+    // Member Functions
+
+        //- Return the initial points for the conformalVoronoiMesh
+        virtual std::vector<Vb::Point> initialPoints() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/conformalVoronoiMesh/initialPointsMethod/initialPointsMethod/initialPointsMethod.C

@@ -50,11 +50,14 @@ initialPointsMethod::initialPointsMethod
     dictionary(initialPointsDict),
     cvMesh_(cvMesh),
     detailsDict_(subDict(type + "Details")),
-    minimumSurfaceDistance_
+    minimumSurfaceDistanceCoeffSqr_
     (
-        readScalar
+        sqr
         (
-            initialPointsDict.lookup("minimumSurfaceDistance")
+            readScalar
+            (
+                initialPointsDict.lookup("minimumSurfaceDistanceCoeff")
+            )
         )
     )
 {}

src/conformalVoronoiMesh/initialPointsMethod/initialPointsMethod/initialPointsMethod.H

@@ -67,8 +67,9 @@ protected:
         dictionary detailsDict_;
 
         //- Only allow the placement of initial points that are within the
-        //  surfaces to be meshed by minimumSurfaceDistance
+        //  surfaces to be meshed by minimumSurfaceDistanceCoeff multiplied by
-        scalar minimumSurfaceDistance_;
+        //  the local target cell size.  Store square of value.
+        scalar minimumSurfaceDistanceCoeffSqr_;
 
 private:
 

src/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C

@@ -73,7 +73,12 @@ std::vector<Vb::Point> pointFile::initialPoints() const
     Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside
     (
         points,
-        sqr(minimumSurfaceDistance_)
+        minimumSurfaceDistanceCoeffSqr_
+       *cvMesh_.cellSizeControl().cellSize
+        (
+            points,
+            List<bool>(points.size(), false)
+        )
     );
 
     forAll(insidePoints, i)

src/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C

@@ -82,6 +82,8 @@ std::vector<Vb::Point> uniformGrid::initialPoints() const
 
     std::vector<Vb::Point> initialPoints;
 
+    List<bool> isSurfacePoint(nk, false);
+
     for (int i = 0; i < ni; i++)
     {
         for (int j = 0; j < nj; j++)
@@ -116,7 +118,8 @@ std::vector<Vb::Point> uniformGrid::initialPoints() const
             Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside
             (
                 points,
-                sqr(minimumSurfaceDistance_)
+                minimumSurfaceDistanceCoeffSqr_
+               *cvMesh_.cellSizeControl().cellSize(points, isSurfacePoint)
             );
 
             forAll(insidePoints, i)