femwind/wind_streamlines.m

   1 function wind_speed = wind_streamline(X, CX, W, params)
   2
   3 %disp('wind_streamlines not done yet')
   4 %Extracting min and max values from mesh grid to create a valid n-d grid
   5 %for streamlines function
   6 % Most work goes into creating a grid space with
   7 %Constructing the wind field interpolant
   8 %FX,
   9
  10 [F1, F2, F3] = wind_interp(W,CX);
  11  [n(1),n(2),n(3)]=size(X{1});
  12 F1.Method = 'natural';
  13 F2.Method = 'natural';
  14 F3.Method = 'natural';
  15 level = 1:n(3);
  16
  17 if params.st_contour == 1
  18     xmin = min(X{1}(:));
  19     xmax = max(X{1}(:));
  20     ymin = min(X{2}(:));
  21     ymax = max(X{2}(:));
  22     zmin = min(X{3}(:));
  23     zmax = max(X{3}(:));
  24     %Creating nd-meshgrid
  25
  26
  27     %The slice function requires a meshgrid to work, while X is an ndgrid, so all of the coordinates are transposed.
  28     %Create a meshgrid with the same array and physical dimensions as the domain arrays
  29     %CX. Evaluate each grid point at the scattered interpolant
  30     [CXG_X,CXG_Y,CXG_Z] = meshgrid(linspace(min(X{2}(:)),max(X{2}(:)), n(2)),...
  31         linspace(min(X{1}(:)),max(X{1}(:)), n(1)), ...
  32         linspace(min(X{3}(:)),max(X{3}(:)), n(3)));
  33
  34     WX = F1(CXG_Y(:), CXG_X(:), CXG_Z(:));
  35     WY = F2(CXG_Y(:), CXG_X(:), CXG_Z(:));
  36     WZ = F3(CXG_Y(:), CXG_X(:), CXG_Z(:));
  37
  38
  39     %Convert individual wind vectors (x,y,z) into a full grid with 3-D wind
  40     %components at each coordinate in the 3-D grid space.
  41     WX = reshape(WX,n);
  42     WY = reshape(WY,n);
  43     WZ = reshape(WZ,n);
  44
  45
  46     wind_speed = sqrt(WX.^2 + WY.^2 + WZ.^2);
  47     %Normalizing wind speed to assist in interpolation
  48     %wind_speed = wind_speed/max(wind_speed(:));
  49
  50     %Creating contour surfaces and colormaps of the magnitude of the
  51     %wind-field at 3 different locations
  52     hsurfaces = slice(CXG_X,CXG_Y,CXG_Z,wind_speed, [ymin, (ymax-ymin)/2, ymax], xmax,zmin);
  53     set(hsurfaces,'FaceColor','interp','EdgeColor','interp');
  54     colormap(turbo)
  55 %     shading interp
  56 %     view(3); axis vis3d; camlight;
  57     %else
  58     %colormap jet
  59     %end
  60     hcont = ...
  61         contourslice(CXG_X,CXG_Y,CXG_Z,wind_speed,[ymin, (ymax-ymin)/2, ymax], xmax,zmin);
  62     set(hcont,'EdgeColor',[0.7 0.7 0.7],'LineWidth',0.5);
  63     % % drawing on current figure
  64     %Setting up colorbar
  65
  66     c = colorbar;
  67     c.TickLabelInterpreter = 'tex';
  68     c.Location = 'southoutside';
  69     c.Label.String = 'Wind Speed [m/s]';
  70 %     c.Limits = [0,300];
  71 %     c.Ticks = [0, 150, 300];
  72 %     c.TickLabels = {'0', '5', '10'};
  73
  74
  75
  76 end
  77     %Note: To use the scattered interpolant function arrays
  78     %into column-vector format
  79     if ~exist('scale_t','var') || isempty(scale_t)
  80         scale_t = 1.1;
  81     end
  82
  83     if ~exist('t_final','var') || isempty(t_final)
  84         t_final = scale_t*max(X{1}(:))/mean(W{1}(:));
  85     end
  86
  87      tspan = [0 t_final];
  88
  89      x0 = 0;
  90      y0 = X{2}(1,5,1): (n(2) - 10) :X{2}(1,n(2) - 5,1);
  91      z0 = params.in_height_stream;
  92
  93
  94     %plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,2,2])
  95 if ~exist('scale','var')
  96     scale = 0.9;
  97 end
  98
  99 for i = 1:length(y0)
 100     for j = 1:length(z0)
 101         [t,y] = ode45(@(t,y) odefun(t,y, F1, F2, F3), tspan, [x0,y0(i),z0(j)]);
 102         % Using stiff ODE23 solver
 103         % [t,y] = ode23s(@(t,y) odefun(t,y, F1, F2, F3), tspan, [x0,y0(i),z0(j)]);
 104         % Find indices that constrain the streamline to stay inside the
 105         % domain of the mesh domain.
 106         count = 2;
 107         if max(y(:,1)) < max(X{1}(:))
 108
 109 %             t0 = 0;
 110 %             tf = count*scale_t*max(X{1}(:))/mean(W{1}(:));
 111 %             t = [t0, tf];
 112             L = length(y(:,1));
 113             [t,y_new] = ode23s(@(t,y_new) odefun(t,y_new, F1, F2, F3), tspan, [y(L,1),y(L,2),y(L,3)]);
 114
 115             y = cat(1,y,y_new);
 116
 117             count = count+ 1
 118         end
 119         count  = 0;
 120         ind1 = find(y(:,1)< max(X{1}(:)));
 121         ind1 = ind1(1: length(ind1) );
 122 %         ind2 = find(y(:,2)< max(X{2}(:)));
 123 %         ind2 = ind2(1: length(ind2) - 1);
 124 %         if length(ind1) < length(ind2)
 125 %             ind1 = ind1;
 126 %         else
 127 %             ind1 = ind2;
 128 %         end
 129
 130
 131         ind1 = find(y(:,1)< max(X{1}(:)));
 132         ind1 = ind1(1: length(ind1) - 1);
 133
 134
 135         hold on
 136
 137         plot3(y(ind1,1),y(ind1,2), y(ind1,3));
 138
 139         %Plot vectors until second to last index inside the bound to keep
 140         %vectors in grid
 141         if params.st_quiver == 1
 142         quiver3(y(ind1,1),y(ind1,2), y(ind1,3), F1(y(ind1,1),y(ind1,2), y(ind1,3)),...
 143              F2(y(ind1,1),y(ind1,2), y(ind1,3)),F3(y(ind1,1),y(ind1,2), y(ind1,3)), scale);
 144         end
 145
 146     end
 147 end
 148 % xlabel('x-Coordinate Axis [m]')
 149 % ylabel('y-Coordinate Axis [m]')
 150 % zlabel('z-Coordinate Axis [m]')
 151 axis tight
 152
 153
 154
 155
 156
 157 %%
 158     function [F1,F2,F3] = wind_interp(W, CX)
 159
 160         F1 = scatteredInterpolant(CX{1}(:), CX{2}(:), CX{3}(:),W{1}(:));
 161         F2 = scatteredInterpolant(CX{1}(:), CX{2}(:), CX{3}(:),W{2}(:));
 162         F3 = scatteredInterpolant(CX{1}(:), CX{2}(:), CX{3}(:),W{3}(:));
 163     end
 164
 165
 166
 167     function dXdt = odefun(t,y, F1, F2, F3)
 168         dXdt = zeros(3,1);
 169         dXdt(1) = F1(y(1),y(2),y(3));
 170         dXdt(2) = F2(y(1),y(2),y(3));
 171         dXdt(3) = F3(y(1),y(2),y(3));
 172     end
 173 end