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/functors/analytical/indicator/indicatorF2D.hh | 256 ++++++++++++++++++++++ 1 file changed, 256 insertions(+) create mode 100644 src/functors/analytical/indicator/indicatorF2D.hh (limited to 'src/functors/analytical/indicator/indicatorF2D.hh') diff --git a/src/functors/analytical/indicator/indicatorF2D.hh b/src/functors/analytical/indicator/indicatorF2D.hh new file mode 100644 index 0000000..c244cb1 --- /dev/null +++ b/src/functors/analytical/indicator/indicatorF2D.hh @@ -0,0 +1,256 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2016 Thomas Henn, Cyril Masquelier, Jan Marquardt, 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. +*/ + +#ifndef INDICATOR_F_2D_HH +#define INDICATOR_F_2D_HH + +#include + +#include "indicatorF2D.h" +#include "utilities/vectorHelpers.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif +#define M_PI2 1.57079632679489661923 + +namespace olb { + +// Warning : the cuboid is only defined parallel to the plans x=0 and y=0 !!! +// For other cuboids, please use the parallelepiped version +template +IndicatorF2DfromIndicatorF3D::IndicatorF2DfromIndicatorF3D(IndicatorF3D& indicator3D) + : _indicator3D(indicator3D) +{ + this->_myMin[0] = _indicator3D.getMin()[0]; + this->_myMin[1] = _indicator3D.getMin()[1]; + this->_myMax[0] = _indicator3D.getMax()[0]; + this->_myMax[1] = _indicator3D.getMax()[1]; +} + +template +bool IndicatorF2DfromIndicatorF3D::operator()(bool output[], const S input[]) +{ + S input3D[3]; + input3D[0] = input[0]; + input3D[1] = input[1]; + input3D[2] = (_indicator3D.getMax()[2] - _indicator3D.getMin()[2]) * 0.5 + _indicator3D.getMin()[2]; + _indicator3D(output, input3D); + return true; +} + + + +// Warning : the cuboid is only defined parallel to the plans x=0 and y=0 !!! +// For other cuboids, please use the parallelepiped version +template +IndicatorCuboid2D::IndicatorCuboid2D(Vector extend, Vector origin, S theta) + : _theta(theta) +{ + this->_myMin = origin; + this->_myMax = origin + extend; + _center = origin + S(.5)*extend; + _xLength = extend[0]; + _yLength = extend[1]; +} + +template +IndicatorCuboid2D::IndicatorCuboid2D(S xLength, S yLength, Vector center, S theta ) + : _center(center), _xLength(xLength), _yLength(yLength), _theta(-theta) +{ + this->_myMin = {_center[0] - _xLength/2., _center[1] - _yLength/2.}; + this->_myMax = {_center[0] + _xLength/2., _center[1] + _yLength/2.}; +} + + +// returns true if x is inside the cuboid +template +bool IndicatorCuboid2D::operator()(bool output[], const S input[]) +{ + S x, y; + if ( !util::nearZero(_theta) ) { + x = _center[0] + (input[0] - _center[0])*std::cos(_theta) - (input[1] - _center[1])*std::sin(_theta); + y = _center[1] + (input[0] - _center[0])*std::sin(_theta) + (input[1] - _center[1])*std::cos(_theta); + } else { + x = input[0]; + y = input[1]; + } + + output[0] = ( (fabs(_center[0] - x) < _xLength/2. || util::approxEqual(fabs(_center[0] - x),_xLength/2.) ) + && (fabs(_center[1] - y) < _yLength/2. || util::approxEqual(fabs(_center[1] - y), _yLength/2.) ) ); + return true; +} + +// creator function +template +IndicatorCuboid2D* createIndicatorCuboid2D(XMLreader const& params, bool verbose) +{ + OstreamManager clout(std::cout,"createIndicatorCuboid2D"); + params.setWarningsOn(verbose); + + Vector center; + S xLength; + S yLength; + + std::stringstream xmlCenter( params.getAttribute("center") ); + xmlCenter >> center[0] >> center[1]; + std::stringstream xmlRadius( params.getAttribute("length") ); + xmlRadius >> xLength >> yLength; + + return new IndicatorCuboid2D(xLength, yLength, center); +} + + +template +IndicatorCircle2D::IndicatorCircle2D(Vector center, S radius) + : _center(center), + _radius2(radius*radius) +{ + this->_myMin = _center - radius; + this->_myMax = _center + radius; +} + +// returns true if x is inside the circle +template +bool IndicatorCircle2D::operator()(bool output[], const S input[]) +{ + output[0] = ( std::pow(_center[0] - input[0],2) + std::pow(_center[1] - input[1], 2) <= _radius2 ); + return output[0]; +} + + +template +bool IndicatorCircle2D::distance(S& distance, const Vector& origin, const Vector& direction, int iC) +{ + S a = direction[0]*direction[0] + direction[1]*direction[1]; + + // returns 0 if point is at the boundary of the sphere + if ( util::approxEqual(a,_radius2) ) { + distance = S(0); + return true; + } + // norm of direction + a = sqrt(a); + + S b = 2.*((origin[0] - _center[0])*direction[0] + + (origin[1] - _center[1])*direction[1])/a; + S c = -_radius2 + (origin[0] - _center[0])*(origin[0] - _center[0]) + + (origin[1] - _center[1])*(origin[1] - _center[1]); + + // discriminant + S d = b*b - 4.*c; + if (d < 0) { + return false; + } + + S x1 = (- b + sqrt(d)) *0.5; + S x2 = (- b - sqrt(d)) *0.5; + + // case if origin is inside the sphere + if ((x1<0.) || (x2<0.)) { + if (x1>0.) { + distance = x1; + return true; + } + if (x2>0.) { + distance = x2; + return true; + } + } + // case if origin is ouside the sphere + else { + distance = std::min(x1,x2); + return true; + } + + return false; +} + + +template +bool IndicatorCircle2D::normal(Vector& normal, const Vector& origin, const Vector& direction, int iC) +{ + S dist; + if (!(distance(dist, origin, direction, iC)) ) { + return false; + } + + Vector intresection(origin + dist*direction); + + normal = intresection - _center; + + return true; +} + + + +template +IndicatorCircle2D* createIndicatorCircle2D(XMLreader const& params, bool verbose) +{ + OstreamManager clout(std::cout,"createIndicatorCircle2D"); + params.setWarningsOn(verbose); + + Vector center; + S radius = 1; + + std::stringstream xmlCenter( params.getAttribute("center") ); + xmlCenter >> center[0] >> center[1]; + std::stringstream xmlRadius( params.getAttribute("radius") ); + xmlRadius >> radius; + + return new IndicatorCircle2D(center, radius); +} + +template +IndicatorLayer2D::IndicatorLayer2D(IndicatorF2D& indicatorF, S layerSize) + : _indicatorF(indicatorF), _layerSize(layerSize) +{ + this->_myMin = indicatorF.getMin() - layerSize; + this->_myMax = indicatorF.getMax() + layerSize; + OLB_ASSERT( (this->_myMax[0]-this->_myMin[0]) <= std::numeric_limits::epsilon() ,"Indicator reduced to zero-set in x direction"); + OLB_ASSERT( (this->_myMax[1]-this->_myMin[1]) <= std::numeric_limits::epsilon() ,"Indicator reduced to zero-set in y direction"); + _isPositive = std::signbit(layerSize); +} + +// returns true if x is inside the layer +template +bool IndicatorLayer2D::operator()(bool output[], const S input[]) +{ + output[0] = !_isPositive; + S r[2]; + for (int iX =- 1; iX < 2; ++iX) { + for (int iY =- 1; iY < 2; ++iY) { + r[0] = input[0] + iX*_layerSize; + r[1] = input[1] + iY*_layerSize; + _indicatorF(output,r); + if (output[0] == !_isPositive) { + return true; + } + } + } + return true; +} + +} // namespace olb + +#endif -- cgit v1.2.3