# HG changeset patch
# User Daniel Kraft <address@hidden>
# Date 1393931450 -3600
# Node ID 76829daacb1251f674922f572ca99ea1e69db5cf
# Parent  eaba4d4eaf1536bb5a85c6140855a67f0dacbe60
Extend feval function to support also passing the coordinates separately
as multiple arguments and evaluating the function at multiple
points at once.

diff -r eaba4d4eaf15 -r 76829daacb12 src/feval.cc
--- a/src/feval.cc	Tue Mar 04 12:10:11 2014 +0100
+++ b/src/feval.cc	Tue Mar 04 12:10:50 2014 +0100
@@ -1,5 +1,5 @@
 /*
- Copyright (C) 2013 Marco Vassallo <address@hidden>
+ 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
@@ -17,20 +17,48 @@
 
 #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\
address@hidden {Function File} address@hidden = \
+feval (@var{function_name}, @var{x}, @var{y})\n\
address@hidden {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\
-The input parameters are the function and the point where it has to\
-be evaluated.\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\
address@hidden 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;
+  octave_value retval = 0;
   
-  if (nargin < 2 || nargin > 2)
+  if (nargin < 2)
     print_usage ();
   else
     {
@@ -45,25 +73,68 @@
         {
           const function & fspo =
             static_cast<const function&> (args(0).get_rep ());
-          Array<double> point= args(1).array_value ();
+          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 (!error_state)
+          if (nargin == 2)
             {
-              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);
-              dolfin::Array<double> x(point.length (), point.fortran_vec ());
-              dolfin::Array<double> values(res.length (), res.fortran_vec ());
+              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);
 
-              f->eval (values, x);
+              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;
+                }
 
-              retval = octave_value (res);
+              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;
 }