Mercurial > fem-fenics-eugenio
view src/feval.cc @ 261:f22588ae37af
Improve template meshfunction implementation
* inst/import_meshfunction_type.m: provide to the auxiliary functions a type name and
a valid identifier
* inst/private/generate_mf_*.m: add space before closing angle bracket to avoid parse
errors
author | Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> |
---|---|
date | Thu, 07 Aug 2014 11:13:54 +0200 |
parents | 5e9b5bbdc56b |
children | 61830a4f9ab9 |
line wrap: on
line source
/* Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net> Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, see <http://www.gnu.org/licenses/>. */ #include "function.h" #include <stdexcept> #include <vector> #include "dolfin_compat.h" DEFUN_DLD (feval, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Function File} {[@var{fx}, @var{fy}, @var{fz}]} = \ feval (@var{function_name}, @var{x}, @var{y}, @var{z})\n\ Evaluate a function at a specific point of the domain and return the value. \n\ The input parameters are the function and the coordinates of the point where \n\ it has to be evaluated.\n\n\ Be aware that the number of input arguments and outputs depends on the \n\ dimensions of its domain and codomain.\n\ @example\n\ value = feval (p, x , y)\n\ @end example\n\ This is the expected call for a scalar field on a bidimensional domain.\n\ @seealso{Function}\n\ @end deftypefn") { int nargin = args.length (); octave_value_list retval; if (! function_type_loaded) { function::register_type (); function_type_loaded = true; mlock (); } if (args(0).type_id () == function::static_type_id ()) { const function & fspo = static_cast<const function&> (args(0).get_rep ()); if (!error_state) { const SHARED_PTR <const dolfin::Function> & f = fspo.get_pfun (); octave_idx_type pdim = f->geometric_dimension (); if (nargin != pdim + 1) print_usage (); else { std::vector <Matrix> coordinates; dim_vector dims; for (octave_idx_type in = 1; in <= pdim; ++in) { if (! args(in).is_real_type ()) error ("invalid argument"); else { Matrix aux = args(in).matrix_value (); if (in == 1) dims = aux.dims (); else { dim_vector newdims = aux.dims (); if (dims != newdims) { std::string msg = "all the input matrices should"; msg += " have the same size"; error (msg.c_str ()); break; } } coordinates.push_back (aux); } } if (! error_state) { octave_idx_type vdim = f->value_dimension (0); if (nargout != vdim) error ("wrong number of output arguments"); else { std::vector <Matrix> evaluations; for (octave_idx_type out = 0; out < vdim; ++out) evaluations.push_back (Matrix (dims)); for (octave_idx_type k = 0; k < dims.numel (); ++k) { Array<double> point (dim_vector (pdim, 1)); for (octave_idx_type el = 0; el < pdim; ++el) point (el) = coordinates[el] (k); dolfin::Array<double> x (point.length (), point.fortran_vec ()); Array<double> res (dim_vector (vdim, 1)); dolfin::Array<double> values (res.length (), res.fortran_vec ()); try { f->eval (values, x); } catch (std::runtime_error & err) { std::string msg = "cannot evaluate a function"; msg += " outside of its domain"; error (msg.c_str ()); break; } for (octave_idx_type el = 0; el < vdim; ++el) evaluations[el] (k) = res (el); } if (! error_state) { for (std::vector<Matrix>::iterator it = evaluations.begin (); it != evaluations.end (); ++it) retval.append (octave_value (*it)); } } } } } } return retval; }