/*
 Copyright (C) 2013-2014 Marco Vassallo <address@hidden>

 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"

/**
 * Internal helper routine to evaluate the given function at a single
 * point (as Array<double>).  Returned is the function value also as
 * Array<double>, which gets stored into the given array.
 * too, so that they don't have to be constructed.
 * @param f Dolfin function handle to evaluate.
 * @param pt Point where we want to evaluate f as Octave array of coordinates.
 * @param res Store the value here, should already be of the correct size.
 */
static void
evaluate (const dolfin::Function& f, const Array<double>& pt,
          Array<double>& res)
{
  Array<double>& ptNonconst = const_cast<Array<double>&> (pt);
  dolfin::Array<double> x(ptNonconst.length (), ptNonconst.fortran_vec ());
  dolfin::Array<double> values(res.length (), res.fortran_vec ());

  f.eval (values, x);
}

DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\
@deftypefn {Function File} address@hidden = \
feval (@var{function_name}, @var{Coordinate})\n\
@deftypefnx {Function File} address@hidden = \
feval (@var{function_name}, @var{x}, @var{y})\n\
@deftypefnx {Function File} address@hidden = \
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\
With only two arguments, @var{Coordinate} is assumed to be an array holding \n\
the coordinates of the point where the function should be evaluated.  With \n\
three or more arguments, the coordinates are passed separately in @var{x}, \n\
@var{y} and optionally @var{z}.  In the latter case and if the function \n\
returns scalar values, @var{x}, @var{y} and @var{z} may also be arrays to \n\
evaluate the function at multiple points at once.\n\
@seealso{Function}\n\
@end deftypefn")
{

  int nargin = args.length ();
  octave_value retval = 0;
  
  if (nargin < 2)
    print_usage ();
  else
    {
      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 ());
          const boost::shared_ptr<const dolfin::Function> 
            & f = fspo.get_pfun ();
          dim_vector dims;
          dims.resize (2);
          dims(0) = 1;
          dims(1) = f->value_dimension (0);
          Array<double> res(dims);

          if (nargin == 2)
            {
              Array<double> point = args(1).array_value ();
              if (!error_state)
                {
                  evaluate (*f, point, res);
                  retval = octave_value (res);
                }
            }
          else
            {
              dim_vector inDims;
              inDims.resize (2);
              inDims(0) = 1;
              inDims(1) = nargin - 1;
              Array<double> point(inDims);

              std::vector<Array<double> > coords;
              coords.reserve (inDims(1));
              dim_vector argDims;
              for (unsigned i = 1; i < nargin; ++i)
                {
                  coords.push_back (args(i).array_value ());
                  const dim_vector curDims = coords.back ().dims ();
                  if (i > 1 && argDims != curDims)
                    {
                      error ("feval: coordinate dimensions mismatch");
                      break;
                    }
                  argDims = curDims;
                }

              if (res.nelem () != 1)
                error ("feval: separate coordinate version only supported"
                       "for scalar functions");

              if (!error_state)
                {
                  Array<double> retValArray(argDims);

                  for (size_t i = 0; i < retValArray.nelem (); ++i)
                    {
                      for (unsigned j = 0; j < coords.size (); ++j)
                        point(j) = coords[j].fortran_vec ()[i];

                      evaluate (*f, point, res);
                      retValArray.fortran_vec ()[i] = res(0);
                    }

                  retval = octave_value (retValArray);
                }
            }
        }
    }

  return retval;
}