Mercurial > fem-fenics-eugenio

annotate 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
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>
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
20 #include <sstream>
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
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
22 DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
23 @deftypefn {Function File} {[@var{value}]} = \
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
24 feval (@var{function_name}, @var{Coordinate})\n\
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\
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
26 The input parameters are the function and the point where it has to\
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
27 be evaluated.\n\
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
28 The point can be either a vector or a matrix. If it is a matrix, each column\
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
29 represents a different point at which we evaluates the 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
30 @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
31 @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
32 {
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
33
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
34 int nargin = args.length ();
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
35 octave_value retval=0;
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 if (nargin < 2 || nargin > 2)
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
38 print_usage ();
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
39 else
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
40 {
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
41 if (! function_type_loaded)
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
42 {
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
43 function::register_type ();
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
44 function_type_loaded = true;
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
45 mlock ();
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
46 }
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
47
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
48 if (args(0).type_id () == function::static_type_id ())
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
49 {
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
50 const function & fspo =
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
51 static_cast<const function&> (args(0).get_rep ());
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
52 Matrix point= args(1).matrix_value ();
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
53
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
54 if (!error_state)
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
55 {
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
56 const boost::shared_ptr<const dolfin::Function>
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
57 & f = fspo.get_pfun ();
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
58
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
59 if (point.rows () == 1)
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
60 point = point.transpose ();
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
61
233
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
62 if (point.rows () != f->geometric_dimension ())
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
63 {
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
64 error ("feval: wrong coordinates dimension");
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
65 }
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
66 else
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
67 {
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
68 dim_vector dims;
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
69 dims.resize (2);
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
70 dims(0) = f->value_dimension (0);
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
71 dims(1) = point.cols ();
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
72 Matrix res (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
73
233
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
74 dim_vector dims_tmp;
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
75 dims_tmp.resize (2);
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
76 dims_tmp(0) = f->value_dimension (0);
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
77 dims_tmp(1) = 1;
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
78 Array<double> res_tmp (dims_tmp);
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
79
233
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
80 for (uint i = 0; i < point.cols (); ++i)
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
81 {
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
82 Array<double> point_tmp = point.column (i);
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
83 dolfin::Array<double>
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
84 x (point_tmp.length (), point_tmp.fortran_vec ());
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
85 dolfin::Array<double>
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
86 values (res_tmp.length (), res_tmp.fortran_vec ());
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
87 try
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
88 {
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
89 f->eval (values, x);
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
90 }
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
91 catch (std::runtime_error & err)
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
92 {
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
93 std::string msg = "feval: cannot evaluate a function";
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
94 msg += " outside of its domain";
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
95 error (msg.c_str ());
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
96 }
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
97 for (uint j = 0; j < res_tmp.length (); ++j)
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
98 res(i, j) = res_tmp (j);
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
99 }
b07ec30369ab Add checks to feval
Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
parents: 216
diff changeset
100 retval = octave_value (res);
216
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
101 }
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
102 }
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
103 }
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
104 }
a28b50969020 New version of feval which accept a matrix or an array as points
gedeone-octave <marcovass89@hotmail.it>
parents: 212
diff changeset
105 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
106 }