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/cuboidGeometry3D.hh | 1233 ++++++++++++++++++++++++++++++++++++++
1 file changed, 1233 insertions(+)
create mode 100644 src/geometry/cuboidGeometry3D.hh
(limited to 'src/geometry/cuboidGeometry3D.hh')
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
+ *
+ *
+ * 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
+#include
+#include
+#include
+#include
+
+#include "geometry/cuboidGeometry3D.h"
+#include "functors/analytical/indicator/indicatorF3D.h"
+#include "communication/loadBalancer.h"
+
+
+namespace olb {
+
+template
+CuboidGeometry3D::CuboidGeometry3D()
+ : _motherCuboid(0,0,0,0,0,0,0), _periodicityOn(false), clout(std::cout, "CuboidGeometry3D")
+{
+ _cuboids.reserve(2);
+ add(_motherCuboid);
+ split(0, 1);
+}
+
+template
+CuboidGeometry3D::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
+CuboidGeometry3D::CuboidGeometry3D(std::vector origin, T deltaR,
+ std::vector extent, int nC)
+ : CuboidGeometry3D(origin[0], origin[1], origin[2], deltaR,
+ extent[0], extent[1], extent[2], nC)
+{
+ _cuboids.reserve(nC+2);
+}
+
+template
+CuboidGeometry3D::CuboidGeometry3D(IndicatorF3D& 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
+CuboidGeometry3D::CuboidGeometry3D(std::shared_ptr> indicator_sharedPtrF, T voxelSize, int nC)
+ : CuboidGeometry3D(*indicator_sharedPtrF, voxelSize, nC)
+{
+}
+
+template
+CuboidGeometry3D::CuboidGeometry3D(IndicatorF3D& 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 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
+CuboidGeometry3D::CuboidGeometry3D(std::shared_ptr> indicator_sharedPtrF, T voxelSize, int nC, std::string minimizeBy)
+ : CuboidGeometry3D(*indicator_sharedPtrF, voxelSize, nC, minimizeBy)
+{
+}
+
+template
+CuboidGeometry3D::~CuboidGeometry3D() {};
+
+
+template
+Cuboid3D& CuboidGeometry3D::get(int iC)
+{
+ return _cuboids[iC];
+}
+
+template
+Cuboid3D const& CuboidGeometry3D::get(int iC) const
+{
+ return _cuboids[iC];
+}
+
+template
+Cuboid3D CuboidGeometry3D::getMotherCuboid()
+{
+ return _motherCuboid;
+}
+
+template
+Cuboid3D const& CuboidGeometry3D::getMotherCuboid() const
+{
+ return _motherCuboid;
+}
+
+template
+void CuboidGeometry3D::setPeriodicity(bool periodicityX, bool periodicityY, bool periodicityZ)
+{
+ _periodicityOn[0] = periodicityX;
+ _periodicityOn[1] = periodicityY;
+ _periodicityOn[2] = periodicityZ;
+}
+
+
+template
+int CuboidGeometry3D::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
+int CuboidGeometry3D::get_iC(Vector coords, int offset) const
+{
+ return get_iC(coords[0], coords[1], coords[2], offset);
+}
+
+template
+int CuboidGeometry3D::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
+bool CuboidGeometry3D::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
+bool CuboidGeometry3D::getC(std::vector physR, int& iC) const
+{
+ const int iCtmp = get_iC(physR);
+ if (iCtmp < getNc()) {
+ iC = iCtmp;
+ return true;
+ }
+ else {
+ return false;
+ }
+}
+
+template
+bool CuboidGeometry3D::getC(const Vector& physR, int& iC) const
+{
+ const int iCtmp = get_iC(physR);
+ if (iCtmp < getNc()) {
+ iC = iCtmp;
+ return true;
+ }
+ else {
+ return false;
+ }
+}
+
+template
+bool CuboidGeometry3D::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
+bool CuboidGeometry3D::getFloorLatticeR(const std::vector& physR, std::vector& 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
+bool CuboidGeometry3D::getFloorLatticeR(
+ const Vector& physR, Vector& 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
+void CuboidGeometry3D::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
+void CuboidGeometry3D::getPhysR(T physR[3], const int latticeR[4]) const
+{
+ getPhysR(physR, latticeR[0], latticeR[1], latticeR[2], latticeR[3]);
+}
+
+
+template
+void CuboidGeometry3D::getNeighbourhood(int cuboid, std::vector& neighbours, int overlap)
+{
+ neighbours.clear();
+
+ std::set 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 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 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 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 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 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 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::iterator it = dummy.begin();
+ for (; it != dummy.end(); ++it) {
+ neighbours.push_back(*it);
+ }
+}
+
+template
+int CuboidGeometry3D::getNc() const
+{
+ return _cuboids.size();
+}
+
+template
+T CuboidGeometry3D::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
+T CuboidGeometry3D::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
+Vector CuboidGeometry3D::getMinPhysR() const
+{
+ Vector 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
+Vector CuboidGeometry3D::getMaxPhysR() const
+{
+ Vector 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
+T CuboidGeometry3D::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
+T CuboidGeometry3D::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
+size_t CuboidGeometry3D::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
+size_t CuboidGeometry3D::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
+size_t CuboidGeometry3D::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
+size_t CuboidGeometry3D::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
+T CuboidGeometry3D::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
+T CuboidGeometry3D::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
+bool CuboidGeometry3D::operator==(CuboidGeometry3D& rhs)
+{
+ return _motherCuboid == rhs._motherCuboid
+ && _periodicityOn == rhs._periodicityOn
+ && _cuboids == rhs._cuboids;
+}
+
+
+
+template
+void CuboidGeometry3D::add(Cuboid3D cuboid)
+{
+
+ _cuboids.push_back(cuboid);
+}
+
+template
+void CuboidGeometry3D::remove(int iC)
+{
+
+ _cuboids.erase(_cuboids.begin() + iC);
+}
+
+
+template
+void CuboidGeometry3D::remove(IndicatorF3D& indicatorF)
+{
+
+ //IndicatorIdentity3D tmpIndicatorF(indicatorF);
+
+ std::vector 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
+void CuboidGeometry3D::shrink(int iC, IndicatorF3D& 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
+void CuboidGeometry3D::shrink(IndicatorF3D& indicatorF)
+{
+ for (int iC = getNc() - 1; iC >= 0; iC--) {
+ shrink(iC, indicatorF);
+ }
+ // shrink mother cuboid
+ Vector minPhysR = getMinPhysR();
+ Vector maxPhysR = getMaxPhysR();
+ T minDelataR = getMinDeltaR();
+ _motherCuboid = Cuboid3D(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
+void CuboidGeometry3D::split(int iC, int p)
+{
+
+ Cuboid3D 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
+void CuboidGeometry3D::splitByWeight(int iC, int p, IndicatorF3D& indicatorF)
+{
+ T averageWeight = get(iC).getWeight() / (T) p;
+ // clout << "Mother " << get(iC).getWeight() << " " << averageWeight << std::endl;
+ Cuboid3D 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 globPos_child = get(iC).getOrigin();
+ std::vector 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 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 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 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 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 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 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
+void CuboidGeometry3D::swap(CuboidGeometry3D& 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(this->clout, rhs.clout);
+}
+
+template
+void CuboidGeometry3D::swapCuboids(std::vector< Cuboid3D >& cuboids)
+{
+ _cuboids.swap(cuboids);
+}
+
+template
+void CuboidGeometry3D::replaceCuboids(std::vector< Cuboid3D >& cuboids)
+{
+ this->_cuboids.clear();
+ for ( unsigned iC = 0; iC < cuboids.size(); iC++) {
+ add(cuboids[iC]);
+ }
+}
+
+
+template
+void CuboidGeometry3D::setWeights(IndicatorF3D& indicatorF)
+{
+
+ //IndicatorIdentity3D tmpIndicatorF(indicatorF);
+ int nC = getNc();
+ int latticeR[4];
+ T physR[3];
+ for ( int iC = nC - 1; iC >= 0; iC--) { // assemble neighbourhood information
+ latticeR[0] = iC;
+ int xN = get(iC).getNx();
+ int yN = get(iC).getNy();
+ int zN = get(iC).getNz();
+ size_t fullCells = 0;
+ 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);
+ bool inside[1];
+ indicatorF(inside,physR);
+ if (inside[0]) {
+ fullCells++;
+ }
+ }
+ }
+ }
+ if (fullCells > 0) {
+ get(iC).setWeight(fullCells);
+ }
+ else {
+ remove(iC);
+ }
+ }
+}
+
+
+template
+size_t CuboidGeometry3D::getNblock() const
+{
+ return 1 // _periodicityOn
+ + _motherCuboid.getNblock() // _motherCuboid
+ + _cuboids.size() > 0 ? 1 + _cuboids.size() * _cuboids[0].getNblock() : 0; // _cuboids;
+}
+
+
+template
+size_t CuboidGeometry3D::getSerializableSize() const
+{
+ return 3 * sizeof(bool) // _periodicityOn
+ + _motherCuboid.getSerializableSize() // _motherCuboid
+ + (_cuboids.size() > 0 ?
+ sizeof(size_t) + _cuboids.size() * _cuboids[0].getSerializableSize() :
+ 0); // _cuboids;
+}
+
+template
+bool* CuboidGeometry3D::getBlock(std::size_t iBlock, std::size_t& sizeBlock, bool loadingMode)
+{
+ std::size_t currentBlock = 0;
+ size_t sizeBufferIndex = 0;
+ bool* dataPtr = nullptr;
+
+ registerVar (iBlock, sizeBlock, currentBlock, dataPtr, _periodicityOn[0], 3);
+ registerSerializableOfConstSize (iBlock, sizeBlock, currentBlock, dataPtr, _motherCuboid, loadingMode);
+ registerStdVectorOfSerializablesOfConstSize (iBlock, sizeBlock, currentBlock, sizeBufferIndex, dataPtr,
+ _cuboids, loadingMode);
+
+ return dataPtr;
+}
+
+
+template
+void CuboidGeometry3D::print() const
+{
+ clout << "---Cuboid Stucture Statistics---" << std::endl;
+ clout << " Number of Cuboids: " << "\t" << getNc() << std::endl;
+ clout << " Delta (min): " << "\t" << "\t" << getMinDeltaR() << std::endl;
+ clout << " (max): " << "\t" << "\t" << getMaxDeltaR() << std::endl;
+ clout << " Ratio (min): " << "\t" << "\t" << getMinRatio() << std::endl;
+ clout << " (max): " << "\t" << "\t" << getMaxRatio() << std::endl;
+ clout << " Nodes (min): " << "\t" << "\t" << getMinLatticeVolume() << std::endl;
+ clout << " (max): " << "\t" << "\t" << getMaxLatticeVolume() << std::endl;
+ clout << " Weight (min): " << "\t" << "\t" << getMinLatticeWeight() << std::endl;
+ clout << " (max): " << "\t" << "\t" << getMaxLatticeWeight() << std::endl;
+ clout << "--------------------------------" << std::endl;
+}
+
+template
+void CuboidGeometry3D::printExtended()
+{
+ clout << "Mothercuboid :" << std::endl;
+ getMotherCuboid().print();
+
+ for (int iC = 0; iC < getNc(); iC++) {
+ clout << "Cuboid #" << iC << ": " << std::endl;
+ get(iC).print();
+ }
+}
+
+template
+void CuboidGeometry3D::writeToExistingFile(std::string completeFileName, LoadBalancer& loadBalancer)
+{
+ std::ofstream fout;
+ if ( singleton::mpi().isMainProcessor() ) {
+
+ // Open File
+ fout.open(completeFileName.c_str(), std::ios::app);
+ if (!fout) {
+ clout << "Error: could not open " << completeFileName << std::endl;
+ }
+
+ // --- Preamble --- //
+ fout << "\n";
+
+ // TODO: Move Cuboid XML Serialization to Cuboid3D class
+ for (int iC = 0; iC < getNc(); ++iC) {
+ fout << "\n";
+ }
+
+ fout << "\n";
+
+ // Close File
+ fout.close();
+ }
+}
+
+
+template
+void CuboidGeometry3D::writeToFile(std::string fileName, LoadBalancer& loadBalancer)
+{
+ std::string fname = singleton::directories().getLogOutDir() + fileName + ".xml";
+ std::ofstream fout;
+ if (singleton::mpi().isMainProcessor()) {
+ fout.open(fname.c_str(), std::ios::trunc);
+ fout << "\n";
+ fout << "\n";
+ fout.close();
+ fout.clear();
+ }
+
+ writeToExistingFile(fname, loadBalancer);
+
+ if (singleton::mpi().isMainProcessor()) {
+ fout.open(fname.c_str(), std::ios::app);
+ fout << "\n";
+ fout.close();
+ }
+}
+
+
+// TODO: Move this method to Cuboid3D class
+/// Helper Function to create cuboid parameters for XML tag
+template
+std::string CuboidGeometry3D::_cuboidParameters(Cuboid3D const& cub)
+{
+ std::stringstream ss;
+ ss.flags(std::ios::scientific);
+ ss.precision (std::numeric_limits::digits10 + 1);
+ ss << " extent=\"";
+ for (int i = 0; i<3; i++) {
+ ss << cub.getExtend()[i] << " ";
+ }
+
+ ss << "\" origin=\"";
+ for (int i = 0; i<3; i++) {
+ ss << cub.getOrigin()[i] << " ";
+ }
+
+ ss << "\" deltaR=\"" << cub.getDeltaR();
+ ss << "\" weight=\"" << cub.getWeightValue();
+ ss << "\" refinementLevel=\"" << cub.getRefinementLevel() << "\"";
+ return ss.str();
+}
+
+
+} // namespace olb
+
+#endif
--
cgit v1.2.3