tributarymatrix.m

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

                    % first purlin: all load transferred to attachment point at s_max:
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a2 = ds.*dy.*(dY-y_p-.5*dy)/dY;
                    a1 = zeros(size(a2));
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a2_1 = ds.*dy.*(y_p+.5*dy)/dY;
                    a1_1 = zeros(size(a2_1));
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j+1),tap_ind_i) = A_mx(attach_ind_i(j+1),tap_ind_i) + a2;
                elseif j == length(attach_ind_i)
                    % last purlin: all load transferred to attachment point at s_min:
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a1 = ds.*dy.*(dY-y_p-.5*dy)/dY;
                    a2 = zeros(size(a1));
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a1_1 = ds.*dy.*(y_p+.5*dy)/dY;
                    a2_1 = zeros(size(a1_1));
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j),tap_ind_i) = A_mx(attach_ind_i(j),tap_ind_i) + a1;
                else
                    % tributary area vectors for 2 attachment points at w_min (1) and w_max (2):
                    a1 = ds.*dy.*(dS-s_p-.5*ds).*(dY-y_p-.5*dy)/(dS*dY);
                    a2 = ds.*dy.*(s_p+.5*ds).*(dY-y_p-.5*dy)/(dS*dY);
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a1_1 = ds.*dy.*(dS-s_p-.5*ds).*(y_p+.5*dy)/(dS*dY);
                    a2_1 = ds.*dy.*(s_p+.5*ds).*(y_p+.5*dy)/(dS*dY);
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j),tap_ind_i) = A_mx(attach_ind_i(j),tap_ind_i) + a1;
                    A_mx(attach_ind_i(j+1),tap_ind_i) = A_mx(attach_ind_i(j+1),tap_ind_i) + a2;
                end

                if plot_on==2
                    ind_nz = find(a1 | a2);
                    for k = ind_nz
                        figure(h3);
                        clf;
                        plot(lc_i(:,3), lc_i(:,4), '+'); % plot all taps on current face
                        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:'); % plot frames
                        plot([s_min s_min], [y_f(1) y_f(2)],'k-',...
                            [s_max s_max], [y_f(1) y_f(2)],'k-'); % plot girts/purlins
                        plot([s_min s_max],[y_f(2) y_f(2)],'kx'); % mark attachment points to current frame
                        plot([s_min s_max],[y_f(1) y_f(1)],'kx'); % mark attachment points to previous frame
                        plot(lc_i(k,3), lc_i(k,4),'o'); % circle current tap
                        % plot bounds of tributary area for current tap on current span:
                        plot(s_min+s_p(k)+[0 ds(k) ds(k) 0 0],y_f(2)-y_p(k)-[0 0 dy(k) dy(k) 0],'-');
                        % label total tributary area of current tap:
                        text(lc_i(k,3), lc_i(k,4),...
                            ['Area = ' num2str(ds(k)*dy(k)) ' (sum = ' num2str(a1(k)+a2(k)+a1_1(k)+a2_1(k)) ')']); 
                        % label portion of tributary area applied to each attachment point:
                        text(s_min,y_f(2),['Area = ' num2str(a1(k))]);
                        text(s_max,y_f(2),['Area = ' num2str(a2(k))]);
                        text(s_min,y_f(1),['Area = ' num2str(a1_1(k))]);
                        text(s_max,y_f(1),['Area = ' num2str(a2_1(k))]);
                        title(['Application of tributary area for tap # ' num2str(lc_i(k,1)) ' on face ' num2str(i)]);
                        hold off;
                        pause;
                    end
                end

            end
            if y_f(3)~=y_f(2)
                % treat contributions from span B, beyond frame of interest:
                dY = y_f(3)-y_f(2);
                tap_bnds_ij_B = tap_bnds_ij;
                tap_bnds_ij_B(find(tap_bnds_ij_B(:,5)<y_f(2)),5) = y_f(2);
                tap_bnds_ij_B(find(tap_bnds_ij_B(:,6)<y_f(2)),6) = y_f(2);
                s_p = tap_bnds_ij_B(:,3)'-s_min;
                ds = tap_bnds_ij_B(:,4)'-tap_bnds_ij_B(:,3)';
                y_p = tap_bnds_ij_B(:,5)'-y_f(2);
                dy = tap_bnds_ij_B(:,6)'-tap_bnds_ij_B(:,5)';

                if j == 0
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a2 = ds.*dy.*(dY-y_p-.5*dy)/dY;
                    a1 = zeros(size(a2));
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a2_3 = ds.*dy.*(y_p+.5*dy)/dY;
                    a1_3 = zeros(size(a2_3));
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j+1),tap_ind_i) = A_mx(attach_ind_i(j+1),tap_ind_i) + a2;
                elseif j == length(attach_ind_i)
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a1 = ds.*dy.*(dY-y_p-.5*dy)/dY;
                    a2 = zeros(size(a1));
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a1_3 = ds.*dy.*(y_p+.5*dy)/dY;
                    a2_3 = zeros(size(a1_3));
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j),tap_ind_i) = A_mx(attach_ind_i(j),tap_ind_i) + a1;
                else
                    % tributary area vectors for 2 attachment points at s_min (1) and s_max (2):
                    a1 = ds.*dy.*(dS-s_p-.5*ds).*(dY-y_p-.5*dy)/(dS*dY);
                    a2 = ds.*dy.*(s_p+.5*ds).*(dY-y_p-.5*dy)/(dS*dY);
                    % tributary area vectors for attachment points on previous frame (for checking):
                    a1_3 = ds.*dy.*(dS-s_p-.5*ds).*(y_p+.5*dy)/(dS*dY);
                    a2_3 = ds.*dy.*(s_p+.5*ds).*(y_p+.5*dy)/(dS*dY);
                    % add contributions to tributary area matrix:
                    A_mx(attach_ind_i(j),tap_ind_i) = A_mx(attach_ind_i(j),tap_ind_i) + a1;
                    A_mx(attach_ind_i(j+1),tap_ind_i) = A_mx(attach_ind_i(j+1),tap_ind_i) + a2;
                end

                if plot_on==2
                    ind_nz = find(a1 | a2);
                    for k = ind_nz
                        figure(h3);
                        clf;
                        plot(lc_i(:,3), lc_i(:,4), '+'); % plot all taps on current face
                        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:'); % plot frames
                        plot([s_min s_min], [y_f(2) y_f(3)],'k-',...
                            [s_max s_max], [y_f(2) y_f(3)],'k-'); % plot girts/purlins
                        plot([s_min s_max],[y_f(2) y_f(2)],'kx'); % mark attachment points to current frame
                        plot([s_min s_max],[y_f(3) y_f(3)],'kx'); % mark attachment points to next frame
                        plot(lc_i(k,3), lc_i(k,4),'o'); % circle current tap
                        % plot bounds of tributary area for current tap on current span:
                        plot(s_min+s_p(k)+[0 ds(k) ds(k) 0 0],y_f(2)+y_p(k)+[0 0 dy(k) dy(k) 0],'-');
                        % label total tributary area of current tap:
                        text(lc_i(k,3), lc_i(k,4),...
                            ['Area = ' num2str(ds(k)*dy(k)) ' (sum = ' num2str(a1(k)+a2(k)+a1_3(k)+a2_3(k)) ')']); 
                        % label portion of tributary area applied to each attachment point:
                        text(s_min,y_f(2),['Area = ' num2str(a1(k))]);
                        text(s_max,y_f(2),['Area = ' num2str(a2(k))]);
                        text(s_min,y_f(3),['Area = ' num2str(a1_3(k))]);
                        text(s_max,y_f(3),['Area = ' num2str(a2_3(k))]);
                        title(['Application of tributary area for tap # ' num2str(lc_i(k,1)) ' on face ' num2str(i)]);
                        hold off;
                        pause;
                    end
                end

            end

        end

    end

end



function tap_trib_bnds_i = taptribbnds( tap_lc_i, lc_bnds_i )

    tap_row_y = unique(tap_lc_i(:,4));
    tap_col_w = unique(tap_lc_i(:,3));
    n_tap_row_y = zeros(size(tap_row_y));
    n_tap_col_w = zeros(size(tap_col_w));
    for I = 1:length(tap_row_y)
        n_tap_row_y(I) = length(find(tap_lc_i(:,4) == tap_row_y(I))); % number of taps per row
    end
    for I = 1:length(tap_col_w)
        n_tap_col_w(I) = length(find(tap_lc_i(:,3) == tap_col_w(I))); % number of taps per column
    end
    d_n_row = diff(n_tap_row_y); % change in number of taps per row
    d_n_col = diff(n_tap_col_w); % change in number of taps per column
    
    % if tap array is a simple grid: constant number of taps per row and per column
    if max(abs(d_n_row))==0 & max(abs(d_n_col))==0 

        tap_trib_bnds_i = simplegrid( tap_lc_i, lc_bnds_i );
    
    % if the number of taps per row changes only once on the face:
    elseif length(find(d_n_row))==1
        % slice the current face horizontally to yeild two simple grids:
        y_slice = mean([tap_row_y(find(d_n_row)) tap_row_y(find(d_n_row)+1)]);
        ind_A = find(tap_lc_i(:,4)<y_slice);
        ind_B = find(tap_lc_i(:,4)>y_slice);
        tap_lc_i_A = tap_lc_i(ind_A,:);
        tap_lc_i_B = tap_lc_i(ind_B,:);
        lc_bnds_i_A = lc_bnds_i; lc_bnds_i_A(4) = y_slice;
        lc_bnds_i_B = lc_bnds_i; lc_bnds_i_B(3) = y_slice;

        tap_trib_bnds_i_A = simplegrid( tap_lc_i_A, lc_bnds_i_A );
        tap_trib_bnds_i_B = simplegrid( tap_lc_i_B, lc_bnds_i_B );
        
        tap_trib_bnds_i = zeros(size(tap_lc_i,1),6);
        tap_trib_bnds_i(ind_A, :) = tap_trib_bnds_i_A;
        tap_trib_bnds_i(ind_B, :) = tap_trib_bnds_i_B;

    else
        error('Script can currently handle only one change in the number of taps per row on any face.')
    end

end

function tap_trib_bnds_i = simplegrid( tap_lc_i, lc_bnds_i )

    tap_trib_bnds_i = zeros(size(tap_lc_i,1),6);
    tap_trib_bnds_i(:,1:2) = tap_lc_i(:,1:2);
    tap_row_y = unique(tap_lc_i(:,4));
    tap_col_w = unique(tap_lc_i(:,3));

    % assign the left edge of the current face as the lower bound for taps in the first column:
    tap_trib_bnds_i(find(tap_lc_i(:,3)==tap_col_w(1)),3) = lc_bnds_i(1);
    % assign the right edge of the current face as the upper bound for taps in the last column:
    tap_trib_bnds_i(find(tap_lc_i(:,3)==tap_col_w(end)),4) = lc_bnds_i(2);
    % loop through intermediate tap rows and assign the midpoints as upper and lower bounds:
    for J = 1:length(tap_col_w)-1
        tap_trib_bnds_i(find(tap_lc_i(:,3)==tap_col_w(J)),4) = tap_col_w(J)+0.5*(tap_col_w(J+1)-tap_col_w(J));
        tap_trib_bnds_i(find(tap_lc_i(:,3)==tap_col_w(J+1)),3) = tap_col_w(J)+0.5*(tap_col_w(J+1)-tap_col_w(J));
    end
    % assign the lower edge of the current face as the lower bound for taps in the first row:
    tap_trib_bnds_i(find(tap_lc_i(:,4)==tap_row_y(1)),5) = lc_bnds_i(3);
    % assign the upper edge of the current face as the uper bound for taps in the last row:
    tap_trib_bnds_i(find(tap_lc_i(:,4)==tap_row_y(end)),6) = lc_bnds_i(4);
    % loop through intermediate tap rows and assign the midpoints as upper and lower bounds:
    for J = 1:length(tap_row_y)-1
        tap_trib_bnds_i(find(tap_lc_i(:,4)==tap_row_y(J)),6) = tap_row_y(J)+0.5*(tap_row_y(J+1)-tap_row_y(J));
        tap_trib_bnds_i(find(tap_lc_i(:,4)==tap_row_y(J+1)),5) = tap_row_y(J)+0.5*(tap_row_y(J+1)-tap_row_y(J));
    end
    
end


end