1## Copyright (C) 2011, 2012 Carlo de Falco 2## 3## This program is free software; you can redistribute it and/or modify 4## it under the terms of the GNU General Public License as published by 5## the Free Software Foundation; either version 3 of the License, or 6## (at your option) any later version. 7## 8## This program is distributed in the hope that it will be useful, 9## but WITHOUT ANY WARRANTY; without even the implied warranty of 10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11## GNU General Public License for more details. 12## 13## You should have received a copy of the GNU General Public License 14## along with Octave; see the file COPYING. If not, see 15## <http://www.gnu.org/licenses/>. 16 17## -*- texinfo -*- 18## 19## @deftypefn {Function File} {@var{data}} = bim3c_intrp (@var{msh}, @var{n_data}, @var{e_data}, @var{points}) 20## 21## Compute interpolated values of node centered multicomponent node centered field @var{n_data} and 22## cell centered field @var{n_data} at an arbitrary set of points whos coordinates are given in the 23## n_by_3 matrix @var{points}. 24## 25## @end deftypefn 26 27 28## Author: Carlo de Falco <carlo@guglielmo.local> 29## Created: 2012-10-01 30 31function data = bim3c_intrp (msh, n_data, e_data, p) 32 33 %% for each point, find the enclosing tetrahedron 34 [t_list, b_list] = tsearchn (msh.p.', msh.t(1:4, :)', p); 35 36 %% only keep points within tetrahedra 37 invalid = isnan (t_list); 38 t_list = t_list (! invalid); 39 ntl = numel (t_list); 40 b_list = b_list(! invalid, :); 41 points(invalid,:) = []; 42 43 data = []; 44 if (! isempty (n_data)) 45 data = cat (1, data, squeeze ( 46 sum (reshape (n_data(msh.t(1:4, t_list), :), 47 [4, ntl, (columns (n_data))]) .* 48 repmat (b_list.', [1, 1, (columns (n_data))]), 1))); 49 endif 50 51 if (! isempty (e_data)) 52 data = cat (1, data, e_data(t_list, :)); 53 endif 54 55endfunction 56 57%!test 58%! msh = bim3c_mesh_properties (msh3m_structured_mesh (linspace (0, 1, 11), linspace (0, 1, 9), linspace (0, 1, 13), 1, 1:6)); 59%! x = y = z = linspace (0, 1, 100).'; 60%! u = msh.p(1, :).'; 61%! ui = bim3c_intrp (msh, u, [], [x, y, z]); 62%! assert (ui, linspace (0, 1, 100), 10*eps);