#include "Laplace.h"
using namespace Laplace;

#include <dolfin.h>
#include "functionspace.h"
#include "boundarycondition.h"
#include "coefficient.h"

DEFUN_DLD (test_coeff, args, , "test_bc: functionspace V, boundarycondition bc")
{

  int nargin = args.length ();
  octave_value retval=0;
  
  if (nargin < 2)
    print_usage ();
  else
    {
      if (args(0).type_id () == functionspace::static_type_id () &&
          args(1).type_id () == boundarycondition::static_type_id ())
        {
          const functionspace & fspo =
            static_cast<const functionspace&> (args(0).get_rep ());
          const boundarycondition & bc =
            static_cast<const boundarycondition&> (args(1).get_rep ());

          if (!error_state)
            {
              const dolfin::FunctionSpace V = fspo.get_fsp ();
              LinearForm L (V);
              std::size_t ncoef = L.num_coefficients ();
              if (nargin < 2 + ncoef || nargin > 2 + ncoef)
                error ("Wrong number of coefficient");
              else
                {
                  for (octave_idx_type i = 0; i < ncoef; ++i)
                    {
                      const coefficient & cf =static_cast<const coefficient&> (args(i+2).get_rep ());
                      std::size_t n = L.coefficient_number (cf.get_str ());
                      const boost::shared_ptr<const expression> & pexp = cf.get_expr ();
                      L.set_coefficient (n, pexp);
                    }

                  const std::vector<boost::shared_ptr <const dolfin::DirichletBC> > &
                  pbc = bc.get_bc ();
                  std::vector<const dolfin::BoundaryCondition*> bcu;

                  for (octave_idx_type i = 0; i < pbc.size (); ++i)
                    bcu.push_back (& (* (pbc[i])));

                  //Now use dolfin, fem-fenics not yet available
                  dolfin::Vector b;
                  dolfin::assemble (b, L);

                  for (std::size_t i = 0; i < bcu.size(); i++)
                    bcu[i]->apply(b);

                  BilinearForm a (V, V);
                  dolfin::Matrix A;
                  dolfin::assemble (A, a);

                  for (std::size_t i = 0; i < bcu.size(); i++)
                    bcu[i]->apply(A);

                  dolfin::Function u(V);
                  dolfin::solve(A, *u.vector(), b, "gmres", "default");

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

                  dolfin::plot (u);
                  dolfin::interactive ();

                }
            }
        }
    }
  return retval;

}