Collision check and particlePosition -> random

- A new class is added for simple collision check
- position particles in utility is upgraded
- morton sorting is not active yet for particlesPhasicFlow
This commit is contained in:
Hamidreza Norouzi 2024-04-13 07:07:36 -07:00
parent 89d7e1f0ba
commit e395c379cb
15 changed files with 353 additions and 399 deletions

View File

@ -9,6 +9,7 @@ particles/regularParticleIdHandler/regularParticleIdHandler.cpp
SphereParticles/sphereShape/sphereShape.cpp
SphereParticles/sphereParticles/sphereParticles.cpp
SphereParticles/sphereParticles/sphereParticlesKernels.cpp
Insertion/collisionCheck/collisionCheck.cpp
Insertion/insertionRegion/insertionRegion.cpp
Insertion/insertion/insertion.cpp
Insertion/shapeMixture/shapeMixture.cpp

View File

@ -0,0 +1,91 @@
#include "collisionCheck.hpp"
bool
pFlow::collisionCheck::build()
{
fill(head_, static_cast<uint32>(-1));
fill(next_, static_cast<uint32>(-1));
for (auto i = 0uL; i < position_.size(); i++)
{
if(!searchBox_.isInside(position_[i]))
{
fatalErrorInFunction<<
"Point "<< position_[i]<< "is not in search box"<<endl;
}
const auto ind = pointIndex(position_[i]);
next_[i] = head_(ind.x(), ind.y(), ind.z());
head_(ind.x(), ind.y(), ind.z()) = static_cast<uint32>(i);
}
return true;
}
pFlow::collisionCheck::collisionCheck(
box sBox,
real dx,
const realx3Vector& pos,
const realVector& diam
)
: searchBox_(sBox),
dx_(dx),
nCells_((sBox.maxPoint() - sBox.minPoint()) / dx + realx3(1.0)),
position_(pos),
diameters_(diam),
next_("next", pos.size()),
head_("head", nCells_.x(), nCells_.y(), nCells_.z())
{
if(!build())
{
fatalExit;
}
}
bool
pFlow::collisionCheck::checkPoint(const realx3& p, const real d) const
{
if(!searchBox_.isInside(p)) return false;
const auto ind = pointIndex(p);
uint32 n = head_(ind.x(), ind.y(), ind.z());
while( n != -1)
{
if( (position_[n]-p).length() - 0.5*(diameters_[n]+d )<0.0 )
{
return false;
}
n = next_[n];
}
return true;
}
bool
pFlow::collisionCheck::mapLastAddedParticle()
{
size_t n = position_.size()-1;
if( next_.size() != n )
{
fatalErrorInFunction<<
"size mismatch of next and position"<<endl;
return false;
}
next_.push_back(-1);
const auto& p = position_[n];
if(!searchBox_.isInside(p))
{
fatalErrorInFunction<<
"Point "<< p <<" is not inside the search box"<<endl;
return false;
}
const auto ind = pointIndex(p);
next_[n] = head_(ind.x(), ind.y(), ind.z());
head_(ind.x(), ind.y(), ind.z()) = static_cast<uint32>(n);
return true;
}

View File

@ -0,0 +1,54 @@
#ifndef __collisionCheck_hpp__
#define __collisionCheck_hpp__
#include "Vectors.hpp"
#include "VectorSingles.hpp"
#include "box.hpp"
namespace pFlow
{
class collisionCheck
{
private:
box searchBox_;
real dx_;
int32x3 nCells_;
const realx3Vector& position_;
const realVector& diameters_;
uint32Vector next_;
ViewType3D<uint32, HostSpace> head_;
int32x3 pointIndex(const realx3& p)const
{
return int32x3( (p - searchBox_.minPoint())/dx_ );
}
bool build();
public:
collisionCheck(
box sBox,
real dx,
const realx3Vector& pos,
const realVector& diam
);
bool checkPoint(const realx3& p, const real d)const;
bool mapLastAddedParticle();
};
}
#endif //__collisionCheck_hpp__

View File

@ -65,9 +65,19 @@ public:
//// - Methods
virtual bool isInside(const realx3& point) const override;
bool isInside(const realx3& point) const override;
virtual realx3 peek() const override;
realx3 peek() const override;
const realx3& minPoint()const override
{
return region_.minPoint();
}
const realx3& maxPoint()const override
{
return region_.maxPoint();
}
//// - IO operatoins

View File

@ -62,6 +62,16 @@ public:
realx3 peek() const;
const auto& minPoint()const
{
return minPoint_;
}
const auto& maxPoint()const
{
return maxPoint_;
}
//// IO operation
bool read(const dictionary& dict);

View File

@ -57,6 +57,10 @@ public:
virtual realx3 peek() const = 0;
virtual const realx3& minPoint()const = 0;
virtual const realx3& maxPoint()const = 0;
//// - IO operatoins
virtual bool read(const dictionary& dict) = 0;

View File

@ -3,7 +3,7 @@ set(source_files
particlesPhasicFlow.cpp
positionParticles/positionParticles.cpp
positionOrdered/positionOrdered.cpp
#positionRandom/positionRandom.cpp
positionRandom/positionRandom.cpp
empty/empty.cpp
)
#set(link_lib phasicFlow Kokkos::kokkos Interaction Utilities)

View File

@ -30,7 +30,7 @@ pFlow::empty::empty(
positionParticles(control, dict),
position_
(
"empty",maxNumberOfParticles_, 0, RESERVE()
"empty",maxNumberOfParticles(), 0, RESERVE()
)
{
}

View File

@ -22,7 +22,6 @@ Licence:
#define __empty_hpp__
#include "positionParticles.hpp"
#include "box.hpp"
namespace pFlow
{
@ -32,10 +31,7 @@ class empty
:
public positionParticles
{
protected:
dictionary emptyDict_;
private:
realx3Vector position_;
@ -54,33 +50,33 @@ public:
empty,
dictionary);
virtual ~empty() = default;
~empty() final = default;
//// - Methods
virtual uint64 numPoints()const
uint32 numPoints()const final
{
return 0;
}
virtual uint64 size()const
uint32 size()const final
{
return 0;
}
real maxDiameter() const override
real maxDiameter() const final
{
return 1.0;
}
// - const access to position
virtual const realx3Vector& position()const
const realx3Vector& position()const final
{
return position_;
}
// - access to position
virtual realx3Vector& position()
realx3Vector& position() final
{
return position_;
}

View File

@ -18,11 +18,9 @@ Licence:
-----------------------------------------------------------------------------*/
#include "positionOrdered.hpp"
#include "error.hpp"
#include "streams.hpp"
#include "dictionary.hpp"
#include "positionOrdered.hpp"
bool pFlow::positionOrdered::findAxisIndex()
@ -45,9 +43,9 @@ bool pFlow::positionOrdered::findAxisIndex()
return false;
}
realx3 uV[3];
std::array<realx3,3> uV;
size_t i=0;
for(auto& ca: axisOrder_)
for(const auto& ca: axisOrder_)
{
if(ca == "x")
{
@ -82,15 +80,16 @@ bool pFlow::positionOrdered::positionPointsOrdered()
position_.clear();
realx3 dl(diameter_);
auto minP = region_->minPoint();
auto maxP = region_->maxPoint();
const auto& region = pRegion();
auto minP = region.minPoint();
auto maxP = region.maxPoint();
auto cntr = minP;
size_t n = 0;
while( n < numPoints_ )
{
if(region_->isInside(cntr))
if(region.isInside(cntr))
{
position_.push_back(cntr);
n++;
@ -130,7 +129,7 @@ pFlow::positionOrdered::positionOrdered
positionParticles(control, dict),
poDict_
(
dict.subDict("positionOrderedInfo")
dict.subDict("orderedInfo")
),
diameter_
(
@ -146,7 +145,10 @@ pFlow::positionOrdered::positionOrdered
),
position_
(
"positionOrdered", maxNumberOfParticles_, numPoints_ ,RESERVE()
"positionOrdered",
max(maxNumberOfParticles(), numPoints_),
numPoints_ ,
RESERVE()
)
{
@ -155,13 +157,6 @@ pFlow::positionOrdered::positionOrdered
fatalExit;
}
if(!region_)
{
fatalErrorInFunction<<"You must provided a region (box, cylinder, ...) for positioning particles in dictionary "<<
dict.globalName()<<endl;
fatalExit;
}
if(!positionPointsOrdered())
{
fatalExit;

View File

@ -31,13 +31,13 @@ class positionOrdered
:
public positionParticles
{
protected:
private:
dictionary poDict_;
real diameter_;
size_t numPoints_;
uint32 numPoints_;
wordList axisOrder_;
@ -60,7 +60,7 @@ protected:
public:
// - type Info
TypeInfo("positionOrdered");
TypeInfo("ordered");
positionOrdered(
systemControl& control,
@ -72,33 +72,33 @@ public:
positionOrdered,
dictionary);
virtual ~positionOrdered() = default;
~positionOrdered() final = default;
//// - Methods
virtual uint64 numPoints()const
uint32 numPoints()const final
{
return position_.size();
return static_cast<uint32>(position_.size());
}
virtual uint64 size()const
uint32 size()const final
{
return position_.size();
return static_cast<uint32>(position_.size());
}
real maxDiameter() const override
real maxDiameter() const final
{
return diameter_;
}
// - const access to position
virtual const realx3Vector& position()const
const realx3Vector& position()const final
{
return position_;
}
// - access to position
virtual realx3Vector& position()
realx3Vector& position() final
{
return position_;
}

View File

@ -20,16 +20,11 @@ Licence:
#include "positionParticles.hpp"
#include "vocabs.hpp"
#include "box.hpp"
#include "cylinder.hpp"
#include "sphere.hpp"
#include "cells.hpp"
//#include "contactSearchFunctions.hpp"
#include "dictionary.hpp"
#include "systemControl.hpp"
#include "streams.hpp"
pFlow::realx3Vector pFlow::positionParticles::sortByMortonCode(realx3Vector& position)const
pFlow::realx3Vector pFlow::positionParticles::sortByMortonCode(
const realx3Vector& position)const
{
struct indexMorton
{
@ -81,25 +76,20 @@ pFlow::positionParticles::positionParticles
systemControl& control,
const dictionary& dict
)
:
regionType_(dict.getValOrSet<word>("regionType", "domain")),
maxNumberOfParticles_(dict.getValOrSet(
"maxNumberOfParticles",
static_cast<uint32>(10000))),
mortonSorting_(dict.getValOrSet("mortonSorting", Logical("Yes")))
{
maxNumberOfParticles_ = dict.getValOrSet("maxNumberOfParticles", static_cast<uint64>(10000));
mortonSorting_ = dict.getValOrSet("mortonSorting", Logical("Yes"));
if( dict.containsDictionay("box") )
if( regionType_ != "domain" )
{
region_ = makeUnique<region<box>>(dict.subDict("box"));
}
else if(dict.containsDictionay("cylinder"))
{
WARNING<<"cylinder region is not set!"<<END_WARNING;
//region_ = makeUnique<region<cylinder>>(dict.subDict("cylinder"));
}
else if(dict.containsDictionay("sphere"))
{
WARNING<<"sphere region is not set!"<<END_WARNING;
//region_ = makeUnique<region<sphere>>(dict.subDict("sphere"));
}
pRegion_ = peakableRegion::create(
regionType_,
dict.subDict(regionType_+"Info"));
}
else
{
fileDictionary domainDict
@ -113,7 +103,7 @@ pFlow::positionParticles::positionParticles
},
&control.settings()
);
region_ = makeUnique<region<box>>( domainDict.subDict("globalBox"));
pRegion_ = peakableRegion::create(regionType_,domainDict.subDict("globalBox"));
}
}
@ -128,9 +118,8 @@ pFlow::realx3Vector pFlow::positionParticles::getFinalPosition()
else
{
realx3Vector vec("final position",position().capacity(), 0 , RESERVE());
vec.assign( position().begin(), position().end());
return std::move(vec);
vec.assign( position().begin(), position().end());
return vec;
}
}

View File

@ -23,96 +23,39 @@ Licence:
#include "virtualConstructor.hpp"
#include "Vectors.hpp"
#include "dictionary.hpp"
#include "systemControl.hpp"
#include "peakableRegion.hpp"
namespace pFlow
{
class regionBase
{
public:
regionBase() = default;
regionBase(const regionBase&) = default;
class dictionary;
class systemControl;
regionBase& operator =(const regionBase&) = default;
virtual ~regionBase() = default;
virtual bool isInside(const realx3 point)const = 0;
virtual realx3 minPoint()const =0;
virtual realx3 maxPoint()const =0;
virtual word name()const =0;
};
template<typename T>
class region
:
public regionBase
{
protected:
T region_;
public:
region(const T& rgn)
:
region_(rgn)
{}
region(const dictionary& dict)
:
region_(dict)
{}
region(const region&) = default;
region& operator =(const region&) = default;
virtual ~region()=default;
bool isInside(const realx3 point) const override
{
return region_.isInside(point);
}
realx3 minPoint()const override
{
return region_.minPoint();
}
realx3 maxPoint()const override
{
return region_.maxPoint();
}
word name()const override
{
return region_.typeName();
}
};
class positionParticles
{
protected:
private:
uniquePtr<regionBase> region_ = nullptr;
uniquePtr<peakableRegion> pRegion_ = nullptr;
size_t maxNumberOfParticles_ = 10000;
word regionType_;
uint32 maxNumberOfParticles_ = 10000;
Logical mortonSorting_;
static const size_t numReports_ = 40;
realx3Vector sortByMortonCode(realx3Vector& position)const;
realx3Vector sortByMortonCode(const realx3Vector& position)const;
protected:
static const uint32 numReports_ = 40;
const auto& pRegion()const
{
return pRegion_();
}
public:
@ -133,11 +76,22 @@ public:
virtual ~positionParticles() = default;
//// - Methods
//// - Methods
virtual uint64 numPoints()const = 0;
bool mortonSorting()const
{
return mortonSorting_();
}
virtual uint64 size()const = 0;
inline
auto maxNumberOfParticles()const
{
return maxNumberOfParticles_;
}
virtual uint32 numPoints()const = 0;
virtual uint32 size()const = 0;
virtual real maxDiameter() const = 0;

View File

@ -20,141 +20,81 @@ Licence:
#include "positionRandom.hpp"
#include "uniformRandomReal.hpp"
#include "NBSLevel0.hpp"
#include "unsortedPairs.hpp"
#include "box.hpp"
#include "collisionCheck.hpp"
namespace pFlow
bool pFlow::positionRandom::positionOnePass(collisionCheck& collCheck)
{
using SearchType = NBSLevel0<DefaultExecutionSpace> ;
using ContainerType = unsortedPairs<DefaultExecutionSpace, int32>;
auto const& region = pRegion();
auto n = static_cast<uint32>(position_.size());
int32 findCollisions(
ContainerType& pairs,
int32Vector_HD& flags);
int32 findCollisions(int32 num, realx3Vector_HD& points, real diam)
{
int32 res =0;
for(auto i=0; i<num;i++)
for(uint32 iter = 0; iter<numPoints_; iter++)
{
for(auto j=i+1; j<num; j++)
realx3 p = region.peek();
if( collCheck.checkPoint(p, diameter_) )
{
if(sphereSphereCheck(points[i],points[j],diam,diam))res++;
}
}
return res;
}
}
bool pFlow::positionRandom::positionOnePass(int32 pass, int32 startNum)
{
realVector_D diameter(startNum , diameter_);
int32Vector_HD flagHD(startNum, 0);
realx3Vector_HD positionHD(startNum);
auto minP = region_->minPoint();
auto maxP = region_->maxPoint();
SearchType search(
box(minP, maxP),
diameter_,
positionHD.deviceViewAll(),
diameter.deviceViewAll());
ContainerType pairs(3*startNum);
REPORT(1)<< "Positioning "<<
greenText("(Pass #"<< pass+1<<")")<<
": started with "<< startNum <<" points."<<endREPORT;
fillPoints(startNum, positionHD, flagHD);
search.broadSearch(pairs, range(0, startNum));
int32 numCollisions = findCollisions(pairs, flagHD);
REPORT(2)<< "Positioned " << cyanText(startNum - numCollisions) <<
" without collision \n"<<endREPORT;
if(startNum-numCollisions >= numPoints_ )
{
REPORT(1)<<"Selected "<< cyanText(numPoints_)<< " for the final field.\n"<<endREPORT;
positionHD.syncViews();
position_.clear();
int32 n=0;
for(int32 i=0; i<startNum; i++)
{
if(flagHD[i] == 0 )
{
position_.push_back( positionHD[i]);
n++;
if(n==numPoints_)break;
}
position_.push_back(p);
diameters_.push_back(diameter_);
if(!collCheck.mapLastAddedParticle())
{
fatalErrorInFunction;
return false;
}
n++;
if(n == numPoints_) break;
}
return true;
}
return false;
return true;
}
bool pFlow::positionRandom::positionPointsRandom()
{
position_.clear();
diameters_.clear();
if(numPoints_ == 0)return true;
size_t pass = 0;
int32 startNum = numPoints_;
uint32 pass = 0;
collisionCheck collCheck(
box(pRegion().minPoint(), pRegion().maxPoint()),
diameter_,
position_,
diameters_);
while ( pass <maxIterations_)
{
if( positionOnePass(pass, startNum) )return true;
startNum = 1.1*startNum+1;
if( !positionOnePass(collCheck) )
{
return false;
}
pass++;
REPORT(1)<<"Positioning "<< Green_Text("(Pass #"<< pass<<")")<<
": number of non-colliding spheres is "<<
Yellow_Text(position_.size())<<END_REPORT;
if( position_.size() == numPoints_ )
{
REPORT(0)<<END_REPORT;
return true;
}
}
fatalErrorInFunction<<
" cannot position "<< numPoints_ << " in the domain in " << maxIterations_ << " iterations.\n" <<
" Cannot position "<< numPoints_ << " in the domain in " <<
maxIterations_ << " iterations.\n" <<
" you may increase maxIterations for positioning points.\n";
return false;
}
bool pFlow::positionRandom::inCollision
(
const realx3 &cntr,
real diam
)
{
for(const auto& cp: position_)
{
if( length(cp-cntr) <= diam ) return true;
}
return false;
}
pFlow::positionRandom::positionRandom
(
@ -165,7 +105,7 @@ pFlow::positionRandom::positionRandom
positionParticles(control, dict),
prDict_
(
dict.subDict("positionRandomInfo")
dict.subDict("randomInfo")
),
diameter_
(
@ -173,107 +113,35 @@ pFlow::positionRandom::positionRandom
),
numPoints_
(
prDict_.getVal<size_t>("numPoints")
prDict_.getVal<uint32>("numPoints")
),
maxIterations_
(
prDict_.getValOrSet("maxIterations", 10)
prDict_.getValOrSet("maxIterations", 10u)
),
position_
(
maxNumberOfParticles_, RESERVE()
"position",
maxNumberOfParticles(),
0,
RESERVE()
),
diameters_
(
"diameters",
maxNumberOfParticles(),
0,
RESERVE()
)
{
reportInterval_ = max(numPoints_/numReports_, static_cast<size_t>(2));
reportInterval_ = max(numPoints_/numReports_, static_cast<uint32>(2));
if( !positionPointsRandom() )
{
fatalExit;
}
if(!region_)
{
fatalErrorInFunction<<"You must provided a region (box, cylinder, ...) for positioning particles in dictionary "<<
dict.globalName()<<endl;
fatalExit;
}
}
void pFlow::positionRandom::fillPoints(
uint numPoints,
realx3Vector_HD& points,
int32Vector_HD& flags )
{
uniformRandomReal rand;
auto minP = region_().minPoint();
auto maxP = region_().maxPoint();
for(size_t i=0; i<numPoints; i++)
{
if(flags[i] == 0)
{
bool loop=true;
size_t n=0;
while (loop)
{
auto pos = rand(minP, maxP);
if( region_().isInside(pos))
{
points[i] =pos;
loop = false;
}
n++;
if(n>100)
{
fatalErrorInFunction<<
"could not find a point inside region"<<region_->name()<<endl;
fatalExit;
}
}
}
}
points.modifyOnHost();
points.syncViews();
}
pFlow::int32 pFlow::findCollisions(
ContainerType& pairs,
int32Vector_HD& flags)
{
auto allPairs = pairs.getPairs();
auto num = pairs.capacity();
auto dFlags = flags.deviceView();
int32 numCollisions = 0;
Kokkos::parallel_reduce(
"positionRandom::findCollisions",
num,
LAMBDA_HD(int32 i, int32& valueToUpdate){
if(allPairs.isValid(i))
{
auto pair = allPairs.getPair(i);
if( dFlags[pair.first] ==0 )
{
dFlags[pair.first] = 1;
valueToUpdate++;
}
}
}, numCollisions);
flags.modifyOnDevice();
flags.syncViews();
return numCollisions;
}

View File

@ -2,90 +2,79 @@
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __positionOrdered_hpp__
#define __positionOrdered_hpp__
#ifndef __positionRandom_hpp__
#define __positionRandom_hpp__
#include "positionParticles.hpp"
#include "VectorSingles.hpp"
#include "VectorDuals.hpp"
#include "dictionary.hpp"
namespace pFlow
{
class collisionCheck;
class positionRandom
:
public positionParticles
class positionRandom : public positionParticles
{
protected:
private:
dictionary prDict_;
dictionary prDict_;
real diameter_;
real diameter_;
size_t numPoints_;
uint32 numPoints_;
size_t maxIterations_;
realx3Vector position_;
uint32 maxIterations_;
size_t reportInterval_;
realx3Vector position_;
bool positionOnePass(int32 pass, int32 startNum);
bool positionPointsRandom();
realVector diameters_;
bool inCollision(const realx3 &cntr, real diam);
uint32 reportInterval_;
void fillPoints(
uint numPoints,
realx3Vector_HD& points,
int32Vector_HD& flags );
bool positionOnePass(collisionCheck& collCheck);
bool positionPointsRandom();
public:
// - type Info
TypeInfo("positionRandom");
TypeInfo("random");
positionRandom(
systemControl& control,
const dictionary& dict);
positionRandom(systemControl& control, const dictionary& dict);
// - add this class to vCtor selection table
add_vCtor(
positionParticles,
positionRandom,
dictionary);
// - add this class to vCtor selection table
add_vCtor
(
positionParticles,
positionRandom,
dictionary
);
virtual ~positionRandom() = default;
~positionRandom() final = default;
//// - Methods
//// - Methods
virtual uint64 numPoints()const
uint32 numPoints() const final
{
return position_.size();
}
virtual uint64 size()const
uint32 size() const final
{
return position_.size();
}
@ -96,25 +85,18 @@ public:
}
// - const access to position
virtual const realx3Vector& position()const
const realx3Vector& position() const final
{
return position_;
}
// - access to position
virtual realx3Vector& position()
// - access to position
realx3Vector& position() final
{
return position_;
}
};
};
}
#endif // __positionOrdered_hpp__