numeric.h

Go to the documentation of this file.

#ifndef NUMERIC_H
#define NUMERIC_H
#include "dimensionality.h"
#include <fstream>
#include <cmath>
#include <cassert>
#include <valarray>
#include <iostream>
//---------------------------------------------------------------------------------------------------------------------

inline double frand() {return rand()/double(RAND_MAX);}
inline double srand() {    // Give a random number of 1 or -1
  double dr;
  double fr = frand();
  if (fr < 0.5) dr = -1.0; else dr = 1.0;
  return dr;
}
template<typename T>
inline T sqr(T val) { 
  return val*val;
}
template<typename T>
inline T max(T a, T b) { 
  return a>b?a:b;
}

//---------------------------------------------------------------------------------------------------------------------
#ifndef ONE_DIMENSIONAL 
template<class T>
class grid {
 protected:
  T *G;
  GridPosition size;
  int gsize;
 public:
  grid() : G(0), size(Position(0)), gsize(0) { }
  grid(GridPosition size_) : G(0), size(size_), gsize((size_+UnitGrid).volume()) {newData();}
  // Copy constructor
  grid(const grid &agrid) : G(0), size(agrid.size), gsize(agrid.gsize) {
    newData();
    for (int i=0; i<gsize; ++i) G[i] = agrid.G[i];
  }
  ~grid() {
    static int ic = 0;
    static int im = 3;
    ic += 1;
    if (ic==im) {
      //std::cout << "Hier kommt jetzt der Destructor!  " << ic << std::endl;
      im += 1;
      if (G) delete[] G;
      //std::cout << "Fertig!" << std::endl;
    }
  }
  const GridPosition& Size () const {return size;}
        
  void Resize (const GridPosition &asize) {
    size = asize; 
    gsize = (size+UnitGrid).volume();
    newData();
  }
        
  const grid& operator=(const std::valarray<T> arr) {
    assert (arr.size()==G.size()); 
    for (int i=0; i<gsize; ++i) G[i] = arr[i];
    return *this;
  }
  const grid& operator=(const T val) {
    for (int i=0; i<gsize; ++i) G[i] = val;
    return *this;
  }

  const grid& operator+=(const std::valarray<T> arr) {
    assert (arr.size()==G.size());
    for (int i=0; i<gsize; ++i) G[i] += arr[i];
    return *this;
  }
  const grid& operator+=(const grid& arr) {
    assert (arr.size==size);
    for (int i=0; i<gsize; ++i) G[i] += arr.G[i];
    return *this;
  }
  const grid& operator-=(const std::valarray<T> arr) {
    assert (arr.size()==G.size());
    for (int i=0; i<gsize; ++i) G[i] -= arr[i];
    return *this;
  }
  const grid& operator-=(const grid& arr) {
    assert (arr.size==size);
    for (int i=0; i<gsize; ++i) G[i] -= arr.G[i];
    return *this;
  }
  const grid& operator*=(const T val) {
    for (int i = 0; i < gsize; ++i) G[i] *= val;
    return *this;
  }
  const grid& operator/=(const T val) {
    for (int i=0; i<gsize; ++i) G[i] /= val;
    return *this;
  }
  //operator std::valarray<T>& () {    // Get a reference of the valarray (dirty!)
  //  return G;
  //}
  T& operator[] (const GridPosition &i) {
    int pos = i.arraypos(size);
    //if ((pos >= gsize) || (pos < 0)) {
    //  std::cerr << i << " Index out of bounds (1)" << std::endl;
    //  exit(-1);
    //}
    return G[pos];
  }
  // Returns reference/value at grid position
  T operator[] (const GridPosition &i) const {
    int pos = i.arraypos(size);
    //if ((pos>=gsize) || (pos<0)) {
    //  std::cerr << i << " Index out of bounds (2)" << std::endl;
    //  exit(-1);
    //}
    return G[pos];
  }
  // The linearly interpolated value at a position between the grid points
  T operator() (const Position &X) const {
    GridPosition p = GridPosition(X);
    p.limitto(size - 1);
    GridPosition pp = p + 1;
    Position a = Position(pp) - X;
    Position b = X - Position(p);
    return WeightedSum(*this,p,pp,a,b);
  }
  // Returns the derivative. Assumes that type T can be cast to a double
  grid<Position> derive(Position ddx);
  // Sets the values in the grid with the use of a function that takes no arguments, eg frand()
  void apply(T fkt());
 private:
  void newData() {
    if (G) delete[] G;
    G = new T[gsize];
  }
};
template<class T>
void grid<T>::apply(T fkt()) {
  for (int i = 0; i < gsize; i++) {
    T tmp = fkt();
    G[i] = tmp;
  }
}
#endif
//---------------------------------------------------------------------------------------------------------------------
#ifdef ONE_DIMENSIONAL 
template<class T>
class grid {
 protected:
  std::valarray<T> G;
  GridPosition size;
 public:
  grid() : G(), size(Position(0)) { }
  grid(GridPosition size_) : G(0.0, size_.volume()), size(size_) { }
  grid(const grid &agrid) : G(agrid.G), size(agrid.size) { }
  ~grid() { }
  const GridPosition& Size () const {return size;}

  void Resize (const GridPosition &asize) {
    size = asize;
    G.resize((size+UnitGrid).volume());
  }
     
  const grid& operator=(const std::valarray<T> arr) {
    assert (arr.size()==G.size());
    G = arr;
    return *this;
  }
  const grid& operator=(const T val) {
    G = val;
    return *this;
  }
        
  const grid& operator+=(const std::valarray<T> arr) {
    assert (arr.size()==G.size());
    G += arr;
    return *this;
  }       
  const grid& operator-=(const std::valarray<T> arr) {
    assert (arr.size()==G.size());
    G -= arr;
    return *this;
  }             
  const grid& operator*=(const T val) {
    G *= val;
    return *this;
  }     
  const grid& operator/=(const T val) {
    G /= val;
    return *this;
  }

                        
 operator std::valarray<T>& () {    // Get a reference of the valarray (dirty!)
   return G;
 }
  
 T& operator[] (const GridPosition &i) {
   int pos = i.arraypos(size);
   return G[i.arraypos(pos)];
 }
 // Returns reference/value at grid position
 T operator[] (const GridPosition &i) const {
   int pos = i.arraypos(size);
   return G[i.arraypos(pos)];
 }      

 // The linearly interpolated value at a position between the grid
 T operator() (const Position &X) const {
   GridPosition p = GridPosition(X);
   p.limitto(size - 1);
   GridPosition pp = p + 1;
   Position a = Position(pp) - X;
   Position b = X - Position(p);
   return WeightedSum(*this,p,pp,a,b);
 }
 // Returns the derivative. Assumes that type T can be cast to a double
 grid<Position> derive(Position ddx);
 // Sets the values in the grid with the use of a function that takes no arguments, eg frand()
 void apply(T fkt());
};
template<class T>
void grid<T>::apply(T fkt()) {
  for (int i=0; i<G.size(); i++) {
    T tmp = fkt();
    G[i] = tmp;
  }
}
template<class T>
grid<Position> grid<T>::derive(Position ddx) {
  int i;
  grid<Position> D;
  GridPosition NG = Size();
  D.Resize(NG);
  ddx *= 2.0;
  grid<T> &G = *this;
  // Derivative in  the volume
  #pragma omp parallel for private(i) default(shared)
  for (i=NG[0]-2; i>0; i--) D[i] = (G[i-1] - G[i+1]) / ddx[0];
  // Derivative at the boundaries
  D[0]       = 2.0 * (G[0] - G[1]) / ddx[0];                // First cell
  D[NG[0]-1] = 2.0 * (G[NG[0]-2] - G[NG[0]-1]) / ddx[0];    // Last cell
  return D;
}

void TriMatrixSolver(std::valarray<double> &a, std::valarray<double> &b, std::valarray<double> &c, std::valarray<double> &r, std::valarray<double> &u, int nt);
template<class T>
T WeightedSum(const grid<T> &g, GridPosition p1, GridPosition p2, const Position& wa, const Position& wb) {
  return (g[p1]*wa.P[0] + g[p2]*wb.P[0]);
}
#endif
//---------------------------------------------------------------------------------------------------------------------
#ifdef TWO_DIMENSIONAL 
template<class T>
T WeightedSum(const grid<T> &g,GridPosition p1, GridPosition p2, const Position& wa, const Position& wb) {
  T sum = g[p1]*wa.P[0]*wa.P[1] + g[p2]*wb.P[0]*wb.P[1];
  ++p1.P[0];
  --p2.P[0];
  return (sum + g[p1]*wb.P[0]*wa.P[1] + g[p2]*wa.P[0]*wb.P[1]);
}
#endif
//---------------------------------------------------------------------------------------------------------------------
#ifdef THREE_DIMENSIONAL 
template<class T>
T WeightedSum(const grid<T> &g,GridPosition p1, GridPosition p2,const Position& wa, const Position& wb) {
  T sum = g[p1]*wa.P[0]*wa.P[1]*wa.P[2] + g[p2]*wb.P[0]*wb.P[1]*wb.P[2];
  ++p1.P[0];
  --p2.P[0]; 
  sum += g[p1]*wb.P[0]*wa.P[1]*wa.P[2] + g[p2]*wa.P[0]*wb.P[1]*wb.P[2];
  ++p1.P[1];
  --p2.P[1];
  sum += g[p1]*wb.P[0]*wb.P[1]*wa.P[2] + g[p2]*wa.P[0]*wa.P[1]*wb.P[2];
  --p1.P[0];
  ++p2.P[0];
  return sum + g[p1]*wa.P[0]*wb.P[1]*wa.P[2] + g[p2]*wb.P[0]*wa.P[1]*wb.P[2];
}
#endif
//---------------------------------------------------------------------------------------------------------------------
#endif