tributarymatrix.m

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

function Astruct = tributarymatrix( hdf_input, y_frame, attach_pts, varargin )

% [ A,  ind_ts ] = tributarymatrix( hdf_filename , y_frame, attach_pts, varargin ) 
% [ A,  ind_ts ] = tributarymatrix( hdf_struct , y_frame, attach_pts, varargin ) 
% [ A,  ind_ts ] = tributarymatrix( [] , y_frame, attach_pts, varargin ) --> prompts user to select HDF file

plot_on = 0; % (0 ==> no plots, 1 ==> some plots, 2 ==> all plots) default can be overridden by optional user input

n_frames = size(y_frame,1);

if ischar( hdf_input )
    hdf_filename = hdf_input;
    hdf_struct = load_hdf_coords( hdf_filename );
elseif isstruct( hdf_input )
    hdf_struct = hdf_input;
else
    [filename pathname]= uigetfile({'*.hdf','HDF-files (*.hdf)'},'Select file for plotting taps');
    hdf_full_filename = [pathname filename]
    hdf_struct = load_hdf_coords( hdf_full_filename );
    assignin('base','hdf_full_filename', hdf_full_filename); % make filename available in base workspace for convenience
end

% Extract model dimensions (full-scale) from HDF data structure:
L = hdf_struct.Building_Length_Ft;
W = hdf_struct.Building_Width_Ft;
H = hdf_struct.Building_Height_Ft; % eave height
R = hdf_struct.Roof_Slope_in_12*W/2; % slope is given as rise (in inches) per foot of horizontal distance

% Initialize dimensions of structure with dimensions of model (default):
L0 = L; W0 = W; H0 = H; R0 = R;

% Process any optional inputs: assign alternative dimensions of structure, activate plotting
if nargin>3
    if mod(nargin-3,2)~=0
        error('Optional input arguments must occur in pairs: ...''PropertyName'',PropertyValue... ')
    else
        for i = 1:2:nargin-4
            if ischar(varargin{i}) & max(strcmp(varargin{i},{'plot_on','L0','W0','H0','R0'}))
                if isnumeric(varargin{i+1})
                    eval([varargin{i} '=' num2str(varargin{i+1}) ';']);
                else
                    error(['Non-numeric value specified for optional input >> ' varargin{i} ' <<']);
                end
            else
                error(['The specified optional input >> ' num2str(varargin{i}) ' << is an invalid variable name.' ]);
            end
        end
    end
end

if plot_on, plot3Dtaps( hdf_struct ); h1 = gcf; pause; end;

tap_c0_3d = hdf_struct.Tap_Coordinates_3D; % tap #, face #, x-coord, y-coord, z-coord
 
% tol = 10^-4; % tolerance for rounding
% tap_c0_3d(find(abs(tap_c0_3d)<tol))=0;
% tap_c0_3d(find(abs(tap_c0_3d(:,3)-W/2)<tol),3)=W/2;
% tap_c0_3d(find(abs(tap_c0_3d(:,3)-W)<tol),3)=W;
% tap_c0_3d(find(abs(tap_c0_3d(:,4)-L)<tol),3)=L;


% Face numbering and coordinate system in UWO HDF files:
%   Face 1 is the leeward face for zero wind direction, Face 3 is windward
%   (z-axis is upwards, walls are folded out here for display)
%   Internal pressure taps are labeled Face 0
%                                          
%                                          |x-axis
%                                          |
%                __________________________|
%               |                          |
%               |            4             |
%          _____|__________________________|_____
%         |     |                          |     |
%         |     |            6             |     |
%  0 deg  |     |                          |     |
%   ----> |  3  |----------ridge-----------|  1  |
%         |     |                          |     |
%         |     |            5             |     |
%  _______|_____|__________________________|_____|
%  y-axis       |                          |Origin
%               |            2             |
%               |__________________________|
%
%   Shift origin to opposite corner (on windward face for zero wind dir),
%   eliminate internal taps and end walls, and renumber faces as follows: 
%                __________________________
%               |                          |
%               |            1             |
%               |__________________________|_________________
%               |Origin                    |          y-axis
%               |            2             |     
%     0 deg     |                          |     
%       ---->   |--- ridge: face = 0 ------|    
%               |                          |     
%               |            3             |     
%               |__________________________|
%               |                          |
%               |            4             |
%               |__________________________|
%               |
%               |
%               |x-axis

% shift origin to opposite corner of model:
tap_c_3d = tap_c0_3d;
tap_c_3d(:,3) = W-tap_c0_3d(:,3);
tap_c_3d(:,4) = L-tap_c0_3d(:,4);
% eliminate internal pressure taps (labeled face 0):
ext_ind = find(tap_c_3d(:,2)~=0); 
tap_c_3d = tap_c_3d(ext_ind,:);
% eliminate taps on end walls:
ind_no_ends = intersect(find(tap_c_3d(:,4)~=0), find(tap_c_3d(:,4)~=L));
tap_c_3d = tap_c_3d(ind_no_ends,:);
% assign new face numbers:
tap_c_3d(find(tap_c_3d(:,3)==0),2) = 1; % left wall
tap_c_3d(find(tap_c_3d(:,3)>0 & tap_c_3d(:,3)<W/2 ),2) = 2; % left roof slope
tap_c_3d(find(tap_c_3d(:,3)==W/2 ),2) = 0; % ridge line
tap_c_3d(find(tap_c_3d(:,3)>W/2 & tap_c_3d(:,3)<W),2) = 3; % right roof slope
tap_c_3d(find(tap_c_3d(:,3)==W),2) = 4; % right wall

if plot_on
    color_list = {'g','k','b','r','c'};
    figure(h1);
    for i = 0:4
        ind_face = find(tap_c_3d(:,2)==i);
        if length(ind_face)
            plot3(tap_c_3d(ind_face,3),tap_c_3d(ind_face,4),tap_c_3d(ind_face,5),[color_list{i+1} '.'])
            hold on;
            x_face = mean(tap_c_3d(ind_face,3));
            y_face = mean(tap_c_3d(ind_face,4));
            z_face = mean(tap_c_3d(ind_face,5));
            text(x_face, y_face, z_face, ['Face ' num2str(i)],'Color',color_list{i+1});
        end
    end
    xlabel('x-axis'); ylabel('y-axis'); zlabel('z-axis');
    title('Tap coordinates with origin shifted to opposite corner, end walls removed, and faces renumbered');
    axis equal
    axis vis3d
    hold off;
    pause;
end

% Scale tap coordinates to match dimensions of structure of interest:
scale = 0;
if H0~=H | W0~=W | L0~=L | R0~=R
    scale =1;
    % scale x-coordinates:
    tap_c_3d(:,3) = tap_c_3d(:,3)*W0/W;
    % scale y-coordinates:
    tap_c_3d(:,4) = tap_c_3d(:,4)*L0/L;
    % scale z-coordinates:
    roof_ind = find(tap_c_3d(:,5)>H);
    wall_ind = find(tap_c_3d(:,5)<=H);
    tap_c_3d(roof_ind,5) = H0+(tap_c_3d(roof_ind,5)-H)*R0/R;
    tap_c_3d(wall_ind,5) = tap_c_3d(wall_ind,5)*H0/H;
end

if plot_on && scale
    color_list = {'g','k','b','r','c'};
    figure(h1);
    for i = 0:4
        ind_face = find(tap_c_3d(:,2)==i);
        if length(ind_face)
            plot3(tap_c_3d(ind_face,3),tap_c_3d(ind_face,4),tap_c_3d(ind_face,5),[color_list{i+1} '.'])
            hold on;
            x_face = mean(tap_c_3d(ind_face,3));
            y_face = mean(tap_c_3d(ind_face,4));
            z_face = mean(tap_c_3d(ind_face,5));
            text(x_face, y_face, z_face, ['Face ' num2str(i)],'Color',color_list{i+1});
        end
    end
    xlabel('x-axis'); ylabel('y-axis'); zlabel('z-axis');
    title('Tap coordinates scaled to dimensions of structure of interest');
    axis equal
    axis vis3d
    hold off;
    pause;
end
    
% Define 2D local coordinates for each face:
% (y-axis coincides with global y-axis, w-axis explained below)
%               _____________________________________________
%               |Origin                    |          y-axis
%     0 deg     |                          |     
%       ---->   |                          |    
%               |__________________________|
%               |
%               |
%               |w-axis
%               
% On the side walls (faces 1 & 4) the local w-axis is in the vertical plane;
% On the roof (faces 2 and 3), the local w-axis is in the horizontal plane.
% The w-coordinate increases upwards on face 1 and downwards on face 4;
% The w-coordinate increases from left to right on faces 2 and 3.
%                                        .
%                                   .         .
%                              .    2           3  .
%                         .---------->   .------->      .
%                                                       |
%                         .                             |
%                       1 ^                             v 4
%                         |                             .
%                         |                             .
%                         .                             .
%          

tap_lc = zeros(size(tap_c_3d,1),4);
tap_lc(:,1:2) = tap_c_3d(:,1:2); % tap #, face #
% y-coordinates:
tap_lc(:,4) = tap_c_3d(:,4); 
% s-coordinates:
tap_lc(find(tap_lc(:,2)==1),3)= tap_c_3d(find(tap_lc(:,2)==1),5);
tap_lc(find(tap_lc(:,2)==4),3)= tap_c_3d(find(tap_lc(:,2)==4),5);
hyp0 = sqrt((W0/2)^2+R0^2); % hypotenuse of triangle formed by W/2 and R
tap_lc(find(tap_lc(:,2)==2),3)= tap_c_3d(find(tap_lc(:,2)==2),3)*2*hyp0/W0;
tap_lc(find(tap_lc(:,2)==3),3)= (W0-tap_c_3d(find(tap_lc(:,2)==3),3))*2*hyp0/W0;

% duplicate the taps along the ridge line to be included in both faces 2 and 3:
tap_lc(find(tap_lc(:,2)==0),3)= W0/2;
tap_ridge_dup = tap_lc(find(tap_lc(:,2)==0),:);
tap_ridge_dup(:,2) = 3;