// SpaFEDTe, a Template based C++ library for creating 
// Discontinuous Finite Element Spaces,
// Copyright (C) 2012 Lorenzo Alessio Botti

/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */ 
/* License as published by the Free Software Foundation either */ 
/* version 3.0 of the License, or (at your option) any later version. */

/* This software 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 */
/* Lesser General Public License for more details. */

/* You should have received a copy of the GNU Lesser General Public */
/* License along with this software; if not, a copy of the full */
/* GNU Lesser General Public License can be found at */
/* http://www.gnu.org/licenses/. */

#ifndef __exact_solution_hpp__
#define __exact_solution_hpp__

#include "FELib.hpp"

struct ExactSolution2DRegular
{
  static constexpr double pi = 3.1415926535897932384626433832795029L;

  template<typename Derived> static double exact_solution (const Eigen::MatrixBase<Derived>& p) 
  {
      const double x = p(0);
      const double y = p(1);
      return sin(pi*x)*sin(pi*y);
  }
  
  template<typename Derived> static Eigen::Matrix<double,2,1> exact_gradient (const Eigen::MatrixBase<Derived>& p) 
  {
    const double x = p(0);
    const double y = p(1);
    Eigen::Matrix<double, 2, 1> res;
    res <<
    pi * cos(pi*x)*sin(pi*y),
    pi * sin(pi*x)*cos(pi*y);
    return res;
  }

  template<typename Derived> static double forcingTerm (const Eigen::MatrixBase<Derived>& p) 
  {
    const double x = p(0);
    const double y = p(1);
    return 2.*pi*pi*(sin(pi*x)*sin(pi*y));
  }
};

struct ExactSolution3DRegular
{
  template<typename Derived> static double exact_solution (const Eigen::MatrixBase<Derived>& p) 
  {
      const double x = p(0);
      const double y = p(1);
      const double z = p(2);
      return cos(x) * exp(y) * (1. - z);
  }
      
  template<typename Derived> static Eigen::Matrix<double,2,1> exact_gradient (const Eigen::MatrixBase<Derived>& p) 
  {
    const double x = p(0);
    const double y = p(1);
    const double z = p(2);
    Eigen::Matrix<double, 3, 1> res;
    res <<
      -sin(x) * exp(y) * (1. - z),
      cos(x) * exp(y) * (1. - z),
      -cos(x) * exp(y);
    return res;
  }
  
  template<typename Derived> static double forcingTerm (const Eigen::MatrixBase<Derived>& p) 
  {
    const double x = p(0);
    const double y = p(1);
    const double z = p(2);
    return 0;
  }
};

struct ExactSolution3DSingular
{
  static constexpr double pi = 3.1415926535897932384626433832795029L;
  template<typename Derived> static double exact_solution (const Eigen::MatrixBase<Derived>& p) 
  {
      const double x = p(0);
      const double y = p(1);
      const double z = p(2);
      double theta = atan2(y,x);
      if (theta < 0)
        theta += 2. * pi;
      return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
  }
      
  template<typename Derived> static Eigen::Matrix<double,2,1> exact_gradient (const Eigen::MatrixBase<Derived>& p) 
  {
    const double x = p(0);
    const double y = p(1);
    const double z = p(2);
    const double tt = 2./3.;
    const double ot = 1./3.;
    const double r2 = x*x + y*y;
    double theta = atan2(y,x);
    if (theta < 0)
      theta += 2. * pi;
    Eigen::Matrix<double, 3, 1> res;
    res <<
      tt*x*pow(r2,-tt)*sin(tt*theta) - pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*y/x/x,
      tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x,
      1;
    return res;
  }
  
  template<typename Derived> static double forcingTerm (const Eigen::MatrixBase<Derived>& p) 
  {
    const double x = p(0);
    const double y = p(1);
    const double z = p(2);
    return 0;
  }
};
#endif
