/*
 Copyright (C) 2013 Marco Vassallo  <gedeone-octave@users.sourceforge.net> 

 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 3 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 "form.h"
#include "boundarycondition.h"

DEFUN_DLD (assemble_system, args, nargout, "-*- texinfo -*-\n\
@deftypefn {Function File} {[@var{A}], [@var{b}], [@var{x}(Optional)]} = \
assemble_system (@var{form a}, @var{form L}, @var{DirichletBC}(Optional), @var{...}) \n\
The input arguments are\n\
@itemize @bullet\n\
@item @var{form a} the bilinear form to assemble.\n\
@item @var{form a} the linear form to assemble.\n\
@item @var{DirichletBC} represents the optional BC that you wish to apply to\n\
the system. If more than one BC has to be applied, just list them.\n\
@end itemize \n\
The output @var{A} is a discretized representation of the system:\n\
@itemize @bullet\n\
@item @var{A} is the sparse Matrix corresponding to the @var{form a}\n\
@item @var{A} is the Vector corresponding to the @var{form L}\n\
@end itemize \n\
If you need to apply boundary condition to a system for a nonlinear problem \n\
then you should provide as 3rd argument the vector and you will receive it back\n\
as the third output argument. For an example of this situation, you can look\n\
the example HyperElasticity.m\n\
@seealso{BilinearForm, LinearForm, ResidualForm, JacobianForm}\n\
@end deftypefn")
{
  int nargin = args.length ();
  octave_value_list retval;
  if (! form_type_loaded)
    {
      form::register_type ();
      form_type_loaded = true;
      mlock ();
    }

  if (! boundarycondition_type_loaded)
    {
      boundarycondition::register_type ();
      boundarycondition_type_loaded = true;
      mlock ();
    }

  if (nargout == 2)
    {
      if (nargin < 2)
        print_usage ();
      else
        {

          if (args(0).type_id () == form::static_type_id ()
              && args(1).type_id () == form::static_type_id ())
            {
              const form & frm1 = static_cast<const form&> (args(0).get_rep ());
              const form & frm2 = static_cast<const form&> (args(1).get_rep ());

              if (! error_state)
                {
                  const dolfin::Form & a = frm1.get_form ();
                  const dolfin::Form & b = frm2.get_form ();
                  a.check ();
                  b.check ();

                  if (a.rank () == 2 && b.rank () == 1)
                    {
                      dolfin::Matrix A;
                      dolfin::assemble (A, a);
                      dolfin::Vector B;
                      dolfin::assemble (B, b);

                      for (std::size_t i = 2; 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(A, B);
                            }
                          else
                            error ("assemble_system: unknown argument type");
                        }

//                    It provides an estimation for excess of nnz elements
                      boost::tuples::tuple<const std::size_t*, const std::size_t*, const double*, int> aa = A.data ();
                      int nnz = aa.get<3> ();

                      std::cout << "Assembling Matrix from the bilinear form..."
                        << std::endl;
                      std::size_t nr = A.size (0), nc = A.size (1);
                      std::vector<double> data_tmp;
                      std::vector<std::size_t> cidx_tmp;

                      octave_idx_type dims = nnz, nz = 0, ii = 0;
                      ColumnVector ridx (dims), cidx (dims), data (dims);

                      for (std::size_t i = 0; i < nr; ++i)
                       {
                         A.getrow (i, cidx_tmp, data_tmp);
                         nz += cidx_tmp.size ();
                         if (dims < nz)
                           {
                             dims = 1.2 * ((nr * nz) / (i + 1));;
                             ridx.resize (dims);
                             cidx.resize (dims);
                             data.resize (dims);
                           }

                         for (octave_idx_type j = 0; j < cidx_tmp.size (); ++j)
                           {
                             ridx.xelem (ii + j) = i + 1;
                             cidx.xelem (ii + j) = cidx_tmp [j] + 1;
                             data.xelem (ii + j) = data_tmp [j];
                           }
                         ii = nz;
                       }

                      ridx.resize (ii);
                      cidx.resize (ii);
                      data.resize (ii);

                      SparseMatrix sm (data, ridx, cidx, nr, nc);
                      retval(0) = sm;

                      std::cout << "Assembling Vector from the linear form..."
                        << std::endl;
                      dim_vector dim;
                      dim.resize (2);
                      dim(0) = B.size ();
                      dim(1) = 1;
                      Array<double> myb (dim);

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

                      retval(1) = myb;
                    }
                }
              else
                error ("assemble_system: unknown size");
            }
        }
    }
  else if (nargout == 3)
    {
      std::cout << "Assemble_system: apply boundary condition to a vector for a nonlinear problem..."
        << std::endl;
      if (nargin < 3)
        print_usage ();
      else
        {
          if (args(0).type_id () == form::static_type_id ()
              && args(1).type_id () == form::static_type_id ())
            {
              const form & frm1 = static_cast<const form&> (args(0).get_rep ());
              const form & frm2 = static_cast<const form&> (args(1).get_rep ());
              const Array<double> myx = args(2).array_value ();

              if (! error_state)
                {
                  const dolfin::Form & a = frm1.get_form ();
                  const dolfin::Form & b = frm2.get_form ();
                  a.check ();
                  b.check ();

                  if (a.rank () == 2 && b.rank () == 1)
                    {
                      dolfin::Matrix A;
                      dolfin::assemble (A, a);
                      dolfin::Vector B;
                      dolfin::assemble (B, b);
                      dolfin::Vector x (myx.length ());

                      for (std::size_t i = 0; i < myx.length (); ++i)
                        x.setitem (i, myx.xelem (i));

                      for (std::size_t i = 3; 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(A, B, x);
                            }
                          else
                            error ("assemble_system: unknown argument type");
                        }

//                    It provides an estimation for excess of nnz elements
                      boost::tuples::tuple<const std::size_t*, const std::size_t*, const double*, int> aa = A.data ();
                      int nnz = aa.get<3> ();

                      std::cout << "Assembling Matrix from the bilinear form..."
                        << std::endl;
                      std::size_t nr = A.size (0), nc = A.size (1);
                      std::vector<double> data_tmp;
                      std::vector<std::size_t> cidx_tmp;

                      octave_idx_type dims = nnz, nz = 0, ii = 0;
                      ColumnVector ridx (dims), cidx (dims), data (dims);

                      for (std::size_t i = 0; i < nr; ++i)
                       {
                         A.getrow (i, cidx_tmp, data_tmp);
                         nz += cidx_tmp.size ();
                         if (dims < nz)
                           {
                             dims = 1.2 * ((nr * nz) / (i + 1));;
                             ridx.resize (dims);
                             cidx.resize (dims);
                             data.resize (dims);
                           }

                         for (octave_idx_type j = 0; j < cidx_tmp.size (); ++j)
                           {
                             ridx.xelem (ii + j) = i + 1;
                             cidx.xelem (ii + j) = cidx_tmp [j] + 1;
                             data.xelem (ii + j) = data_tmp [j];
                           }

                         ii = nz;
                       }

                      ridx.resize (ii);
                      cidx.resize (ii);
                      data.resize (ii);

                      SparseMatrix sm (data, ridx, cidx, nr, nc);
                      retval(0) = sm;

                      std::cout << "Assembling Vector from the linear form..."
                        << std::endl;
                      dim_vector dim;
                      dim.resize (2);
                      dim(0) = B.size ();
                      dim(1) = 1;
                      Array<double> myb (dim), myc (dim);

                      for (std::size_t i = 0; i < B.size (); ++i)
                        {
                          myb.xelem (i) = B[i];
                          myc.xelem (i) = x[i];
                        }

                      retval(1) = myb;
                      retval(2) = myc;
                    }
                }
              else
                error ("assemble_system: unknown size");
            }
        }
    }
  return retval;
}