diff options
Diffstat (limited to 'src/functors/analytical/frameChangeF2D.hh')
| -rw-r--r-- | src/functors/analytical/frameChangeF2D.hh | 297 |
1 files changed, 297 insertions, 0 deletions
diff --git a/src/functors/analytical/frameChangeF2D.hh b/src/functors/analytical/frameChangeF2D.hh new file mode 100644 index 0000000..7b9811d --- /dev/null +++ b/src/functors/analytical/frameChangeF2D.hh @@ -0,0 +1,297 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2014-2015 Mathias J. Krause, Marie-Luise Maier + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * <http://www.openlb.net/> + * + * 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 ANALYTICAL_FRAME_CHANGE_F_2D_HH +#define ANALYTICAL_FRAME_CHANGE_F_2D_HH + +#include<vector> +#include<cmath> +#include<string> + +#include "frameChangeF2D.h" +#include "frameChangeF3D.h" +#include "functors/genericF.h" +#include "analyticalF.h" +#include "functors/lattice/superBaseF2D.h" +#include "geometry/superGeometry2D.h" +#include "functors/lattice/superLatticeLocalF2D.h" +#include "core/superLattice2D.h" +#include "dynamics/lbHelpers.h" // for computation of lattice rho and velocity +#include "utilities/vectorHelpers.h" // for normalize + +using namespace olb::util; + +namespace olb { + +template <typename T> +PowerLaw2D<T>::PowerLaw2D(std::vector<T> axisPoint, std::vector<T> axisDirection, T maxVelocity, T radius, T exponent) : AnalyticalF2D<T,T>(2) +{ + this->getName() = "PowerLaw2D"; + _axisPoint.resize(2); + _axisDirection.resize(2); + for (int i = 0; i < 2; ++i) { + _axisDirection[i] = axisDirection[i]; + _axisPoint[i] = axisPoint[i]; + } + _maxVelocity = maxVelocity; + _radius = radius; + _exponent = exponent; +} + +template <typename T> +PowerLaw2D<T>::PowerLaw2D(SuperGeometry2D<T>& superGeometry, int material, T maxVelocity, T distance2Wall, T exponent) : AnalyticalF2D<T,T>(2) +{ + this->getName() = "PowerLaw2D"; + _axisPoint.resize(2); + _axisDirection.resize(2); + _axisPoint = superGeometry.getStatistics().getCenterPhysR(material); + std::vector<int> discreteNormal = superGeometry.getStatistics().computeDiscreteNormal(material); + for (int i = 0; i < 2; ++i) { + _axisDirection[i] = discreteNormal[i]; + } + + _radius = T(distance2Wall); + for (int iD = 0; iD < 2; iD++) { + _radius += (superGeometry.getStatistics().getPhysRadius(material)[iD]); + } + _maxVelocity = maxVelocity; + _exponent = exponent; +} + + +template <typename T> +bool PowerLaw2D<T>::operator()(T output[], const T x[]) +{ + T d = fabs(_axisDirection[1]*(x[0] - _axisPoint[0]) - _axisDirection[0]*(x[1] - _axisPoint[1])); + output[0] = _maxVelocity*_axisDirection[0]*(1. - pow(d/_radius,_exponent)); + output[1] = _maxVelocity*_axisDirection[1]*(1. - pow(d/_radius,_exponent)); + if ( 1. - pow(d/_radius,_exponent) < 0.) { + output[0] = T(); + output[1] = T(); + } + return true; +} + + +template <typename T> +Poiseuille2D<T>::Poiseuille2D(std::vector<T> axisPoint, std::vector<T> axisDirection, T maxVelocity, T radius) : PowerLaw2D<T>(axisPoint, axisDirection, maxVelocity, radius, 2) +{ + this->getName() = "Poiseuille2D"; +} + + +template <typename T> +Poiseuille2D<T>::Poiseuille2D(SuperGeometry2D<T>& superGeometry, int material, T maxVelocity, T distance2Wall) : PowerLaw2D<T>(superGeometry, material, maxVelocity, distance2Wall, 2) +{ + this->getName() = "Poiseuille2D"; +} + + +template <typename T, typename S, typename DESCRIPTOR> +PoiseuilleStrainRate2D<T,S,DESCRIPTOR>::PoiseuilleStrainRate2D(UnitConverter<T,DESCRIPTOR> const& converter, T ly) + : AnalyticalF2D<T,S>(4), _lengthY(ly), _maxVelocity(converter.getCharPhysVelocity()) +{ + this->getName() = "PoiseuilleStrainRate2D"; +} + + +template <typename T, typename S, typename DESCRIPTOR> +bool PoiseuilleStrainRate2D<T,S,DESCRIPTOR>::operator()(T output[], const S input[]) +{ + T y = input[1]; + + T DuDx = T(); + T DuDy = (T) _maxVelocity*(-2.*(y-(_lengthY/2.))/((_lengthY/2.)*(_lengthY/2.))); + T DvDx = T(); + T DvDy = T(); + + output[0] = (DuDx + DuDx)/2.; + output[1] = (DuDy + DvDx)/2.; + output[2] = (DvDx + DuDy)/2.; + output[3] = (DvDy + DvDy)/2.; + + return true; +}; + +template <typename T> +AnalyticalPorousVelocity2D<T>::AnalyticalPorousVelocity2D(std::vector<T> axisDirection_, T K_, T mu_, T gradP_, T radius_, T eps_) + : AnalyticalF2D<T,T>(2), axisDirection(axisDirection_), K(K_), mu(mu_), gradP(gradP_),radius(radius_), eps(eps_) +{ + this->getName() = "AnalyticalPorousVelocity2D"; +}; + + +template <typename T> +T AnalyticalPorousVelocity2D<T>::getPeakVelocity() +{ + T uMax = K / mu*gradP*(1. - 1./(cosh((sqrt(1./K))*radius))); + + return uMax/eps; +}; + + +template <typename T> +bool AnalyticalPorousVelocity2D<T>::operator()(T output[], const T input[]) +{ + output[0] = K / mu*gradP*(1. - (cosh((sqrt(1./K))*(input[1] - radius)))/(cosh((sqrt(1./K))*radius))); + output[1] = K / mu*gradP*(1. - (cosh((sqrt(1./K))*(input[0] - radius)))/(cosh((sqrt(1./K))*radius))); + + output[0] *= axisDirection[0]/eps; + output[1] *= axisDirection[1]/eps; + + return true; +}; + + +////////// Polar & Cartesian /////////////// + + +// constructor to obtain Cartesian coordinates of polar coordinates, +// with _polarOrigin, in x-y-plane +template<typename T, typename S> +PolarToCartesian2D<T, S>::PolarToCartesian2D(std::vector<T> polarOrigin) + : AnalyticalF2D<T, S>(2), + _polarOrigin(polarOrigin) +{ + this->getName() = "const"; +} + +// operator returns Cartesian coordinates of polar coordinates, +// input is x[0] = radius , x[1] = phi +template<typename T, typename S> +bool PolarToCartesian2D<T, S>::operator()(T output[], const S x[]) +{ + output[0] = x[0]*cos(x[1]) + _polarOrigin[0]; + output[1] = x[0]*sin(x[1]) + _polarOrigin[1]; + return true; +} + + +// constructor to obtain polar coordinates of Cartesian coordinates +template<typename T, typename S> +CartesianToPolar2D<T, S>::CartesianToPolar2D(T cartesianOriginX, + T cartesianOriginY, T cartesianOriginZ, + T axisDirectionX, T axisDirectionY, T axisDirectionZ, + T orientationX, T orientationY, T orientationZ) + : AnalyticalF2D<T, S>(2) +{ + _cartesianOrigin.push_back(cartesianOriginX); + _cartesianOrigin.push_back(cartesianOriginY); + _cartesianOrigin.push_back(cartesianOriginZ); + + _axisDirection.push_back(axisDirectionX); + _axisDirection.push_back(axisDirectionY); + _axisDirection.push_back(axisDirectionZ); + + _orientation.push_back(orientationX); + _orientation.push_back(orientationY); + _orientation.push_back(orientationZ); +} + +template<typename T, typename S> +CartesianToPolar2D<T, S>::CartesianToPolar2D(std::vector<T> cartesianOrigin, + std::vector<T> axisDirection, + std::vector<T> orientation) + : AnalyticalF2D<T, S>(2), + _cartesianOrigin(cartesianOrigin), + _axisDirection(axisDirection), + _orientation(orientation) +{ + this->getName() = "const"; +} + +// operator returns polar coordinates, output[0] = radius(>= 0), +// output[1] = phi in [0, 2Pi) +template<typename T, typename S> +bool CartesianToPolar2D<T, S>::operator()(T output[], const S x[]) +{ + //Vector<T, 3> x_rot(x[this->getSourceDim()]); //// CAUSES ERROR + int dim = this->getSourceDim(); + //std::cout<<"Source dim return: "<<dim<<std::endl; //sourceDim = 2 + Vector<T, 3> x_rot; + for (int i = 0; i < dim; ++i) { + x_rot[i] = x[i]; + } + Vector<T, 3> e3(T(), T(), 1); + Vector<T, 3> axisDirection(_axisDirection); + + Vector<T, 3> normal = crossProduct3D(axisDirection, e3); + Vector<T, 3> normalAxisDir(axisDirection); + normalAxisDir.normalize(); + + // if axis has to be rotated + if (!( util::nearZero(normalAxisDir[0]) && util::nearZero(normalAxisDir[1]) && util::nearZero(normalAxisDir[2]-1) ) ) { + + if ( !util::nearZero(util::norm(_orientation)) ) { + normal = _orientation; + } + + // normal is orientation direction + AngleBetweenVectors3D<T, S> angle(fromVector3(e3), fromVector3(normal)); + T tmp[3] = {_axisDirection[0],_axisDirection[1],_axisDirection[2]}; + T alpha[1] = {}; + angle(alpha, tmp); + + // cross is rotation axis + //Vector<T, 3> cross = crossProduct3D(e3, axisDirection); + // rotation with angle alpha to rotAxisDir + RotationRoundAxis3D<T, S> rotRAxis(_cartesianOrigin, fromVector3(normal), -alpha[0]); + T x_tmp[3] = {}; + x_tmp[0] = x[0]; + x_tmp[1] = x[1]; + x_tmp[2] = x[2]; + T tmp2[3] = {}; + rotRAxis(tmp2,x_tmp); + x_rot[0] = tmp2[0]; + x_rot[1] = tmp2[1]; + x_rot[2] = tmp2[2]; + } + + // calculates r, phi related to cartesianAxisDirection and their origin + Vector<T, 2> difference; + difference[0] = x_rot[0] - _cartesianOrigin[0]; + difference[1] = x_rot[1] - _cartesianOrigin[1]; + T distance = T(); + + for (int i = 0; i < 2; ++i) { + distance += difference[i]*difference[i]; + } + + distance = sqrt(distance); + T phi[1] = {}; + + if (distance > T()) { + Vector<T, 3> e1(T(1), T(), T()); + AngleBetweenVectors3D<T, S> angle(fromVector3(e1), fromVector3(e3)); + + T x_help[3] = {difference[0],difference[1],0.}; + angle(phi, x_help); + } + + output[0] = distance; + output[1] = phi[0]; + return true; +} + +} // end namespace olb +#endif |