%Example code for reading and plotting both the instantaneous and mean flow %fields for a given file, make sure this script is atleast the same %directory level or at most a directory level above. % %NaNMax - This variable is the cause of the difference between the %MeanfromMean and MeanfromInstantaneous. Not all vectors through time are %'valid' and are marked by -9999 in the csv files. NaNMax specifies the %ratio of valid vectors to be able to display the mean value. %% House Keeping clear all % close all clc %% Define file to be read MeanFileName = 'SquareBack_NW_5_0yaw_Z=0194m_FlowField'; % don't include the extension, will be added in later %% Find the files to load in Files = [dir('*'); dir('*\*'); dir('*\*\*')]; FilesCell = struct2cell(Files); FileLocation = contains(FilesCell(1,:),MeanFileName); InstLocation = contains(FilesCell(1,:),'Instantaneous'); %% Percentage to valid vectors required to calculate a mean in the instantaneous files NaNMax = 0.995; %% Mean data, in a try loop as this code may be run with either or both types of data %%Mean try %% Load in the files MeanPath = [Files(FileLocation & ~InstLocation).folder filesep Files(FileLocation & ~InstLocation).name]; MeanData = readtable(MeanPath); %Variables names will look different to what is presented in the .csv file, MATLAB doesn't play nicely with special characters (){}=\/,. in variable names %% Arrange the data for plotting UniqueX = uniquetol(MeanData.x_m_); UniqueY = uniquetol(MeanData.y_m_); UniqueZ = uniquetol(MeanData.z_m_); XSamples = length(UniqueX); YSamples = length(UniqueY); ZSamples = length(UniqueZ); if any([XSamples YSamples ZSamples] == 1) disp('Data requested is planar PIV data') if ZSamples == 1 %This section is required for the Z=0.194m plane as the data %was taken in tunnel coordinates and post processed into model %coordinates. % %As a result, the information on the matrix size has been %'lost' when writing a dataset, so the original matrix size has %been hard coded into this plotting code. It is understood this %isn't good practice. disp('Data Requested is Z=*m') reshapeMatrix = [247 355]; else %Sorting the Planar Data for plotting as MATLAB can be picky MeanData = sortrows(MeanData,[3 2 1]); reshapeMatrix = [XSamples YSamples ZSamples]; reshapeMatrix(reshapeMatrix == 1) = []; %remove singleton dimensions, saves doing squeeze() later. end else disp('Data requested is tomographic PIV data') reshapeMatrix = [XSamples YSamples ZSamples]; %Sorting the Tomographic Data for plotting as MATLAB can be picky MeanData = sortrows(MeanData,[3 2 1]); end PlotX = reshape(MeanData.x_m_,reshapeMatrix); PlotY = reshape(MeanData.y_m_,reshapeMatrix); PlotZ = reshape(MeanData.z_m_,reshapeMatrix); PlotU = reshape(MeanData.u_mean_m_s_,reshapeMatrix); %Reshape the variables for plotting PlotV = reshape(MeanData.v_mean_m_s_,reshapeMatrix); PlotW = reshape(MeanData.w_mean_m_s_,reshapeMatrix); PlotU(PlotU == -9999) = NaN; %Setting the non-values to NaN prevents them from being plotted in figures PlotV(PlotV == -9999) = NaN; PlotW(PlotW == -9999) = NaN; if ~any(~isnan(PlotV(:))) %Catching the 2D2C PIV Planes in XZ and XY VelMag = sqrt(PlotU.^2 + PlotW.^2); elseif ~any(~isnan(PlotW(:))) VelMag = sqrt(PlotU.^2 + PlotV.^2); else VelMag = sqrt(PlotU.^2 + PlotV.^2 + PlotW.^2); end %% Plotting normalised velocity if length(size(PlotX)) == 2 %to plot planar data figure('name','MeanFromMeanFile') surf(PlotX,PlotY,PlotZ,... VelMag,... 'EdgeColor','none',... 'LineStyle','none',... 'FaceColor','interp',... 'FaceLighting','none') xlabel('x (m)') ylabel('y (m)') zlabel('z (m)') colorbar('southoutside') set(gca,'DataAspectRatio',[1 1 1]) xlim([0.4 1.25]) ylim([-0.25 0.25]) zlim([0 0.4]) elseif length(size(PlotX)) == 3 %to plot tomographic data figure('name','MeanFromMeanFile') hold on %plot approximately mid vertical and mid horizontal planes mid_vert_pos = find(abs(UniqueY) == min(abs(UniqueY))); mid_horz_pos = find(abs(UniqueZ - 0.1945) == min(abs(UniqueZ - 0.1945))); surf(squeeze(PlotX(:,mid_vert_pos,:)),... squeeze(PlotY(:,mid_vert_pos,:)),... squeeze(PlotZ(:,mid_vert_pos,:)),... squeeze(VelMag(:,mid_vert_pos,:)),... 'EdgeColor','none',... 'LineStyle','none',... 'FaceColor','interp',... 'FaceLighting','none') surf(squeeze(PlotX(:,:,mid_horz_pos)),... squeeze(PlotY(:,:,mid_horz_pos)),... squeeze(PlotZ(:,:,mid_horz_pos)),... squeeze(VelMag(:,:,mid_horz_pos)),... 'EdgeColor','none',... 'LineStyle','none',... 'FaceColor','interp',... 'FaceLighting','none') view([-37.5 30]) xlabel('x (m)') ylabel('y (m)') zlabel('z (m)') colorbar('southoutside') set(gca,'DataAspectRatio',[1 1 1]) xlim([0.4 1.25]) ylim([-0.25 0.25]) zlim([0 0.4]) end catch disp('No mean data found') end %%Instantaneous try %% Load in the files InstPath = [Files(FileLocation & InstLocation).folder filesep Files(FileLocation & InstLocation).name]; InstData = readtable(InstPath); %% Arrange the data for plotting UniqueX = unique(InstData.x_m_); UniqueY = unique(InstData.y_m_); UniqueZ = unique(InstData.z_m_); XSamples = length(UniqueX); YSamples = length(UniqueY); ZSamples = length(UniqueZ); [~,TSamples] = size(InstData); TSamples = (TSamples-3)/3; %-3 for the xyz, /3 as data includes u, v and w components if any([XSamples YSamples ZSamples] == 1) if ZSamples == 1 %This section is required for the Z=0.194m plane as the data %was taken in tunnel coordinates and post processed into model %coordinates. % %As a result, the information on the matrix size has been %'lost' when writing a dataset, so the original matrix size has %been hard coded into this plotting code. It is understood this %isn't good practice. reshapeMatrix = [247 355 TSamples]; else %%Sorting the InstData as MATLAB is picky for plotting InstData = sortrows(InstData,[3 2 1]); reshapeMatrix = [XSamples YSamples ZSamples TSamples]; reshapeMatrix(reshapeMatrix == 1) = []; %remove singleton dimensions, saves doing squeeze() later. end else reshapeMatrix = [XSamples YSamples ZSamples TSamples]; %%Sorting the InstData as MATLAB is picky for plotting InstData = sortrows(InstData,[3 2 1]); end PlotX = reshape(InstData.x_m_,reshapeMatrix(1:end-1)); %don't need the time portion of reshapeMatrix PlotY = reshape(InstData.y_m_,reshapeMatrix(1:end-1)); PlotZ = reshape(InstData.z_m_,reshapeMatrix(1:end-1)); PlotU = reshape(table2array(InstData(:,[4:3:3+TSamples*3])),reshapeMatrix); %Reshape the variables for plotting PlotV = reshape(table2array(InstData(:,[5:3:3+TSamples*3])),reshapeMatrix); PlotW = reshape(table2array(InstData(:,[6:3:3+TSamples*3])),reshapeMatrix); ndims = length(size(PlotU)); PlotU(PlotU == -9999) = NaN; %Setting the non-values to NaN prevents them from being plotted in figures PlotV(PlotV == -9999) = NaN; PlotW(PlotW == -9999) = NaN; NaNPos = sum(isnan(PlotU),ndims); %U is present on all planes and the vector will be NaN for UVW if it is NaN for either U, V or W. NaNPos = NaNPos > TSamples*(1-NaNMax); PlotU = mean(PlotU,ndims,'omitnan'); %Take a mean of the variables along the 3rd (time) dimension PlotV = mean(PlotV,ndims,'omitnan'); PlotW = mean(PlotW,ndims,'omitnan'); %Removing those with bad NaN data PlotU(NaNPos) = NaN; PlotV(NaNPos) = NaN; PlotW(NaNPos) = NaN; if ~any(~isnan(PlotV(:))) %Catching the 2D2C PIV Planes in XZ and XY VelMag = sqrt(PlotU.^2 + PlotW.^2); elseif ~any(~isnan(PlotW(:))) VelMag = sqrt(PlotU.^2 + PlotV.^2); else VelMag = sqrt(PlotU.^2 + PlotV.^2 + PlotW.^2); end %% Plotting mean from instantaneous velocity if length(size(PlotX)) == 2 %to plot planar data figure('name','MeanFromInstantaneousFile') surf(PlotX,PlotY,PlotZ,... VelMag,... 'EdgeColor','none',... 'LineStyle','none',... 'FaceColor','interp',... 'FaceLighting','none') xlabel('x (m)') ylabel('y (m)') zlabel('z (m)') colorbar('southoutside') set(gca,'DataAspectRatio',[1 1 1]) xlim([0.4 1.25]) ylim([-0.25 0.25]) zlim([0 0.4]) elseif length(size(PlotX)) == 3 %to plot tomographic data figure('name','MeanFromInstantaneousFile') hold on %plot approximately mid vertical and mid horizontal planes mid_vert_pos = find(abs(UniqueY) == min(abs(UniqueY))); mid_horz_pos = find(abs(UniqueZ - 0.1945) == min(abs(UniqueZ - 0.1945))); surf(squeeze(PlotX(:,mid_vert_pos,:)),... squeeze(PlotY(:,mid_vert_pos,:)),... squeeze(PlotZ(:,mid_vert_pos,:)),... squeeze(VelMag(:,mid_vert_pos,:)),... 'EdgeColor','none',... 'LineStyle','none',... 'FaceColor','interp',... 'FaceLighting','none') surf(squeeze(PlotX(:,:,mid_horz_pos)),... squeeze(PlotY(:,:,mid_horz_pos)),... squeeze(PlotZ(:,:,mid_horz_pos)),... squeeze(VelMag(:,:,mid_horz_pos)),... 'EdgeColor','none',... 'LineStyle','none',... 'FaceColor','interp',... 'FaceLighting','none') view([-37.5 30]) xlabel('x (m)') ylabel('y (m)') zlabel('z (m)') colorbar('southoutside') set(gca,'DataAspectRatio',[1 1 1]) xlim([0.4 1.25]) ylim([-0.25 0.25]) zlim([0 0.4]) end catch disp('No instantaneous data found') end