rectMesh postProcess revisited

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp

@@ -178,4 +178,26 @@ bool PostprocessOperationAverage::write(const fileSystem &parDir) const
     return true;
 }
 
+bool PostprocessOperationAverage::write(iOstream &os) const
+{
+    if(! postprocessOperation::write(os))
+    {
+        return false;
+    }
+    if(!calculateFluctuation2_())
+    {
+        return true;
+    }
+    
+    return 
+        std::visit
+        (
+            [&](auto&& arg)->bool
+            {
+                return arg.writeFieldToVtk(os);
+            },
+            fluctuation2FieldPtr_()
+        );
+}
+
 } // namespace pFlow::postprocessData

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp

@@ -195,6 +195,8 @@ public:
     /// write to os stream
     bool write(const fileSystem &parDir)const override;
 
+    bool write(iOstream& os)const override;
+
 
     /// @brief Execute average operation on field values
     /// @param weights Weight factors for particles

src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp

@@ -142,7 +142,7 @@ regionField<T> executeFluctuation2Operation
 )
 {
     const auto& regPoints = fieldAvg.regPoints();
-    regionField<T> processedField(regFieldName, regPoints, T{});
+    regionField<T> processedField(regFieldName+"_fluctuation2", regPoints, T{});
     auto vols = regPoints.volumes();
 
     for(uint32 reg =0; reg<regPoints.size(); reg++)

src/PostprocessData/operation/postprocessOperation/postprocessOperation.cpp

@@ -124,6 +124,25 @@ bool postprocessOperation::write(const fileSystem &parDir) const
     return true;
 }
 
+bool postprocessOperation::write(iOstream& os)const
+{
+    if(!regPoints().writeToSameTimeFile())
+    {
+        const auto& field = processedField();
+
+        return 
+        std::visit
+        (
+            [&](auto&& arg)->bool
+            {
+                return arg.writeFieldToVtk(os);
+            },
+            field
+        );
+    }
+    return false;
+}
+
 uniquePtr<postprocessOperation> postprocessOperation::create
 (
     const dictionary &opDict,

src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp

@@ -255,7 +255,7 @@ public:
     /// write the result to output stream (possibly a file)
     /// @param os Output stream to write the result.
     virtual 
-    bool write(iOstream& os)const {return true;}
+    bool write(iOstream& os)const;
 
     /// Create the polymorphic object using the virtual constructor.
     /// @param opDict Dictionary containing operation-specific parameters.

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp

@@ -152,11 +152,25 @@ bool pFlow::postprocessData::PostprocessComponent<RegionType, ProcessMethodType>
     }
     else
     {
-        notImplementedFunction;
+        word chNum = real2FixedStripZeros(database().time().currentTime() *1000000, 0);
+        fileSystem  file = parDir + (name() +"-"+chNum+".vtk");
+
+        auto osPtr = makeUnique<oFstream>(file);
+
+        regPoints().write(osPtr());
+
+        for(auto& operation:operatios_)
+        {
+            if(!operation->write(osPtr()))
+            {
+                fatalErrorInFunction
+                    <<"Error occurred in writing operation defined in dict "
+                    << operation->operationDict()
+                    <<endl;
                 return false;
             }
-    
-
+        }   
+    }
 
     return true;
 }

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.hpp

@@ -122,8 +122,6 @@ public:
 
 };
 
-
-
 }
 
 #include "PostprocessComponent.cpp"

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp

@@ -47,8 +47,8 @@ public:
     :
         PostprocessComponent<RegionType,GaussianDistribution>(dict, fieldsDB, defaultTimeControl)
     {
-        /// initializes the Gaussian distribution for all elements of region
-        //const uint32 n = this->regPoints().size();
+
+        this->regPoints().setRegionExtension(2);
         auto d = this->regPoints().eqDiameters();
         auto c = this->regPoints().centers();
         auto& regs = this->regionProecessMethod();

src/PostprocessData/region/regionFields/regionField.hpp

@@ -24,6 +24,7 @@ Licence:
 #include "types.hpp"
 #include "regionPoints.hpp"
 #include "Field.hpp"
+#include "cellMapper.hpp"
 
 namespace pFlow::postprocessData
 {
@@ -101,6 +102,11 @@ public:
         return field_.size();
     }
 
+    uint32x3 shape()const
+    {
+        return regionPoints_.shape();
+    }
+
     bool empty()const
     {
         return field_.empty();

src/PostprocessData/region/regionFields/regionFieldTemplate.cpp

@@ -1,3 +1,6 @@
+
+
+
 namespace pFlow::postprocessData
 {
 
@@ -9,6 +12,86 @@ regionField<T>::regionField(
 :
     field_(name, "regionFieldValue", rPoints.size(), rPoints.size(), defaultVal),
     regionPoints_(rPoints)
+
 {}
 
+template<typename T>
+inline
+bool regionField<T>::writeFieldToVtk(iOstream& os)const
+{
+    fatalErrorInFunction<< "This type is not supported for vtk conversion:"<<
+	field_.typeName()<<endl;
+	fatalExit;
+    return false;
+}
+
+template<>
+inline
+bool regionField<real>::writeFieldToVtk(iOstream& os)const
+{
+    os<<"FIELD FieldData 1 " << field_.name() << " 1 "<< field_.size() << " float\n";
+    
+    
+    const auto mapper = cellMapper{shape()};
+
+    for(uint32 k=0; k<mapper.nz(); k++)
+    {
+        for(uint32 j=0; j<mapper.ny(); j++)
+        {
+            for(uint32 i=0; i<mapper.nx(); i++)
+            {
+                os<< field_[ mapper(i,j,k) ]<<'\n';
+            }
+        }
+    }
+    os<<endl;
+    return true;
+}
+
+template<>
+inline
+bool regionField<realx3>::writeFieldToVtk(iOstream& os)const
+{
+    os<<"FIELD FieldData 1 " << field_.name() << " 3 "<< field_.size() << " float\n";
+    
+    
+    const auto mapper = cellMapper{shape()}; 
+
+    for(uint32 k=0; k<mapper.nz(); k++)
+    {
+        for(uint32 j=0; j<mapper.ny(); j++)
+        {
+            for(uint32 i=0; i<mapper.nx(); i++)
+            {
+                os<<field_[mapper(i,j,k)].x()<<' '<<field_[mapper(i,j,k)].y()<<' '<<field_[mapper(i,j,k)].z()<<'\n';
+            }
+        }
+    }
+    os<<endl;
+    return true;
+}
+
+template<>
+inline
+bool regionField<uint32>::writeFieldToVtk(iOstream& os)const
+{
+    os<<"FIELD FieldData 1 " << field_.name() << " 1 "<< field_.size() << " int\n";
+    
+    
+    const auto mapper = cellMapper{shape()} ; 
+
+    for(uint32 k=0; k<mapper.nz(); k++)
+    {
+        for(uint32 j=0; j<mapper.ny(); j++)
+        {
+            for(uint32 i=0; i<mapper.nx(); i++)
+            {
+                os<< field_[ mapper(i,j,k) ]<<'\n';
+            }
+        }
+    }
+    os<<endl;
+    return true;
+}
+
 } // End namespace pFlow::postprocessData

src/PostprocessData/region/regionPoints/rectMeshRegionPoints/cellMapper.hpp

@@ -0,0 +1,56 @@
+
+#ifndef __cellMapper_hpp__
+#define __cellMapper_hpp__
+
+#include "types.hpp"
+
+namespace pFlow::postprocessData
+{
+
+struct cellMapper
+{
+    uint32x3 cells_;
+
+    cellMapper()
+    :
+        cells_()
+    {}
+
+    explicit cellMapper(uint32x3 cells)
+    :
+        cells_(cells)
+    {}
+
+    cellMapper(const cellMapper&) = default;
+    cellMapper(cellMapper&&) = default;
+    cellMapper& operator=(const cellMapper&) = default;
+    cellMapper& operator=(cellMapper&&) = default;
+    ~cellMapper() = default;
+
+    inline
+    uint32 operator()(uint32 i, uint32 j, uint32 k)const
+    {
+        return (k*(cells_.y()*cells_.x()))+j*cells_.x() + i;
+    }
+
+    uint32 nx()const
+    {
+        return cells_.x();
+    }
+
+    uint32 ny()const
+    {
+        return cells_.y();
+    }
+
+    uint32 nz()const
+    {
+        return cells_.z();
+    }
+
+};
+
+
+} //pFlow::postprocessData
+
+#endif //__cellMapper_hpp__

src/PostprocessData/region/regionPoints/rectMeshRegionPoints/rectMeshRegionPoints.cpp

@@ -2,6 +2,80 @@
 #include "fieldsDataBase.hpp"
 #include "numericConstants.hpp"
 
+void pFlow::postprocessData::rectMeshRegionPoints::findPointsBeyoundCells()
+{
+    // check if pointsBeyoundCells_ is initialized
+    if(!pointsBeyoundCells_)
+    {
+        pointsBeyoundCells_ = makeUnique<decltype(pointsOnCells_)>
+        (
+            "selectedPoints2",
+            this->size()
+        );
+    }
+
+    // get the reference to pointsBeyoundCells_ and clear it
+    auto& selectedPoints = pointsBeyoundCells_();
+
+    // point positions are obtained from the database
+    const auto points = database().updatePoints();
+
+    // iterate through all cells to find points that are within the search radius
+    for(int32 i=0; i<mapper_.nx(); i++)
+    {
+        for(int32 j=0; j<mapper_.ny(); j++)
+        {
+            for(int32 k=0; k<mapper_.nz(); k++)
+            {
+                uint32 cellIndex = mapper_(i,j,k);
+
+                // copy the points in the center cell
+                auto& cellIndices = selectedPoints[cellIndex];
+                cellIndices.clear();
+                
+                
+                if(pointsOnCells_[cellIndex].empty())
+                    continue;
+
+                const auto cellCenter = centerPoints_[cellIndex];
+                const auto rad = 0.5*diameter_[cellIndex];
+                
+                for(int32 ii=-2; ii <= 2; ++ii)
+                {
+                    for(int32 jj=-2; jj <= 2; ++jj)
+                    {
+                        for(int32 kk=-2; kk <= 2; ++kk)
+                        {
+
+                            int32 ni = i + ii;
+                            int32 nj = j + jj;
+                            int32 nk = k + kk;
+                            if(ni < 0 || nj < 0 || nk < 0)
+                                continue;
+                            
+                            if(ni >= mapper_.nx() || nj >= mapper_.ny() || nk >= mapper_.nz())
+                                continue;
+                            
+                            uint32 neighborIndex = mapper_(ni, nj, nk);
+
+                            const auto& neighborPoints = pointsOnCells_[neighborIndex];
+
+                            for(auto nIndx : neighborPoints)
+                            {
+                                if( (points[nIndx]-cellCenter).length() < rad )
+                                {
+                                    cellIndices.push_back(nIndx);
+                                }
+                            }
+
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
 
 pFlow::postprocessData::rectMeshRegionPoints::rectMeshRegionPoints
 (
@@ -9,23 +83,28 @@ pFlow::postprocessData::rectMeshRegionPoints::rectMeshRegionPoints
     fieldsDataBase &fieldsDataBase
 )
 :
-    regionPoints(dict, fieldsDataBase)
+    regionPoints(dict, fieldsDataBase),
+    boxRegion_(dict.subDict("rectMeshInfo")),
+    pointsOnCells_("selectedPoints"),
+    selectedPoints_(pointsOnCells_)
 {
     const auto& rectMeshInfo = dict.subDict("rectMeshInfo");
-    boxRegion_ = box(rectMeshInfo.subDict("boxInfo"));
-    nx = rectMeshInfo.getValMax<uint32>("nx", 1);
-    ny = rectMeshInfo.getValMax<uint32>("ny", 1);
-    nz = rectMeshInfo.getValMax<uint32>("nz", 1);
     
-    // boxRegion_ = box(minP, maxP);
-    cells_ = uint32x3(nx, ny, nz);
-    uint32 nCells = nx * ny * nz;
+    auto nx = rectMeshInfo.getValMax<uint32>("nx", 1);
+    auto ny = rectMeshInfo.getValMax<uint32>("ny", 1);
+    auto nz = rectMeshInfo.getValMax<uint32>("nz", 1);
     
-    volumes_.resize(nCells, boxRegion_.volume() / nCells);
-    diameter_.resize(nCells, 2 * pow(3 * boxRegion_.volume() / nCells / 4.0 / Pi, 1.0 / 3.0));
-    selectedPoints_.resize(nCells);
+    mapper_ = cellMapper(uint32x3(nx, ny, nz));
+    uint32 nCells = mapper_.nx() * mapper_.ny() * mapper_.nz();
+    real vol = boxRegion_.volume() / nCells;
+    volumes_.resize(nCells, vol);
+    diameter_.resize(nCells, 2 * pow(3 * vol / 4.0 / Pi, 0.3333333));
+    pointsOnCells_.resize(nCells);
     centerPoints_.resize(nCells);
 
+    real dx = (boxRegion_.maxPoint().x() - boxRegion_.minPoint().x()) / mapper_.nx();
+    real dy = (boxRegion_.maxPoint().y() - boxRegion_.minPoint().y()) / mapper_.ny();
+    real dz = (boxRegion_.maxPoint().z() - boxRegion_.minPoint().z()) / mapper_.nz();
     
     for(uint32 i = 0; i < nx; ++i)
     {
@@ -33,57 +112,86 @@ pFlow::postprocessData::rectMeshRegionPoints::rectMeshRegionPoints
         {
             for(uint32 k = 0; k < nz; ++k)
             {
-                // calculate the center point of each cell
-                uint32 index = (i * ny + j) * nz + k;
                 realx3 center = boxRegion_.minPoint() + 
                     realx3(
-                        (i + 0.5) * (boxRegion_.maxPoint().x() - boxRegion_.minPoint().x()) / nx,
-                        (j + 0.5) * (boxRegion_.maxPoint().y() - boxRegion_.minPoint().y()) / ny,
-                        (k + 0.5) * (boxRegion_.maxPoint().z() - boxRegion_.minPoint().z()) / nz
+                        ( static_cast<real>(i) + 0.5) * dx,
+                        ( static_cast<real>(j) + 0.5) * dy,
+                        ( static_cast<real>(k) + 0.5) * dz
                     );
-                centerPoints_[index] = center;
+                centerPoints_[mapper_(i, j, k)] = center;
             }
         }
     }
 }
 
+void pFlow::postprocessData::rectMeshRegionPoints::setRegionExtension(real ext)
+{
+    regionExtension_ = max(one, ext);
+    real vf = pow(regionExtension_, 3);
+    for(auto& v:volumes_)
+    {
+        v *= vf;
+    }
+
+    for(auto& d:diameter_)
+    {
+        d *= regionExtension_;
+    }
+}
+
 bool pFlow::postprocessData::rectMeshRegionPoints::update()
 {
     const auto points = database().updatePoints();
-    for (auto& elem : selectedPoints_)
+    for (auto& elem : pointsOnCells_)
     {
         elem.clear();
     }
  
-    real dx = (boxRegion_.maxPoint().x() - boxRegion_.minPoint().x()) / cells_.x();
-    real dy = (boxRegion_.maxPoint().y() - boxRegion_.minPoint().y()) / cells_.y();
-    real dz = (boxRegion_.maxPoint().z() - boxRegion_.minPoint().z()) / cells_.z();
+    real dx = (boxRegion_.maxPoint().x() - boxRegion_.minPoint().x()) / mapper_.nx();
+    real dy = (boxRegion_.maxPoint().y() - boxRegion_.minPoint().y()) / mapper_.ny();
+    real dz = (boxRegion_.maxPoint().z() - boxRegion_.minPoint().z()) / mapper_.nz();
 
     for (uint32 i = 0; i < points.size(); ++i)
     {
         if(boxRegion_.isInside(points[i]))
         {
-            uint indexX = (points[i] - boxRegion_.minPoint()).x() / dx;
-            uint indexY = (points[i] - boxRegion_.minPoint()).y() / dy;
-            uint indexZ = (points[i] - boxRegion_.minPoint()).z() / dz;
-            uint cellIndex = (indexX * cells_.y() + indexY) * cells_.z() + indexZ;
-            selectedPoints_[cellIndex].push_back(i);
+            uint32 indexX = (points[i] - boxRegion_.minPoint()).x() / dx;
+            uint32 indexY = (points[i] - boxRegion_.minPoint()).y() / dy;
+            uint32 indexZ = (points[i] - boxRegion_.minPoint()).z() / dz;
+            pointsOnCells_[mapper_(indexX, indexY, indexZ)].push_back(i);
         }
     }
+
+    // search beyound cells is not required
+    if( equal(regionExtension_,one))
+    {
+        selectedPoints_ = pointsOnCells_;
+        return true;
+    }
+
+    // search beyound cells is required
+    findPointsBeyoundCells();
+    selectedPoints_ = pointsBeyoundCells_();
+
     return true;
 }
 
 bool pFlow::postprocessData::rectMeshRegionPoints::write(iOstream &os) const
 {
+
+    auto [x, y , z] = boxRegion_.minPoint();
+    auto [nx, ny, nz] = mapper_.cells_;
+    real dx = (boxRegion_.maxPoint().x() - boxRegion_.minPoint().x()) / mapper_.nx();
+    real dy = (boxRegion_.maxPoint().y() - boxRegion_.minPoint().y()) / mapper_.ny();
+    real dz = (boxRegion_.maxPoint().z() - boxRegion_.minPoint().z()) / mapper_.nz();
+
+
     os << "# vtk DataFile Version 3.0" << endl;
     os << "postProcessData" << endl;
     os << "ASCII" << endl;
     os << "DATASET RECTILINEAR_GRID" << endl;
 	os << "DIMENSIONS " << nx + 1 << " " << ny + 1 << " " << nz + 1 << endl;
 	
-    auto [x, y , z] = boxRegion_.minPoint();
-    auto [dx, dy, dz] = (boxRegion_.maxPoint() - boxRegion_.minPoint()) / realx3(nx, ny,  nz);
-
     os << "X_COORDINATES " << nx + 1 << " float\n";
     for(int32 i = 0; i < nx + 1; i++)
     {

src/PostprocessData/region/regionPoints/rectMeshRegionPoints/rectMeshRegionPoints.hpp

@@ -27,7 +27,6 @@ Licence:
  * It inherits from regionPoints and implements all required virtual methods.
  * 
  * @see regionPoints
- * @see rectMesh
  * @see fieldsDataBase
  */
 
@@ -37,11 +36,13 @@ Licence:
 #include "regionPoints.hpp"
 #include "box.hpp"
 #include "Vectors.hpp"
-#include "cells.hpp"
+#include "cellMapper.hpp"
+
 
 namespace pFlow::postprocessData
 {
 
+
 class rectMeshRegionPoints
 :
     public regionPoints
@@ -51,11 +52,8 @@ private:
     /// box object defining the region for point selection
     box  boxRegion_;
 
-    /// Number of cells in each direction
-    uint32 nx, ny, nz; 
-
     /// store the cells that are inside the box region
-    uint32x3   cells_;
+    cellMapper      mapper_;
 
     /// Center points of each cell in the rectMesh region
     realx3Vector    centerPoints_;
@@ -66,8 +64,15 @@ private:
     /// Diameter of each cell in the rectMesh region
     realVector    diameter_;
 
+    Vector<uint32Vector>            pointsOnCells_;
+
+    uniquePtr<Vector<uint32Vector>> pointsBeyoundCells_;
+
     /// Indices of points that are selected by this region
-    Vector<uint32Vector>    selectedPoints_;
+    Vector<uint32Vector>&           selectedPoints_;
+
+
+    void findPointsBeyoundCells();
 
 public:
 
@@ -95,14 +100,29 @@ public:
     }
 
     /**
-     * @brief Check if the region is empty
-     * @return Always returns false
+     * return the shape of the field 
+     */
+    uint32x3 shape()const override
+    {
+        return mapper_.cells_;
+    }
+
+    const cellMapper&  mapper()const
+    {
+        return mapper_;
+    }
+
+    /**
+     * @brief Update the points selected by this region
+     * @return True if update was successful
      */
     bool empty()const override
     {
         return volumes_.empty();
     }
 
+    void setRegionExtension(real ext) override;
+
     /**
      * @brief Get the volume of the rectMesh region
      * @return A span containing the volume of the region
@@ -174,11 +194,10 @@ public:
 
     /**
      * @brief Determine if data should be written to the same time file
-     * @return Always returns true
      */
     bool writeToSameTimeFile()const override
     {
-        return true;
+        return false;
     }
 
     /**

src/PostprocessData/region/regionPoints/regionPoints/regionPoints.hpp

@@ -54,6 +54,11 @@ class regionPoints
     /// Reference to the fields database containing simulation data
     fieldsDataBase& fieldsDataBase_;
 
+protected:
+
+    /// extends the search radius to a distance farther than the region
+    real            regionExtension_ = 1.0;
+
 public:
 
     TypeInfo("regionPoints");
@@ -79,10 +84,22 @@ public:
     virtual 
     uint32 size()const = 0;
 
+    virtual 
+    uint32x3 shape()const
+    {
+      return uint32x3(size(), 1u, 1u);
+    }
+    
     /// @brief  check if the region is empty
     virtual 
     bool empty()const = 0;  
     
+    /// by default it does nothing
+    /// But, it can be used for the methods that needs to search for 
+    /// particles which are beyound the region 
+    virtual void setRegionExtension(real ext)
+    {}
+
     /// @brief  volume of elements
     /// @return sapn for accessing the volume of elements 
     virtual