#include <dolfin.h>
#include "Laplace.h"
#include "expression.h"
#include "functionspace.h"


class Source : public dolfin::Expression
{
  void eval(dolfin::Array<double>& values, const dolfin::Array<double>& x) const
  {
    double dx = x[0] - 0.5;
    double dy = x[1] - 0.5;
    values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
  }
};

class dUdN : public dolfin::Expression
{
  void eval(dolfin::Array<double>& values, const dolfin::Array<double>& x) const
  {
    values[0] = sin(5*x[0]);
  }
};


DEFUN_DLD (test_expr, args, , "test_expr: functionspace V, fcn_handle f, mesh_label label")
{

  int nargin = args.length ();
  octave_value retval=0;

  if (nargin < 3 || nargin > 3)
    print_usage ();
  else
    {
      if (args(0).type_id () == functionspace::static_type_id ())
        {
          const functionspace & fspo = static_cast<const functionspace&> (args(0).get_rep ());
          octave_fcn_handle * fh = args(1).fcn_handle_value ();
          Array<int> side = args(2).array_value ();

          if (!error_state)
            {
              const dolfin::FunctionSpace V = fspo.get_fsp ();
              dolfin::Mesh mesh = *(V.mesh());
              expression ff (*fh);

              //We use dolfin now, fem-fenics not yet available

              std::vector<const dolfin::BoundaryCondition*> bcu;
              Source f;
              dUdN g;

              for (octave_idx_type i = 0; i < side.length (); ++i)
                bcu.push_back (new dolfin::DirichletBC (V, ff, side(i)));


              Laplace::BilinearForm a (V, V);
              Laplace::LinearForm L (V);
              L.f = f;
              L.g = g;

              dolfin::Function u(V);
              dolfin::solve(a == L, u, bcu);

              dolfin::File file("fem-fenics.pvd");
              file << u;

              dolfin::plot(u);
              dolfin::interactive();
            }
        }
    }
  return retval;

}