Mercurial > fem-fenics-eugenio

comparison src/feval.cc @ 233:b07ec30369ab

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 a28b50969020
children a51c09492e30
comparison
equal deleted inserted replaced
232:f1c717e8a971 233:b07ec30369ab
14 You should have received a copy of the GNU General Public License along with 14 You should have received a copy of the GNU General Public License along with
15 this program; if not, see <http://www.gnu.org/licenses/>. 15 this program; if not, see <http://www.gnu.org/licenses/>.
16 */ 16 */
17 17
18 #include "function.h" 18 #include "function.h"
19 #include <stdexcept>
20 #include <sstream>
19 21
20 DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\ 22 DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\
21 @deftypefn {Function File} {[@var{value}]} = \ 23 @deftypefn {Function File} {[@var{value}]} = \
22 feval (@var{function_name}, @var{Coordinate})\n\ 24 feval (@var{function_name}, @var{Coordinate})\n\
23 Evaluate a function at a specific point of the domain and return the value. \n\ 25 Evaluate a function at a specific point of the domain and return the value. \n\
55 & f = fspo.get_pfun (); 57 & f = fspo.get_pfun ();
56 58
57 if (point.rows () == 1) 59 if (point.rows () == 1)
58 point = point.transpose (); 60 point = point.transpose ();
59 61
60 dim_vector dims; 62 if (point.rows () != f->geometric_dimension ())
61 dims.resize (2); 63 {
62 dims(0) = f->value_dimension (0); 64 error ("feval: wrong coordinates dimension");
63 dims(1) = point.cols (); 65 }
64 Matrix res (dims); 66 else
67 {
68 dim_vector dims;
69 dims.resize (2);
70 dims(0) = f->value_dimension (0);
71 dims(1) = point.cols ();
72 Matrix res (dims);
65 73
66 dim_vector dims_tmp; 74 dim_vector dims_tmp;
67 dims_tmp.resize (2); 75 dims_tmp.resize (2);
68 dims_tmp(0) = f->value_dimension (0); 76 dims_tmp(0) = f->value_dimension (0);
69 dims_tmp(1) = 1; 77 dims_tmp(1) = 1;
70 Array<double> res_tmp (dims_tmp); 78 Array<double> res_tmp (dims_tmp);
71 79
72 for (uint i = 0; i < point.cols (); ++i) 80 for (uint i = 0; i < point.cols (); ++i)
73 { 81 {
74 Array<double> point_tmp = point.column (i); 82 Array<double> point_tmp = point.column (i);
75 dolfin::Array<double> 83 dolfin::Array<double>
76 x (point_tmp.length (), point_tmp.fortran_vec ()); 84 x (point_tmp.length (), point_tmp.fortran_vec ());
77 dolfin::Array<double> 85 dolfin::Array<double>
78 values (res_tmp.length (), res_tmp.fortran_vec ()); 86 values (res_tmp.length (), res_tmp.fortran_vec ());
79 f->eval (values, x); 87 try
80 for (uint j = 0; j < res_tmp.length (); ++j) 88 {
81 res(i, j) = res_tmp (j); 89 f->eval (values, x);
90 }
91 catch (std::runtime_error & err)
92 {
93 std::string msg = "feval: cannot evaluate a function";
94 msg += " outside of its domain";
95 error (msg.c_str ());
96 }
97 for (uint j = 0; j < res_tmp.length (); ++j)
98 res(i, j) = res_tmp (j);
99 }
100 retval = octave_value (res);
82 } 101 }
83 retval = octave_value (res);
84 } 102 }
85 } 103 }
86 } 104 }
87 return retval; 105 return retval;
88 } 106 }