tributarymatrix.m

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

tap_lc(find(tap_lc(:,2)==0),2)= 2;
tap_lc = [tap_lc; tap_ridge_dup];

% define bounds of each face in s-y plane:
lc_bnds = [0 H0 0 L0; 0 hyp0 0 L0; 0 hyp0 0 L0; 0 H0 0 L0];

% Loop through faces, evaluate bounds of tributary area for each tap:
tap_trib_bnds = zeros(size(tap_lc,1),6); % initialize matrix for storing results: 
                                         %  tap #, face #, s_min, s_max, y_min, y_max

if plot_on, h2 = figure; end;
for i = 1:4
    ind_i = find(tap_lc(:,2)==i);
    n_tap_i = length(ind_i);
    tap_lc_i = tap_lc(ind_i,:);
    lc_bnds_i = lc_bnds(i,:);

    tap_trib_bnds_i = taptribbnds( tap_lc_i, lc_bnds_i );
    
    if plot_on
        figure(h2);
        plot(tap_lc_i(:,3), tap_lc_i(:,4), '+');
        hold on;
        plot([lc_bnds_i(1:2) lc_bnds_i(2:-1:1) lc_bnds_i(1)],...
            [lc_bnds_i(3) lc_bnds_i(3) lc_bnds_i(4) lc_bnds_i(4) lc_bnds_i(3)],'-');
        hold off;
        title(['Pressure tap layout on face ' num2str(i)]);
        pause;
        for j = 1:n_tap_i
            figure(h2);
            hold on;
            plot(tap_lc_i(j,3), tap_lc_i(j,4),'o');
            plot([tap_trib_bnds_i(j,3:4) tap_trib_bnds_i(j,4:-1:3) tap_trib_bnds_i(j,3)],...
                [tap_trib_bnds_i(j,5) tap_trib_bnds_i(j,5) tap_trib_bnds_i(j,6) tap_trib_bnds_i(j,6) tap_trib_bnds_i(j,5)], '-');
            hold off;
        end
        title(['Pressure tap layout on face ' num2str(i) ' with tributary areas indicated']);
        pause;
    end
    tap_trib_bnds(ind_i,:) = tap_trib_bnds_i;
end

% Reflect local tap coordinates and tributary bounds about axes of symmetry:
% (b) Reflect about x = W0/2 (ridge axis) by reflecting face numbers:
tap_lc_b = tap_lc;
tap_lc_b(find(tap_trib_bnds(:,2)==1),2)=4;
tap_lc_b(find(tap_trib_bnds(:,2)==2),2)=3;
tap_lc_b(find(tap_trib_bnds(:,2)==3),2)=2;
tap_lc_b(find(tap_trib_bnds(:,2)==4),2)=1;
tap_trib_bnds_b = tap_trib_bnds; 
tap_trib_bnds_b(find(tap_trib_bnds(:,2)==1),2)=4;
tap_trib_bnds_b(find(tap_trib_bnds(:,2)==2),2)=3;
tap_trib_bnds_b(find(tap_trib_bnds(:,2)==3),2)=2;
tap_trib_bnds_b(find(tap_trib_bnds(:,2)==4),2)=1;
% (c) Reflect about y = L0/2 by reflecting y-coordinates:
tap_lc_c = tap_lc;
tap_lc_c(:,4) = L0-tap_lc(:,4);
tap_trib_bnds_c = tap_trib_bnds;
tap_trib_bnds_c(:,5)=L0-tap_trib_bnds(:,6); % y_min becomes y_max after reflection
tap_trib_bnds_c(:,6)=L0-tap_trib_bnds(:,5); % y_max becomes y_min after reflection
% (d) Reflect about both axis by reflecting face numbers and y-coordinates:
tap_lc_d = tap_lc;
tap_lc_d(find(tap_trib_bnds(:,2)==1),2)=4;
tap_lc_d(find(tap_trib_bnds(:,2)==2),2)=3;
tap_lc_d(find(tap_trib_bnds(:,2)==3),2)=2;
tap_lc_d(find(tap_trib_bnds(:,2)==4),2)=1;
tap_lc_d(:,4) = L0-tap_lc(:,4);
tap_trib_bnds_d = tap_trib_bnds;
tap_trib_bnds_d(find(tap_trib_bnds(:,2)==1),2)=4;
tap_trib_bnds_d(find(tap_trib_bnds(:,2)==2),2)=3;
tap_trib_bnds_d(find(tap_trib_bnds(:,2)==3),2)=2;
tap_trib_bnds_d(find(tap_trib_bnds(:,2)==4),2)=1;
tap_trib_bnds_d(:,5)=L0-tap_trib_bnds(:,6); % y_min becomes y_max after reflection
tap_trib_bnds_d(:,6)=L0-tap_trib_bnds(:,5); % y_max becomes y_min after reflection

% Select the taps whose tributary area overlaps the tributary area of the specified frames:
ind_trib = [];
for f = 1:n_frames
    ind_trib_ab = intersect(find(tap_trib_bnds(:,6)>y_frame(f,1)),find(tap_trib_bnds(:,5)<y_frame(f,3)));
    ind_trib_cd = intersect(find(tap_trib_bnds_c(:,6)>y_frame(f,1)),find(tap_trib_bnds_c(:,5)<y_frame(f,3)));
    ind_trib = union( ind_trib_ab, ind_trib_cd );
%end

    tap_lc = tap_lc(ind_trib,:);
    tap_lc_b = tap_lc_b(ind_trib,:);
    tap_lc_c = tap_lc_c(ind_trib,:);
    tap_lc_d = tap_lc_d(ind_trib,:);
    tap_trib_bnds = tap_trib_bnds(ind_trib,:);
    tap_trib_bnds_b = tap_trib_bnds_b(ind_trib,:);
    tap_trib_bnds_c = tap_trib_bnds_c(ind_trib,:);
    tap_trib_bnds_d = tap_trib_bnds_d(ind_trib,:);
    tap_num = tap_lc(:,1);
    [tf, tap_ind] = ismember(tap_num, hdf_struct.Tap_Position_List);
    Astruct(f).tap_ind = tap_ind;

    % error checking (can be removed eventually)
    % min(tap_lc(:,1)==tap_trib_bnds(:,1));
    if ismember(0, tf)
        error('Tap number not found in "Tap_Position_List" field in HDF file.');
    elseif ~isempty(setxor( tap_num, hdf_struct.Tap_Position_List(tap_ind)))
        error('Error in matching tap number with "Tap_Position_List" in HDF file.');
    end

    if plot_on
        % Find indices in tap_c_3d of taps that are tributary to current frame:
        [ tf2, ind_3d ] = ismember( tap_num, tap_c_3d(:,1) );
        color_list = {'b','r','g','c','m','y'};

        figure(h1);
        plot3([0 W W W/2 0 0],[0 0 0 0 0 0],[0 0 H H+R H 0],'k-', ...
        [0 W W W/2 0 0],[L L L L L L],[0 0 H H+R H 0],'k-', ...
        [0 0],[0 L],[0 0],'k-',[0 0],[0 L],[H H],'k-',...
        [W W],[0 L],[0 0],'k-',[W W],[0 L],[H H],'k-',...
        [W/2 W/2],[0 L],[H+R H+R],'k-'); hold on;
        plot3(tap_c_3d(:,3),tap_c_3d(:,4),tap_c_3d(:,5),'k.');
        plot3(tap_c_3d(ind_3d,3),tap_c_3d(ind_3d,4),tap_c_3d(ind_3d,5),'ro');
%    for i = 1:n_frames
        f_clr = mod(f,6); f_clr(find(i_clr==0))=6;
        clr = color_list{f_clr};
        plot3([0 0 W0/2 W0 W0], y_frame(f,1)*ones(1,5),[0 H0 H0+R0 H0 0],'k:','LineWidth',1);
        plot3([0 0 W0/2 W0 W0], y_frame(f,2)*ones(1,5),[0 H0 H0+R0 H0 0],[clr '-'],'LineWidth',2);
        plot3([0 0 W0/2 W0 W0], y_frame(f,3)*ones(1,5),[0 H0 H0+R0 H0 0],'k:','LineWidth',1);
        plot3([0 0 W0/2 W0 W0], (L0-y_frame(f,1))*ones(1,5),[0 H0 H0+R0 H0 0],'k:','LineWidth',1);
        plot3([0 0 W0/2 W0 W0], (L0-y_frame(f,2))*ones(1,5),[0 H0 H0+R0 H0 0],[clr '-'],'LineWidth',2);
        plot3([0 0 W0/2 W0 W0], (L0-y_frame(f,3))*ones(1,5),[0 H0 H0+R0 H0 0],'k:','LineWidth',1);
%    end
        xlabel('x-axis'); ylabel('y-axis'); zlabel('z-axis');
        title('Tap coordinates plotted with specified frame and tributary taps circled in red.');
        axis equal
        axis vis3d
        rotate3d on;
        hold off;
        pause;
    end

%for f = 1:n_frames
    Astruct(f).Aa = trib_mx( y_frame(f,:), tap_trib_bnds, tap_lc );
    Astruct(f).Ab = trib_mx( y_frame(f,:), tap_trib_bnds_b, tap_lc_b );
    Astruct(f).Ac = trib_mx( y_frame(f,:), tap_trib_bnds_c, tap_lc_c );
    Astruct(f).Ad = trib_mx( y_frame(f,:), tap_trib_bnds_d, tap_lc_d );
end

function A_mx = trib_mx( y_f, tap_bnds, lc )

    % find indices of taps that overlap current frame:
    trib_ind = intersect(find(tap_bnds(:,6)>y_f(1)),find(tap_bnds(:,5)<y_f(3)));
    % Adjust tributary bounds to lie within the tributary area of the current frame:
    tap_bnds(find(tap_bnds(:,6)>y_f(3)),6) = y_f(3);
    tap_bnds(find(tap_bnds(:,5)<y_f(1)),5) = y_f(1);

    A_mx = zeros(size(attach_pts,1), length(tap_ind)); % initialize tributary area matrix
    if plot_on==2, h3 = figure; end;
    for i = 1:4
        % find taps on current face that overlap current frame:
        tap_ind_i = intersect( trib_ind, find(lc(:,2)==i));
        n_tap_i = length(tap_ind_i);
        lc_i = lc(tap_ind_i,:);
        lc_bnds_i = lc_bnds(i,:);
        tap_bnds_i = tap_bnds(tap_ind_i, :);

        if plot_on==2
            figure(h3);
            plot(lc_i(:,3), lc_i(:,4), '+');
            hold on;
            plot(lc_bnds_i(1:2),[y_f(1) y_f(1)],'k:',...
                lc_bnds_i(1:2),[y_f(2) y_f(2)],'k-',...
                lc_bnds_i(1:2),[y_f(3) y_f(3)],'k:');
            hold off;
            title(['Pressure taps on face ' num2str(i) ' that are tributary to specified frame.']);
            pause;
            for j = 1:n_tap_i
                figure(h3);
                hold on;
                plot(lc_i(j,3), lc_i(j,4),'o');
                plot([tap_bnds_i(j,3:4) tap_bnds_i(j,4:-1:3) tap_bnds_i(j,3)],...
                    [tap_bnds_i(j,5) tap_bnds_i(j,5) tap_bnds_i(j,6) tap_bnds_i(j,6) tap_bnds_i(j,5)], '-');
                hold off;
            end
            title(['Pressure tap layout on face ' num2str(i) ' with tributary areas indicated']);
            pause;
        end
        % find attachment points on current face:
        attach_ind_i = find(attach_pts(:,1)==i);
        attach_s_i = attach_pts(attach_ind_i,2); % s-coordinates

        for j = 0:length(attach_ind_i)

            if j == 0
                s_min = lc_bnds_i(1);
                s_max = attach_s_i(j+1);
            elseif j == length(attach_ind_i)
                s_min = attach_s_i(j);
                s_max = lc_bnds_i(2);
            else
                s_min = attach_s_i(j);
                s_max = attach_s_i(j+1);
            end
            dS = s_max - s_min;

            if dS<=0
                continue
            end
            % adjust tributary bounds to lie within current panel:
            tap_bnds_ij = tap_bnds_i;
            tap_bnds_ij(find(tap_bnds_ij(:,3)>s_max),3) = s_max;
            tap_bnds_ij(find(tap_bnds_ij(:,4)>s_max),4) = s_max;
            tap_bnds_ij(find(tap_bnds_ij(:,3)<s_min),3) = s_min;
            tap_bnds_ij(find(tap_bnds_ij(:,4)<s_min),4) = s_min;

            if y_f(1)~=y_f(2)
                % treat contributions from span A preceeding frame of interest:
                dY = y_f(2)-y_f(1);
                tap_bnds_ij_A = tap_bnds_ij;
                tap_bnds_ij_A(find(tap_bnds_ij_A(:,5)>y_f(2)),5) = y_f(2);
                tap_bnds_ij_A(find(tap_bnds_ij_A(:,6)>y_f(2)),6) = y_f(2);
                s_p = tap_bnds_ij_A(:,3)'-s_min; % "s prime" (lever arm from panel origin to corner of tributary area)
                ds = tap_bnds_ij_A(:,4)'-tap_bnds_ij_A(:,3)'; 
                y_p = y_f(2)-tap_bnds_ij_A(:,6)'; % "y prime" (lever arm from panel origin to corner of tributary area)
                dy = tap_bnds_ij_A(:,6)'-tap_bnds_ij_A(:,5)';

                if j == 0