fem二维（一）

 

（1） 

  （2）

方程如下

Ne=8;Np=9;
alpha_x=ones(Ne,1)*1;
alpha_y=ones(Ne,1)*1;
beta=zeros(Ne,1);
f=ones(Ne,1)*(-1);

nie=[4 2 1;4 5 2;5 3 2;5 6 3;7 5 4;7 8 5;8 6 5;8 9 6];
x1=2;x2=1;x3=0;
y1=0;y2=1;y3=2;
pos=[x1,y1;x2,y1;x3,y1; ...  %postion of points
    x1,y2;x2,y2;x3,y2; ...
    x1,y3;x2,y3;x3,y3];

b1=pos(nie(:,2),2)-pos(nie(:,3),2);
b2=pos(nie(:,3),2)-pos(nie(:,1),2);
b3=pos(nie(:,1),2)-pos(nie(:,2),2);
c1=pos(nie(:,3),1)-pos(nie(:,2),1);
c2=pos(nie(:,1),1)-pos(nie(:,3),1);
c3=pos(nie(:,2),1)-pos(nie(:,1),1);
s=1/2*(b1.*c2-b2.*c1);
%get K and b
K=zeros(Np,Np);
b=zeros(Np,1);
for ie=1:Ne
    Ke=1/4/s(ie)*(alpha_x(ie)*[b1(ie);b2(ie);b3(ie)]*[b1(ie),b2(ie),b3(ie)]+ ...
        alpha_y(ie)*[c1(ie);c2(ie);c3(ie)]*[c1(ie),c2(ie),c3(ie)])+...
        s(ie)/12*beta(ie)*(1+eye(3));
    K(nie(ie,:),nie(ie,:))=K(nie(ie,:),nie(ie,:))+Ke;    
    be=s(ie)*f(ie)/3;
    b(nie(ie,:))= b(nie(ie,:))+be;
end


 %apply boundary condition on AB and CD
 %phi(1:3)=1;phi(7:9)=0;
K(1:3,:)=0;K(1,1)=1;K(2,2)=1;K(3,3)=1;b(1:3)=1;
K(7:9,:)=0;K(7,7)=1;K(8,8)=1;K(9,9)=1;b(7:9)=0;

%resolve
phi=K^(-1)*b;