diff options
Diffstat (limited to 'src/geometry/cuboidGeometry3D.hh')
-rw-r--r-- | src/geometry/cuboidGeometry3D.hh | 1233 |
1 files changed, 1233 insertions, 0 deletions
diff --git a/src/geometry/cuboidGeometry3D.hh b/src/geometry/cuboidGeometry3D.hh new file mode 100644 index 0000000..a39a6cd --- /dev/null +++ b/src/geometry/cuboidGeometry3D.hh @@ -0,0 +1,1233 @@ +/* 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 + * <http://www.openlb.net/> + * + * 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 vector of 3D cuboid -- generic implementation. + */ + + +#ifndef CUBOID_GEOMETRY_3D_HH +#define CUBOID_GEOMETRY_3D_HH + + +#include <iostream> +#include <math.h> +#include <algorithm> +#include <set> +#include <limits> + +#include "geometry/cuboidGeometry3D.h" +#include "functors/analytical/indicator/indicatorF3D.h" +#include "communication/loadBalancer.h" + + +namespace olb { + +template<typename T> +CuboidGeometry3D<T>::CuboidGeometry3D() + : _motherCuboid(0,0,0,0,0,0,0), _periodicityOn(false), clout(std::cout, "CuboidGeometry3D") +{ + _cuboids.reserve(2); + add(_motherCuboid); + split(0, 1); +} + +template<typename T> +CuboidGeometry3D<T>::CuboidGeometry3D(T originX, T originY, T originZ, T deltaR, + int nX, int nY, int nZ, int nC) + : _motherCuboid(originX, originY, originZ, deltaR, nX, nY, nZ), + _periodicityOn(false), clout(std::cout, "CuboidGeometry3D") +{ + _cuboids.reserve(nC+2); + add(_motherCuboid); + split(0, nC); +} + +template<typename T> +CuboidGeometry3D<T>::CuboidGeometry3D(std::vector<T> origin, T deltaR, + std::vector<int> extent, int nC) + : CuboidGeometry3D(origin[0], origin[1], origin[2], deltaR, + extent[0], extent[1], extent[2], nC) +{ + _cuboids.reserve(nC+2); +} + +template<typename T> +CuboidGeometry3D<T>::CuboidGeometry3D(IndicatorF3D<T>& indicatorF, T voxelSize, int nC) + : _motherCuboid( 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)), + _periodicityOn(false), clout(std::cout, "CuboidGeometry3D") +{ + _cuboids.reserve(nC+2); + add(_motherCuboid); + split(0, nC); + shrink(indicatorF); +} +template<typename T> +CuboidGeometry3D<T>::CuboidGeometry3D(std::shared_ptr<IndicatorF3D<T>> indicator_sharedPtrF, T voxelSize, int nC) + : CuboidGeometry3D<T>(*indicator_sharedPtrF, voxelSize, nC) +{ +} + +template<typename T> +CuboidGeometry3D<T>::CuboidGeometry3D(IndicatorF3D<T>& indicatorF, T voxelSize, int nC, std::string minimizeBy) + : _motherCuboid( 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)), + _periodicityOn(false), clout(std::cout, "CuboidGeometry3D") +{ + _cuboids.reserve(nC+2); + add(_motherCuboid); + + if ( minimizeBy == "volume" ) { + // Search for the largest multiplier not dividable by two + int initalNc = nC; + while ( initalNc % 2 == 0 ) { + initalNc /= 2; + } + + // Split evenly in initalNc many cuboids and shrink all + split(0, initalNc); + shrink(indicatorF); + + while ( getNc() < nC ) { + // clout << "Starting with " << getNc() << " Cuboids." << std::endl; + // Search for the largest child cuboid + T iCVolume[getNc()]; + int maxiC = 0; + for (int iC = 0; iC < getNc(); iC++) { + iCVolume[iC] = get(iC).getLatticeVolume(); + if ( iCVolume[iC] > iCVolume[maxiC] ) { + maxiC = iC; + } + // clout << "Cuboid #" << iC << " with volume " << iCVolume[iC] << " and weight " << get(iC).getWeight() << std::endl; + } + // clout << "Found maximum volume cuboid #" << maxiC << " with volume " << iCVolume[maxiC] << std::endl; + + // looking for largest extend, because halfing the cuboid by its largest extend will result in the smallest surface and therfore in the least comminication cells + if ( get(maxiC).getNx() >= get(maxiC).getNy() && get(maxiC).getNx() >= get(maxiC).getNz()) { + // clout << "Cut in x direction!" << std::endl; + get(maxiC).divide(2,1,1,_cuboids); + } + else if ( get(maxiC).getNy() >= get(maxiC).getNx() && get(maxiC).getNy() >= get(maxiC).getNz()) { + // clout << "Cut in y direction!" << std::endl; + get(maxiC).divide(1,2,1,_cuboids); + } + else { + // clout << "Cut in z direction!" << std::endl; + get(maxiC).divide(1,1,2,_cuboids); + } + remove(maxiC); + // shrink the two new cuboids + shrink(_cuboids.size()-2,indicatorF); + shrink(_cuboids.size()-1,indicatorF); + } + + } + else if ( minimizeBy == "weight" ) { + + setWeights(indicatorF); + + // conduct a prime factorisation for the number of cuboids nC + std::vector<int> factors; + int initalNc = nC; + // iterate over the prime numbes from 0 to 100 (may have to be extended) + for ( int i : { + 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71 + } ) { + while ( initalNc % i == 0 ) { + initalNc /= i; + factors.push_back(i); + } + } + + // recursively split the current cuboids by each prime factor + for ( int i = factors.size() - 1; i >= 0; i-- ) { + int currentNc = _cuboids.size(); + for ( int iC = 0; iC < currentNc; iC++ ) { + // clout << "split cuboid number #" << iC << " in " << factors[i] << " parts" << std::endl; + splitByWeight(iC, factors[i], indicatorF); + shrink(indicatorF); + } + for ( int iC = 0; iC < currentNc; iC++ ) { + // clout << "delet cuboid number #" << iC << std::endl; + remove(0); + } + } + + } + +} +template<typename T> +CuboidGeometry3D<T>::CuboidGeometry3D(std::shared_ptr<IndicatorF3D<T>> indicator_sharedPtrF, T voxelSize, int nC, std::string minimizeBy) + : CuboidGeometry3D<T>(*indicator_sharedPtrF, voxelSize, nC, minimizeBy) +{ +} + +template<typename T> +CuboidGeometry3D<T>::~CuboidGeometry3D() {}; + + +template<typename T> +Cuboid3D<T>& CuboidGeometry3D<T>::get(int iC) +{ + return _cuboids[iC]; +} + +template<typename T> +Cuboid3D<T> const& CuboidGeometry3D<T>::get(int iC) const +{ + return _cuboids[iC]; +} + +template<typename T> +Cuboid3D<T> CuboidGeometry3D<T>::getMotherCuboid() +{ + return _motherCuboid; +} + +template<typename T> +Cuboid3D<T> const& CuboidGeometry3D<T>::getMotherCuboid() const +{ + return _motherCuboid; +} + +template<typename T> +void CuboidGeometry3D<T>::setPeriodicity(bool periodicityX, bool periodicityY, bool periodicityZ) +{ + _periodicityOn[0] = periodicityX; + _periodicityOn[1] = periodicityY; + _periodicityOn[2] = periodicityZ; +} + + +template<typename T> +int CuboidGeometry3D<T>::get_iC(T x, T y, T z, int offset) const +{ + unsigned i; + for (i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].checkPoint(x, y, z, offset)) { + return (int)i; + } + } + return (int)i; +} + + +template<typename T> +int CuboidGeometry3D<T>::get_iC(Vector<T,3> coords, int offset) const +{ + return get_iC(coords[0], coords[1], coords[2], offset); +} + +template<typename T> +int CuboidGeometry3D<T>::get_iC(T x, T y, T z, int orientationX, int orientationY, + int orientationZ) const +{ + unsigned i; + for (i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].checkPoint(x, y, z) && + _cuboids[i].checkPoint(x + orientationX / _cuboids[i].getDeltaR(), + y + orientationY / _cuboids[i].getDeltaR(), + z + orientationZ / _cuboids[i].getDeltaR())) { + return (int)i; + } + } + return (int)i; +} + +template<typename T> +bool CuboidGeometry3D<T>::getC(T physR[3], int& iC) const +{ + const int iCtmp = get_iC(physR[0], physR[1], physR[2]); + if (iCtmp < getNc()) { + iC = iCtmp; + return true; + } + else { + return false; + } +} + +template<typename T> +bool CuboidGeometry3D<T>::getC(std::vector<T> physR, int& iC) const +{ + const int iCtmp = get_iC(physR); + if (iCtmp < getNc()) { + iC = iCtmp; + return true; + } + else { + return false; + } +} + +template<typename T> +bool CuboidGeometry3D<T>::getC(const Vector<T,3>& physR, int& iC) const +{ + const int iCtmp = get_iC(physR); + if (iCtmp < getNc()) { + iC = iCtmp; + return true; + } + else { + return false; + } +} + +template<typename T> +bool CuboidGeometry3D<T>::getLatticeR(int latticeR[4], const T physR[3]) const +{ + int iCtmp = get_iC(physR[0], physR[1], physR[2]); + if (iCtmp < getNc()) { + latticeR[0] = iCtmp; + latticeR[1] = (int)floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() + .5); + latticeR[2] = (int)floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() + .5); + latticeR[3] = (int)floor( (physR[2] - _cuboids[latticeR[0]].getOrigin()[2] ) / _cuboids[latticeR[0]].getDeltaR() + .5); + return true; + } + else { + return false; + } +} + +template<typename T> +bool CuboidGeometry3D<T>::getFloorLatticeR(const std::vector<T>& physR, std::vector<int>& latticeR) const +{ + int iCtmp = get_iC(physR[0], physR[1], physR[2]); + if (iCtmp < getNc()) { + latticeR[0] = iCtmp; + latticeR[1] = (int)floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() ); + latticeR[2] = (int)floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() ); + latticeR[3] = (int)floor( (physR[2] - _cuboids[latticeR[0]].getOrigin()[2] ) / _cuboids[latticeR[0]].getDeltaR() ); + return true; + } + else { + return false; + } +} + +template<typename T> +bool CuboidGeometry3D<T>::getFloorLatticeR( + const Vector<T,3>& physR, Vector<int,4>& latticeR) const +{ + int iCtmp = get_iC(physR[0], physR[1], physR[2]); + if (iCtmp < getNc()) { + latticeR[0] = iCtmp; + latticeR[1] = (int)floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() ); + latticeR[2] = (int)floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() ); + latticeR[3] = (int)floor( (physR[2] - _cuboids[latticeR[0]].getOrigin()[2] ) / _cuboids[latticeR[0]].getDeltaR() ); + return true; + } + else { + return false; + } +} + +template<typename T> +void CuboidGeometry3D<T>::getPhysR(T physR[3], const int& iCglob, const int& iX, const int& iY, const int& iZ) const +{ + _cuboids[iCglob].getPhysR(physR, iX, iY, iZ); + for (int iDim = 0; iDim < 3; iDim++) { + if (_periodicityOn[iDim]) { + //std::cout << "!!! " << iDim << _periodicityOn[iDim] <<":"<< _motherCuboid.getDeltaR()*(_motherCuboid.getExtend()[iDim]) << std::endl; + physR[iDim] = remainder( physR[iDim] - _motherCuboid.getOrigin()[iDim] + + _motherCuboid.getDeltaR() * (_motherCuboid.getExtend()[iDim]), + _motherCuboid.getDeltaR() * (_motherCuboid.getExtend()[iDim])); + // solving the rounding error problem for double + if ( physR[iDim]*physR[iDim] < 0.001 * _motherCuboid.getDeltaR()*_motherCuboid.getDeltaR() ) { + physR[iDim] = T(); + } + // make it to mod instead remainer + if ( physR[iDim] < 0 ) { + physR[iDim] += _motherCuboid.getDeltaR() *( _motherCuboid.getExtend()[iDim]); + } + // add origin + physR[iDim] += _motherCuboid.getOrigin()[iDim]; + } + } + return; +} + +template<typename T> +void CuboidGeometry3D<T>::getPhysR(T physR[3], const int latticeR[4]) const +{ + getPhysR(physR, latticeR[0], latticeR[1], latticeR[2], latticeR[3]); +} + + +template<typename T> +void CuboidGeometry3D<T>::getNeighbourhood(int cuboid, std::vector<int>& neighbours, int overlap) +{ + neighbours.clear(); + + std::set<int> dummy; + + for (int iC = 0; iC < getNc(); iC++) { + if (cuboid == iC) { + continue; + } + T globX = get(iC).getOrigin()[0]; + T globY = get(iC).getOrigin()[1]; + T globZ = get(iC).getOrigin()[2]; + T nX = get(iC).getNx(); + T nY = get(iC).getNy(); + T nZ = get(iC).getNz(); + T deltaR = get(iC).getDeltaR(); + if (get(cuboid).checkInters(globX - overlap * deltaR, + globX + (nX + overlap - 1)*deltaR, + globY - overlap * deltaR, + globY + (nY + overlap - 1)*deltaR, + globZ - overlap * deltaR, + globZ + (nZ + overlap - 1)*deltaR, overlap)) { + //neighbours.push_back(iC); + dummy.insert(iC); + } + + if (_periodicityOn[0]) { + if (get(cuboid).getOrigin()[0] + (get(cuboid).getNx() + overlap - 1)*get(cuboid).getDeltaR() > getMaxPhysR()[0]) { + Cuboid3D<T> cub(get(cuboid).getOrigin()[0]-getMaxPhysR()[0], + get(cuboid).getOrigin()[1], + get(cuboid).getOrigin()[2], + get(cuboid).getDeltaR(), + get(cuboid).getNx(), + get(cuboid).getNy(), + get(cuboid).getNz()); + if (cub.checkInters(globX - overlap * deltaR, + globX + (nX + overlap - 1)*deltaR, + globY - overlap * deltaR, + globY + (nY + overlap - 1)*deltaR, + globZ - overlap * deltaR, + globZ + (nZ + overlap - 1)*deltaR, overlap)) { + dummy.insert(iC); + } + } + if (get(cuboid).getOrigin()[0] - overlap*get(cuboid).getDeltaR() < getMinPhysR()[0]) { + Cuboid3D<T> cub(get(cuboid).getOrigin()[0]+getMaxPhysR()[0], + get(cuboid).getOrigin()[1], + get(cuboid).getOrigin()[2], + get(cuboid).getDeltaR(), + get(cuboid).getNx(), + get(cuboid).getNy(), + get(cuboid).getNz()); + if (cub.checkInters(globX - overlap * deltaR, + globX + (nX + overlap - 1)*deltaR, + globY - overlap * deltaR, + globY + (nY + overlap - 1)*deltaR, + globZ - overlap * deltaR, + globZ + (nZ + overlap - 1)*deltaR, overlap)) { + dummy.insert(iC); + } + } + } + + if (_periodicityOn[1]) { + if (get(cuboid).getOrigin()[1] + (get(cuboid).getNy() + overlap - 1)*get(cuboid).getDeltaR() > getMaxPhysR()[1]) { + Cuboid3D<T> cub(get(cuboid).getOrigin()[0], + get(cuboid).getOrigin()[1]-getMaxPhysR()[1], + get(cuboid).getOrigin()[2], + get(cuboid).getDeltaR(), + get(cuboid).getNx(), + get(cuboid).getNy(), + get(cuboid).getNz()); + if (cub.checkInters(globX - overlap * deltaR, + globX + (nX + overlap - 1)*deltaR, + globY - overlap * deltaR, + globY + (nY + overlap - 1)*deltaR, + globZ - overlap * deltaR, + globZ + (nZ + overlap - 1)*deltaR, overlap)) { + dummy.insert(iC); + } + } + if (get(cuboid).getOrigin()[1] - overlap*get(cuboid).getDeltaR() < getMinPhysR()[1]) { + Cuboid3D<T> cub(get(cuboid).getOrigin()[0], + get(cuboid).getOrigin()[1]+getMaxPhysR()[1], + get(cuboid).getOrigin()[2], + get(cuboid).getDeltaR(), + get(cuboid).getNx(), + get(cuboid).getNy(), + get(cuboid).getNz()); + if (cub.checkInters(globX - overlap * deltaR, + globX + (nX + overlap - 1)*deltaR, + globY - overlap * deltaR, + globY + (nY + overlap - 1)*deltaR, + globZ - overlap * deltaR, + globZ + (nZ + overlap - 1)*deltaR, overlap)) { + dummy.insert(iC); + } + } + } + + if (_periodicityOn[2]) { + if (get(cuboid).getOrigin()[2] + (get(cuboid).getNz() + overlap - 1)*get(cuboid).getDeltaR() > getMaxPhysR()[2]) { + Cuboid3D<T> cub(get(cuboid).getOrigin()[0], + get(cuboid).getOrigin()[1], + get(cuboid).getOrigin()[2]-getMaxPhysR()[2], + get(cuboid).getDeltaR(), + get(cuboid).getNx(), + get(cuboid).getNy(), + get(cuboid).getNz()); + if (cub.checkInters(globX - overlap * deltaR, + globX + (nX + overlap - 1)*deltaR, + globY - overlap * deltaR, + globY + (nY + overlap - 1)*deltaR, + globZ - overlap * deltaR, + globZ + (nZ + overlap - 1)*deltaR, overlap)) { + dummy.insert(iC); + } + } + if (get(cuboid).getOrigin()[2] - overlap*get(cuboid).getDeltaR() < getMinPhysR()[2]) { + Cuboid3D<T> cub(get(cuboid).getOrigin()[0], + get(cuboid).getOrigin()[1], + get(cuboid).getOrigin()[2]+getMaxPhysR()[2], + get(cuboid).getDeltaR(), + get(cuboid).getNx(), + get(cuboid).getNy(), + get(cuboid).getNz()); + if (cub.checkInters(globX - overlap * deltaR, + globX + (nX + overlap - 1)*deltaR, + globY - overlap * deltaR, + globY + (nY + overlap - 1)*deltaR, + globZ - overlap * deltaR, + globZ + (nZ + overlap - 1)*deltaR, overlap)) { + dummy.insert(iC); + } + } + } + + } + std::set<int>::iterator it = dummy.begin(); + for (; it != dummy.end(); ++it) { + neighbours.push_back(*it); + } +} + +template<typename T> +int CuboidGeometry3D<T>::getNc() const +{ + return _cuboids.size(); +} + +template<typename T> +T CuboidGeometry3D<T>::getMinRatio() const +{ + T minRatio = 1.; + for (unsigned i = 0; i < _cuboids.size(); i++) { + if ((T)_cuboids[i].getNx() / (T)_cuboids[i].getNy() < minRatio) { + minRatio = (T)_cuboids[i].getNx() / (T)_cuboids[i].getNy(); + } + if ((T)_cuboids[i].getNy() / (T)_cuboids[i].getNz() < minRatio) { + minRatio = (T)_cuboids[i].getNy() / (T)_cuboids[i].getNz(); + } + if ((T)_cuboids[i].getNz() / (T)_cuboids[i].getNx() < minRatio) { + minRatio = (T)_cuboids[i].getNz() / (T)_cuboids[i].getNx(); + } + } + return minRatio; +} + +template<typename T> +T CuboidGeometry3D<T>::getMaxRatio() const +{ + T maxRatio = 1.; + for (unsigned i = 0; i < _cuboids.size(); i++) { + if ((T)_cuboids[i].getNx() / (T)_cuboids[i].getNy() > maxRatio) { + maxRatio = (T)_cuboids[i].getNx() / (T)_cuboids[i].getNy(); + } + if ((T)_cuboids[i].getNy() / (T)_cuboids[i].getNz() > maxRatio) { + maxRatio = (T)_cuboids[i].getNy() / (T)_cuboids[i].getNz(); + } + if ((T)_cuboids[i].getNz() / (T)_cuboids[i].getNx() > maxRatio) { + maxRatio = (T)_cuboids[i].getNz() / (T)_cuboids[i].getNx(); + } + } + return maxRatio; +} + +template<typename T> +Vector<T,3> CuboidGeometry3D<T>::getMinPhysR() const +{ + Vector<T,3> output (_cuboids[0].getOrigin()); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getOrigin()[0] < output[0]) { + output[0] = _cuboids[i].getOrigin()[0]; + } + if (_cuboids[i].getOrigin()[1] < output[1]) { + output[1] = _cuboids[i].getOrigin()[1]; + } + if (_cuboids[i].getOrigin()[2] < output[2]) { + output[2] = _cuboids[i].getOrigin()[2]; + } + } + return output; +} + +template<typename T> +Vector<T,3> CuboidGeometry3D<T>::getMaxPhysR() const +{ + Vector<T,3> output (_cuboids[0].getOrigin()); + output[0] += _cuboids[0].getNx()*_cuboids[0].getDeltaR(); + output[1] += _cuboids[0].getNy()*_cuboids[0].getDeltaR(); + output[2] += _cuboids[0].getNz()*_cuboids[0].getDeltaR(); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getOrigin()[0] + _cuboids[i].getNx()*_cuboids[i].getDeltaR() > output[0]) { + output[0] = _cuboids[i].getOrigin()[0] + _cuboids[i].getNx()*_cuboids[i].getDeltaR(); + } + if (_cuboids[i].getOrigin()[1] + _cuboids[i].getNy()*_cuboids[i].getDeltaR() > output[1]) { + output[1] = _cuboids[i].getOrigin()[1] + _cuboids[i].getNy()*_cuboids[i].getDeltaR(); + } + if (_cuboids[i].getOrigin()[2] + _cuboids[i].getNz()*_cuboids[i].getDeltaR() > output[2]) { + output[2] = _cuboids[i].getOrigin()[2] + _cuboids[i].getNz()*_cuboids[i].getDeltaR(); + } + } + return output; +} + +template<typename T> +T CuboidGeometry3D<T>::getMinPhysVolume() const +{ + T minVolume = _cuboids[0].getPhysVolume(); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getPhysVolume() < minVolume) { + minVolume = _cuboids[i].getPhysVolume(); + } + } + return minVolume; +} + +template<typename T> +T CuboidGeometry3D<T>::getMaxPhysVolume() const +{ + T maxVolume = _cuboids[0].getPhysVolume(); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getPhysVolume() > maxVolume) { + maxVolume = _cuboids[i].getPhysVolume(); + } + } + return maxVolume; +} + +template<typename T> +size_t CuboidGeometry3D<T>::getMinLatticeVolume() const +{ + size_t minNodes = _cuboids[0].getLatticeVolume(); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getLatticeVolume() < minNodes) { + minNodes = _cuboids[i].getLatticeVolume(); + } + } + return minNodes; +} + +template<typename T> +size_t CuboidGeometry3D<T>::getMaxLatticeVolume() const +{ + size_t maxNodes = _cuboids[0].getLatticeVolume(); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getLatticeVolume() > maxNodes) { + maxNodes = _cuboids[i].getLatticeVolume(); + } + } + return maxNodes; +} + +template<typename T> +size_t CuboidGeometry3D<T>::getMinLatticeWeight() const +{ + size_t minNodes = _cuboids[0].getWeight(); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getWeight() < minNodes) { + minNodes = _cuboids[i].getWeight(); + } + } + return minNodes; +} + +template<typename T> +size_t CuboidGeometry3D<T>::getMaxLatticeWeight() const +{ + size_t maxNodes = _cuboids[0].getWeight(); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getWeight() > maxNodes) { + maxNodes = _cuboids[i].getWeight(); + } + } + return maxNodes; +} + +template<typename T> +T CuboidGeometry3D<T>::getMinDeltaR() const +{ + T minDelta = _cuboids[0].getDeltaR(); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getDeltaR() < minDelta) { + minDelta = _cuboids[i].getDeltaR(); + } + } + return minDelta; +} + +template<typename T> +T CuboidGeometry3D<T>::getMaxDeltaR() const +{ + T maxDelta = _cuboids[0].getDeltaR(); + for (unsigned i = 0; i < _cuboids.size(); i++) { + if (_cuboids[i].getDeltaR() > maxDelta) { + maxDelta = _cuboids[i].getDeltaR(); + } + } + return maxDelta; +} + + +template<typename T> +bool CuboidGeometry3D<T>::operator==(CuboidGeometry3D<T>& rhs) +{ + return _motherCuboid == rhs._motherCuboid + && _periodicityOn == rhs._periodicityOn + && _cuboids == rhs._cuboids; +} + + + +template<typename T> +void CuboidGeometry3D<T>::add(Cuboid3D<T> cuboid) +{ + + _cuboids.push_back(cuboid); +} + +template<typename T> +void CuboidGeometry3D<T>::remove(int iC) +{ + + _cuboids.erase(_cuboids.begin() + iC); +} + + +template<typename T> +void CuboidGeometry3D<T>::remove(IndicatorF3D<T>& indicatorF) +{ + + //IndicatorIdentity3D<T> tmpIndicatorF(indicatorF); + + std::vector<bool> allZero; + int latticeR[4]; + T physR[3]; + for (unsigned iC = 0; iC < _cuboids.size(); iC++) { + latticeR[0] = iC; + allZero.push_back(true); + for (int iX = 0; iX < _cuboids[iC].getNx(); iX++) { + for (int iY = 0; iY < _cuboids[iC].getNy(); iY++) { + for (int iZ = 0; iZ < _cuboids[iC].getNz(); iZ++) { + latticeR[1] = iX; + latticeR[2] = iY; + latticeR[3] = iZ; + getPhysR(physR,latticeR); + bool inside[1]; + indicatorF(inside,physR); + if (inside[0]) { + allZero[iC] = 0; + } + } + } + } + } + for (int iC = _cuboids.size() - 1; iC >= 0; iC--) { + if (allZero[iC] ) { + remove(iC); + } + } +} + +template<typename T> +void CuboidGeometry3D<T>::shrink(int iC, IndicatorF3D<T>& indicatorF) +{ + int latticeR[4]; + T physR[3]; + bool inside[1]; + + latticeR[0] = iC; + size_t fullCells = 0; + int xN = get(iC).getNx(); + int yN = get(iC).getNy(); + int zN = get(iC).getNz(); + int maxX = 0; + int maxY = 0; + int maxZ = 0; + int newX = xN - 1; + int newY = yN - 1; + int newZ = zN - 1; + for (int iX = 0; iX < xN; iX++) { + for (int iY = 0; iY < yN; iY++) { + for (int iZ = 0; iZ < zN; iZ++) { + latticeR[1] = iX; + latticeR[2] = iY; + latticeR[3] = iZ; + getPhysR(physR,latticeR); + indicatorF(inside,physR); + if (inside[0]) { + fullCells++; + maxX = std::max(maxX, iX); + maxY = std::max(maxY, iY); + maxZ = std::max(maxZ, iZ); + newX = std::min(newX, iX); + newY = std::min(newY, iY); + newZ = std::min(newZ, iZ); + } + } + } + } + // if (maxX+2 < xN) maxX+=2; else if (maxX+1 < xN) maxX+=1; + // if (maxY+2 < yN) maxY+=2; else if (maxY+1 < yN) maxY+=1; + // if (maxZ+2 < zN) maxZ+=2; else if (maxZ+1 < zN) maxZ+=1; + // + // if (newX-2 >= 0) newX-=2; else if (newX-1 >= 0) newX-=1; + // if (newY-2 >= 0) newY-=2; else if (newY-1 >= 0) newY-=1; + // if (newZ-2 >= 0) newZ-=2; else if (newZ-1 >= 0) newZ-=1; + + if (fullCells > 0) { + get(iC).setWeight(fullCells); + _cuboids[iC].resize(newX, newY, newZ, maxX - newX + 1, maxY - newY + 1, maxZ - newZ + 1); + } + else { + remove(iC); + } +} + +template<typename T> +void CuboidGeometry3D<T>::shrink(IndicatorF3D<T>& indicatorF) +{ + for (int iC = getNc() - 1; iC >= 0; iC--) { + shrink(iC, indicatorF); + } + // shrink mother cuboid + Vector<T,3> minPhysR = getMinPhysR(); + Vector<T,3> maxPhysR = getMaxPhysR(); + T minDelataR = getMinDeltaR(); + _motherCuboid = Cuboid3D<T>(minPhysR[0], minPhysR[1], minPhysR[2], minDelataR, + (int)((maxPhysR[0]-minPhysR[0])/minDelataR + 0.5), + (int)((maxPhysR[1]-minPhysR[1])/minDelataR + 0.5), + (int)((maxPhysR[2]-minPhysR[2])/minDelataR + 0.5)); +} + + +template<typename T> +void CuboidGeometry3D<T>::split(int iC, int p) +{ + + Cuboid3D<T> temp(_cuboids[iC].getOrigin()[0], _cuboids[iC].getOrigin()[1], + _cuboids[iC].getOrigin()[2], _cuboids[iC].getDeltaR(), + _cuboids[iC].getNx(), _cuboids[iC].getNy(), _cuboids[iC].getNz()); + temp.divide(p, _cuboids); + remove(iC); +} + + +template<typename T> +void CuboidGeometry3D<T>::splitByWeight(int iC, int p, IndicatorF3D<T>& indicatorF) +{ + T averageWeight = get(iC).getWeight() / (T) p; + // clout << "Mother " << get(iC).getWeight() << " " << averageWeight << std::endl; + Cuboid3D<T> temp(_cuboids[iC].getOrigin()[0], _cuboids[iC].getOrigin()[1], + _cuboids[iC].getOrigin()[2], _cuboids[iC].getDeltaR(), + _cuboids[iC].getNx(), _cuboids[iC].getNy(), _cuboids[iC].getNz()); + + int latticeR[4]; + T physR[3]; + latticeR[0] = iC; + int xN = get(iC).getNx(); + int yN = get(iC).getNy(); + int zN = get(iC).getNz(); + T deltaR = get(iC).getDeltaR(); + int fullCells = 0; + + Vector<T, 3> globPos_child = get(iC).getOrigin(); + std::vector<int> extend_child = {xN, yN, zN}; + int localPos_child = 0; + + // looking for largest extend, because halfing the cuboid by its largest extend will result in the smallest surface and therfore in the least comminication cells + if ( get(iC).getNx() >= get(iC).getNy() && get(iC).getNx() >= get(iC).getNz()) { + // clout << "Cut in x direction!" << std::endl; + + // for each child cuboid search for the optimal cutting plane + for ( int iChild = 0; iChild < p - 1; iChild++) { + fullCells = 0; + int fullCells_minusOne = 0; + + for (int iX = localPos_child; iX < xN; iX++) { + fullCells_minusOne = fullCells; + for (int iY = 0; iY < yN; iY++) { + for (int iZ = 0; iZ < zN; iZ++) { + latticeR[1] = iX; + latticeR[2] = iY; + latticeR[3] = iZ; + getPhysR(physR,latticeR); + bool inside[1]; + indicatorF(inside,physR); + if (inside[0]) { + fullCells++; + } + } + } + // the optimal cutting plane is determined, so that the child cuboid's cells inside the indicator are the closest to the total cells inside the indicator per number of children + if ( fullCells >= averageWeight ) { + if ( (fullCells - averageWeight) > (averageWeight - fullCells_minusOne) ) { + iX--; + } + // clout << "found optimal iX = " << iX << std::endl; + extend_child[0] = iX - localPos_child + 1; + + Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]); + _cuboids.push_back(child); + + globPos_child[0] += extend_child[0]*deltaR; + localPos_child += extend_child[0] + 1; + // clout << "added child " << iChild << " of " << p << std::endl; + + break; + } + } + + } + + extend_child[0] = xN - localPos_child + p - 1; + + Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]); + _cuboids.push_back(child); + + // clout << "added last child of " << p << std::endl; + + } + else if ( get(iC).getNy() >= get(iC).getNx() && get(iC).getNy() >= get(iC).getNz()) { + // clout << "Cut in y direction!" << std::endl; + + for ( int iChild = 0; iChild < p - 1; iChild++) { + fullCells = 0; + int fullCells_minusOne = 0; + + for (int iY = localPos_child; iY < yN; iY++) { + fullCells_minusOne = fullCells; + for (int iX = 0; iX < xN; iX++) { + for (int iZ = 0; iZ < zN; iZ++) { + latticeR[1] = iX; + latticeR[2] = iY; + latticeR[3] = iZ; + getPhysR(physR,latticeR); + bool inside[1]; + indicatorF(inside,physR); + if (inside[0]) { + fullCells++; + } + } + } + if ( fullCells >= averageWeight ) { + if ( (fullCells - averageWeight) > (averageWeight - fullCells_minusOne) ) { + iY--; + } + // clout << "found optimal iY = " << iY << std::endl; + extend_child[1] = iY - localPos_child + 1; + + Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]); + _cuboids.push_back(child); + + globPos_child[1] += extend_child[1]*deltaR; + localPos_child += extend_child[1] + 1; + // clout << "added child " << iChild << " of " << p << std::endl; + + break; + } + } + + } + + extend_child[1] = yN - localPos_child + p - 1; + + Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]); + _cuboids.push_back(child); + + // clout << "added last child of " << p << std::endl; + } + else { + // clout << "Cut in z direction!" << std::endl; + + for ( int iChild = 0; iChild < p - 1; iChild++) { + fullCells = 0; + int fullCells_minusOne = 0; + + for (int iZ = localPos_child; iZ < zN; iZ++) { + fullCells_minusOne = fullCells; + for (int iY = 0; iY < yN; iY++) { + for (int iX = 0; iX < xN; iX++) { + latticeR[1] = iX; + latticeR[2] = iY; + latticeR[3] = iZ; + getPhysR(physR,latticeR); + bool inside[1]; + indicatorF(inside,physR); + if (inside[0]) { + fullCells++; + } + } + } + if ( fullCells >= averageWeight ) { + if ( (fullCells - averageWeight) > (averageWeight - fullCells_minusOne) ) { + iZ--; + } + // clout << "found optimal iZ = " << iZ << std::endl; + extend_child[2] = iZ - localPos_child + 1; + + Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]); + _cuboids.push_back(child); + + globPos_child[2] += extend_child[2]*deltaR; + localPos_child += extend_child[2] + 1; + // clout << "added child " << iChild << " of " << p << std::endl; + + break; + } + } + + } + + extend_child[2] = zN - localPos_child + p - 1; + + Cuboid3D<T> child(globPos_child[0], globPos_child[1], globPos_child[2], deltaR, extend_child[0], extend_child[1], extend_child[2]); + _cuboids.push_back(child); + + // clout << "added last child of " << p << std::endl; + } +} + + +template<typename T> +void CuboidGeometry3D<T>::swap(CuboidGeometry3D<T>& rhs) +{ + std::swap(this->_cuboids, rhs._cuboids); + std::swap(this->_motherCuboid, rhs._motherCuboid); + std::swap(this->_periodicityOn[0], rhs._periodicityOn[0]); + std::swap(this->_periodicityOn[1], rhs._periodicityOn[1]); + std::swap(this->_periodicityOn[2], rhs._periodicityOn[2]); + std::swap(t |