/*
 Copyright (C) 2013 Marco Vassallo

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

 This program 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 General Public License for more
 details.

 You should have received a copy of the GNU General Public License along with
 this program; if not, see <http://www.gnu.org/licenses/>.
*/

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

DEFUN_DLD (fem_lhs, args, , "-*- texinfo -*-\n\
@deftypefn {Function File} {[@var{bc}]} = \
fem_rhs (@var{Functional Space}, @var{Coefficient}, \
@var{Boundary Condition}) \n\
The input parameters are\n\
@itemize @bullet \n\
@item @var{Functional Space} is the fem-fenics functional space where\
the bilinear form is defined\n\
@item @var{Boundary Condition} contains the value of the essential bc that\
we want to apply to the Bilinear Form\n\
@item @var{Coefficient} is a variable of type coefficient which contains\
the value of the coefficient for the bilinear form\n\
@end itemize\n\
The output @var{A} is a sparse matrix which represents the bilinear form\n\
@seealso{fem_init_mesh, fem_fs}\n\
@end deftypefn")
{

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

  if (nargin < 1)
    print_usage ();
  else
    {

      if (! functionspace_type_loaded)
        {
          functionspace::register_type ();
          functionspace_type_loaded = true;
          mlock ();
        }
      if (args(0).type_id () == functionspace::static_type_id ())
        {
          const functionspace & fspo
            = static_cast<const functionspace&> (args(0).get_rep ());

          if (! error_state)
            {
              const dolfin::FunctionSpace V = fspo.get_fsp ();
              Laplace::LinearForm L (V);
              std::size_t ncoef = L.num_coefficients (), nc = 0;

              if (! coefficient_type_loaded)
                {
                  coefficient::register_type ();
                  coefficient_type_loaded = true;
                  mlock ();
                }
              for (std::size_t i = 1; i < nargin; ++i)
                {
                  if (args(i).type_id () == coefficient::static_type_id ())
                    {
                      const coefficient & cf
                        = static_cast <const coefficient&> (args(i).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);
                      ++nc;
                    }
                 }

              if (nc != ncoef)
                error ("Wrong number of coefficient");
              else
                {

                  dolfin::Vector b;
                  dolfin::assemble (b, L);

                  if (! boundarycondition_type_loaded)
                    {
                      boundarycondition::register_type ();
                      boundarycondition_type_loaded = true;
                      mlock ();
                    }
                  for (std::size_t i = 1; i < nargin; ++i)
                    {
                      if (args(i).type_id () == boundarycondition::static_type_id ())
                        {
                          const boundarycondition & bc 
                            = static_cast<const boundarycondition&> (args(i).get_rep ());

                          const std::vector<boost::shared_ptr <const dolfin::DirichletBC> > & pbc
                                = bc.get_bc ();
                              for (std::size_t j = 0; j < pbc.size (); ++j)
                                pbc[j]->apply(b);
                        }
                    }

                  dim_vector dims;
                  dims.resize (2);
                  dims(0) = b.size ();
                  dims(1) = 1;
                  Array<double> myb (dims);

                  for (std::size_t i = 0; i < b.size (); ++i)
                    myb(i) = b[i];

                  retval = myb;
                }
            }
        }
    }
  return retval;
}