|
11171
|
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}} = bim2c_intrp (@var{msh}, @var{n_data}, @var{e_data}, @var{points}) |
|
|
20 ## |
|
12628
|
21 ## Compute interpolated values of multicomponent node centered field @var{n_data} and/or |
|
|
22 ## cell centered field @var{n_data} at an arbitrary set of points whose coordinates are given in the |
|
|
23 ## n_by_2 matrix @var{points}. |
|
11171
|
24 ## |
|
|
25 ## @end deftypefn |
|
|
26 |
|
|
27 |
|
|
28 ## Author: Carlo de Falco <carlo@guglielmo.local> |
|
|
29 ## Created: 2012-10-01 |
|
|
30 |
|
|
31 function data = bim2c_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:3, :)', 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:3, t_list), :), |
|
|
47 [3, 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 |
|
|
55 endfunction |
|
|
56 |
|
|
57 %!test |
|
|
58 %! msh = bim2c_mesh_properties (msh2m_structured_mesh (linspace (0, 1, 11), linspace (0, 1, 13), 1, 1:4)); |
|
|
59 %! x = y = linspace (0, 1, 100).'; |
|
|
60 %! u = msh.p(1, :).'; |
|
|
61 %! ui = bim2c_intrp (msh, u, [], [x, y]); |
|
12628
|
62 %! assert (ui, linspace (0, 1, 100), 10*eps); |
