firtrinfem.m

possion方程的fem求解

% 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