clear all;clc;close all;
%标准国际单位制
Length=2;%2m
min_H=0.2;max_H=0.6;%截面高度
Width=1;%截面宽度
E=3e10;%Pa 3e4MPa
u=0.2;%  泊松比
q0=1000;%均布力；1kN/m
F_jizhong=2e3;%集中力2/kN

%定义节点和单元
Ne_x=16;%单元列数4 
Ne_y=8;%单元行数2
D=[E/(1-u^2),u*E/(1-u^2),0;
       u*E/(1-u^2),E/(1-u^2),0;
       0,0,E/(2*(1+u))];
       

Ne=Ne_x*Ne_y;          % 单元总数 
Nnod=(Ne_x+1)*(Ne_y+1);;    % 节点总数 
%节点坐标 
x=linspace(-1,1,Ne_x+1); %节点从左到右的横坐标
rt=[repmat(Ne_y+1,1,Ne_x+1)]; % 每列节点数 5 3 5 3 5 3  ……
x=repelem(x,rt);       % x坐标扩展为所有节点的坐标
y=[];
dy=(max_H-min_H)/Ne_x/4;
for i=1:Ne_x
    y=[y,linspace(-max_H/2+(i-1)*2*dy,max_H/2-(i-1)*2*dy,Ne_y+1)];
end
y=[y,linspace(-min_H/2,min_H/2,Ne_y+1)];
Nodes=[x;y]';

% 单元节点编号
num=0;
for i=1:Ne_x
    for j=1:Ne_y
        num=num+1;
        Elements(num,:)=[(i-1)*(Ne_y+1)+j,(i-1)*(Ne_y+1)+j+1,i*(Ne_y+1)+j+1,i*(Ne_y+1)+j];
    end
end
%定义分布荷载
E_ID_p1 = Ne_y:Ne_y:Ne_x*Ne_y;      % 分布荷载作用单元编号

N_ID_p1 = zeros(size(E_ID_p1,2),size(Elements,2));%分布荷载所在单元节点编号
j = 1; 
for i_q = E_ID_p1
 N_ID_p1(j,:) = Elements(i_q,:);              %分布荷载所在单元节点编号
 j = j+1;
end

%施加集中荷载的节点编号
E_ID_p2 = max(Ne_y:Ne_y:Ne);            % 集中荷载作用单元编号***
Nod_1=Elements(E_ID_p2,4);                % 集中荷载节点编号***

%构件跨中节点编号
Elt_p3=max(Ne_y:Ne_y:Ne/2);            % 构件跨中单元编号
Nod_mid=Elements(Elt_p3,3) ;               % 构件跨中节点编号

%构件左端节点编号
Elt_p4=Ne_y           ;              % 构件左端单元编号
Nod_left=Elements(Elt_p4,2)    ;            % 构件左端节点编号

%边界条件：type = 1 -- fixed disp [ 1--节点固定; 0--节点自由 ]; 
%          type = 2 -- point load 
% node, x, y, z, rx, ry, rz, type 

Boundary =  [1 1.0 1.0 0.0 0.0 0.0 0.0 1; 
Nod_left 1.0 1.0 0.0 0.0 0.0 0.0 1; 
Nod_1 0.0 F_jizhong 0.0 0.0 0.0 0.0 2 ];

% 初始化结果存储空间
K = zeros( 2 * Nnod , 2 * Nnod );
F = zeros( 2 * Nnod , 1 );

for i=1:Ne
    ke=ele_stiff_matrix(i,Nodes,Elements,E,u,Width);
    ENodes( : , 1 ) = Elements( i , : ); 
    %组装为整体刚度阵
    iEnum = size( ENodes , 1 );      % 单元的结点总数 
    iDOFnum = 2 * iEnum;              % 单元的总自由度数 
    iDof = zeros( 2 * iEnum , 1 );   % 各结点自由度编号 

    for il=1:iEnum 
      inod = ENodes( il ); 
      iDof( 2 * il-1 : 2 * il) = 2 * inod-1 : 2 * inod; 
    end
    % 拼装：ir为行局部编号, irg为对应整体编号 
    for ir = 1 : iDOFnum 
      irg = iDof( ir ); 
      % ic为列局部编号, icg为对应整体编号 
      for ic = 1 : iDOFnum 
        icg = iDof( ic ); 
        K( irg , icg ) = K( irg , icg ) + ke( ir , ic ); 
      end 
    end
end
%施加荷载边界条件
%均布荷载
iLenElt = size(N_ID_p1,1);          %均布荷载所在单元编号
for ie = 1:iLenElt
  for ind = 1:2                   % ind结点自由度编号
    iNode = N_ID_p1(ie,:);           %获取单元结点
    irow = (iNode-1)*2+ind;
    Pe = UniLoad(ie,N_ID_p1,q0,Nodes,Elements,Width);    % 3级子程序：EquivalentLoad(ie)
    F(irow) = F(irow)+Pe(ind:2:end,:);
  end
end

% 集中荷载
iBct = size( Boundary , 1 ); 
for ib = 1 : iBct 
  iNode = Boundary( ib , 1 ); 
  iType = Boundary( ib , 8 ); 
  if iType == 2                       % 点荷载边界条件 
    for ind = 1 : 2                   
      irow = ( iNode-1 ) * 2 + ind; 
      F( irow ) = F( irow ) + Boundary( ib , ind + 1 ); 
    end 
  end 
end 
%施加支座边界条件
for ib = 1 : iBct
    iNode = Boundary( ib , 1 );       % 结点编号
    iType = Boundary( ib , 8 );
    if iType == 1             % 固定点边界条件
        for ind = 1 : 2
            irow = ( iNode-1 ) * 2 + ind;
            if ( Boundary( ib , ind + 1 ) == 1 )
                K( irow , : ) = 0.0;
                K( : , irow ) = 0.0;
                K( irow , irow ) = 1.0;       % 划零置一法
                F(irow)=0.0;
                %         gKA( irow , irow ) = 1.0e+20;  % 乘大数法
            end
        end
    end
end
Disp = K \ F;
% 存储结点位移
gNTu = zeros( Nnod , 2 );
for id = 1 : Nnod
    gNTu( id , 1 : 2 ) = Disp( 2 * id - 1 : 2 * id );
end
% 绘制变形图
SF=1000;
PlotResults2(Nodes,Elements,gNTu,SF);
fprintf('自由端最大挠度%fmm\n',max(abs(gNTu(Nnod,2)))*1000);
ElementNodeCount=4;
[NodeStrain,NodeStress]=CalculateStrainAndStress2(Disp,D,Nodes,Elements);
 PlotContour(Nodes,Elements,Disp,NodeStress(1,:))