// -*- C++ -*-
#ifndef PHOSOLUTION_H
#define PHOSOLUTION_H

#include "CommonTemplateTypes.hpp"
#include "InputData.hpp"


namespace ho
{
  struct ExactSolution
  {
    typedef CommonTemplateTypes::ExactSolutionType ExactSolutionType;
    typedef CommonTemplateTypes::DiffusivityValueType DiffusivityValueType;
    typedef CommonTemplateTypes::DiffusivityType DiffusivityType;
    typedef CommonTemplateTypes::LoadType LoadType;
    typedef CommonTemplateTypes::ExactGradientType ExactGradientType;
    typedef CommonTemplateTypes::Point Point;
  
    DiffusivityType nu;
    int degreeNu;
    LoadType f;
  
    ExactSolutionType u;
    ExactGradientType G;
  };
  
  //------------------------------------------------------------------------------
  //------------------------------------------------------------------------------
  
  struct Regular : public ExactSolution
  {
  
    Regular(const double & _nu = 1.)
    {
      degreeNu = 0;
  
      nu = [_nu](const Point & x) -> DiffusivityValueType {
        return _nu * DiffusivityValueType::Identity();
      };
  
      u = [](const Point & x) -> double {
        return sin(InputData::pi*x(0))*sin(InputData::pi*x(1));
      };
  
      G = [](const Point & x) -> Eigen::Matrix<double, 2, 1> {
        Eigen::Matrix<double, 2, 1> res;
        res <<
        InputData::pi * cos(InputData::pi*x(0))*sin(InputData::pi*x(1)),
        InputData::pi * sin(InputData::pi*x(0))*cos(InputData::pi*x(1));
        return res;
      };
  
      f = [this,_nu](const Point & x) -> double {
        return 2.*_nu*std::pow(InputData::pi, 2)*sin(InputData::pi*x(0))*sin(InputData::pi*x(1));
      };
    }
  
  };
  
  //------------------------------------------------------------------------------
  //------------------------------------------------------------------------------
  
  struct ExactNeumannOmo : public ExactSolution
  {
  
    ExactNeumannOmo(const double & _nu = 1.)
    {
      degreeNu = 0;
  
      nu = [_nu](const Point & x) -> DiffusivityValueType {
        return _nu * DiffusivityValueType::Identity();
      };
  
      u = [](const Point & x) -> double {
        return pow(x(0),2) * (3. - 2. * x(0)) - pow(x(1),2) * (3. - 2. * x(1));
      };
  
      G = [](const Point & x) -> Eigen::Matrix<double, 2, 1> {
        Eigen::Matrix<double, 2, 1> res;
        res <<
        6. * (x(0) - pow(x(0),2)),
        6. * (pow(x(1),2) - x(1));
        return res;
      };
  
      f = [this,_nu](const Point & x) -> double {
        return 12. * (x(0) - x(1));
      };
    }
  
  };
  
  //------------------------------------------------------------------------------
  //------------------------------------------------------------------------------
  
  struct ExactNeumannAff : public ExactSolution
  {
  
    ExactNeumannAff(const double & _nu = 1.)
    {
      degreeNu = 0;
  
      nu = [_nu](const Point & x) -> DiffusivityValueType {
        return _nu * DiffusivityValueType::Identity();
      };
  
      u = [](const Point & x) -> double {
        return x(0) - x(1);
      };
  
      G = [](const Point & x) -> Eigen::Matrix<double, 2, 1> {
        Eigen::Matrix<double, 2, 1> res;
        res <<
        1.,
        - 1.;
        return res;
      };
  
      f = [this,_nu](const Point & x) -> double {
        return 0.;
      };
    }
  
  };
  
  //------------------------------------------------------------------------------
  //------------------------------------------------------------------------------
  
  struct Anisotropic : public ExactSolution
  {
  
    Anisotropic(const double & _epsilon = 1.)
    {
      degreeNu = 0;
  
      nu = [_epsilon](const Point & x) -> DiffusivityValueType {
        DiffusivityValueType res;
        res << 1., 0., 0., _epsilon;
        return res;
      };
  
      u = [](const Point & x) -> double {
        return sin(InputData::pi*x(0))*sin(InputData::pi*x(1));
      };
  
      G = [](const Point & x) -> Eigen::Matrix<double, 2, 1> {
        Eigen::Matrix<double, 2, 1> res;
        res <<
        InputData::pi * cos(InputData::pi*x(0))*sin(InputData::pi*x(1)),
        InputData::pi * sin(InputData::pi*x(0))*cos(InputData::pi*x(1));
        return res;
      };
  
      f = [this,_epsilon](const Point & x) -> double {
        return (1.+_epsilon)*std::pow(InputData::pi, 2)*sin(InputData::pi*x(0))*sin(InputData::pi*x(1));
      };
    }
  
  };
  
  //------------------------------------------------------------------------------
  //------------------------------------------------------------------------------
  
  struct LePotier : public ExactSolution
  {
  
    LePotier(const double & _epsilon = 1.)
    {
      degreeNu = 2;
  
      nu = [_epsilon](const Point & x) -> DiffusivityValueType {
        DiffusivityValueType res;
        res << 
        std::pow(x(1)+0.1,2) + _epsilon * std::pow(x(0)+0.1,2), -(1.-_epsilon)*(x(0)+0.1)*(x(1)+0.1), 
        -(1.-_epsilon)*(x(0)+0.1)*(x(1)+0.1), std::pow(x(0)+0.1,2) + _epsilon * std::pow(x(1)+0.1,2);
        return res;
      };
  
      u = [](const Point & x) -> double {
        return sin(InputData::pi*x(0))*sin(InputData::pi*x(1));
      };
  
      G = [](const Point & x) -> Eigen::Matrix<double, 2, 1> {
        Eigen::Matrix<double, 2, 1> res;
        res <<
        InputData::pi * cos(InputData::pi*x(0))*sin(InputData::pi*x(1)),
        InputData::pi * sin(InputData::pi*x(0))*cos(InputData::pi*x(1));
        return res;
      };
  
      f = [this,_epsilon](const Point & x) -> double {
        return 
        std::pow(InputData::pi,2)*(1.+_epsilon)*sin(InputData::pi*x(0))*sin(InputData::pi*x(1))*(std::pow(x(0)+0.1,2)+std::pow(x(1)+0.1,2))
        + InputData::pi*(1.-3.*_epsilon)*cos(InputData::pi*x(0))*sin(InputData::pi*x(1))*(x(0)+0.1)
        + InputData::pi*(1.-3.*_epsilon)*sin(InputData::pi*x(0))*cos(InputData::pi*x(1))*(x(1)+0.1)
        + 2.*std::pow(InputData::pi,2)*(1.-_epsilon)*cos(InputData::pi*x(0))*cos(InputData::pi*x(1))*(x(0)+0.1)*(x(1)+0.1);
      };
    }
  };
  
  //------------------------------------------------------------------------------
  //------------------------------------------------------------------------------
  
  struct MixedSecondOrder : public ExactSolution
  {
  
    MixedSecondOrder(const double & _nu = 1.)
    {
      degreeNu = 0;
  
      nu = [_nu](const Point & x) -> DiffusivityValueType {
        return _nu * DiffusivityValueType::Identity();
      };
  
      u = [](const Point & x) -> double {
        return x(0)*x(1);
      };
  
      G = [](const Point & x) -> Eigen::Matrix<double, 2, 1> {
        Eigen::Matrix<double, 2, 1> res;
        res << x(1), x(0);
  
        return res;
      };
  
      f = [this,_nu](const Point & x) -> double {
        return 0.;
      };
    }
  
  };

}  
  #endif
