Mercurial > fem-fenics-eugenio

comparison src/feval.cc @ 237:418a5119047b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
New interface for function evaluation
author Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
date Mon, 23 Jun 2014 20:28:28 +0200
parents a51c09492e30
children 0f14cdbcaed3
comparison
equal deleted inserted replaced
236:a51c09492e30 237:418a5119047b
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> 19 #include <stdexcept>
20 #include <vector>
20 21
21 DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\ 22 DEFUN_DLD (feval, args, nargout, "-*- texinfo -*-\n\
22 @deftypefn {Function File} {[@var{value}]} = \ 23 @deftypefn {Function File} {[@var{fx}, @var{fy}, @var{fz}]} = \
23 feval (@var{function_name}, @var{Coordinate})\n\ 24 feval (@var{function_name}, @var{x}, @var{y}, @var{z})\n\
24 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\
25 The input parameters are the function and the point where it has to\ 26 The input parameters are the function and the coordinates of the point where \n\
26 be evaluated.\n\ 27 it has to be evaluated.\n\n\
27 The point can be either a vector or a matrix. If it is a matrix, each column\ 28 Be aware that the number of input arguments and outputs depends on the \n\
28 represents a different point at which we evaluates the function.\n\ 29 dimensions of its domain and codomain.\n\
30 @example\n\
31 value = feval (p, x , y)\n\
32 @end example\n\
33 This is the expected call for a scalar field on a bidimensional domain.\n\
29 @seealso{Function}\n\ 34 @seealso{Function}\n\
30 @end deftypefn") 35 @end deftypefn")
31 { 36 {
32 37
33 int nargin = args.length (); 38 int nargin = args.length ();
34 octave_value retval=0; 39 octave_value_list retval;
35 40
36 if (nargin < 2 || nargin > 2) 41 if (! function_type_loaded)
37 print_usage ();
38 else
39 { 42 {
40 if (! function_type_loaded) 43 function::register_type ();
44 function_type_loaded = true;
45 mlock ();
46 }
47
48 if (args(0).type_id () == function::static_type_id ())
49 {
50 const function & fspo =
51 static_cast<const function&> (args(0).get_rep ());
52
53 if (!error_state)
41 { 54 {
42 function::register_type (); 55 const boost::shared_ptr<const dolfin::Function>
43 function_type_loaded = true; 56 & f = fspo.get_pfun ();
44 mlock ();
45 }
46 57
47 if (args(0).type_id () == function::static_type_id ()) 58 octave_idx_type pdim = f->geometric_dimension ();
48 { 59 if (nargin != pdim + 1)
49 const function & fspo = 60 print_usage ();
50 static_cast<const function&> (args(0).get_rep ()); 61 else
51 Matrix point= args(1).matrix_value (); 62 {
63 std::vector <Matrix> coordinates;
64 dim_vector dims;
52 65
53 if (!error_state) 66 for (octave_idx_type in = 1; in <= pdim; ++in)
54 { 67 {
55 const boost::shared_ptr<const dolfin::Function> 68 if (! args(in).is_real_type ())
56 & f = fspo.get_pfun (); 69 error ("invalid argument");
70 else
71 {
72 Matrix aux = args(in).matrix_value ();
73 if (in == 1)
74 dims = aux.dims ();
75 else
76 {
77 dim_vector newdims = aux.dims ();
78 if (dims != newdims)
79 {
80 std::string msg = "all the input matrices should";
81 msg += " have the same size";
82 error (msg.c_str ());
83 }
84 }
85 coordinates.push_back (aux);
86 }
87 }
57 88
58 if (point.rows () == 1) 89 octave_idx_type vdim = f->value_dimension (0);
59 point = point.transpose (); 90 if (nargout != vdim)
91 error ("wrong number of output arguments");
60 92
61 if (point.rows () != f->geometric_dimension ()) 93 std::vector <Matrix> evaluations;
94 for (octave_idx_type out = 0; out < vdim; ++out)
95 evaluations.push_back (Matrix (dims));
96
97 for (octave_idx_type k = 0; k < dims.numel (); ++k)
62 { 98 {
63 error ("feval: wrong coordinates dimension"); 99 Array<double> point (dim_vector (pdim, 1));
100 for (octave_idx_type el = 0; el < pdim; ++el)
101 point (el) = coordinates[el] (k);
102 dolfin::Array<double>
103 x (point.length (), point.fortran_vec ());
104
105 Array<double> res (dim_vector (vdim, 1));
106 dolfin::Array<double>
107 values (res.length (), res.fortran_vec ());
108 try
109 {
110 f->eval (values, x);
111 }
112 catch (std::runtime_error & err)
113 {
114 std::string msg = "cannot evaluate a function";
115 msg += " outside of its domain";
116 error (msg.c_str ());
117 }
118
119 for (octave_idx_type el = 0; el < vdim; ++el)
120 evaluations[el] (k) = res (el);
64 } 121 }
65 else
66 {
67 dim_vector dims;
68 dims.resize (2);
69 dims(0) = f->value_dimension (0);
70 dims(1) = point.cols ();
71 Matrix res (dims);
72 122
73 dim_vector dims_tmp; 123 for (std::vector<Matrix>::iterator it = evaluations.begin ();
74 dims_tmp.resize (2); 124 it != evaluations.end (); ++it)
75 dims_tmp(0) = f->value_dimension (0); 125 retval.append (octave_value (*it));
76 dims_tmp(1) = 1;
77 Array<double> res_tmp (dims_tmp);
78
79 for (uint i = 0; i < point.cols (); ++i)
80 {
81 Array<double> point_tmp = point.column (i);
82 dolfin::Array<double>
83 x (point_tmp.length (), point_tmp.fortran_vec ());
84 dolfin::Array<double>
85 values (res_tmp.length (), res_tmp.fortran_vec ());
86 try
87 {
88 f->eval (values, x);
89 }
90 catch (std::runtime_error & err)
91 {
92 std::string msg = "feval: cannot evaluate a function";
93 msg += " outside of its domain";
94 error (msg.c_str ());
95 }
96 for (uint j = 0; j < res_tmp.length (); ++j)
97 res(i, j) = res_tmp (j);
98 }
99 retval = octave_value (res);
100 }
101 } 126 }
102 } 127 }
103 } 128 }
129
104 return retval; 130 return retval;
105 } 131 }