diff --git a/src/phasicFlow/triSurface/multiTriSurface.cpp b/src/phasicFlow/triSurface/multiTriSurface.cpp index b2c04995..e73f4382 100644 --- a/src/phasicFlow/triSurface/multiTriSurface.cpp +++ b/src/phasicFlow/triSurface/multiTriSurface.cpp @@ -21,139 +21,6 @@ Licence: #include "multiTriSurface.hpp" -/*void pFlow::multiTriSurface::calculateVars() -{ - numSurfaces_ = surfaceNames_.size(); - - // make sure the host and device are sync - lastPointIndex_.syncViews(); - - surfaceNumPoints_.clear(); - - int32 last = 0; - for(auto& pi:lastPointIndex_) - { - surfaceNumPoints_.push_back(pi+1-last); - last = pi+1; - } - - // update views - surfaceNumPoints_.syncViews(); - - surfaceNumVertices_.clear(); - last = 0; - - lastVertexIndex_.syncViews(); - - - for(auto& vi:lastVertexIndex_) - { - surfaceNumVertices_.push_back(vi+1-last); - last = vi+1; - } - surfaceNumVertices_.syncViews(); - - pointsStartPos_.reallocate(surfaceNumPoints_.capacity()); - pointsStartPos_.clear(); - - pointsStartPos_.push_back(0); - - ForAll( i, surfaceNumPoints_) - { - if(i==0)continue; - - pointsStartPos_.push_back( pointsStartPos_[i-1] + surfaceNumPoints_[i-1] ); - } - - pointsStartPos_.syncViews(); - - - verticesStartPos_.reallocate(surfaceNumVertices_.capacity()); - verticesStartPos_.clear(); - - verticesStartPos_.push_back(0); - ForAll(i, surfaceNumVertices_) - { - if(i==0)continue; - verticesStartPos_.push_back(verticesStartPos_[i-1] + surfaceNumVertices_[i-1]); - } - verticesStartPos_.syncViews(); - - -}*/ - -/*pFlow::multiTriSurface::multiTriSurface() -: - triSurface(), - lastPointIndex_("lastPointIndex", "lastPointIndex"), - lastVertexIndex_("lastVertexIndex", "lastVertexIndex"), - surfaceNames_("surfaceNames", "surfaceNames") -{ - calculateVars(); -}*/ - -/*bool pFlow::multiTriSurface::addTriSurface -( - const word& name, - const triSurface& tSurf -) -{ - - const auto& newPoints = tSurf.points(); - const auto& newVertices = tSurf.vertices(); - const auto& newAreas = tSurf.area(); - - // - - - points_.append(newPoints); - - - // add new vertices to the existing one - auto vOldSize = vertices_.size(); - auto vNewSize = vOldSize + newVertices.size(); - vertices_.resize(vNewSize); - area_.resize(vNewSize); - - auto verVec = vertices_.deviceViewAll(); - auto areaVec = area_.deviceViewAll(); - - auto newVerVec = newVertices.deviceViewAll(); - auto newArea = newAreas.deviceViewAll(); - - auto maxIdx = maxIndex(); - - Kokkos::parallel_for( - "multiTriSurface::addTriSurface", - newVertices.size(), - LAMBDA_HD(int32 i){ - verVec[vOldSize+i] = newVerVec[i]+(maxIdx+1); - areaVec[vOldSize+i] = newArea[i]; - } - ); - Kokkos::fence(); - - if( !check() ) - { - fatalErrorInFunction<< - "the indices and number of points do not match. \n"; - return false; - } - - lastPointIndex_.push_back(points_.size()-1); - lastPointIndex_.syncViews(); - - lastVertexIndex_.push_back(vertices_.size()-1); - lastVertexIndex_.syncViews(); - - surfaceNames_.push_back(name); - - calculateVars(); - - return true; - -}*/ - pFlow::multiTriSurface::multiTriSurface ( const objectFile &obj, @@ -185,7 +52,7 @@ pFlow::multiTriSurface::multiTriSurface } -bool pFlow::multiTriSurface::appendTriSurface +bool pFlow::multiTriSurface::appendSurface ( const word &name, const realx3x3Vector &triangles @@ -194,7 +61,7 @@ bool pFlow::multiTriSurface::appendTriSurface uint32 start = size(); uint32 pointStart = numPoints(); - if(!triSurface::appendTriSurface(triangles)) + if(!triSurface::append(triangles)) { fatalErrorInFunction; return false; @@ -236,45 +103,3 @@ bool pFlow::multiTriSurface::write return triSurface::write(os,iop); } - -/*pFlow::real3* - pFlow::multiTriSurface::beginSurfacePoint -( - unit i -) -{ - if(i== 0) return points_.data(); - if(i>=numSurfaces())return points_.data()+numPoints(); - return points_.data()+lastPointIndex_[i-1]+1; -} - -const pFlow::real3* - pFlow::multiTriSurface::beginSurfacePoint -( - unit i -)const -{ - if(i== 0) return points_.data(); - if(i>=numSurfaces())return points_.data()+numPoints(); - return points_.data()+lastPointIndex_[i-1]+1; -} - -pFlow::real3* - pFlow::multiTriSurface::endSurfacePoint -( - unit i -) -{ - if(i>=numSurfaces())return points_.data()+numPoints(); - return points_.data()+lastPointIndex_[i]+1; -} - -const pFlow::real3* - pFlow::multiTriSurface::endSurfacePoint -( - unit i -)const -{ - if(i>=numSurfaces())return points_.data()+numPoints(); - return points_.data()+lastPointIndex_[i]+1; -}*/ \ No newline at end of file diff --git a/src/phasicFlow/triSurface/multiTriSurface.hpp b/src/phasicFlow/triSurface/multiTriSurface.hpp index d00130a7..b95977d3 100644 --- a/src/phasicFlow/triSurface/multiTriSurface.hpp +++ b/src/phasicFlow/triSurface/multiTriSurface.hpp @@ -29,7 +29,6 @@ Licence: namespace pFlow { - class multiTriSurface : public triSurface, @@ -39,10 +38,6 @@ private: subSurfaceList subSurfaces_; - // - the index of last point of each triSurface - - //void calculateVars(); - public: // - type info @@ -50,33 +45,42 @@ public: //// - Constructors - // + /// @brief Construct from objectFile and owner repository. + /// This is mainly used for reading from file. + /// @param obj object file + /// @param owner owner repository multiTriSurface(const objectFile& obj, repository* owner); + /// @brief Construct from another multiTriSurface multiTriSurface( const objectFile& objf, repository* owner, const multiTriSurface& surf); + /// @brief Copy construct (default) multiTriSurface(const multiTriSurface&) = default; + /// @brief Copy assignment (default) multiTriSurface& operator = (const multiTriSurface&) = default; + /// @brief No move construct multiTriSurface(multiTriSurface&&) = delete; + /// @brief No move assignment multiTriSurface& operator = (multiTriSurface&&) = delete; + /// @brief default destructor ~multiTriSurface() override = default; //// - Methods - //bool addTriSurface(const word& name, const triSurface& tSurf); - - bool appendTriSurface(const word& name, const realx3x3Vector& vertices); + bool appendSurface( + const word& name, + const realx3x3Vector& vertices); uint32 numSurfaces()const { - return subSurfaces_.size(); + return static_cast(subSurfaces_.size()); } const subSurfaceList& subSurfaces()const @@ -84,7 +88,7 @@ public: return subSurfaces_; } - rangeU32 subSurfaceRange(uint32 nSub) + rangeU32 subSurfaceRange(uint32 nSub)const { if( !(nSub < numSurfaces() ) ) { @@ -106,55 +110,50 @@ public: subSurfaces_[nSub].pointEnd()}; } - /*void clear() + /// @brief Clear the content of object + void clear() { triSurface::clear(); - - lastPointIndex_.clear(); - surfaceNames_.clear(); - }*/ + subSurfaces_.clear(); + } - /*const auto& pointsStartPos()const + uint32 subSurfaceSize(uint32 nSub)const { - return pointsStartPos_; + if(nSub> (iIstream & is, multiTriSurface & tri ) -{ - if(!tri.readMultiTriSurface(is)) - { - ioErrorInFile(is.name(), is.lineNumber())<< - " error in reading multiTriSurface from file.\n"; - fatalExit; - } - return is; -}*/ inline iOstream& operator << (iOstream& os, const multiTriSurface& tri) { @@ -189,7 +178,7 @@ inline iOstream& operator << (iOstream& os, const multiTriSurface& tri) return os; } -} +} // pFlow #endif diff --git a/src/phasicFlow/triSurface/triSurface.cpp b/src/phasicFlow/triSurface/triSurface.cpp index 6d3fea97..4daf1387 100644 --- a/src/phasicFlow/triSurface/triSurface.cpp +++ b/src/phasicFlow/triSurface/triSurface.cpp @@ -20,12 +20,12 @@ Licence: #include "triSurface.hpp" -#include "Vectors.hpp" -#include "triangleFunctions.hpp" #include "error.hpp" #include "iOstream.hpp" +#include "Vectors.hpp" +#include "triSurfaceKernels.hpp" + -//#include "triSurfaceKernels.hpp" namespace pFlow { @@ -104,75 +104,23 @@ bool convertToTriSurfaceComponents } - -pFlow::triSurface::triSurface -( - const objectFile &obj, - repository *owner -) -: - IOobject - ( - obj, - IOPattern::AllProcessorsSimilar, - owner - ), - points_("points", "points"), - vertices_("vertices", "vertices"), - area_("area", "area"), - normals_("normals","normals") +bool pFlow::triSurface::calculateNormals() { + return pFlow::triSurfaceKernels::calculateNormals( + points_, + vertices_, + normals_); } -pFlow::triSurface::triSurface -( - const objectFile &objf, - repository* owner, - const triSurface &surf -) -: - IOobject - ( - objectFile - ( - objf.name(), - objf.localPath(), - objectFile::READ_NEVER, - objf.wFlag() - ), - IOPattern::AllProcessorsSimilar, - owner - ), - points_(surf.points_), - vertices_(surf.vertices_), - area_(surf.area_), - normals_(surf.normals_) -{} - -pFlow::triSurface::triSurface( - const realx3x3Field_H &triangles, - repository *owner) - : IOobject( - objectFile( - triangles.name(), - "", - IOobject::READ_NEVER, - IOobject::WRITE_ALWAYS), - IOPattern::AllProcessorsSimilar, - owner), - points_("points", "points"), - vertices_("vertices", "vertices"), - area_("area", "area"), - normals_("normals", "normals") +bool pFlow::triSurface::calculateArea() { - if( !appendTriSurface(triangles) ) - { - fatalExit; - } - + return pFlow::triSurfaceKernels::calculateArea( + points_, + vertices_, + area_); } -bool pFlow::triSurface::appendTriSurface +bool pFlow::triSurface::append ( const realx3x3Field_H &triangles ) @@ -207,7 +155,7 @@ bool pFlow::triSurface::appendTriSurface return true; } -bool pFlow::triSurface::appendTriSurface(const realx3x3Vector &triangles) +bool pFlow::triSurface::append(const realx3x3Vector &triangles) { uint32 basePointIndex = numPoints(); @@ -239,6 +187,73 @@ bool pFlow::triSurface::appendTriSurface(const realx3x3Vector &triangles) return true; } +pFlow::triSurface::triSurface +( + const objectFile &obj, + repository *owner +) +: + IOobject + ( + obj, + IOPattern::AllProcessorsSimilar, + owner + ) +{ + // this constructor is used by multiTrisurface and read operation is + // done in that class +} + +pFlow::triSurface::triSurface +( + const objectFile &objf, + repository* owner, + const triSurface &surf +) +: + IOobject + ( + objectFile + ( + objf.name(), + objf.localPath(), + objectFile::READ_NEVER, + objf.wFlag() + ), + IOPattern::AllProcessorsSimilar, + owner + ), + points_(surf.points_), + vertices_(surf.vertices_), + area_(surf.area_), + normals_(surf.normals_) +{} + +pFlow::triSurface::triSurface +( + const realx3x3Field_H &triangles, + repository *owner +) +: + IOobject + ( + objectFile( + triangles.name(), + "", + IOobject::READ_NEVER, + IOobject::WRITE_ALWAYS), + IOPattern::AllProcessorsSimilar, + owner + ) +{ + if( !append(triangles) ) + { + fatalExit; + } +} + + + bool pFlow::triSurface::read(iIstream &is, const IOPattern &iop) { points_.clear(); @@ -257,8 +272,10 @@ bool pFlow::triSurface::read(iIstream &is, const IOPattern &iop) return false; } - WARNING<<"You should calculate area and normal after reading "< dPoints_; + + deviceViewType1D dVectices_; +public: + + INLINE_FUNCTION_H + triangleAccessor( + uint32 numPoints, + deviceViewType1D points, + uint32 numTrianlges, + deviceViewType1D vertices ) + : + numPoints_(numPoints), + numTriangles_(numTrianlges), + dPoints_(points), + dVectices_(vertices) + {} + + INLINE_FUNCTION_HD + triangleAccessor(const triangleAccessor&)= default; + + INLINE_FUNCTION_HD + triangleAccessor& operator=(const triangleAccessor&)= default; + + INLINE_FUNCTION_HD + triangleAccessor(triangleAccessor&&)= default; + + INLINE_FUNCTION_HD + triangleAccessor& operator=(triangleAccessor&&)= default; + + INLINE_FUNCTION_HD + ~triangleAccessor()=default; + + INLINE_FUNCTION_HD + realx3x3 triangle(uint32 i)const { + auto v = dVectices_[i]; + return realx3x3( + dPoints_[v.x_], + dPoints_[v.y_], + dPoints_[v.z_]); + } + + INLINE_FUNCTION_HD + realx3x3 operator()(uint32 i)const { return triangle(i); } + + INLINE_FUNCTION_HD + realx3x3 operator[](uint32 i)const { return triangle(i); } + + INLINE_FUNCTION_HD + uint32 numPoints()const { return numPoints_; } + + INLINE_FUNCTION_HD + uint32 numTrianlges()const { return numTriangles_;} +}; + + + class triSurface : public IOobject @@ -42,19 +107,29 @@ class triSurface private: /// points of triangles - realx3Field_D points_; + realx3Field_D points_{"points", "points"}; /// vectices indices of triangles - uint32x3Field_D vertices_; + uint32x3Field_D vertices_{"vertices", "vertices"}; /// area of each triangle - realField_D area_; + realField_D area_{"area", "area"}; /// normal vector of triangles - realx3Field_D normals_; + realx3Field_D normals_{"normals", "normals"}; protected: + + bool calculateNormals(); + bool calculateArea(); + + bool append(const realx3x3Field_H& triangles); + + bool append(const realx3x3Vector& triangles); + + + // to be used by multiTriSurface triSurface(const objectFile& obj, repository* owner); public: @@ -77,9 +152,6 @@ public: ~triSurface() override = default; - bool appendTriSurface(const realx3x3Field_H& triangles); - - bool appendTriSurface(const realx3x3Vector& triangles); //// - Methods uint32 numPoints() const @@ -144,6 +216,16 @@ public: area_.clear(); normals_.clear(); } + + /// Obtain an object for accessing triangles + auto getTriangleAccessor()const + { + return triangleAccessor( + points_.size(), + points_.deviceViewAll(), + vertices_.size(), + vertices_.deviceViewAll()); + } //// - IO operations diff --git a/src/phasicFlow/triSurface/triSurfaceKernels.hpp b/src/phasicFlow/triSurface/triSurfaceKernels.hpp index fbb1d7d6..7f79d06a 100644 --- a/src/phasicFlow/triSurface/triSurfaceKernels.hpp +++ b/src/phasicFlow/triSurface/triSurfaceKernels.hpp @@ -29,19 +29,47 @@ namespace pFlow::triSurfaceKernels { INLINE_FUNCTION_H -bool calculateArea(const realx3Field_D& points, const int32x3Field_D& vertices, realField_D& area) +bool calculateArea( + const realx3Field_D& points, + const uint32x3Field_D& vertices, + realField_D& area) { auto numTri = vertices.size(); - auto areaD = area.deviceViewAll(); - auto pointsD = points.deviceViewAll(); - auto verticesD = vertices.deviceViewAll(); + auto& areaD = area.deviceViewAll(); + auto& pointsD = points.deviceViewAll(); + auto& verticesD = vertices.deviceViewAll(); Kokkos::parallel_for( "pFlow::triSurfaceKernels::calculateArea", numTri, - LAMBDA_HD(int32 i){ + LAMBDA_HD(uint32 i){ auto v = verticesD[i]; - areaD[i] = pFlow::triangleFunctions::triangleSurface( + areaD[i] = pFlow::triangle::surface( + pointsD[v.x()], + pointsD[v.y()], + pointsD[v.z()]); + }); + + return true; +} + +INLINE_FUNCTION_H +bool calculateNormals( + const realx3Field_D& points, + const uint32x3Field_D& vertices, + realx3Field_D& normals) +{ + auto numTri = vertices.size(); + auto& normalsD = normals.deviceViewAll(); + auto& pointsD = points.deviceViewAll(); + auto& verticesD = vertices.deviceViewAll(); + + Kokkos::parallel_for( + "pFlow::triSurfaceKernels::calculateNormals", + numTri, + LAMBDA_HD(uint32 i){ + auto v = verticesD[i]; + normalsD[i] = pFlow::triangle::normal( pointsD[v.x()], pointsD[v.y()], pointsD[v.z()]);