/*
 Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
 Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>

 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 "function.h"
#include <stdexcept>
#include <vector>
#include "dolfin_compat.h"

DEFUN_DLD (feval, args, nargout, "-*- texinfo -*-\n\
@deftypefn {Function File} {[@var{fx}, @var{fy}, @var{fz}]} = \
feval (@var{function_name}, @var{x}, @var{y}, @var{z})\n\
Evaluate a function at a specific point of the domain and return the value. \n\
The input parameters are the function and the coordinates of the point where \n\
it has to be evaluated.\n\n\
Be aware that the number of input arguments and outputs depends on the \n\
dimensions of its domain and codomain.\n\
@example\n\
value = feval (p, x , y)\n\
@end example\n\
This is the expected call for a scalar field on a bidimensional domain.\n\
@seealso{Function}\n\
@end deftypefn")
{

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

  if (! function_type_loaded)
    {
      function::register_type ();
      function_type_loaded = true;
      mlock ();
    }

  if (args(0).type_id () == function::static_type_id ())
    {
      const function & fspo =
        static_cast<const function &> (args(0).get_rep ());

      if (!error_state)
        {
          const SHARED_PTR <const dolfin::Function>
            & f = fspo.get_pfun ();

          octave_idx_type pdim = f->geometric_dimension ();
          if (nargin != pdim + 1)
            { print_usage (); }
          else
            {
              std::vector <Matrix> coordinates;
              dim_vector dims;

              for (octave_idx_type in = 1; in <= pdim; ++in)
                {
                  if (! args(in).is_real_type ())
                    { error ("invalid argument"); }
                  else
                    {
                      Matrix aux = args(in).matrix_value ();
                      if (in == 1)
                        { dims = aux.dims (); }
                      else
                        {
                          dim_vector newdims = aux.dims ();
                          if (dims != newdims)
                            {
                              std::string msg = "all the input matrices should";
                              msg += " have the same size";
                              error (msg.c_str ());
                              break;
                            }
                        }
                      coordinates.push_back (aux);
                    }
                }

              if (! error_state)
                {
                  octave_idx_type vdim = f->value_dimension (0);
                  if (nargout != vdim)
                    { error ("wrong number of output arguments"); }
                  else
                    {
                      std::vector <Matrix> evaluations;
                      for (octave_idx_type out = 0; out < vdim; ++out)
                        { evaluations.push_back (Matrix (dims)); }

                      for (octave_idx_type k = 0; k < dims.numel (); ++k)
                        {
                          Array<double> point (dim_vector (pdim, 1));
                          for (octave_idx_type el = 0; el < pdim; ++el)
                            { point (el) = coordinates[el] (k); }
                          dolfin::Array<double>
                            x (point.length (), point.fortran_vec ());

                          Array<double> res (dim_vector (vdim, 1));
                          dolfin::Array<double>
                            values (res.length (), res.fortran_vec ());
                          try
                            {
                              f->eval (values, x);
                            }
                          catch (std::runtime_error & err)
                            {
                              std::string msg = "cannot evaluate a function";
                              msg += " outside of its domain";
                              error (msg.c_str ());
                              break;
                            }

                          for (octave_idx_type el = 0; el < vdim; ++el)
                            { evaluations[el] (k) = res (el); }
                        }

                      if (! error_state)
                        {
                          for (std::vector<Matrix>::iterator it =
                                 evaluations.begin ();
                               it != evaluations.end (); ++it)
                            { retval.append (octave_value (*it)); }
                        }
                    }
                }
            }
        }
    }

  return retval;
}