diff options
Initialize at openlb-1-3
Diffstat (limited to 'src/geometry/blockGeometryStatistics3D.hh')
-rw-r--r-- | src/geometry/blockGeometryStatistics3D.hh | 1471 |
1 files changed, 1471 insertions, 0 deletions
diff --git a/src/geometry/blockGeometryStatistics3D.hh b/src/geometry/blockGeometryStatistics3D.hh new file mode 100644 index 0000000..86ed1ec --- /dev/null +++ b/src/geometry/blockGeometryStatistics3D.hh @@ -0,0 +1,1471 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2011, 2014 Mathias J. Krause, Simon Zimny + * 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 + * Representation of a statistic for a 3D geometry -- generic implementation. + */ + +#ifndef BLOCK_GEOMETRY_STATISTICS_3D_HH +#define BLOCK_GEOMETRY_STATISTICS_3D_HH + +#include <iostream> +#include <math.h> +#include <fstream> +#include <sstream> +#include <cmath> + +#include "geometry/blockGeometry3D.h" +#include "geometry/blockGeometryStatistics3D.h" + +namespace olb { + +template<typename T> +BlockGeometryStatistics3D<T>::BlockGeometryStatistics3D( BlockGeometryStructure3D<T>* blockGeometry) + : _blockGeometry(blockGeometry), clout(std::cout,"BlockGeometryStatistics3D") +{ + _statisticsUpdateNeeded=true; +} + + +template<typename T> +bool& BlockGeometryStatistics3D<T>::getStatisticsStatus() +{ + return _statisticsUpdateNeeded; +} + +template<typename T> +bool const & BlockGeometryStatistics3D<T>::getStatisticsStatus() const +{ + return _statisticsUpdateNeeded; +} + + +template<typename T> +void BlockGeometryStatistics3D<T>::update(bool verbose) +{ + + if (getStatisticsStatus() ) { + _material2n.clear(); + + _nX = _blockGeometry->getNx(); + _nY = _blockGeometry->getNy(); + _nZ = _blockGeometry->getNz(); + _h = _blockGeometry->getDeltaR(); + + for (int iX = 0; iX < _nX; ++iX) { + for (int iY = 0; iY < _nY; ++iY) { + for (int iZ = 0; iZ < _nZ; ++iZ) { + takeStatistics(iX, iY, iZ); + } + } + } + + _nMaterials=int(); + std::map<int, int>::iterator iter; + for (iter = _material2n.begin(); iter != _material2n.end(); iter++) { + _nMaterials++; + } + + if (verbose) { + clout << "updated" << std::endl; + } + getStatisticsStatus()=false; + } +} + + +template<typename T> +int BlockGeometryStatistics3D<T>::getNmaterials() +{ + update(); + return _nMaterials; +} + +template<typename T> +int BlockGeometryStatistics3D<T>::getNvoxel(int material) +{ + update(); + return _material2n[material]; +} + +template<typename T> +std::map<int, int> BlockGeometryStatistics3D<T>::getMaterial2n() +{ + update(); + return _material2n; +} + +template<typename T> +int BlockGeometryStatistics3D<T>::getNvoxel() +{ + update(); + int total = 0; + std::map<int, int>::iterator iter; + for (iter = _material2n.begin(); iter != _material2n.end(); iter++) { + if (iter->first!=0) { + total+=iter->second; + } + } + return total; +} + +template<typename T> +std::vector<int> BlockGeometryStatistics3D<T>::getMinLatticeR(int material) +{ + update(); + return _material2min[material]; +} + +template<typename T> +std::vector<int> BlockGeometryStatistics3D<T>::getMaxLatticeR(int material) +{ + update(); + return _material2max[material]; +} + +template<typename T> +std::vector<T> BlockGeometryStatistics3D<T>::getMinPhysR(int material) +{ + std::vector<T> tmp(3,T()); + _blockGeometry->getPhysR(&(tmp[0]), &(getMinLatticeR(material)[0])); + return tmp; +} + +template<typename T> +std::vector<T> BlockGeometryStatistics3D<T>::getMaxPhysR(int material) +{ + std::vector<T> tmp(3,T()); + _blockGeometry->getPhysR(&(tmp[0]), &(getMaxLatticeR(material)[0])); + return tmp; +} + +template<typename T> +std::vector<T> BlockGeometryStatistics3D<T>::getLatticeExtend(int material) +{ + update(); + std::vector<T> extend; + for (int iDim = 0; iDim < 3; iDim++) { + extend.push_back(_material2max[material][iDim] - _material2min[material][iDim]); + } + return extend; +} + +template<typename T> +std::vector<T> BlockGeometryStatistics3D<T>::getPhysExtend(int material) +{ + update(); + std::vector<T> extend; + for (int iDim = 0; iDim < 3; iDim++) { + extend.push_back(getMaxPhysR(material)[iDim] - getMinPhysR(material)[iDim]); + } + return extend; +} + +template<typename T> +std::vector<T> BlockGeometryStatistics3D<T>::getPhysRadius(int material) +{ + update(); + std::vector<T> radius; + for (int iDim=0; iDim<3; iDim++) { + radius.push_back((getMaxPhysR(material)[iDim] - getMinPhysR(material)[iDim])/2.); + } + return radius; +} + +template<typename T> +std::vector<T> BlockGeometryStatistics3D<T>::getCenterPhysR(int material) +{ + update(); + std::vector<T> center; + for (int iDim=0; iDim<3; iDim++) { + center.push_back(getMinPhysR(material)[iDim] + getPhysRadius(material)[iDim]); + } + return center; +} + +template<typename T> +std::vector<int> BlockGeometryStatistics3D<T>::getType(int iX, int iY, int iZ, bool anyNormal) +{ + + update(); + std::vector<int> discreteNormal(4, 0); + std::vector<int> discreteNormal2(4, 0); + std::vector<int> nullVector(4, 0); + + if (_blockGeometry->getMaterial(iX, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ) != 0) { + + //boundary0N and boundary 0P + if (_blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0) { + + if (_blockGeometry->getMaterial(iX + 1, iY, iZ) == 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 0; + discreteNormal[1] = -1; + discreteNormal[2] = 0; + discreteNormal[3] = 0; + } else { + discreteNormal2[0] = 0; + discreteNormal2[1] = -1; + discreteNormal2[2] = 0; + discreteNormal2[3] = 0; + } + } + + if (_blockGeometry->getMaterial(iX - 1, iY, iZ) == 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 0; + discreteNormal[1] = 1; + discreteNormal[2] = 0; + discreteNormal[3] = 0; + } else { + discreteNormal2[0] = 0; + discreteNormal2[1] = 1; + discreteNormal2[2] = 0; + discreteNormal2[3] = 0; + } + } + } + + // boundary1N and boundary1P + if (_blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0) { + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) == 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 0; + discreteNormal[1] = 0; + discreteNormal[2] = -1; + discreteNormal[3] = 0; + } else { + discreteNormal2[0] = 0; + discreteNormal2[1] = 0; + discreteNormal2[2] = -1; + discreteNormal2[3] = 0; + } + } + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ) == 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 0; + discreteNormal[1] = 0; + discreteNormal[2] = 1; + discreteNormal[3] = 0; + } else { + discreteNormal2[0] = 0; + discreteNormal2[1] = 0; + discreteNormal2[2] = 1; + discreteNormal2[3] = 0; + } + } + } + + // boundary2N and boundary2P + if (_blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0) { + + if (_blockGeometry->getMaterial(iX, iY, iZ + 1) == 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 0; + discreteNormal[1] = 0; + discreteNormal[2] = 0; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 0; + discreteNormal2[1] = 0; + discreteNormal2[2] = 0; + discreteNormal2[3] = -1; + } + } + + if (_blockGeometry->getMaterial(iX, iY, iZ - 1) == 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 0; + discreteNormal[1] = 0; + discreteNormal[2] = 0; + discreteNormal[3] = 1; + } else { + discreteNormal2[0] = 0; + discreteNormal2[1] = 0; + discreteNormal2[2] = 0; + discreteNormal2[3] = 1; + } + } + } + + // externalCornerNNN and externalCornerNPN + if (_blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ + 1) != 0) { + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX + 1, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY + 1, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY + 1, iZ + 1) == 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 1; + discreteNormal[1] = -1; + discreteNormal[2] = -1; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 1; + discreteNormal2[1] = -1; + discreteNormal2[2] = -1; + discreteNormal2[3] = -1; + } + } + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX + 1, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY - 1, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY - 1, iZ + 1) == 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 1; + discreteNormal[1] = -1; + discreteNormal[2] = 1; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 1; + discreteNormal2[1] = -1; + discreteNormal2[2] = 1; + discreteNormal2[3] = -1; + } + } + } + + // externalCornerNPP and externalCornerNNP + if (_blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ - 1) != 0) { + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX + 1, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY - 1, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY - 1, iZ - 1) == 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 1; + discreteNormal[1] = -1; + discreteNormal[2] = 1; + discreteNormal[3] = 1; + } else { + discreteNormal2[0] = 1; + discreteNormal2[1] = -1; + discreteNormal2[2] = 1; + discreteNormal2[3] = 1; + } + } + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX + 1, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY + 1, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY + 1, iZ - 1) == 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 1; + discreteNormal[1] = -1; + discreteNormal[2] = -1; + discreteNormal[3] = 1; + } + + else { + discreteNormal2[0] = 1; + discreteNormal2[1] = -1; + discreteNormal2[2] = -1; + discreteNormal2[3] = 1; + } + } + } + + // externalCornerPPP and externalCornerPNP + if (_blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ - 1) != 0) { + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY - 1, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX - 1, iY - 1, iZ - 1) == 1) { + + if (discreteNormal == nullVector) { + + discreteNormal[0] = 1; + discreteNormal[1] = 1; + discreteNormal[2] = 1; + discreteNormal[3] = 1; + } + + else { + discreteNormal2[0] = 1; + discreteNormal2[1] = 1; + discreteNormal2[2] = 1; + discreteNormal2[3] = 1; + } + } + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY + 1, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX - 1, iY + 1, iZ - 1) == 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 1; + discreteNormal[1] = 1; + discreteNormal[2] = -1; + discreteNormal[3] = 1; + } else { + discreteNormal2[0] = 1; + discreteNormal2[1] = 1; + discreteNormal2[2] = -1; + discreteNormal2[3] = 1; + } + } + } + + // externalCornerPNN and externalCornerPPN + if (_blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ + 1) != 0) { + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY + 1, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX - 1, iY + 1, iZ + 1) == 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 1; + discreteNormal[1] = 1; + discreteNormal[2] = -1; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 1; + discreteNormal2[1] = 1; + discreteNormal2[2] = -1; + discreteNormal2[3] = -1; + } + } + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY - 1, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX - 1, iY - 1, iZ + 1) == 1) { + + if (discreteNormal == nullVector) { + + discreteNormal[0] = 1; + discreteNormal[1] = 1; + discreteNormal[2] = 1; + discreteNormal[3] = -1; + } + + else { + discreteNormal2[0] = 1; + discreteNormal2[1] = 1; + discreteNormal2[2] = 1; + discreteNormal2[3] = -1; + } + } + } + + // internalCornerPPP and internalCornerPNP + if (_blockGeometry->getMaterial(iX - 1, iY, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) == 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0) { + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 2; + discreteNormal[1] = 1; + discreteNormal[2] = 1; + discreteNormal[3] = 1; + } else { + discreteNormal2[0] = 2; + discreteNormal2[1] = 1; + discreteNormal2[2] = 1; + discreteNormal2[3] = 1; + } + } + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 2; + discreteNormal[1] = 1; + discreteNormal[2] = -1; + discreteNormal[3] = 1; + } + + else { + discreteNormal2[0] = 2; + discreteNormal2[1] = 1; + discreteNormal2[2] = -1; + discreteNormal2[3] = 1; + } + } + } + + // internalCornerPNN and InternalCornerPPN + if (_blockGeometry->getMaterial(iX - 1, iY, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) == 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0) { + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 2; + discreteNormal[1] = 1; + discreteNormal[2] = -1; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 2; + discreteNormal2[1] = 1; + discreteNormal2[2] = -1; + discreteNormal2[3] = -1; + } + } + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 2; + discreteNormal[1] = 1; + discreteNormal[2] = 1; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 2; + discreteNormal2[1] = 1; + discreteNormal2[2] = 1; + discreteNormal2[3] = -1; + } + } + } + + // internalCornerNPP and internalCornerNNP + if (_blockGeometry->getMaterial(iX + 1, iY, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) == 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0) { + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 2; + discreteNormal[1] = -1; + discreteNormal[2] = 1; + discreteNormal[3] = 1; + } else { + discreteNormal2[0] = 2; + discreteNormal2[1] = -1; + discreteNormal2[2] = 1; + discreteNormal2[3] = 1; + } + } + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 2; + discreteNormal[1] = -1; + discreteNormal[2] = -1; + discreteNormal[3] = 1; + } + + else { + discreteNormal2[0] = 2; + discreteNormal2[1] = -1; + discreteNormal2[2] = -1; + discreteNormal2[3] = 1; + } + } + } + + // internalCornerNPN and internalCornerNNN + if (_blockGeometry->getMaterial(iX + 1, iY, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) == 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0) { + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + + discreteNormal[0] = 2; + discreteNormal[1] = -1; + discreteNormal[2] = 1; + discreteNormal[3] = -1; + } + + else { + discreteNormal2[0] = 2; + discreteNormal2[1] = -1; + discreteNormal2[2] = 1; + discreteNormal2[3] = -1; + } + } + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) == 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 2; + discreteNormal[1] = -1; + discreteNormal[2] = -1; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 2; + discreteNormal2[1] = -1; + discreteNormal2[2] = -1; + discreteNormal2[3] = -1; + } + } + } + + // externalEdge0PN and externalEdge0NN + if (_blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ + 1) != 1) { + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ + 1) == 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = 0; + discreteNormal[2] = 1; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = 0; + discreteNormal2[2] = 1; + discreteNormal2[3] = -1; + } + } + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ + 1) == 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = 0; + discreteNormal[2] = -1; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = 0; + discreteNormal2[2] = -1; + discreteNormal2[3] = -1; + } + } + } + + // externalEdge0NP and externalEdge0PP + if (_blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ - 1) != 1) { + + if (_blockGeometry->getMaterial(iX, iY + 1, iZ - 1) == 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = 0; + discreteNormal[2] = -1; + discreteNormal[3] = 1; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = 0; + discreteNormal2[2] = -1; + discreteNormal2[3] = 1; + } + } + + if (_blockGeometry->getMaterial(iX, iY - 1, iZ - 1) == 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = 0; + discreteNormal[2] = 1; + discreteNormal[3] = 1; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = 0; + discreteNormal2[2] = 1; + discreteNormal2[3] = 1; + } + } + } + + // externalEdge1NN and externalEdge1NP + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0) { + + if (_blockGeometry->getMaterial(iX + 1, iY, iZ + 1) == 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = -1; + discreteNormal[2] = 0; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = -1; + discreteNormal2[2] = 0; + discreteNormal2[3] = -1; + } + } + + if (_blockGeometry->getMaterial(iX - 1, iY, iZ + 1) == 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = 1; + discreteNormal[2] = 0; + discreteNormal[3] = -1; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = 1; + discreteNormal2[2] = 0; + discreteNormal2[3] = -1; + } + } + } + + // externalEdge1PN and externalEdge1PP + if (_blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0) { + + if (_blockGeometry->getMaterial(iX + 1, iY, iZ - 1) == 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = -1; + discreteNormal[2] = 0; + discreteNormal[3] = 1; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = -1; + discreteNormal2[2] = 0; + discreteNormal2[3] = 1; + } + } + + if (_blockGeometry->getMaterial(iX - 1, iY, iZ - 1) == 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY, iZ) != 0 + && _blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 1) { + + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = 1; + discreteNormal[2] = 0; + discreteNormal[3] = 1; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = 1; + discreteNormal2[2] = 0; + discreteNormal2[3] = 1; + } + } + } + + // externalEdge2NN and externalEdge2PN + if (_blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY + 1, iZ) != 0) { + + if (_blockGeometry->getMaterial(iX + 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX + 1, iY + 1, iZ) == 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = -1; + discreteNormal[2] = -1; + discreteNormal[3] = 0; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = -1; + discreteNormal2[2] = -1; + discreteNormal2[3] = 0; + } + } + + if (_blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY + 1, iZ) == 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = 1; + discreteNormal[2] = -1; + discreteNormal[3] = 0; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = 1; + discreteNormal2[2] = -1; + discreteNormal2[3] = 0; + } + } + } + + // externalEdge2PP and externalEdge2NP + if (_blockGeometry->getMaterial(iX, iY, iZ + 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ + 1) != 0 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 1 + && _blockGeometry->getMaterial(iX, iY, iZ - 1) != 0 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 1 + && _blockGeometry->getMaterial(iX, iY - 1, iZ) != 0) { + + if (_blockGeometry->getMaterial(iX - 1, iY, iZ) != 1 + && _blockGeometry->getMaterial(iX - 1, iY - 1, iZ) == 1) { + if (discreteNormal == nullVector) { + discreteNormal[0] = 3; + discreteNormal[1] = 1; + discreteNormal[2] = 1; + discreteNormal[3] = 0; + } else { + discreteNormal2[0] = 3; + discreteNormal2[1] = 1; + discreteNormal2[2] = 1; + discreteNormal2[3] = 0; + } |