557 lines
12 KiB
C++
557 lines
12 KiB
C++
/*------------------------------- phasicFlow ---------------------------------
|
|
O C enter of
|
|
O O E ngineering and
|
|
O O M ultiscale modeling of
|
|
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
|
|
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
|
|
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
|
|
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
|
|
|
-----------------------------------------------------------------------------*/
|
|
|
|
|
|
#include "pointStructure.hpp"
|
|
#include "pointStructureKernels.hpp"
|
|
#include "setFieldList.hpp"
|
|
#include "error.hpp"
|
|
#include "iOstream.hpp"
|
|
//#include "Time.hpp"
|
|
#include "mortonIndexing.hpp"
|
|
|
|
FUNCTION_H
|
|
bool pFlow::pointStructure::evaluatePointStructure()
|
|
{
|
|
if(pointFlag_.size() != pointPosition_.size())
|
|
{
|
|
fatalErrorInFunction<<
|
|
"number of elements in pointFlag and pointPosition is not equal \n";
|
|
return false;
|
|
}
|
|
|
|
setNumMaxPoints();
|
|
|
|
int32 minActive, maxActive;
|
|
numActivePoints_ = pFlow::pointStructureKernels::scanPointFlag(
|
|
0,
|
|
numPoints_,
|
|
static_cast<int8>(pointStructure::ACTIVE),
|
|
pointFlag_.deviceVectorAll(),
|
|
minActive,
|
|
maxActive
|
|
);
|
|
|
|
activeRange_ = {minActive, maxActive};
|
|
|
|
return true;
|
|
}
|
|
|
|
FUNCTION_H
|
|
void pFlow::pointStructure::setNumMaxPoints()
|
|
{
|
|
maxPoints_ = pointFlag_.capacity();
|
|
numPoints_ = pointFlag_.size();
|
|
}
|
|
|
|
FUNCTION_H
|
|
pFlow::realx3Field_D& pFlow::pointStructure::pointPosition()
|
|
{
|
|
return pointPosition_;
|
|
}
|
|
|
|
FUNCTION_H
|
|
pFlow::int8Field_HD& pFlow::pointStructure::pointFlag()
|
|
{
|
|
return pointFlag_;
|
|
}
|
|
|
|
pFlow::uniquePtr<pFlow::int32IndexContainer>
|
|
pFlow::pointStructure::getNewPointsIndices(int32 numNewPoints)const
|
|
{
|
|
|
|
|
|
if( capacity() - activeRange_.second >= numNewPoints )
|
|
{
|
|
// fill the sequence starting from activeRange_.second-1
|
|
return makeUnique<int32IndexContainer>(
|
|
activeRange_.second,
|
|
activeRange_.second+numNewPoints);
|
|
|
|
}
|
|
|
|
// second, check if there is space at the beginning
|
|
if( activeRange_.first >= numNewPoints)
|
|
{
|
|
return makeUnique<int32IndexContainer>(
|
|
0,
|
|
numNewPoints);
|
|
}
|
|
|
|
// otherwise scan the points from first to the end to find empty spaces
|
|
int32Vector newPoints(
|
|
numNewPoints,
|
|
RESERVE());
|
|
|
|
newPoints.clear();
|
|
int32 numAdded = 0;
|
|
ForAll(i, pointFlag_)
|
|
{
|
|
if(!isActive(i))
|
|
{
|
|
newPoints.push_back(static_cast<int32>(i));
|
|
numAdded++;
|
|
}
|
|
|
|
if(numAdded == numNewPoints)
|
|
{
|
|
return makeUnique<int32IndexContainer>(
|
|
newPoints.data(),
|
|
numNewPoints);
|
|
}
|
|
}
|
|
|
|
// check if there is space at the end for the remaining of points
|
|
if( numAdded <numNewPoints && capacity() - size() >= numNewPoints - numAdded )
|
|
{
|
|
int32 ind = size();
|
|
for(int32 i=numAdded; i<numNewPoints; i++)
|
|
{
|
|
newPoints.push_back(ind);
|
|
ind++;
|
|
}
|
|
|
|
return makeUnique<int32IndexContainer>(
|
|
newPoints.data(),
|
|
numNewPoints);
|
|
}
|
|
else
|
|
{
|
|
fatalErrorInFunction<<"not enough capacity for inserting particles into the point structure\n";
|
|
return nullptr;
|
|
}
|
|
|
|
return nullptr;
|
|
}
|
|
|
|
pFlow::pointStructure::pointStructure()
|
|
:
|
|
eventSubscriber(),
|
|
pointFlag_("pointFlag", "pointFlag", maxPoints_, 0 , RESERVE()),
|
|
pointPosition_("pointPosition", "pointPosition", maxPoints_, 0, RESERVE()),
|
|
activeRange_(0,0)
|
|
{}
|
|
|
|
pFlow::pointStructure::pointStructure
|
|
(
|
|
const int8Vector& flgVec,
|
|
const realx3Vector& posVec
|
|
)
|
|
:
|
|
eventSubscriber(),
|
|
maxPoints_(posVec.capacity()),
|
|
pointFlag_("pointFlag", "pointFlag", maxPoints_, 0 , RESERVE()),
|
|
pointPosition_("pointPosition", "pointPosition", maxPoints_, 0, RESERVE()),
|
|
activeRange_(0,0)
|
|
{
|
|
|
|
pointFlag_.assign(flgVec);
|
|
|
|
pointPosition_.assign(posVec);
|
|
|
|
if( !evaluatePointStructure() )
|
|
{
|
|
fatalExit;
|
|
}
|
|
}
|
|
|
|
pFlow::pointStructure::pointStructure
|
|
(
|
|
const realx3Vector& posVec
|
|
)
|
|
:
|
|
eventSubscriber(),
|
|
maxPoints_(posVec.capacity()),
|
|
pointFlag_("pointFlag", "pointFlag", maxPoints_, 0 , RESERVE()),
|
|
pointPosition_("pointPosition", "pointPosition", maxPoints_, 0, RESERVE()),
|
|
activeRange_(0,0)
|
|
{
|
|
|
|
pointPosition_.assign(posVec);
|
|
|
|
pointFlag_.resize(pointPosition_.size(), static_cast<int8>(pointStructure::ACTIVE) );
|
|
//pointFlag_.syncViews();
|
|
|
|
if( !evaluatePointStructure() )
|
|
{
|
|
fatalExit;
|
|
}
|
|
}
|
|
|
|
FUNCTION_H
|
|
const pFlow::realx3Field_D& pFlow::pointStructure::pointPosition()const
|
|
{
|
|
return pointPosition_;
|
|
}
|
|
|
|
FUNCTION_H
|
|
const pFlow::int8Field_HD& pFlow::pointStructure::pointFlag()const
|
|
{
|
|
return pointFlag_;
|
|
}
|
|
|
|
// - size of data structure
|
|
FUNCTION_H
|
|
pFlow::label pFlow::pointStructure::size()const
|
|
{
|
|
return numPoints_;
|
|
}
|
|
|
|
FUNCTION_H
|
|
pFlow::label pFlow::pointStructure::capacity()const
|
|
{
|
|
return maxPoints_;
|
|
}
|
|
|
|
FUNCTION_H
|
|
pFlow::label pFlow::pointStructure::numActive() const
|
|
{
|
|
return numActivePoints_;
|
|
}
|
|
|
|
FUNCTION_H
|
|
bool pFlow::pointStructure::allActive()const
|
|
{
|
|
return numActivePoints_ == numPoints_;
|
|
}
|
|
|
|
FUNCTION_H
|
|
bool pFlow::pointStructure::mortonSortPoints(const box& domain, real dx)
|
|
{
|
|
if( !getSortedIndex(
|
|
domain,
|
|
dx,
|
|
activeRange_,
|
|
pointPosition_.deviceVectorAll(),
|
|
pointFlag_.deviceVectorAll(),
|
|
mortonSortedIndex_) )
|
|
{
|
|
fatalErrorInFunction<<"failed to perform morton sorting!"<<endl;
|
|
return false;
|
|
}
|
|
|
|
pointPosition_.sortItems(mortonSortedIndex_);
|
|
pointFlag_.sortItems(mortonSortedIndex_);
|
|
|
|
auto oldSize = size();
|
|
auto oldCapacity = capacity();
|
|
auto oldRange = activeRange();
|
|
|
|
// update size, range, capacity
|
|
setNumMaxPoints();
|
|
activeRange_ = {0, static_cast<int>(mortonSortedIndex_.size())};
|
|
numActivePoints_ = mortonSortedIndex_.size();
|
|
|
|
eventMessage msg(eventMessage::REARRANGE);
|
|
|
|
if(oldSize != size() )
|
|
{
|
|
msg.add(eventMessage::SIZE_CHANGED);
|
|
}
|
|
|
|
if(oldCapacity != capacity())
|
|
{
|
|
msg.add(eventMessage::CAP_CHANGED);
|
|
}
|
|
|
|
if( oldRange != activeRange_)
|
|
{
|
|
msg.add(eventMessage::RANGE_CHANGED);
|
|
}
|
|
|
|
// notify all the registered objects except the exclusionList
|
|
if( !this->notify(msg) ) return false;
|
|
|
|
return true;
|
|
}
|
|
|
|
FUNCTION_H
|
|
size_t pFlow::pointStructure::markDeleteOutOfBox(const box& domain)
|
|
{
|
|
if(numPoints_==0) return 0;
|
|
|
|
int32 minRange, maxRange;
|
|
int32 numMarked =
|
|
pFlow::pointStructureKernels::markDeleteOutOfBox(
|
|
domain,
|
|
activeRange_.first,
|
|
activeRange_.second,
|
|
pointStructure::DELETED,
|
|
pointPosition_.deviceVectorAll(),
|
|
pointFlag_.deviceVectorAll(),
|
|
activePointsMaskD(),
|
|
minRange, maxRange );
|
|
|
|
if(numMarked)
|
|
{
|
|
pointFlag_.modifyOnDevice();
|
|
pointFlag_.syncViews();
|
|
}
|
|
|
|
if( numMarked<=numActivePoints_)
|
|
{
|
|
numActivePoints_ -= numMarked;
|
|
}
|
|
else
|
|
{
|
|
fatalErrorInFunction<<
|
|
"number of points marked as delete ("<<numMarked<<
|
|
") is greater than the activePoints ("<<numActivePoints_<<
|
|
").\n";
|
|
fatalExit;
|
|
}
|
|
|
|
range newRange = {minRange, maxRange};
|
|
|
|
if( activeRange_ != newRange )
|
|
{
|
|
activeRange_ = newRange;
|
|
eventMessage msg(eventMessage::RANGE_CHANGED);
|
|
|
|
// notify all the registered objects about active range change
|
|
if( !this->notify(msg) )
|
|
{
|
|
fatalExit<<
|
|
"something went wrong in notifying registered object about range change. \n";
|
|
fatalExit;
|
|
}
|
|
}
|
|
|
|
return numMarked;
|
|
|
|
}
|
|
|
|
FUNCTION_H
|
|
bool pFlow::pointStructure::updateForDelete()
|
|
{
|
|
notImplementedFunction;
|
|
|
|
return true;
|
|
}
|
|
|
|
FUNCTION_H
|
|
pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
|
|
(
|
|
const realx3Vector& pos,
|
|
const setFieldList& setField,
|
|
repository& owner,
|
|
const List<eventObserver*>& exclusionList
|
|
)
|
|
{
|
|
|
|
auto numNew = pos.size();
|
|
if( numNew==0)
|
|
{
|
|
return makeUnique<int32IndexContainer>();
|
|
}
|
|
|
|
auto newPointsPtr = getNewPointsIndices( numNew );
|
|
|
|
if(!newPointsPtr)return nullptr;
|
|
|
|
auto oldSize = size();
|
|
auto oldCapacity = capacity();
|
|
auto oldRange = activeRange();
|
|
|
|
tobeInsertedIndex_ = newPointsPtr();
|
|
|
|
// set the position of new points
|
|
|
|
if(!pointPosition_.insertSetElement(
|
|
newPointsPtr(),
|
|
pos)
|
|
)return nullptr;
|
|
|
|
|
|
if(!pointFlag_.insertSetElement(
|
|
newPointsPtr(),
|
|
static_cast<int8>(PointFlag::ACTIVE))
|
|
)return nullptr;
|
|
|
|
|
|
setNumMaxPoints();
|
|
auto minInd = newPointsPtr().min();
|
|
auto maxInd = newPointsPtr().max();
|
|
|
|
|
|
List<eventObserver*> localExlusion(exclusionList);
|
|
|
|
for(auto sfEntry:setField)
|
|
{
|
|
|
|
if(void* fieldPtr =
|
|
sfEntry.setPointFieldSelectedAll(
|
|
owner,
|
|
newPointsPtr(),
|
|
false );
|
|
fieldPtr)
|
|
|
|
localExlusion.push_back(
|
|
static_cast<eventObserver*>(fieldPtr)
|
|
);
|
|
else
|
|
return nullptr;
|
|
}
|
|
|
|
// changes the active rage based on the new inserted points
|
|
activeRange_ = { static_cast<int>(min(activeRange_.first, minInd )),
|
|
static_cast<int>(max(activeRange_.second, maxInd+1))};
|
|
|
|
numActivePoints_ += numNew;
|
|
|
|
eventMessage msg(eventMessage::INSERT);
|
|
|
|
if( oldRange != activeRange_ )
|
|
msg.add(eventMessage::RANGE_CHANGED);
|
|
|
|
if( oldSize != size() )
|
|
msg.add(eventMessage::SIZE_CHANGED);
|
|
|
|
if( oldCapacity != capacity() )
|
|
msg.add(eventMessage::CAP_CHANGED);
|
|
|
|
// notify all the registered objects except the exclusionList
|
|
if( !this->notify(msg, localExlusion) ) return nullptr;
|
|
|
|
return newPointsPtr;
|
|
}
|
|
|
|
|
|
|
|
FUNCTION_H
|
|
bool pFlow::pointStructure::readPointStructure
|
|
(
|
|
iIstream& is
|
|
)
|
|
{
|
|
auto psCapacity = is.lookupDataOrSet("pointStructureCapacity", maxSizeDefault_);
|
|
pointPosition_.reallocate(psCapacity);
|
|
pointFlag_.reallocate(psCapacity);
|
|
|
|
if( !pointPosition_.read(is))
|
|
{
|
|
ioErrorInFile(is.name(), is.lineNumber())<<
|
|
"Error in reading pointPosition in pointStructure \n";
|
|
return false;
|
|
}
|
|
|
|
if(! pointFlag_.read(is, true))
|
|
{
|
|
ioErrorInFile(is.name(), is.lineNumber())<<
|
|
"Error in reading pointFlag in pointStructure \n";
|
|
return false;
|
|
}
|
|
|
|
|
|
return evaluatePointStructure();
|
|
}
|
|
|
|
FUNCTION_H
|
|
bool pFlow::pointStructure::writePointStructure
|
|
(
|
|
iOstream& os
|
|
)const
|
|
{
|
|
os.writeWordEntry("pointStructureCapacity", maxPoints_);
|
|
|
|
if(!pointPosition_.write(os))
|
|
{
|
|
ioErrorInFile(os.name(), os.lineNumber())<<
|
|
"error in writing pointPosition to file \n";
|
|
return false;
|
|
}
|
|
|
|
if(!pointFlag_.write(os))
|
|
{
|
|
ioErrorInFile(os.name(), os.lineNumber())<<
|
|
"error in writing pointFlag to file \n";
|
|
return false;
|
|
}
|
|
return true;
|
|
}
|
|
|
|
|
|
/*pFlow::uniquePtr<pFlow::int32Vector>
|
|
pFlow::pointStructure::newPointsIndices(
|
|
int32 numNewPoints
|
|
)const
|
|
{
|
|
|
|
auto newPointsPtr = makeUnique<int32Vector>(
|
|
"pointStructure::newPointsPtr",
|
|
numNewPoints);
|
|
|
|
auto& newPoints = newPointsPtr();
|
|
|
|
|
|
// first, check if there is space at the end
|
|
if( capacity() - activeRange_.second >= numNewPoints )
|
|
{
|
|
// fill the sequence starting from activeRange_.second-1
|
|
fillSequence(newPoints, activeRange_.second-1);
|
|
return newPointsPtr;
|
|
}
|
|
|
|
// second, check if there is space at the beggining
|
|
if( activeRange_.first >= numNewPoints)
|
|
{
|
|
// fill the sequence starting from 0
|
|
fillSequence(newPoints, 0);
|
|
return newPointsPtr;
|
|
}
|
|
|
|
// otherwise scan the points from first to the end to find empty spaces
|
|
newPoints.clear();
|
|
int32 numAdded = 0;
|
|
ForAll(i, pointFlag_)
|
|
{
|
|
if(!isActive(i))
|
|
{
|
|
newPoints.push_back(static_cast<int32>(i));
|
|
numAdded++;
|
|
}
|
|
|
|
if(numAdded == numNewPoints)
|
|
{
|
|
|
|
return newPointsPtr;
|
|
}
|
|
}
|
|
|
|
// check if there is space at the end for the remaining of points
|
|
if( capacity() - size() >= numNewPoints - numAdded )
|
|
{
|
|
int32 ind = size();
|
|
for(int32 i=numAdded; i<numNewPoints; i++)
|
|
{
|
|
newPoints.push_back(ind);
|
|
ind++;
|
|
}
|
|
|
|
return newPointsPtr;
|
|
}
|
|
else
|
|
{
|
|
fatalErrorInFunction<<"not enough capacity for inserting particles into the point structure\n";
|
|
return nullptr;
|
|
}
|
|
|
|
return nullptr;
|
|
}*/ |