二维三角形元有限元的MATLAB实现

%% 主程序
%% 自定的各变量参数
%%
Lx = 2;
Ly = 2;
numnodx =10;
numnody =10;
numel = (numnodx-1)*(numnody-1)*2; %总的三角形单元数
X = linspace(0,Lx,numnodx);
Y = linspace(0,Ly,numnody);
Ke = sparse(numnodx*numnody,numnodx*numnody);
Fe = sparse(numnodx*numnody,1);
%% 边界条件
bpoint = [(1:1:numnodx),(numnodx*numnody-numnodx+1:1:numnodx*numnody)];   % 强制性边界点的编号，本例子中是两端
bvalue = [50.*ones(1,numnodx),100.*ones(1,numnodx)];  %强制性边界条件的取值，在第一个点的地方为u0,最后一个点为u1
%% 组装整体刚度矩阵和整体荷载向量
NOD = connect_m(numnodx, numnody);
for k = 1 : numel
    ke = Stiffness_matrix(k,numnodx, numnody, X, Y,NOD);
    fe = load_vector(k,numnodx, numnody, X, Y,NOD);
    sctr = NOD(:,k);
    Ke(sctr,sctr) = Ke(sctr,sctr) + ke;
    Fe(sctr) = Fe(sctr) + fe;
end 
 K0 = Ke;
%% 约束条件处理
%% 处理已知部分
for i = 1:length(bpoint)
   n = bpoint(i);
   Ke(n,:) = 0;
   Ke(n,n) = 1;
   Fe(n) = bvalue(i);
   
end
  u_coeff = Ke\Fe; 

%% 绘图，可视化
u_cal = u_coeff;
u_cal_re = reshape(u_coeff,numnody,numnodx);
u_cal_re = full(u_cal_re);
h = mesh(linspace(0,Lx,numnodx)',linspace(0,Ly,numnody)',u_cal_re);
title([' FE Solutions ','(',num2str(numnodx),'×',num2str(numnody),')']);%标题

%%  调用的函数 
%% 关联矩阵
function NOD = connect_m(numnodx, numnody) %nel可以省略，因为在二维有限元中三角形元的nel是3
  Snumnelx = numnodx - 1; %横向矩形形单元数
  Snumnely = numnody - 1; %纵向矩形单元数
  n = 1 : numnodx * numnody; %总结点数
  A = reshape(n, numnodx, numnody);  %节点同形状的矩阵
  NOD = zeros(3, 2*(Snumnelx * Snumnely));
  for k = 1:2*(Snumnelx * Snumnely)
  k0 = ceil(k/2); %向上取整，算出i在第几个矩形元
    xnum = rem(k0, Snumnelx ); %往左数第几个矩形单元
    if xnum == 0
        xnum = Snumnelx;%余为零时，即为最右边那一个
    end
    ynum = ceil(k0/Snumnelx ); %ceil向上取整得到单元i是从下往上数第几个
    a = A(xnum:xnum + 1,ynum:ynum + 1); %一个矩形元对应四个节点
    a1 = a(:);
    erem2 = rem(k,2);
    if erem2 == 1 %找出这个单元在矩形元中是哪种类型的三角元
        NOD0 = a1([3 1 4])';
    else
        NOD0 = a1([2 4 1])';
    end
    NOD(:,k) = [NOD0(1),NOD0(2),NOD0(3)];
  end
end

%% 节点坐标
function pi = Node_coordinate(i,numnodx,numnody,X,Y)
    xnum = rem(i, numnodx ); %往左数第几个点
    if xnum == 0
        xnum = numnodx; %余为零时，即为最右边那一个
    end
    ynum = ceil(i/numnodx ); %ceil向上取整得到单元i是从下往上数第几个
    pi=[X(xnum),Y(ynum)];
end


function f = fun(x0,y0)
f = 0;
end

%% 单元荷载向量程序
function Fe = load_vector(i,numnodx, numnody, X, Y,NOD)
%% 单元各节点所对应的坐标
 Lnod = NOD(:,i);%节点的编号

p(1,:) = Node_coordinate(Lnod(1),numnodx,numnody,X,Y);%编号为i的节点的横纵坐标
p(2,:) = Node_coordinate(Lnod(2),numnodx,numnody,X,Y);
p(3,:) = Node_coordinate(Lnod(3),numnodx,numnody,X,Y);

%% 计算单元荷载向量所需要的a(1,2,3),b(1,2,3),c(1,2,3)
E = [p(1,1), p(1,2),1;p(2,1), p(2,2),1;p(3,1), p(3,2),1];
detE = det(E);
E0 = (detE)^(-1);
a1 = det([p(2,1),1;p(3,1),1]); b1 = (-1) * det([p(2,1),1;p(3,1),1]);
c1 = det([p(2,1),p(2,1);p(3,1),p(3,2)]);
a2 = det([p(3,1),1;p(1,1),1]); b2 = (-1) * det([p(3,1),1;p(1,1),1]); 
c2 = det([p(3,1),p(3,2);p(1,1),p(1,2)]);
a3 = det([p(1,2),1;p(2,1),1]); b3 = (-1) * det([p(1,2),1;p(2,2),1]); 
c3 = det([p(1,1),p(1,2);p(2,1),p(2,1)]);

%% 计算N（i,j,m）
N1=@(x,y)E0 * (a1 * x + b1 * y + c1) * fun(x,y);
N2=@(x,y)E0 * (a2 * x + b2 * y + c2) * fun(x,y);
N3=@(x,y)E0 * (a3 * x + b3 * y + c3) * fun(x,y);

%% 积分 
Fe = zeros(3,1);
t = @(x)2-x;
Fe(1) = integral2(N1,0,2,0,t); 
Fe(2) = integral2(N2,0,2,0,t); 
Fe(3) = integral2(N3,0,2,0,t); 
end

%% 单元刚度矩阵
function Ke = Stiffness_matrix(i,numnodx, numnody, X, Y,NOD)
%% 单位各节点所对应的坐标
Lnod = NOD(:,i); %节点的编号
p1 = Node_coordinate(Lnod(1),numnodx,numnody,X,Y);%编号为i的节点的横纵坐标
p2 = Node_coordinate(Lnod(2),numnodx,numnody,X,Y);
p3 = Node_coordinate(Lnod(3),numnodx,numnody,X,Y);
x1 = p1(1); y1 = p1(2); 
x2 = p2(1); y2 = p2(2);
x3 = p3(1); y3 = p3(2); 
%% 计算单元刚度矩阵所需要的a(i,j,m),b(i,j,m),c(i,j,m)
E = [x1, y1,1;x2, y2,1;x3, y3,1];
detE = det(E);
a1 = det([y2,1;y3,1]); a2 = det([y3,1;y1,1]); a3 = det([y1,1;y2,1]);
b1 = (-1) * det([x2,1;x3,1]); 
b2 = (-1) * det([x3,1;x1,1]); 
b3 = (-1) * det([x1,1;x2,1]);
% c1 = det([x2,y2;x3,y3]); c2 = det([x3,x3;x1,y1]); c3 = det([x1,y1;x2,y2]);
B = [a1, b1; a2, b2; a3,b3];
%% 计算准单元刚度矩阵
Ke = (2 * detE)^(-1) .* (B * B');
end