Mercurial > octave
view scripts/plot/draw/streamribbon.m @ 32632:2e484f9f1f18 stable
maint: update Octave Project Developers copyright for the new year
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 22 Dec 2023 12:08:17 -0500 |
parents | 597f3ee61a48 |
children |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2020-2024 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {} streamribbon (@var{x}, @var{y}, @var{z}, @var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz}) ## @deftypefnx {} {} streamribbon (@var{u}, @var{v}, @var{w}, @var{sx}, @var{sy}, @var{sz}) ## @deftypefnx {} {} streamribbon (@var{xyz}, @var{x}, @var{y}, @var{z}, @var{anlr_spd}, @var{lin_spd}) ## @deftypefnx {} {} streamribbon (@var{xyz}, @var{anlr_spd}, @var{lin_spd}) ## @deftypefnx {} {} streamribbon (@var{xyz}, @var{anlr_rot}) ## @deftypefnx {} {} streamribbon (@dots{}, @var{width}) ## @deftypefnx {} {} streamribbon (@var{hax}, @dots{}) ## @deftypefnx {} {@var{h} =} streamribbon (@dots{}) ## Calculate and display streamribbons. ## ## The streamribbon is constructed by rotating a normal vector around a ## streamline according to the angular rotation of the vector field. ## ## The vector field is given by @code{[@var{u}, @var{v}, @var{w}]} and is ## defined over a rectangular grid given by @code{[@var{x}, @var{y}, @var{z}]}. ## The streamribbons start at the seed points ## @code{[@var{sx}, @var{sy}, @var{sz}]}. ## ## @code{streamribbon} can be called with a cell array that contains ## pre-computed streamline data. To do this, @var{xyz} must be created with ## the @code{stream3} function. @var{lin_spd} is the linear speed of the ## vector field and can be calculated from @code{[@var{u}, @var{v}, @var{w}]} ## by the square root of the sum of the squares. The angular speed ## @var{anlr_spd} is the projection of the angular velocity onto the velocity ## of the normalized vector field and can be calculated with the @code{curl} ## command. This option is useful if you need to alter the integrator step ## size or the maximum number of streamline vertices. ## ## Alternatively, ribbons can be created from an array of vertices @var{xyz} of ## a path curve. @var{anlr_rot} contains the angles of rotation around the ## edges between adjacent vertices of the path curve. ## ## The input parameter @var{width} sets the width of the streamribbons. ## ## Streamribbons are colored according to the total angle of rotation along the ## ribbon. ## ## If the first argument @var{hax} is an axes handle, then plot into this axes, ## rather than the current axes returned by @code{gca}. ## ## The optional return value @var{h} is a graphics handle to the plot objects ## created for each streamribbon. ## ## Example: ## ## @example ## @group ## [x, y, z] = meshgrid (0:0.2:4, -1:0.2:1, -1:0.2:1); ## u = - x + 10; ## v = 10 * z.*x; ## w = - 10 * y.*x; ## streamribbon (x, y, z, u, v, w, [0, 0], [0, 0.6], [0, 0]); ## view (3); ## @end group ## @end example ## ## @seealso{streamline, stream3, streamtube, ostreamtube} ## ## @end deftypefn ## References: ## ## @inproceedings{ ## title = {Feature Detection from Vector Quantities in a Numerically Simulated Hypersonic Flow Field in Combination with Experimental Flow Visualization}, ## author = {Pagendarm, Hans-Georg and Walter, Birgit}, ## year = {1994}, ## publisher = {IEEE Computer Society Press}, ## booktitle = {Proceedings of the Conference on Visualization ’94}, ## pages = {117–123}, ## } ## ## @article{ ## title = {Efficient streamline, streamribbon, and streamtube constructions on unstructured grids}, ## author = {Ueng, Shyh-Kuang and Sikorski, C. and Ma, Kwan-Liu}, ## year = {1996}, ## month = {June}, ## publisher = {IEEE Transactions on Visualization and Computer Graphics}, ## } ## ## @inproceedings{ ## title = {Visualization of 3-D vector fields - Variations on a stream}, ## author = {Dave Darmofal and Robert Haimes}, ## year = {1992} ## } ## ## @techreport{ ## title = {Parallel Transport Approach to Curve Framing}, ## author = {Andrew J. Hanson and Hui Ma}, ## year = {1995} ## } ## ## @article{ ## title = {There is More than One Way to Frame a Curve}, ## author = {Bishop, Richard}, ## year = {1975}, ## month = {03}, ## volume = {82}, ## publisher = {The American Mathematical Monthly} ## } function h = streamribbon (varargin) [hax, varargin, nargin] = __plt_get_axis_arg__ ("streamribbon", varargin{:}); width = []; xyz = []; anlr_spd = []; lin_spd = []; anlr_rot = []; switch (nargin) case 2 [xyz, anlr_rot] = varargin{:}; case 3 if (numel (varargin{3}) == 1) [xyz, anlr_rot, width] = varargin{:}; else [xyz, anlr_spd, lin_spd] = varargin{:}; [m, n, p] = size (anlr_spd); [x, y, z] = meshgrid (1:n, 1:m, 1:p); endif case 4 [xyz, anlr_spd, lin_spd, width] = varargin{:}; [m, n, p] = size (anlr_spd); [x, y, z] = meshgrid (1:n, 1:m, 1:p); case 6 if (iscell (varargin{1})) [xyz, x, y, z, anlr_spd, lin_spd] = varargin{:}; else [u, v, w, spx, spy, spz] = varargin{:}; [m, n, p] = size (u); [x, y, z] = meshgrid (1:n, 1:m, 1:p); endif case 7 if (iscell (varargin{1})) [xyz, x, y, z, anlr_spd, lin_spd, width] = varargin{:}; else [u, v, w, spx, spy, spz, width] = varargin{:}; [m, n, p] = size (u); [x, y, z] = meshgrid (1:n, 1:m, 1:p); endif case 9 [x, y, z, u, v, w, spx, spy, spz] = varargin{:}; case 10 [x, y, z, u, v, w, spx, spy, spz, width] = varargin{:}; otherwise print_usage (); endswitch if (isempty (xyz)) xyz = stream3 (x, y, z, u, v, w, spx, spy, spz); anlr_spd = curl (x, y, z, u, v, w); lin_spd = sqrt (u.*u + v.*v + w.*w); endif ## Derive scale factor from the bounding box diagonal if (isempty (width)) mxx = mnx = mxy = mny = mxz = mnz = []; j = 1; for i = 1 : length (xyz) sl = xyz{i}; if (! isempty (sl)) slx = sl(:,1); sly = sl(:,2); slz = sl(:,3); mxx(j) = max (slx); mnx(j) = min (slx); mxy(j) = max (sly); mny(j) = min (sly); mxz(j) = max (slz); mnz(j) = min (slz); j += 1; endif endfor dx = max (mxx) - min (mnx); dy = max (mxy) - min (mny); dz = max (mxz) - min (mnz); width = sqrt (dx*dx + dy*dy + dz*dz) / 25; elseif (! isreal (width) || width <= 0) error ("streamribbon: WIDTH must be a real scalar > 0"); endif if (! isempty (anlr_rot)) for i = 1 : length (xyz) if (rows (anlr_rot{i}) != rows (xyz{i})) error ("streamribbon: ANLR_ROT must have same length as XYZ"); endif endfor endif if (isempty (hax)) hax = gca (); else hax = hax(1); endif ## Angular speed of a paddle wheel spinning around a streamline in a fluid ## flow "V": ## dtheta/dt = 0.5 * <curl(V), V/norm(V)> ## ## Integration along a streamline segment with the length "h" yields the ## rotation angle: ## theta = 0.25 * h * <curl(V), V(0)/norm(V(0))^2) + V(h)/norm(V(h))^2)> ## ## Alternative approach using the curl angular speed "c = curl()": ## theta = 0.5 * h * (c(0)/norm(V(0)) + c(h)/norm(V(h))) ## ## Hints: ## i. ) For integration use trapezoidal rule ## ii.) "V" can be assumend to be piecewise linear and curl(V) to be ## piecewise constant because of the used linear interpolation h = []; for i = 1 : length (xyz) sl = xyz{i}; num_vertices = rows (sl); if (! isempty (sl) && num_vertices > 1) if (isempty (anlr_rot)) ## Plot from vector field ## Interpolate onto streamline vertices [lin_spd_sl, anlr_spd_sl, max_vertices] = ... interp_sl (x, y, z, lin_spd, anlr_spd, sl); if (max_vertices > 1) ## Euclidean distance between two adjacent vertices stepsize = vecnorm (diff (sl(1:max_vertices, :)), 2, 2); ## Angular rotation around edges between two adjacent sl-vertices ## Note: Potential "division by zero" is checked in interp_sl() anlr_rot_sl = 0.5 * stepsize.*(anlr_spd_sl(1:max_vertices - 1)./ ... lin_spd_sl(1:max_vertices - 1) + ... anlr_spd_sl(2:max_vertices)./ ... lin_spd_sl(2:max_vertices)); htmp = plotribbon (hax, sl, anlr_rot_sl, max_vertices, 0.5 * width); h = [h; htmp]; endif else ## Plot from vertice array anlr_rot_sl = anlr_rot{i}; htmp = plotribbon (hax, sl, anlr_rot_sl, num_vertices, 0.5 * width); h = [h; htmp]; endif endif endfor endfunction function h = plotribbon (hax, sl, anlr_rot_sl, max_vertices, width2) total_angle = cumsum (anlr_rot_sl); total_angle = [0; total_angle]; ## 1st streamline segment X0 = sl(1,:); X1 = sl(2,:); R = X1 - X0; RE = R / norm (R); ## Initial vector KE which is to be transported along the vertice array KE = get_normal2 (RE); XS10 = - width2 * KE + X0; XS20 = width2 * KE + X0; ## Apply angular rotation cp = cos (anlr_rot_sl(1)); sp = sin (anlr_rot_sl(1)); KE = rotation (KE, RE, cp, sp).'; XS1 = - width2 * KE + X1; XS2 = width2 * KE + X1; px = zeros (2, max_vertices); py = zeros (2, max_vertices); pz = zeros (2, max_vertices); pc = zeros (2, max_vertices); px(:,1) = [XS10(1); XS20(1)]; py(:,1) = [XS10(2); XS20(2)]; pz(:,1) = [XS10(3); XS20(3)]; pc(:,1) = total_angle(1) * [1; 1]; px(:,2) = [XS1(1); XS2(1)]; py(:,2) = [XS1(2); XS2(2)]; pz(:,2) = [XS1(3); XS2(3)]; pc(:,2) = total_angle(2) * [1; 1]; for i = 3 : max_vertices ## Next streamline segment X0 = X1; X1 = sl(i,:); R = X1 - X0; RE = R / norm (R); ## Project KE onto RE and get the difference in order to transport ## the normal vector KE along the vertex array Kp = KE - RE * dot (KE, RE); KE = Kp / norm (Kp); ## Apply angular rotation to KE cp = cos (anlr_rot_sl(i - 1)); sp = sin (anlr_rot_sl(i - 1)); KE = rotation (KE, RE, cp, sp).'; XS1 = - width2 * KE + X1; XS2 = width2 * KE + X1; px(:,i) = [XS1(1); XS2(1)]; py(:,i) = [XS1(2); XS2(2)]; pz(:,i) = [XS1(3); XS2(3)]; pc(:,i) = total_angle(i) * [1; 1]; endfor h = surface (hax, px, py, pz, pc); endfunction ## Interpolate speed and divergence onto the streamline vertices and ## return the first chunck of valid samples until a singularity / ## zero is hit or the streamline vertex array "sl" ends function [lin_spd_sl, anlr_spd_sl, max_vertices] = ... interp_sl (x, y, z, lin_spd, anlr_spd, sl) anlr_spd_sl = interp3 (x, y, z, anlr_spd, sl(:,1), sl(:,2), sl(:,3)); lin_spd_sl = interp3 (x, y, z, lin_spd, sl(:,1), sl(:,2), sl(:,3)); is_singular_anlr_spd = find (isnan (anlr_spd_sl), 1, "first"); is_zero_lin_spd = find (lin_spd_sl == 0, 1, "first"); max_vertices = rows (sl); if (! isempty (is_singular_anlr_spd)) max_vertices = min (max_vertices, is_singular_anlr_spd - 1); endif if (! isempty (is_zero_lin_spd)) max_vertices = min (max_vertices, is_zero_lin_spd - 1); endif endfunction ## N normal to X, so that N is in span ([0 0 1], X) ## If not possible then span ([1 0 0], X) function N = get_normal2 (X) if ((X(1) == 0) && (X(2) == 0)) A = [1, 0, 0]; else A = [0, 0, 1]; endif ## Project A onto X and get the difference N = A - X * dot (A, X) / (norm (X)^2); N /= norm (N); endfunction ## Rotate X around U where |U| = 1 ## cp = cos (angle), sp = sin (angle) function Y = rotation (X, U, cp, sp) ux = U(1); uy = U(2); uz = U(3); Y(1,:) = X(1) * (cp + ux * ux * (1 - cp)) + ... X(2) * (ux * uy * (1 - cp) - uz * sp) + ... X(3) * (ux * uz * (1 - cp) + uy * sp); Y(2,:) = X(1) * (uy * ux * (1 - cp) + uz * sp) + ... X(2) * (cp + uy * uy * (1 - cp)) + ... X(3) * (uy * uz * (1 - cp) - ux * sp); Y(3,:) = X(1) * (uz * ux * (1 - cp) - uy * sp) + ... X(2) * (uz * uy * (1 - cp) + ux * sp) + ... X(3) * (cp + uz * uz * (1 - cp)); endfunction %!demo %! clf; %! [x, y, z] = meshgrid (0:0.2:4, -1:0.2:1, -1:0.2:1); %! u = - x + 10; %! v = 10 * z.*x; %! w = - 10 * y.*x; %! sx = [0, 0]; %! sy = [0, 0.6]; %! sz = [0, 0]; %! streamribbon (x, y, z, u, v, w, sx, sy, sz); %! hold on; %! quiver3 (x, y, z, u, v, w); %! colormap (jet); %! shading interp; %! camlight ("headlight"); %! view (3); %! axis tight equal off; %! set (gca, "cameraviewanglemode", "manual"); %! hcb = colorbar; %! title (hcb, "Angle"); %! title ("Streamribbon"); %!demo %! clf; %! t = (0:pi/50:2*pi).'; %! xyz{1} = [cos(t), sin(t), 0*t]; %! twist{1} = ones (numel (t), 1) * pi / (numel (t) - 1); %! streamribbon (xyz, twist, 0.5); %! colormap (jet); %! view (3); %! camlight ("headlight"); %! axis tight equal off; %! title ("Moebius Strip"); ## Test input validation %!error <Invalid call> streamribbon () %!error <Invalid call> streamribbon (1) %!error <Invalid call> streamribbon (1,2,3,4,5) %!error <Invalid call> streamribbon (1,2,3,4,5,6,7,8) %!error <Invalid call> streamribbon (1,2,3,4,5,6,7,8,9,10,11) %!error <WIDTH must be a real scalar . 0> streamribbon (1,2,3,1i) %!error <WIDTH must be a real scalar . 0> streamribbon (1,2,3,0) %!error <WIDTH must be a real scalar . 0> streamribbon (1,2,3,-1) %!error <ANLR_ROT must have same length as XYZ> streamribbon ({[1,1,1;2,2,2]},{[1,1,1]})