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 --- .../analytical/indicator/indicatorBaseF3D.hh | 321 +++++++++++++++++++++ 1 file changed, 321 insertions(+) create mode 100644 src/functors/analytical/indicator/indicatorBaseF3D.hh (limited to 'src/functors/analytical/indicator/indicatorBaseF3D.hh') diff --git a/src/functors/analytical/indicator/indicatorBaseF3D.hh b/src/functors/analytical/indicator/indicatorBaseF3D.hh new file mode 100644 index 0000000..bc10532 --- /dev/null +++ b/src/functors/analytical/indicator/indicatorBaseF3D.hh @@ -0,0 +1,321 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2014-2016 Cyril Masquelier, Albert Mink, Mathias J. Krause, Benjamin Förster + * 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_BASE_F_3D_HH +#define INDICATOR_BASE_F_3D_HH + + +#include +#include "indicatorBaseF3D.h" +#include "utilities/vectorHelpers.h" +#include "math.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif +#define M_PI2 1.57079632679489661923 + +namespace olb { + +template +IndicatorF3D::IndicatorF3D() : GenericF(1, 3) +{} + +template +Vector& IndicatorF3D::getMin() +{ + return _myMin; +} + +template +Vector& IndicatorF3D::getMax() +{ + return _myMax; +} + +template +bool IndicatorF3D::distance(S& distance, const Vector& origin, const Vector& direction, int iC) +{ + bool originValue; + bool currentValue; + S precision = .0001; + S pitch = 0.5; + + // start at origin and move into given direction + Vector currentPoint(origin); + + (*this)(&originValue, origin.data); + (*this)(¤tValue, currentPoint.data); + + while (currentValue == originValue && isInsideBox(currentPoint)) { + currentPoint += direction; + // update currentValue until the first point on the other side (inside/outside) is found + (*this)(¤tValue, currentPoint.data); + } + + // return false if no point was found in given direction + if (!isInsideBox(currentPoint) && !originValue) { + return false; + } + + + while (pitch >= precision) { + if (!isInsideBox(currentPoint) && originValue) { + currentPoint -= pitch * direction; + pitch /= 2.; + } else { + (*this)(¤tValue, currentPoint.data); + if (currentValue == originValue) { + currentPoint += pitch * direction; + pitch /= 2.; + } else { + currentPoint -= pitch * direction; + pitch /= 2.; + } + } + } + + distance = (currentPoint - origin).norm(); + return true; +} + +template +bool IndicatorF3D::distance(S& distance, const Vector& origin) +{ + std::cout << "you should not be here!" << std::endl; + return true; +} + +template +bool IndicatorF3D::rotOnAxis(Vector& vec_rot, const Vector& vec, const Vector& axis, S& theta) +{ + /// http://mathworld.wolfram.com/RodriguesRotationFormula.html https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula + + //Vector axisN(axis*(1/(axis).norm())); // normalize rotation axis + Vector axisN(axis*(1/const_cast&> (axis).norm())); // normalize rotation axis + + vec_rot = vec ; + + Vector crossProd; + + crossProd[0] = axisN[1]*vec[2] - axisN[2]*vec[1]; + crossProd[1] = axisN[2]*vec[0] - axisN[0]*vec[2]; + crossProd[2] = axisN[0]*vec[1] - axisN[1]*vec[0]; + + S dotProd = axisN[0]*vec[0] + axisN[1]*vec[1] + axisN[2]*vec[2]; + + //v_rot = std::cos(theta)*vec + (crossProd)*std::sin(theta) + axisN*(dotProduct3D(axisN,vec))*(1 - std::cos(theta)); + vec_rot = std::cos(theta)*vec + (crossProd)*std::sin(theta) + axisN*(dotProd)*(1 - std::cos(theta)); + + return true; + +} + +template +bool IndicatorF3D::normal(Vector& normal, const Vector& origin, const Vector& direction, int iC) +{ + //OstreamManager clout(std::cout,"IndicatorF3D"); +#ifdef OLB_DEBUG + std::cout << "Calculating IndicatorF3D Normal" << std::endl; +#endif + + bool originValue; + (*this)(&originValue, origin.data); + Vector currentPoint(origin); + + S precision = .0001; + + S dist; + distance(dist, origin, direction, iC); + + S dirMag = const_cast&> (direction).norm(); +#ifdef OLB_DEBUG + std::cout << "magnitude = " << dirMag << std::endl; +#endif + + //Vector POS(origin + dist*direction*(1/const_cast&> (direction).norm())); //Point on Surface + //direction = direction*(1/const_cast&> (direction).norm()); + Vector directionN(direction*(1/dirMag)); + Vector POS(origin + dist*directionN); //Point on Surface + + /// find perpendicular vector to direction + Vector directionPerp; + if ( (util::nearZero(directionN[0]) && util::nearZero(directionN[1]) && !util::nearZero(directionN[2])) + || (util::nearZero(directionN[0]) && !util::nearZero(directionN[1]) && util::nearZero(directionN[2])) ) { + directionPerp[0] = 1; + directionPerp[1] = 0; + directionPerp[2] = 0; + } else if ( !util::nearZero(directionN[0]) && util::nearZero(directionN[1]) && util::nearZero(directionN[2]) ) { + directionPerp[0] = 0; + directionPerp[1] = 0; + directionPerp[2] = 1; + } else if ( ( !util::nearZero(directionN[0]) || !util::nearZero(directionN[1]) ) && !util::nearZero(directionN[2]) ) { + directionPerp[0] = directionN[0]; + directionPerp[1] = directionN[1]; + directionPerp[2] = -(directionN[0] + directionN[1])/directionN[2]; + } else { + std::cout << "Error: unknown case for perpendicular check" << std::endl; + return false; + } + + Vector directionPerpN(directionPerp*(dirMag/(directionPerp).norm())); + + + Vector point1; + Vector point2; + Vector point3; + + bool currentValue; + + /// Loop 3 times to find three points on the surface to use for normal calc. + /// orthogonal to direction vector 120 degree to each other + + + for (int n: { + 0,120,240 + }) { + S thetaMain = n*M_PI/180.; + + /// rotate directionPerpN through 3 angles {0,120,240} + Vector perp; + rotOnAxis(perp, directionPerpN, directionN, thetaMain); + Vector perpPoint(POS + perp); + + //std::cout << "perp = [" << perp[0] << "," << perp[1] << "," << perp[2] << "]" << std::endl; + + //S rotate = 90.; + S rotate = 179.; + S pitch = rotate/2.; + + Vector vec(perp); + + Vector rotAxis; + + //rotAxis = cross(perp,directionN); + rotAxis[0] = perp[1]*directionN[2] - perp[2]*directionN[1]; + rotAxis[1] = perp[2]*directionN[0] - perp[0]*directionN[2]; + rotAxis[2] = perp[0]*directionN[1] - perp[1]*directionN[0]; + + //normal = rotAxis; + //return true; + //std::cout << "rotAxis = [" << rotAxis[0] << "," << rotAxis[1] << "," << rotAxis[2] << "]" << std::endl; + + /// Find 'positive' angle + Vector testPOS; + S testAngle(45.*M_PI/180.); + rotOnAxis(testPOS, perp, rotAxis, testAngle); + Vector testPoint( POS + testPOS); + //std::cout << "testPOS = [" << testPOS[0] << "," << testPOS[1] << "," << testPOS[2] << "]" << std::endl; + //std::cout << "testPoint = [" << testPoint[0] << "," << testPoint[1] << "," << testPoint[2] << "]" << std::endl; + + S distTestPoint = (testPoint - origin).norm(); + S distPerpPoint = (perpPoint - origin).norm(); + + S mod = 0; + if (distTestPoint < distPerpPoint) { // pos. angle rotates towards + mod = -1; + } else { + mod = 1; + } + + while (std::abs(pitch) >= precision) { + + S theta(pitch*M_PI/180); + + currentPoint = POS + vec; + (*this)(¤tValue, currentPoint.data); + + S temp; + if (currentValue == originValue) { + temp = mod*theta; + rotOnAxis(vec, vec, rotAxis, temp); + } else { + temp = -mod*theta; + rotOnAxis(vec, vec, rotAxis, temp); + } + pitch /= 2.; + + + + } + + if (n == 0) { + point1 = currentPoint; + } else if (n == 120) { + point2 = currentPoint; + } else if (n == 240) { + point3 = currentPoint; + } else { + std::cout << "Something broke" << std::endl; + return false; + } + + } + + /// Calculate Normal + Vector vec1 (point1 - point2); + Vector vec2 (point1 - point3); + + normal[0] = -(vec1[1]*vec2[2] - vec1[2]*vec2[1]); + normal[1] = -(vec1[2]*vec2[0] - vec1[0]*vec2[2]); + normal[2] = -(vec1[0]*vec2[1] - vec1[1]*vec2[0]); + + normalize(normal); + + + //S dist; + //Vector dist; + //distance(dist, origin, direction, iC); + //normal = Vector(dist); + //normal = POS; + //normal = directionPerpN; + //normal = directionN; + return true; + +} + +template +bool IndicatorF3D::isInsideBox(Vector point) +{ + return point >= _myMin && point <= _myMax; +} + + +template +IndicatorIdentity3D::IndicatorIdentity3D(std::shared_ptr> f) + : _f(f) +{ + this->_myMin = _f->getMin(); + this->_myMax = _f->getMax(); +} + +template +bool IndicatorIdentity3D::operator() (bool output[1], const S input[3]) +{ + return (_f)->operator()(output, input); +} + + +} // namespace olb + +#endif -- cgit v1.2.3