稳定渗流场的有限元

k=1; a=ones（83）; b=ones（83）; c=ones（83）; A=ones（81）; 
Ke=zeros（998）;Kq=zeros（98）;K3=zeros（998）;K1=zeros（998）;
G=zeros（398）;Ce=zeros（328）;
Q=1;%单位面积产量（井壁处流速）
Pe=10;%定压边界
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%单元集合阵
G（111）=1;G（251）=1;G（321）=1;
G（152）=1;G（232）=1;G（322）=1;
G（113）=1;G（243）=1;G（353）=1;
G（154）=1;G（264）=1;G（334）=1;
G（145）=1;G（275）=1;G（355）=1;
G（176）=1;G（286）=1;G（356）=1;
G（187）=1;G（297）=1;G（357）=1;
G（158）=1;G（298）=1;G（368）=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%节点坐标阵
C=[0 6;5 6;10 6;0 3;5 3;10 3;0 0;5 0;10 0];%所有点的坐标
for i=1:8
   Ce（::i）=G（::i）*C;  %单元的节点坐标阵
   a（1i）=Ce（21i）*Ce（32i）-Ce（22i）*Ce（31i）;b（1i）=Ce（22i）-Ce（32i）;c（1i）=Ce（31i）-Ce（21i）;
   a（2i）=Ce（31i）*Ce（12i）-Ce（32i）*Ce（11i）;b（2i）=Ce（32i）-Ce（12i）;c（2i）=Ce（11i）-Ce（31i）;
   a（3i）=Ce（11i）*Ce（22i）-Ce（12i）*Ce（21i）;b（3i）=Ce（12i）-Ce（22i）;c（3i）=Ce（21i）-Ce（11i）;
   A（i）=（a（1i）+a（2i）+a（3i））/2;
   N_x=[b（1i）b（2i）b（3i）]/2/A（i）;%partial N-x
   N_y=[c（1i）c（2i）c（3i）]/2/A（i）;%partial N-y
   Nxx=N_x‘*N_x;
   Nyy=N_y‘*N_y;
   %Ke（::i）=G（::i）‘*（Nxx+Nyy）*G（::i）*k;%k 是流度单元面积分；
   Ke（::i）=G（::i）‘*（Nxx+Nyy）*G（::i）*k*A（i）;
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%边界1上的积分
    n=0;x=zeros（21）;N_i