mortong sorting added, contactList left
This commit is contained in:
parent
b5572d3f7f
commit
93945516cb
|
@ -150,6 +150,35 @@ pFlow::particles::particles
|
|||
|
||||
}
|
||||
|
||||
bool pFlow::particles::beforeIteration()
|
||||
{
|
||||
auto domain = this->control().domain();
|
||||
|
||||
auto numMarked = dynPointStruct_.markDeleteOutOfBox(domain);
|
||||
|
||||
if(time_.sortTime())
|
||||
{
|
||||
real min_dx, max_dx;
|
||||
boundingSphereMinMax(min_dx, max_dx);
|
||||
Timer t;
|
||||
t.start();
|
||||
REPORT(0)<<"Performing morton sorting on particles ...."<<endREPORT;
|
||||
if(!pStruct().mortonSortPoints(domain, max_dx))
|
||||
{
|
||||
fatalErrorInFunction<<"Morton sorting was not successful"<<endl;
|
||||
return false;
|
||||
}
|
||||
t.end();
|
||||
REPORT(1)<<"Active range is "<< pStruct().activeRange()<<endREPORT;
|
||||
REPORT(1)<<"It took "<< yellowText(t.totalTime())<<" seconds."<<endREPORT;
|
||||
}
|
||||
|
||||
this->zeroForce();
|
||||
this->zeroTorque();
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
pFlow::uniquePtr<pFlow::List<pFlow::eventObserver*>>
|
||||
pFlow::particles::getFieldObjectList()const
|
||||
{
|
||||
|
|
|
@ -241,18 +241,7 @@ public:
|
|||
return shapeName_;
|
||||
}
|
||||
|
||||
bool beforeIteration() override
|
||||
{
|
||||
auto domain = this->control().domain();
|
||||
|
||||
auto numMarked = dynPointStruct_.markDeleteOutOfBox(domain);
|
||||
|
||||
|
||||
this->zeroForce();
|
||||
this->zeroTorque();
|
||||
|
||||
return true;
|
||||
}
|
||||
bool beforeIteration() override;
|
||||
|
||||
virtual
|
||||
bool insertParticles
|
||||
|
|
|
@ -45,12 +45,14 @@ repository/IOobject/IOobject.cpp
|
|||
repository/IOobject/IOfileHeader.cpp
|
||||
|
||||
structuredData/box/box.cpp
|
||||
structuredData/cells/cells.cpp
|
||||
structuredData/cylinder/cylinder.cpp
|
||||
structuredData/sphere/sphere.cpp
|
||||
structuredData/iBox/iBoxs.cpp
|
||||
structuredData/line/line.cpp
|
||||
structuredData/zAxis/zAxis.cpp
|
||||
structuredData/pointStructure/pointStructure.cpp
|
||||
structuredData/pointStructure/mortonIndexing.cpp
|
||||
structuredData/pointStructure/selectors/pStructSelector/pStructSelector.cpp
|
||||
structuredData/pointStructure/selectors/selectBox/selectBox.cpp
|
||||
structuredData/pointStructure/selectors/selectRange/selectRange.cpp
|
||||
|
|
|
@ -26,6 +26,8 @@ Licence:
|
|||
#include <Kokkos_DualView.hpp>
|
||||
#include <Kokkos_UnorderedMap.hpp>
|
||||
|
||||
#include "iOstream.hpp"
|
||||
|
||||
|
||||
namespace pFlow
|
||||
{
|
||||
|
@ -51,9 +53,12 @@ using DefaultExecutionSpace = Kokkos::DefaultExecutionSpace;
|
|||
template<typename T1, typename T2>
|
||||
using kPair = Kokkos::pair<T1,T2>;
|
||||
|
||||
using range = kPair<int,int>;
|
||||
template<typename T>
|
||||
using kRange = kPair<T,T>;
|
||||
|
||||
using range64 = kPair<int long,int long>;
|
||||
using range = kRange<int>;
|
||||
|
||||
using range64 = kRange<int long>;
|
||||
|
||||
template<typename T, typename... properties>
|
||||
using ViewTypeScalar = Kokkos::View<T,properties...>;
|
||||
|
@ -132,6 +137,13 @@ using deviceAtomicViewType3D =
|
|||
T***,
|
||||
Kokkos::MemoryTraits<std::is_same<DefaultExecutionSpace,Serial>::value?0:Kokkos::Atomic>>;
|
||||
|
||||
template<typename T>
|
||||
iOstream& operator <<(iOstream& os, const kRange<T>& rng)
|
||||
{
|
||||
os<<"["<<rng.first<<" "<<rng.second<<")";
|
||||
return os;
|
||||
}
|
||||
|
||||
|
||||
} // pFlow
|
||||
|
||||
|
|
|
@ -195,8 +195,8 @@ void permuteSort(const Type* first, PermuteType* pFirst, int32 numElems)
|
|||
};
|
||||
|
||||
compOperator compare{first};
|
||||
fillSequence<Type, useParallel>(pFirst, numElems, static_cast<PermuteType>(0));
|
||||
sort<Type, compOperator, useParallel>(pFirst, numElems, compare);
|
||||
fillSequence<PermuteType, useParallel>(pFirst, numElems, static_cast<PermuteType>(0));
|
||||
sort<PermuteType, compOperator, useParallel>(pFirst, numElems, compare);
|
||||
}
|
||||
|
||||
|
||||
|
|
|
@ -197,6 +197,26 @@ bool pFlow::Vector<T, Allocator>::deleteElement
|
|||
return false;
|
||||
}
|
||||
|
||||
template<typename T, typename Allocator>
|
||||
void pFlow::Vector<T, Allocator>::sortItems(
|
||||
const int32IndexContainer& indices)
|
||||
{
|
||||
if(indices.size() == 0)
|
||||
{
|
||||
this->resize(0);
|
||||
return;
|
||||
}
|
||||
size_t newSize = indices.size();
|
||||
auto hIndices = indices.hostView();
|
||||
VectorType sortedVec(name(), capacity(), newSize, RESERVE());
|
||||
|
||||
ForAll(i, hIndices)
|
||||
{
|
||||
sortedVec[i] = vectorType::operator[](i);
|
||||
}
|
||||
*this = std::move(sortedVec);
|
||||
}
|
||||
|
||||
template<typename T, typename Allocator>
|
||||
bool pFlow::Vector<T, Allocator>::insertSetElement(
|
||||
const int32IndexContainer& indices,
|
||||
|
|
|
@ -323,6 +323,9 @@ public:
|
|||
// return false if out of range
|
||||
bool deleteElement(label index);
|
||||
|
||||
/// Sort elements based on the indices
|
||||
void sortItems(const int32IndexContainer& indices);
|
||||
|
||||
// - set or insert new elements into the vector
|
||||
// return false if it fails
|
||||
bool insertSetElement(const int32IndexContainer& indices, const T& val);
|
||||
|
|
|
@ -545,6 +545,35 @@ public:
|
|||
return true;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_H
|
||||
void sortItems(const int32IndexContainer& indices)
|
||||
{
|
||||
if(indices.size() == 0)
|
||||
{
|
||||
setSize(0);
|
||||
return;
|
||||
}
|
||||
size_t newSize = indices.size();
|
||||
|
||||
deviceViewType sortedView("sortedView", newSize);
|
||||
auto dVec = deviceVectorAll();
|
||||
|
||||
auto d_indices = indices.deviceView();
|
||||
|
||||
Kokkos::parallel_for(
|
||||
"sortItems",
|
||||
newSize,
|
||||
LAMBDA_HD(int32 i){
|
||||
sortedView[i] = dVec[d_indices[i]];
|
||||
}
|
||||
);
|
||||
Kokkos::fence();
|
||||
setSize(newSize);
|
||||
copy(deviceVector(), sortedView);
|
||||
modifyOnDevice();
|
||||
syncViews();
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_H
|
||||
bool insertSetElement(const int32IndexContainer& indices, const T& val)
|
||||
{
|
||||
|
|
|
@ -558,6 +558,46 @@ public:
|
|||
return false;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_H
|
||||
void sortItems(const int32IndexContainer& indices)
|
||||
{
|
||||
if(indices.size() == 0)
|
||||
{
|
||||
setSize(0);
|
||||
return;
|
||||
}
|
||||
|
||||
size_t newSize = indices.size();
|
||||
viewType sortedView("sortedView", newSize);
|
||||
|
||||
if constexpr (isHostAccessible_)
|
||||
{
|
||||
auto h_indices = indices.hostView();
|
||||
Kokkos::parallel_for(
|
||||
"sortItems",
|
||||
newSize,
|
||||
LAMBDA_HD(int32 i){
|
||||
sortedView[i] = view_[h_indices[i]];
|
||||
});
|
||||
}else
|
||||
{
|
||||
auto d_indices = indices.deviceView();
|
||||
Kokkos::parallel_for(
|
||||
"sortItems",
|
||||
newSize,
|
||||
LAMBDA_HD(int32 i){
|
||||
sortedView[i] = view_[d_indices[i]];
|
||||
});
|
||||
}
|
||||
|
||||
Kokkos::fence();
|
||||
setSize(newSize);
|
||||
|
||||
copy(deviceVector(), sortedView);
|
||||
return;
|
||||
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_H
|
||||
bool insertSetElement(const int32IndexContainer& indices, const Vector<T>& vals)
|
||||
{
|
||||
|
|
|
@ -19,3 +19,6 @@ Licence:
|
|||
-----------------------------------------------------------------------------*/
|
||||
|
||||
#include "indexContainer.hpp"
|
||||
|
||||
|
||||
template class pFlow::indexContainer<pFlow::int32>;
|
|
@ -43,6 +43,10 @@ public:
|
|||
// - viewType of data on host
|
||||
using HostViewType = typename DualViewType::t_host;
|
||||
|
||||
using HostType = typename HostViewType::device_type;
|
||||
|
||||
using DeviceType = typename DeviceViewType::device_type;
|
||||
|
||||
template<typename ViewType>
|
||||
class IndexAccessor
|
||||
{
|
||||
|
@ -101,6 +105,10 @@ public:
|
|||
|
||||
indexContainer& operator = (const indexContainer&) = default;
|
||||
|
||||
indexContainer(indexContainer&&) = default;
|
||||
|
||||
indexContainer& operator = (indexContainer&&) = default;
|
||||
|
||||
~indexContainer() = default;
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
|
@ -150,6 +158,16 @@ public:
|
|||
return views_.d_view;
|
||||
}
|
||||
|
||||
HostViewType& hostView()
|
||||
{
|
||||
return views_.h_view;
|
||||
}
|
||||
|
||||
DeviceViewType& deviceView()
|
||||
{
|
||||
return views_.d_view;
|
||||
}
|
||||
|
||||
auto indicesHost()const
|
||||
{
|
||||
return IndexAccessor<HostViewType>(views_.h_view);
|
||||
|
@ -157,7 +175,45 @@ public:
|
|||
|
||||
auto indicesDevice()const
|
||||
{
|
||||
return IndexAccessor<DeviceViewType>(views_.d_veiw);
|
||||
return IndexAccessor<DeviceViewType>(views_.d_view);
|
||||
}
|
||||
|
||||
void modifyOnHost()
|
||||
{
|
||||
views_.modify_host();
|
||||
}
|
||||
|
||||
void modifyOnDevice()
|
||||
{
|
||||
views_.modify_device();
|
||||
}
|
||||
|
||||
void syncViews()
|
||||
{
|
||||
bool findMinMax = false;
|
||||
if(views_.template need_sync<HostType>())
|
||||
{
|
||||
Kokkos::deep_copy(views_.d_view, views_.h_view);
|
||||
findMinMax = true;
|
||||
}
|
||||
else if(views_.template need_sync<DeviceType>())
|
||||
{
|
||||
Kokkos::deep_copy(views_.h_view, views_.d_view);
|
||||
findMinMax = true;
|
||||
}
|
||||
|
||||
if(findMinMax)
|
||||
{
|
||||
min_ = pFlow::min(views_.d_view, 0, size_);
|
||||
max_ = pFlow::max(views_.d_view, 0, size_);
|
||||
}
|
||||
}
|
||||
|
||||
size_t setSize(size_t ns)
|
||||
{
|
||||
auto tmp = size_;
|
||||
size_ = ns;
|
||||
return tmp;
|
||||
}
|
||||
};
|
||||
|
||||
|
|
|
@ -124,6 +124,13 @@ bool pFlow::pointField<VectorField, T, MemorySpace>::update(const eventMessage&
|
|||
return this->insertSetElement(newElems, defaultValue_);
|
||||
}
|
||||
|
||||
if(msg.isRearranged())
|
||||
{
|
||||
auto sortedIndex = pStruct().mortonSortedIndex();
|
||||
this->sortItems(sortedIndex);
|
||||
return true;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
|
|
@ -61,8 +61,16 @@ pFlow::timeControl::timeControl
|
|||
(
|
||||
startTime_,
|
||||
dict.getValOrSet("timersReportInterval", 0.04)
|
||||
),
|
||||
performSorting_
|
||||
(
|
||||
dict.getValOrSet("performSorting", Logical("No"))
|
||||
),
|
||||
sortingInterval_
|
||||
(
|
||||
startTime_,
|
||||
dict.getValOrSet("sortingInterval", static_cast<real>(1.0))
|
||||
)
|
||||
|
||||
{
|
||||
checkForOutputToFile();
|
||||
}
|
||||
|
@ -95,6 +103,15 @@ pFlow::timeControl::timeControl(
|
|||
(
|
||||
startTime_,
|
||||
dict.getValOrSet("timersReportInterval", 0.04)
|
||||
),
|
||||
performSorting_
|
||||
(
|
||||
dict.getValOrSet("performSorting", Logical("No"))
|
||||
),
|
||||
sortingInterval_
|
||||
(
|
||||
startTime_,
|
||||
dict.getValOrSet("sortingInterval", static_cast<real>(1.0))
|
||||
)
|
||||
{
|
||||
checkForOutputToFile();
|
||||
|
@ -160,6 +177,11 @@ bool pFlow::timeControl::timersReportTime()const
|
|||
return timersReportInterval_.isMember(currentTime_, dt_);
|
||||
}
|
||||
|
||||
bool pFlow::timeControl::sortTime()const
|
||||
{
|
||||
return performSorting_()&&sortingInterval_.isMember(currentTime_,dt_);
|
||||
}
|
||||
|
||||
void pFlow::timeControl::setSaveTimeFolder(
|
||||
bool saveToFile,
|
||||
const word& timeName)
|
||||
|
|
|
@ -76,6 +76,10 @@ protected:
|
|||
|
||||
realStridedRange timersReportInterval_;
|
||||
|
||||
Logical performSorting_;
|
||||
|
||||
realStridedRange sortingInterval_;
|
||||
|
||||
int32StridedRagne screenReportInterval_ ={0,100};
|
||||
|
||||
bool outputToFile_ = false;
|
||||
|
@ -165,6 +169,8 @@ public:
|
|||
|
||||
bool timersReportTime()const;
|
||||
|
||||
bool sortTime()const;
|
||||
|
||||
bool setOutputToFile(real writeTime, const word& timeName)
|
||||
{
|
||||
if(managedExternaly_)
|
||||
|
|
|
@ -0,0 +1,21 @@
|
|||
/*------------------------------- 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 "cells.hpp"
|
|
@ -0,0 +1,259 @@
|
|||
/*------------------------------- 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.
|
||||
|
||||
-----------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef __cells_hpp__
|
||||
#define __cells_hpp__
|
||||
|
||||
|
||||
#include "types.hpp"
|
||||
#include "box.hpp"
|
||||
|
||||
namespace pFlow
|
||||
{
|
||||
|
||||
template<typename indexType>
|
||||
class cells
|
||||
{
|
||||
public:
|
||||
|
||||
using CellType = triple<indexType>;
|
||||
|
||||
protected:
|
||||
|
||||
// - domain
|
||||
box domain_{realx3(0.0), realx3(1.0)};
|
||||
|
||||
// - cell size
|
||||
realx3 cellSize_{1,1,1};
|
||||
|
||||
CellType numCells_{1,1,1};
|
||||
|
||||
|
||||
// - protected methods
|
||||
INLINE_FUNCTION_H
|
||||
void calculate()
|
||||
{
|
||||
numCells_ = (domain_.maxPoint()-domain_.minPoint())/cellSize_ + realx3(1.0);
|
||||
numCells_ = max( numCells_ , CellType(static_cast<indexType>(1)) );
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
cells()
|
||||
{}
|
||||
|
||||
INLINE_FUNCTION_H
|
||||
cells(const box& domain, real cellSize)
|
||||
:
|
||||
domain_(domain),
|
||||
cellSize_(cellSize)
|
||||
{
|
||||
calculate();
|
||||
}
|
||||
|
||||
|
||||
INLINE_FUNCTION_H
|
||||
cells(const box& domain, int32 nx, int32 ny, int32 nz)
|
||||
:
|
||||
domain_(domain),
|
||||
cellSize_(
|
||||
(domain_.maxPoint() - domain_.minPoint())/realx3(nx, ny, nz)
|
||||
),
|
||||
numCells_(nx, ny, nz)
|
||||
{}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
cells(const cells&) = default;
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
cells& operator = (const cells&) = default;
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
cells(cells &&) = default;
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
cells& operator=(cells&&) = default;
|
||||
|
||||
cells getCells()const
|
||||
{
|
||||
return *this;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_H
|
||||
void setCellSize(real cellSize)
|
||||
{
|
||||
cellSize_ = cellSize;
|
||||
calculate();
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_H
|
||||
void setCellSize(realx3 cellSize)
|
||||
{
|
||||
cellSize_ = cellSize;
|
||||
calculate();
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
realx3 cellSize()const
|
||||
{
|
||||
return cellSize_;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
const CellType& numCells()const
|
||||
{
|
||||
return numCells_;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
indexType nx()const
|
||||
{
|
||||
return numCells_.x();
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
indexType ny()const
|
||||
{
|
||||
return numCells_.y();
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
indexType nz()const
|
||||
{
|
||||
return numCells_.z();
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
int64 totalCells()const
|
||||
{
|
||||
return static_cast<int64>(numCells_.x())*
|
||||
static_cast<int64>(numCells_.y())*
|
||||
static_cast<int64>(numCells_.z());
|
||||
}
|
||||
|
||||
const auto& domain()const
|
||||
{
|
||||
return domain_;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
CellType pointIndex(const realx3& p)const
|
||||
{
|
||||
return CellType( (p - domain_.minPoint())/cellSize_ );
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
bool pointIndexInDomain(const realx3 p, CellType& index)const
|
||||
{
|
||||
if( !domain_.isInside(p) ) return false;
|
||||
|
||||
index = this->pointIndex(p);
|
||||
return true;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
bool inDomain(const realx3& p)const
|
||||
{
|
||||
return domain_.isInside(p);
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
bool isInRange(const CellType& cell)const
|
||||
{
|
||||
if(cell.x()<0)return false;
|
||||
if(cell.y()<0)return false;
|
||||
if(cell.z()<0)return false;
|
||||
if(cell.x()>numCells_.x()-1) return false;
|
||||
if(cell.y()>numCells_.y()-1) return false;
|
||||
if(cell.z()>numCells_.z()-1) return false;
|
||||
return true;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
bool isInRange(indexType i, indexType j, indexType k)const
|
||||
{
|
||||
if(i<0)return false;
|
||||
if(j<0)return false;
|
||||
if(k<0)return false;
|
||||
if(i>numCells_.x()-1) return false;
|
||||
if(j>numCells_.y()-1) return false;
|
||||
if(k>numCells_.z()-1) return false;
|
||||
return true;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
void extendBox(
|
||||
const CellType& p1,
|
||||
const CellType& p2,
|
||||
const CellType& p3,
|
||||
indexType extent,
|
||||
CellType& minP,
|
||||
CellType& maxP)const
|
||||
{
|
||||
minP = min( min( p1, p2), p3)-extent;
|
||||
maxP = max( max( p1, p2), p3)+extent;
|
||||
|
||||
minP = bound(minP);
|
||||
maxP = bound(maxP);
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
void extendBox(
|
||||
const realx3& p1,
|
||||
const realx3& p2,
|
||||
const realx3& p3,
|
||||
real extent,
|
||||
realx3& minP,
|
||||
realx3& maxP)const
|
||||
{
|
||||
minP = min(min(p1,p2),p3) - extent*cellSize_ ;
|
||||
maxP = max(max(p1,p2),p3) + extent*cellSize_ ;
|
||||
|
||||
minP = bound(minP);
|
||||
maxP = bound(maxP);
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
CellType bound(CellType p)const
|
||||
{
|
||||
return CellType(
|
||||
min( numCells_.x()-1, max(0,p.x())),
|
||||
min( numCells_.y()-1, max(0,p.y())),
|
||||
min( numCells_.z()-1, max(0,p.z()))
|
||||
);
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
realx3 bound(realx3 p)const
|
||||
{
|
||||
return realx3(
|
||||
min( domain_.maxPoint().x(), max(domain_.minPoint().x(),p.x())),
|
||||
min( domain_.maxPoint().y(), max(domain_.minPoint().y(),p.y())),
|
||||
min( domain_.maxPoint().z(), max(domain_.minPoint().z(),p.z()))
|
||||
);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif
|
|
@ -0,0 +1,75 @@
|
|||
/*------------------------------- 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 "mortonIndexing.hpp"
|
||||
#include "cells.hpp"
|
||||
|
||||
|
||||
bool pFlow::getSortedIndex(
|
||||
box boundingBox,
|
||||
real dx,
|
||||
range activeRange,
|
||||
ViewType1D<realx3> pos,
|
||||
ViewType1D<int8> flag,
|
||||
int32IndexContainer& sortedIndex)
|
||||
{
|
||||
|
||||
// obtain the morton code of the particles
|
||||
cells<size_t> allCells( boundingBox, dx);
|
||||
int32IndexContainer index(activeRange.first, activeRange.second);
|
||||
|
||||
ViewType1D<uint64_t> mortonCode("mortonCode", activeRange.second);
|
||||
|
||||
using rpMorton =
|
||||
Kokkos::RangePolicy<Kokkos::IndexType<int32>>;
|
||||
int32 numActive = 0;
|
||||
Kokkos::parallel_reduce
|
||||
(
|
||||
"mortonIndexing::getIndex::morton",
|
||||
rpMorton(activeRange.first, activeRange.second),
|
||||
LAMBDA_HD(int32 i, int32& sumToUpdate){
|
||||
if( flag[i] == 1 ) // active point
|
||||
{
|
||||
auto cellInd = allCells.pointIndex(pos[i]);
|
||||
mortonCode[i] = xyzToMortonCode64(cellInd.x(), cellInd.y(), cellInd.z());
|
||||
sumToUpdate++;
|
||||
}else
|
||||
{
|
||||
mortonCode[i] = xyzToMortonCode64(-1,-1,-1);
|
||||
}
|
||||
},
|
||||
numActive
|
||||
);
|
||||
|
||||
|
||||
permuteSort(
|
||||
mortonCode,
|
||||
activeRange.first,
|
||||
activeRange.second,
|
||||
index.deviceView(),
|
||||
0 );
|
||||
index.modifyOnDevice();
|
||||
index.setSize(numActive);
|
||||
index.syncViews();
|
||||
|
||||
sortedIndex = index;
|
||||
|
||||
return true;
|
||||
}
|
|
@ -0,0 +1,88 @@
|
|||
/*------------------------------- 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.
|
||||
|
||||
-----------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef __mortonIndexing_hpp__
|
||||
#define __mortonIndexing_hpp__
|
||||
|
||||
#include "types.hpp"
|
||||
#include "box.hpp"
|
||||
#include "indexContainer.hpp"
|
||||
|
||||
namespace pFlow
|
||||
{
|
||||
|
||||
bool getSortedIndex(
|
||||
box boundingBox,
|
||||
real dx,
|
||||
range activeRange,
|
||||
ViewType1D<realx3> pos,
|
||||
ViewType1D<int8> flag,
|
||||
int32IndexContainer& sortedIndex);
|
||||
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
uint64_t splitBy3(const uint64_t val){
|
||||
uint64_t x = val;
|
||||
x = (x | x << 32) & 0x1f00000000ffff;
|
||||
x = (x | x << 16) & 0x1f0000ff0000ff;
|
||||
x = (x | x << 8) & 0x100f00f00f00f00f;
|
||||
x = (x | x << 4) & 0x10c30c30c30c30c3;
|
||||
x = (x | x << 2) & 0x1249249249249249;
|
||||
return x;
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
uint64_t xyzToMortonCode64(uint64_t x, uint64_t y, uint64_t z)
|
||||
{
|
||||
return splitBy3(x) | (splitBy3(y) << 1) | (splitBy3(z) << 2);
|
||||
}
|
||||
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
uint64_t getThirdBits(uint64_t x)
|
||||
{
|
||||
x = x & 0x9249249249249249;
|
||||
x = (x | (x >> 2)) & 0x30c30c30c30c30c3;
|
||||
x = (x | (x >> 4)) & 0xf00f00f00f00f00f;
|
||||
x = (x | (x >> 8)) & 0x00ff0000ff0000ff;
|
||||
x = (x | (x >> 16)) & 0xffff00000000ffff;
|
||||
x = (x | (x >> 32)) & 0x00000000ffffffff;
|
||||
return x;
|
||||
|
||||
}
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
void mortonCode64Toxyz(uint64_t morton, uint64_t& x, uint64_t& y, uint64_t& z)
|
||||
{
|
||||
x = getThirdBits(morton);
|
||||
y = getThirdBits(morton >> 1);
|
||||
z = getThirdBits(morton >> 2);
|
||||
}
|
||||
|
||||
struct indexMorton
|
||||
{
|
||||
size_t morton;
|
||||
size_t index;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
#endif //__mortonIndexing_hpp__
|
|
@ -24,7 +24,8 @@ Licence:
|
|||
#include "setFieldList.hpp"
|
||||
#include "error.hpp"
|
||||
#include "iOstream.hpp"
|
||||
#include "Time.hpp"
|
||||
//#include "Time.hpp"
|
||||
#include "mortonIndexing.hpp"
|
||||
|
||||
FUNCTION_H
|
||||
bool pFlow::pointStructure::evaluatePointStructure()
|
||||
|
@ -231,6 +232,55 @@ 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)
|
||||
|
@ -306,7 +356,6 @@ pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
|
|||
)
|
||||
{
|
||||
|
||||
|
||||
auto numNew = pos.size();
|
||||
if( numNew==0)
|
||||
{
|
||||
|
@ -359,8 +408,8 @@ pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
|
|||
}
|
||||
|
||||
// changes the active rage based on the new inserted points
|
||||
activeRange_ = { min(activeRange_.first, minInd ),
|
||||
max(activeRange_.second, maxInd+1)};
|
||||
activeRange_ = { static_cast<int>(min(activeRange_.first, minInd )),
|
||||
static_cast<int>(max(activeRange_.second, maxInd+1))};
|
||||
|
||||
numActivePoints_ += numNew;
|
||||
|
||||
|
|
|
@ -163,9 +163,12 @@ protected:
|
|||
// index range of active points (half-open range)
|
||||
range activeRange_;
|
||||
|
||||
// - index vector for points to be inserted
|
||||
/// Index vector for points to be inserted
|
||||
int32IndexContainer tobeInsertedIndex_;
|
||||
|
||||
/// Sorted index of particles based on morton code
|
||||
int32IndexContainer mortonSortedIndex_;
|
||||
|
||||
|
||||
//// - protected methods
|
||||
FUNCTION_H
|
||||
|
@ -298,6 +301,10 @@ public:
|
|||
FUNCTION_H
|
||||
virtual bool updateForDelete();
|
||||
|
||||
|
||||
FUNCTION_H
|
||||
virtual bool mortonSortPoints(const box& domain, real dx);
|
||||
|
||||
///////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
// - const access to points to be newly inserted
|
||||
|
@ -320,6 +327,13 @@ public:
|
|||
}
|
||||
|
||||
|
||||
FUNCTION_H
|
||||
auto mortonSortedIndex()const
|
||||
{
|
||||
return mortonSortedIndex_;
|
||||
}
|
||||
|
||||
|
||||
// - update data structure by inserting/setting new points
|
||||
// Notifies all the fields in the registered list of data structure
|
||||
// and exclude the fields that re in the exclusionList
|
||||
|
|
Loading…
Reference in New Issue