/** \file navierStokes.h */

#ifndef NAVIER_STOKES_H
#define NAVIER_STOKES_H

#include "AMDiS.h"
#include "Views.h"
#include "GeometryTools.h"

struct InflowBC : AbstractFunction<double, WorldVector<double> >
{
  InflowBC(double H_=4.1, double Um_=1.5) : H(H_), Um(Um_) {}
  double operator()(const WorldVector<double> &x) const {
    return 4.0 * Um * x[1] * (H - x[1]) / sqr(H);
  }
protected:
  double H;
  double Um;
};


class Polygon : public AbstractFunction<double, WorldVector<double> >
{
public:

  Polygon(WorldVector<double> x0_, WorldVector<double> x1_, WorldVector<double> x2_, WorldVector<double> x3_)
  {
    vertices.push_back(x0_);
    vertices.push_back(x1_);
    vertices.push_back(x2_);
    vertices.push_back(x3_);
    vertices.push_back(x0_);
  }

  Polygon(std::vector<WorldVector<double> > xi_) : vertices(xi_) { }

  double operator()(const WorldVector<double>& x) const
  {
    double result = 1.e15;
    for (size_t i = 0; i < vertices.size()-1; i++)
      result = std::min(result, meshconv::distance_point_line_2d(x.begin(), vertices[i].begin(), vertices[i+1].begin()));
    return result * (meshconv::point_in_polygon(x.begin(), vertices) ? -1.0 : 1.0);
  }
  
  void refine(int np)
  {
    std::vector<WorldVector<double> > newVertices;
    
    for (size_t i = 0; i < vertices.size()-1; i++) {
      for (size_t j = 0; j < np-1; j++) {
	double lambda = static_cast<double>(j)/static_cast<double>(np - 1.0);
	WorldVector<double> p = lambda*vertices[i+1] + (1.0-lambda)*vertices[i];
	newVertices.push_back(p);
      }
    }
    swap(vertices, newVertices);
  }
  
  void move(const DOFVector<WorldVector<double> >* velocity)
  {
    for (size_t i = 0; i < vertices.size()-1; i++) {
      WorldVector<double> shift = evalAtPoint(*velocity, vertices[i]);
      vertices[i] += shift;
    }
    vertices[vertices.size()-1] = vertices[0];
  }
  
  void move(AbstractFunction<WorldVector<double>, WorldVector<double> >* velocity)
  {
    for (size_t i = 0; i < vertices.size()-1; i++) {
      WorldVector<double> shift = (*velocity)(vertices[i]);
      vertices[i] += shift;
    }
    vertices[vertices.size()-1] = vertices[0];
  }

private:
  std::vector<WorldVector<double> > vertices;
};

#endif