Mercurial > fem-fenics-eugenio

comparison src/feval.cc @ 216:a28b50969020

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
New version of feval which accept a matrix or an array as points where the evaluation is done.
author gedeone-octave <marcovass89@hotmail.it>
date Wed, 05 Mar 2014 11:59:02 +0000
parents ca7eb016cf64
children b07ec30369ab
comparison
equal deleted inserted replaced
215:cdbf0263863b 216:a28b50969020
1 /* 1 /*
2 Copyright (C) 2013-2014 Marco Vassallo <gedeone-octave@users.sourceforge.net> 2 Copyright (C) 2013 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
40 DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\ 20 DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\
41 @deftypefn {Function File} {[@var{value}]} = \ 21 @deftypefn {Function File} {[@var{value}]} = \
42 feval (@var{function_name}, @var{Coordinate})\n\ 22 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\
47 Evaluate a function at a specific point of the domain and return the value. \n\ 23 Evaluate a function at a specific point of the domain and return the value. \n\
48 With only two arguments, @var{Coordinate} is assumed to be an array holding \n\ 24 The input parameters are the function and the point where it has to\
49 the coordinates of the point where the function should be evaluated. With \n\ 25 be evaluated.\n\
50 three or more arguments, the coordinates are passed separately in @var{x}, \n\ 26 The point can be either a vector or a matrix. If it is a matrix, each column\
51 @var{y} and optionally @var{z}. In the latter case and if the function \n\ 27 represents a different point at which we evaluates 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\
54 @seealso{Function}\n\ 28 @seealso{Function}\n\
55 @end deftypefn") 29 @end deftypefn")
56 { 30 {
57 31
58 int nargin = args.length (); 32 int nargin = args.length ();
59 octave_value retval = 0; 33 octave_value retval=0;
60 34
61 if (nargin < 2) 35 if (nargin < 2 || nargin > 2)
62 print_usage (); 36 print_usage ();
63 else 37 else
64 { 38 {
65 if (! function_type_loaded) 39 if (! function_type_loaded)
66 { 40 {
71 45
72 if (args(0).type_id () == function::static_type_id ()) 46 if (args(0).type_id () == function::static_type_id ())
73 { 47 {
74 const function & fspo = 48 const function & fspo =
75 static_cast<const function&> (args(0).get_rep ()); 49 static_cast<const function&> (args(0).get_rep ());
76 const boost::shared_ptr<const dolfin::Function> 50 Matrix point= args(1).matrix_value ();
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);
83 51
84 if (nargin == 2) 52 if (!error_state)
85 { 53 {
86 Array<double> point = args(1).array_value (); 54 const boost::shared_ptr<const dolfin::Function>
87 if (!error_state) 55 & f = fspo.get_pfun ();
56
57 if (point.rows () == 1)
58 point = point.transpose ();
59
60 dim_vector dims;
61 dims.resize (2);
62 dims(0) = f->value_dimension (0);
63 dims(1) = point.cols ();
64 Matrix res (dims);
65
66 dim_vector dims_tmp;
67 dims_tmp.resize (2);
68 dims_tmp(0) = f->value_dimension (0);
69 dims_tmp(1) = 1;
70 Array<double> res_tmp (dims_tmp);
71
72 for (uint i = 0; i < point.cols (); ++i)
88 { 73 {
89 evaluate (*f, point, res); 74 Array<double> point_tmp = point.column (i);
90 retval = octave_value (res); 75 dolfin::Array<double>
76 x (point_tmp.length (), point_tmp.fortran_vec ());
77 dolfin::Array<double>
78 values (res_tmp.length (), res_tmp.fortran_vec ());
79 f->eval (values, x);
80 for (uint j = 0; j < res_tmp.length (); ++j)
81 res(i, j) = res_tmp (j);
91 } 82 }
92 } 83 retval = octave_value (res);
93 else
94 {
95 dim_vector inDims;
96 inDims.resize (2);
97 inDims(0) = 1;
98 inDims(1) = nargin - 1;
99 Array<double> point(inDims);
100
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 }
115
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 }
135 } 84 }
136 } 85 }
137 } 86 }
138
139 return retval; 87 return retval;
140 } 88 }