Mercurial > fem-fenics-eugenio

comparison src/feval.cc @ 212:ca7eb016cf64

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Extend feval function to support also passing the coordinates separately as multiple arguments and evaluating the function at multiple points at once. From Daniel Kraft.
author gedeone-octave <marcovass89@hotmail.it>
date Tue, 04 Mar 2014 12:09:26 +0000
parents 66071811eef8
children a28b50969020
comparison
equal deleted inserted replaced
209:37bcbf7095be 212:ca7eb016cf64
1 /* 1 /*
2 Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net> 2 Copyright (C) 2013-2014 Marco Vassallo <gedeone-octave@users.sourceforge.net>
3 3
4 This program is free software; you can redistribute it and/or modify it under 4 This program is free software; you can redistribute it and/or modify it under
5 the terms of the GNU General Public License as published by the Free Software 5 the terms of the GNU General Public License as published by the Free Software
6 Foundation; either version 3 of the License, or (at your option) any later 6 Foundation; either version 3 of the License, or (at your option) any later
7 version. 7 version.
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 19
20 /**
21 * Internal helper routine to evaluate the given function at a single
22 * point (as Array<double>). Returned is the function value also as
23 * Array<double>, which gets stored into the given array.
24 * too, so that they don't have to be constructed.
25 * @param f Dolfin function handle to evaluate.
26 * @param pt Point where we want to evaluate f as Octave array of coordinates.
27 * @param res Store the value here, should already be of the correct size.
28 */
29 static void
30 evaluate (const dolfin::Function& f, const Array<double>& pt,
31 Array<double>& res)
32 {
33 Array<double>& ptNonconst = const_cast<Array<double>&> (pt);
34 dolfin::Array<double> x(ptNonconst.length (), ptNonconst.fortran_vec ());
35 dolfin::Array<double> values(res.length (), res.fortran_vec ());
36
37 f.eval (values, x);
38 }
39
20 DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\ 40 DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\
21 @deftypefn {Function File} {[@var{value}]} = \ 41 @deftypefn {Function File} {[@var{value}]} = \
22 feval (@var{function_name}, @var{Coordinate})\n\ 42 feval (@var{function_name}, @var{Coordinate})\n\
43 @deftypefnx {Function File} {[@var{value}]} = \
44 feval (@var{function_name}, @var{x}, @var{y})\n\
45 @deftypefnx {Function File} {[@var{value}]} = \
46 feval (@var{function_name}, @var{x}, @var{y}, @var{z})\n\
23 Evaluate a function at a specific point of the domain and return the value. \n\ 47 Evaluate a function at a specific point of the domain and return the value. \n\
24 The input parameters are the function and the point where it has to\ 48 With only two arguments, @var{Coordinate} is assumed to be an array holding \n\
25 be evaluated.\n\ 49 the coordinates of the point where the function should be evaluated. With \n\
50 three or more arguments, the coordinates are passed separately in @var{x}, \n\
51 @var{y} and optionally @var{z}. In the latter case and if the function \n\
52 returns scalar values, @var{x}, @var{y} and @var{z} may also be arrays to \n\
53 evaluate the function at multiple points at once.\n\
26 @seealso{Function}\n\ 54 @seealso{Function}\n\
27 @end deftypefn") 55 @end deftypefn")
28 { 56 {
29 57
30 int nargin = args.length (); 58 int nargin = args.length ();
31 octave_value retval=0; 59 octave_value retval = 0;
32 60
33 if (nargin < 2 || nargin > 2) 61 if (nargin < 2)
34 print_usage (); 62 print_usage ();
35 else 63 else
36 { 64 {
37 if (! function_type_loaded) 65 if (! function_type_loaded)
38 { 66 {
43 71
44 if (args(0).type_id () == function::static_type_id ()) 72 if (args(0).type_id () == function::static_type_id ())
45 { 73 {
46 const function & fspo = 74 const function & fspo =
47 static_cast<const function&> (args(0).get_rep ()); 75 static_cast<const function&> (args(0).get_rep ());
48 Array<double> point= args(1).array_value (); 76 const boost::shared_ptr<const dolfin::Function>
77 & f = fspo.get_pfun ();
78 dim_vector dims;
79 dims.resize (2);
80 dims(0) = 1;
81 dims(1) = f->value_dimension (0);
82 Array<double> res(dims);
49 83
50 if (!error_state) 84 if (nargin == 2)
51 { 85 {
52 const boost::shared_ptr<const dolfin::Function> 86 Array<double> point = args(1).array_value ();
53 & f = fspo.get_pfun (); 87 if (!error_state)
54 dim_vector dims; 88 {
55 dims.resize (2); 89 evaluate (*f, point, res);
56 dims(0) = 1; 90 retval = octave_value (res);
57 dims(1) = f->value_dimension (0); 91 }
58 Array<double> res (dims); 92 }
59 dolfin::Array<double> x(point.length (), point.fortran_vec ()); 93 else
60 dolfin::Array<double> values(res.length (), res.fortran_vec ()); 94 {
95 dim_vector inDims;
96 inDims.resize (2);
97 inDims(0) = 1;
98 inDims(1) = nargin - 1;
99 Array<double> point(inDims);
61 100
62 f->eval (values, x); 101 std::vector<Array<double> > coords;
102 coords.reserve (inDims(1));
103 dim_vector argDims;
104 for (unsigned i = 1; i < nargin; ++i)
105 {
106 coords.push_back (args(i).array_value ());
107 const dim_vector curDims = coords.back ().dims ();
108 if (i > 1 && argDims != curDims)
109 {
110 error ("feval: coordinate dimensions mismatch");
111 break;
112 }
113 argDims = curDims;
114 }
63 115
64 retval = octave_value (res); 116 if (res.nelem () != 1)
117 error ("feval: separate coordinate version only supported"
118 "for scalar functions");
119
120 if (!error_state)
121 {
122 Array<double> retValArray(argDims);
123
124 for (size_t i = 0; i < retValArray.nelem (); ++i)
125 {
126 for (unsigned j = 0; j < coords.size (); ++j)
127 point(j) = coords[j].fortran_vec ()[i];
128
129 evaluate (*f, point, res);
130 retValArray.fortran_vec ()[i] = res(0);
131 }
132
133 retval = octave_value (retValArray);
134 }
65 } 135 }
66 } 136 }
67 } 137 }
138
68 return retval; 139 return retval;
69 } 140 }