femwind/sl_cs.m

   1 function stream_cs(X, CX, W, init_z, n_lines, end_z, y0)
   2 %X - Wind mesh grid
   3 %CX Center coordinates of each grid element
   4 %W Wind Array Used For Streamlines
   5 %Parameters set by user
   6 % init_z: Initial height for plotting first streamline
   7 % n_lines: number of streamlines
   8 % end_z:
   9 %W Final wind, or gradient matrix at center points of the grid elements
  10 %disp('wind_streamlines not done yet')
  11 %Extracting min and max values from mesh grid to create a valid n-d grid
  12 %for streamlines function
  13 % Most work goes into creating a grid space with
  14     xmin = min(X{1}(:));
  15     xmax = max(X{1}(:));
  16     ymin = min(X{2}(:));
  17     ymax = max(X{2}(:));
  18     zmin = min(X{3}(:));
  19     zmax = max(X{3}(:));
  20     %Creating nd-meshgrid
  21     [n(1),n(2),n(3)]=size(X{1});
  22
  23     level = 1:n(3);
  24     [F1, F2, F3] = wind_interp(W,CX);
  25     %The slice function requires a meshgrid to work, while X is an ndgrid, so all of the coordinates are transposed.
  26     %Create a meshgrid with the same array and physical dimensions as the domain arrays
  27     %CX. Evaluate each grid point at the scattered interpolant
  28 %     [CXG_X,CXG_Y,CXG_Z] = meshgrid(linspace(min(X{2}(:)),max(X{2}(:)), n(2)),...
  29 %         linspace(min(X{1}(:)),max(X{1}(:)), n(1)), ...
  30 %         linspace(min(X{3}(:)),max(X{3}(:)), n(3)));
  31 %
  32 %     WX = F1(CXG_Y(:), CXG_X(:), CXG_Z(:));
  33 %     WY = F2(CXG_Y(:), CXG_X(:), CXG_Z(:));
  34 %     WZ = F3(CXG_Y(:), CXG_X(:), CXG_Z(:));
  35 %
  36 %
  37 %     %Convert individual wind vectors (x,y,z) into a full grid with 3-D wind
  38 %     %components at each coordinate in the 3-D grid space.
  39 %     WX = reshape(WX,n);
  40 %     WY = reshape(WY,n);
  41 %     WZ = reshape(WZ,n);
  42 %
  43 %
  44 %     wind_speed = sqrt(WX.^2 + WY.^2 + WZ.^2);
  45 %
  46 %     %Creating contour surfaces and colormaps of the magnitude of the
  47 %     %wind-field at 3 different locations
  48 %     hsurfaces = slice(CXG_X,CXG_Y,CXG_Z,wind_speed, [ymin, (ymax-ymin)/2, ymax], xmax,zmin);
  49 %     set(hsurfaces,'FaceColor','interp','EdgeColor','none')
  50 %     colormap turbo
  51 %     hcont = ...
  52 %         contourslice(CXG_X,CXG_Y,CXG_Z,wind_speed,[ymin, (ymax-ymin)/2, ymax], xmax,zmin);
  53 %     set(hcont,'EdgeColor',[0.7 0.7 0.7],'LineWidth',0.5)
  54 %     % % drawing on current figure
  55 %
  56 %     colorbar
  57 %     axis tight
  58 %     title('Streamline Plot with Meshgrid and Wind Color and Contour Map')
  59
  60
  61     %Note: To use the scattered interpolant function arrays
  62     %into column-vector format
  63
  64      tspan = [0 100000];
  65      x0 = 0;
  66      y0 = y0; %(max(X{2}(:)) + min(X{2}(:)))/2 + 1;
  67      %z0 = init_z:n_lines: end_z;
  68      z0 = 2:25: max(X{3}(:));
  69
  70
  71     %plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,2,2])
  72 if ~exist('scale','var')
  73     scale = 0.9;
  74 end
  75 for i = 1:length(y0)
  76     for j = 1:length(z0)
  77         [t,y] = ode45(@(t,y) odefun(t,y, F1, F2, F3), tspan, [x0,y0(i),z0(j)]);
  78         % Find indices that constrain the streamline to stay inside the
  79         % domain of the mesh domain.
  80         k1 = find(y(:,1)< max(X{1}(:)));
  81         k1 = k1(1: length(k1));
  82
  83         hold on
  84         plot(y(k1,1), y(k1,3))
  85
  86         %Plot vectors until second to last index inside the bound to keep
  87         %vectors in grid
  88 %         quiver(y(k1,1), y(k1,3),...
  89 %              F1(y(k1,1),y(k1,2), y(k1,3)),F3(y(k1,1),y(k1,2), y(k1,3)))
  90
  91     title(['XZ-Planar Cross Sectional Profile of Streamlines Positioned at Y=', num2str(y0)])
  92     xlabel('X')
  93     ylabel('Z')
  94     end
  95 end
  96 axis tight
  97
  98
  99 %%
 100     function [F1,F2,F3] = wind_interp(W, CX)
 101
 102         F1 = scatteredInterpolant(CX{1}(:), CX{2}(:), CX{3}(:),W{1}(:));
 103         F2 = scatteredInterpolant(CX{1}(:), CX{2}(:), CX{3}(:),W{2}(:));
 104         F3 = scatteredInterpolant(CX{1}(:), CX{2}(:), CX{3}(:),W{3}(:));
 105     end
 106
 107
 108
 109     function dXdt = odefun(t,y, F1, F2, F3)
 110         dXdt = zeros(3,1);
 111         dXdt(1) = F1(y(1),y(2),y(3));
 112         dXdt(2) = F2(y(1),y(2),y(3));
 113         dXdt(3) = F3(y(1),y(2),y(3));
 114     end
 115 end