Particle insertion (domain check)

- Domain check is added to prevent points that are outside of internalDomainBox
This commit is contained in:
Hamidreza Norouzi 2024-04-13 01:41:11 -07:00
parent 9c2a9a81b0
commit 89d7e1f0ba
16 changed files with 109 additions and 27 deletions

View File

@ -68,6 +68,16 @@ bool pFlow::InsertionRegion<ShapeType>::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:"<<internalBox.maxPoint()<<END_WARNING;
}
names.reserve(max(newNum,names.capacity()));
pos.reserve(max(newNum,pos.capacity()));
names.clear();
@ -75,12 +85,6 @@ bool pFlow::InsertionRegion<ShapeType>::insertParticles
realVector diams("diams", newNum, 0, RESERVE());
for(uint32 i=0; i<newNum; i++)
{
uint32 idx;
shapes_.shapeNameToIndex(names[i], idx);
diams[i] = shapes_.boundingDiameter(idx);
}
size_t n = 0;
uint32 idx;
@ -94,6 +98,10 @@ bool pFlow::InsertionRegion<ShapeType>::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);

View File

@ -22,6 +22,7 @@ Licence:
#define __InsertionRegion_hpp__
#include "dictionary.hpp"
#include "pointStructure.hpp"
#include "insertionRegion.hpp"
namespace pFlow

View File

@ -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)
{

View File

@ -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_;

View File

@ -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)
{

View File

@ -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);

View File

@ -23,7 +23,7 @@ template <typename Type>
bool pFlow::setFieldEntry::checkForType()const
{
word typeName( entry_.firstPart() );
return basicTypeName<Type>() == typeName;
return getTypeName<Type>() == typeName;
};
template <typename Type>

View File

@ -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_;

View File

@ -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))

View File

@ -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);

View File

@ -79,10 +79,12 @@ public:
//// - Methods
INLINE_FUNCTION_HD
bool isInside(const realx3& point)const
{
return point > min_ && point <max_;
// box planes are also included as the inside of the box
return point >= min_ && point <=max_;
}
INLINE_FUNCTION_HD

View File

@ -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)

View File

@ -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

View File

@ -543,7 +543,7 @@ pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::changeCapacity
while( emptySpots < reqEmptySpots )
{
newCap = newCap*pFlow::gSettings::vectorGrowthFactor__+1;
uint32 emptySpots = newCap - numActive_;
emptySpots = newCap - numActive_;
}
viewType newFlags(flags_.label(), newCap);

View File

@ -154,6 +154,12 @@ public:
return boundaries_.extendedDomain();
}
inline
auto internalDomainBox()const
{
return boundaries_.internalDomainBox();
}
// - IO methods
/// @brief read the point structure from the input stream.

View File

@ -76,7 +76,7 @@ public:
{
auto cPoint = point-center_;
auto dist2 = dot(cPoint,cPoint);
return dist2 < radius2_;
return dist2 <= radius2_;
}