%计算每个三角单元
L = cell(600,1);
%b = zeros(600,1);
%c = zeros(600,1);
a = double(zeros(331,331));
for i = 1:1:600
m = B(i,1);
n = B(i,2);
l = B(i,3);
b(m) = (A(B(i,2),2)-A(B(i,3),2))*3/pi;
b(n) = (A(B(i,3),2)-A(B(i,1),2))*3/pi;
b(l) = (A(B(i,1),2)-A(B(i,2),2))*3/pi;
c(m) = (A(B(i,1),1)-A(B(i,3),1))*3/pi;
c(n) = (A(B(i,1),1)-A(B(i,3),1))*3/pi;
c(l) = (A(B(i,2),1)-A(B(i,1),1))*3/pi;
L(m,1) = {(A(B(i,2),1)*A(B(i,3),2)-A(B(i,3),1)*A(B(i,2),2)+(A(B(i,2),2)-A(B(i,3),2))*x+(A(B(i,1),1)-A(B(i,3),1))*y)*3/pi};
L(n,1) = {(A(B(i,3),1)*A(B(i,2),2)-A(B(i,1),1)*A(B(i,3),2)+(A(B(i,3),2))-A(B(i,1),2)*x+(A(B(i,1),1)-A(B(i,3),1))*y)*3/pi};
L(l,1) = {(A(B(i,1),1)*A(B(i,2),2)-A(B(i,2),1)*A(B(i,1),2)+(A(B(i,1),2))-A(B(i,2),2)*x+(A(B(i,2),1)-A(B(i,1),1))*y)*3/pi};
%三角形计算9次
for j = 1:1:3
M = B(i,j);
for k = 1:1:3
N = B(i,k);
%二重积分计算
syms x y;
I = int(L{M,1}.*L{N,1},y,0,1-x);
P = int(I,x,0,1);
P = eval(P);%类型转换
%计算刚度矩阵
a(M,N) = a(M,N)+((b(M).*b(N)+c(M).*c(N)));
end
end
end |