# HG changeset patch
# User Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
# Date 1403357287 -7200
# Node ID b07ec30369ab4faebec0ed9c11fd7b3365d2d13a
# Parent  f1c717e8a9713d5cfb6f8733ba429563d6c98d80
Add checks to feval

  * src/feval.cc: check the spatial dimensions of the coordinates and that
                  they belong to the domain

diff -r f1c717e8a971 -r b07ec30369ab src/feval.cc
--- a/src/feval.cc	Sun Jun 15 14:49:46 2014 +0200
+++ b/src/feval.cc	Sat Jun 21 15:28:07 2014 +0200
@@ -16,6 +16,8 @@
 */
 
 #include "function.h"
+#include <stdexcept>
+#include <sstream>
 
 DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\
 @deftypefn {Function File} {[@var{value}]} = \
@@ -57,30 +59,46 @@
               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);
+              if (point.rows () != f->geometric_dimension ())
+                {
+                  error ("feval: wrong coordinates dimension");
+                }
+              else
+                {
+                  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);
+                  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);
+                  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 ());
+                      try
+                        {
+                          f->eval (values, x);
+                        }
+                      catch (std::runtime_error & err)
+                        {
+                          std::string msg = "feval: cannot evaluate a function";
+                          msg += " outside of its domain";
+                          error (msg.c_str ());
+                        }
+                      for (uint j = 0; j < res_tmp.length (); ++j)
+                        res(i, j) = res_tmp (j);
+                    }
+                  retval = octave_value (res);
                 }
-              retval = octave_value (res);
             }
         }
     }