summaryrefslogtreecommitdiff
path: root/src/geometry/blockGeometryStatistics3D.hh
diff options
context:
space:
mode:
authorAdrian Kummerlaender2019-06-24 14:43:36 +0200
committerAdrian Kummerlaender2019-06-24 14:43:36 +0200
commit94d3e79a8617f88dc0219cfdeedfa3147833719d (patch)
treec1a6894679563e271f5c6ea7a17fa3462f7212a3 /src/geometry/blockGeometryStatistics3D.hh
downloadgrid_refinement_openlb-94d3e79a8617f88dc0219cfdeedfa3147833719d.tar
grid_refinement_openlb-94d3e79a8617f88dc0219cfdeedfa3147833719d.tar.gz
grid_refinement_openlb-94d3e79a8617f88dc0219cfdeedfa3147833719d.tar.bz2
grid_refinement_openlb-94d3e79a8617f88dc0219cfdeedfa3147833719d.tar.lz
grid_refinement_openlb-94d3e79a8617f88dc0219cfdeedfa3147833719d.tar.xz
grid_refinement_openlb-94d3e79a8617f88dc0219cfdeedfa3147833719d.tar.zst
grid_refinement_openlb-94d3e79a8617f88dc0219cfdeedfa3147833719d.zip
Initialize at openlb-1-3
Diffstat (limited to 'src/geometry/blockGeometryStatistics3D.hh')
-rw-r--r--src/geometry/blockGeometryStatistics3D.hh1471
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;
+ }