ENH: particleDistribution function object - updated output

src/functionObjects/field/particleDistribution/particleDistribution.C

@@ -46,6 +46,29 @@ namespace functionObjects
 }
 }
 
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::particleDistribution::writeFileHeader
+(
+    Ostream& os
+) const
+{
+    writeHeader(os, "Particle distribution");
+    writeHeaderValue(os, "Cloud", cloudName_);
+    forAll(nameVsBinWidth_, i)
+    {
+        writeHeaderValue
+        (
+            os,
+            "Field:" + nameVsBinWidth_[i].first(),
+            nameVsBinWidth_[i].second()
+        );
+    }
+    os << endl;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::functionObjects::particleDistribution::particleDistribution
@@ -58,12 +81,12 @@ Foam::functionObjects::particleDistribution::particleDistribution
     fvMeshFunctionObject(name, runTime, dict),
     writeFile(runTime, name, typeName, dict),
     cloudName_("unknown-cloudName"),
-    fieldNames_(),
+    nameVsBinWidth_(),
-    tagFieldName_("unknown-tagFieldName"),
+    tagFieldName_("none"),
-    distributionBinWidth_(0),
     rndGen_(1234, -1)
 {
     read(dict);
+    writeFileHeader(file());
 }
 
 
@@ -80,11 +103,10 @@ bool Foam::functionObjects::particleDistribution::read(const dictionary& dict)
     if (fvMeshFunctionObject::read(dict) && writeFile::read(dict))
     {
         dict.lookup("cloud") >> cloudName_;
-        dict.lookup("fields") >> fieldNames_;
+        dict.lookup("nameVsBinWidth") >> nameVsBinWidth_;
-        dict.lookup("tagField") >> tagFieldName_;
+        dict.readIfPresent("tagField", tagFieldName_);
-        dict.lookup("distributionBinWidth") >> distributionBinWidth_;
 
-        Info<< type() << " " << name() << " output:"
+        Info<< type() << " " << name() << " output:" << nl
             << "    Processing cloud : " << cloudName_ << nl
             << endl;
 
@@ -107,9 +129,12 @@ bool Foam::functionObjects::particleDistribution::write()
 
     if (!mesh_.foundObject<cloud>(cloudName_))
     {
+        wordList cloudNames(mesh_.names<cloud>());
+
         WarningInFunction
             << "Unable to find cloud " << cloudName_
-            << " in the mesh database" << endl;
+            << " in the mesh database.  Available clouds include:"
+            << cloudNames << endl;
 
         return false;
     }
@@ -152,25 +177,31 @@ bool Foam::functionObjects::particleDistribution::write()
         }
     }
 
+    file() << "# Time: " << mesh_.time().timeName() << nl;
+
     bool ok = false;
-    forAll(fieldNames_, i)
+    forAll(nameVsBinWidth_, i)
     {
-        const word fName = fieldNames_[i];
+        ok = false;
-        ok = ok || processField<scalar>(cloudObr, fName, tagAddr);
+        ok = ok || processField<scalar>(cloudObr, i, tagAddr);
-        ok = ok || processField<vector>(cloudObr, fName, tagAddr);
+        ok = ok || processField<vector>(cloudObr, i, tagAddr);
-        ok = ok || processField<tensor>(cloudObr, fName, tagAddr);
+        ok = ok || processField<tensor>(cloudObr, i, tagAddr);
-        ok = ok || processField<sphericalTensor>(cloudObr, fName, tagAddr);
+        ok = ok || processField<sphericalTensor>(cloudObr, i, tagAddr);
-        ok = ok || processField<symmTensor>(cloudObr, fName, tagAddr);
+        ok = ok || processField<symmTensor>(cloudObr, i, tagAddr);
-        ok = ok || processField<tensor>(cloudObr, fName, tagAddr);
+        ok = ok || processField<tensor>(cloudObr, i, tagAddr);
 
         if (log && !ok)
         {
             WarningInFunction
-                << "Unable to find field " << fName << " in the "
+                << "Unable to find field " << nameVsBinWidth_[i].first()
-                << cloudName_ << " cloud database" << endl;
+                << " in the " << cloudName_ << " cloud database" << endl;
         }
     }
 
+    if (ok)
+    {
+        file() << nl;
+    }
 
     return true;
 }
@@ -180,9 +211,15 @@ void Foam::functionObjects::particleDistribution::generateDistribution
 (
     const word& fieldName,
     const scalarField& field,
+    const scalar binWidth,
     const label tag
 )
 {
+    if (field.empty())
+    {
+        return;
+    }
+
     Ostream& os = file();
 
     if (tag != -1)
@@ -193,13 +230,11 @@ void Foam::functionObjects::particleDistribution::generateDistribution
     distributionModels::general distribution
     (
         field,
-        distributionBinWidth_,
+        binWidth,
         rndGen_
     );
 
-    distribution.writeData(os);
+    os  << fieldName << distribution.writeDict(mesh_.time().timeName());
-
-    os  << endl;
 }
 
 

src/functionObjects/field/particleDistribution/particleDistribution.H

@@ -28,7 +28,7 @@ Group
     grpFieldFunctionObjects
 
 Description
-    Generates a particle distrbution for lagrangian data.
+    Generates a particle distrbution for lagrangian data at a given time.
 
 Usage
     \verbatim
@@ -38,8 +38,11 @@ Usage
         libs        ("libfieldFunctionObjects.so");
         ...
         cloud       "myCloud";
-        field       "d"
+        nameVsBinWidth
-        distributionBinWidth 0.1;
+        (
+            (d 0.1)
+            (U 10)
+        );
     }
     \endverbatim
 
@@ -48,8 +51,8 @@ Usage
         Property    | Description              |  Required  | Default value
         type        | Type name: particleDistribution | yes |
         cloud       | Name of cloud to process | Yes        |
-        field       | Name of cloud field to process | Yes  |
+        nameVsBinWidth | List of cloud field vs bin width | Yes |
-        distributionBinWidth | Width of distribution bins | yes |
+        tagField    | Name of cloud field to use to group particles | no | none
     \endtable
 
 See also
@@ -69,6 +72,7 @@ SourceFiles
 #include "writeFile.H"
 #include "scalarField.H"
 #include "cachedRandom.H"
+#include "Tuple2.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -93,26 +97,24 @@ protected:
         //- Cloud name
         word cloudName_;
 
-        //- Field name
+        //- List of field name vs. bin width
-        wordList fieldNames_;
+        List<Tuple2<word, scalar>> nameVsBinWidth_;
 
         //- Tag field name - used to filter the particles into groups
         word tagFieldName_;
 
-        //- Distribution bin width
-        scalar distributionBinWidth_;
-
         //- Random number generator - used by distribution models
         cachedRandom rndGen_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
 
         //- Generate the distribution
         void generateDistribution
         (
             const word& fieldName,
             const scalarField& field,
+            const scalar binWidth,
             const label tag = -1
         );
 
@@ -121,7 +123,7 @@ protected:
         bool processField
         (
             const objectRegistry& obr,
-            const word& fieldName,
+            const label fieldi,
             const List<DynamicList<label>>& addr
         );
 
@@ -131,6 +133,9 @@ protected:
         //- Disallow default bitwise assignment
         void operator=(const particleDistribution&) = delete;
 
+        //- Output file header information
+        virtual void writeFileHeader(Ostream& os) const;
+
 
 public:
 

src/functionObjects/field/particleDistribution/particleDistributionTemplates.C

@@ -29,10 +29,13 @@ template<class Type>
 bool Foam::functionObjects::particleDistribution::processField
 (
     const objectRegistry& obr,
-    const word& fieldName,
+    const label fieldi,
     const List<DynamicList<label>>& addr
 )
 {
+    const word& fieldName = nameVsBinWidth_[fieldi].first();
+    const scalar binWidth = nameVsBinWidth_[fieldi].second();
+
     if (obr.foundObject<IOField<Type>>(fieldName))
     {
         const IOField<Type>& field =
@@ -45,7 +48,13 @@ bool Foam::functionObjects::particleDistribution::processField
                 const Field<Type> subField(field, addr[i]);
                 for (direction d = 0; d < pTraits<Type>::nComponents; ++d)
                 {
-                    generateDistribution(fieldName, subField.component(d), i);
+                    generateDistribution
+                    (
+                        fieldName,
+                        subField.component(d),
+                        binWidth,
+                        i
+                    );
                 }
             }
         }
@@ -53,7 +62,8 @@ bool Foam::functionObjects::particleDistribution::processField
         {
             for (direction d = 0; d < pTraits<Type>::nComponents; ++d)
             {
-                generateDistribution(fieldName, field.component(d));
+                const word fName = fieldName + pTraits<Type>::componentNames[d];
+                generateDistribution(fName, field.component(d), binWidth);
             }
         }
 

src/lagrangian/distributionModels/general/general.C

@@ -42,6 +42,8 @@ namespace distributionModels
 
 void Foam::distributionModels::general::initialise()
 {
+    static scalar eps = ROOTVSMALL;
+
     const label nEntries = xy_.size();
 
     integral_.setSize(nEntries);
@@ -50,7 +52,7 @@ void Foam::distributionModels::general::initialise()
     integral_[0] = 0.0;
     for (label i = 1; i < nEntries; i++)
     {
-        scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0]);
+        scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0] + eps);
         scalar d = xy_[i-1][1] - k*xy_[i-1][0];
         scalar y1 = xy_[i][0]*(0.5*k*xy_[i][0] + d);
         scalar y0 = xy_[i-1][0]*(0.5*k*xy_[i-1][0] + d);
@@ -61,12 +63,12 @@ void Foam::distributionModels::general::initialise()
 
     scalar sumArea = integral_.last();
 
-    meanValue_ = sumArea/(maxValue() - minValue());
+    meanValue_ = sumArea/(maxValue() - minValue() + eps);
 
     for (label i=0; i < nEntries; i++)
     {
-        xy_[i][1] /= sumArea;
+        xy_[i][1] /= sumArea + eps;
-        integral_[i] /= sumArea;
+        integral_[i] /= sumArea + eps;
     }
 }
 
@@ -245,7 +247,17 @@ Foam::dictionary Foam::distributionModels::general::writeDict
 {
     // dictionary dict = distributionModel::writeDict(dictName);
     dictionary dict(dictName);
-    dict.add("distribution", xy_);
+    List<scalar> data(xy_.size());
+    forAll(data, i)
+    {
+        data[i] = xy_[i][0];
+    }
+    dict.add("x", data);
+    forAll(data, i)
+    {
+        data[i] = xy_[i][1];
+    }
+    dict.add("y", data);
 
     return dict;
 }
@@ -254,7 +266,16 @@ Foam::dictionary Foam::distributionModels::general::writeDict
 void Foam::distributionModels::general::readDict(const dictionary& dict)
 {
     // distributionModel::readDict(dict);
-    dict.lookup("distribution") >> xy_;
+    List<scalar> x(dict.lookup("x"));
+    List<scalar> y(dict.lookup("y"));
+
+    xy_.setSize(x.size());
+    forAll(xy_, i)
+    {
+        xy_[i][0] = x[i];
+        xy_[i][1] = y[i];
+    }
+
     initialise();
 }