computedif.m

低矮房屋风压系数、风荷载计算分析matlab程序

function DIF_struct = computeDIF( bldg_struct, HDF_path, out_opt_struct )

% Input arguments:

% bldg_struct           a structure array containing 
% HDF_path
%   
clr = get(0,'DefaultUicontrolBackgroundColor'); % Background color for figures and dialog boxes

% if input arguments are not specified, prompt user to select:
% Building input data:
if ~exist('bldg_struct') || isempty(bldg_struct)
    bldg_struct = readBLDGfile( );
    if isfield(bldg_struct,'error')
        DIF_struct.error = 1; return;
    end
end
% Location of HDF files:
if ~exist('HDF_path') || isempty(HDF_path)
    % Load pathname for HDF files used in last analysis, if available, otherwise create empty string:
    if exist('tmp_HDF_path.mat','file'), load('tmp_HDF_path'); end
    if ~exist('HDF_path','var') || ~isstr(HDF_path), HDF_path = ''; end
    % Prompt use to select pressure database file:
    [HDF_filename HDF_path] = uigetfile({'*.hdf','HDF files (*.hdf)'},'Select an HDF pressure database file in the directory of interest',HDF_path);
    if any(HDF_filename~=0)
        save('tmp_HDF_path','HDF_path');
    else
        err = errordlg('HDF pressure database file selection cancelled', 'Analysis cancelled');
        uiwait(err); DIF_struct.error = 1; return;
    end
end
% Output options:
field_names = {'ts_plot','obs_pks','est_pks','mult_pks','pk_loads'};
if ~exist('out_opt_struct','var') || ~isstruct(out_opt_struct) || ~isempty(setdiff(field_names,fieldnames(out_opt_struct)))
    % Open dialog box to prompt user to select:
    [out_opt_struct, ok] = output_options;
    if ok~=1
        err = errordlg('Output options selection cancelled', 'Analysis cancelled');
        uiwait(err); DIF_struct.error = 1; return;
    end
end


% Extract building model information from 'bldg_struct':
d0 = bldg_struct.d0;
W0 = d0(1); L0 = d0(2); H0 = d0(3); R0 = d0(4);
length_units = bldg_struct.length_units;
attach_pts = bldg_struct.attach_pts;
frame_coords = bldg_struct.frame_coords;
n_f = size(frame_coords,1); % number of frames
resp_names = bldg_struct.resp_names; % descriptive names of responses of interest
resp_units = bldg_struct.resp_units; % units of responses of interest
N = bldg_struct.inf_coeff; % influence coefficients
[n_a, n_r] = size(N); % number of attachment points, number of response quantities
% Frame labels, for time series plots:
for i = 1:n_f
    frame_label{i} = ['Frame ' num2str(i) ': y = ' num2str(frame_coords(i,2)) ' ' length_units ];
end
% Response labels, for time series plots:
for i = 1:n_r
    resp_label{i} = ['Response ' num2str(i) ': ' resp_names{i}];
end

% Obtain list of which files in the pressure database folder to process:
dir_listing = dir(HDF_path);
file_list = dir_listing(find(cat(1,dir_listing.isdir)==0)); % eliminate directories from list
dir_ind = 1; % initialize index for counting the number of wind directions
for i = 1:length(file_list)
    filename = file_list(i).name; 
    hdf_ext_ind = strfind( filename, '.HDF' );
    if ~isempty(hdf_ext_ind)
        hdf_filename_root(dir_ind,:) = filename(1:hdf_ext_ind-5);
        theta_a_str(dir_ind,:) = filename(hdf_ext_ind-4:hdf_ext_ind-1);
        theta_a(dir_ind,1) = str2num( theta_a_str(dir_ind,:) )/10;
        dir_ind = dir_ind+1;
    end
end
hdf_filename_root = unique(hdf_filename_root,'rows');
if size(hdf_filename_root,1)~=1
    errordlg(['All HDF file names in the selected directory must be identical '...
        'except for the wind direction, specified by the last 4 characters before the file extension.'], ...
        'Error in pressure database folder contents');
    DIF_struct.error = 1; return;
end

% Display list of wind directions available in selected directory, allow user to select directions to use in analysis:
theta_a_cell = mat2cell(theta_a_str(:,1:3),ones(size(theta_a)));
[theta_a_sel,ok] = listdlg('ListString',theta_a_cell, 'InitialValue', [1:length(theta_a_cell)], ...
    'PromptString', {'Select which wind directions to','use in the analysis from the','following list for which pressure','databases are available:'}, ...
    'Name','Wind direction selection.');
if ok
    theta_a = theta_a(theta_a_sel); theta_a_str = theta_a_str(theta_a_sel,:); p = length(theta_a);
else
    errordlg('Wind direction selection cancelled by user.','Analysis cancelled');
    DIF_struct.error = 1; return;
end

theta_a(find(theta_a==360))=0; % ensure that >=0 & <360
theta_b = 360 - theta_a;
theta_b(find(theta_b==360))=0; % ensure that >=0 & <360
theta_c = 180 - theta_a;
theta_c(find(theta_c<0))=theta_c(find(theta_c<0))+360; % make >=0 & <360
theta_d = 180 + theta_a;
theta_d(find(theta_d>=360))=theta_d(find(theta_d>=360))-360; % make >=0 & <360

% Load pressure database model information from the first file in the list:
HDF_struct = load_hdf_coords( fullfile( HDF_path, [hdf_filename_root theta_a_str(1,:) '.HDF'] ) );

% extract full-scale dimensions of building:
W = HDF_struct.Building_Width_Ft; L = HDF_struct.Building_Length_Ft; H = HDF_struct.Building_Height_Ft;
R = HDF_struct.Roof_Slope_in_12*W/2; 
d = [W L H R];
if HDF_struct.A_Tap_Coordinates_3D_Units ~= 'FS_Feet'
    errordlg(['Unrecognized length units in HDF pressure database file:  >> ' HDF_struct.A_Tap_Coordinates_3D_Units ' << ("FS_Feet" expected)'],...
        'Error reading HDF file');
    DIF_struct.error = 1; return;
elseif ~strcmp(bldg_struct.terrain,'Directional') && ~strcmp(HDF_struct.Exposure_Name, bldg_struct.terrain)
    errordlg(['The terrain conditions specified in the HDF file ("' HDF_struct.Exposure_Name '") do not match the terrain' ...
        ' specified in the building input file ("' bldg_struct.terrain '").'], 'Error: terrain mismatch');
    DIF_struct.error = 1; return;
end
DIF_struct.terrain = HDF_struct.Exposure_Name; % 'Open_Country' or 'Suburban'

% store quantities for output in a structure array:
DIF_struct.resp_names = resp_names;
DIF_struct.resp_units = bldg_struct.resp_units;
% output dimensions in same units as building input file
DIF_struct.length_units = bldg_struct.length_units; 
DIF_struct.d0 = bldg_struct.d0;
DIF_struct.d = d;
DIF_struct.frame_coords = bldg_struct.frame_coords;
DIF_struct.attach_pts = bldg_struct.attach_pts;
DIF_struct.force_units = bldg_struct.force_units;
DIF_struct.ws_units = bldg_struct.ws_units;


% tol = 10^-3;
% % Compare dimensions in HDF file with dimensions in building input file and prompt user to confirm scaling of tap coordinates:
% button = 'Yes';
% if abs(H-H0)/H>tol || abs(W-W0)/W>tol || abs(L-L0)/L>tol || abs(R-R0)/R>tol
%     button = questdlg(['The building dimensions in the pressure database file ' ...
%         '(W=' num2str(W) ', L=' num2str(L) ', H=' num2str(H) ',  R=' num2str(R) ') ' ...
%         'do not match the dimensions specified in the building input file ' ...
%         '(W=' num2str(W0) ', L=' num2str(L0) ', H=' num2str(H0) ',  R=' num2str(R0) ') ' ...
%         'Proceed with scaling the tap coordinates to match the dimensions ' ...
%         'in the building input file?'], 'Mismatch of building dimensions: scale tap coordinates?','Yes','Cancel','Yes');
% end
% if ~strcmp(button,'Yes')
%     error('Mismatch of building dimensions: analysis cancelled by user.');
% end

% Assemble tributary matrices for specified frames:
Astruct = tributarymatrix( HDF_struct, frame_coords, attach_pts, 'plot_on', 0, ...
    'H0', H0, 'W0', W0, 'L0', L0, 'R0', R0);

ro = 0.0023769; % Air density (lb*s^2/ft^4) at sea level, 15 degrees C

% Initialize arrays for storing peak response quantities:
X_max_obs_a = zeros(n_r,p,n_f);
X_max_obs_b = zeros(n_r,p,n_f);
X_max_obs_c = zeros(n_r,p,n_f);
X_max_obs_d = zeros(n_r,p,n_f);
X_min_obs_a = zeros(n_r,p,n_f);
X_min_obs_b = zeros(n_r,p,n_f);
X_min_obs_c = zeros(n_r,p,n_f);
X_min_obs_d = zeros(n_r,p,n_f);
if out_opt_struct.est_pks==1
    X_max_est_a = zeros(n_r,p,n_f);
    X_max_est_b = zeros(n_r,p,n_f);
    X_max_est_c = zeros(n_r,p,n_f);
    X_max_est_d = zeros(n_r,p,n_f);
    X_min_est_a = zeros(n_r,p,n_f);
    X_min_est_b = zeros(n_r,p,n_f);
    X_min_est_c = zeros(n_r,p,n_f);
    X_min_est_d = zeros(n_r,p,n_f);
end    
if out_opt_struct.pk_loads == 1
    load_max_a = zeros(n_a, n_r, p, n_f);
    load_max_b = zeros(n_a, n_r, p, n_f);
    load_max_c = zeros(n_a, n_r, p, n_f);
    load_max_d = zeros(n_a, n_r, p, n_f);
    load_min_a = zeros(n_a, n_r, p, n_f);
    load_min_b = zeros(n_a, n_r, p, n_f);
    load_min_c = zeros(n_a, n_r, p, n_f);
    load_min_d = zeros(n_a, n_r, p, n_f);
end

h1 = waitbar(0, ['Please wait... Analyzing HDF file ' num2str(1) ' of ' num2str(p)],'Color',clr,'WindowStyle','modal',...
    'CreateCancelBtn',@WaitBarCancel,'Name','Analysis in progress...','CloseRequestFcn',@WaitBarCancel); % cancel button

% loop through list of wind directions
for i = 1:p
    
    HDF_filename_i = [hdf_filename_root theta_a_str(i,:) '.HDF'];
    for f = 1:n_f
        drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
        waitbar((i-1)/p, h1, ['Please wait... Analyzing HDF file ' num2str(i) ' of ' num2str(p)]); drawnow; % update waitbar

    	Cp = cpmatrix( fullfile(HDF_path, HDF_filename_i) , Astruct(f).tap_ind ); % Pressure coefficient time series
        TS_struct(f).n_s = size(Cp,2); % number of samples
        drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
        
        % Compute time series of response quantities of interest using influence vectors:
        TS_struct(f).ts_a = 0.5*ro*N'*Astruct(f).Aa*Cp;
        TS_struct(f).ts_b = 0.5*ro*N'*Astruct(f).Ab*Cp;
        TS_struct(f).ts_c = 0.5*ro*N'*Astruct(f).Ac*Cp; 
        TS_struct(f).ts_d = 0.5*ro*N'*Astruct(f).Ad*Cp;
    
        % Compute peak value (largest or smallest) of each internal force quantity;
        % Identify time sample at which peak occurs for each quantity and 
        % store the corresponding load distribution:
        [X_max_obs_a(:,i,f), TS_struct(f).max_sample_a] = max(TS_struct(f).ts_a,[],2);
        [X_max_obs_b(:,i,f), TS_struct(f).max_sample_b] = max(TS_struct(f).ts_b,[],2);
        [X_max_obs_c(:,i,f), TS_struct(f).max_sample_c] = max(TS_struct(f).ts_c,[],2);
        [X_max_obs_d(:,i,f), TS_struct(f).max_sample_d] = max(TS_struct(f).ts_d,[],2);

        [X_min_obs_a(:,i,f), TS_struct(f).min_sample_a] = min(TS_struct(f).ts_a,[],2);
        [X_min_obs_b(:,i,f), TS_struct(f).min_sample_b] = min(TS_struct(f).ts_b,[],2);
        [X_min_obs_c(:,i,f), TS_struct(f).min_sample_c] = min(TS_struct(f).ts_c,[],2);
        [X_min_obs_d(:,i,f), TS_struct(f).min_sample_d] = min(TS_struct(f).ts_d,[],2);

        if out_opt_struct.pk_loads == 1
            load_max_a( :, :, i, f ) = 0.5*ro*Astruct(f).Aa*Cp(:, TS_struct(f).max_sample_a );
            load_max_b( :, :, i, f ) = 0.5*ro*Astruct(f).Ab*Cp(:, TS_struct(f).max_sample_b );
            load_max_c( :, :, i, f ) = 0.5*ro*Astruct(f).Ac*Cp(:, TS_struct(f).max_sample_c );
            load_max_d( :, :, i, f ) = 0.5*ro*Astruct(f).Ad*Cp(:, TS_struct(f).max_sample_d );
            load_min_a( :, :, i, f ) = 0.5*ro*Astruct(f).Ad*Cp(:, TS_struct(f).min_sample_a );
            load_min_b( :, :, i, f ) = 0.5*ro*Astruct(f).Ab*Cp(:, TS_struct(f).min_sample_b );
            load_min_c( :, :, i, f ) = 0.5*ro*Astruct(f).Ac*Cp(:, TS_struct(f).min_sample_c );
            load_min_d( :, :, i, f ) = 0.5*ro*Astruct(f).Ad*Cp(:, TS_struct(f).min_sample_d );
        end
        drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
        if out_opt_struct.est_pks==1
            [X_max_est_a(:,i,f), X_min_est_a(:,i,f)] = maxminest( TS_struct(f).ts_a );
            [X_max_est_b(:,i,f), X_min_est_b(:,i,f)] = maxminest( TS_struct(f).ts_b );
            drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
            [X_max_est_c(:,i,f), X_min_est_c(:,i,f)] = maxminest( TS_struct(f).ts_c );
            [X_max_est_d(:,i,f), X_min_est_d(:,i,f)] = maxminest( TS_struct(f).ts_d );
        end
    end % end loop through frames
    if out_opt_struct.ts_plot
        drawnow; if ~ishandle(h1), return; end % check if user pressed cancel
        set(h1,'Visible','off'); % hide waitbar
        h2 = TSplot;
        uiwait(h2);    
        if isfield(DIF_struct,'error'), return; end
        set(h1,'Visible','on');
    end % end if plot_on
end % end loop through wind directions

if ishandle(h1), delete(h1); end % close waitbar

theta_all = [theta_a; theta_b; theta_c; theta_d];
theta = unique(theta_all);
n_q = length(theta);
DIF_struct.theta_all = theta_all(:);
DIF_struct.theta = theta(:);
if out_opt_struct.obs_pks
    X_max_obs_all = cat(2,X_max_obs_a, X_max_obs_b, X_max_obs_c, X_max_obs_d);
    X_min_obs_all = cat(2,X_min_obs_a, X_min_obs_b, X_min_obs_c, X_min_obs_d);
    X_max_obs = zeros(n_r,n_q,n_f);
    X_min_obs = zeros(n_r,n_q,n_f);
    for i = 1:n_q
        X_max_obs(:,i,:)=mean(X_max_obs_all(:,find(theta_all==theta(i)),:),2);
        X_min_obs(:,i,:)=mean(X_min_obs_all(:,find(theta_all==theta(i)),:),2);