Mercurial > fem-fenics-eugenio

annotate 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
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
1 /*
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
2 Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
3
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
4 This program is free software; you can redistribute it and/or modify it under
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
5 the terms of the GNU General Public License as published by the Free Software
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
6 Foundation; either version 3 of the License, or (at your option) any later
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
7 version.
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
8
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
9 This program is distributed in the hope that it will be useful, but WITHOUT
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
10 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
11 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
12 details.
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
13
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
14 You should have received a copy of the GNU General Public License along with
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
15 this program; if not, see <http://www.gnu.org/licenses/>.
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
16 */
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
17
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
18 #include "function.h"
233
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
19 #include <stdexcept>
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
20 #include <vector>
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
21
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
22 DEFUN_DLD (feval, args, nargout, "-*- texinfo -*-\n\
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
23 @deftypefn {Function File} {[@var{fx}, @var{fy}, @var{fz}]} = \
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
24 feval (@var{function_name}, @var{x}, @var{y}, @var{z})\n\
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
25 Evaluate a function at a specific point of the domain and return the value. \n\
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
26 The input parameters are the function and the coordinates of the point where \n\
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
27 it has to be evaluated.\n\n\
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
28 Be aware that the number of input arguments and outputs depends on the \n\
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
29 dimensions of its domain and codomain.\n\
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
30 @example\n\
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
31 value = feval (p, x , y)\n\
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
32 @end example\n\
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
33 This is the expected call for a scalar field on a bidimensional domain.\n\
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
34 @seealso{Function}\n\
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
35 @end deftypefn")
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
36 {
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
37
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
38 int nargin = args.length ();
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
39 octave_value_list retval;
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
40
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
41 if (! function_type_loaded)
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
42 {
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
43 function::register_type ();
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
44 function_type_loaded = true;
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
45 mlock ();
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
46 }
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
47
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
48 if (args(0).type_id () == function::static_type_id ())
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
49 {
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
50 const function & fspo =
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
51 static_cast<const function&> (args(0).get_rep ());
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
52
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
53 if (!error_state)
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
54 {
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
55 const boost::shared_ptr<const dolfin::Function>
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
56 & f = fspo.get_pfun ();
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
57
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
58 octave_idx_type pdim = f->geometric_dimension ();
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
59 if (nargin != pdim + 1)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
60 print_usage ();
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
61 else
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
62 {
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
63 std::vector <Matrix> coordinates;
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
64 dim_vector dims;
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
65
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
66 for (octave_idx_type in = 1; in <= pdim; ++in)
233
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
67 {
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
68 if (! args(in).is_real_type ())
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
69 error ("invalid argument");
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
70 else
233
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
71 {
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
72 Matrix aux = args(in).matrix_value ();
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
73 if (in == 1)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
74 dims = aux.dims ();
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
75 else
233
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
76 {
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
77 dim_vector newdims = aux.dims ();
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
78 if (dims != newdims)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
79 {
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
80 std::string msg = "all the input matrices should";
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
81 msg += " have the same size";
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
82 error (msg.c_str ());
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
83 }
233
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
84 }
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
85 coordinates.push_back (aux);
233
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
86 }
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
87 }
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
88
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
89 octave_idx_type vdim = f->value_dimension (0);
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
90 if (nargout != vdim)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
91 error ("wrong number of output arguments");
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
92
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
93 std::vector <Matrix> evaluations;
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
94 for (octave_idx_type out = 0; out < vdim; ++out)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
95 evaluations.push_back (Matrix (dims));
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
96
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
97 for (octave_idx_type k = 0; k < dims.numel (); ++k)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
98 {
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
99 Array<double> point (dim_vector (pdim, 1));
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
100 for (octave_idx_type el = 0; el < pdim; ++el)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
101 point (el) = coordinates[el] (k);
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
102 dolfin::Array<double>
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
103 x (point.length (), point.fortran_vec ());
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
104
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
105 Array<double> res (dim_vector (vdim, 1));
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
106 dolfin::Array<double>
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
107 values (res.length (), res.fortran_vec ());
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
108 try
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
109 {
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
110 f->eval (values, x);
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
111 }
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
112 catch (std::runtime_error & err)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
113 {
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
114 std::string msg = "cannot evaluate a function";
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
115 msg += " outside of its domain";
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
116 error (msg.c_str ());
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
117 }
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
118
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
119 for (octave_idx_type el = 0; el < vdim; ++el)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
120 evaluations[el] (k) = res (el);
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
121 }
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
122
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
123 for (std::vector<Matrix>::iterator it = evaluations.begin ();
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
124 it != evaluations.end (); ++it)
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
125 retval.append (octave_value (*it));
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
126 }
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
127 }
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
128 }
237
418a5119047b New interface for function evaluation
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 236
diff changeset
129
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
130 return retval;
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
131 }