From c90f775156c48a0c7f861cbb0d85efc0fec3c487 Mon Sep 17 00:00:00 2001 From: Hamidreza Date: Tue, 1 Jul 2025 18:18:53 +0330 Subject: [PATCH] rectMesh postProcess revisited --- .../PostprocessOperationAverage.cpp | 22 +++ .../PostprocessOperationAverage.hpp | 2 + .../operationFunctions.hpp | 2 +- .../postprocessOperation.cpp | 19 ++ .../postprocessOperation.hpp | 4 +- .../PostprocessComponent.cpp | 22 ++- .../PostprocessComponent.hpp | 4 +- .../PostprocessComponentGaussian.hpp | 4 +- .../region/regionFields/regionField.hpp | 6 + .../regionFields/regionFieldTemplate.cpp | 83 +++++++++ .../rectMeshRegionPoints/cellMapper.hpp | 56 ++++++ .../rectMeshRegionPoints.cpp | 170 ++++++++++++++---- .../rectMeshRegionPoints.hpp | 43 +++-- .../regionPoints/regionPoints.hpp | 17 ++ 14 files changed, 399 insertions(+), 55 deletions(-) create mode 100644 src/PostprocessData/region/regionPoints/rectMeshRegionPoints/cellMapper.hpp diff --git a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp index 07a4411e..46983bfb 100644 --- a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp +++ b/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 \ No newline at end of file diff --git a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp index 1840c0c4..de779d01 100644 --- a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp +++ b/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 diff --git a/src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp b/src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp index d9c23552..fae9f581 100644 --- a/src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp +++ b/src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp @@ -142,7 +142,7 @@ regionField executeFluctuation2Operation ) { const auto& regPoints = fieldAvg.regPoints(); - regionField processedField(regFieldName, regPoints, T{}); + regionField processedField(regFieldName+"_fluctuation2", regPoints, T{}); auto vols = regPoints.volumes(); for(uint32 reg =0; regbool + { + return arg.writeFieldToVtk(os); + }, + field + ); + } + return false; +} + uniquePtr postprocessOperation::create ( const dictionary &opDict, diff --git a/src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp b/src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp index 73aa3112..c932a7e3 100644 --- a/src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp +++ b/src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp @@ -99,7 +99,7 @@ public: private: /// Dictionary containing operation-specific parameters. - pFlow::dictionary operationDict_; + pFlow::dictionary operationDict_; /// This Threshold is used to exclude the regions which contain /// fewer than this value. @@ -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. diff --git a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp index 90dcff16..0bc83ec5 100644 --- a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp +++ b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp @@ -152,11 +152,25 @@ bool pFlow::postprocessData::PostprocessComponent } else { - notImplementedFunction; - return false; - } - + word chNum = real2FixedStripZeros(database().time().currentTime() *1000000, 0); + fileSystem file = parDir + (name() +"-"+chNum+".vtk"); + auto osPtr = makeUnique(file); + + regPoints().write(osPtr()); + + for(auto& operation:operatios_) + { + if(!operation->write(osPtr())) + { + fatalErrorInFunction + <<"Error occurred in writing operation defined in dict " + << operation->operationDict() + < volumeFactor_; - bool executed_{false}; + bool executed_{false}; dictionaryList operationDicts_; @@ -122,8 +122,6 @@ public: }; - - } #include "PostprocessComponent.cpp" diff --git a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp index 2f5cbd3e..ac8b33fa 100644 --- a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp +++ b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp @@ -47,8 +47,8 @@ public: : PostprocessComponent(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(); diff --git a/src/PostprocessData/region/regionFields/regionField.hpp b/src/PostprocessData/region/regionFields/regionField.hpp index 9094cd34..177dddef 100644 --- a/src/PostprocessData/region/regionFields/regionField.hpp +++ b/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(); diff --git a/src/PostprocessData/region/regionFields/regionFieldTemplate.cpp b/src/PostprocessData/region/regionFields/regionFieldTemplate.cpp index 408107c6..42e61d48 100644 --- a/src/PostprocessData/region/regionFields/regionFieldTemplate.cpp +++ b/src/PostprocessData/region/regionFields/regionFieldTemplate.cpp @@ -1,3 +1,6 @@ + + + namespace pFlow::postprocessData { @@ -9,6 +12,86 @@ regionField::regionField( : field_(name, "regionFieldValue", rPoints.size(), rPoints.size(), defaultVal), regionPoints_(rPoints) + {} +template +inline +bool regionField::writeFieldToVtk(iOstream& os)const +{ + fatalErrorInFunction<< "This type is not supported for vtk conversion:"<< + field_.typeName()< +inline +bool regionField::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 +inline +bool regionField::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 +inline +bool regionField::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 + ( + "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() || 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("nx", 1); - ny = rectMeshInfo.getValMax("ny", 1); - nz = rectMeshInfo.getValMax("nz", 1); - // boxRegion_ = box(minP, maxP); - cells_ = uint32x3(nx, ny, nz); - uint32 nCells = nx * ny * nz; + auto nx = rectMeshInfo.getValMax("nx", 1); + auto ny = rectMeshInfo.getValMax("ny", 1); + auto nz = rectMeshInfo.getValMax("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(i) + 0.5) * dx, + ( static_cast(j) + 0.5) * dy, + ( static_cast(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++) { diff --git a/src/PostprocessData/region/regionPoints/rectMeshRegionPoints/rectMeshRegionPoints.hpp b/src/PostprocessData/region/regionPoints/rectMeshRegionPoints/rectMeshRegionPoints.hpp index f8bbaa8b..ee5c96da 100644 --- a/src/PostprocessData/region/regionPoints/rectMeshRegionPoints/rectMeshRegionPoints.hpp +++ b/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 pointsOnCells_; + + uniquePtr> pointsBeyoundCells_; + /// Indices of points that are selected by this region - Vector selectedPoints_; + Vector& 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 @@ -163,7 +183,7 @@ public: fatalExit; } - return span(selectedPoints_[elem].data(), selectedPoints_[elem].size()); + return span(selectedPoints_[elem].data(), selectedPoints_[elem].size()); } /** @@ -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; } /** diff --git a/src/PostprocessData/region/regionPoints/regionPoints/regionPoints.hpp b/src/PostprocessData/region/regionPoints/regionPoints/regionPoints.hpp index 54cca825..1e8874bb 100644 --- a/src/PostprocessData/region/regionPoints/regionPoints/regionPoints.hpp +++ b/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"); @@ -78,11 +83,23 @@ public: /// @brief size of elements 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