有问题请联系QQ：291873404
function [u,err,zuiERR,Ujing]=possion(p,e,t,possionf,fjing)
%本程序为解决－（u''xx+u''yy）=f,

%输入p,e,t为几何图形分解产生而来的
%输入f为句柄函数
%g为边界条件
N=size(t' ,1);
M=size(p,2);
edge_Num=length(e);

%以下程序为找边界点


nodebian=1;
for h=1:edge_Num
    i=e(1,h);
    j=e(2,h);
    nodebian=[nodebian j]; 
end
nodebian=sort(nodebian);
nodebian(1)=[];


mm=length(nodebian);
K=zeros(M);
f=zeros(M,1);
%以下程序为计算矩阵刚度

for h=1:N
    i=t(1,h);
    j=t(2,h);
    k=t(3,h);
    sm=[1,p(1,i),p(2,i);1,p(1,j),p(2,j);1,p(1,k),p(2,k)];
    s(h)=det(sm)/2;
    
    Kii=[(p(2,j)-p(2,k))^2+(p(1,k)-p(1,j))^2]/(4*s(h));
    K(i,i)=K(i,i)+Kii;
    
      Kjj=[(p(2,k)-p(2,i))^2+(p(1,k)-p(1,i))^2]/(4*s(h));
    K(j,j)=K(j,j)+Kjj;
    
      Kkk=[(p(2,j)-p(2,i))^2+(p(1,i)-p(1,j))^2]/(4*s(h));                                                                                                                                   
    K(k,k)=K(k,k)+Kkk;
         
    Kij=[(p(2,j)-p(2,k))*(p(2,k)-p(2,i))+(p(1,k)-p(1,j))*(p(1,i)-p(1,k))]/(4*s(h));
    K(i,j)=K(i,j)+Kij;
     K(j,i)=K(j,i)+Kij;
    
     Kjk=[(p(2,k)-p(2,i))*(p(2,i)-p(2,j))+(p(1,i)-p(1,k))*(p(1,j)-p(1,i))]/(4*s(h));
    K(j,k)=K(j,k)+Kjk;
     K(k,j)=K(k,j)+Kjk;
     
      Kki=[(p(2,i)-p(2,j))*(p(2,j)-p(2,k))+(p(1,j)-p(1,i))*(p(1,k)-p(1,j))]/(4*s(h));
    K(k,i)=K(k,i)+Kki;
     K(i,k)=K(i,k)+Kki; 
     
     B=[p(1,i),p(2,i);p(1,j),p(2,j);p(1,k),p(2,k)];
x=sanjiaojifen('possionf',B);
  f(i)=f(i)+x(1);
  f(j)=f(j)+x(2);
  f(k)=f(k)+x(3);
end


neibian=0;
for i=1:M
    k=length(find(i==nodebian));
    if(k==0)
        neibian=[neibian i];
    end
end
neibian(1)=[];

 K(nodebian,:)=[];
  K(:,nodebian)=[];
f(nodebian)=[];
% ux=liezhuxiaoyuan(K,f);
v=sparse(K);
ux=v\f;
u(neibian)=ux;
u(nodebian)=0;
%v=K*u-f;


for i=1:M
    Ujing(i)=feval(fjing,p(1,i),p(2,i));
end
err=abs(u-Ujing);
zuiERR=max(abs(u-Ujing));