From ceb3e0596c6609026f33bb2fde83c587b2c5e122 Mon Sep 17 00:00:00 2001 From: hamidrezanorouzi Date: Wed, 7 Sep 2022 22:22:23 +0430 Subject: [PATCH] positionParticles-ordered modified to accept cylinder&sphere region --- .gitignore | 1 + CMakeLists.txt | 2 +- .../rotatingDrum_1/caseSetup/interaction | 64 ++++++ .../caseSetup/particleInsertion | 14 ++ .../rotatingDrum_1/caseSetup/sphereShape | 11 ++ benchmarks/rotatingDrum_1/cleanThisCase | 7 + benchmarks/rotatingDrum_1/runThisCase | 21 ++ .../rotatingDrum_1/settings/geometryDict | 63 ++++++ .../rotatingDrum_1/settings/particlesDict | 44 +++++ .../rotatingDrum_1/settings/settingsDict | 36 ++++ src/phasicFlow/CMakeLists.txt | 3 + src/phasicFlow/structuredData/box/box.C | 4 +- src/phasicFlow/structuredData/box/box.H | 8 +- .../structuredData/cylinder/cylinder.C | 187 ++++++++++++++++++ .../structuredData/cylinder/cylinder.H | 158 +++++++++++++++ .../cylinderRegion/cylinderRegion.C | 73 +++++++ .../cylinderRegion/cylinderRegion.H | 62 ++++++ .../sphereRegion/sphereRegion.C | 26 +-- .../sphereRegion/sphereRegion.H | 5 +- src/phasicFlow/structuredData/sphere/sphere.C | 143 ++++++++++++++ src/phasicFlow/structuredData/sphere/sphere.H | 132 +++++++++++++ utilities/pFlowToVTK/pointFieldToVTK.H | 3 +- .../positionOrdered/positionOrdered.C | 23 ++- .../positionOrdered/positionOrdered.H | 3 - .../positionParticles/positionParticles.C | 18 +- .../positionParticles/positionParticles.H | 66 +++++++ 26 files changed, 1138 insertions(+), 39 deletions(-) create mode 100755 benchmarks/rotatingDrum_1/caseSetup/interaction create mode 100755 benchmarks/rotatingDrum_1/caseSetup/particleInsertion create mode 100755 benchmarks/rotatingDrum_1/caseSetup/sphereShape create mode 100755 benchmarks/rotatingDrum_1/cleanThisCase create mode 100755 benchmarks/rotatingDrum_1/runThisCase create mode 100644 benchmarks/rotatingDrum_1/settings/geometryDict create mode 100644 benchmarks/rotatingDrum_1/settings/particlesDict create mode 100644 benchmarks/rotatingDrum_1/settings/settingsDict create mode 100644 src/phasicFlow/structuredData/cylinder/cylinder.C create mode 100644 src/phasicFlow/structuredData/cylinder/cylinder.H create mode 100644 src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.C create mode 100644 src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.H create mode 100644 src/phasicFlow/structuredData/sphere/sphere.C create mode 100644 src/phasicFlow/structuredData/sphere/sphere.H diff --git a/.gitignore b/.gitignore index 3dbc24c1..65101b44 100644 --- a/.gitignore +++ b/.gitignore @@ -36,6 +36,7 @@ build/** include/** bin/** lib/** +test*/** **/**notnow # all possible time folders **/[0-9] diff --git a/CMakeLists.txt b/CMakeLists.txt index 6cb8f48d..9d03d0b7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -75,7 +75,7 @@ add_subdirectory(solvers) add_subdirectory(utilities) -#add_subdirectory(test) +add_subdirectory(test_newFeatures) install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H" diff --git a/benchmarks/rotatingDrum_1/caseSetup/interaction b/benchmarks/rotatingDrum_1/caseSetup/interaction new file mode 100755 index 00000000..ca9f7e4f --- /dev/null +++ b/benchmarks/rotatingDrum_1/caseSetup/interaction @@ -0,0 +1,64 @@ +/* -------------------------------*- C++ -*--------------------------------- *\ +| phasicFlow File | +| copyright: www.cemf.ir | +\* ------------------------------------------------------------------------- */ + +objectName interaction; +objectType dicrionary; + +materials (glassMat wallMat); // a list of materials names +densities (2500.0 2500); // density of materials [kg/m3] + +contactListType sortedContactList; + +model +{ + contactForceModel nonLinearLimited; + rollingFrictionModel normal; + + /* + Property (glassMat-glassMat glassMat-wallMat + wallMat-wallMat); + */ + + Yeff (1.0e6 1.0e6 // Young modulus [Pa] + 1.0e6); + + Geff (0.8e6 0.8e6 // Shear modulus [Pa] + 0.8e6); + + nu (0.25 0.25 // Poisson's ratio [-] + 0.25); + + en (0.97 0.85 // coefficient of normal restitution + 1.00); + + et (1.0 1.0 // coefficient of tangential restitution + 1.0); + + mu (0.65 0.65 // dynamic friction + 0.65); + + mur (0.1 0.1 // rolling friction + 0.1); + +} + +contactSearch +{ + method NBS; // method for broad search particle-particle + wallMapping cellsSimple; // method for broad search particle-wall + + NBSInfo + { + updateFrequency 10; // each 20 timesteps, update neighbor list + sizeRatio 1.05; // bounding box size to particle diameter (max) + } + + cellsSimpleInfo + { + updateFrequency 10; // each 20 timesteps, update neighbor list + cellExtent 0.6; // bounding box for particle-wall search (> 0.5) + } + +} \ No newline at end of file diff --git a/benchmarks/rotatingDrum_1/caseSetup/particleInsertion b/benchmarks/rotatingDrum_1/caseSetup/particleInsertion new file mode 100755 index 00000000..eec7b7f9 --- /dev/null +++ b/benchmarks/rotatingDrum_1/caseSetup/particleInsertion @@ -0,0 +1,14 @@ +/* -------------------------------*- C++ -*--------------------------------- *\ +| phasicFlow File | +| copyright: www.cemf.ir | +\* ------------------------------------------------------------------------- */ + +objectName particleInsertion; +objectType dicrionary; + + +active no; // is insertion active? + +collisionCheck No; // not implemented for yes + + diff --git a/benchmarks/rotatingDrum_1/caseSetup/sphereShape b/benchmarks/rotatingDrum_1/caseSetup/sphereShape new file mode 100755 index 00000000..a99e7f9c --- /dev/null +++ b/benchmarks/rotatingDrum_1/caseSetup/sphereShape @@ -0,0 +1,11 @@ +/* -------------------------------*- C++ -*--------------------------------- *\ +| phasicFlow File | +| copyright: www.cemf.ir | +\* ------------------------------------------------------------------------- */ + +objectName sphereDict; +objectType sphereShape; + +names (glassBead); // names of shapes +diameters (0.0055); // diameter of shapes +materials (glassMat); // material names for shapes diff --git a/benchmarks/rotatingDrum_1/cleanThisCase b/benchmarks/rotatingDrum_1/cleanThisCase new file mode 100755 index 00000000..8a0ab919 --- /dev/null +++ b/benchmarks/rotatingDrum_1/cleanThisCase @@ -0,0 +1,7 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +ls | grep -P "^(([0-9]+\.?[0-9]*)|(\.[0-9]+))$" | xargs -d"\n" rm -rf +rm -rf VTK + +#------------------------------------------------------------------------------ diff --git a/benchmarks/rotatingDrum_1/runThisCase b/benchmarks/rotatingDrum_1/runThisCase new file mode 100755 index 00000000..c48d71fe --- /dev/null +++ b/benchmarks/rotatingDrum_1/runThisCase @@ -0,0 +1,21 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory +echo "\n<--------------------------------------------------------------------->" +echo "1) Creating particles" +echo "<--------------------------------------------------------------------->\n" +particlesPhasicFlow + +echo "\n<--------------------------------------------------------------------->" +echo "2) Creating geometry" +echo "<--------------------------------------------------------------------->\n" +geometryPhasicFlow + +echo "\n<--------------------------------------------------------------------->" +echo "3) Running the case" +echo "<--------------------------------------------------------------------->\n" +sphereGranFlow + + + + +#------------------------------------------------------------------------------ diff --git a/benchmarks/rotatingDrum_1/settings/geometryDict b/benchmarks/rotatingDrum_1/settings/geometryDict new file mode 100644 index 00000000..a01bb327 --- /dev/null +++ b/benchmarks/rotatingDrum_1/settings/geometryDict @@ -0,0 +1,63 @@ +/* -------------------------------*- C++ -*--------------------------------- *\ +| phasicFlow File | +| copyright: www.cemf.ir | +\* ------------------------------------------------------------------------- */ + +objectName geometryDict; +objectType dictionary; + +motionModel rotatingAxisMotion; + +surfaces +{ + + cylinder + { + type cylinderWall; + p1 (0.0 0.0 0.0); + p2 (0.0 0.0 0.8); + radius1 0.2; + radius2 0.2; + resolution 24; + material wallMat; + motion rotAxis; + } + + + wall1 + { + type planeWall; + p1 (-0.2 -0.2 0.0); + p2 ( 0.2 -0.2 0.0); + p3 ( 0.2 0.2 0.0); + p4 (-0.2 0.2 0.0); + material wallMat; + motion rotAxis; + } + + /* + This is a plane wall at the front end of cylinder + */ + wall2 + { + type planeWall; + p1 (-0.2 -0.2 0.8); + p2 ( 0.2 -0.2 0.8); + p3 ( 0.2 0.2 0.8); + p4 (-0.2 0.2 0.8); + material wallMat; + motion rotAxis; + } + +} + +// information for rotatingAxisMotion motion model +rotatingAxisMotionInfo +{ + rotAxis + { + p1 (0.0 0.0 0.0); + p2 (0.0 0.0 1.0); + omega 1.256; // rotation speed (rad/s) => 12 rpm + } +} \ No newline at end of file diff --git a/benchmarks/rotatingDrum_1/settings/particlesDict b/benchmarks/rotatingDrum_1/settings/particlesDict new file mode 100644 index 00000000..1bcb9482 --- /dev/null +++ b/benchmarks/rotatingDrum_1/settings/particlesDict @@ -0,0 +1,44 @@ +/* -------------------------------*- C++ -*--------------------------------- *\ +| phasicFlow File | +| copyright: www.cemf.ir | +\* ------------------------------------------------------------------------- */ + +objectName particlesDict; +objectType dictionary; + +setFields +{ + + defaultValue + { + velocity realx3 (0 0 0); // linear velocity (m/s) + acceleration realx3 (0 0 0); // linear acceleration (m/s2) + rotVelocity realx3 (0 0 0); // rotational velocity (rad/s) + shapeName word glassBead; // name of the particle shape + } + + selectors + {} +} + +positionParticles +{ + method positionOrdered; + + maxNumberOfParticles 746002; + mortonSorting Yes; + + cylinder // box for positioning particles + { + p1 ( 0.0 0.0 0.0); // lower corner point of the box + p2 ( 0.0 0.0 0.8); // upper corner point of the box + radius 0.195; + } + + positionOrderedInfo + { + diameter 0.006; // minimum space between centers of particles + numPoints 440000; // number of particles in the simulation + axisOrder (z x y); // axis order for filling the space with particles + } +} diff --git a/benchmarks/rotatingDrum_1/settings/settingsDict b/benchmarks/rotatingDrum_1/settings/settingsDict new file mode 100644 index 00000000..4a67601a --- /dev/null +++ b/benchmarks/rotatingDrum_1/settings/settingsDict @@ -0,0 +1,36 @@ +/* -------------------------------*- C++ -*--------------------------------- *\ +| phasicFlow File | +| copyright: www.cemf.ir | +\* ------------------------------------------------------------------------- */ +objectName settingsDict; +objectType dictionary;; + +run rotatingDrum_1; + +dt 0.00001; // time step for integration (s) + +startTime 0; // start time for simulation + +endTime 10; // end time for simulation + +saveInterval 0.2; // time interval for saving the simulation + +timePrecision 5; // maximum number of digits for time folder + +g (0 -9.8 0); // gravity vector (m/s2) + +/* + Simulation domain + every particles that goes outside this domain is deleted. +*/ +domain +{ + min (-0.21 -0.21 -0.01); + max ( 0.21 0.21 0.81); +} + +integrationMethod AdamsBashforth3; // integration method + +timersReport Yes; // report timers? + +timersReportInterval 0.01; // time interval for reporting timers diff --git a/src/phasicFlow/CMakeLists.txt b/src/phasicFlow/CMakeLists.txt index a5420f70..bf46ae87 100644 --- a/src/phasicFlow/CMakeLists.txt +++ b/src/phasicFlow/CMakeLists.txt @@ -44,6 +44,8 @@ repository/IOobject/IOobject.C repository/IOobject/IOfileHeader.C structuredData/box/box.C +structuredData/cylinder/cylinder.C +structuredData/sphere/sphere.C structuredData/iBox/iBoxs.C structuredData/line/line.C structuredData/pointStructure/pointStructure.C @@ -54,6 +56,7 @@ structuredData/trisurfaceStructure/triSurface.C structuredData/trisurfaceStructure/multiTriSurface.C structuredData/trisurfaceStructure/stlFile.C structuredData/peakableRegion/sphereRegion/sphereRegion.C +structuredData/peakableRegion/cylinderRegion/cylinderRegion.C structuredData/peakableRegion/boxRegion/boxRegion.C structuredData/peakableRegion/peakableRegion/peakableRegion.C structuredData/peakableRegion/peakableRegions.C diff --git a/src/phasicFlow/structuredData/box/box.C b/src/phasicFlow/structuredData/box/box.C index 65846b30..731c8106 100644 --- a/src/phasicFlow/structuredData/box/box.C +++ b/src/phasicFlow/structuredData/box/box.C @@ -114,7 +114,7 @@ bool pFlow::box::write } FUNCTION_H -pFlow::iIstream& pFlow::operator << (iIstream& is, box& b) +pFlow::iIstream& pFlow::operator >>(iIstream& is, box& b) { if(! b.read(is)) { @@ -126,7 +126,7 @@ pFlow::iIstream& pFlow::operator << (iIstream& is, box& b) } FUNCTION_H -pFlow::iOstream& pFlow::operator >> (iOstream& os, const box& b) +pFlow::iOstream& pFlow::operator <<(iOstream& os, const box& b) { if(! b.write(os)) diff --git a/src/phasicFlow/structuredData/box/box.H b/src/phasicFlow/structuredData/box/box.H index 1720c734..85b25df3 100644 --- a/src/phasicFlow/structuredData/box/box.H +++ b/src/phasicFlow/structuredData/box/box.H @@ -77,13 +77,13 @@ public: } INLINE_FUNCTION_HD - const realx3& minPoint()const + realx3 minPoint()const { return min_; } INLINE_FUNCTION_HD - const realx3& maxPoint()const + realx3 maxPoint()const { return max_; } @@ -103,10 +103,10 @@ public: }; FUNCTION_H -iIstream& operator << (iIstream& is, box& b); +iIstream& operator >>(iIstream& is, box& b); FUNCTION_H -iOstream& operator >> (iOstream& os, const box& b); +iOstream& operator << (iOstream& os, const box& b); } diff --git a/src/phasicFlow/structuredData/cylinder/cylinder.C b/src/phasicFlow/structuredData/cylinder/cylinder.C new file mode 100644 index 00000000..02319a7f --- /dev/null +++ b/src/phasicFlow/structuredData/cylinder/cylinder.C @@ -0,0 +1,187 @@ +/*------------------------------- 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 "cylinder.H" + +FUNCTION_H +bool pFlow::cylinder::calculateParams() +{ + + auto p1p2 = p2_ - p1_; + + if( p1p2.length() > smallValue ) + { + axisVector2_ = dot(p1p2,p1p2); + axisVector_ = p1p2; + return true; + }else + { + return false; + } +} + +FUNCTION_H +pFlow::cylinder::cylinder( + const realx3& p1, + const realx3& p2, + const real radius) +: + p1_(p1), + p2_(p2), + radius2_(radius*radius) +{ + if(!calculateParams()) + { + fatalErrorInFunction<< + "error in the input parameters for cylinder"<("p1") + ), + p2_ + ( + dict.getVal("p2") + ) +{ + auto rad = dict.getVal("radius"); + radius2_= rad*rad; + + if(!calculateParams()) + { + fatalErrorInFunction<< + "error in the input parameters for cylinder in dictionary "<< dict.globalName()<("p1", p1_)) return false; + if(!is.nextData("p2", p2_)) return false; + real rad; + if(!is.nextData("radius", rad)) return false; + radius2_ =rad*rad; + return true; +} + +FUNCTION_H +bool pFlow::cylinder::write(iOstream& os)const +{ + os.writeWordEntry("p1", p1_); + os.writeWordEntry("p2", p2_); + os.writeWordEntry("radius", sqrt(radius2_)); + return os.check(FUNCTION_NAME); +} + +FUNCTION_H +bool pFlow::cylinder::read +( + const dictionary& dict +) +{ + p1_ = dict.getVal("p1"); + p2_ = dict.getVal("p2"); + auto rad = dict.getVal("radius"); + radius2_ = rad*rad; + return true; +} + +FUNCTION_H +bool pFlow::cylinder::write +( + dictionary& dict +)const +{ + if(!dict.add("p1", p1_)) + { + fatalErrorInFunction<< + " error in writing p1 to dictionary "<>(iIstream& is, cylinder& b) +{ + if(! b.read(is)) + { + ioErrorInFile(is.name(), is.lineNumber())<< + "error in reading cylinder. \n"; + fatalExit; + } + return is; +} + +FUNCTION_H +pFlow::iOstream& pFlow::operator << (iOstream& os, const cylinder& b) +{ + + if(! b.write(os)) + { + ioErrorInFile(os.name(), os.lineNumber())<< + "error in writing cylinder. \n"; + fatalExit; + } + return os; +} \ No newline at end of file diff --git a/src/phasicFlow/structuredData/cylinder/cylinder.H b/src/phasicFlow/structuredData/cylinder/cylinder.H new file mode 100644 index 00000000..dd93e1ee --- /dev/null +++ b/src/phasicFlow/structuredData/cylinder/cylinder.H @@ -0,0 +1,158 @@ +/*------------------------------- 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 __cylinder_H__ +#define __cylinder_H__ + +#include "types.H" +#include "dictionary.H" +#include "iIstream.H" +#include "iOstream.H" + +namespace pFlow +{ + +class cylinder +{ +protected: + + // - begin point + realx3 p1_; + + // - end point + realx3 p2_; + + // - radius^2 + real radius2_; + + realx3 axisVector_; + + real axisVector2_; + + FUNCTION_H + bool calculateParams(); + +public: + + // - type info + TypeNameNV("cylinder"); + + //// - Constructors + FUNCTION_H + cylinder(const realx3& p1, const realx3& p2, const real radius); + + FUNCTION_H + cylinder(const dictionary& dict); + + FUNCTION_H + cylinder(iIstream& is); + + FUNCTION_HD + cylinder(const cylinder&) = default; + + FUNCTION_HD + cylinder(cylinder&&) = default; + + FUNCTION_HD + cylinder& operator=(const cylinder&) = default; + + FUNCTION_HD + cylinder& operator=(cylinder&&) = default; + + ~cylinder()=default; + + //// - Methods + + INLINE_FUNCTION_HD + bool isInside(const realx3& point)const + { + auto p1Point = point-p1_; + auto H = cross(p1Point , axisVector_); + auto H2 = dot(H,H); + if( H2 < radius2_*axisVector2_) + { + real t = dot(p1Point, axisVector_)/axisVector2_; + if(t >= 0.0 && t <= 1.0) + return true; + else + return false; + } + else + { + return false; + } + + } + + INLINE_FUNCTION_HD + const realx3& p1()const + { + return p1_; + } + + INLINE_FUNCTION_HD + const realx3& p2()const + { + return p2_; + } + + INLINE_FUNCTION_HD + realx3 minPoint()const + { + return min( p1_ - realx3(radius()), p2_ - realx3(radius())); // should be improved + } + + INLINE_FUNCTION_HD + realx3 maxPoint()const + { + return max( p1_ + realx3(radius()), p2_ + realx3(radius())); // should be improved + } + + INLINE_FUNCTION_HD + real radius()const + { + return sqrt(radius2_); + } + + //// - IO operation + FUNCTION_H + bool read(iIstream & is); + + FUNCTION_H + bool write(iOstream& os)const; + + FUNCTION_H + bool read(const dictionary& dict); + + FUNCTION_H + bool write(dictionary& dict)const; +}; + +FUNCTION_H +iIstream& operator >>(iIstream& is, cylinder& b); + +FUNCTION_H +iOstream& operator << (iOstream& os, const cylinder& b); + + +} // pFlow + + +#endif // __cylinder_H__ diff --git a/src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.C b/src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.C new file mode 100644 index 00000000..a664c08a --- /dev/null +++ b/src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.C @@ -0,0 +1,73 @@ +/*------------------------------- 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 "cylinderRegion.H" + + +pFlow::cylinderRegion::cylinderRegion +( + const dictionary& dict +) +: + cylinder_(dict), + random_() +{ + +} + +bool pFlow::cylinderRegion::isInside +( + const realx3& p +) const +{ + return cylinder_.isInside(p); +} + +pFlow::realx3 pFlow::cylinderRegion::peek()const +{ + for(int32 i=0; i<100;i++) + { + auto p = + random_.randomNumber(cylinder_.minPoint(), cylinder_.maxPoint()); + if( cylinder_.isInside(p)) return p; + } + + fatalErrorInFunction<< + "cannot peek a random point from cylinderRegion. \n"; + fatalExit; + return 0; +} + + +bool pFlow::cylinderRegion::read +( + const dictionary& dict +) +{ + return cylinder_.read(dict); +} + +bool pFlow::cylinderRegion::write +( + dictionary& dict +)const +{ + return cylinder_.write(dict); +} \ No newline at end of file diff --git a/src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.H b/src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.H new file mode 100644 index 00000000..bb6a69d2 --- /dev/null +++ b/src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.H @@ -0,0 +1,62 @@ +/*------------------------------- 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 __cylinderRegion_H__ +#define __cylinderRegion_H__ + +#include "cylinder.H" +#include "uniformRandomReal.H" + +namespace pFlow +{ + +class cylinderRegion +{ +protected: + cylinder cylinder_; + + mutable uniformRandomReal random_; +public: + + // - type info + TypeNameNV("cylinderRegion"); + + cylinderRegion(const dictionary& dict); + + ~cylinderRegion() = default; + + //// - methods + bool isInside(const realx3& p) const; + + realx3 peek()const; + + + //// IO operation + bool read(const dictionary& dict); + + bool write(dictionary& dict)const; + + +}; + +} + +#endif diff --git a/src/phasicFlow/structuredData/peakableRegion/sphereRegion/sphereRegion.C b/src/phasicFlow/structuredData/peakableRegion/sphereRegion/sphereRegion.C index edecfe18..c05df7e9 100644 --- a/src/phasicFlow/structuredData/peakableRegion/sphereRegion/sphereRegion.C +++ b/src/phasicFlow/structuredData/peakableRegion/sphereRegion/sphereRegion.C @@ -26,8 +26,7 @@ pFlow::sphereRegion::sphereRegion const dictionary& dict ) : - center_(dict.getVal("center")), - radius_(dict.getVal("radius")), + sphere_(dict), random_() { @@ -38,15 +37,17 @@ bool pFlow::sphereRegion::isInside const realx3& p ) const { - return length(center_-p) <= radius_; + return sphere_.isInside(p); } pFlow::realx3 pFlow::sphereRegion::peek()const { for(int32 i=0; i<100; ++i) { - auto p = random_.randomNumber(center_ - realx3(radius_), center_ + realx3(radius_)); - if(isInside(p)) return p; + auto p = random_.randomNumber( + sphere_.center() - realx3(sphere_.radius()), + sphere_.center() + realx3(sphere_.radius())); + if(sphere_.isInside(p)) return p; } fatalErrorInFunction<< @@ -62,8 +63,8 @@ bool pFlow::sphereRegion::read const dictionary& dict ) { - center_ = dict.getVal("center"); - radius_ = dict.getVal("radius"); + sphere_.read(dict); + return true; } @@ -72,17 +73,10 @@ bool pFlow::sphereRegion::write dictionary& dict )const { - if(!dict.add("center", center_)) + if(!sphere_.write(dict)) { fatalErrorInFunction<< - " error in writing center to dictionary "<("center") + ) +{ + auto rad = dict.getVal("radius"); + radius2_= rad*rad; +} + +FUNCTION_H +pFlow::sphere::sphere +( + iIstream& is +) +{ + if( !read(is)) + { + ioErrorInFile(is.name(), is.lineNumber())<< + "error in reading sphere from file. \n"; + fatalExit; + } +} + + +FUNCTION_H +bool pFlow::sphere::read(iIstream & is) +{ + if(!is.nextData("center", center_)) return false; + real rad; + if(!is.nextData("radius", rad)) return false; + radius2_ =rad*rad; + return true; +} + +FUNCTION_H +bool pFlow::sphere::write(iOstream& os)const +{ + os.writeWordEntry("center", center_); + os.writeWordEntry("radius", sqrt(radius2_)); + return os.check(FUNCTION_NAME); +} + +FUNCTION_H +bool pFlow::sphere::read +( + const dictionary& dict +) +{ + center_ = dict.getVal("center"); + auto rad = dict.getVal("radius"); + radius2_ = rad*rad; + return true; +} + +FUNCTION_H +bool pFlow::sphere::write +( + dictionary& dict +)const +{ + if(!dict.add("center", center_)) + { + fatalErrorInFunction<< + " error in writing center to dictionary "<>(iIstream& is, sphere& b) +{ + if(! b.read(is)) + { + ioErrorInFile(is.name(), is.lineNumber())<< + "error in reading sphere. \n"; + fatalExit; + } + return is; +} + +FUNCTION_H +pFlow::iOstream& pFlow::operator << (iOstream& os, const sphere& b) +{ + + if(! b.write(os)) + { + ioErrorInFile(os.name(), os.lineNumber())<< + "error in writing sphere. \n"; + fatalExit; + } + return os; +} \ No newline at end of file diff --git a/src/phasicFlow/structuredData/sphere/sphere.H b/src/phasicFlow/structuredData/sphere/sphere.H new file mode 100644 index 00000000..cce94370 --- /dev/null +++ b/src/phasicFlow/structuredData/sphere/sphere.H @@ -0,0 +1,132 @@ +/*------------------------------- 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 __sphere_H__ +#define __sphere_H__ + +#include "types.H" +#include "dictionary.H" +#include "iIstream.H" +#include "iOstream.H" + +namespace pFlow +{ + +class sphere +{ +protected: + + // - center + realx3 center_; + + // - radius^2 + real radius2_; + + +public: + + // - type info + TypeNameNV("sphere"); + + //// - Constructors + FUNCTION_H + sphere(const realx3& center, const real radius); + + FUNCTION_H + sphere(const dictionary& dict); + + FUNCTION_H + sphere(iIstream& is); + + FUNCTION_HD + sphere(const sphere&) = default; + + FUNCTION_HD + sphere(sphere&&) = default; + + FUNCTION_HD + sphere& operator=(const sphere&) = default; + + FUNCTION_HD + sphere& operator=(sphere&&) = default; + + ~sphere()=default; + + //// - Methods + + INLINE_FUNCTION_HD + bool isInside(const realx3& point)const + { + auto cPoint = point-center_; + auto dist2 = dot(cPoint,cPoint); + return dist2 < radius2_; + + } + + INLINE_FUNCTION_HD + const realx3& center()const + { + return center_; + } + + + INLINE_FUNCTION_HD + realx3 minPoint()const + { + return center_ - realx3(radius()); + } + + INLINE_FUNCTION_HD + realx3 maxPoint()const + { + return center_ + realx3(radius()); + } + + INLINE_FUNCTION_HD + real radius()const + { + return sqrt(radius2_); + } + + //// - IO operation + FUNCTION_H + bool read(iIstream & is); + + FUNCTION_H + bool write(iOstream& os)const; + + FUNCTION_H + bool read(const dictionary& dict); + + FUNCTION_H + bool write(dictionary& dict)const; +}; + +FUNCTION_H +iIstream& operator >>(iIstream& is, sphere& b); + +FUNCTION_H +iOstream& operator << (iOstream& os, const sphere& b); + + +} // pFlow + + +#endif // __sphere_H__ diff --git a/utilities/pFlowToVTK/pointFieldToVTK.H b/utilities/pFlowToVTK/pointFieldToVTK.H index 826d2898..59e82e6d 100755 --- a/utilities/pFlowToVTK/pointFieldToVTK.H +++ b/utilities/pFlowToVTK/pointFieldToVTK.H @@ -324,7 +324,8 @@ bool convertTimeFolderPointFields( auto posVec = std::as_const(pStruct).pointPosition().hostVectorAll(); auto* pos = posVec.data(); - Report(1)<<"Writing pointStructure to vtk file."<(0.5)*dl + box_.minPoint(); - auto maxP = static_cast(0.5)*dl + box_.maxPoint(); + auto minP = region_->minPoint(); + auto maxP = region_->maxPoint(); auto cntr = minP; size_t n = 0; while( n < numPoints_ ) { - position_.push_back(cntr); + if(region_->isInside(cntr)) + { + position_.push_back(cntr); + n++; + } cntr += dl*uVector1_; @@ -111,7 +115,7 @@ bool pFlow::positionOrdered::positionPointsOrdered() } } } - n++; + } return true; @@ -139,10 +143,6 @@ pFlow::positionOrdered::positionOrdered ( poDict_.getValOrSet("axisOrder", wordList{"x", "y", "z"}) ), - box_ - ( - poDict_.subDict("box") - ), position_ ( maxNumberOfParticles_, RESERVE() @@ -154,6 +154,13 @@ pFlow::positionOrdered::positionOrdered fatalExit; } + if(!region_) + { + fatalErrorInFunction<<"You must provided a region (box, cylinder, ...) for positioning particles in dictionary "<< + dict.globalName()<(10000)); - mortonSorting_ = dict.getValOrSet("mortonSorting", Logical("Yes")); + mortonSorting_ = dict.getValOrSet("mortonSorting", Logical("Yes")); + + if( dict.containsDictionay("box") ) + { + region_ = makeUnique>(dict.subDict("box")); + } + else if(dict.containsDictionay("cylinder")) + { + region_ = makeUnique>(dict.subDict("cylinder")); + } + else if(dict.containsDictionay("sphere")) + { + region_ = makeUnique>(dict.subDict("sphere")); + } } diff --git a/utilities/particlesPhasicFlow/positionParticles/positionParticles.H b/utilities/particlesPhasicFlow/positionParticles/positionParticles.H index 41f8f7d6..cdc82706 100755 --- a/utilities/particlesPhasicFlow/positionParticles/positionParticles.H +++ b/utilities/particlesPhasicFlow/positionParticles/positionParticles.H @@ -28,11 +28,77 @@ Licence: namespace pFlow { +class regionBase +{ +public: + + regionBase() = default; + + regionBase(const regionBase&) = default; + + 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; + + +}; + +template +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(); + } + +}; class positionParticles { protected: + uniquePtr region_ = nullptr; + size_t maxNumberOfParticles_ = 10000; Logical mortonSorting_;