From 89d7e1f0bae182bc464399f6160c29bb4b36f378 Mon Sep 17 00:00:00 2001 From: Hamidreza Norouzi Date: Sat, 13 Apr 2024 01:41:11 -0700 Subject: [PATCH] Particle insertion (domain check) - Domain check is added to prevent points that are outside of internalDomainBox --- .../InsertionRegion/InsertionRegion.cpp | 22 ++++++++++----- .../InsertionRegion/InsertionRegion.hpp | 1 + .../Insertion/insertion/insertion.cpp | 6 +++++ .../Insertion/insertion/insertion.hpp | 3 +++ .../insertionRegion/insertionRegion.cpp | 6 +++++ .../insertionRegion/insertionRegion.hpp | 5 +++- .../setFieldList/setFieldEntryTemplates.cpp | 2 +- .../boundaries/boundaryBase/boundaryBase.hpp | 19 +++++++++++++ .../boundaries/boundaryList.cpp | 27 ++++++++++++++----- .../boundaries/boundaryList.hpp | 12 +++++++++ src/phasicFlow/structuredData/box/box.hpp | 4 ++- .../structuredData/cylinder/cylinder.hpp | 4 ++- .../infinitePlane/infinitePlane.hpp | 15 +++++------ .../internalPoints/pointFlagKernels.hpp | 2 +- .../pointStructure/pointStructure.hpp | 6 +++++ .../structuredData/sphere/sphere.hpp | 2 +- 16 files changed, 109 insertions(+), 27 deletions(-) diff --git a/src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp b/src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp index dda67f80..18f7c23e 100644 --- a/src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp +++ b/src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp @@ -68,6 +68,16 @@ bool pFlow::InsertionRegion::insertParticles if(newNum == 0) return true; + auto internalBox = pStruct().internalDomainBox(); + + if( !(internalBox.minPoint() < internalBox.maxPoint())) + { + WARNING<<"Minimum point of internal point is not lower than "<< + "the maximum point \n"<< + "minimum point: "<< internalBox.minPoint()<< + "\nmaximum point:"<::insertParticles realVector diams("diams", newNum, 0, RESERVE()); - for(uint32 i=0; i::insertParticles { realx3 p = pReg.peek(); + + // check if point is inside internal box + if(!internalBox.isInside(p))continue; + if( !checkForContact(pos, diams, p, d) ) { names.push_back(name); diff --git a/src/Particles/Insertion/InsertionRegion/InsertionRegion.hpp b/src/Particles/Insertion/InsertionRegion/InsertionRegion.hpp index 7efb6191..a1bc383c 100644 --- a/src/Particles/Insertion/InsertionRegion/InsertionRegion.hpp +++ b/src/Particles/Insertion/InsertionRegion/InsertionRegion.hpp @@ -22,6 +22,7 @@ Licence: #define __InsertionRegion_hpp__ #include "dictionary.hpp" +#include "pointStructure.hpp" #include "insertionRegion.hpp" namespace pFlow diff --git a/src/Particles/Insertion/insertion/insertion.cpp b/src/Particles/Insertion/insertion/insertion.cpp index e2f909c4..09f29f98 100644 --- a/src/Particles/Insertion/insertion/insertion.cpp +++ b/src/Particles/Insertion/insertion/insertion.cpp @@ -39,6 +39,12 @@ pFlow::insertion::insertion(particles& prtcl) readInsertionDict(*this); } +const pFlow::pointStructure& +pFlow::insertion::pStruct() const +{ + return particles_.pStruct(); +} + bool pFlow::insertion::read(iIstream& is, const IOPattern& iop) { diff --git a/src/Particles/Insertion/insertion/insertion.hpp b/src/Particles/Insertion/insertion/insertion.hpp index c70ad77a..92d91f3e 100644 --- a/src/Particles/Insertion/insertion/insertion.hpp +++ b/src/Particles/Insertion/insertion/insertion.hpp @@ -27,6 +27,7 @@ namespace pFlow // forward class particles; +class pointStructure; /** * Base class for particle insertion @@ -88,6 +89,8 @@ public: return particles_; } + const pointStructure& pStruct()const; + inline bool readFromFile() const { return readFromFile_; diff --git a/src/Particles/Insertion/insertionRegion/insertionRegion.cpp b/src/Particles/Insertion/insertionRegion/insertionRegion.cpp index 437b3009..364dca8f 100644 --- a/src/Particles/Insertion/insertionRegion/insertionRegion.cpp +++ b/src/Particles/Insertion/insertionRegion/insertionRegion.cpp @@ -158,6 +158,12 @@ pFlow::insertionRegion::insertionRegion( } } +const pFlow::pointStructure& +pFlow::insertionRegion::pStruct() const +{ + return Insertion().pStruct(); +} + pFlow::uint32 pFlow::insertionRegion::numberToBeInserted(uint32 iter, real t, real dt) { diff --git a/src/Particles/Insertion/insertionRegion/insertionRegion.hpp b/src/Particles/Insertion/insertionRegion/insertionRegion.hpp index 2e9c3eb1..89b57884 100644 --- a/src/Particles/Insertion/insertionRegion/insertionRegion.hpp +++ b/src/Particles/Insertion/insertionRegion/insertionRegion.hpp @@ -31,6 +31,7 @@ namespace pFlow class dictionary; class insertion; +class pointStructure; /** * This class defines all the necessary enteties for defining an insertion @@ -92,7 +93,7 @@ private: /// insertion region dictionary const dictionary& dict_; - /// ref to pointStructure + /// ref to insertion const insertion& insertion_; /// @brief time control for insertion events @@ -163,6 +164,8 @@ public: return insertion_; } + const pointStructure& pStruct()const; + inline bool insertionTime(uint32 iter, real t, real dt) const { return tControl_.timeEvent(iter, t, dt); diff --git a/src/phasicFlow/setFieldList/setFieldEntryTemplates.cpp b/src/phasicFlow/setFieldList/setFieldEntryTemplates.cpp index 45489eca..2ee5576c 100644 --- a/src/phasicFlow/setFieldList/setFieldEntryTemplates.cpp +++ b/src/phasicFlow/setFieldList/setFieldEntryTemplates.cpp @@ -23,7 +23,7 @@ template bool pFlow::setFieldEntry::checkForType()const { word typeName( entry_.firstPart() ); - return basicTypeName() == typeName; + return getTypeName() == typeName; }; template diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp index 06766b97..15b8284f 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp @@ -153,12 +153,30 @@ public: (dict, bplane, internal, bndrs, thisIndex) ); + /// The length from boundary plane into the domain + /// where beyond that distance internal points exist. + /// By conventions is it always equal to neighborLength_ + real neighborLengthIntoInternal()const + { + return neighborLength_; + } + + /// The distance length from boundary plane + /// where neighbor particles exist in that distance. + /// This length may be modified in each boundary type + /// as required. In this case the boundaryExtensionLength + /// method should also be modified accordingly. virtual real neighborLength()const { return neighborLength_; } + /// The extention length (in vector form) for the boundary + /// as required by each boundary type. It is allowed for + /// each boundary type to be extended outward to allow + /// particles to stay more in its list before being removed + /// from its list. virtual realx3 boundaryExtensionLength()const { @@ -240,6 +258,7 @@ public: const boundaryBase& mirrorBoundary()const; + /// the actual boundary plane of this boundary const plane& boundaryPlane()const { return boundaryPlane_; diff --git a/src/phasicFlow/structuredData/boundaries/boundaryList.cpp b/src/phasicFlow/structuredData/boundaries/boundaryList.cpp index 81a227a5..9073d9e7 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryList.cpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryList.cpp @@ -143,12 +143,27 @@ bool pFlow::boundaryList::setLists() return true; } -bool pFlow::boundaryList::beforeIteration -( - uint32 iter, - real t, - real dt -) +pFlow::box +pFlow::boundaryList::internalDomainBox() const +{ + const auto& thisBox = pStruct_.thisDomain().domainBox(); + + const realx3 lowerPointDisplacement = { + boundary(0).neighborLengthIntoInternal(), + boundary(2).neighborLengthIntoInternal(), + boundary(4).neighborLengthIntoInternal()}; + + const realx3 upperPointDisplacement = { + boundary(1).neighborLengthIntoInternal(), + boundary(3).neighborLengthIntoInternal(), + boundary(5).neighborLengthIntoInternal()}; + + return {thisBox.minPoint() + lowerPointDisplacement, + thisBox.maxPoint() - upperPointDisplacement}; +} + +bool +pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt) { // it is time to update lists if(timeControl_.timeEvent(iter, t, dt)) diff --git a/src/phasicFlow/structuredData/boundaries/boundaryList.hpp b/src/phasicFlow/structuredData/boundaries/boundaryList.hpp index 6581d6cd..c5e3b25e 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryList.hpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryList.hpp @@ -46,6 +46,8 @@ private: domain extendedDomain_; + box internalDomainBox_; + bool listSet_ = false; void setExtendedDomain(); @@ -105,6 +107,16 @@ public: return extendedDomain_; } + inline + const auto& extendedDomainBox()const + { + return extendedDomain_.domainBox(); + } + + box internalDomainBox()const; + + + bool beforeIteration(uint32 iter, real t, real dt); bool iterate(uint32 iter, real t, real dt); diff --git a/src/phasicFlow/structuredData/box/box.hpp b/src/phasicFlow/structuredData/box/box.hpp index 8a350405..67d0e8dc 100644 --- a/src/phasicFlow/structuredData/box/box.hpp +++ b/src/phasicFlow/structuredData/box/box.hpp @@ -79,10 +79,12 @@ public: //// - Methods + INLINE_FUNCTION_HD bool isInside(const realx3& point)const { - return point > min_ && point = min_ && point <=max_; } INLINE_FUNCTION_HD diff --git a/src/phasicFlow/structuredData/cylinder/cylinder.hpp b/src/phasicFlow/structuredData/cylinder/cylinder.hpp index cc3617ad..ad54b4c1 100644 --- a/src/phasicFlow/structuredData/cylinder/cylinder.hpp +++ b/src/phasicFlow/structuredData/cylinder/cylinder.hpp @@ -90,7 +90,9 @@ public: auto p1Point = point-p1_; auto H = cross(p1Point , axisVector_); auto H2 = dot(H,H); - if( H2 < radius2_*axisVector2_) + + // the shell itslef is considered as inside point + if( H2 <= radius2_*axisVector2_) { real t = dot(p1Point, axisVector_)/axisVector2_; if(t >= 0.0 && t <= 1.0) diff --git a/src/phasicFlow/structuredData/infinitePlane/infinitePlane.hpp b/src/phasicFlow/structuredData/infinitePlane/infinitePlane.hpp index 28c548bc..2b8440df 100644 --- a/src/phasicFlow/structuredData/infinitePlane/infinitePlane.hpp +++ b/src/phasicFlow/structuredData/infinitePlane/infinitePlane.hpp @@ -103,6 +103,12 @@ public: return pointFromPlane(p)>=0; } + INLINE_FUNCTION_HD + bool pointInNegativeSide(const realx3& p)const + { + return pointFromPlane(p)<0; + } + INLINE_FUNCTION_HD bool inPositiveDistance(const realx3& p, real dist)const { @@ -114,14 +120,7 @@ public: bool inNegativeDistance(const realx3& p, real dist)const { real d = pointFromPlane(p); - return d < 0.0 && d <= -dist; - } - - - INLINE_FUNCTION_HD - bool pointInNegativeSide(const realx3& p)const - { - return pointFromPlane(p)<0; + return d < 0.0 && d >= -dist; } INLINE_FUNCTION_HD diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp index 3c5c1540..91764cf6 100644 --- a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp +++ b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp @@ -543,7 +543,7 @@ pFlow::uint32 pFlow::pointFlag::changeCapacity while( emptySpots < reqEmptySpots ) { newCap = newCap*pFlow::gSettings::vectorGrowthFactor__+1; - uint32 emptySpots = newCap - numActive_; + emptySpots = newCap - numActive_; } viewType newFlags(flags_.label(), newCap); diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.hpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.hpp index 27679178..065e58cd 100644 --- a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.hpp +++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.hpp @@ -153,6 +153,12 @@ public: { return boundaries_.extendedDomain(); } + + inline + auto internalDomainBox()const + { + return boundaries_.internalDomainBox(); + } // - IO methods diff --git a/src/phasicFlow/structuredData/sphere/sphere.hpp b/src/phasicFlow/structuredData/sphere/sphere.hpp index 1a8236b1..21b20d1c 100644 --- a/src/phasicFlow/structuredData/sphere/sphere.hpp +++ b/src/phasicFlow/structuredData/sphere/sphere.hpp @@ -76,7 +76,7 @@ public: { auto cPoint = point-center_; auto dist2 = dot(cPoint,cPoint); - return dist2 < radius2_; + return dist2 <= radius2_; }