frameChangeF3D.h

Go to the documentation of this file.

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2013, 2015 Gilles Zahnd, 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 FRAME_CHANGE_F_3D_H
#define FRAME_CHANGE_F_3D_H
 
#include<vector>
#include<string>
#include <random>
 
#include "analyticalF.h"
 
/* To enable simulations in a rotating frame, the axis is set in the
  * constructor with axisPoint and axisDirection. The axisPoint can be the
  * coordinate of any point on the axis. The axisDirection has to be a normed to
  * 1. The pulse w is in rad/s. It determines the pulse speed by its norm while
  * the trigonometric or clockwise direction is determined by its sign: When the
  * axisDirection is pointing "towards you", a positive pulse makes it turn in
  * the trigonometric way. It has to be noticed that putting both axisDirection
  * into -axisDirection and w into -w yields an exactly identical situation.
  */
 
namespace olb {
 
 
template<typename T, unsigned D> class SuperGeometry;
 
// PART 1: /////////////////////////////////////////////////////////////////////
// Functors for rotating the coordinate system (velocity, pressure, force,...)
 
template <typename T>
class RotatingLinear3D final : public AnalyticalF3D<T,T> {
protected:
  std::vector<T> axisPoint;
  std::vector<T> axisDirection;
  T w;
  T scale;
public:
  RotatingLinear3D(std::vector<T> axisPoint_, std::vector<T> axisDirection_, T w_, T scale_=1);
  bool operator()(T output[], const T x[]) override;
};
 
template <typename T>
class RotatingLinearAnnulus3D final : public AnalyticalF3D<T,T> {
protected:
  std::vector<T> axisPoint;
  std::vector<T> axisDirection;
  T w;
  T ri;
  T ro;
  T scale;
public:
  RotatingLinearAnnulus3D(std::vector<T> axisPoint_, std::vector<T> axisDirection_, T w_, T ri_, T ro_, T scale_=1);
  bool operator()(T output[], const T x[]);
};
 
 
template <typename T>
class RotatingQuadratic1D final : public AnalyticalF3D<T,T> {
protected:
  std::vector<T> axisPoint;
  std::vector<T> axisDirection;
  T w;
  T scale;
  T additive;
public:
  RotatingQuadratic1D(std::vector<T> axisPoint_, std::vector<T> axisDirection_,
                      T w_, T scale_=1, T additive_=0);
  bool operator()(T output[], const T x[]) override;
};
 
 
// PART 2: /////////////////////////////////////////////////////////////////////
// Functors for setting velocities on a velocity boundary of a pipe
 
template <typename T>
class CirclePowerLaw3D : public AnalyticalF3D<T,T> {
protected:
  olb::Vector<T, 3> _center;
  std::vector<T> _normal;
  T _radius;
  T _maxVelocity;
  T _n;
  T _scale;
 
public:
  CirclePowerLaw3D(olb::Vector<T, 3> axisPoint, std::vector<T> axisDirection,  T maxVelocity, T radius, T n, T scale = T(1));
  CirclePowerLaw3D(T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T maxVelocity, T n, T scale = T(1));
  CirclePowerLaw3D(const SuperGeometry<T,3>& superGeometry, int material, T maxVelocity, T n, T distance2Wall, T scale = T(1));
 
  CirclePowerLaw3D(bool useMeanVelocity, std::vector<T> axisPoint, std::vector<T> axisDirection,  T Velocity, T radius, T n, T scale = T(1));
  CirclePowerLaw3D(bool useMeanVelocity, T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T Velocity, T n, T scale = T(1));
  CirclePowerLaw3D(bool useMeanVelocity, SuperGeometry<T,3>& superGeometry, int material, T Velocity, T n, T distance2Wall, T scale = T(1));
 
  olb::Vector<T, 3> getCenter()
  {
    return _center;
  };
  std::vector<T> getNormal()
  {
    return _normal;
  };
  T getRadius()
  {
    return _radius;
  };
 
  bool operator()(T output[], const T x[]) override;
};
 
template <typename T>
class CirclePowerLawTurbulent3D : public CirclePowerLaw3D<T> {
private:
  T _turbulenceIntensity;
  std::random_device _rd;
  std::mt19937 _generator;
  std::normal_distribution<T> _dist;
 
public:
  CirclePowerLawTurbulent3D(std::vector<T> axisPoint_, std::vector<T> axisDirection,  T maxVelocity, T radius, T n = 7, T turbulenceIntensity = 0.05, T scale = T(1));
  CirclePowerLawTurbulent3D(T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T maxVelocity, T n = 7, T turbulenceIntensity = 0.05, T scale = T(1));
  CirclePowerLawTurbulent3D(SuperGeometry<T,3>& superGeometry, int material, T maxVelocity, T distance2Wall, T n = 7, T turbulenceIntensity = 0.05, T scale = T(1));
 
  CirclePowerLawTurbulent3D(bool useMeanVelocity, std::vector<T> axisPoint, std::vector<T> axisDirection,  T Velocity, T radius, T n = 7, T turbulenceIntensity = 0.05, T scale = T(1));
  CirclePowerLawTurbulent3D(bool useMeanVelocity, T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T Velocity, T n = 7, T turbulenceIntensity = 0.05, T scale = T(1));
  CirclePowerLawTurbulent3D(bool useMeanVelocity, SuperGeometry<T,3>& superGeometry, int material, T Velocity, T distance2Wall, T n = 7, T turbulenceIntensity = 0.05, T scale = T(1));
 
  bool operator()(T output[], const T x[]) override;
};
 
template <typename T>
class CirclePoiseuille3D final : public CirclePowerLaw3D<T> {
 
public:
  CirclePoiseuille3D(std::vector<T> axisPoint, std::vector<T> axisDirection,  T maxVelocity, T radius, T scale = T(1));
  CirclePoiseuille3D(T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T maxVelocity, T scale = T(1));
  CirclePoiseuille3D(SuperGeometry<T,3>& superGeometry, int material, T maxVelocity, T distance2Wall, T scale = T(1));
 
  CirclePoiseuille3D(bool useMeanVelocity, std::vector<T> axisPoint, std::vector<T> axisDirection,  T Velocity, T radius, T scale = T(1));
  CirclePoiseuille3D(bool useMeanVelocity, T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T Velocity, T scale = T(1));
  CirclePoiseuille3D(bool useMeanVelocity, SuperGeometry<T,3>& superGeometry, int material, T Velocity, T distance2Wall, T scale = T(1));
};
 
template < typename T,typename DESCRIPTOR>
class CirclePoiseuilleStrainRate3D : public AnalyticalF3D<T,T> {
protected:
  T lengthY;
  T lengthZ;
  T maxVelocity;
 
public:
  CirclePoiseuilleStrainRate3D(UnitConverter<T, DESCRIPTOR> const& converter, T ly);
  bool operator()(T output[], const T input[]) override;
};
 
template <typename T>
class RectanglePoiseuille3D final : public AnalyticalF3D<T,T> {
protected:
  OstreamManager clout;
  olb::Vector<T, 3> x0;
  olb::Vector<T, 3> x1;
  olb::Vector<T, 3> x2;
  std::vector<T> maxVelocity;
 
public:
  RectanglePoiseuille3D(olb::Vector<T, 3>& x0_, olb::Vector<T, 3>& x1_, olb::Vector<T, 3>& x2_, std::vector<T>& maxVelocity_);
 
  RectanglePoiseuille3D(SuperGeometry<T,3>& superGeometry_, int material_, std::vector<T>& maxVelocity_, T offsetX, T offsetY, T offsetZ);
  bool operator()(T output[], const T x[]) override;
};
 
//More accurate Version using a trigonometric formulation
template <typename T>
class RectangleTrigonometricPoiseuille3D final : public AnalyticalF3D<T,T> {
protected:
  OstreamManager clout;
  std::vector<T> x0;
  std::vector<T> x1;
  std::vector<T> x2;
  std::vector<T> maxVelocity;
  int numOfPolynoms;
  T profilePeak;
  T profileRatioAvgToPeak;
  void calcPeakAndAvg();
 
public:
  RectangleTrigonometricPoiseuille3D(std::vector<T>& x0_, std::vector<T>& x1_, std::vector<T>& x2_, std::vector<T>& maxVelocity, int numOfPolynoms_ = 1);
 
  RectangleTrigonometricPoiseuille3D(SuperGeometry<T,3>& superGeometry_, int material_, std::vector<T>& maxVelocity_, T offsetX, T offsetY, T offsetZ, int numOfPolynoms_ = 1);
  T getProfilePeak();
  T getRatioAvgToPeak();
  bool operator()(T output[], const T x[]) override;
};
 
template <typename T>
class EllipticPoiseuille3D final : public AnalyticalF3D<T,T> {
protected:
  OstreamManager clout;
  std::vector<T> _center;
  T _a2, _b2;
  T _maxVel;
 
public:
  EllipticPoiseuille3D(std::vector<T> center, T a, T b, T maxVel);
  bool operator()(T output[], const T x[]) override;
};
 
 
template <typename T>
class AnalyticalPorousVelocity3D : public AnalyticalF3D<T,T> {
protected:
  olb::Vector<T, 3> center;
  olb::Vector<T, 3> normal;
  T K, mu, gradP, radius;
  T eps;
public:
  AnalyticalPorousVelocity3D(SuperGeometry<T,3>& superGeometry, int material, T K_, T mu_, T gradP_, T radius_, T eps_=T(1));
  T getPeakVelocity();
  bool operator()(T output[], const T input[]) override;
};
 
 
 
 
template <typename T, typename S>
class AngleBetweenVectors3D final : public AnalyticalF3D<T,S> {
protected:
  std::vector<T> _referenceVector;
  std::vector<T> _orientation;
public:
  AngleBetweenVectors3D(std::vector<T> referenceVector,
                        std::vector<T> orientation);
  bool operator()(T output[], const S x[]) override;
};
 
 
template <typename T, typename S>
class RotationRoundAxis3D final : public AnalyticalF3D<T,S> {
protected:
  olb::Vector<T, 3>_origin;
  olb::Vector<T, 3> _rotAxisDirection;
  T _alpha;
public:
  RotationRoundAxis3D(olb::Vector<T, 3> origin, olb::Vector<T, 3> rotAxisDirection,
                      T alpha);
  bool operator()(T output[], const S x[]) override;
};
 
 
 
 
template <typename T, typename S>
class CylinderToCartesian3D final : public AnalyticalF3D<T,S> {
protected:
  std::vector<T> _cylinderOrigin;
public:
  CylinderToCartesian3D(std::vector<T> cylinderOrigin);
  bool operator()(T output[], const S x[]) override;
};
 
 
template <typename T, typename S>
class CartesianToCylinder3D final : public AnalyticalF3D<T,S> {
protected:
    olb::Vector<T, 3>_cartesianOrigin;
    olb::Vector<T, 3> _axisDirection;
    olb::Vector<T, 3> _orientation;
public:
  CartesianToCylinder3D(olb::Vector<T, 3> cartesianOrigin, T& angle,
                        olb::Vector<T, 3> orientation = {T(1),T(),T()});
  CartesianToCylinder3D(olb::Vector<T, 3> cartesianOrigin,
                        olb::Vector<T, 3> axisDirection,
                        olb::Vector<T, 3> orientation = {T(1),T(),T()});
  CartesianToCylinder3D(T cartesianOriginX, T cartesianOriginY,
                        T cartesianOriginZ,
                        T axisDirectionX, T axisDirectionY,
                        T axisDirectionZ,
                        T orientationX = T(1), T orientationY = T(),
                        T orientationZ = T());
  bool operator()(T output[], const S x[]) override;
  olb::Vector<T, 3> const& getAxisDirection() const;
  olb::Vector<T, 3>& getAxisDirection();
  olb::Vector<T, 3> const& getCartesianOrigin() const;
  olb::Vector<T, 3>& getCartesianOrigin();
  void setAngle(T angle);
};
 
 
 
 
template <typename T, typename S>
class SphericalToCartesian3D final : public AnalyticalF3D<T,S> {
  //protected:
public:
  SphericalToCartesian3D();
  bool operator()(T output[], const S x[]) override;
};
 
 
template <typename T, typename S>
class CartesianToSpherical3D final : public AnalyticalF3D<T,S> {
protected:
  std::vector<T> _cartesianOrigin;
  std::vector<T> _axisDirection;
public:
  CartesianToSpherical3D(std::vector<T> cartesianOrigin,
                         std::vector<T> axisDirection);//, std::vector<T> orientation);
  bool operator()(T output[], const S x[]) override;
};
 
 
 
 
template <typename T, typename S>
class MagneticForceFromCylinder3D final : public AnalyticalF3D<T,S> {
protected:
  CartesianToCylinder3D<T,S>& _car2cyl;
  T _length;
  T _radWire;
  T _cutoff;
  T _Mw;
  T _Mp;
  T _radParticle;
  T _mu0;
  T _factor;
public:
  MagneticForceFromCylinder3D(CartesianToCylinder3D<T,S>& car2cyl, T length,
                              T radWire, T cutoff, T Mw, T Mp = T(1),
                              T radParticle = T(1), T mu0 = 4*3.14159265e-7);
  bool operator()(T output[], const S x[]) override;
};
 
template <typename T, typename S>
class MagneticFieldFromCylinder3D final : public AnalyticalF3D<T, S> {
protected:
  CartesianToCylinder3D<T, S>& _car2cyl;
  T _length;
  T _radWire;
  T _cutoff;
  T _Mw;
  T _factor;
public:
  MagneticFieldFromCylinder3D(CartesianToCylinder3D<T, S>& car2cyl,
                              T length,
                              T radWire,
                              T cutoff,
                              T Mw
                             );
  bool operator()(T output[], const S x[]) override;
};
 
} // end namespace olb
 
#endif