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