MATLAB有限元二维编程（三角单元）

首先，原理我是参考了这位博主，博主讲的很细致，但是有一些细节没有详细的提及，在此把自己在参照该博主的matlab代码后自己编程的过程纪录一下，如果有人看博主的文章就已经懂了，那么我这篇文章可能就对这些游客多余了。https://blog.csdn.net/lusongno1/article/details/81125167
博主是自己画的网格，COMSOL有导出网格的功能，所以我是借助了COMSOL里画了一个简单的模型然后导出网格再用MATLAB处理的，自我感觉网格的大小可以自定义，所以选择这个方法。
处理网格，因为COMSOL导出的网格有自己的格式，且其边界单元有一个缺点，不会区分给出的节点是那条边上的，本问题中由于边界都取０，所以没什么影响，但是如果狄式边界取值为非零值，则COMSOL不会告诉你哪条边为非零值。当然，把边界都设为０，作为有限元二维编程入门足够了。网格的数据格式如下
进入正题，求解的题目为：
此时我有必要说明一下，希望大家不会犯我一样的错误，以前我以为借助编程求解函数，那么很多步骤其实自己只要了解就好，不用从头到尾的把函数求解一遍，感觉那样跟自己手动求解没什么两样。这种想法有缺陷，编程实际上就是把你整个推导过程写成程序，如果你连每一步怎么来的，接下来该怎么走都不清楚，又何谈把它编成程序，所以先把方程整个推导一遍是基础。在我理解看来，编程也就是在算积分，单元总装，求解部分有用，其他都是这几步的铺垫。
推导过程：

%读入网格
clear all;
filename='mm1.txt';
%%节点坐标
[node_num]=textread(filename,'%n',1);%节点个数
n_coor=zeros(node_num,2);%节点坐标
m=1;
[n1,n2]=textread(filename,'%n%n','headerlines',1);
n_coor(:,1)=n1(1:node_num);
n_coor_x=n1(1:node_num);
n_coor_y=n2(1:node_num);
n_coor(:,2)=n2(1:node_num);
%%边界点
bou_num=n1(node_num+1);
boundary=zeros(bou_num,2);
a1=node_num+2;
a2=a1+bou_num-1;
boundary(:,1)=n1(a1:a2);
boundary(:,2)=n2(a1:a2);
%%单元索引
ele_num=n1(node_num+2+bou_num);%单元个数
ele_index=zeros(ele_num,3);%单元索引
a3=node_num+5+bou_num;
[ele_index(:,1),ele_index(:,2),ele_index(:,3)]=textread(filename,'%n%n%n','headerlines',a3);

%%计算单个矩阵
K=sparse(node_num,node_num);
F=sparse(node_num,1);
for i=1:ele_num
    ke=caculate1(i,ele_index,n_coor);
    f=caculate2(i,ele_index,n_coor);   
    q1=ele_index(i,1)+1;
    q2=ele_index(i,2)+1;
    q3=ele_index(i,3)+1;
    q=zeros(1,3);
    q(1)=q1;
    q(2)=q2;
    q(3)=q3;
    K(q,q)=K(q,q)+ke;
    F(q,1)=F(q,1)+f;
end
%施加边界条件
b=zeros(node_num,1);
u_b=0;
K(1,1)=1;
for i=1:bou_num
    w1=boundary(i,1)+1;
    w2=boundary(i,2)+1;
   if(b(w1)==0)
       F(w1,1)=u_b;
       K(w1,:)=0;
       K(:,w1)=0;
       K(w1,w1)=1.0;
       b(w1)=1;
   end
   if(b(w2)==0)
       F(w2,1)=u_b;      
       K(w2,:)=0;
       K(:,w2)=0;
       K(w2,w2)=1;
       b(w2)=1;
   end   
end
u=KF;
subplot(1,2,1);
plot3(n_coor_x,n_coor_y,u);
%scatter3(n_coor(:,1),n_coor(:,2),u);
L =1;
nsamp = 1001;
xsamp = linspace(0,L,nsamp);%100等分区间中间有100个数
[X,Y] = meshgrid(xsamp,xsamp);
uexact = sin(pi*X).*sin(pi*Y);
uexact_re = reshape(uexact,nsamp,nsamp);
subplot(1,2,2);
mmm=mesh(xsamp,xsamp,uexact_re);%%%%%


function [k] = caculate1(i,ele_index,n_coor)
%UNTITLED11 此处显示有关此函数的摘要
%   此处显示详细说明
a1=ele_index(i,1)+1;
a2=ele_index(i,2)+1;
a3=ele_index(i,3)+1;
x1=n_coor(a1,1);
y1=n_coor(a1,2);
x2=n_coor(a2,1);
y2=n_coor(a2,2);
x3=n_coor(a3,1);
y3=n_coor(a3,2);
A=0.5*((x2*y3-x3*y2)-(y3-y2)*x1+y1*(x3-x2));
A=abs(A);
A=1/(2*A);
J=(x3-x1)*(y2-y3)-(y3-y1)*(x2-x3);
J1=A*[(y2-y3) (x3-x2);(y3-y1) (x1-x3)];

a11=J.*([1 0]*J1*(J1'*[1;0])).*0.5;
a12=J.*([0 1]*J1*(J1'*[1;0])).*0.5;
a13=J.*([-1 -1]*J1*(J1'*[1;0])).*(0.5);
a22=J.*([0 1]*J1*(J1'*[0;1])).*0.5;
a23=J.*([-1 -1]*J1*(J1'*[0;1])).*(0.5);
a33=J.*(([-1 -1]*J1)*(J1'*[-1;-1]))*(0.5);
k=[a11 a12 a13; a12 a22 a23;a13 a23 a33];
end

function [f] = caculate2(i,ele_index,n_coor)
%UNTITLED12 此处显示有关此函数的摘要
%   此处显示详细说明
a1=ele_index(i,1)+1;
a2=ele_index(i,2)+1;
a3=ele_index(i,3)+1;
x1=n_coor(a1,1);
y1=n_coor(a1,2);
x2=n_coor(a2,1);
y2=n_coor(a2,2);
x3=n_coor(a3,1);
y3=n_coor(a3,2);
J=(x3-x1)*(y2-y3)-(y3-y1)*(x2-x3);
f1=@(lam1,lam2) fun(x1*lam1+x2*lam2+x3*(1-lam2-lam1),y1*lam1+y2*lam2+y3*(1-lam2-lam1)).*lam1.*J;
f2=@(lam1,lam2) fun(x1*lam1+x2*lam2+x3*(1-lam2-lam1),y1*lam1+y2*lam2+y3*(1-lam2-lam1)).*lam2.*J;
f3=@(lam1,lam2) fun(x1*lam1+x2*lam2+x3*(1-lam2-lam1),y1*lam1+y2*lam2+y3*(1-lam2-lam1)).*(1-lam1-lam2).*J;
lam=@(lam1) 1-lam1;
g1=integral2(f1,0,1,0,lam);
g2=integral2(f2,0,1,0,lam);
g3=integral2(f3,0,1,0,lam);
f=zeros(3,1);
f(1)=g1;
f(2)=g2;
f(3)=g3;
end

function bx = fun(x,y)
bx = (2*pi^2)*sin(pi*x).*sin(pi*y);
end

真解和所求解对比图
If you have any questions, feel free to contact me or write your problems in the comment region.

Higher you climb, more view you will see.