% function [k,b]=Firtrinfem(u,v,a0,b0)
% a 是x軸上的邊界長度 ,b是y軸上的邊界長度
% m是x軸上的等分數，n是y軸上的等分數，Nnode為節點總數
% node全局編碼數組，為一維數組，存放節點全局編號
% Elematrix單元號與節點局部編碼和全局編碼的二維數組
clear
clc
tic,
l=6;
u=2^l;v=2^l;a0=1;b0=1;
hx=a0/u;hy=b0/v;
N_node=(u+1)*(v+1); N_elem=2*u*v;
delta_e=hx*hy*0.5;

%生成節點總體編碼與各單元的節點局部編碼的矩陣
%本程序中的節點編碼及其排序是根據《偏微分方程數值解法》（陸金甫，關治）P240中例子 

for e=1:N_elem
    alphax(e)=1;
    alphay(e)=1;
    beta(e)=0;
    t=2*u;
    if mod(e,t)==0
        t=e/t-1;
        x(e)=ceil(e/2)+t;
        M_ju_zong(e,1)=x(e)+1;
        M_ju_zong(e,2)=x(e)+u+2;
        M_ju_zong(e,3)=x(e);
    else
        x(e)=ceil(e/2)+floor(e/t);
        if mod(e,2)==1
            M_ju_zong(e,1)=x(e)+u+1;
            M_ju_zong(e,2)=x(e);
            M_ju_zong(e,3)=x(e)+u+2;
        else
            M_ju_zong(e,1)=x(e)+1;
            M_ju_zong(e,2)=x(e)+u+2;
            M_ju_zong(e,3)=x(e);
        end
    end
end

%求出每個結點的橫縱坐標

for i=1:N_node
    x(i)=mod(i-1,u+1)*hx;
    y(i)=floor((i-1)/(u+1))*hy;
end

%初始化總體剛度矩陣
k=zeros(N_node,N_node);
b=zeros(N_node,1);

%計算ai,aj,am,bi,bj,bm

for e=1:N_elem
    i=M_ju_zong(e,1);
    j=M_ju_zong(e,2);
    m=M_ju_zong(e,3);
    be(1)=y(j)-y(m);
    be(2)=y(m)-y(i);
    be(3)=y(i)-y(j);
    ce(1)=x(m)-x(j);
    ce(2)=x(i)-x(m);
    ce(3)=x(j)-x(i);
    fe(i)=-2*pi^2*sin(pi*x(i))*sin(pi*y(i));
    fe(j)=-2*pi^2*sin(pi*x(j))*sin(pi*y(j));
    fe(m)=-2*pi^2*sin(pi*x(m))*sin(pi*y(m));
     
    deltae=0.5*(be(1)*ce(2)-be(2)*ce(1));
%計算delta_ij
     for i=1:3
         for j=1:3
             if (i~=j)
                 del_ij=1.0;
             else
                 del_ij=0.0;
             end
         end
     end
%計算單元剛度矩陣
     for i=1:3
         b(1)=deltae*fe(i)/3;b(2)=deltae*fe(j)/3;b(3)=deltae*fe(m)/3;
         for j=1:3
             kk(i,j)=(alphax(e)*be(i)*be(j)+alphay(e)*ce(i)*ce(j))/(4*deltae);
         end
     end
%合成單元剛度矩陣到剛度矩陣
     for i=1:3
        for j=1:3
            k(M_ju_zong(e,i),M_ju_zong(e,j))=k(M_ju_zong(e,i),M_ju_zong(e,j))+kk(i,j);
            b(M_ju_zong(e,i))=b(M_ju_zong(e,i))+b(i);
        end
     end
end

%邊界條件的處理

 k(1:u+1,:)=[];%去掉矩陣前u+1行
 k(:,1:u+1)=[];%去掉矩陣前u+1列
 k=k(1:N_node-2*(u+1),:);%去掉矩陣后u+1行
 k=k(:,1:N_node-2*(u+1));%去掉矩陣后u+1列
 b(1:u+1)=[];%去掉右端向量前u+1行
 b(N_node-2*(u+1)+1:N_node-(u+1))=[];%去掉右端向量后u+1行

 %去掉左邊界結點所對應的行列
 
 g=N_node-2*(u+1);p=0;
 for i=1:(g-p)
     if mod(i+p,u+1)==1&i+p<=g
         k(i,:)=[];
         k(:,i)=[];
         b(i)=[];
         p=p+1;
     end
 end
%去掉右邊界結點所對應的行列

g1=g-v+1;p=0;
 for i=1:(g1-p)
     if mod(i+p,u)==0&i+p<=g1
         k(i,:)=[];
         k(:,i)=[];
         b(i)=[];
         p=p+1;
     end
 end
%{
 spy(k);
 disp(k)
 disp(b);
 k;
 size(k)
 size(b)
 %}


toc