path: root/src/functors/analytical/frameChangeF3D.h

/*  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 "analyticalF.h"

/** \file
  This file contains two different classes of functors, in the FIRST part
  - for simulations in a rotating frame
  - different functors for
      velocity (3d, RotatingLinear3D),
      pressure (1d, RotatingQuadratic1D) and
      force    (3d, RotatingForceField3D)
  The functors return the displacement of a point x in a fixed amount of time.

  The ones in the SECOND part are useful to set Poiseuille velocity profiles on
  - pipes with round cross-section and
  - pipes with square-shaped cross-section.
*/

/** 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> class SuperGeometry3D;

// PART 1: /////////////////////////////////////////////////////////////////////
// Functors for rotating the coordinate system (velocity, pressure, force,...)

/**
  * This functor gives a linar profile for a given point x as it computes
  * the distance between x and the axis.
  *
  * The field in outcome is the velocity field of q rotating solid
  */

/// Functor with a linear profile e.g. for rotating velocity fields.
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;
};

/**
  * This functor gives a linar profile in an annulus for a given point x between the inner and outer radius as it computes
  * the distance between x and the inner and outer radius.
  *
  * The field in outcome is the velocity field of q rotating solid in an annulus
  */

/// Functor with a linear profile e.g. for rotating velocity fields.
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[]);
};


/**
  * This functor gives a parabolic profile for a given point x as it computes
  * the distance between x and the axis.
  *
  * This field is a scalar field, a vector with one component will be used
  */

/// Functor with a parabolic profile e.g. for rotating pressure fields.
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

/**
  * This functor returns a quadratic Poiseuille profile for use with a pipe with
  * round cross-section. It uses cylinder coordinates and is valid for the
  * entire length of the pipe.
  *
  * This functor gives a parabolic velocity profile for a given point x as it
  * computes the distance between x and the axis.
  *
  * The axis is set in the input with axisPoint and axisDirection. The axisPoint
  * can be the coordinate of any point where the axis passes.
  * axisDirection has to be normed to 1.
  * Once the axis is set in the middle of the pipe, the radius of the
  * pipe "radius" and the velocity in the middle of the pipe "maxVelocity"
  * determine the Poisseuille profile entierly.
  */

/// Velocity profile for round pipes and power law fluids: u(r)=u_max*(1-(r/R)^((n+1)/n)). The exponent n characterizes the fluid behavior.
/// n<1: Pseudoplastic, n=1: Newtonian fluid, n>1: Dilatant
template <typename T>
class CirclePowerLaw3D : public AnalyticalF3D<T,T> {
protected:
  std::vector<T> _center;
  std::vector<T> _normal;
  T _radius;
  T _maxVelocity;
  T _n;
  T _scale;

public:
  CirclePowerLaw3D(std::vector<T> 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(SuperGeometry3D<T>& superGeometry, int material, T maxVelocity, T n, 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, SuperGeometry3D<T>& superGeometry, int material, T Velocity, T n, T scale = T(1));

  /// Returns centerpoint vector
  std::vector<T> getCenter()
  {
    return _center;
  };
  /// Returns normal vector
  std::vector<T> getNormal()
  {
    return _normal;
  };
  /// Returns radi
  T getRadius()
  {
    return _radius;
  };

  bool operator()(T output[], const T x[]) override;
};

/// Velocity profile for round pipes and turbulent flows: u(r)=u_max*(1-r/R)^(1/n) The exponent n can be calculated by n = 1.03 * ln(Re) - 3.6
/// n=7 is used for many flow applications
template <typename T>
class CirclePowerLawTurbulent3D : public CirclePowerLaw3D<T> {
public:
  CirclePowerLawTurbulent3D(std::vector<T> axisPoint_, std::vector<T> axisDirection,  T maxVelocity, T radius, T n, T scale = T(1));
  CirclePowerLawTurbulent3D(T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T maxVelocity, T n, T scale = T(1));
  CirclePowerLawTurbulent3D(SuperGeometry3D<T>& superGeometry, int material, T maxVelocity, T n, T scale = T(1));

  CirclePowerLawTurbulent3D(bool useMeanVelocity, std::vector<T> axisPoint, std::vector<T> axisDirection,  T Velocity, T radius, T n, 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, T scale = T(1));
  CirclePowerLawTurbulent3D(bool useMeanVelocity, SuperGeometry3D<T>& superGeometry, int material, T Velocity, T n, T scale = T(1));

  bool operator()(T output[], const T x[]) override;
};

/// Velocity profile for round pipes and a laminar flow of a Newtonian fluid: u(r)=u_max*(1-(r/R)^2)
template <typename T>
class CirclePoiseuille3D final : public CirclePowerLaw3D<T> {

public:
  CirclePoiseuille3D(std::vector<T> axisPoint, std::vector<T> axisDirection,  T maxVelocity, T radiu