# PhasicFlow: src/Geometry/geometry/geometry.cpp Source File
PhasicFlow is an open-source parallel DEM (discrete element method) package for simulating granular flow. It is developed in C++ and can be exectued on both GPU (like CUDA) and CPU.

geometry.cpp
Go to the documentation of this file.
/*------------------------------- 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 "geometry.hpp"
#include "systemControl.hpp"
#include "vocabs.hpp"
bool pFlow::geometry::createPropertyId()
{
if(materialName_.size() != numSurfaces() )
{
fatalErrorInFunction<<
"number of subSurface and material names do not match"<<endl;
return false;
}
uint32Vector propId(
"propId",
capacity(),
size(),
RESERVE());
ForAll(i, materialName_)
{
uint32 pIdx =0;
if( !wallProperty_.nameToIndex(materialName_[i], pIdx) )
{
fatalErrorInFunction<<
"Property/material name is invalid "<<materialName_[i]<<
". A list of valid names are \n"<< wallProperty_.materials()<<endl;
return false;
}
auto triRange = subSurfaceRange(i);
propId.fill(triRange.start(), triRange.end(), pIdx);
}
propertyId_.assign(propId);
return true;
}
void pFlow::geometry::zeroForce()
{
contactForceWall_.fill(zero3);
}
pFlow::geometry::geometry
(
systemControl& control,
const property& prop
)
:
multiTriSurface
(
objectFile
(
triSurfaceFile__,
"",
objectFile::READ_ALWAYS,
objectFile::WRITE_ALWAYS
),
&control.geometry()
),
demComponent
(
"geometry",
control
),
wallProperty_(prop),
propertyId_
(
objectFile
(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER
),
*this,
0u
),
contactForceWall_
(
objectFile
(
"contactForcWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
normalStressWall_
(
objectFile
(
"normalStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
shearStressWall_
(
objectFile
(
"shearStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
)
{
readWholeObject_ = false;
if( !IOobject::readObject() )
{
fatalErrorInFunction<<
"Error in reading from file "<<IOobject::path()<<endl;
fatalExit;
}
readWholeObject_ = true;
if( this->numSurfaces() != motionComponentName_.size() )
{
fatalErrorInFunction<<
"Number of surfaces is not equal to number of motion component names"<<endl;
fatalExit;
}
if(!createPropertyId())
{
fatalExit;
}
}
pFlow::geometry::geometry
(
systemControl &control,
const property &prop,
multiTriSurface &surf,
const wordVector &motionCompName,
const wordVector &materialName,
const dictionary& motionDict
)
:
multiTriSurface
(
objectFile
(
triSurfaceFile__,
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
&control.geometry(),
surf
),
demComponent
(
"geometry",
control
),
wallProperty_
(
prop
),
propertyId_
(
objectFile
(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER
),
*this,
0u
),
contactForceWall_
(
objectFile
(
"contactForcWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
normalStressWall_
(
objectFile
(
"normalStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
),
shearStressWall_
(
objectFile
(
"shearStressWall",
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
),
*this,
zero3
)
{
motionComponentName_.assign(motionCompName);
materialName_.assign(materialName);
if( this->numSurfaces() != motionComponentName_.size() )
{
fatalErrorInFunction<<
"Number of surfaces ("<< this->numSurfaces() <<
") is not equal to number of motion component names("<<
motionComponentName_.size()<<")"<<endl;
fatalExit;
}
if(!createPropertyId())
{
fatalExit;
}
}
bool pFlow::geometry::beforeIteration()
{
zeroForce();
return true;
}
bool pFlow::geometry::iterate()
{
return true;
}
bool pFlow::geometry::afterIteration()
{
auto numTris = size();
auto& normalsD = normals().deviceViewAll();
auto& areaD = area().deviceViewAll();
auto& cForceD = contactForceWall_.deviceViewAll();
auto& sStressD = shearStressWall_.deviceViewAll();
auto& nStressD = normalStressWall_.deviceViewAll();
Kokkos::parallel_for(
"pFlow::geometry::afterIteration",
deviceRPolicyStatic(0, numTris),
LAMBDA_HD(uint32 i)
{
realx3 n = normalsD[i];
real A = max(areaD[i],smallValue);
realx3 nF = dot(cForceD[i],n)*n;
realx3 sF = cForceD[i] - nF;
nStressD[i] = nF/A;
sStressD[i] = sF/A;
});
Kokkos::fence();
return true;
}
bool pFlow::geometry::read(iIstream &is, const IOPattern &iop)
{
motionComponentName_.read(is, iop);
materialName_.read(is, iop);
if( readWholeObject_ )
{
return multiTriSurface::read(is, iop);
}
return true;
}
bool pFlow::geometry::write(iOstream &os, const IOPattern &iop) const
{
if( !motionComponentName_.write(os, iop) )
{
fatalErrorInFunction;
return false;
}
if( !materialName_.write(os,iop))
{
fatalErrorInFunction;
return false;
}
return multiTriSurface::write(os,iop);
}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create
(
systemControl& control,
const property& prop
)
{
//
fileDictionary dict( motionModelFile__, control.time().geometry().path());
word model = dict.getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<END_REPORT;
if( systemControlvCtorSelector_.search(geomModel) )
{
auto objPtr = systemControlvCtorSelector_[geomModel] (control, prop);
REPORT(2)<<"Model "<< Green_Text(geomModel)<<" is created.\n"<<END_REPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< Yellow_Text(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
systemControlvCtorSelector_
);
fatalExit;
}
return nullptr;
}
pFlow::uniquePtr<pFlow::geometry>
pFlow::geometry::create
(
systemControl& control,
const property& prop,
multiTriSurface& surf,
const wordVector& motionCompName,
const wordVector& materialName,
const dictionary& motionDic
)
{
word model = motionDic.getVal<word>("motionModel");
auto geomModel = angleBracketsNames("geometry", model);
REPORT(1)<< "Selecting geometry model . . ."<<END_REPORT;
if( dictionaryvCtorSelector_.search(geomModel) )
{
auto objPtr = dictionaryvCtorSelector_[geomModel]
(
control,
prop,
surf,
motionCompName,
materialName,
motionDic
);
REPORT(2)<<"Model "<< Green_Text(geomModel)<<" is created.\n"<<END_REPORT;
return objPtr;
}
else
{
printKeys
(
fatalError << "Ctor Selector "<< Yellow_Text(geomModel) << " dose not exist. \n"
<<"Avaiable ones are: \n\n"
,
dictionaryvCtorSelector_
);
fatalExit;
}
return nullptr;
}
/*pFlow::geometry::geometry
(
systemControl& control,
const property& prop,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& matName
)
:
demGeometry(control),
wallProperty_(prop),
geometryRepository_(control.geometry()),
triSurface_(
control.geometry().emplaceObject<multiTriSurface>(
objectFile(
"triSurface",
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
triSurface
)
),
motionComponentName_(
control.geometry().emplaceObject<wordField>(
objectFile(
"motionComponentName",
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
"motionNamesList",
motionCompName
)
),
materialName_(
control.geometry().emplaceObject<wordField>(
objectFile(
"materialName",
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS
),
"materialNamesList",
matName
)
),
propertyId_(
control.geometry().emplaceObject<int8TriSurfaceField_D>(
objectFile(
"propertyId",
"",
objectFile::READ_NEVER,
objectFile::WRITE_NEVER),
surface(),
0 ) ),
contactForceWall_(
control.geometry().emplaceObject<realx3TriSurfaceField_D>(
objectFile(
"contactForceWall",
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS),
surface(),
zero3) ),
stressWall_(
control.geometry().emplaceObject<realx3TriSurfaceField_D>(
objectFile(
"stressWall",
"",
objectFile::READ_NEVER,
objectFile::WRITE_ALWAYS),
surface(),
zero3) )
{
if(!findPropertyId())
{
fatalExit;
}
}
pFlow::geometry::geometry
(
systemControl& control,
const property& prop,
const dictionary& dict,
const multiTriSurface& triSurface,
const wordVector& motionCompName,
const wordVector& matName
)
:
geometry(control, prop, triSurface, motionCompName, matName)
{}
bool pFlow::geometry::beforeIteration()
{
this->zeroForce();
return true;
}
bool pFlow::geometry::afterIteration()
{
auto Force = contactForceWall_.deviceVectorAll();
auto area = triSurface_.area().deviceVectorAll();
auto stress = stressWall_.deviceVectorAll();
auto numTri =triSurface_.size();
Kokkos::parallel_for(
"geometry::calculateStress",
numTri,
LAMBDA_HD(int32 i){
stress[i] = Force[i]/area[i];
});
Kokkos::fence();
return true;
}*/
Green_Text
#define Green_Text(text)
Definition: iOstream.hpp:42
pFlow::geometry::materialName_
wordField_H materialName_
Material name of each wall surface
Definition: geometry.hpp:62
pFlow::triSurface::capacity
uint32 capacity() const
Definition: triSurface.hpp:167
pFlow::real
float real
Definition: builtinTypes.hpp:45
pFlow::smallValue
const real smallValue
Definition: numericConstants.hpp:31
fatalExit
#define fatalExit
Fatal exit.
Definition: error.hpp:98
REPORT
#define REPORT(n)
Definition: streams.hpp:39
systemControl.hpp
pFlow::geometry::create
static uniquePtr< geometry > create(systemControl &control, const property &prop)
Definition: geometry.cpp:336
pFlow::uint32
unsigned int uint32
Definition: builtinTypes.hpp:56
pFlow::word
std::string word
Definition: builtinTypes.hpp:64
pFlow::systemControl
Definition: systemControl.hpp:41
pFlow::printKeys
iOstream & printKeys(iOstream &os, const wordHashMap< T > &m)
pFlow::max
T max(const internalField< T, MemorySpace > &iField)
Definition: internalFieldAlgorithms.hpp:79
pFlow::VectorSingle< word, HostSpace >::size
INLINE_FUNCTION_H uint32 size() const
Size of the vector.
Definition: wordVectorHost.hpp:231
pFlow::geometry::createPropertyId
bool createPropertyId()
Find property id of each triangle based on the supplied material name and the surface wall that the t...
Definition: geometry.cpp:26
pFlow::multiTriSurface::write
bool write(iOstream &os, const IOPattern &iop) const override
Definition: multiTriSurface.cpp:91
pFlow::zero3
const realx3 zero3(0.0)
Definition: types.hpp:137
pFlow::objectFile::WRITE_ALWAYS
@ WRITE_ALWAYS
Definition: objectFile.hpp:43
pFlow::motionModelFile__
const char *const motionModelFile__
Definition: vocabs.hpp:45
pFlow::geometry::wallProperty_
const property & wallProperty_
Const reference to physical property of materials.
Definition: geometry.hpp:54
pFlow::Time::geometry
const repository & geometry() const
Definition: Time.hpp:68
pFlow::endl
iOstream & endl(iOstream &os)
Add newline and flush stream.
Definition: iOstream.hpp:341
Yellow_Text
#define Yellow_Text(text)
Definition: iOstream.hpp:40
pFlow::multiTriSurface
Definition: multiTriSurface.hpp:32
pFlow::objectFile::WRITE_NEVER
@ WRITE_NEVER
Definition: objectFile.hpp:44
pFlow::geometry::afterIteration
bool afterIteration() override
This is called in time loop, after iterate.
Definition: geometry.cpp:275
RESERVE
Definition: Vector.hpp:40
dot
INLINE_FUNCTION_HD T dot(const quadruple< T > &oprnd1, const quadruple< T > &oprnd2)
pFlow::iIstream
Interface class for any input stream
Definition: iIstream.hpp:37
pFlow::triSurface::size
uint32 size() const
Definition: triSurface.hpp:162
pFlow::IOobject::readObject
bool readObject(bool rdHdr=true)
Definition: IOobject.cpp:83
fatalErrorInFunction
#define fatalErrorInFunction
Report a fatal error and function name and exit the application.
Definition: error.hpp:77
pFlow::property::materials
const auto & materials() const
Return list of material names.
Definition: property.hpp:117
pFlow::IOfileHeader::path
fileSystem path() const
Definition: IOfileHeader.cpp:63
pFlow::angleBracketsNames
word angleBracketsNames(const word &w1, const word &w2)
Output <w1,w2>
Definition: bTypesFunctions.cpp:156
pFlow::demComponent
A base class for every main component of DEM system.
Definition: demComponent.hpp:42
pFlow::geometry::propertyId_
uint32TriSurfaceField_D propertyId_
Property id of each triangle in the set of wall surfaces.
Definition: geometry.hpp:67
pFlow::systemControl::time
const Time & time() const
Definition: systemControl.hpp:122
pFlow::objectFile::READ_ALWAYS
@ READ_ALWAYS
Definition: objectFile.hpp:36
pFlow::IOPattern
Definition: IOPattern.hpp:32
ForAll
#define ForAll(i, container)
Definition: pFlowMacros.hpp:75
pFlow::multiTriSurface::read
bool read(iIstream &is, const IOPattern &iop) override
Definition: multiTriSurface.cpp:76
geometry.hpp
pFlow::geometry::beforeIteration
bool beforeIteration() override
This is called in time loop, before iterate.
Definition: geometry.cpp:264
pFlow::geometry::read
bool read(iIstream &is, const IOPattern &iop) override
Definition: geometry.cpp:301
pFlow::property::nameToIndex
bool nameToIndex(const word &name, uint32 &idx) const
Get the name of material in index idx Return true, if the name found, otherwise false.
Definition: property.hpp:172
pFlow::geometry::iterate
bool iterate() override
This is called in time loop.
Definition: geometry.cpp:270
END_REPORT
#define END_REPORT
Definition: streams.hpp:40
pFlow::objectFile
Definition: objectFile.hpp:30
pFlow::geometry::zeroForce
void zeroForce()
Initialize contact force to zero.
Definition: geometry.cpp:64
fatalError
#define fatalError
Report a fatal error and exit the applicaiton.
Definition: error.hpp:70
pFlow::triSurfaceFile__
const char *const triSurfaceFile__
Definition: vocabs.hpp:43
pFlow::repository::path
virtual fileSystem path() const
Definition: repository.cpp:70
pFlow::dictionary::getVal
T getVal(const word &keyword) const
get the value of data entry
Definition: dictionary.hpp:379
pFlow::property
property holds the pure properties of materials.
Definition: property.hpp:37
pFlow::objectFile::READ_NEVER
@ READ_NEVER
Definition: objectFile.hpp:37
n
uint32 n
Definition: NBSLoop.hpp:24
pFlow::systemControl::geometry
const repository & geometry() const
Definition: systemControl.hpp:132
pFlow::fileDictionary
Definition: fileDictionary.hpp:29
pFlow::uniquePtr
Definition: uniquePtr.hpp:42
LAMBDA_HD
#define LAMBDA_HD
Definition: pFlowMacros.hpp:58
pFlow::deviceRPolicyStatic
Kokkos::RangePolicy< Kokkos::DefaultExecutionSpace, Kokkos::Schedule< Kokkos::Static >, Kokkos::IndexType< pFlow::uint32 > > deviceRPolicyStatic
Definition: KokkosTypes.hpp:66
pFlow::Vector::fill
void fill(const T &val)
Definition: VectorI.hpp:23
pFlow::triSurfaceField::assign
void assign(const std::vector< T > &vals)
Definition: triSurfaceField.hpp:123
pFlow::objectFile::READ_IF_PRESENT
@ READ_IF_PRESENT
Definition: objectFile.hpp:38
vocabs.hpp
pFlow::triple< real >
pFlow::Vector< uint32 >
pFlow::iOstream
Interface class for any output stream.
Definition: iOstream.hpp:59
pFlow::multiTriSurface::numSurfaces
uint32 numSurfaces() const
Definition: multiTriSurface.hpp:81
pFlow::geometry::write
bool write(iOstream &os, const IOPattern &iop) const override
write
Definition: geometry.cpp:315
pFlow::dictionary
Dictionary holds a set of data entries or sub-dictionaries that are enclosed in a curely braces or ar...
Definition: dictionary.hpp:67
pFlow::geometry::geometry
geometry(systemControl &control, const property &prop)
Construct from controlSystem and property, for reading from file.
Definition: geometry.cpp:70
pFlow::multiTriSurface::subSurfaceRange
rangeU32 subSurfaceRange(uint32 nSub) const
Definition: multiTriSurface.hpp:91