path: root/src/geometry/blockGeometryStructure3D.hh

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2014 Mathias J. Krause
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
 */

/** \file
 * Representation of a 3d block geometry structure -- generic implementation.
 */

#ifndef BLOCK_GEOMETRY_STRUCTURE_3D_HH
#define BLOCK_GEOMETRY_STRUCTURE_3D_HH

#include <math.h>
#include "geometry/blockGeometryStructure3D.h"
#include "functors/analytical/indicator/indicatorF3D.h"


namespace olb {

template<typename T>
BlockGeometryStructure3D<T>::BlockGeometryStructure3D(int iCglob)
  : _iCglob(iCglob), _statistics(this), clout(std::cout,"BlockGeometryStructure3D") { }

template<typename T>
int const& BlockGeometryStructure3D<T>::getIcGlob() const
{
  return _iCglob;
}

template<typename T>
Vector<int,3> const BlockGeometryStructure3D<T>::getExtend() const
{
  return Vector<int,3> (getNx(), getNy(), getNz());
}

template<typename T>
void BlockGeometryStructure3D<T>::getPhysR(T physR[3], const int latticeR[3]) const
{
  getPhysR(physR, latticeR[0], latticeR[1], latticeR[2]);
  return;
}

template<typename T>
bool BlockGeometryStructure3D<T>::isInside(int iX, int iY, int iZ) const
{
  return 0 <= iX && iX < getNx() &&
         0 <= iY && iY < getNy() &&
         0 <= iZ && iZ < getNz();
}

template<typename T>
int&  BlockGeometryStructure3D<T>::get(std::vector<int> latticeR)
{
  return get(latticeR[0], latticeR[1], latticeR[2]);
}

template<typename T>
int const&  BlockGeometryStructure3D<T>::get(std::vector<int> latticeR) const
{
  return get(latticeR[0], latticeR[1], latticeR[2]);
}

template<typename T>
int BlockGeometryStructure3D<T>::clean(bool verbose)
{
  int counter=0;
  for (int iX = 0; iX < getNx(); iX++) {
    for (int iY = 0; iY < getNy(); iY++) {
      for (int iZ = 0; iZ < getNz(); iZ++) {

        if (get(iX, iY, iZ) != 1 && get(iX, iY, iZ) != 0) {

          if (
            getMaterial(iX - 1, iY, iZ) != 1
            && getMaterial(iX, iY - 1, iZ) != 1
            && getMaterial(iX, iY, iZ - 1) != 1
            && getMaterial(iX - 1, iY - 1, iZ) != 1
            && getMaterial(iX, iY - 1, iZ - 1) != 1
            && getMaterial(iX - 1, iY, iZ - 1) != 1
            && getMaterial(iX - 1, iY - 1, iZ - 1) != 1
            && getMaterial(iX + 1, iY, iZ) != 1
            && getMaterial(iX, iY + 1, iZ) != 1
            && getMaterial(iX, iY, iZ + 1) != 1
            && getMaterial(iX + 1, iY + 1, iZ) != 1
            && getMaterial(iX, iY + 1, iZ + 1) != 1
            && getMaterial(iX + 1, iY, iZ + 1) != 1
            && getMaterial(iX + 1, iY + 1, iZ + 1) != 1
            && getMaterial(iX - 1, iY + 1, iZ) != 1
            && getMaterial(iX + 1, iY - 1, iZ) != 1
            && getMaterial(iX, iY - 1, iZ + 1) != 1
            && getMaterial(iX, iY + 1, iZ - 1) != 1
            && getMaterial(iX - 1, iY, iZ + 1) != 1
            && getMaterial(iX + 1, iY, iZ - 1) != 1
            && getMaterial(iX + 1, iY + 1, iZ - 1) != 1
            && getMaterial(iX + 1, iY - 1, iZ - 1) != 1
            && getMaterial(iX + 1, iY - 1, iZ + 1) != 1
            && getMaterial(iX - 1, iY + 1, iZ + 1) != 1
            && getMaterial(iX - 1, iY - 1, iZ +