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);