Helpers.h

/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


#ifndef EXTENSIONS_HELPERS_H
#define EXTENSIONS_HELPERS_H

// #include "AMDiS.h"

#include "FixVec.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "DOFIterator.h"

#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/StdMpi.h"
#endif

#ifdef _WIN32
#include <io.h>
#else
#include <unistd.h>
#endif
#include <ios>
#include <iostream>
#include <fstream>
#include <string>

#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>
#include "boost/lexical_cast.hpp"
#include "boost/date_time/posix_time/posix_time.hpp"
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/special_functions/atanh.hpp>

#include <time.h>
#include <stdarg.h>

#include "VectorOperations.h"
#include "Views.h"
#include "operations/norm.hpp"
#include "operations/product.hpp"

#define TEST_WARNING(test) if ((test));else WARNING

namespace AMDiS { namespace extensions {
  
using namespace std;

static long random_seed_initial_value = 0;

namespace Helpers {

  // math routines
  // ==============
  
  inline double cint(double x) {
    double intpart; 
    if (modf(x, &intpart) >= 0.5)
      return x >=0 ? std::ceil(x) : std::floor(x);
    else
      return x < 0 ? std::ceil(x) : std::floor(x);
  }
  
  inline double round (double r, double places) {
    double off = std::pow(10.0, places);
    return cint(r*off)/off;
  }
  
  inline double signum(double val, double tol = DBL_TOL) { return val > tol ? 1.0 : (val < -tol ? -1.0 : 0.0); }
  inline int signum(int val) { return val > 0 ? 1 : (val < 0 ? -1 : 0); }
  inline int signum(bool val) { return val ? 1 : -1; }
  inline int toInt(bool val) { return val ? 1 : 0; }

  inline bool toBool(double val, double tol = DBL_TOL) { return val > tol || val < -tol ? true : false; }
  inline bool isPositiv(double val, double tol = DBL_TOL) { return val > -tol; }
  inline bool isStrictPositiv(double val, double tol = DBL_TOL) { return val > tol; }

  template<typename T>
  T atanh(T x)
  {
    T result;
    try {
      result = boost::math::atanh(x);
    } catch(...) { // domain_error, or overflow_error
      if (x > 0.0)
	result = 10.0;
      else
	result = -10.0;
    }
    return result;
  }
  
  // routines for conversion to string
  // =================================
  
  template<typename T>
  inline std::string toString(const T &value, 
			      ios_base::fmtflags flag = ios_base::scientific)
  {
    std::stringstream ss;
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    return ss.str();
  }
  
  
  inline std::string toString(const int &value)
  {
    return boost::lexical_cast<std::string>(value);
  }
  
  
  template<typename T>
  inline std::string toString(const WorldVector<T> &value, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
  {
    std::stringstream ss;
    if (brackets) ss<<"[";
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    
    if (brackets) ss<<"]";
    return ss.str();
  }
  
  
  template<typename T>
  inline std::string toString(const std::vector<T> &vec, 
			      bool brackets = true, 
			      std::string separator = ",",
			      ios_base::fmtflags flag = ios_base::scientific,
			      size_t precision = 6
 			    )
  {
    std::stringstream ss;
    if (brackets) ss<<"[";
    ss.setf(flag);
    ss.precision(precision);
    for (size_t i = 0; i < vec.size(); ++i) {
      if (!(ss<<(i>0?(separator+" "):std::string(""))<<vec[i]))
	throw(std::runtime_error("toString() not possible"));
    }
    if (brackets) ss<<"]";
    return ss.str();
  }
  
  
  template<class T>
  inline T fromString(const std::string& s)
  {
    std::istringstream stream(s);
    T t;
    stream >> t;
    return t;
  }

  
  inline std::string fillString(int length, char c, int numValues, ...)
  {
    va_list values;
    int value;
    va_start(values, numValues);

    int len = length;
    for(int i=0; i<numValues; i++) {
      value = va_arg(values, int);
      double tmp = static_cast<double>(value);
      while (true) {
	len --;
	tmp /= 10.0;
	if (tmp < 1.0 - DBL_TOL)
	  break;
      }
    }
    va_end(values);

    return std::string(len, c);
  }

  
  /// produce printable string of memory
  inline std::string memoryToString(double mem, int startIdx = 0) {
    int idx = startIdx;
    double mem_ = mem;
    while (mem_/1024.0 > 1.0) {
      mem_ /= 1024.0;
      idx++;
    }
    std::string result = toString(mem_)+" "+(idx==0?"B":(idx==1?"KB":(idx==2?"MB":"GB")));
    return result;
  }
  
  
  // some mesh routines
  // ===========================
  
  /// rotate scale and shift the mesh coordinates
  void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);

  /// scale a vector of meshes by different values in all directions
  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
  
  /// scale and shift by different values in all directions
  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
  
  /// scale by different values in all directions
  void scaleMesh(Mesh *mesh, WorldVector<double> scale);
  
  /// scale and shift by the same values in all directions
  void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
  
  /**
   * calculate the dimension of a mesh, by mesh traversal.
   * The mthod determines the maximal mesh coordinates and assumes symmetry of the
   * mesh around the origin, i.e. mesh = [-rsult, result]
   **/
  WorldVector<double> getMeshDimension(Mesh *mesh);
  
  /// calculate the dimension of a mesh, by mesh traversal.
  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
  
  inline double calcMeshSizes(Mesh* mesh, double& minH, double& maxH, int& minLevel, int& maxLevel) 
  {
    FixVec<WorldVector<double>, VERTEX> coords(mesh->getDim(), NO_INIT);

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    minH = 1e15; maxH = 0.0;
    int k = 0;
    minLevel = 100;
    maxLevel = 0;
    while (elInfo) {
      maxLevel = std::max(maxLevel,elInfo->getLevel());
      minLevel = std::min(minLevel,elInfo->getLevel());
      coords = elInfo->getCoords();
      double h = 0.0;
      for (int i = 0; i < coords.getSize(); i++) {
        for (int j = 0; j < coords.getSize(); j++) {
	  if (i != j)
            h = std::max(h,norm(coords[i]-coords[j]));
        }
      }
      minH = std::min(h, minH);
      maxH = std::max(h, maxH);
      elInfo = stack.traverseNext(elInfo);
      k++;
    }
    minLevel += mesh->getMacroElementLevel();
    maxLevel += mesh->getMacroElementLevel();
    
    return minH;
  }
  
  /// read DOFVector from AMDiS .dat-files
  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);

  // some linear algebra methods
  // ===========================

  /// calculate the determant of a 3x3 matrix of type dense2D
  inline double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
  {
    if(num_rows(A)==3 && num_cols(A)==3) {
      double det = A(0,0)*A(1,1)*A(2,2)+A(0,1)*A(1,2)*A(2,0)+A(0,2)*A(1,0)*A(2,1);
      return det-(A(0,2)*A(1,1)*A(2,0)+A(0,1)*A(1,0)*A(2,2)+A(0,0)*A(1,2)*A(2,1));
    }
    return 1.0;
  }
  
  
  /// calculate the determant of a 3x3 matrix of type WorldMatrix
  inline double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
  {
    if(A.getNumRows()==3 && A.getNumCols()==3) {
      double det = A[0][0]*A[1][1]*A[2][2]+A[0][1]*A[1][2]*A[2][0]+A[0][2]*A[1][0]*A[2][1];
      return det-(A[0][2]*A[1][1]*A[2][0]+A[0][1]*A[1][0]*A[2][2]+A[0][0]*A[1][2]*A[2][1]);
    }
    return 1.0;
  }
  
  
  /// inverse power method to find minimal eigenvalue of matrix A and appropriate eigenvector x, using m iterations
  template <class LinearSolver, class Matrix, class EigenVector>
  inline double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
  {
    using mtl::two_norm; using AMDiS::two_norm;
    using mtl::dot; using AMDiS::dot;

    EigenVector y( size(x) ), res( size(x) );
    double lambda = 0.0, lambda_old = 0.0, relErr;
    mtl::seed<double> seed;
    random(x, seed);
    x /= two_norm(x);

    solver.setMultipleRhs(true);
    for (int i = 0; i < m; ++i) {
      solver.solveSystem(A, y, x ); // solve Ay=x --> y
      lambda_old = lambda;
      lambda = dot(y,x) / dot(y,y);
      relErr = abs((lambda - lambda_old) / lambda);
      if (relErr < 1.e-5)
	break;
      x = y / two_norm(y);
    }
    solver.setMultipleRhs(false);

    return lambda;
  }

  
  /// calculate approximation to condition number of matrix A using the LinearSolver solver
  template <class LinSolver, class Matrix>
  inline double condition (LinSolver& solver, const Matrix& A, int m=10)
  {
    using mtl::two_norm; using AMDiS::two_norm;
    using mtl::dot; using AMDiS::dot;
    mtl::dense_vector<double>   x(num_rows(A)), y(num_rows(A));
    mtl::seed<double>           seed;
    random(x, seed);

    double lambda = 0.0, lambda_old = 0.0, relErr, result1;
    x /= two_norm(x);
    for (int i = 0; i < m; ++i) {
      y = A * x;
      lambda_old = lambda;
      lambda = dot(y,x);
      relErr = std::abs((lambda-lambda_old)/lambda);
      if (relErr < 1.e-5)
	break;
      x = y / two_norm(y) ;
    }
    result1 = std::abs(lambda / inverse_iteration(solver, A, x, m));
    return result1;
  }
  
  
  /// interpolate values of DOFVector along line from (x1,y1) -> (x2,y1) with nPoints points
  template<typename Container>
  inline void interpolOverLine(const Container &vec, 
			       double x1, double y1, double x2, double y2, 
			       int nPoints,
			       std::vector<WorldVector<double> > &x, std::vector<double> &y)
  { FUNCNAME("Helpers::interpolOverLine()");

    WorldVector<double> p;

    TEST_EXIT(nPoints > 1)("Zu wenig Punkte fuer eine Interpolation!\n");
    
    for (int i = 0; i < nPoints; ++i) {
      double lambda = static_cast<double>(i)/static_cast<double>(nPoints - 1.0);
      p[0] = lambda*x2 + (1.0-lambda)*x1;
      p[1] = lambda*y2 + (1.0-lambda)*y1;
      
      double value = evalAtPoint(vec, p);
      x.push_back(p);
      y.push_back(value);
    }
  }
  

  /// calculate maxima of DOFVector along x-axis, using interpolOverLine
  void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
  
  /// calculate maxima of DOFVector along y-axis, using interpolOverLine
  void calcMaxOnYAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);

  /// calc normal vectors of surface from element normals by averaging
  void getNormalsWeighted(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals);
  
  /// calc normal vectors of surface from element normals by local approximation by quartic
  void getNormals(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals, DOFVector<WorldMatrix<double> > *gradNormals = NULL);
  
  /// calc curvature from given normal vectors
  void getCurvature(DOFVector<WorldVector<double> >* normals, DOFVector<double>* curvature);
  
  // misc routines
  // =============
  
  void plot(std::vector<double> &values, int numRows = 7, int numCols = 20, std::string symbol = "*");
  
  inline long getMicroTimestamp() 
  {
    using namespace boost::posix_time;
    ptime t0(min_date_time);
    ptime now = microsec_clock::local_time();
    time_duration td = now-t0;
    return td.total_microseconds();
  }
  
  
  inline long getRandomSeed(int param = 0) 
  {
    using namespace boost::posix_time;
    ptime t0(min_date_time);
    ptime now = microsec_clock::local_time();
    time_duration td = now-t0;
    long value0 = td.total_microseconds();
    long value1 = clock();
    long value2 = random_seed_initial_value++;

    return value0 + value1 + value2 + param;
  }

  
} // end namespace Helpers

} }

using namespace AMDiS::extensions;

#endif // EXTENSIONS_HELPERS_H