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/blockGeometryStructure2D.hh | 534 +++++++++++++++++++++++++++++++
1 file changed, 534 insertions(+)
create mode 100644 src/geometry/blockGeometryStructure2D.hh
(limited to 'src/geometry/blockGeometryStructure2D.hh')
diff --git a/src/geometry/blockGeometryStructure2D.hh b/src/geometry/blockGeometryStructure2D.hh
new file mode 100644
index 0000000..8ef4c0d
--- /dev/null
+++ b/src/geometry/blockGeometryStructure2D.hh
@@ -0,0 +1,534 @@
+/* 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 a 2d block geometry structure -- generic implementation.
+ */
+
+#ifndef BLOCK_GEOMETRY_STRUCTURE_2D_HH
+#define BLOCK_GEOMETRY_STRUCTURE_2D_HH
+
+#include
+#include "geometry/blockGeometryStructure2D.h"
+#include "functors/lattice/indicator/blockIndicatorBaseF2D.h"
+
+
+namespace olb {
+
+template
+BlockGeometryStructure2D::BlockGeometryStructure2D(int iCglob)
+ : _iCglob(iCglob), _statistics(this), clout(std::cout,"BlockGeometryStructure2D") { }
+
+template
+int const& BlockGeometryStructure2D::getIcGlob() const
+{
+ return _iCglob;
+}
+
+template
+Vector const BlockGeometryStructure2D::getExtend() const
+{
+ return Vector (getNx(), getNy());
+}
+
+template
+void BlockGeometryStructure2D::getPhysR(T physR[2], const int latticeR[2]) const
+{
+ getPhysR(physR, latticeR[0], latticeR[1]);
+ return;
+}
+
+template
+int& BlockGeometryStructure2D::get(std::vector latticeR)
+{
+ return get(latticeR[0], latticeR[1]);
+}
+
+template
+int const& BlockGeometryStructure2D::get(std::vector latticeR) const
+{
+ return get(latticeR[0], latticeR[1]);
+}
+
+template
+int BlockGeometryStructure2D::clean(bool verbose)
+{
+ int counter=0;
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) != 1 && get(iX, iY)!= 0) {
+ if ( getMaterial(iX, iY) != 1
+ && getMaterial(iX + 1, iY) != 1
+ && getMaterial(iX - 1, iY) != 1
+ && getMaterial(iX, iY + 1) != 1
+ && getMaterial(iX + 1, iY + 1) != 1
+ && getMaterial(iX - 1, iY + 1) != 1
+ && getMaterial(iX, iY - 1) != 1
+ && getMaterial(iX + 1, iY - 1) != 1
+ && getMaterial(iX - 1, iY - 1) != 1 ) {
+ get(iX, iY) = 0;
+ counter++;
+ }
+ }
+ }
+ }
+ if (verbose) {
+ clout << "cleaned "<< counter << " outer boundary voxel(s)" << std::endl;
+ }
+ return counter;
+}
+
+template
+int BlockGeometryStructure2D::outerClean(bool verbose)
+{
+ int counter=0;
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) == 1) {
+ if ( getMaterial(iX + 1, iY) == 0
+ || getMaterial(iX - 1, iY) == 0
+ || getMaterial(iX, iY + 1) == 0
+ || getMaterial(iX + 1, iY + 1) == 0
+ || getMaterial(iX - 1, iY + 1) == 0
+ || getMaterial(iX, iY - 1) == 0
+ || getMaterial(iX + 1, iY - 1) == 0
+ || getMaterial(iX - 1, iY - 1) == 0 ) {
+ get(iX, iY) = 0;
+ counter++;
+ }
+ }
+ }
+ }
+ if (verbose) {
+ clout << "cleaned "<< counter << " outer fluid voxel(s)" << std::endl;
+ }
+ return counter;
+}
+
+template
+int BlockGeometryStructure2D::innerClean(bool verbose)
+{
+ int count = 0;
+ int count2 = 0;
+
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) != 1 && get(iX, iY) != 0) {
+ count++;
+
+ if (getMaterial(iX - 1, iY) == 1) {
+ if (getMaterial(iX, iY + 1) == 1) {
+ if (getMaterial(iX, iY - 1) == 1) {
+ get(iX, iY) = 1;
+ count2++;
+ }
+ }
+ }
+ if (getMaterial(iX + 1, iY) == 1) {
+ if (getMaterial(iX, iY + 1) == 1) {
+ if (getMaterial(iX, iY - 1) == 1) {
+ get(iX, iY) = 1;
+ count2++;
+ }
+ }
+ }
+ if (getMaterial(iX, iY + 1) == 1) {
+ if (getMaterial(iX + 1, iY) == 1) {
+ if (getMaterial(iX - 1, iY) == 1) {
+ get(iX, iY) = 1;
+ count2++;
+ }
+ }
+ }
+ if (getMaterial(iX, iY - 1) == 1) {
+ if (getMaterial(iX + 1, iY) == 1) {
+ if (getMaterial(iX - 1, iY) == 1) {
+ get(iX, iY) = 1;
+ count2++;
+ }
+ }
+ }
+ }
+ }
+ }
+ if (verbose) {
+ this->clout << "cleaned "<< count2 << " inner boundary voxel(s)" << std::endl;
+ }
+ return count2;
+}
+
+template
+int BlockGeometryStructure2D::innerClean(int fromM, bool verbose)
+{
+ int count = 0;
+ int count2 = 0;
+
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) != 1 && get(iX, iY)!= 0 && get(iX, iY) == fromM) {
+ count++;
+ if (getMaterial(iX - 1, iY) == 1) {
+ if (getMaterial(iX, iY + 1) == 1) {
+ if (getMaterial(iX, iY - 1) == 1) {
+ get(iX, iY) = 1;
+ count2++;
+ }
+ }
+ }
+ if (getMaterial(iX + 1, iY) == 1) {
+ if (getMaterial(iX, iY + 1) == 1) {
+ if (getMaterial(iX, iY - 1) == 1) {
+ get(iX, iY) = 1;
+ count2++;
+ }
+ }
+ }
+ if (getMaterial(iX, iY + 1) == 1) {
+ if (getMaterial(iX + 1, iY) == 1) {
+ if (getMaterial(iX - 1, iY) == 1) {
+ get(iX, iY) = 1;
+ count2++;
+ }
+ }
+ }
+ if (getMaterial(iX, iY - 1) == 1) {
+ if (getMaterial(iX + 1, iY) == 1) {
+ if (getMaterial(iX - 1, iY) == 1) {
+ get(iX, iY) = 1;
+ count2++;
+ }
+ }
+ }
+ }
+ }
+ }
+ if (verbose)
+ this->clout << "cleaned "<< count2
+ << " inner boundary voxel(s) of Type " << fromM << std::endl;
+ return count2;
+}
+
+template
+void BlockGeometryStructure2D::reset(IndicatorF2D& domain)
+{
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ T physR[2] { };
+ getPhysR(physR, iX, iY);
+ if (domain(physR)) {
+ get(iX, iY) = 0;
+ }
+ }
+ }
+}
+
+template
+template
+bool BlockGeometryStructure2D::findStreamDirections(
+ int iX, int iY,
+ BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator,
+ bool streamDirections[])
+{
+ if (boundaryIndicator(iX, iY)) {
+ bool found = false;
+ streamDirections[0] = false;
+ for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) {
+ streamDirections[iPop] = false;
+ if (bulkIndicator(iX + descriptors::c(iPop,0), iY + descriptors::c(iPop,1))) {
+ streamDirections[iPop] = true;
+ found = true;
+ }
+ }
+ return found;
+ }
+ else {
+ return false;
+ }
+}
+
+template
+template
+bool BlockGeometryStructure2D::findStreamDirections(int iX, int iY, int material, std::list bulkMaterials, bool streamDirections[])
+{
+ bool found = false;
+ if (getMaterial(iX, iY) != material) {
+ return false;
+ }
+ else {
+ std::list::iterator mat;
+ streamDirections[0] = false;
+ for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) {
+ streamDirections[iPop] = false;
+ for (mat=bulkMaterials.begin(); !streamDirections[iPop] && mat!=bulkMaterials.end(); ++mat) {
+ if (getMaterial(iX + descriptors::c(iPop,0), iY + descriptors::c(iPop,1)) == *mat ) {
+ streamDirections[iPop] = true;
+ found = true;
+ }
+ }
+ }
+ return found;
+ }
+}
+
+template
+bool BlockGeometryStructure2D::find(int material, unsigned offsetX, unsigned offsetY,
+ int& foundX, int& foundY)
+{
+
+ bool found = false;
+ for (foundX = 0; foundX < getNx(); foundX++) {
+ for (foundY = 0; foundY < getNy(); foundY++) {
+ found = check(material, foundX, foundY, offsetX, offsetY);
+ if (found) {
+ return found;
+ }
+ }
+ }
+ return found;
+}
+
+template
+bool BlockGeometryStructure2D::check(int material, int iX, int iY,
+ unsigned offsetX, unsigned offsetY)
+{
+ bool found = true;
+ for (int iOffsetX = -offsetX; iOffsetX <= (int) offsetX; ++iOffsetX) {
+ for (int iOffsetY = -offsetY; iOffsetY <= (int) offsetY; ++iOffsetY) {
+ if (getMaterial(iX + iOffsetX, iY + iOffsetY) != material) {
+ found = false;
+ }
+ }
+ }
+ return found;
+}
+
+template
+bool BlockGeometryStructure2D::checkForErrors(bool verbose) const
+{
+ bool error = false;
+
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) == 0) {
+ if ( getMaterial(iX + 1, iY) == 1
+ || getMaterial(iX - 1, iY) == 1
+ || getMaterial(iX, iY + 1) == 1
+ || getMaterial(iX + 1, iY + 1) == 1
+ || getMaterial(iX - 1, iY + 1) == 1
+ || getMaterial(iX, iY - 1) == 1
+ || getMaterial(iX + 1, iY - 1) == 1
+ || getMaterial(iX - 1, iY - 1) == 1 ) {
+ error = true;
+ }
+ }
+ }
+ }
+ if (verbose) {
+ if (error) {
+ this->clout << "error!" << std::endl;
+ }
+ else {
+ this->clout << "the model is correct!" << std::endl;
+ }
+ }
+ return error;
+}
+
+template
+void BlockGeometryStructure2D::rename(int fromM, int toM)
+{
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) == fromM) {
+ get(iX, iY) = toM;
+ }
+ }
+ }
+}
+
+template
+void BlockGeometryStructure2D::rename(int fromM, int toM, IndicatorF2D& condition)
+{
+ T physR[2];
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) == fromM) {
+ getPhysR(physR, iX,iY);
+ bool inside[1];
+ condition(inside, physR);
+ if (inside[0]) {
+ get(iX, iY) = toM;
+ }
+ }
+ }
+ }
+}
+
+template
+void BlockGeometryStructure2D::rename(int fromM, int toM, unsigned offsetX,
+ unsigned offsetY)
+{
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) == fromM) {
+ bool found = true;
+ for (int iOffsetX = -offsetX; iOffsetX <= (int) offsetX; ++iOffsetX) {
+ for (int iOffsetY = -offsetY; iOffsetY <= (int) offsetY; ++iOffsetY) {
+ if (getMaterial(iX + iOffsetX, iY + iOffsetY) != fromM) {
+ if (getMaterial(iX + iOffsetX, iY + iOffsetY) != 1245) {
+ found = false;
+ }
+ }
+ }
+ }
+ if (found) {
+ get(iX, iY) = 1245;
+ }
+ }
+ }
+ }
+ rename(1245,toM);
+}
+
+template
+void BlockGeometryStructure2D::rename(int fromM, int toM, int testM,
+ std::vector testDirection)
+{
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) == fromM) {
+
+ // flag that indicates the renaming of the current voxel, valid voxels are not renamed
+ bool isValid = true;
+ for (int iOffsetX = std::min(testDirection[0],0); iOffsetX <= std::max(testDirection[0],0); ++iOffsetX) {
+ for (int iOffsetY = std::min(testDirection[1],0); iOffsetY <= std::max(testDirection[1],0); ++iOffsetY) {
+ if (iOffsetX!=0 || iOffsetY!=0) {
+ if (getMaterial(iX + iOffsetX, iY + iOffsetY) != testM) {
+ isValid = false;
+ }
+ }
+ }
+ }
+ if (!isValid) {
+ get(iX, iY) = toM;
+ }
+ }
+ }
+ }
+}
+
+
+template
+void BlockGeometryStructure2D::rename(int fromM, int toM, int fluidM,
+ IndicatorF2D& condition, std::vector discreteNormal)
+{
+ rename(fromM, toM, condition);
+ std::vector testDirection(discreteNormal);
+ T physR[2];
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) == toM) {
+ getPhysR(physR, iX,iY);
+ bool inside[1];
+ condition(inside, physR);
+ if (inside[0]) {
+ if (getMaterial(iX+testDirection[0],iY+testDirection[1])!=fluidM ||
+ getMaterial(iX+2*testDirection[0],iY+2*testDirection[1])!=fluidM ||
+ getMaterial(iX-testDirection[0],iY-testDirection[1])!=0 ) {
+ get(iX, iY) = fromM;
+ }
+ }
+ }
+ }
+ }
+}
+
+template
+void BlockGeometryStructure2D::rename(int fromM, int toM, int fluidM,
+ IndicatorF2D& condition)
+{
+ rename(fromM, toM, condition);
+ std::vector testDirection = getStatistics().computeDiscreteNormal(toM);
+ T physR[3];
+ for (int iX = 0; iX < getNx(); iX++) {
+ for (int iY = 0; iY < getNy(); iY++) {
+ if (get(iX, iY) == toM) {
+ getPhysR(physR, iX,iY);
+ bool inside[1];
+ condition(inside, physR);
+ if (inside[0]) {
+ if (getMaterial(iX+testDirection[0],iY+testDirection[1])!=fluidM ||
+ getMaterial(iX+2*testDirection[0],iY+2*testDirection[1])!=fluidM ||
+ getMaterial(iX-testDirection[0],iY-testDirection[1])!=0 ) {
+ get(iX, iY) = fromM;
+ }
+ }
+ }
+ }
+ }
+}
+
+template
+void BlockGeometryStructure2D::regionGrowing(int fromM, int toM, int seedX, int seedY,
+ int offsetX, int offsetY, std::map, int>* tmp)
+{
+ std::map, int> tmp2;
+ bool firstCall = false;
+ if (tmp == nullptr) {
+ tmp = &tmp2;
+ firstCall = true;
+ }
+
+ if (getMaterial(seedX, seedY) == fromM) {
+ std::vector found;
+ found.push_back(seedX);
+ found.push_back(seedY);
+ if (tmp->count(found) == 0) {
+ (*tmp)[found] = 2;
+ if (offsetX != 0) {
+ regionGrowing(fromM, toM, seedX + 1, seedY, offsetX,
+ offsetY, tmp);
+ regionGrowing(fromM, toM, seedX - 1, seedY, offsetX,
+ offsetY, tmp);
+ }
+ if (offsetY != 0) {
+ regionGrowing(fromM, toM, seedX, seedY + 1, offsetX,
+ offsetY, tmp);
+ regionGrowing(fromM, toM, seedX, seedY - 1, offsetX,
+ offsetY, tmp);
+ }
+ }
+ }
+ if (firstCall) {
+ std::map, int>::iterator iter;
+ for (iter = tmp->begin(); iter != tmp->end(); iter++) {
+ get((iter->first)[0],(iter->first)[1]) = toM;
+ }
+ }
+ return;
+}
+
+
+} // namespace olb
+
+#endif
--
cgit v1.2.3