%有限单元法计算的场值
X=[-100 -100 -100 -50 -50 -50 0 0 0 50 50 50 100 100 100];%节点X坐标
Y=[0 20 40 0 20 40 0 20 40 0 20 40 0 20 40];%节点Y坐标
SZ=[4 5 5 5 5 5 5 8 8 8 8 11 11 11 11 12;... % i
2 2 3 6 4 7 8 9 7 10 11 12 10 13 14 14;... % j
1 4 2 3 7 8 6 6 10 11 9 9 13 14 12 15;]; % m 单元顶点i j m对应节点号
ND=15;%节点总数
NE=16;%单元总数
KE=zeros(ND);
K=zeros(ND);
%单元矩阵KE相加合成总矩阵系数K
for I=1:NE
x=X(SZ(:,I));
y=Y(SZ(:,I));
a(1)=y(2)-y(3);
a(2)=y(3)-y(1);
a(3)=y(1)-y(2);
b(1)=x(3)-x(2);
b(2)=x(1)-x(3);
b(3)=x(2)-x(1);
s=2*(a(1)*b(2)-a(2)*b(1));
for i=1:3
for j=1:3
KE(i,j)=(a(i)*a(j)+b(i)*b(j))/s;
kni=SZ(i,I);
knj=SZ(j,I);
K(kni,knj)=K(kni,knj)+KE(i,j);
end
end
end
%--------------------------------------------------------------------
%代入边界条件，将线性方程组KU=0转换成KU=B，解方程得到区域内节点场值
NDBJ=[1 2 3 4 6 7 9 10 12 13 14 15];%边界节点号
m=1;
G=6.67e-11;
R=50;
v=2.5*pi*R.^2;
D=100;
X=[-100:50:100 ];%节点X坐标
Y=[0:20:40];%节点Y坐标
[X,Y]=meshgrid(X,Y);
dg=2*G*v*D./(X.^2+(D+Y).^2);
surf(X,Y,dg);
Ul=dg(:)
for z=1:ND
if z==NDBJ(m)
B(z)=dg(z)*10e10*K(z,z);
m=m+1;
K(z,z)=10e10*K(z,z);
else
B(z)=0;
end
B=B(:);
end
Uj=inv(K)*B;
%列出计算场值U与理论场值dg，并求出其误差
error=[Uj,Ul,Uj-Ul]
