# HG changeset patch
# User gedeone-octave <marcovass89@hotmail.it>
# Date 1394020742 0
# Node ID a28b50969020c7b40b39df37ee6584b0fa619a7a
# Parent  cdbf0263863b347e85aa40f897eac87f90fc38ae
New version of feval which accept a matrix or an array as points
where the evaluation is done.

diff -r cdbf0263863b -r a28b50969020 src/feval.cc
--- a/src/feval.cc	Wed Mar 05 10:32:56 2014 +0000
+++ b/src/feval.cc	Wed Mar 05 11:59:02 2014 +0000
@@ -1,140 +1,88 @@
-/*
- Copyright (C) 2013-2014 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 "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} {[@var{value}]} = \
-feval (@var{function_name}, @var{Coordinate})\n\
-@deftypefnx {Function File} {[@var{value}]} = \
-feval (@var{function_name}, @var{x}, @var{y})\n\
-@deftypefnx {Function File} {[@var{value}]} = \
-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;
-}
+/*
+ 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 "function.h"
+
+DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\
+@deftypefn {Function File} {[@var{value}]} = \
+feval (@var{function_name}, @var{Coordinate})\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\
+The point can be either a vector or a matrix. If it is a matrix, each column\
+represents a different point at which we evaluates the function.\n\
+@seealso{Function}\n\
+@end deftypefn")
+{
+
+  int nargin = args.length ();
+  octave_value retval=0;
+  
+  if (nargin < 2 || 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 ());
+          Matrix point= args(1).matrix_value ();
+
+          if (!error_state)
+            {
+              const boost::shared_ptr<const dolfin::Function> 
+                & f = fspo.get_pfun ();
+
+              if (point.rows () == 1)
+                point = point.transpose ();
+
+              dim_vector dims;
+              dims.resize (2);
+              dims(0) = f->value_dimension (0);
+              dims(1) = point.cols ();
+              Matrix res (dims);
+
+              dim_vector dims_tmp;
+              dims_tmp.resize (2);
+              dims_tmp(0) = f->value_dimension (0);
+              dims_tmp(1) = 1;
+              Array<double> res_tmp (dims_tmp);
+
+              for (uint i = 0; i < point.cols (); ++i)
+                {
+                  Array<double> point_tmp = point.column (i);
+                  dolfin::Array<double> 
+                    x (point_tmp.length (), point_tmp.fortran_vec ());
+                  dolfin::Array<double> 
+                    values (res_tmp.length (), res_tmp.fortran_vec ());
+                  f->eval (values, x);
+                  for (uint j = 0; j < res_tmp.length (); ++j)
+                    res(i, j) = res_tmp (j);
+                }
+              retval = octave_value (res);
+            }
+        }
+    }
+  return retval;
+}