poisson_fem.m

poisson问题有限元求解通用程序(线性元). 压缩包内exp.m包含了一个实际问题,matlab下直接运行即可,如效率不高,可适当修改网格参数h.

%%%%%%%%%%%%%%%%%%%有限元方法(三角形元)求解二维Poisson混合边值问题%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-delta(c*u)+au     =f      in [x1,x2]*[y1,y2]
%u                  =f_D1   on T_D1 is (d1x1,d1y1) to (d1x2,d1y2)
%                   =f_D2   on T_D2 is (d2x1,d2y1) to (d2x2,d2y2)
%                   ...
%laplace(u)*n+bu    =f_N1   on T_N1 is (n1x1,n1y1) to (n1x2,n1y2)
%                   =f_N2   on T_N2 is (n2x1,n2y1) to (n2x2,n2y2)
%                   ...
%Eq=[a;b;c];                    ---constant
%F=[f]                          ---sym object
%Omega=[x1;x2;y1;y2];           ---constant
%Ed1=[  d1x1,d1y1,d1x2,d1y2;    ---constant
%       d2x1,d2y1;d2x2,d2y2;
%                   ...     ];                       
%F1=[f_D1;f_D2;...];            ---constant
%Ed2=[  n1x1,n1y1,n1x2,n1y2;
%       n2x1,n2y1,n2x2,n2y2;
%                   ...     ];
%F2=[f_D1;f_D2;...];            ---sym object
%%%%%%%%%%%%%%%%%%%%%Warn%%%%%%%%%%%%%%%%%%%%%%%%
%1  Ed must be constant,edge point,normal order.
%2  point in Ed1 or Ed2 must be in Ed.
%function [A,B,U,t]=Poisson_FEM(Eq,F,Omega,Ed1,F1,Ed2,F2,h);
%Function Introduction
%1  adapt to any 2D-area;
%2  adapt to D or N boundary condition
%To be continued to deal equation with function coefficients Eq. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序段开始%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A,B,U,t,pData]=Poisson_FEM(Eq,F,Omega,Ed,Ed1,F1,Ed2,F2,h);
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;                                      %全局变量，子函数调用需要
%变量定义%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gw=[1/3,1/3,1/3];
ap=[2/3,1/6,1/6;1/6,1/6,2/3;1/6,2/3,1/6];
a=Eq(1);b=Eq(2);c=Eq(3);                                                                            %方程参数---------------------------引用
f=F;

%以下过程需要提供边界点数据--坐标形式--边界正向
pData=[Ed,zeros(length(Ed),3)];
edData=[1:length(pData)]';

PreTreat(Ed1,Ed2,F1,F2);                                                                            %两类边界点预处理
%此处考虑N,D界点存在，故需特殊处理，即把界点定义为D类点
pData(ed2Data,3)        =2*ones(length(ed2Data),1);                                                 %N类点类型赋值
pData(ed1Data,3)        =ones(length(ed1Data),1);                                                   %D类点类型赋值  

MyDelaunay(h);                                                                                      %三角剖分
A               =zeros(length(pData));                                                              %拟刚度矩阵
B               =zeros(length(pData),1);                                                            %拟荷载矩阵
U               =ones(length(pData),1);                                                             %向量形式点数据(函数值)，未知点处为1，已知点处为函数值                                    
U(ed1Data)      =pData(ed1Data,4);                                                                  %边界点函数值赋值 
[unTag,temp]    =find(pData(:,3)~=1);                                                               %未知点确定

%%%%%%%%%%%%%%开始计算%%%%%%%%%%%%%%%%%%%%%%%%
tic;
y=MyStiff(0);
t=toc;
temp_A              =A;
for i=1:length(pData)
    temp_A(i,:)=temp_A(i,:).*U';
end
temp_B              =B-sum(temp_A(:,ed1Data)')';
temp_A              =A(unTag,unTag);
temp_B(ed1Data,:)   =[];
temp_U              =inv(temp_A)*temp_B;
U(unTag)            =temp_U;

%%%%%%%%%%%结果显示%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
trimesh(eData,pData(:,1),pData(:,2),U);
xlabel('X');ylabel('Y');zlabel('Z');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序段结束%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算拟刚度矩阵A和B%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%f      基函数索引值，可以是向量  
function y=MyStiff(x)
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
y=0;
for k=1:length(eData)
    
    temp_m              =eData(k,1:3);
    A(temp_m,temp_m)    =A(temp_m,temp_m)+MyEStiff(k);                                              %累计单元刚度   
    temp_eload=MyELoad(k);
    B(temp_m,1)         =B(temp_m,1)+temp_eload;                                                    %累计单元荷载
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算单元刚度矩阵
%e      三角形单元索引值     
function y=MyEStiff(e)
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
y=zeros(3,3);
temp_t      =eData(e,:);
temp_v      =EFunction(e,temp_t);                                                                   %单元上基函数向量u
temp_V      =temp_v*temp_v';                                                                        %单元上基函数乘积矩阵u*v
temp_dv     =[diff(temp_v,'x'),diff(temp_v,'y')];                                                   %单元上基函数梯度乘积向量du
temp_Dv     =temp_dv*temp_dv';                                                                      %单元上基函数梯度乘积矩阵du*dv

y           =y+TInt(c*temp_Dv+a*temp_V,e);                                                          %单元刚度贡献1,2

temp_side   =[temp_t',[temp_t(2:3)';temp_t(1)]];
temp_tag    =SideOnEdge(temp_side(:,1),temp_side(:,2),edData);
temp_s1t    =find(temp_tag>=1 & temp_tag<=2);
if length(temp_s1t)>0
    temp_s1     =temp_side(temp_s1t,:);
    temp_n      =Nor(pData(temp_s1(1),1:2),pData(temp_s1(2),1:2));                                  %边界T_D上单元法向
    y           =y+c*LInt(temp_dv*temp_n'*temp_v',pData(temp_s1(1),1:2),pData(temp_s1(2),1:2));     %边界T_D上荷载贡献
end
temp_s2t    =find(temp_tag>2);
if length(temp_s2t)>0
    temp_s2     =temp_side(temp_s2t,:);
    y           =y+b*c*LInt(temp_V,pData(temp_s2(1),1:2),pData(temp_s2(2),1:2));                    %边界T_N上荷载贡献
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算单元荷载向量%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%e      三角形单元索引值     
function y=MyELoad(e)
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
y           =zeros(3,1);
temp_t      =eData(e,1:3);
temp_v      =EFunction(e,temp_t);

y           =y+TInt(f*temp_v,e);                                                                    %三角形单元荷载贡献

temp_side   =[temp_t',[temp_t(2:3)';temp_t(1)]];
temp_tag    =SideOnEdge(temp_side(:,1),temp_side(:,2),edData);
temp_s2t    =find(temp_tag>2);
if length(temp_s2t)>0
    temp_s2     =temp_side(temp_s2t,:);
    temp_fN     =(pData(temp_s2(1),5)+pData(temp_s2(2),5))/2;                                           %边界T_N上f_N函数值
    y           =y+c*temp_fN*LInt(temp_v,pData(temp_s2(1),1:2),pData(temp_s2(2),1:2));                  %边界T_N上荷载贡献
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算三角形面积%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%e      三角形单元索引值
function y=TArea(e)
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
temp_m  =[1,1,1;pData(eData(e,:),1:2)']';
y       =abs(det(temp_m))/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算三角形中心%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%e      三角形单元索引值
function y=Center(e)
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
y       =sum(pData(eData(e,:),1:2))/3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算线性单元函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%y      单元函数形式，object对象
%e      三角形单元索引值
%p      点索引值，可以是向量
function y=EFunction(e,p)
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
y       =[];
temp_s  =TArea(e);
for i=1:length(p)
    temp_tag    =find(eData(e,:)==p(i));
    if length(temp_tag)==0
        temp_f      =0;
    else
        temp_m                  =sym([1,1,1;pData(eData(e,:),1:2)']');
        temp_m(temp_tag(1),:)   =sym('[1,x,y]');
        temp_f                  =det(temp_m)/(2*temp_s);
    end
    y           =[y;temp_f];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算三角形单元上的高斯积分%%%%%%%%%%%%%%%%%%%%%%
function y=TInt(F,e)
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
temp_t          =eData(e,:);
temp_s          =TArea(e);
temp_p          =ap*pData(temp_t,1:2);
[temp_r,temp_c] =size(F);
y               =zeros(temp_r,temp_c);
for i=1:temp_r
    for j=1:temp_c
        temp_F  =[];
        for k=1:3
            temp_F  =[temp_F;double(subs(F(i,j),{'x','y'},{temp_p(k,1),temp_p(k,2)}))];
        end
        y(i,j)  =temp_s*gw*temp_F;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%判断边ij是否在边界上%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=SideOnEdge(s,e,ed);
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
[temp_r,temp_c] =size(s);
y               =zeros(temp_r,1);
for i=1:length(s)
    temp_s      =find(ed==s(i));
    temp_e      =find(ed==e(i));
    if length(temp_s)*length(temp_e)>0 & (temp_e-temp_s==1 | temp_e-temp_s==1-length(ed))
        y(i)    =pData(s(i),3)*pData(e(i),3);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%三角剖分%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=MyDelaunay(h);
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
temp_M              =max(pData);
temp_m              =min(pData);
temp_x1             =temp_m(1);     temp_x2             =temp_M(1);
temp_y1             =temp_m(2);     temp_y2             =temp_M(2);
[temp_xM,temp_yM]   =meshgrid(temp_x1:h:temp_x2,temp_y1:h:temp_y2);
temp_select         =[];
for i=1:((temp_x2-temp_x1)/h+1)*((temp_y2-temp_y1)/h+1)
    temp_p      =[temp_xM(i),temp_yM(i)];
    if PtInCurve(temp_p,pData(:,1:2))==1
        temp_select     =[temp_select;i];
    end
end
pData               =[pData;[temp_xM(temp_select),temp_yM(temp_select),zeros(length(temp_select),3)]];
eData               =delaunay(pData(:,1),pData(:,2));
%正向三角形
for i=1:length(eData)
    tV=pData(eData(i,2:3),1:2)-pData(eData(i,1:2),1:2);
    if det(tV)<0
        eData(i,1:3)=eData(i,3:-1:1);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%两类边界点数据初始化(索引)%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=PreTreat(Ed1,Ed2,F1,F2);
global a b c A B U pData eData edData ed1Data ed2Data f gw ap;
%D边界点数据初始化(索引)开始
ed1Data         =[];
ed2Data         =[];
[temp_r,temp_c] =size(Ed1);
for i=1:temp_r
    temp_tag    =find(pData(:,1)==Ed1(i,1) & pData(:,2)==Ed1(i,2));
    temp_start  =find(edData(:,1)==temp_tag);               %start在边界中的序号
    temp_tag    =find(pData(:,1)==Ed1(i,3) & pData(:,2)==Ed1(i,4));
    temp_end    =find(edData(:,1)==temp_tag);               %end在边界中的序号
    temp_ed1    =[];
    if length(temp_start)+length(temp_end)==2
        %D类边界赋值(索引)
        if temp_end>1
            temp_ed1        =edData(temp_start:temp_end);
        else
            temp_ed1        =edData(temp_start:length(edData));
            temp_ed1        =[temp_ed1;edData(1)];
        end
        %D类点函数值赋值
        pData(temp_ed1,4)   =double(subs(F1(i,:),{'x','y'},{pData(temp_ed1,1),pData(temp_ed1,2)}))+zeros(length(temp_ed1),1);
        ed1Data             =[ed1Data;temp_ed1];
    else
        msgbox('Ed1 error!');
    end
end
%D边界点数据初始化(索引)结束
 %N边界点数据初始化(索引)开始
[temp_r,temp_c] =size(Ed2);
for i=1:temp_r
    temp_tag    =find(pData(:,1)==Ed2(i,1) & pData(:,2)==Ed2(i,2));
    temp_start  =find(edData(:,1)==temp_tag);
    temp_tag    =find(pData(:,1)==Ed2(i,3) & pData(:,2)==Ed2(i,4));
    temp_end    =find(edData(:,1)==temp_tag);
    temp_ed2    =[];
    if length(temp_start)+length(temp_end)==2
        %N类点边界赋值(索引)
        if temp_end>1
            ed2Data         =[ed2Data;edData(temp_start:temp_end)];
        else
            ed2Data         =[ed2Data;edData(temp_start:length(edData))];
            ed2Data         =[ed2Data;edData(1)];
        end
        %N类点函数值赋值
        pData(temp_ed2,5)   =double(subs(F2(i,:),{'x','y'},{pData(temp_ed2,1),pData(temp_ed2,2)}))+zeros(length(temp_ed2),1);
        ed2Data             =[ed2Data;temp_ed2];
    else
        msgbox('Ed2 error!');
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算两点距离%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=LLength(p1,p2)
y=sqrt((p2-p1)*(p2-p1)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算中点%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=Mid(p1,p2)
y=(p1+p2)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算两点法向%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=Nor(p1,p2)
temp_v=[p2-p1];
y=[p2(2)-p1(2),p2(1)-p1(1)]/LLength(p1,p2);
if det([temp_v;y])>0
    y=-y;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算两点积分%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%F---必须是sym object对象
function y=LInt(F,p1,p2)
temp_tag        =p1==p2;
temp_f          =0;
[temp_r,temp_c] =size(F);
y               =zeros(temp_r,temp_c);
for i=1:temp_r
    for j=1:temp_c
        switch sum(temp_tag')
            case 0
                temp_f          =subs(F(i,j),'y',(p2(2)-p1(2))/(p2(1)-p1(1))*(sym('x')-1)+p1(2));
            case 1
                if temp_tag(1)==0
                    temp_f      =subs(F(i,j),'y',p1(2));
                else
                    temp_f      =subs(F(i,j),'x',p1(1));
                end
            otherwise
                temp_f          =0;
        end
    end
    if temp_f==0
        temp_y  =0;
    else
        temp_y  =int(temp_f,p1(2),p2(2));
    end
    y(i,j)      =double(temp_y);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%判断点point是否在封闭曲线curve内部%%%%%%%%%%%%%%
function y=PtInCurve(point,curve);
temp_delta  =curve-ones(length(curve),1)*point;
temp_d      =temp_delta(:,1).^2+temp_delta(:,2).^2;
temp_m      =min(temp_d);
temp_tag    =find(temp_d==temp_m);
if temp_tag(1)==length(curve)
    temp_v1     =curve(1,1:2)-curve(temp_tag(1),1:2);
    temp_v2     =curve(temp_tag(1),1:2)-curve(temp_tag(1)-1,1:2);
elseif temp_tag(1)==1
    temp_v1     =curve(temp_tag(1)+1,1:2)-curve(temp_tag(1),1:2);
    temp_v2     =curve(temp_tag(1),1:2)-curve(length(curve),1:2);
else
    temp_v1     =curve(temp_tag(1)+1,1:2)-curve(temp_tag(1),1:2);
    temp_v2     =curve(temp_tag(1),1:2)-curve(temp_tag(1)-1,1:2);
end
    temp_v      =point-curve(temp_tag(1),1:2);
if det([temp_v1;temp_v])>0 & det([temp_v2;temp_v])>0
    y       =1;
else
    y       =0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%