Mercurial > fem-fenics-eugenio
changeset 233:b07ec30369ab
Add checks to feval
* src/feval.cc: check the spatial dimensions of the coordinates and that
they belong to the domain
| author | Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> |
|---|---|
| date | Sat, 21 Jun 2014 15:28:07 +0200 |
| parents | f1c717e8a971 |
| children | 94895618701f |
| files | src/feval.cc |
| diffstat | 1 files changed, 39 insertions(+), 21 deletions(-) [+] |
line wrap: on
line diff
--- 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); } } }