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/MakeHeader | 37 ++ src/functors/analytical/indicator/indicCalc2D.cpp | 56 ++ src/functors/analytical/indicator/indicCalc2D.h | 127 ++++ src/functors/analytical/indicator/indicCalc2D.hh | 202 ++++++ src/functors/analytical/indicator/indicCalc3D.cpp | 50 ++ src/functors/analytical/indicator/indicCalc3D.h | 92 +++ src/functors/analytical/indicator/indicCalc3D.hh | 110 ++++ src/functors/analytical/indicator/indicator2D.h | 34 + src/functors/analytical/indicator/indicator2D.hh | 34 + src/functors/analytical/indicator/indicator3D.h | 34 + src/functors/analytical/indicator/indicator3D.hh | 35 + .../analytical/indicator/indicatorBaseF2D.cpp | 35 + .../analytical/indicator/indicatorBaseF2D.h | 106 +++ .../analytical/indicator/indicatorBaseF2D.hh | 210 ++++++ .../analytical/indicator/indicatorBaseF3D.cpp | 34 + .../analytical/indicator/indicatorBaseF3D.h | 77 +++ .../analytical/indicator/indicatorBaseF3D.hh | 321 +++++++++ src/functors/analytical/indicator/indicatorF2D.cpp | 41 ++ src/functors/analytical/indicator/indicatorF2D.h | 123 ++++ src/functors/analytical/indicator/indicatorF2D.hh | 256 ++++++++ src/functors/analytical/indicator/indicatorF3D.cpp | 57 ++ src/functors/analytical/indicator/indicatorF3D.h | 241 +++++++ src/functors/analytical/indicator/indicatorF3D.hh | 731 +++++++++++++++++++++ src/functors/analytical/indicator/module.mk | 27 + .../indicator/smoothIndicatorBaseF2D.cpp | 34 + .../analytical/indicator/smoothIndicatorBaseF2D.h | 144 ++++ .../analytical/indicator/smoothIndicatorBaseF2D.hh | 302 +++++++++ .../indicator/smoothIndicatorBaseF3D.cpp | 34 + .../analytical/indicator/smoothIndicatorBaseF3D.h | 146 ++++ .../analytical/indicator/smoothIndicatorBaseF3D.hh | 342 ++++++++++ .../indicator/smoothIndicatorCalcF2D.cpp | 35 + .../analytical/indicator/smoothIndicatorCalcF2D.h | 55 ++ .../analytical/indicator/smoothIndicatorCalcF2D.hh | 75 +++ .../indicator/smoothIndicatorCalcF3D.cpp | 34 + .../analytical/indicator/smoothIndicatorCalcF3D.h | 55 ++ .../analytical/indicator/smoothIndicatorCalcF3D.hh | 73 ++ .../analytical/indicator/smoothIndicatorF2D.cpp | 39 ++ .../analytical/indicator/smoothIndicatorF2D.h | 79 +++ .../analytical/indicator/smoothIndicatorF2D.hh | 239 +++++++ .../analytical/indicator/smoothIndicatorF3D.cpp | 40 ++ .../analytical/indicator/smoothIndicatorF3D.h | 95 +++ .../analytical/indicator/smoothIndicatorF3D.hh | 439 +++++++++++++ 42 files changed, 5330 insertions(+) create mode 100644 src/functors/analytical/indicator/MakeHeader create mode 100644 src/functors/analytical/indicator/indicCalc2D.cpp create mode 100644 src/functors/analytical/indicator/indicCalc2D.h create mode 100644 src/functors/analytical/indicator/indicCalc2D.hh create mode 100644 src/functors/analytical/indicator/indicCalc3D.cpp create mode 100644 src/functors/analytical/indicator/indicCalc3D.h create mode 100644 src/functors/analytical/indicator/indicCalc3D.hh create mode 100644 src/functors/analytical/indicator/indicator2D.h create mode 100644 src/functors/analytical/indicator/indicator2D.hh create mode 100644 src/functors/analytical/indicator/indicator3D.h create mode 100644 src/functors/analytical/indicator/indicator3D.hh create mode 100644 src/functors/analytical/indicator/indicatorBaseF2D.cpp create mode 100644 src/functors/analytical/indicator/indicatorBaseF2D.h create mode 100644 src/functors/analytical/indicator/indicatorBaseF2D.hh create mode 100644 src/functors/analytical/indicator/indicatorBaseF3D.cpp create mode 100644 src/functors/analytical/indicator/indicatorBaseF3D.h create mode 100644 src/functors/analytical/indicator/indicatorBaseF3D.hh create mode 100644 src/functors/analytical/indicator/indicatorF2D.cpp create mode 100644 src/functors/analytical/indicator/indicatorF2D.h create mode 100644 src/functors/analytical/indicator/indicatorF2D.hh create mode 100644 src/functors/analytical/indicator/indicatorF3D.cpp create mode 100644 src/functors/analytical/indicator/indicatorF3D.h create mode 100644 src/functors/analytical/indicator/indicatorF3D.hh create mode 100644 src/functors/analytical/indicator/module.mk create mode 100644 src/functors/analytical/indicator/smoothIndicatorBaseF2D.cpp create mode 100644 src/functors/analytical/indicator/smoothIndicatorBaseF2D.h create mode 100644 src/functors/analytical/indicator/smoothIndicatorBaseF2D.hh create mode 100644 src/functors/analytical/indicator/smoothIndicatorBaseF3D.cpp create mode 100644 src/functors/analytical/indicator/smoothIndicatorBaseF3D.h create mode 100644 src/functors/analytical/indicator/smoothIndicatorBaseF3D.hh create mode 100644 src/functors/analytical/indicator/smoothIndicatorCalcF2D.cpp create mode 100644 src/functors/analytical/indicator/smoothIndicatorCalcF2D.h create mode 100644 src/functors/analytical/indicator/smoothIndicatorCalcF2D.hh create mode 100644 src/functors/analytical/indicator/smoothIndicatorCalcF3D.cpp create mode 100644 src/functors/analytical/indicator/smoothIndicatorCalcF3D.h create mode 100644 src/functors/analytical/indicator/smoothIndicatorCalcF3D.hh create mode 100644 src/functors/analytical/indicator/smoothIndicatorF2D.cpp create mode 100644 src/functors/analytical/indicator/smoothIndicatorF2D.h create mode 100644 src/functors/analytical/indicator/smoothIndicatorF2D.hh create mode 100644 src/functors/analytical/indicator/smoothIndicatorF3D.cpp create mode 100644 src/functors/analytical/indicator/smoothIndicatorF3D.h create mode 100644 src/functors/analytical/indicator/smoothIndicatorF3D.hh (limited to 'src/functors/analytical/indicator') diff --git a/src/functors/analytical/indicator/MakeHeader b/src/functors/analytical/indicator/MakeHeader new file mode 100644 index 0000000..66496d7 --- /dev/null +++ b/src/functors/analytical/indicator/MakeHeader @@ -0,0 +1,37 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias 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. + + +generic := + +precompiled := indicatorBaseF2D \ + indicatorBaseF3D \ + indicatorF2D \ + indicatorF3D \ + indicCalc2D \ + indicCalc3D \ + smoothIndicatorF2D \ + smoothIndicatorF3D \ + smoothIndicatorBaseF2D \ + smoothIndicatorBaseF3D \ + smoothIndicatorCalcF2D \ + smoothIndicatorCalcF3D diff --git a/src/functors/analytical/indicator/indicCalc2D.cpp b/src/functors/analytical/indicator/indicCalc2D.cpp new file mode 100644 index 0000000..ba54bde --- /dev/null +++ b/src/functors/analytical/indicator/indicCalc2D.cpp @@ -0,0 +1,56 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2014-2016 Mathias J. Krause, Cyril Masquelier, + * Benjamin Förster, Albert Mink + * 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. +*/ + +#include "indicCalc2D.h" +#include "indicCalc2D.hh" + + +namespace olb { + +// arithmetic helper class for Indical 1d functors +template class IndicCalc1D; +template class IndicPlus1D; +template class IndicMinus1D; +template class IndicMultiplication1D; + +// arithmetic helper class for Indical 2d functors +template class IndicCalc2D; +template class IndicCalc2D; +template class IndicCalc2D; + +template std::shared_ptr> operator+(std::shared_ptr> lhs, + std::shared_ptr> rhs); +template std::shared_ptr> operator-(std::shared_ptr> lhs, + std::shared_ptr> rhs); +template std::shared_ptr> operator*(std::shared_ptr> lhs, + std::shared_ptr> rhs); + +template std::shared_ptr> operator+(IndicatorIdentity2D & lhs, + std::shared_ptr>); +template std::shared_ptr> operator-(IndicatorIdentity2D & lhs, + std::shared_ptr>); +template std::shared_ptr> operator*(IndicatorIdentity2D & lhs, + std::shared_ptr>); + +} diff --git a/src/functors/analytical/indicator/indicCalc2D.h b/src/functors/analytical/indicator/indicCalc2D.h new file mode 100644 index 0000000..a5f6449 --- /dev/null +++ b/src/functors/analytical/indicator/indicCalc2D.h @@ -0,0 +1,127 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2014-2016 Mathias J. Krause, Cyril Masquelier, + * Benjamin Förster, Albert Mink + * 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 INDIC_CALC_2D_H +#define INDIC_CALC_2D_H + +#include "indicatorBaseF2D.h" +#include "utilities/arithmetic.h" + +namespace olb { + + +/* + * arithmetic helper classes for IndicatorF1D, IndicatorF2D, smoothIndicator2D + * UNION + + * WITHOUT - + * INTERSECTION * +*/ + +//////////////////////////////// IndicCalc1D //////////////////////////////// +/// arithmetic helper class for Indicator 1d functors +template +class IndicCalc1D : public IndicatorF1D { +protected: + IndicatorF1D& _f; + IndicatorF1D& _g; +public: + // set image/target dimensions of IndicCalc1D as well + IndicCalc1D(IndicatorF1D& f, IndicatorF1D& g); +}; + +/// addition functor acts as union +template +class IndicPlus1D : public IndicCalc1D { +public: + IndicPlus1D(IndicatorF1D& f, IndicatorF1D& g); + bool operator() (bool output[], const S input[]) override; +}; + +/// subtraction functor acts as without +template +class IndicMinus1D : public IndicCalc1D { +public: + IndicMinus1D(IndicatorF1D& f, IndicatorF1D& g); + bool operator() (bool output[], const S input[]) override; +}; + +/// multiplication functor acts as intersection +template +class IndicMultiplication1D : public IndicCalc1D { +public: + IndicMultiplication1D(IndicatorF1D& f, IndicatorF1D& g); + bool operator() (bool output[], const S input[]) override; +}; + + + +//////////////////////////////// indicCalc2D //////////////////////////////// +/// arithmetic helper class for Indicator 2D functors +template class F> +class IndicCalc2D : public IndicatorF2D { +protected: + std::shared_ptr> _f; + std::shared_ptr> _g; +public: + IndicCalc2D( std::shared_ptr> f, std::shared_ptr> g ); + + bool operator() (bool output[], const S input[2]) override; +}; + +/// Addition functor (W==bool: Union) +template +using IndicPlus2D = IndicCalc2D; +template +using IndicMinus2D = IndicCalc2D; +template +using IndicMultiplication2D = IndicCalc2D; + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator+(std::shared_ptr> lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator-(std::shared_ptr> lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator*(std::shared_ptr> lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator+(F1 & lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator-(F1 & lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator*(F1 & lhs, std::shared_ptr> rhs); + + +} // end namespace olb + +#endif diff --git a/src/functors/analytical/indicator/indicCalc2D.hh b/src/functors/analytical/indicator/indicCalc2D.hh new file mode 100644 index 0000000..f325bc7 --- /dev/null +++ b/src/functors/analytical/indicator/indicCalc2D.hh @@ -0,0 +1,202 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2014-2016 Mathias J. Krause, Cyril Masquelier, + * Benjamin Förster, Albert Mink + * 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 INDIC_CALC_2D_HH +#define INDIC_CALC_2D_HH + +#include "indicCalc2D.h" + +namespace olb { + + +//////////////////////////////// IndicCalc1D //////////////////////////////// +template +IndicCalc1D::IndicCalc1D(IndicatorF1D& f, IndicatorF1D& g) + : _f(f), _g(g) +{ + this->_myMin[0] = std::min(f.getMin()[0], g.getMin()[0]); + this->_myMax[0] = std::max(f.getMax()[0], g.getMax()[0]); + std::swap(f._ptrCalcC, this->_ptrCalcC); +} + + +template +IndicPlus1D::IndicPlus1D(IndicatorF1D& f, IndicatorF1D& g) + : IndicCalc1D(f, g) +{} + +// returns 1 if( f==1 || g==1 ) UNION +template +bool IndicPlus1D::operator() (bool output[], const S input[]) +{ + this->_f(output, input); + bool tmp; + this->_g(&tmp, input); + output[0] |= tmp; + return true; +} + + +template +IndicMinus1D::IndicMinus1D(IndicatorF1D& f, IndicatorF1D& g) + : IndicCalc1D(f, g) +{} + +// returns 1 if( f==1 && g==0 ) WITHOUT +template +bool IndicMinus1D::operator()(bool output[], const S input[]) +{ + this->_f(output, input); + bool tmp; + this->_g(&tmp, input); + output[0] &= !tmp; + return true; +} + + + +template +IndicMultiplication1D::IndicMultiplication1D(IndicatorF1D& f, IndicatorF1D& g) + : IndicCalc1D(f, g) +{} + +// returns 1 if( f==1 && g==1 ) INTERSECTION +template +bool IndicMultiplication1D::operator() (bool output[], const S input[]) +{ + this->_f(output, input); + bool tmp; + this->_g(&tmp, input); + output[0] &= tmp; + return true; +} + + + +template +IndicatorF1D& IndicatorF1D::operator+(IndicatorF1D& rhs) +{ + auto tmp = std::make_shared< IndicPlus1D >(*this,rhs); + this->_ptrCalcC = tmp; + return *tmp; +} + +template +IndicatorF1D& IndicatorF1D::operator-(IndicatorF1D& rhs) +{ + auto tmp = std::make_shared< IndicMinus1D >(*this,rhs); + this->_ptrCalcC = tmp; + return *tmp; +} + +template +IndicatorF1D& IndicatorF1D::operator*(IndicatorF1D& rhs) +{ + auto tmp = std::make_shared< IndicMultiplication1D >(*this,rhs); + this->_ptrCalcC = tmp; + return *tmp; +} + + + + +//////////////////////////////// IndicCalc2D //////////////////////////////// +template class F> +IndicCalc2D::IndicCalc2D(std::shared_ptr> f, std::shared_ptr> g) + : _f(f), _g(g) +{ + for ( int i=0; i<2; i++) { + this->_myMin[i] = std::min(_f->getMin()[i], _g->getMin()[i]); + this->_myMax[i] = std::max(_f->getMax()[i], _g->getMax()[i]); + } +} + +template class F> +bool IndicCalc2D::operator()( bool output[], const S input[2]) +{ + // componentwise operation on equidimensional functors + bool* outputF = output; + _f->operator()(outputF, input); + + bool outputG[this->getTargetDim()]; + _g->operator()(outputG, input); + + for (int i = 0; i < this->getTargetDim(); i++) { + output[i] = F()(outputF[i], outputG[i]); + } + return output; +} + + + + + +//// no association to a operator+ from a class is needed, thus we have these free functions +template class F1, template class F2, + typename> +std::shared_ptr> operator+(std::shared_ptr> lhs, std::shared_ptr> rhs) +{ + return std::make_shared>(lhs, rhs); +} + +template class F1, template class F2, + typename> +std::shared_ptr> operator-(std::shared_ptr> lhs, std::shared_ptr> rhs) +{ + return std::make_shared>(lhs, rhs); +} + +template class F1, template class F2, + typename> +std::shared_ptr> operator*(std::shared_ptr> lhs, std::shared_ptr> rhs) +{ + return std::make_shared>(lhs, rhs); +} + +// template specialization for indicatorIdentity +template class F1, template class F2, + typename> +std::shared_ptr> operator+(F1 & lhs, std::shared_ptr> rhs) +{ + return lhs._f + rhs; +} + +template class F1, template class F2, + typename> +std::shared_ptr> operator-(F1 & lhs, std::shared_ptr> rhs) +{ + return lhs._f - rhs; +} + +template class F1, template class F2, + typename> +std::shared_ptr> operator*(F1 & lhs, std::shared_ptr> rhs) +{ + return lhs._f * rhs; +} + + +} // end namespace olb + +#endif diff --git a/src/functors/analytical/indicator/indicCalc3D.cpp b/src/functors/analytical/indicator/indicCalc3D.cpp new file mode 100644 index 0000000..8ae86cf --- /dev/null +++ b/src/functors/analytical/indicator/indicCalc3D.cpp @@ -0,0 +1,50 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Albert Mink + * 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. +*/ + +#include "indicCalc3D.h" +#include "indicCalc3D.hh" + + +namespace olb { + +// arithmetic helper class for indicator 3d functors +template class IndicCalc3D; +template class IndicCalc3D; +template class IndicCalc3D; + +template std::shared_ptr> operator+(std::shared_ptr> lhs, + std::shared_ptr> rhs); +template std::shared_ptr> operator-(std::shared_ptr> lhs, + std::shared_ptr> rhs); +template std::shared_ptr> operator*(std::shared_ptr> lhs, + std::shared_ptr> rhs); + +template std::shared_ptr> operator+(IndicatorIdentity3D & lhs, + std::shared_ptr>); +template std::shared_ptr> operator-(IndicatorIdentity3D & lhs, + std::shared_ptr>); +template std::shared_ptr> operator*(IndicatorIdentity3D & lhs, + std::shared_ptr>); + +} + diff --git a/src/functors/analytical/indicator/indicCalc3D.h b/src/functors/analytical/indicator/indicCalc3D.h new file mode 100644 index 0000000..8d7f667 --- /dev/null +++ b/src/functors/analytical/indicator/indicCalc3D.h @@ -0,0 +1,92 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Albert Mink + * 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 INDIC_CALC_3D_H +#define INDIC_CALC_3D_H + +#include "utilities/arithmetic.h" +#include "indicatorBaseF3D.h" + +namespace olb { + + +/* + * arithmetic helper classes for IndicatorF3D, smoothIndicator3D + * UNION + + * WITHOUT - + * INTERSECTION * +*/ + +//////////////////////////////// indicCalc3D //////////////////////////////// +/// arithmetic helper class for Indicator 3d functors +template class F> +class IndicCalc3D : public IndicatorF3D { +protected: + std::shared_ptr> _f; + std::shared_ptr> _g; +public: + IndicCalc3D( std::shared_ptr> f, std::shared_ptr> g ); + + bool operator() (bool output[], const S input[3]) override; +}; + +/// Addition functor (W==bool: Union) +template +using IndicPlus3D = IndicCalc3D; +template +using IndicMinus3D = IndicCalc3D; +template +using IndicMultiplication3D = IndicCalc3D; + +/** Free function implements lhs+rhs, only for IndicaotrsF3D types through enable_if and is_base_of + * + * \tparam S usual type for source dimension of the functor + * \tparam F1 lhs has to be derived from IndicatorF3D, otherwise function is disabled + * \tparam F2 rhs + */ +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator+(std::shared_ptr> lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator-(std::shared_ptr> lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator*(std::shared_ptr> lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator+(F1 & lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator-(F1 & lhs, std::shared_ptr> rhs); + +template class F1, template class F2, + typename=typename std::enable_if, F1>::value>::type> +std::shared_ptr> operator*(F1 & lhs, std::shared_ptr> rhs); +} // end namespace olb + +#endif diff --git a/src/functors/analytical/indicator/indicCalc3D.hh b/src/functors/analytical/indicator/indicCalc3D.hh new file mode 100644 index 0000000..d087ed9 --- /dev/null +++ b/src/functors/analytical/indicator/indicCalc3D.hh @@ -0,0 +1,110 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Albert Mink + * 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 INDIC_CALC_3D_HH +#define INDIC_CALC_3D_HH + +#include "indicCalc3D.h" + +namespace olb { + + +template class F> +IndicCalc3D::IndicCalc3D(std::shared_ptr> f, std::shared_ptr> g) + : _f(f), _g(g) +{ + for ( int i=0; i<3; i++) { + this->_myMin[i] = std::min(_f->getMin()[i], _g->getMin()[i]); + this->_myMax[i] = std::max(_f->getMax()[i], _g->getMax()[i]); + } +} + + + +template class F> +bool IndicCalc3D::operator()( bool output[], const S input[3]) +{ + // componentwise operation on equidimensional functors + bool* outputF = output; + _f->operator()(outputF, input); + + bool outputG[this->getTargetDim()]; + _g->operator()(outputG, input); + + for (int i = 0; i < this->getTargetDim(); i++) { + output[i] = F()(outputF[i], outputG[i]); + } + return output; +} + + + + +//// no association to a operator+ from a class is needed, thus we have these free functions +template class F1, template class F2, + typename> +std::shared_ptr> operator+(std::shared_ptr> lhs, std::shared_ptr> rhs) +{ + return std::make_shared>(lhs, rhs); +} + +template class F1, template class F2, + typename> +std::shared_ptr> operator-(std::shared_ptr> lhs, std::shared_ptr> rhs) +{ + return std::make_shared>(lhs, rhs); +} + +template class F1, template class F2, + typename> +std::shared_ptr> operator*(std::shared_ptr> lhs, std::shared_ptr> rhs) +{ + return std::make_shared>(lhs, rhs); +} + +// template specialization for indicatorIdentity +template class F1, template class F2, + typename> +std::shared_ptr> operator+(F1 & lhs, std::shared_ptr> rhs) +{ + return lhs._f + rhs; +} + +template class F1, template class F2, + typename> +std::shared_ptr> operator-(F1 & lhs, std::shared_ptr> rhs) +{ + return lhs._f - rhs; +} + +template class F1, template class F2, + typename> +std::shared_ptr> operator*(F1 & lhs, std::shared_ptr> rhs) +{ + return lhs._f * rhs; +} + + +} // end namespace olb + +#endif diff --git a/src/functors/analytical/indicator/indicator2D.h b/src/functors/analytical/indicator/indicator2D.h new file mode 100644 index 0000000..9da168e --- /dev/null +++ b/src/functors/analytical/indicator/indicator2D.h @@ -0,0 +1,34 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012 Lukas Baron, 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 + * Groups all include files for the directory genericFunctions. + */ + +#include "indicatorBaseF2D.h" +#include "indicatorF2D.h" +#include "indicCalc2D.h" + +#include "smoothIndicatorBaseF2D.h" +#include "smoothIndicatorF2D.h" +#include "smoothIndicatorCalcF2D.h" diff --git a/src/functors/analytical/indicator/indicator2D.hh b/src/functors/analytical/indicator/indicator2D.hh new file mode 100644 index 0000000..1ea92d9 --- /dev/null +++ b/src/functors/analytical/indicator/indicator2D.hh @@ -0,0 +1,34 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012 Lukas Baron, 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 + * Groups all include files for the directory genericFunctions. + */ + +#include "indicatorBaseF2D.hh" +#include "indicatorF2D.hh" +#include "indicCalc2D.hh" + +#include "smoothIndicatorBaseF2D.hh" +#include "smoothIndicatorF2D.hh" +#include "smoothIndicatorCalcF2D.hh" diff --git a/src/functors/analytical/indicator/indicator3D.h b/src/functors/analytical/indicator/indicator3D.h new file mode 100644 index 0000000..73d390e --- /dev/null +++ b/src/functors/analytical/indicator/indicator3D.h @@ -0,0 +1,34 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012-2016 Lukas Baron, Mathias J. Krause, Albert Mink, 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. +*/ + +/** \file + * Groups all include files for the directory genericFunctions. + */ + +#include "indicatorBaseF3D.h" +#include "indicatorF3D.h" +#include "indicCalc3D.h" + +#include "smoothIndicatorBaseF3D.h" +#include "smoothIndicatorF3D.h" +#include "smoothIndicatorCalcF3D.h" diff --git a/src/functors/analytical/indicator/indicator3D.hh b/src/functors/analytical/indicator/indicator3D.hh new file mode 100644 index 0000000..8be9c86 --- /dev/null +++ b/src/functors/analytical/indicator/indicator3D.hh @@ -0,0 +1,35 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012-2016 Lukas Baron, Mathias J. Krause, + * Albert Mink, 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. +*/ + +/** \file + * Groups all include files for the directory genericFunctions. + */ + +#include "indicatorBaseF3D.hh" +#include "indicatorF3D.hh" +#include "indicCalc3D.hh" + +#include "smoothIndicatorBaseF3D.hh" +#include "smoothIndicatorF3D.hh" +#include "smoothIndicatorCalcF3D.hh" diff --git a/src/functors/analytical/indicator/indicatorBaseF2D.cpp b/src/functors/analytical/indicator/indicatorBaseF2D.cpp new file mode 100644 index 0000000..578aa82 --- /dev/null +++ b/src/functors/analytical/indicator/indicatorBaseF2D.cpp @@ -0,0 +1,35 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2015-2016 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. +*/ + + +#include "indicatorBaseF2D.h" +#include "indicatorBaseF2D.hh" + +namespace olb { + +template class IndicatorF1D; + +template class IndicatorF2D; +template class IndicatorIdentity2D; + +} diff --git a/src/functors/analytical/indicator/indicatorBaseF2D.h b/src/functors/analytical/indicator/indicatorBaseF2D.h new file mode 100644 index 0000000..4f76b28 --- /dev/null +++ b/src/functors/analytical/indicator/indicatorBaseF2D.h @@ -0,0 +1,106 @@ +/* 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_H +#define INDICATOR_BASE_F_2D_H + +#include + +#include "core/vector.h" +#include "functors/genericF.h" + +namespace olb { + + +/** IndicatorF1D is an application from \f$ \Omega \subset R^3 \to {0,1} \f$. + * \param _myMin holds minimal(component wise) vector of the domain \f$ \Omega \f$. + * \param _myMax holds maximal(component wise) vector of the domain \f$ \Omega \f$. + */ +template +class IndicatorF1D : public GenericF { +protected: + IndicatorF1D(); + Vector _myMin; + Vector _myMax; +public: + virtual Vector& getMin(); + virtual Vector& getMax(); + IndicatorF1D& operator+(IndicatorF1D& rhs); + IndicatorF1D& operator-(IndicatorF1D& rhs); + IndicatorF1D& operator*(IndicatorF1D& rhs); +}; + + +/** IndicatorF2D is an application from \f$ \Omega \subset R^3 \to {0,1} \f$. + * \param _myMin holds minimal(component wise) vector of the domain \f$ \Omega \f$. + * \param _myMax holds maximal(component wise) vector of the domain \f$ \Omega \f$. + */ +template +class IndicatorF2D : public GenericF { +protected: + IndicatorF2D(); + Vector _myMin; + Vector _myMax; +public: + using GenericF::operator(); + + virtual Vector& getMin(); + virtual Vector& getMax(); + + /// returns false or true and pos. distance if there was one found for an given origin and direction + /** + * (mind that the default computation is done by a numerical approximation which searches .. [TODO]) + */ + virtual bool distance(S& distance, const Vector& origin, const Vector& direction, int iC=-1); + + /// returns true and the normal if there was one found for an given origin and direction + /** + * (mind that the default computation is done by a numerical approximation which searches .. [TODO]) + */ + virtual bool normal(Vector& normal, const Vector& origin, const Vector& direction, int iC=-1); + + /// Returns true if `point` is inside a cube with corners `_myMin` and `_myMax` + bool isInsideBox(Vector point); + + /// Indicator specific function operator overload + /** + * \return Domain indicator i.e. `true` iff the input lies within the described domain. + **/ + virtual bool operator() (const S input[]); +}; + + +/////////////////////////////////////IdentityF////////////////////////////////// +template +class IndicatorIdentity2D : public IndicatorF2D { +public: + std::shared_ptr> _f; + + IndicatorIdentity2D(std::shared_ptr> f); + bool operator() (bool output[1], const S input[2]) override; +}; + + +} + +#endif 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 diff --git a/src/functors/analytical/indicator/indicatorBaseF3D.cpp b/src/functors/analytical/indicator/indicatorBaseF3D.cpp new file mode 100644 index 0000000..09936f0 --- /dev/null +++ b/src/functors/analytical/indicator/indicatorBaseF3D.cpp @@ -0,0 +1,34 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2015-2016 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. +*/ + + +#include "indicatorBaseF3D.h" +#include "indicatorBaseF3D.hh" +#include "indicCalc3D.hh" + +namespace olb { + +template class IndicatorF3D; +template class IndicatorIdentity3D; + +} diff --git a/src/functors/analytical/indicator/indicatorBaseF3D.h b/src/functors/analytical/indicator/indicatorBaseF3D.h new file mode 100644 index 0000000..d7eced8 --- /dev/null +++ b/src/functors/analytical/indicator/indicatorBaseF3D.h @@ -0,0 +1,77 @@ +/* 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_H +#define INDICATOR_BASE_F_3D_H + +#include + +#include "core/vector.h" +#include "functors/genericF.h" + +namespace olb { + + +/** IndicatorF3D is an application from \f$ \Omega \subset R^3 \to \{0,1\} \f$. + * \param _myMin holds minimal(component wise) vector of the domain \f$ \Omega \f$. + * \param _myMax holds maximal(component wise) vector of the domain \f$ \Omega \f$. + */ +template +class IndicatorF3D : public GenericF { +protected: + IndicatorF3D(); + Vector _myMin; + Vector _myMax; +public: + virtual Vector& getMin(); + virtual Vector& getMax(); + /** \returns false or true and pos. distance if there was one found for a given origin and direction. + * Mind that the default computation is done by a numerical approximation which searches .. [TODO: CYRIL] + */ + virtual bool distance(S& distance, const Vector& origin, const Vector& direction, int iC=-1); + virtual bool distance(S& distance, const Vector& origin); + /// returns true and the normal if there was one found for an given origin and direction + /** + * (mind that the default computation is done by a numerical approximation which searches .. [TODO]) + */ + virtual bool normal(Vector& normal, const Vector& origin, const Vector& direction, int iC=-1); + ///Rotate vector around axis by angle theta + virtual bool rotOnAxis(Vector& vec_rot, const Vector& vec, const Vector& axis, S& theta); + /// Returns true if `point` is inside a cube with corners `_myMin` and `_myMax` + bool isInsideBox(Vector point); +}; + + +template +class IndicatorIdentity3D : public IndicatorF3D { +public: + std::shared_ptr> _f; + + IndicatorIdentity3D(std::shared_ptr> f); + bool operator() (bool output[1], const S input[3]) override; +}; + + +} + +#endif 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 diff --git a/src/functors/analytical/indicator/indicatorF2D.cpp b/src/functors/analytical/indicator/indicatorF2D.cpp new file mode 100644 index 0000000..dc57e54 --- /dev/null +++ b/src/functors/analytical/indicator/indicatorF2D.cpp @@ -0,0 +1,41 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2013 Lukas Baron, 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. +*/ + +#include "io/stlReader.h" +#include "indicatorF2D.h" +#include "indicatorF2D.hh" + +namespace olb { + +// st