From 58dae2de43ab706ee8725208414d2c4af7614a58 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Sat, 24 Nov 2018 13:14:57 +0100
Subject: [PATCH] ENH: improvements for dataCloud function object (issue #1044)

- now supports a parcel selection mechanism like vtkCloud,
  giving the ability to select a subset of parcels.
  For example, a given stride, or removal of parcels with a small
  diameter.

  Eg,
      dataCloud output Time: 3.2
      Applying parcel filtering to 994 parcels
      - add stride 4
      - subtract field U : (less 0.2)
      After filtering using 214/994 parcels

- add output precision control for dataCloud
---
 src/functionObjects/lagrangian/Make/files     |   2 +-
 .../lagrangian/cloudInfo/cloudInfo.C          |  51 +++--
 .../lagrangian/cloudInfo/cloudInfo.H          |  14 +-
 .../parcelSelectionDetail.C                   |   0
 .../parcelSelectionDetail.H                   |   0
 .../lagrangian/dataCloud/dataCloud.C          |  28 ++-
 .../lagrangian/dataCloud/dataCloud.H          |  74 ++++++-
 .../lagrangian/dataCloud/dataCloudTemplates.C | 187 +++++++++++++-----
 .../simplifiedSiwek/system/dataCloud          |   1 +
 .../filter/system/controlDict                 |   1 +
 .../filter/system/dataCloud                   |  42 ++++
 .../reactingParcelFoam/filter/system/vtkCloud |   1 +
 12 files changed, 302 insertions(+), 99 deletions(-)
 rename src/functionObjects/lagrangian/{vtkCloud => common}/parcelSelectionDetail.C (100%)
 rename src/functionObjects/lagrangian/{vtkCloud => common}/parcelSelectionDetail.H (100%)
 create mode 100644 tutorials/lagrangian/reactingParcelFoam/filter/system/dataCloud

diff --git a/src/functionObjects/lagrangian/Make/files b/src/functionObjects/lagrangian/Make/files
index 8a65315b9a..e349d4f545 100644
--- a/src/functionObjects/lagrangian/Make/files
+++ b/src/functionObjects/lagrangian/Make/files
@@ -1,9 +1,9 @@
+common/parcelSelectionDetail.C
 dataCloud/dataCloud.C
 cloudInfo/cloudInfo.C
 icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.C
 dsmcFields/dsmcFields.C
 
 vtkCloud/vtkCloud.C
-vtkCloud/parcelSelectionDetail.C
 
 LIB = $(FOAM_LIBBIN)/liblagrangianFunctionObjects
diff --git a/src/functionObjects/lagrangian/cloudInfo/cloudInfo.C b/src/functionObjects/lagrangian/cloudInfo/cloudInfo.C
index d1d379b817..3bbe4a132f 100644
--- a/src/functionObjects/lagrangian/cloudInfo/cloudInfo.C
+++ b/src/functionObjects/lagrangian/cloudInfo/cloudInfo.C
@@ -78,12 +78,6 @@ Foam::functionObjects::cloudInfo::cloudInfo
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::functionObjects::cloudInfo::~cloudInfo()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::functionObjects::cloudInfo::read(const dictionary& dict)
@@ -96,10 +90,10 @@ bool Foam::functionObjects::cloudInfo::read(const dictionary& dict)
         if (writeToFile() && names().size())
         {
             Info<< "applying to clouds:" << nl;
-            forAll(names(), i)
+            forAll(names(), cloudi)
             {
-                Info<< "    " << names()[i] << nl;
-                writeFileHeader(files(i));
+                Info<< "    " << names()[cloudi] << nl;
+                writeFileHeader(files(cloudi));
             }
             Info<< endl;
         }
@@ -107,8 +101,6 @@ bool Foam::functionObjects::cloudInfo::read(const dictionary& dict)
         {
             Info<< "no clouds to be processed" << nl << endl;
         }
-
-        return true;
     }
 
     return true;
@@ -123,24 +115,26 @@ bool Foam::functionObjects::cloudInfo::execute()
 
 bool Foam::functionObjects::cloudInfo::write()
 {
-    forAll(names(), i)
+    forAll(names(), cloudi)
     {
-        const word& cloudName = names()[i];
+        const word& cloudName = names()[cloudi];
 
         const kinematicCloud& cloud =
             obr_.lookupObject<kinematicCloud>(cloudName);
 
-        label nParcels = returnReduce(cloud.nParcels(), sumOp<label>());
-        scalar massInSystem =
+        const label nTotParcels =
+            returnReduce(cloud.nParcels(), sumOp<label>());
+
+        const scalar totMass =
             returnReduce(cloud.massInSystem(), sumOp<scalar>());
 
-        scalar Dmax = cloud.Dmax();
-        scalar D10 = cloud.Dij(1, 0);
-        scalar D32 = cloud.Dij(3, 2);
+        const scalar Dmax = cloud.Dmax();
+        const scalar D10 = cloud.Dij(1, 0);
+        const scalar D32 = cloud.Dij(3, 2);
 
         Log << type() << " " << name() <<  " write:" << nl
-            << "    number of parcels : " << nParcels << nl
-            << "    mass in system    : " << massInSystem << nl
+            << "    number of parcels : " << nTotParcels << nl
+            << "    mass in system    : " << totMass << nl
             << "    maximum diameter  : " << Dmax << nl
             << "    D10 diameter      : " << D10 << nl
             << "    D32 diameter      : " << D32 << nl
@@ -148,14 +142,15 @@ bool Foam::functionObjects::cloudInfo::write()
 
         if (writeToFile())
         {
-            writeTime(files(i));
-            files(i)
-                << token::TAB
-                << nParcels << token::TAB
-                << massInSystem << token::TAB
-                << Dmax << token::TAB
-                << D10 << token::TAB
-                << D32 << token::TAB
+            auto& os = files(cloudi);
+
+            writeTime(os);
+            os
+                << token::TAB << nTotParcels
+                << token::TAB << totMass
+                << token::TAB << Dmax
+                << token::TAB << D10
+                << token::TAB << D32
                 << endl;
         }
     }
diff --git a/src/functionObjects/lagrangian/cloudInfo/cloudInfo.H b/src/functionObjects/lagrangian/cloudInfo/cloudInfo.H
index d4dad2889c..e01b1d4ffb 100644
--- a/src/functionObjects/lagrangian/cloudInfo/cloudInfo.H
+++ b/src/functionObjects/lagrangian/cloudInfo/cloudInfo.H
@@ -52,9 +52,9 @@ Usage
 
     Where the entries comprise:
     \table
-        Property     | Description             | Required    | Default value
-        type         | type name: cloudInfo    | yes         |
-        clouds       | list of clouds names to process |yes  |
+        Property     | Description                          | Required | Default
+        type         | type name: cloudInfo                 | yes |
+        clouds       | list of clouds names to process      | yes |
     \endtable
 
     The output data of each cloud is written to a file named \<cloudName\>.dat
@@ -108,10 +108,6 @@ protected:
         virtual void writeFileHeader(Ostream& os) const;
 
 
-private:
-
-    // Private member functions
-
         //- No copy construct
         cloudInfo(const cloudInfo&) = delete;
 
@@ -137,13 +133,13 @@ public:
 
 
     //- Destructor
-    virtual ~cloudInfo();
+    virtual ~cloudInfo() = default;
 
 
     // Member Functions
 
         //- Read the controls
-        virtual bool read(const dictionary&);
+        virtual bool read(const dictionary& dict);
 
         //- Execute, currently does nothing
         virtual bool execute();
diff --git a/src/functionObjects/lagrangian/vtkCloud/parcelSelectionDetail.C b/src/functionObjects/lagrangian/common/parcelSelectionDetail.C
similarity index 100%
rename from src/functionObjects/lagrangian/vtkCloud/parcelSelectionDetail.C
rename to src/functionObjects/lagrangian/common/parcelSelectionDetail.C
diff --git a/src/functionObjects/lagrangian/vtkCloud/parcelSelectionDetail.H b/src/functionObjects/lagrangian/common/parcelSelectionDetail.H
similarity index 100%
rename from src/functionObjects/lagrangian/vtkCloud/parcelSelectionDetail.H
rename to src/functionObjects/lagrangian/common/parcelSelectionDetail.H
diff --git a/src/functionObjects/lagrangian/dataCloud/dataCloud.C b/src/functionObjects/lagrangian/dataCloud/dataCloud.C
index 53267ab025..bf303c92a3 100644
--- a/src/functionObjects/lagrangian/dataCloud/dataCloud.C
+++ b/src/functionObjects/lagrangian/dataCloud/dataCloud.C
@@ -79,16 +79,29 @@ bool Foam::functionObjects::dataCloud::writeCloud
         return false;
     }
 
+    applyFilter_ = calculateFilter(obrTmp, log);
+    reduce(applyFilter_, orOp<bool>());
+
+
+    // Number of parcels (locally)
+    label nParcels = (applyFilter_ ? parcelAddr_.count() : pointsPtr->size());
+
     // Total number of parcels on all processes
-    label nTotParcels = pointsPtr->size();
-    reduce(nTotParcels, sumOp<label>());
+    const label nTotParcels = returnReduce(nParcels, sumOp<label>());
+
+    if (applyFilter_ && log)
+    {
+        // Report filtered/unfiltered count
+        Log << "After filtering using " << nTotParcels << '/'
+            << (returnReduce(pointsPtr->size(), sumOp<label>()))
+            << " parcels" << nl;
+    }
 
     if (!nTotParcels)
     {
         return false;
     }
 
-
     if (Pstream::master())
     {
         mkDir(outputName.path());
@@ -113,6 +126,9 @@ Foam::functionObjects::dataCloud::dataCloud
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
+    printf_(),
+    precision_(IOstream::defaultPrecision()),
+    applyFilter_(false),
     selectClouds_(),
     fieldName_(),
     directory_()
@@ -139,6 +155,9 @@ bool Foam::functionObjects::dataCloud::read(const dictionary& dict)
         printf_ = "%0" + std::to_string(padWidth) + "d";
     }
 
+    precision_ =
+        dict.lookupOrDefault("precision", IOstream::defaultPrecision());
+
 
     selectClouds_.clear();
     dict.readIfPresent("clouds", selectClouds_);
@@ -152,6 +171,9 @@ bool Foam::functionObjects::dataCloud::read(const dictionary& dict)
 
     dict.readEntry("field", fieldName_);
 
+    // Actions to define selection
+    parcelSelect_ = dict.subOrEmptyDict("selection");
+
     // Output directory
 
     directory_.clear();
diff --git a/src/functionObjects/lagrangian/dataCloud/dataCloud.H b/src/functionObjects/lagrangian/dataCloud/dataCloud.H
index a1ef28b587..002dc54c26 100644
--- a/src/functionObjects/lagrangian/dataCloud/dataCloud.H
+++ b/src/functionObjects/lagrangian/dataCloud/dataCloud.H
@@ -43,19 +43,34 @@ Description
     }
     \endverbatim
 
-Usage
+    \heading Basic Usage
     \table
         Property     | Description                      | Required | Default
         type         | Type name: dataCloud             | yes   |
-        cloud        |                                  | no    | defaultCloud
-        clouds       | wordRe list of clouds            | no    |
+        clouds       | List of clouds (name or regex)   | no  |
+        cloud        | Cloud name                       | no  | defaultCloud
         field        | Name of the field                | yes   |
-        directory    | The output directory name     | no | postProcessing/NAME
-        width        | Padding width for file name      | no    | 8
+        selection    | Parcel selection control         | no  | empty-dict
     \endtable
 
+    \heading Output Options
+    \table
+        Property     | Description                      | Required | Default
+        precision    | The write precision              | no  | same as IOstream
+        directory    | The output directory name    | no | postProcessing/NAME
+        width        | Padding width for file name      | no  | 8
+        writeControl | Output control                   | recommended | timeStep
+    \endtable
+
+Note
+    See Foam::functionObjects::vtkCloud and Foam::Detail::parcelSelection
+    for more details about the parcel selection mechanism.
+
 See also
+    Foam::Detail::parcelSelection
     Foam::functionObjects::vtkCloud
+    Foam::functionObjects::fvMeshFunctionObject
+    Foam::functionObjects::timeControl
 
 SourceFiles
     dataCloud.C
@@ -67,6 +82,7 @@ SourceFiles
 #define functionObjects_dataCloud_H
 
 #include "fvMeshFunctionObject.H"
+#include "parcelSelectionDetail.H"
 #include "vectorField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -82,13 +98,20 @@ namespace functionObjects
 
 class dataCloud
 :
-    public fvMeshFunctionObject
+    public fvMeshFunctionObject,
+    public Foam::Detail::parcelSelection
 {
     // Private data
 
         //- The printf format for zero-padding names
         string printf_;
 
+        //- The output precision
+        unsigned precision_;
+
+        //- Apply output filter (for the current cloud)
+        bool applyFilter_;
+
         //- Requested names of clouds to process
         wordRes selectClouds_;
 
@@ -101,12 +124,47 @@ class dataCloud
 
     // Private Member Functions
 
+        //- Output (point,value) combination on a single line
         template<class Type>
-        static void writeField
+        static void writePointValue
+        (
+            Ostream& os,
+            const vector& pt,
+            const Type& val
+        );
+
+        template<class Type>
+        static void writeList
         (
             Ostream& os,
             const vectorField& points,
-            const Field<Type>& field
+            const List<Type>& field
+        );
+
+        template<class Type>
+        static void writeListParallel
+        (
+            Ostream& os,
+            const vectorField& points,
+            const List<Type>& field
+        );
+
+        template<class Type>
+        static void writeList
+        (
+            Ostream& os,
+            const vectorField& points,
+            const List<Type>& field,
+            const bitSet& selected
+        );
+
+        template<class Type>
+        static void writeListParallel
+        (
+            Ostream& os,
+            const vectorField& points,
+            const List<Type>& field,
+            const bitSet& selected
         );
 
         //- Write to disk
diff --git a/src/functionObjects/lagrangian/dataCloud/dataCloudTemplates.C b/src/functionObjects/lagrangian/dataCloud/dataCloudTemplates.C
index c51033e9e5..2741203320 100644
--- a/src/functionObjects/lagrangian/dataCloud/dataCloudTemplates.C
+++ b/src/functionObjects/lagrangian/dataCloud/dataCloudTemplates.C
@@ -25,33 +25,141 @@ License
 
 #include "IOField.H"
 #include "OFstream.H"
+#include "ListOps.H"
 #include "pointField.H"
 #include "vectorField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class Type>
-void Foam::functionObjects::dataCloud::writeField
+void Foam::functionObjects::dataCloud::writePointValue
+(
+    Ostream& os,
+    const vector& pt,
+    const Type& val
+)
+{
+    os << pt.x() << ' ' << pt.y() << ' ' << pt.z();
+
+    for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
+    {
+        os << ' ' << component(val, cmpt);
+    }
+    os << nl;
+}
+
+
+template<class Type>
+void Foam::functionObjects::dataCloud::writeList
 (
     Ostream& os,
     const vectorField& points,
-    const Field<Type>& field
+    const List<Type>& field
 )
 {
     const label len = field.size();
 
     for (label pointi=0; pointi<len; ++pointi)
     {
-        const point& pt = points[pointi];
-        const Type& val = field[pointi];
+        writePointValue(os, points[pointi], field[pointi]);
+    }
+}
 
-        os << pt.x() << ' ' << pt.y() << ' ' << pt.z();
 
-        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
+template<class Type>
+void Foam::functionObjects::dataCloud::writeListParallel
+(
+    Ostream& os,
+    const vectorField& points,
+    const List<Type>& field
+)
+{
+    if (Pstream::master())
+    {
+        writeList(os, points, field);
+
+        vectorField recvPoints;
+        Field<Type> recvField;
+
+        // Receive and write
+        for (int slave=1; slave<Pstream::nProcs(); ++slave)
         {
-            os << ' ' << component(val, cmpt);
+            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
+
+            fromSlave >> recvPoints >> recvField;
+
+            writeList(os, recvPoints, recvField);
         }
-        os << nl;
+    }
+    else
+    {
+        // Send to master
+        OPstream toMaster
+        (
+            Pstream::commsTypes::blocking,
+            Pstream::masterNo()
+        );
+
+        toMaster
+            << points << field;
+    }
+}
+
+
+template<class Type>
+void Foam::functionObjects::dataCloud::writeList
+(
+    Ostream& os,
+    const vectorField& points,
+    const List<Type>& field,
+    const bitSet& selected
+)
+{
+    for (const label pointi : selected)
+    {
+        writePointValue(os, points[pointi], field[pointi]);
+    }
+}
+
+
+template<class Type>
+void Foam::functionObjects::dataCloud::writeListParallel
+(
+    Ostream& os,
+    const vectorField& points,
+    const List<Type>& field,
+    const bitSet& selected
+)
+{
+    if (Pstream::master())
+    {
+        writeList(os, points, field, selected);
+
+        vectorField recvPoints;
+        Field<Type> recvField;
+
+        // Receive and write
+        for (int slave=1; slave<Pstream::nProcs(); ++slave)
+        {
+            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
+
+            fromSlave >> recvPoints >> recvField;
+
+            writeList(os, recvPoints, recvField);
+        }
+    }
+    else
+    {
+        // Send to master
+        OPstream toMaster
+        (
+            Pstream::commsTypes::blocking,
+            Pstream::masterNo()
+        );
+
+        toMaster
+            << subset(selected, points)
+            << subset(selected, field);
     }
 }
 
@@ -63,16 +171,6 @@ bool Foam::functionObjects::dataCloud::writeField
     const objectRegistry& obrTmp
 ) const
 {
-    // Fields are not always on all processors (eg, multi-component parcels).
-    // Thus need to resolve between all processors.
-
-    const auto* fldPtr = obrTmp.findObject<IOField<Type>>(fieldName_);
-
-    if (!returnReduce((fldPtr != nullptr), orOp<bool>()))
-    {
-        return false;
-    }
-
     const auto* pointsPtr = obrTmp.findObject<vectorField>("position");
 
     if (!pointsPtr)
@@ -81,46 +179,35 @@ bool Foam::functionObjects::dataCloud::writeField
         return false;
     }
 
+    // Fields are not always on all processors (eg, multi-component parcels).
+    // Thus need to resolve between all processors.
+
+    const List<Type>* fldPtr = obrTmp.findObject<IOField<Type>>(fieldName_);
+    const List<Type>& values = (fldPtr ? *fldPtr : List<Type>());
+
+    if (!returnReduce((fldPtr != nullptr), orOp<bool>()))
+    {
+        return false;
+    }
+
+    autoPtr<OFstream> osPtr;
+
     if (Pstream::master())
     {
-        OFstream os(outputName);
+        osPtr.reset(new OFstream(outputName));
+        osPtr->precision(precision_);
 
-        os << "# x y z " << fieldName_ << nl;
+        *(osPtr) << "# x y z " << fieldName_ << nl;
+    }
 
-        // Master
-        if (fldPtr)
-        {
-            writeField(os, *pointsPtr, *fldPtr);
-        }
 
-        // Slaves - recv
-        for (int slave=1; slave<Pstream::nProcs(); ++slave)
-        {
-            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
-            vectorField points(fromSlave);
-            Field<Type> fld(fromSlave);
-
-            writeField(os, points, fld);
-        }
+    if (applyFilter_)
+    {
+        writeListParallel(osPtr.ref(), *pointsPtr, values, parcelAddr_);
     }
     else
     {
-        // Slaves
-
-        OPstream toMaster(Pstream::commsTypes::blocking, Pstream::masterNo());
-
-        if (fldPtr)
-        {
-            toMaster
-                << *pointsPtr
-                << *fldPtr;
-        }
-        else
-        {
-            toMaster
-                << vectorField()
-                << Field<Type>();
-        }
+        writeListParallel(osPtr.ref(), *pointsPtr, values);
     }
 
     return true;
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/dataCloud b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/dataCloud
index 127d1ccf07..132e87de5e 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/dataCloud
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/dataCloud
@@ -1,5 +1,6 @@
 // -*- C++ -*-
 // Minimal example of using the dataCloud function object.
+
 dataCloud
 {
     type    dataCloud;
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/system/controlDict b/tutorials/lagrangian/reactingParcelFoam/filter/system/controlDict
index e54c5d788b..5a2262f241 100644
--- a/tutorials/lagrangian/reactingParcelFoam/filter/system/controlDict
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/system/controlDict
@@ -53,6 +53,7 @@ maxDeltaT       1;
 
 functions
 {
+    #include "dataCloud"
     #include "vtkCloud"
     #include "vtkWrite"
 }
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/system/dataCloud b/tutorials/lagrangian/reactingParcelFoam/filter/system/dataCloud
new file mode 100644
index 0000000000..51a99a2332
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/system/dataCloud
@@ -0,0 +1,42 @@
+// -*- C++ -*-
+
+dataCloud
+{
+    type    dataCloud;
+    libs    ("liblagrangianFunctionObjects.so");
+    log     true;
+
+    // Nothing happens before this anyhow
+    timeStart 0.5;
+
+    writeControl    writeTime;
+
+    // cloud   reactingCloud1;
+    clouds  ( ".*" );
+
+    // Field to output
+    field   d;
+
+    // Optional selection mechanism
+    selection
+    {
+        // Reduced number of output parcels
+        stride
+        {
+            action  add;
+            source  stride;
+            stride  4;
+        }
+
+        Umin
+        {
+            action  subtract;
+            source  field;
+            field   U;
+            accept  (less 0.2);
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/system/vtkCloud b/tutorials/lagrangian/reactingParcelFoam/filter/system/vtkCloud
index ce153de26e..451312c593 100644
--- a/tutorials/lagrangian/reactingParcelFoam/filter/system/vtkCloud
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/system/vtkCloud
@@ -28,6 +28,7 @@ cloudWrite
     //- Output directory name - Default postProcessing
     // directory       "VTK";
 
+    // Optional selection mechanism
     selection
     {
         all