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/MakeHeader | 41 + src/geometry/blockGeometry2D.cpp | 38 + src/geometry/blockGeometry2D.h | 121 +++ src/geometry/blockGeometry2D.hh | 288 ++++++ src/geometry/blockGeometry3D.cpp | 36 + src/geometry/blockGeometry3D.h | 127 +++ src/geometry/blockGeometry3D.hh | 331 +++++++ src/geometry/blockGeometryStatistics2D.cpp | 34 + src/geometry/blockGeometryStatistics2D.h | 142 +++ src/geometry/blockGeometryStatistics2D.hh | 517 ++++++++++ src/geometry/blockGeometryStatistics3D.cpp | 34 + src/geometry/blockGeometryStatistics3D.h | 150 +++ src/geometry/blockGeometryStatistics3D.hh | 1471 ++++++++++++++++++++++++++++ src/geometry/blockGeometryStructure2D.cpp | 40 + src/geometry/blockGeometryStructure2D.h | 171 ++++ src/geometry/blockGeometryStructure2D.hh | 534 ++++++++++ src/geometry/blockGeometryStructure3D.cpp | 35 + src/geometry/blockGeometryStructure3D.h | 168 ++++ src/geometry/blockGeometryStructure3D.hh | 783 +++++++++++++++ src/geometry/blockGeometryView2D.cpp | 35 + src/geometry/blockGeometryView2D.h | 109 +++ src/geometry/blockGeometryView2D.hh | 171 ++++ src/geometry/blockGeometryView3D.cpp | 35 + src/geometry/blockGeometryView3D.h | 113 +++ src/geometry/blockGeometryView3D.hh | 185 ++++ src/geometry/cuboid2D.cpp | 35 + src/geometry/cuboid2D.h | 166 ++++ src/geometry/cuboid2D.hh | 425 ++++++++ src/geometry/cuboid3D.cpp | 35 + src/geometry/cuboid3D.h | 175 ++++ src/geometry/cuboid3D.hh | 758 ++++++++++++++ src/geometry/cuboidGeometry2D.cpp | 35 + src/geometry/cuboidGeometry2D.h | 165 ++++ src/geometry/cuboidGeometry2D.hh | 677 +++++++++++++ src/geometry/cuboidGeometry3D.cpp | 35 + src/geometry/cuboidGeometry3D.h | 273 ++++++ src/geometry/cuboidGeometry3D.hh | 1233 +++++++++++++++++++++++ src/geometry/geometry2D.h | 38 + src/geometry/geometry2D.hh | 38 + src/geometry/geometry3D.h | 38 + src/geometry/geometry3D.hh | 40 + src/geometry/module.mk | 27 + src/geometry/superGeometry2D.cpp | 37 + src/geometry/superGeometry2D.h | 198 ++++ src/geometry/superGeometry2D.hh | 504 ++++++++++ src/geometry/superGeometry3D.cpp | 37 + src/geometry/superGeometry3D.h | 203 ++++ src/geometry/superGeometry3D.hh | 557 +++++++++++ src/geometry/superGeometryStatistics2D.cpp | 34 + src/geometry/superGeometryStatistics2D.h | 128 +++ src/geometry/superGeometryStatistics2D.hh | 392 ++++++++ src/geometry/superGeometryStatistics3D.cpp | 34 + src/geometry/superGeometryStatistics3D.h | 140 +++ src/geometry/superGeometryStatistics3D.hh | 467 +++++++++ 54 files changed, 12633 insertions(+) create mode 100644 src/geometry/MakeHeader create mode 100644 src/geometry/blockGeometry2D.cpp create mode 100644 src/geometry/blockGeometry2D.h create mode 100644 src/geometry/blockGeometry2D.hh create mode 100644 src/geometry/blockGeometry3D.cpp create mode 100644 src/geometry/blockGeometry3D.h create mode 100644 src/geometry/blockGeometry3D.hh create mode 100644 src/geometry/blockGeometryStatistics2D.cpp create mode 100644 src/geometry/blockGeometryStatistics2D.h create mode 100644 src/geometry/blockGeometryStatistics2D.hh create mode 100644 src/geometry/blockGeometryStatistics3D.cpp create mode 100644 src/geometry/blockGeometryStatistics3D.h create mode 100644 src/geometry/blockGeometryStatistics3D.hh create mode 100644 src/geometry/blockGeometryStructure2D.cpp create mode 100644 src/geometry/blockGeometryStructure2D.h create mode 100644 src/geometry/blockGeometryStructure2D.hh create mode 100644 src/geometry/blockGeometryStructure3D.cpp create mode 100644 src/geometry/blockGeometryStructure3D.h create mode 100644 src/geometry/blockGeometryStructure3D.hh create mode 100644 src/geometry/blockGeometryView2D.cpp create mode 100644 src/geometry/blockGeometryView2D.h create mode 100644 src/geometry/blockGeometryView2D.hh create mode 100644 src/geometry/blockGeometryView3D.cpp create mode 100644 src/geometry/blockGeometryView3D.h create mode 100644 src/geometry/blockGeometryView3D.hh create mode 100644 src/geometry/cuboid2D.cpp create mode 100644 src/geometry/cuboid2D.h create mode 100644 src/geometry/cuboid2D.hh create mode 100644 src/geometry/cuboid3D.cpp create mode 100644 src/geometry/cuboid3D.h create mode 100644 src/geometry/cuboid3D.hh create mode 100644 src/geometry/cuboidGeometry2D.cpp create mode 100644 src/geometry/cuboidGeometry2D.h create mode 100644 src/geometry/cuboidGeometry2D.hh create mode 100644 src/geometry/cuboidGeometry3D.cpp create mode 100644 src/geometry/cuboidGeometry3D.h create mode 100644 src/geometry/cuboidGeometry3D.hh create mode 100644 src/geometry/geometry2D.h create mode 100644 src/geometry/geometry2D.hh create mode 100644 src/geometry/geometry3D.h create mode 100644 src/geometry/geometry3D.hh create mode 100644 src/geometry/module.mk create mode 100644 src/geometry/superGeometry2D.cpp create mode 100644 src/geometry/superGeometry2D.h create mode 100644 src/geometry/superGeometry2D.hh create mode 100644 src/geometry/superGeometry3D.cpp create mode 100644 src/geometry/superGeometry3D.h create mode 100644 src/geometry/superGeometry3D.hh create mode 100644 src/geometry/superGeometryStatistics2D.cpp create mode 100644 src/geometry/superGeometryStatistics2D.h create mode 100644 src/geometry/superGeometryStatistics2D.hh create mode 100644 src/geometry/superGeometryStatistics3D.cpp create mode 100644 src/geometry/superGeometryStatistics3D.h create mode 100644 src/geometry/superGeometryStatistics3D.hh (limited to 'src/geometry') diff --git a/src/geometry/MakeHeader b/src/geometry/MakeHeader new file mode 100644 index 0000000..25442ca --- /dev/null +++ b/src/geometry/MakeHeader @@ -0,0 +1,41 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 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. + + +generic := + +precompiled := blockGeometry2D \ + blockGeometry3D \ + blockGeometryStatistics2D \ + blockGeometryStatistics3D \ + blockGeometryStructure2D \ + blockGeometryStructure3D \ + blockGeometryView2D \ + blockGeometryView3D \ + cuboid2D \ + cuboid3D \ + cuboidGeometry2D \ + cuboidGeometry3D \ + superGeometry2D \ + superGeometry3D \ + superGeometryStatistics2D \ + superGeometryStatistics3D diff --git a/src/geometry/blockGeometry2D.cpp b/src/geometry/blockGeometry2D.cpp new file mode 100644 index 0000000..7bcadf3 --- /dev/null +++ b/src/geometry/blockGeometry2D.cpp @@ -0,0 +1,38 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2011 Mathias J. Krause, Simon Zimny + * 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 + * Representation of the 2D block geometry -- template instantiation. + */ + +#include "dynamics/latticeDescriptors.h" + +#include "blockGeometry2D.h" +#include "blockGeometry2D.hh" + +namespace olb { + +template class BlockGeometry2D; + +} + diff --git a/src/geometry/blockGeometry2D.h b/src/geometry/blockGeometry2D.h new file mode 100644 index 0000000..fd94523 --- /dev/null +++ b/src/geometry/blockGeometry2D.h @@ -0,0 +1,121 @@ +/* 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 + * + * + * 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 the 2D block geometry view -- header file. + */ + +#ifndef BLOCK_GEOMETRY_2D_H +#define BLOCK_GEOMETRY_2D_H + +#include +#include +#include "core/blockData2D.h" +#include "geometry/blockGeometryStatistics2D.h" +#include "geometry/blockGeometryStructure2D.h" +#include "geometry/cuboid2D.h" + +/// All OpenLB code is contained in this namespace. +namespace olb { + +/// Representation of a block geometry +/** This class is derived from block geometry structure. It + * holds the actual data with the materials. It stores pointers + * to all dependent block geometry views. + * It presents a volume of voxels where different types are + * given my material numbers which is imporant e.g. to work + * with different boundaries (like for inflow/output regions). + * + */ +template +class BlockGeometry2D final : public BlockData2D, public BlockGeometryStructure2D { + +private: + /// Cuboid which charaterizes the block geometry + Cuboid2D _cuboid; + /// List to all depending statistic status objects + std::list _statisticsUpdateNeeded; + +public: + /// Constructor + BlockGeometry2D(T x0, T y0, T h, int nX, int nY, int iCglob=-1); + /// Constructor + BlockGeometry2D(Cuboid2D& cuboid, int iCglob=-1); + /// Copy constructor + BlockGeometry2D(BlockGeometry2D const& rhs); + /// Copy assignment + BlockGeometry2D& operator=(BlockGeometry2D const& rhs); + + /// Returns the underlying block structure + BlockStructure2D& getBlockStructure() override; + + /// Write access to the associated block statistic + BlockGeometryStatistics2D& getStatistics(bool verbose=true) override; + /// Read only access to the associated block statistic + BlockGeometryStatistics2D const& getStatistics(bool verbose=true) const override; + + /// Read only access to the origin position given in SI units (meter) + Vector getOrigin() const override; + /// Read only access to the voxel size given in SI units (meter) + const T getDeltaR() const override; + /// Returns the extend (number of voxels) in X-direction + int getNx() const override; + /// Returns the extend (number of voxels) in Y-direction + int getNy() const override; + + /// Write access to a material number + int& get(int iX, int iY) override; // override BlockGeometryStructure2D::get() by linking to BlockData2D::get() + /// Read only access to a material number + int const& get(int iX, int iY) const override; // override BlockGeometryStructure2D::get() by linking to BlockData2D::get() + /// returns the (iX,iY) entry in the 2D scalar field + int getMaterial(int iX, int iY) const override; // TODO old + + /// Transforms lattice to physical coordinates (wrapped from cuboid geometry) + void getPhysR(T physR[2], const int& iX, const int& iY) const override; + + // returns the raw data field + // olb::ScalarField2D* getRawData(); + // void resize(int X, int Y, int Z, int nX, int nY, int nZ); + // refine Mesh + // void refineMesh(int level = 2); + + /// Adds a pointer to the list of dependent statistic classes + void addToStatisticsList(bool* statisticStatus) override; + /// Removes a pointer from the list of dependent statistic classes if existing + void removeFromStatisticsList(bool* statisticStatus) override; + + /// Prints a chosen part of the block geometry + void printLayer(int x0, int x1, int y0, int y1, bool linenumber = false); + /// Prints a chosen part of the block geometry + void printLayer(int direction, int layer, bool linenumber = false); + /// Prints a chosen node and its neighbourhood + void printNode(int x0, int y0); + +private: + /// Resets all depending statistic flags + void resetStatistics(); +}; + +} // namespace olb + +#endif diff --git a/src/geometry/blockGeometry2D.hh b/src/geometry/blockGeometry2D.hh new file mode 100644 index 0000000..9b46cdd --- /dev/null +++ b/src/geometry/blockGeometry2D.hh @@ -0,0 +1,288 @@ +/* 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 + * + * + * 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 the 2D block geometry -- generic implementation. + */ + +#ifndef BLOCK_GEOMETRY_2D_HH +#define BLOCK_GEOMETRY_2D_HH + + +#include "geometry/blockGeometry2D.h" + +namespace olb { + +template +BlockGeometry2D::BlockGeometry2D(T x0, T y0, T h, int nX, int nY, int iCglob) + : BlockData2D(nX,nY), BlockGeometryStructure2D(iCglob), _cuboid(x0, y0, h, nX, nY) +{ + this->_statistics = BlockGeometryStatistics2D(this); + addToStatisticsList( &(this->_statistics.getStatisticsStatus()) ); +} + +template +BlockGeometry2D::BlockGeometry2D(Cuboid2D& cuboid, int iCglob) + : BlockData2D(cuboid.getNx(),cuboid.getNy()), BlockGeometryStructure2D(iCglob), + _cuboid(cuboid) +{ + this->_statistics = BlockGeometryStatistics2D(this); + addToStatisticsList( &(this->_statistics.getStatisticsStatus()) ); +} + +template +BlockGeometry2D::BlockGeometry2D(BlockGeometry2D const& rhs) + : BlockData2D(rhs.getNx(), rhs.getNy()), BlockGeometryStructure2D(rhs._iCglob), + _cuboid(rhs._cuboid) +{ + this->_statistics = BlockGeometryStatistics2D(this); + addToStatisticsList( &(this->_statistics.getStatisticsStatus()) ); +} + +template +BlockGeometry2D& BlockGeometry2D::operator=(BlockGeometry2D const& rhs) +{ + this->_nx = rhs._nx; + this->_ny = rhs._ny; + this->_iCglob = rhs._iCglob; + _cuboid = rhs._cuboid; + this->_statistics = BlockGeometryStatistics2D(this); + addToStatisticsList( &(this->_statistics.getStatisticsStatus()) ); + return *this; +} + +template +BlockStructure2D& BlockGeometry2D::getBlockStructure() +{ + return *this; +} + +template +BlockGeometryStatistics2D& BlockGeometry2D::getStatistics(bool verbose) +{ + return this->_statistics; +} + +template +BlockGeometryStatistics2D const& BlockGeometry2D::getStatistics(bool verbose) const +{ + return this->_statistics; +} + +template +Vector BlockGeometry2D::getOrigin() const +{ + return _cuboid.getOrigin(); +} + +template +const T BlockGeometry2D::getDeltaR() const +{ + return _cuboid.getDeltaR(); +} + +template +int BlockGeometry2D::getNx() const +{ + return _cuboid.getNx(); +} + +template +int BlockGeometry2D::getNy() const +{ + return _cuboid.getNy(); +} + + +template +int& BlockGeometry2D::get(int iX, int iY) +{ + resetStatistics(); + return BlockData2D::get(iX,iY); +} + +template +int const& BlockGeometry2D::get(int iX, int iY) const +{ + return BlockData2D::get(iX,iY); +} + +template +int BlockGeometry2D::getMaterial(int iX, int iY) const +{ + int material; + if (iX < 0 || iX + 1 > getNx() || iY < 0 || iY + 1 > getNy() ) { + material = 0; + } else { + material = BlockData2D::get(iX,iY); + } + return material; +} + +template +void BlockGeometry2D::getPhysR(T physR[2], const int& iX, const int& iY) const +{ + _cuboid.getPhysR(physR, iX, iY); + return; +} + +///TODO up, DONE down +/* +template +olb::ScalarField3D* BlockGeometry3D::getRawData() { + return &_geometryData; +}*/ +/* +template +void BlockGeometry3D::resize(int X, int Y, int Z, int nX, int nY, int nZ) { + olb::ScalarField3D geometryDataTemp(nX, nY, nZ); + geometryDataTemp.construct(); + geometryDataTemp.reset(); + for(int iX=0; iX +void BlockGeometry3D::refineMesh(int level) { + + // new number of Voxels in X-,Y-,and Z-direction and new spacing + + int _nXnew = level * getNx(); + int _nYnew = level * getNy(); + int _nZnew = level * getNz(); + T _hnew = getDeltaR() / ((T) level); + + olb::ScalarField3D _refinedGeometryData(_nXnew, _nYnew, + _nZnew); + _refinedGeometryData.construct(); + + for (int iX = 0; iX < getNx(); iX++) { + for (int iY = 0; iY < getNy(); iY++) { + for (int iZ = 0; iZ < getNz(); iZ++) { + for (int li = 0; li < level; li++) { + for (int lj = 0; lj < level; lj++) { + for (int lk = 0; lk < level; lk++) { + + _refinedGeometryData.get(level * iX + li, level + * iY + lj, level * iZ + lk) = getMaterial( + iX, iY, iZ); + + } + } + } + } + } + } + + for (int iX = 0; iX < _nXnew; iX++) { + for (int iY = 0; iY < _nYnew; iY++) { + for (int iZ = 0; iZ < _nZnew; iZ++) { + + if (iX == 0 || iY == 0 || iZ == 0 || iX == _nXnew - 1 || iY + == _nYnew - 1 || iZ == _nZnew - 1) + _refinedGeometryData.get(iX, iY, iZ) = 0; + } + } + } + reInit(getOrigin()[0], getOrigin()[1], getOrigin()[2], _hnew, _nXnew - 2, _nYnew - 2, _nZnew - 2, this->_iCglob, + &_refinedGeometryData); +}*/ + + +template +void BlockGeometry2D::addToStatisticsList(bool* statisticStatus) +{ + _statisticsUpdateNeeded.push_back(statisticStatus); +} + +template +void BlockGeometry2D::removeFromStatisticsList(bool* statisticStatus) +{ + _statisticsUpdateNeeded.remove(statisticStatus); +} + + +template +void BlockGeometry2D::printLayer(int x0, int x1, int y0, int y1, bool linenumber) +{ + for (int x = x0; x <= x1; x++) { + if (linenumber) { + this->clout << x << ": "; + } + for (int y = y0; y <= y1; y++) { + this->clout << getMaterial(x, y) << " "; + } + if (x1 - x0 != 0) { + this->clout << std::endl; + } + } + this->clout << std::endl; +} + +template +void BlockGeometry2D::printLayer(int direction, int layer, bool linenumber) +{ + assert(direction >= 0 && direction <= 2); + switch (direction) { + case 0: + printLayer(layer, layer, 0, getNy() - 1, linenumber); + break; + case 1: + printLayer(0, getNx() - 1, layer, layer, linenumber); + break; + } +} + +template +void BlockGeometry2D::printNode(int x0, int y0) +{ + for (int x = x0 - 1; x <= x0 + 1; x++) { + this->clout << "x=" << x << std::endl; + for (int y = y0 - 1; y <= y0 + 1; y++) { + this->clout << getMaterial(x, y) << " "; + } + this->clout << std::endl; + } + this->clout << std::endl; +} + +template +void BlockGeometry2D::resetStatistics() +{ + for (std::list::iterator it = _statisticsUpdateNeeded.begin(); it != _statisticsUpdateNeeded.end(); it++) { + **it = true; + } +} + +} // namespace olb + +#endif diff --git a/src/geometry/blockGeometry3D.cpp b/src/geometry/blockGeometry3D.cpp new file mode 100644 index 0000000..7f3ca2f --- /dev/null +++ b/src/geometry/blockGeometry3D.cpp @@ -0,0 +1,36 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2011 Mathias J. Krause, Simon Zimny + * 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 + * Representation of the 3D block geometry -- template instantiation. + */ + +#include "blockGeometry3D.h" +#include "blockGeometry3D.hh" + +namespace olb { + +template class BlockGeometry3D; + +} + diff --git a/src/geometry/blockGeometry3D.h b/src/geometry/blockGeometry3D.h new file mode 100644 index 0000000..3d0a401 --- /dev/null +++ b/src/geometry/blockGeometry3D.h @@ -0,0 +1,127 @@ +/* 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 + * + * + * 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 the 3D block geometry view -- header file. + */ + +#ifndef BLOCK_GEOMETRY_3D_H +#define BLOCK_GEOMETRY_3D_H + +#include +#include +#include "core/blockData3D.h" +#include "geometry/blockGeometryStatistics3D.h" +#include "geometry/blockGeometryStructure3D.h" +#include "geometry/cuboid3D.h" + +/// All OpenLB code is contained in this namespace. +namespace olb { + +/// Representation of a block geometry +/** This class is derived from block geometry structure. It + * holds the actual data with the materials. It stores pointers + * to all dependent block geometry views. + * It presents a volume of voxels where different types are + * given by material numbers which is important e.g. to work + * with different boundaries (like for inflow/output regions). + * + * This class is not intended to be derived from. + */ +template +class BlockGeometry3D final : public BlockData3D, public BlockGeometryStructure3D { + +private: + /// Cuboid which charaterizes the block geometry + Cuboid3D _cuboid; + /// List to all depending statistic status objects + std::list _statisticsUpdateNeeded; + +public: + /// Constructor + BlockGeometry3D(T x0, T y0, T z0, T h, int nX, int nY, int nZ, int iCglob=-1); + /// Constructor + BlockGeometry3D(Cuboid3D& cuboid, int iCglob=-1); + /// Copy constructor + BlockGeometry3D(BlockGeometry3D const& rhs); + /// Copy assignment + BlockGeometry3D& operator=(BlockGeometry3D const& rhs); + + /// Returns the underlying block structure + BlockStructure3D& getBlockStructure() override; + + /// Write access to the associated block statistic + BlockGeometryStatistics3D& getStatistics(bool verbose=true) override; + /// Read only access to the associated block statistic + BlockGeometryStatistics3D const& getStatistics(bool verbose=true) const override; + + /// Read only access to the origin position given in SI units (meter) + Vector getOrigin() const override; + /// Read only access to the voxel size given in SI units (meter) + const T getDeltaR() const override; + /// Returns the extend (number of voxels) in X-direction + int getNx() const override; + /// Returns the extend (number of voxels) in Y-direction + int getNy() const override; + /// Returns the extend (number of voxels) in Z-direction + int getNz() const override; + /// Write access to a material number + int& get(int iX, int iY, int iZ) override; // override BlockGeometryStructure3D::get() by linking to BlockData3D::get() + /// Read only access to a material number + int const& get(int iX, int iY, int iZ) const override; // override BlockGeometryStructure3D::get() by linking to BlockData3D::get() + /// returns the (iX,iY,iZ) entry in the 3D scalar field + int getMaterial(int iX, int iY, int iZ) const override; // TODO old + + /// Transforms lattice to physical coordinates (wrapped from cuboid geometry) + void getPhysR(T physR[3], const int& iX, const int& iY, const int& iZ) const override; + + // returns the raw data field + // olb::ScalarField3D* getRawData(); + // void resize(int X, int Y, int Z, int nX, int nY, int nZ); + // refine Mesh + // void refineMesh(int level = 2); + + /// Adds a pointer to the list of dependent statistic classes + void addToStatisticsList(bool* statisticStatus) override; + /// Removes a pointer from the list of dependent statistic classes if existing + void removeFromStatisticsList(bool* statisticStatus) override; + + /// Prints a chosen part of the block geometry + void printLayer(int x0, int x1, int y0, int y1, int z0, int z1, bool linenumber = false); + /// Prints a chosen part of the block geometry + void printLayer(int direction, int layer, bool linenumber = false); + /// Prints a chosen node and its neighbourhood + void printNode(int x0, int y0, int z0); + +private: + + /// Re-initialization of the block geometry + /*void reInit(T originX, T originY, T originZ, T deltaR, int nX, int nY, int nZ, int iCglob = -1, + olb::ScalarField3D* geometryData = nullptr);*/ + /// Resets all depending statistic flags + void resetStatistics(); +}; + +} // namespace olb + +#endif diff --git a/src/geometry/blockGeometry3D.hh b/src/geometry/blockGeometry3D.hh new file mode 100644 index 0000000..6264c8b --- /dev/null +++ b/src/geometry/blockGeometry3D.hh @@ -0,0 +1,331 @@ +/* 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 + * + * + * 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 the 3D block geometry -- generic implementation. + */ + +#ifndef BLOCK_GEOMETRY_3D_HH +#define BLOCK_GEOMETRY_3D_HH + + +#include "geometry/blockGeometry3D.h" + + +namespace olb { + +template +BlockGeometry3D::BlockGeometry3D(T x0, T y0, T z0, T h, int nX, int nY, int nZ, int iCglob) + : BlockData3D(nX,nY,nZ), BlockGeometryStructure3D(iCglob), + _cuboid(x0, y0, z0, h, nX, nY, nZ) +{ + this->_statistics = BlockGeometryStatistics3D(this); + addToStatisticsList( &(this->_statistics.getStatisticsStatus()) ); +} + +template +BlockGeometry3D::BlockGeometry3D(Cuboid3D& cuboid, int iCglob) + : BlockData3D(cuboid.getNx(),cuboid.getNy(),cuboid.getNz()), + BlockGeometryStructure3D(iCglob), _cuboid(cuboid) +{ + this->_statistics = BlockGeometryStatistics3D(this); + addToStatisticsList( &(this->_statistics.getStatisticsStatus()) ); +} + +template +BlockGeometry3D::BlockGeometry3D(BlockGeometry3D const& rhs) + : BlockData3D(rhs.getNx(), rhs.getNy(), rhs.getNz()), + BlockGeometryStructure3D(rhs._iCglob), _cuboid(rhs._cuboid) +{ + this->_statistics = BlockGeometryStatistics3D(this); + addToStatisticsList( &(this->_statistics.getStatisticsStatus()) ); +} + +template +BlockGeometry3D& BlockGeometry3D::operator=(BlockGeometry3D const& rhs) +{ + this->_nx = rhs._nx; + this->_ny = rhs._ny; + this->_nz = rhs._nz; + this->_iCglob = rhs._iCglob; + _cuboid = rhs._cuboid; + this->_statistics = BlockGeometryStatistics3D(this); + addToStatisticsList( &(this->_statistics.getStatisticsStatus()) ); + return *this; +} + +template +BlockStructure3D& BlockGeometry3D::getBlockStructure() +{ + return *this; +} + +template +BlockGeometryStatistics3D& BlockGeometry3D::getStatistics(bool verbose) +{ + return this->_statistics; +} + +template +BlockGeometryStatistics3D const& BlockGeometry3D::getStatistics(bool verbose) const +{ + return this->_statistics; +} + +template +Vector BlockGeometry3D::getOrigin() const +{ + return _cuboid.getOrigin(); +} + +template +const T BlockGeometry3D::getDeltaR() const +{ + return _cuboid.getDeltaR(); +} + +template +int BlockGeometry3D::getNx() const +{ + return _cuboid.getNx(); +} + +template +int BlockGeometry3D::getNy() const +{ + return _cuboid.getNy(); +} + +template +int BlockGeometry3D::getNz() const +{ + return _cuboid.getNz(); +} + + +template +int& BlockGeometry3D::get(int iX, int iY, int iZ) +{ + resetStatistics(); + return BlockData3D::get(iX,iY,iZ,0); +} + +template +int const& BlockGeometry3D::get(int iX, int iY, int iZ) const +{ + return BlockData3D::get(iX,iY,iZ,0); +} + +template +int BlockGeometry3D::getMaterial(int iX, int iY, int iZ) const +{ + int material; + if (iX < 0 || iX + 1 > getNx() || iY < 0 || iY + + 1 > getNy() || iZ < 0 || iZ + 1 > getNz()) { + material = 0; + } else { + material = BlockData3D::get(iX,iY,iZ,0); + } + return material; +} + +template +void BlockGeometry3D::getPhysR(T physR[3], const int& iX, const int& iY, + const int& iZ) const +{ + _cuboid.getPhysR(physR, iX, iY, iZ); + return; +} + + +///TODO up, DONE down +/* +template +olb::ScalarField3D* BlockGeometry3D::getRawData() { + return &_geometryData; +}*/ +/* +template +void BlockGeometry3D::resize(int X, int Y, int Z, int nX, int nY, int nZ) { + olb::ScalarField3D geometryDataTemp(nX, nY, nZ); + geometryDataTemp.construct(); + geometryDataTemp.reset(); + for(int iX=0; iX +void BlockGeometry3D::refineMesh(int level) { + + // new number of Voxels in X-,Y-,and Z-direction and new spacing + + int _nXnew = level * getNx(); + int _nYnew = level * getNy(); + int _nZnew = level * getNz(); + T _hnew = getDeltaR() / ((T) level); + + olb::ScalarField3D _refinedGeometryData(_nXnew, _nYnew, + _nZnew); + _refinedGeometryData.construct(); + + for (int iX = 0; iX < getNx(); iX++) { + for (int iY = 0; iY < getNy(); iY++) { + for (int iZ = 0; iZ < getNz(); iZ++) { + for (int li = 0; li < level; li++) { + for (int lj = 0; lj < level; lj++) { + for (int lk = 0; lk < level; lk++) { + + _refinedGeometryData.get(level * iX + li, level + * iY + lj, level * iZ + lk) = getMaterial( + iX, iY, iZ); + + } + } + } + } + } + } + + for (int iX = 0; iX < _nXnew; iX++) { + for (int iY = 0; iY < _nYnew; iY++) { + for (int iZ = 0; iZ < _nZnew; iZ++) { + + if (iX == 0 || iY == 0 || iZ == 0 || iX == _nXnew - 1 || iY + == _nYnew - 1 || iZ == _nZnew - 1) + _refinedGeometryData.get(iX, iY, iZ) = 0; + } + } + } + reInit(getOrigin()[0], getOrigin()[1], getOrigin()[2], _hnew, _nXnew - 2, _nYnew - 2, _nZnew - 2, this->_iCglob, + &_refinedGeometryData); +}*/ + + +template +void BlockGeometry3D::addToStatisticsList(bool* statisticStatus) +{ + _statisticsUpdateNeeded.push_back(statisticStatus); +} + +template +void BlockGeometry3D::removeFromStatisticsList(bool* statisticStatus) +{ + _statisticsUpdateNeeded.remove(statisticStatus); +} + + +template +void BlockGeometry3D::printLayer(int x0, int x1, int y0, int y1, int z0, + int z1, bool linenumber) +{ + for (int x = x0; x <= x1; x++) { + if (linenumber) { + this->clout << x << ": "; + } + for (int y = y0; y <= y1; y++) { + for (int z = z0; z <= z1; z++) { + this->clout << getMaterial(x, y, z) << " "; + } + if (y1 - y0 != 0 && z1 - z0 != 0) { + this->clout << std::endl; + } + } + if (x1 - x0 != 0) { + this->clout << std::endl; + } + } + this->clout << std::endl; +} + +template +void BlockGeometry3D::printLayer(int direction, int layer, bool linenumber) +{ + assert(direction >= 0 && direction <= 2); + switch (direction) { + case 0: + printLayer(layer, layer, 0, getNy() - 1, 0, getNz() - 1, linenumber); + break; + case 1: + printLayer(0, getNx() - 1, layer, layer, 0, getNz() - 1, linenumber); + break; + case 2: + printLayer(0, getNx() - 1, 0, getNy() - 1, layer, layer, linenumber); + break; + } +} + +template +void BlockGeometry3D::printNode(int x0, int y0, int z0) +{ + for (int x = x0 - 1; x <= x0 + 1; x++) { + this->clout << "x=" << x << std::endl; + for (int y = y0 - 1; y <= y0 + 1; y++) { + for (int z = z0 - 1; z <= z0 + 1; z++) { + this->clout << getMaterial(x, y, z) << " "; + } + this->clout << std::endl; + } + this->clout << std::endl; + } + this->clout << std::endl; +} + +/* +template +void BlockGeometry3D::reInit(T x0, T y0, T z0, T h, int nX, int nY, int nZ, + int iCglob, olb::ScalarField3D* geometryData) { + + _cuboid = Cuboid3D(x0, y0, z0, h, nX, nY, nZ); + if (geometryData != nullptr) { + _geometryData.swap(*geometryData); + } + else { + olb::ScalarField3D geometryDataTemp(nX, nY, nZ); + _geometryData.swap(geometryDataTemp); + _geometryData.construct(); + _geometryData.reset(); + geometryDataTemp.deConstruct(); + } + this->_iCglob=iCglob; +}*/ + +template +void BlockGeometry3D::resetStatistics() +{ + for (std::list::iterator it = _statisticsUpdateNeeded.begin(); it != _statisticsUpdateNeeded.end(); ++it) { + **it = true; + } +} + +} // namespace olb + +#endif diff --git a/src/geometry/blockGeometryStatistics2D.cpp b/src/geometry/blockGeometryStatistics2D.cpp new file mode 100644 index 0000000..9d894d1 --- /dev/null +++ b/src/geometry/blockGeometryStatistics2D.cpp @@ -0,0 +1,34 @@ +/* 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 + * + * + * 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 2D geometry -- template instantiation. + */ + +#include "blockGeometryStatistics2D.h" +#include "blockGeometryStatistics2D.hh" + +namespace olb { +template class BlockGeometryStatistics2D; +} + diff --git a/src/geometry/blockGeometryStatistics2D.h b/src/geometry/blockGeometryStatistics2D.h new file mode 100644 index 0000000..8d682ac --- /dev/null +++ b/src/geometry/blockGeometryStatistics2D.h @@ -0,0 +1,142 @@ +/* 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 + * + * + * 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 2D geometry -- header file. + */ + +#ifndef BLOCK_GEOMETRY_STATISTICS_2D_H +#define BLOCK_GEOMETRY_STATISTICS_2D_H + +#include +#include +#include + +#include "io/ostreamManager.h" + +/// All OpenLB code is contained in this namespace. +namespace olb { + +/// Representation of a statistic for a 2D geometry +/** A block geomety statistic computes different integral + * values, like total number of different materials, + * materials of any kind, min./max. physical position, of an + * underlying block geoemtry structure. + * + * This class is not intended to be derived from. + */ + +template +class BlockGeometryStructure2D; + +template +class BlockGeometryStatistics2D { + +private: + + /// Points to the underlying data from which the statistics is taken + BlockGeometryStructure2D* _blockGeometry; + /// Specifies if an update is needed + bool _statisticsUpdateNeeded; + /// Number of voxels in each direction + int _nX, _nY; + /// Spacing + T _h; + + /// Number of different material numbers + int _nMaterials; + /// Mapping a material number to the number of this kind found in the super geometry + std::map _material2n; + /// Mapping a material number to the min. lattice position in each space direction + std::map > _material2min; + /// Mapping a material number to the max. lattice position in each space direction + std::map > _material2max; + + /// class specific cout + mutable OstreamManager clout; + +public: + + /// Constructor + BlockGeometryStatistics2D(BlockGeometryStructure2D* blockGeometry); + + /// Read and write access to a flag, which indicates if an uptate is needed (=true) + bool& getStatisticsStatus(); + /// Read only access to a flag, which indicates if an uptate is needed (=true) + bool const & getStatisticsStatus() const; + /// Returns the map with the numbers of voxels for each material + std::map getMaterial2n(); + + /// Updates the statistics if it is really needed + void update(bool verbose=true); + + /// Returns the number of different materials + int getNmaterials(); + /// Returns the number of voxels for a given material number + int getNvoxel(int material); + /// Returns the number of voxels with material!=0 + int getNvoxel(); + /// Returns the min. lattice position in each direction + std::vector getMinLatticeR(int material); + /// Returns the max. lattice position in each direction + std::vector getMaxLatticeR(int material); + /// Returns the min. phys position in each direction + std::vector getMinPhysR(int material); + /// Returns the max. phys position in each direction + std::vector getMaxPhysR(int material); + /// Returns the lattice extend as length in each direction + std::vector getLatticeExtend(int material); + /// Returns the phys extend as length in each direction + std::vector getPhysExtend(int material); + /// Returns the phys radius as length in each direction + std::vector getPhysRadius(int material); + /// Returns the center position + std::vector getCenterPhysR(int material); + /// Returns the boundary type which is characterized by a discrete normal (c.f. Zimny) + std::vector getType(int iX, int iY); + + /// Returns normal that points into the fluid for paraxial surfaces + std::vector computeNormal(int iX, int iY); + /// Returns normal that points into the fluid for paraxial surfaces + std::vector computeNormal (int material); + /// Returns discrete normal with norm maxNorm that points into the fluid for paraxial surfaces + /// maxNorm=1.1 implies only normals parallel to the axises + std::vector computeDiscreteNormal (int material, T maxNorm = 1.1); + + // Returns true if at position (iX,iY,iZ) and in a neighbourhood of size (offsetX,offsetY,offsetZ) only voxels of the given material are found + bool check(int material, int iX, int iY, unsigned offsetX, unsigned offsetY); + // Returns true and a position (iX,iY,iZ) if there is a neighbourhood of size (offsetX,offsetY,offsetZ) around (iX,iY,iZ) with only voxels of the given material + bool find(int material, unsigned offsetX, unsigned offsetY, int& iX, int& iY); + + /// Prints some statistic information, i.e. the number of voxels and min. max. physical position for each different material + void print(); + +private: + + /// Helper function to simplify the implementation + void takeStatistics(int iX, int iY); +}; + +} // namespace olb + +#endif diff --git a/src/geometry/blockGeometryStatistics2D.hh b/src/geometry/blockGeometryStatistics2D.hh new file mode 100644 index 0000000..ffd56e7 --- /dev/null +++ b/src/geometry/blockGeometryStatistics2D.hh @@ -0,0 +1,517 @@ +/* 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 + * + * + * 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 2D geometry -- generic implementation. + */ + +#ifndef BLOCK_GEOMETRY_STATISTICS_2D_HH +#define BLOCK_GEOMETRY_STATISTICS_2D_HH + +#include +#include +#include +#include +#include + +#include "geometry/blockGeometry2D.h" +#include "geometry/blockGeometryStatistics2D.h" + +namespace olb { + +template +BlockGeometryStatistics2D::BlockGeometryStatistics2D( BlockGeometryStructure2D* blockGeometry) + : _blockGeometry(blockGeometry), clout(std::cout,"BlockGeometryStatistics2D") +{ + _statisticsUpdateNeeded=true; +} + + +template +bool& BlockGeometryStatistics2D::getStatisticsStatus() +{ + return _statisticsUpdateNeeded; +} + +template +bool const & BlockGeometryStatistics2D::getStatisticsStatus() const +{ + return _statisticsUpdateNeeded; +} + + +template +void BlockGeometryStatistics2D::update(bool verbose) +{ + + if (getStatisticsStatus() ) { + _material2n.clear(); + + _nX = _blockGeometry->getNx(); + _nY = _blockGeometry->getNy(); + _h = _blockGeometry->getDeltaR(); + + for (int iX = 0; iX < _nX; ++iX) { + for (int iY = 0; iY < _nY; ++iY) { + takeStatistics(iX, iY); + } + } + + _nMaterials=int(); + std::map::iterator iter; + for (iter = _material2n.begin(); iter != _material2n.end(); iter++) { + _nMaterials++; + } + + if (verbose) { + clout << "updated" << std::endl; + } + getStatisticsStatus()=false; + } +} + + +template +int BlockGeometryStatistics2D::getNmaterials() +{ + update(); + return _nMaterials; +} + +template +int BlockGeometryStatistics2D::getNvoxel(int material) +{ + update(); + return _material2n[material]; +} + +template +std::map BlockGeometryStatistics2D::getMaterial2n() +{ + update(); + return _material2n; +} + +template +int BlockGeometryStatistics2D::getNvoxel() +{ + update(); + int total = 0; + std::map::iterator iter; + for (iter = _material2n.begin(); iter != _material2n.end(); iter++) { + if (iter->first!=0) { + total+=iter->second; + } + } + return total; +} + +template +std::vector BlockGeometryStatistics2D::getMinLatticeR(int material) +{ + update(); + return _material2min[material]; +} + +template +std::vector BlockGeometryStatistics2D::getMaxLatticeR(int material) +{ + update(); + return _material2max[material]; +} + +template +std::vector BlockGeometryStatistics2D::getMinPhysR(int material) +{ + std::vector tmp(2,T()); + _blockGeometry->getPhysR(&(tmp[0]), &(getMinLatticeR(material)[0])); + return tmp; +} + +template +std::vector BlockGeometryStatistics2D::getMaxPhysR(int material) +{ + std::vector tmp(2,T()); + _blockGeometry->getPhysR(&(tmp[0]), &(getMaxLatticeR(material)[0])); + return tmp; +} + +template +std::vector BlockGeometryStatistics2D::getLatticeExtend(int material) +{ + update(); + std::vector extend; + for (int iDim = 0; iDim < 2; iDim++) { + extend.push_back(_material2max[material][iDim] - _material2min[material][iDim]); + } + return extend; +} + +template +std::vector BlockGeometryStatistics2D::getPhysExtend(int material) +{ + update(); + std::vector extend; + for (int iDim = 0; iDim < 2; iDim++) { + extend.push_back(getMaxPhysR(material)[iDim] - getMinPhysR(material)[iDim]); + } + return extend; +} + +template +std::vector BlockGeometryStatistics2D::getPhysRadius(int material) +{ + update(); + std::vector radius; + for (int iDim=0; iDim<2; iDim++) { + radius.push_back((getMaxPhysR(material)[iDim] - getMinPhysR(material)[iDim])/2.); + } + return radius; +} + +template +std::vector BlockGeometryStatistics2D::getCenterPhysR(int material) +{ + update(); + std::vector center; + for (int iDim=0; iDim<2; iDim++) { + center.push_back(getMinPhysR(material)[iDim] + getPhysRadius(material)[iDim]); + } + return center; +} + +template +std::vector BlockGeometryStatistics2D::getType(int iX, int iY) +{ + + std::vector discreteNormal(3, 0); + + if (_blockGeometry->getMaterial(iX, iY) != 1 + && _blockGeometry->getMaterial(iX, iY) != 0) { + + ///boundary0N and boundary 0P + if (_blockGeometry->getMaterial(iX, iY + 1) != 1 + && _blockGeometry->getMaterial(iX, iY + 1) != 0 + && _blockGeometry->getMaterial(iX, iY - 1) != 1 + && _blockGeometry->getMaterial(iX, iY - 1) != 0) { + + if (_blockGeometry->getMaterial(iX + 1, iY) == 1) { + discreteNormal[0] = 0; + discreteNormal[1] = -1; + discreteNormal[2] = 0; + } + + if (_blockGeometry->getMaterial(iX - 1, iY) == 1) { + discreteNormal[0] = 0; + discreteNormal[1] = 1; + discreteNormal[2] = 0; + } + } + + /// boundary1N and boundary1P + if (_blockGeometry->getMaterial(iX + 1, iY) != 1 + && _blockGeometry->getMaterial(iX + 1, iY) != 0 + && _blockGeometry->getMaterial(iX - 1, iY) != 1 + && _blockGeometry->getMaterial(iX - 1, iY) != 0) { + + if (_blockGeometry->getMaterial(iX, iY + 1) == 1) { + discreteNormal[0] = 0; + discreteNormal[1] = 0; + discreteNormal[2] = -1; + } + + if (_blockGeometry->getMaterial(iX, iY - 1) == 1) { + discreteNormal[0] = 0; + discreteNormal[1] = 0; + discreteNormal[2] = 1; + } + } + + /// externalCornerNN and externalCornerNP + if (_blockGeometry->getMaterial(iX + 1, iY) != 1 + && _blockGeometry->getMaterial(iX + 1, iY) != 0) { + + if (_blockGeometry->getMaterial(iX, iY + 1) != 1 + && _blockGeometry->getMaterial(iX, iY + 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY + 1) == 1) { + discreteNormal[0] = 1; + discreteNormal[1] = -1; + discreteNormal[2] = -1; + } + + if (_blockGeometry->getMaterial(iX, iY - 1) != 1 + && _blockGeometry->getMaterial(iX, iY - 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY - 1) == 1) { + discreteNormal[0] = 1; + discreteNormal[1] = -1; + discreteNormal[2] = 1; + } + } + + /// externalCornerPN and externalCornerPP + if (_blockGeometry->getMaterial(iX - 1, iY) != 1 + && _blockGeometry->getMaterial(iX - 1, iY) != 0) { + + if (_blockGeometry->getMaterial(iX, iY + 1) != 1 + && _blockGeometry->getMaterial(iX, iY + 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY + 1) == 1) { + discreteNormal[0] = 1; + discreteNormal[1] = 1; + discreteNormal[2] = -1; + } + + if (_blockGeometry->getMaterial(iX, iY - 1) != 1 + && _blockGeometry->getMaterial(iX, iY - 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY - 1) == 1) { + discreteNormal[0] = 1; + discreteNormal[1] = 1; + discreteNormal[2] = 1; + } + } + + /// internalCornerNN and internalCornerNP + if (_blockGeometry->getMaterial(iX - 1, iY) != 1 + && _blockGeometry->getMaterial(iX - 1, iY) != 0) { + + if (_blockGeometry->getMaterial(iX, iY - 1) != 1 + && _blockGeometry->getMaterial(iX, iY - 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY - 1) == 0) { + discreteNormal[0] = 2; + discreteNormal[1] = -1; + discreteNormal[2] = -1; + } + + if (_blockGeometry->getMaterial(iX, iY + 1) != 1 + && _blockGeometry->getMaterial(iX, iY + 1) != 0 + && _blockGeometry->getMaterial(iX - 1, iY + 1) == 0) { + discreteNormal[0] = 2; + discreteNormal[1] = -1; + discreteNormal[2] = 1; + } + } + + /// internalCornerPN and internalCornerPP + if (_blockGeometry->getMaterial(iX + 1, iY) != 1 + && _blockGeometry->getMaterial(iX + 1, iY) != 0) { + + if (_blockGeometry->getMaterial(iX, iY - 1) != 1 + && _blockGeometry->getMaterial(iX, iY - 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY - 1) == 0) { + discreteNormal[0] = 2; + discreteNormal[1] = 1; + discreteNormal[2] = -1; + } + + if (_blockGeometry->getMaterial(iX, iY + 1) != 1 + && _blockGeometry->getMaterial(iX, iY + 1) != 0 + && _blockGeometry->getMaterial(iX + 1, iY + 1) == 0) { + discreteNormal[0] = 2; + discreteNormal[1] = 1; + discreteNormal[2] = 1; + } + } + } + return discreteNormal; +} + +template +std::vector BlockGeometryStatistics2D::computeNormal(int iX, int iY) +{ + update(); + std::vector normal (2,int(0)); + + if (iX != 0) { + if (_blockGeometry->getMaterial(iX - 1, iY) == 1) { + normal[0] = -1; + } + } + if (iX != _nX - 1) { + if (_blockGeometry->getMaterial(iX + 1, iY) == 1) { + normal[0] = 1; + } + } + if (iY != 0) { + if (_blockGeometry->getMaterial(iX, iY - 1) == 1) { + normal[1] = -1; + } + } + if (iY != _nY - 1) { + if (_blockGeometry->getMaterial(iX, iY + 1) == 1) { + normal[1] = 1; + } + } + return normal; +} + +template +std::vector BlockGeometryStatistics2D::computeNormal(int material) +{ + + update(); + std::vector normal (2,int(0)); + std::vector minC = getMinLatticeR(material); + std::vector maxC = getMaxLatticeR(material); + for (int iX = minC[0]; iX<=maxC[0]; iX++) { + for (int iY = minC[1]; iY<=maxC[1]; iY++) { + if (_blockGeometry->getMaterial(iX,iY) == material) { + normal[0]+=computeNormal(iX,iY)[0]; + normal[1]+=computeNormal(iX,iY)[1]; + } + } + } + T norm = sqrt(normal[0]*normal[0]+normal[1]*normal[1]); + if (norm>0.) { + normal[0]/=norm; + normal[1]/=norm; + } + return normal; +} + +template +std::vector BlockGeometryStatistics2D::computeDiscreteNormal(int material, T maxNorm) +{ + + update(); + std::vector normal = computeNormal(material); + std::vector discreteNormal(2,int(0)); + + T smallestAngle = T(0); + for (int iX = -1; iX<=1; iX++) { + for (int iY = -1; iY<=1; iY++) { + T norm = sqrt(iX*iX+iY*iY); + if (norm>0.&& norm=smallestAngle) { + smallestAngle=angle; + discreteNormal[0] = iX; + discreteNormal[1] = iY; + } + } + } + } + return discreteNormal; +} + + +template +bool BlockGeometryStatistics2D::check(int material, int iX, int iY, + unsigned offsetX, unsigned offsetY) +{ + update(); + bool found = true; + for (int iOffsetX = -offsetX; iOffsetX <= (int) offsetX; ++iOffsetX) { + for (int iOffsetY = -offsetY; iOffsetY <= (int) offsetY; ++iOffsetY) { + if (_blockGeometry->getMaterial(iX + iOffsetX, iY + iOffsetY) != material) { + found = false; + } + } + } + return found; +} + +template +bool BlockGeometryStatistics2D::find(int material, unsigned offsetX, + unsigned offsetY, int& foundX, int& foundY) +{ + update(); + bool found = false; + for (foundX = 0; foundX < _nX; foundX++) { + for (foundY = 0; foundY < _nY; foundY++) { + found = check(material, foundX, foundY, offsetX, offsetY); + if (found) { + return found; + } + } + } + return found; +} + + +template +void BlockGeometryStatistics2D::print() +{ + + update(); + std::map::iterator iter; + for (iter = _material2n.begin(); iter != _material2n.end(); iter++) { + clout << "materialNumber=" << iter->first + << "; count=" << iter->second + << "; minLatticeR=(" << _material2min[iter->first][0] <<","<< _material2min[iter->first][1] <<")" + << "; maxLatticeR=(" << _material2max[iter->first][0] <<","<< _material2max[iter->first][1] <<")" + << std::endl; + } +} + + +template +void BlockGeometryStatistics2D::takeStatistics(int iX, int iY) +{ + + int type = _blockGeometry->getMaterial(iX, iY); + if (_material2n.count(type) == 0) { + _material2n[type] = 1; + std::vector minCo; + std::vector maxCo; + minCo.push_back(iX); + minCo.push_back(iY); + _material2min[type] = minCo; + maxCo.push_back(iX); + maxCo.push_back(iY); + _material2max[type] = maxCo; + + } else { + _material2n[type]++; + if (iX < _material2min[type][0]) { + _material2min[type][0] = iX; + } + if (iY < _material2min[type][1]) { + _material2min[type][1] = iY; + } + if (iX > _material2max[type][0]) { + _material2max[type][0] = iX; + } + if (iY > _material2max[type][1]) { + _material2max[type][1] = iY; + } + } +} + +// This function compares two discrete normals (discreteNormal, discreteNormal2) in case of a duplicate assignment of boundary types. +// The goal of this function is to combine these special boundaryVoxels to an existing one (in this case boundary or externalEdge) according to +// the x-, y- and z-values of their discrete normals. +// In the following the algorithm is declared only for the x value, but it is also valid for the y and z values. +// +// for x1 = x2, the new value of x is x1 (1) +// for x1*x2 = -1, the new value of x is 0 (2) +// for x1*x2 = 0, the new value is 0 (3) +// +// It may be possible that all three values equal 0. To avoid that the values are tested again (x²+y²+z²==0) and the loosest assumption (3) is +// redefined to. +// +// If x,y and z == 0 --> find those x,y or z which are 0 because of (3) and redefine them to the value !=0 +// +// Additionally the calculated entries are multiplied with (-1) to get the right existing boundary. + +} // namespace olb + +#endif diff --git a/src/geometry/blockGeometryStatistics3D.cpp b/src/geometry/blockGeometryStatistics3D.cpp new file mode 100644 index 0000000..5783afe --- /dev/null +++ b/src/geometry/blockGeometryStatistics3D.cpp @@ -0,0 +1,34 @@ +/* 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 + * + * + * 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 -- template instantiation. + */ + +#include "blockGeometryStatistics3D.h" +#include "blockGeometryStatistics3D.hh" + +namespace olb { +template class BlockGeometryStatistics3D; +} + diff --git a/src/geometry/blockGeometryStatistics3D.h b/src/geometry/blockGeometryStatistics3D.h new file mode 100644 index 0000000..5d4411b --- /dev/null +++ b/src/geometry/blockGeometryStatistics3D.h @@ -0,0 +1,150 @@ +/* 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 + * + * + * 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 -- header file. + */ + +#ifndef BLOCK_GEOMETRY_STATISTICS_3D_H +#define BLOCK_GEOMETRY_STATISTICS_3D_H + +#include +#include +#include + +#include "io/ostreamManager.h" + +/// All OpenLB code is contained in this namespace. +namespace olb { + +/// Representation of a statistic for a 3D geometry +/** A block geomety statistic computes different integral + * values, like total number of different materials, + * materials of any kind, min./max. physical position, of an + * underlying block geoemtry structure. + * + * This class is not intended to be derived from. + */ + +template +class BlockGeometryStructure3D; + +template +class BlockGeometryStatistics3D { + +private: + + /// Points to the underlying data from which the statistics is taken + BlockGeometryStructure3D* _blockGeometry; + /// Specifies if an update is needed + bool _statisticsUpdateNeeded; + /// Number of voxels in each direction + int _nX, _nY, _nZ; + /// Spacing + T _h; + + /// Number of different material numbers + int _nMaterials; + /// Mapping a material number to the number of this kind found in the super geometry + std::map _material2n; + /// Mapping a material number to the min. lattice position in each space direction + std::map > _material2min; + /// Mapping a material number to the max. lattice positio