Mercurial > octave

comparison scripts/plot/draw/ostreamtube.m @ 28135:695bb31e565b stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename "streamtube" to "ostreamtube" (bug #57471). * ostreamtube.m: Rename from "streamtube.m" ahead of adding a Matlab compatible "streamtube" function. * stream3.m, streamline.m, module.mk, plot.txi, NEWS: Change all references.
author Markus Meisinger <chloros2@gmx.de>
date Wed, 19 Feb 2020 07:53:10 +0100
parents scripts/plot/draw/streamtube.m@465be7c652f1
children 23f667483fab
comparison
equal deleted inserted replaced
28131:27c99ff83b99 28135:695bb31e565b
1 ########################################################################
2 ##
3 ## Copyright (C) 2019-2020 The Octave Project Developers
4 ##
5 ## See the file COPYRIGHT.md in the top-level directory of this
6 ## distribution or <https://octave.org/copyright/>.
7 ##
8 ## This file is part of Octave.
9 ##
10 ## Octave is free software: you can redistribute it and/or modify it
11 ## under the terms of the GNU General Public License as published by
12 ## the Free Software Foundation, either version 3 of the License, or
13 ## (at your option) any later version.
14 ##
15 ## Octave is distributed in the hope that it will be useful, but
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ## GNU General Public License for more details.
19 ##
20 ## You should have received a copy of the GNU General Public License
21 ## along with Octave; see the file COPYING. If not, see
22 ## <https://www.gnu.org/licenses/>.
23 ##
24 ########################################################################
25
26 ## -*- texinfo -*-
27 ## @deftypefn {} {} ostreamtube (@var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
28 ## @deftypefnx {} {} ostreamtube (@var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz})
29 ## @deftypefnx {} {} ostreamtube (@var{xyz}, @var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w})
30 ## @deftypefnx {} {} ostreamtube (@dots{}, @var{options})
31 ## @deftypefnx {} {} ostreamtube (@var{hax}, @dots{})
32 ## @deftypefnx {} {@var{h} =} ostreamtube (@dots{})
33 ## Calculate and display streamtubes.
34 ##
35 ## Streamtubes are approximated by connecting circular crossflow areas
36 ## along a streamline. The expansion of the flow is determined by the local
37 ## crossflow divergence.
38 ##
39 ## The vector field is given by @code{[@var{u}, @var{v}, @var{w}]} and is
40 ## defined over a rectangular grid given by @code{[@var{x}, @var{y}, @var{z}]}.
41 ## The streamtubes start at the seed points
42 ## @code{[@var{sx}, @var{sy}, @var{sz}]}.
43 ##
44 ## The tubes are colored based on the local vector field strength.
45 ##
46 ## The input parameter @var{options} is a 2-D vector of the form
47 ## @code{[@var{scale}, @var{n}]}. The first parameter scales the start radius
48 ## of the streamtubes (default 1). The second parameter specifies the number
49 ## of vertices that are used to construct the tube circumference (default 20).
50 ##
51 ## @code{ostreamtube} can be called with a cell array containing pre-computed
52 ## streamline data. To do this, @var{xyz} must be created with the
53 ## @code{stream3} function. This option is useful if you need to alter the
54 ## integrator step size or the maximum number of vertices of the streamline.
55 ##
56 ## If the first argument @var{hax} is an axes handle, then plot into this axes,
57 ## rather than the current axes returned by @code{gca}.
58 ##
59 ## The optional return value @var{h} is a graphics handle to the plot
60 ## objects created for each streamtube.
61 ##
62 ## Example:
63 ##
64 ## @example
65 ## @group
66 ## [x, y, z] = meshgrid (-1:0.1:1, -1:0.1:1, -3:0.1:0);
67 ## u = -x / 10 - y;
68 ## v = x - y / 10;
69 ## w = - ones (size (x)) / 10;
70 ## ostreamtube (x, y, z, u, v, w, 1, 0, 0);
71 ## @end group
72 ## @end example
73 ##
74 ## @seealso{stream3, streamline}
75 ## @end deftypefn
76
77 ## References:
78 ##
79 ## @inproceedings{
80 ## title = {Visualization of 3-D vector fields - Variations on a stream},
81 ## author = {Dave Darmofal and Robert Haimes},
82 ## year = {1992}
83 ## }
84 ##
85 ## @article{
86 ## title = {Efficient streamline, streamribbon, and streamtube constructions on unstructured grids},
87 ## author = {Ueng, Shyh-Kuang and Sikorski, C. and Ma, Kwan-Liu},
88 ## year = {1996},
89 ## month = {June},
90 ## publisher = {IEEE Transactions on Visualization and Computer Graphics},
91 ## }
92
93 function h = ostreamtube (varargin)
94
95 [hax, varargin, nargin] = __plt_get_axis_arg__ ("ostreamtube", varargin{:});
96
97 options = [];
98 xyz = [];
99 switch (nargin)
100 case 0
101 print_usage ();
102 case 6
103 [u, v, w, spx, spy, spz] = varargin{:};
104 [m, n, p] = size (u);
105 [x, y, z] = meshgrid (1:n, 1:m, 1:p);
106 case 7
107 if (iscell (varargin{1}))
108 [xyz, x, y, z, u, v, w] = varargin{:};
109 else
110 [u, v, w, spx, spy, spz, options] = varargin{:};
111 [m, n, p] = size (u);
112 [x, y, z] = meshgrid (1:n, 1:m, 1:p);
113 endif
114 case 8
115 [xyz, x, y, z, u, v, w, options] = varargin{:};
116 case 9
117 [x, y, z, u, v, w, spx, spy, spz] = varargin{:};
118 case 10
119 [x, y, z, u, v, w, spx, spy, spz, options] = varargin{:};
120 otherwise
121 error ("ostreamtube: invalid number of inputs");
122 endswitch
123
124 scale = 1;
125 num_circum = 20;
126 if (! isempty (options))
127 switch (numel (options))
128 case 1
129 scale = options(1);
130 case 2
131 scale = options(1);
132 num_circum = options(2);
133 otherwise
134 error ("ostreamtube: invalid number of OPTIONS elements");
135 endswitch
136
137 if (! isreal (scale) || scale <= 0)
138 error ("ostreamtube: SCALE must be a real scalar > 0");
139 endif
140 if (! isreal (num_circum) || num_circum < 3)
141 error ("ostreamtube: number of tube vertices N must be greater than 2");
142 endif
143 num_circum = fix (num_circum);
144 endif
145
146 if (isempty (hax))
147 hax = gca ();
148 else
149 hax = hax(1);
150 endif
151
152 if (isempty (xyz))
153 xyz = stream3 (x, y, z, u, v, w, spx, spy, spz, 0.2);
154 endif
155
156 div = divergence (x, y, z, u, v, w);
157
158 ## Use the bounding box diagonal to determine the starting radius
159 mxx = mnx = mxy = mny = mxz = mnz = [];
160 j = 1;
161 for i = 1 : length (xyz)
162 sl = xyz{i};
163 if (! isempty (sl))
164 slx = sl(:, 1); sly = sl(:, 2); slz = sl(:, 3);
165 mxx(j) = max (slx); mnx(j) = min (slx);
166 mxy(j) = max (sly); mny(j) = min (sly);
167 mxz(j) = max (slz); mnz(j) = min (slz);
168 j += 1;
169 endif
170 endfor
171 dx = max (mxx) - min (mnx);
172 dy = max (mxy) - min (mny);
173 dz = max (mxz) - min (mnz);
174 rstart = scale * sqrt (dx*dx + dy*dy + dz*dz) / 25;
175
176 h = [];
177 for i = 1 : length (xyz)
178 sl = xyz{i};
179 num_vertices = rows (sl);
180 if (! isempty (sl) && num_vertices > 2)
181
182 usl = interp3 (x, y, z, u, sl(:, 1), sl(:, 2), sl(:, 3));
183 vsl = interp3 (x, y, z, v, sl(:, 1), sl(:, 2), sl(:, 3));
184 wsl = interp3 (x, y, z, w, sl(:, 1), sl(:, 2), sl(:, 3));
185 vv = sqrt (usl.*usl + vsl.*vsl + wsl.*wsl);
186
187 div_sl = interp3 (x, y, z, div, sl(:, 1), sl(:, 2), sl(:, 3));
188 is_singular_div = find (isnan (div_sl), 1, "first");
189
190 if (! isempty (is_singular_div))
191 max_vertices = is_singular_div - 1;
192 else
193 max_vertices = num_vertices;
194 endif
195
196 if (max_vertices > 2)
197
198 htmp = plottube (hax, sl, div_sl, vv, max_vertices, ...
199 rstart, num_circum);
200 h = [h; htmp];
201
202 endif
203 endif
204 endfor
205
206 endfunction
207
208 function h = plottube (hax, sl, div_sl, vv, max_vertices, rstart, num_circum)
209
210 phi = linspace (0, 2*pi, num_circum);
211 cp = cos (phi);
212 sp = sin (phi);
213
214 X0 = sl(1, :);
215 X1 = sl(2, :);
216
217 ## 1st rotation axis
218 R = X1 - X0;
219 RE = R / norm (R);
220
221 ## Initial radius
222 vold = vv(1);
223 vact = vv(2);
224 ract = rstart * exp (0.5 * div_sl(2) * norm (R) / vact) * sqrt (vold / vact);
225 vold = vact;
226 rold = ract;
227
228 ## Guide point and its rotation to create a segment
229 N = get_normal1 (R);
230 K = ract * N;
231 XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_circum);
232
233 px = zeros (num_circum, max_vertices - 1);
234 py = zeros (num_circum, max_vertices - 1);
235 pz = zeros (num_circum, max_vertices - 1);
236 pc = zeros (num_circum, max_vertices - 1);
237
238 px(:, 1) = XS(1, :).';
239 py(:, 1) = XS(2, :).';
240 pz(:, 1) = XS(3, :).';
241 pc(:, 1) = vact * ones (num_circum, 1);
242
243 for i = 3 : max_vertices
244
245 KK = K;
246 X0 = X1;
247 X1 = sl(i, :);
248 R = X1 - X0;
249 RE = R / norm (R);
250
251 ## Tube radius
252 vact = vv(i);
253 ract = rold * exp (0.5 * div_sl(i) * norm (R) / vact) * sqrt (vold / vact);
254 vold = vact;
255 rold = ract;
256
257 ## Project KK onto RE and get the difference in order to calculate the next
258 ## guiding point
259 Kp = KK - RE * dot (KK, RE);
260 K = ract * Kp / norm (Kp);
261
262 ## Rotate around RE and collect surface patches
263 XS = rotation (K, RE, cp, sp) + repmat (X1.', 1, num_circum);
264
265 px(:, i - 1) = XS(1, :).';
266 py(:, i - 1) = XS(2, :).';
267 pz(:, i - 1) = XS(3, :).';
268 pc(:, i - 1) = vact * ones (num_circum, 1);
269
270 endfor
271
272 h = surface (hax, px, py, pz, pc);
273
274 endfunction
275
276 ## Arbitrary N normal to X
277 function N = get_normal1 (X)
278
279 if ((X(3) == 0) && (X(1) == -X(2)))
280 N = [- X(2) - X(3), X(1), X(1)];
281 else
282 N = [X(3), X(3), - X(1) - X(2)];
283 endif
284
285 N /= norm (N);
286
287 endfunction
288
289 ## Rotate X around U where |U| = 1
290 ## cp = cos (angle), sp = sin (angle)
291 function Y = rotation (X, U, cp, sp)
292
293 ux = U(1);
294 uy = U(2);
295 uz = U(3);
296
297 Y(1, :) = X(1) * (cp + ux * ux * (1 - cp)) + ...
298 X(2) * (ux * uy * (1 - cp) - uz * sp) + ...
299 X(3) * (ux * uz * (1 - cp) + uy * sp);
300
301 Y(2, :) = X(1) * (uy * ux * (1 - cp) + uz * sp) + ...
302 X(2) * (cp + uy * uy * (1 - cp)) + ...
303 X(3) * (uy * uz * (1 - cp) - ux * sp);
304
305 Y(3, :) = X(1) * (uz * ux * (1 - cp) - uy * sp) + ...
306 X(2) * (uz * uy * (1 - cp) + ux * sp) + ...
307 X(3) * (cp + uz * uz * (1 - cp));
308
309 endfunction
310
311 %!demo
312 %! clf;
313 %! [x, y, z] = meshgrid (-1:0.1:1, -1:0.1:1, -3.5:0.1:0);
314 %! a = 0.1;
315 %! b = 0.1;
316 %! u = - a * x - y;
317 %! v = x - a * y;
318 %! w = - b * ones (size (x));
319 %! sx = 1.0;
320 %! sy = 0.0;
321 %! sz = 0.0;
322 %! ostreamtube (x, y, z, u, v, w, sx, sy, sz, [1.2, 30]);
323 %! colormap (jet);
324 %! shading interp;
325 %! view ([-47, 24]);
326 %! camlight ();
327 %! lighting gouraud;
328 %! grid on;
329 %! view (3);
330 %! axis equal;
331 %! set (gca, "cameraviewanglemode", "manual");
332 %! title ("Spiral Sink");
333
334 %!demo
335 %! clf;
336 %! [x, y, z] = meshgrid (-2:0.5:2);
337 %! t = sqrt (1.0./(x.^2 + y.^2 + z.^2)).^3;
338 %! u = - x.*t;
339 %! v = - y.*t;
340 %! w = - z.*t;
341 %! [sx, sy, sz] = meshgrid (-2:4:2);
342 %! xyz = stream3 (x, y, z, u, v, w, sx, sy, sz, [0.1, 60]);
343 %! ostreamtube (xyz, x, y, z, u, v, w, [2, 50]);
344 %! colormap (jet);
345 %! shading interp;
346 %! view ([-47, 24]);
347 %! camlight ();
348 %! lighting gouraud;
349 %! grid on;
350 %! view (3);
351 %! axis equal;
352 %! set (gca, "cameraviewanglemode", "manual");
353 %! title ("Integration Towards Sink");
354
355 ## Test input validation
356 %!error ostreamtube ()
357 %!error <invalid number of inputs> ostreamtube (1)
358 %!error <invalid number of inputs> ostreamtube (1,2)
359 %!error <invalid number of inputs> ostreamtube (1,2,3)
360 %!error <invalid number of inputs> ostreamtube (1,2,3,4)
361 %!error <invalid number of inputs> ostreamtube (1,2,3,4,5)
362 %!error <invalid number of OPTIONS> ostreamtube (1,2,3,4,5,6, [1,2,3])
363 %!error <SCALE must be a real scalar . 0> ostreamtube (1,2,3,4,5,6, [1i])
364 %!error <SCALE must be a real scalar . 0> ostreamtube (1,2,3,4,5,6, [0])
365 %!error <N must be greater than 2> ostreamtube (1,2,3,4,5,6, [1, 1i])
366 %!error <N must be greater than 2> ostreamtube (1,2,3,4,5,6, [1, 2])