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/indicatorBaseF2D.hh | 210 +++++++++++++++++++++
1 file changed, 210 insertions(+)
create mode 100644 src/functors/analytical/indicator/indicatorBaseF2D.hh
(limited to 'src/functors/analytical/indicator/indicatorBaseF2D.hh')
diff --git a/src/functors/analytical/indicator/indicatorBaseF2D.hh b/src/functors/analytical/indicator/indicatorBaseF2D.hh
new file mode 100644
index 0000000..0127b96
--- /dev/null
+++ b/src/functors/analytical/indicator/indicatorBaseF2D.hh
@@ -0,0 +1,210 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2014-2016 Cyril Masquelier, 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_2D_HH
+#define INDICATOR_BASE_F_2D_HH
+
+#include
+#include "indicatorBaseF2D.h"
+#include "math.h"
+
+namespace olb {
+
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+template
+IndicatorF1D::IndicatorF1D()
+ : GenericF(1, 1)
+{ }
+
+template
+Vector& IndicatorF1D::getMin()
+{
+ return _myMin;
+}
+
+template
+Vector& IndicatorF1D::getMax()
+{
+ return _myMax;
+}
+
+
+template
+IndicatorF2D::IndicatorF2D()
+ : GenericF(1, 2)
+{ }
+
+template
+Vector& IndicatorF2D::getMin()
+{
+ return _myMin;
+}
+
+template
+Vector& IndicatorF2D::getMax()
+{
+ return _myMax;
+}
+
+template
+bool IndicatorF2D::distance(S& distance, const Vector& origin, const Vector& direction, int iC)
+{
+ bool originValue;
+ (*this)(&originValue, origin.data);
+ Vector currentPoint(origin);
+
+ S precision = .0001;
+ S pitch = 0.5;
+
+ bool currentValue;
+ (*this)(¤tValue, currentPoint.data);
+ while (currentValue == originValue && isInsideBox(currentPoint)) {
+ currentPoint += direction;
+ (*this)(¤tValue, currentPoint.data);//changed until first point on the other side (inside/outside) is found
+ }
+
+ 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;
+}
+
+// given origin (inside) and direction first calculate distance to surface
+// go -90� to direction using POS as origin and check if inside/outside
+// if inside rotate outward, if outside rotate inward
+// iterate until close enough to surface, then store point1
+// repeat for 90� and store point2
+// use point1 and point2 on surface to calculate normal
+template
+bool IndicatorF2D::normal(Vector& normal, const Vector& origin, const Vector& direction, int iC)
+{
+ //OstreamManager clout(std::cout,"normal");
+ //clout << "Calculating IndicatorF2D Normal " << endl;
+ bool originValue;
+ (*this)(&originValue, origin.data);
+ Vector currentPoint(origin);
+
+ S precision = .0001;
+
+ S dist;
+ distance(dist, origin, direction, iC);
+
+
+ Vector POS(origin + dist*direction*(1/const_cast&> (direction).norm())); //Point on Surface
+
+ Vector point1;
+ Vector point2;
+
+ bool currentValue;
+
+ for (int n: {
+ -90,90
+ }) { //std::range n = {-90, 90};
+ S rotate(n);
+ S pitch(rotate/2.);
+ while (std::abs(pitch) >= precision) {
+ S theta(rotate*M_PI/180.);
+
+ Vector vec(std::cos(theta)*direction[0]+std::sin(theta)*direction[1],-std::sin(theta)*direction[0]+std::cos(theta)*direction[1]);
+ currentPoint = POS + vec;
+ (*this)(¤tValue, currentPoint.data);
+
+ if (currentValue == originValue) {
+ rotate -= pitch;
+ } else {
+ rotate += pitch;
+ }
+ pitch /= 2.;
+ }
+
+ if (n == -90) {
+ point1 = currentPoint;
+ } else if (n == 90) {
+ point2 = currentPoint;
+ }
+ }
+ // Calculate Normal from point1 and point2
+ normal = Vector((point2[1] - point1[1]), (-1)*(point2[0] - point1[0]));
+
+
+ //S dist;
+ //Vector dist;
+ //distance(dist, origin, direction, iC);
+ //normal = Vector(dist);
+ return true;
+}
+
+template
+bool IndicatorF2D::isInsideBox(Vector point)
+{
+ return point >= _myMin && point <= _myMax;
+}
+
+template
+bool IndicatorF2D::operator() (const S input[])
+{
+ bool output;
+ this->operator()(&output, input);
+ return output;
+}
+
+template
+IndicatorIdentity2D::IndicatorIdentity2D(std::shared_ptr> f)
+ : _f(f)
+{
+ this->_myMin = _f->getMin();
+ this->_myMax = _f->getMax();
+}
+
+template
+bool IndicatorIdentity2D::operator() (bool output[1], const S input[2])
+{
+ return (_f)->operator()(output, input);
+}
+
+} // namespace olb
+
+#endif
--
cgit v1.2.3