From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 24 Jun 2019 14:43:36 +0200 Subject: Initialize at openlb-1-3 --- src/geometry/cuboid3D.hh | 758 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 758 insertions(+) create mode 100644 src/geometry/cuboid3D.hh (limited to 'src/geometry/cuboid3D.hh') diff --git a/src/geometry/cuboid3D.hh b/src/geometry/cuboid3D.hh new file mode 100644 index 0000000..71b6123 --- /dev/null +++ b/src/geometry/cuboid3D.hh @@ -0,0 +1,758 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2007, 2014 Mathias J. Krause + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program; if not, write to the Free + * Software Foundation, Inc., 51 Franklin Street, Fifth Floor, + * Boston, MA 02110-1301, USA. +*/ + +/** \file + * The description of a single 3D cuboid -- generic implementation. + */ + + +#ifndef CUBOID_3D_HH +#define CUBOID_3D_HH + + +#include +#include "geometry/cuboidGeometry2D.h" +#include "cuboid3D.h" +#include "dynamics/lbHelpers.h" +#include "utilities/vectorHelpers.h" + + +namespace olb { + +////////////////////// Class Cuboid3D ///////////////////////// + +/*template +Cuboid3D::Cuboid3D() : clout(std::cout,"Cuboid3D") { + _weight = -1; +}*/ + +template +Cuboid3D::Cuboid3D() : Cuboid3D(0, 0, 0, 0, 0, 0, 0, 0) +{ +} + +template +Cuboid3D::Cuboid3D(T globPosX, T globPosY, T globPosZ, T delta , int nX, int nY, + int nZ, int refinementLevel) + : _weight(std::numeric_limits::max()), clout(std::cout,"Cuboid3D") +{ + init(globPosX, globPosY, globPosZ, delta, nX, nY, nZ, refinementLevel); +} + +template +Cuboid3D::Cuboid3D(std::vector origin, T delta , std::vector extend, + int refinementLevel) + : clout(std::cout,"Cuboid3D") +{ + _weight = std::numeric_limits::max(); + init(origin[0], origin[1], origin[2], delta, extend[0], extend[1], extend[2], refinementLevel); +} + +template +Cuboid3D::Cuboid3D(IndicatorF3D& indicatorF, T voxelSize, int refinementLevel) + : clout(std::cout,"Cuboid3D") +{ + _weight = std::numeric_limits::max(); + init(indicatorF.getMin()[0], indicatorF.getMin()[1], indicatorF.getMin()[2], voxelSize, + (int)((indicatorF.getMax()[0]-indicatorF.getMin()[0])/voxelSize+1.5), + (int)((indicatorF.getMax()[1]-indicatorF.getMin()[1])/voxelSize+1.5), + (int)((indicatorF.getMax()[2]-indicatorF.getMin()[2])/voxelSize+1.5), refinementLevel); +} + +// copy constructor +template +Cuboid3D::Cuboid3D(Cuboid3D const& rhs, int overlap) : clout(std::cout,"Cuboid3D") +{ + init( rhs._globPosX-rhs._delta*overlap, rhs._globPosY-rhs._delta*overlap, + rhs._globPosZ-rhs._delta*overlap, rhs._delta, rhs._nX+2*overlap, + rhs._nY+2*overlap, rhs._nZ+2*overlap); + _weight = rhs._weight; + _refinementLevel = rhs._refinementLevel; +} + + +// copy assignment +template +Cuboid3D& Cuboid3D::operator=(Cuboid3D const& rhs) +{ + init( rhs._globPosX, rhs._globPosY, rhs._globPosZ, rhs._delta, rhs._nX, + rhs._nY, rhs._nZ); + _weight = rhs._weight; + _refinementLevel = rhs._refinementLevel; + return *this; +} + +template +void Cuboid3D::init(T globPosX, T globPosY, T globPosZ, T delta, int nX, int nY, + int nZ, int refinementLevel) +{ + _globPosX = globPosX; + _globPosY = globPosY; + _globPosZ = globPosZ; + _delta = delta; + _nX = nX; + _nY = nY; + _nZ = nZ; + _refinementLevel = refinementLevel; +} + +template +std::size_t Cuboid3D::getNblock() const +{ + return 9; +} + +template +std::size_t Cuboid3D::getSerializableSize() const +{ + return ( 4 * sizeof(T) ) + ( 5 * sizeof(int) ); +} + +template +Vector const Cuboid3D::getOrigin() const +{ + return Vector (_globPosX, _globPosY, _globPosZ); +} + +template +T Cuboid3D::getDeltaR() const +{ + return _delta; +} + +template +int Cuboid3D::getNx() const +{ + return _nX; +} + +template +int Cuboid3D::getNy() const +{ + return _nY; +} + +template +int Cuboid3D::getNz() const +{ + return _nZ; +} + +template +Vector const Cuboid3D::getExtend() const +{ + return Vector (_nX, _nY, _nZ); +} + +template +T Cuboid3D::getPhysVolume() const +{ + return _nY*_nX*_nZ*_delta*_delta*_delta; +} + +template +int Cuboid3D::getWeightValue() const +{ + return _weight; +} + +template +size_t Cuboid3D::getWeight() const +{ + if (_weight == std::numeric_limits::max()) { + return getLatticeVolume(); + } else { + return _weight; + } +} + +template +void Cuboid3D::setWeight(size_t fullCells) +{ + _weight = fullCells; +} + +template +int Cuboid3D::getRefinementLevel() const +{ + return _refinementLevel; +} + +template +void Cuboid3D::setRefinementLevel(int refLevel) +{ + _refinementLevel = refLevel; +} + +template +size_t Cuboid3D::getLatticeVolume() const +{ + return static_cast(_nY)*static_cast(_nX)*static_cast(_nZ); +} + +template +T Cuboid3D::getPhysPerimeter() const +{ + return 2*_delta*_delta*(_nX*_nY + _nY*_nZ + _nZ*_nX); +} + +template +int Cuboid3D::getLatticePerimeter() const +{ + return 2*((_nX-1)*(_nY-1) + (_nY-1)*(_nZ-1) + (_nZ-1)*(_nX-1)); +} + + +template +bool Cuboid3D::operator==(const Cuboid3D& rhs) const +{ + return util::nearZero(_globPosX - rhs._globPosX) && + util::nearZero(_globPosY - rhs._globPosY) && + util::nearZero(_globPosZ - rhs._globPosZ) && + util::nearZero(_delta - rhs._delta) && + _nX == rhs._nX && + _nY == rhs._nY && + _nZ == rhs._nZ && + _weight == rhs._weight && + _refinementLevel == rhs._refinementLevel; +} + + +template +bool* Cuboid3D::getBlock(std::size_t iBlock, std::size_t& sizeBlock, bool loadingMode) +{ + std::size_t currentBlock = 0; + bool* dataPtr = nullptr; + + registerVar(iBlock, sizeBlock, currentBlock, dataPtr, _globPosX); + registerVar(iBlock, sizeBlock, currentBlock, dataPtr, _globPosY); + registerVar(iBlock, sizeBlock, currentBlock, dataPtr, _globPosZ); + registerVar(iBlock, sizeBlock, currentBlock, dataPtr, _delta); + registerVar(iBlock, sizeBlock, currentBlock, dataPtr, _nX); + registerVar(iBlock, sizeBlock, currentBlock, dataPtr, _nY); + registerVar(iBlock, sizeBlock, currentBlock, dataPtr, _nZ); + registerVar(iBlock, sizeBlock, currentBlock, dataPtr, _weight); + registerVar(iBlock, sizeBlock, currentBlock, dataPtr, _refinementLevel); + + return dataPtr; +} + +template +void Cuboid3D::print() const +{ + clout << "--------Cuboid Details----------" << std::endl; + clout << " Corner (x/y/z): " << "\t" << "(" + << _globPosX-_delta/2. << "/" << _globPosY - _delta/2. + << "/" << _globPosZ - _delta/2. << ")" << std::endl; + clout << " Delta: " << "\t" << "\t" << getDeltaR() << std::endl; + clout << " Perimeter: " << "\t" << "\t" << getPhysPerimeter() << std::endl; + clout << " Volume: " << "\t" << "\t" << getPhysVolume() << std::endl; + clout << " Extent (x/y/z): " << "\t" << "(" + << getNx() << "/" << getNy() << "/" + << getNz() << ")" << std::endl; + clout << " Nodes at Perimeter: " << "\t" << getLatticePerimeter() << std::endl; + clout << " Nodes in Volume: " << "\t" << getLatticeVolume() << std::endl; + clout << " Nodes in Indicator: " << "\t" << getWeight() << std::endl; + clout << " Other Corner: " << "\t"<< "(" << _globPosX + T(_nX-0.5)*_delta << "/" + << _globPosY + T(_nY-0.5)*_delta << "/" + << _globPosZ + T(_nZ-0.5)*_delta << ")" << std::endl; + clout << "--------------------------------" << std::endl; + +} + +template +void Cuboid3D::getPhysR(T physR[3], const int latticeR[3]) const +{ + physR[0] = _globPosX + latticeR[0]*_delta; + physR[1] = _globPosY + latticeR[1]*_delta; + physR[2] = _globPosZ + latticeR[2]*_delta; +} + +template +void Cuboid3D::getPhysR(T physR[3], const int& iX, const int& iY, const int& iZ) const +{ + physR[0] = _globPosX + iX*_delta; + physR[1] = _globPosY + iY*_delta; + physR[2] = _globPosZ + iZ*_delta; +} + +template +void Cuboid3D::getLatticeR(int latticeR[3], const T physR[3]) const +{ + latticeR[0] = (int)floor( (physR[0] - _globPosX )/_delta +.5); + latticeR[1] = (int)floor( (physR[1] - _globPosY )/_delta +.5); + latticeR[2] = (int)floor( (physR[2] - _globPosZ )/_delta +.5); +} + +template +void Cuboid3D::getLatticeR(int latticeR[3], const Vector& physR) const +{ + latticeR[0] = (int)floor( (physR[0] - _globPosX )/_delta +.5); + latticeR[1] = (int)floor( (physR[1] - _globPosY )/_delta +.5); + latticeR[2] = (int)floor( (physR[2] - _globPosZ )/_delta +.5); +} + +template +void Cuboid3D::getFloorLatticeR(const std::vector& physR, std::vector& latticeR) const +{ + getFloorLatticeR(&latticeR[0], &physR[0]); +} + +template +void Cuboid3D::getFloorLatticeR(int latticeR[3], const T physR[3]) const +{ + latticeR[0] = (int)floor( (physR[0] - _globPosX)/_delta); + latticeR[1] = (int)floor( (physR[1] - _globPosY)/_delta); + latticeR[2] = (int)floor( (physR[2] - _globPosZ)/_delta); +} + +template +bool Cuboid3D::checkPoint(T globX, T globY, T globZ, int overlap) const +{ + return _globPosX - _delta / 2. <= globX + overlap * _delta && + _globPosX + T(_nX-0.5+overlap) * _delta > globX && + _globPosY - _delta/2. <= globY + overlap * _delta && + _globPosY + T(_nY-0.5+overlap) * _delta > globY && + _globPosZ - _delta/2. <= globZ + overlap * _delta && + _globPosZ + T(_nZ-0.5+overlap) * _delta > globZ; +} + +template +bool Cuboid3D::physCheckPoint(T globX, T globY, T globZ, double overlap) const +{ + return _globPosX - T(0.5 + overlap) * _delta <= globX && + _globPosX + T(_nX-0.5+overlap)*_delta > globX && + _globPosY - T(0.5 + overlap)*_delta <= globY && + _globPosY + T(_nY-0.5+overlap)*_delta > globY && + _globPosZ - T(0.5 + overlap)*_delta <= globZ && + _globPosZ + T(_nZ-0.5+overlap)*_delta > globZ; +} + +template +bool Cuboid3D::checkPoint(T globX, T globY, T globZ, int &locX, int &locY, + int &locZ, int overlap) const +{ + if (overlap!=0) { + Cuboid3D tmp(_globPosX - overlap*_delta, _globPosY - overlap*_delta, _globPosZ - overlap*_delta, + _delta , _nX + overlap*2, _nY + overlap*2, _nZ + overlap*2); + return tmp.checkPoint(globX, globY, globZ, locX, locY, locZ); + } else if (!checkPoint(globX, globY, globZ)) { + return false; + } else { + locX = (int)floor((globX - (T)_globPosX)/_delta + .5); + locY = (int)floor((globY - (T)_globPosY)/_delta + .5); + locZ = (int)floor((globZ - (T)_globPosZ)/_delta + .5); + return true; + } +} + +template +bool Cuboid3D::checkInters(T globX0, T globX1, T globY0, T globY1, T globZ0, + T globZ1, int overlap) const +{ + T locX0d = std::max(_globPosX-overlap*_delta,globX0); + T locY0d = std::max(_globPosY-overlap*_delta,globY0); + T locZ0d = std::max(_globPosZ-overlap*_delta,globZ0); + T locX1d = std::min(_globPosX+(_nX+overlap-1)*_delta,globX1); + T locY1d = std::min(_globPosY+(_nY+overlap-1)*_delta,globY1); + T locZ1d = std::min(_globPosZ+(_nZ+overlap-1)*_delta,globZ1); + + return locX1d >= locX0d + && locY1d >= locY0d + && locZ1d >= locZ0d; +} + +template +bool Cuboid3D::checkInters(T globX, T globY, T globZ, int overlap) const +{ + return checkInters(globX, globX, globY, globY, globZ, globZ, overlap); +} + +template +bool Cuboid3D::checkInters(T globX0, T globX1, T globY0, T globY1, T globZ0, T globZ1, + int &locX0, int &locX1, int &locY0, int &locY1, int &locZ0, int &locZ1, int overlap) const +{ + if (overlap!=0) { + Cuboid3D tmp(_globPosX - overlap*_delta, _globPosY - overlap*_delta, _globPosZ - overlap*_delta, + _delta , _nX + overlap*2, _nY + overlap*2, _nZ + overlap*2); + return tmp.checkInters(globX0, globX1, globY0, globY1, globZ0, globZ1, + locX0, locX1, locY0, locY1, locZ0, locZ1); + } else if (!checkInters(globX0, globX1, globY0, globY1, globZ0, globZ1)) { + locX0 = 1; + locX1 = 0; + locY0 = 1; + locY1 = 0; + locZ0 = 1; + locZ1 = 0; + return false; + } else { + locX0 = 0; + for (int i=0; _globPosX + i*_delta < globX0; i++) { + locX0 = i+1; + } + locX1 = _nX-1; + for (int i=_nX-1; _globPosX + i*_delta > globX1; i--) { + locX1 = i-1; + } + locY0 = 0; + for (int i=0; _globPosY + i*_delta < globY0; i++) { + locY0 = i+1; + } + locY1 = _nY-1; + for (int i=_nY-1; _globPosY + i*_delta > globY1; i--) { + locY1 = i-1; + } + locZ0 = 0; + for (int i=0; _globPosZ + i*_delta < globZ0; i++) { + locZ0 = i+1; + } + locZ1 = _nZ-1; + for (int i=_nZ-1; _globPosZ + i*_delta > globZ1; i--) { + locZ1 = i-1; + } + return true; + } +} + +template +void Cuboid3D::divide(int nX, int nY, int nZ, std::vector > &childrenC) const +{ + int xN_child = 0; + int yN_child = 0; + int zN_child = 0; + + T globPosX_child = _globPosX; + T globPosY_child = _globPosY; + T globPosZ_child = _globPosZ; + + for (int iX=0; iX child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + globPosZ_child += zN_child*_delta; + } + globPosZ_child = _globPosZ; + globPosY_child += yN_child*_delta; + } + globPosY_child = _globPosY; + globPosX_child += xN_child*_delta; + } +} + + +template +void Cuboid3D::resize(int iX, int iY, int iZ, int nX, int nY, int nZ) +{ + _globPosX = _globPosX+iX*_delta; + _globPosY = _globPosY+iY*_delta; + _globPosZ = _globPosZ+iZ*_delta; + + _nX = nX; + _nY = nY; + _nZ = nZ; +} + + +template +void Cuboid3D::divide(int p, std::vector > &childrenC) const +{ + + OLB_PRECONDITION(p>0); + + int iXX = 1; + int iYY = 1; + int iZZ = p; + int nX = _nX/iXX; + int bestIx = iXX; + int nY = _nY/iYY; + int bestIy = iYY; + int nZ = _nZ/iZZ; + int bestIz = iZZ; + T bestRatio = ((T)(_nX/iXX)/(T)(_nY/iYY)-1)*((T)(_nX/iXX)/(T)(_nY/iYY)-1) + + ((T)(_nY/iYY)/(T)(_nZ/iZZ)-1)*((T)(_nY/iYY)/(T)(_nZ/iZZ)-1) + + ((T)(_nZ/iZZ)/(T)(_nX/iXX)-1)*((T)(_nZ/iZZ)/(T)(_nX/iXX)-1); + + for (int iX=1; iX<=p; iX++) { + for (int iY=1; iY*iX<=p; iY++) { + for (int iZ=p/(iX*iY); iZ*iY*iX<=p; iZ++) { + if ((iX+1)*iY*iZ>p && iX*(iY+1)*iZ>p ) { + T ratio = ((T)(_nX/iX)/(T)(_nY/iY)-1)*((T)(_nX/iX)/(T)(_nY/iY)-1) + + ((T)(_nY/iY)/(T)(_nZ/iZ)-1)*((T)(_nY/iY)/(T)(_nZ/iZ)-1) + + ((T)(_nZ/iZ)/(T)(_nX/iX)-1)*((T)(_nZ/iZ)/(T)(_nX/iX)-1); + if (rationY && nZ>nX) { + + int restY = rest%bestIy; + // split in two cuboid + if (restY==0) { + int restX = rest/bestIy; + CuboidGeometry2D helpG(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restX); + + int yN_child = 0; + T globPosY_child = _globPosY; + + for (int iY=0; iY child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + + } + globPosY_child += yN_child*_delta; + } + return; + } + + // split in four cuboid + + int restX = rest/bestIy+1; + int yN_child = 0; + T globPosY_child = _globPosY; + int splited_nY = (int) (_nY * (T)((bestIx*bestIz+restX)*restY)/(T)p); + CuboidGeometry2D helpG0(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restX); + + for (int iY=0; iY child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + } + globPosY_child += yN_child*_delta; + } + + splited_nY = _nY - splited_nY; + restX = rest/bestIy; + CuboidGeometry2D helpG1(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restX); + yN_child = 0; + + for (int iY=0; iY child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + } + globPosY_child += yN_child*_delta; + } + return; + } + + // add in x than in y direction + else if (nX>nY && nX>nZ) { + int restY = rest%bestIy; + // split in two cuboid + if (restY==0) { + int restZ = rest/bestIy; + CuboidGeometry2D helpG(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restZ); + + int yN_child = 0; + T globPosY_child = _globPosY; + + for (int iY=0; iY child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + + } + globPosY_child += yN_child*_delta; + } + return; + } + + // split in four cuboid + + int restZ = rest/bestIy+1; + + int yN_child = 0; + T globPosY_child = _globPosY; + int splited_nY = (int) (_nY * (T)((bestIx*bestIz+restZ)*restY)/(T)p); + CuboidGeometry2D helpG0(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restZ); + + for (int iY=0; iY child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + } + globPosY_child += yN_child*_delta; + } + + splited_nY = _nY - splited_nY; + restZ = rest/bestIy; + + CuboidGeometry2D helpG1(_globPosX, _globPosZ, _delta, _nX, _nZ, bestIx*bestIz+restZ); + yN_child = 0; + + for (int iY=0; iY child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + } + globPosY_child += yN_child*_delta; + } + return; + } + + // add in y than in x direction + else { + int restX = rest%bestIx; + // split in two cuboid + if (restX==0) { + int restZ = rest/bestIx; + CuboidGeometry2D helpG(_globPosZ, _globPosY, _delta, _nZ, _nY, bestIz*bestIy+restZ); + + + int xN_child = 0; + T globPosX_child = _globPosX; + + for (int iX=0; iX child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + } + globPosX_child += xN_child*_delta; + } + return; + } + + // split in four cuboid + + int restZ = rest/bestIx+1; + int xN_child = 0; + T globPosX_child = _globPosX; + int splited_nX = (int) (_nX * (T)((bestIz*bestIy+restZ)*restX)/(T)p); + CuboidGeometry2D helpG0(_globPosZ, _globPosY, _delta, _nZ, _nY, bestIz*bestIy+restZ); + + for (int iX=0; iX child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + } + globPosX_child += xN_child*_delta; + } + + splited_nX = _nX - splited_nX; + restZ = rest/bestIx; + CuboidGeometry2D helpG1(_globPosZ, _globPosY, _delta, _nZ, _nY, bestIz*bestIy+restZ); + xN_child = 0; + + for (int iX=0; iX child(globPosX_child, globPosY_child, globPosZ_child, + _delta, xN_child, yN_child, zN_child); + childrenC.push_back(child); + } + globPosX_child += xN_child*_delta; + } + return; + } + } +} + +} // namespace olb + +#endif -- cgit v1.2.3