二维有限元，线性插值

设置u=-（x*x+y*y），c=(x+y),可得f=6*（x+y），设置所有边界条件为dirichlet边界条件，其他条件应该也不复杂。boundaryedge矩阵是自己对着生成网格给出来的。感觉最难的地方就是在计算单元刚度矩阵的时候，因为使用了坐标变换，变成平面的标准三角形。（xi,yi),(xj,yj),(xm,ym)分别对应到（0，0），（1，0），（0，1）。

clc
clear all
n=10;
m=10;
%%计算一个二维网格，长度，宽度都为1，均分成10份，一共200个单元，121格节点，采用线性插值
%%P矩阵
for i=1:n+1
    for j=1:m+1
        p(1,(i-1)*(m+1)+j)=1/m*(i-1);
        p(2,(i-1)*(m+1)+j)=(j-1)*1/m;
    end
end
%%T矩阵
a=1;
b=2;
for i=1:n
    for j=1:m
        T(:,2*((i-1)*n+j)-1)=[a,a+11,a+12];
        T(:,2*((i-1)*n+j))=[b,b-1,b+11];
        a=a+1;
        b=b+1;
    end
    a=a+1;
    b=b+1;
end
%%边界边矩阵
    boundaryedge=[1,1,1,12;1,21,12,23;1,41,23,34;1,61,34,45;1,81,45,56;1,101,56,67;1,121,67,78;
        1,141,78,89;1,161,89,100;1,181,100,111;1,181,111,112;1,183,112,113;1,185,113,114;1,187,114,115;
        1,189,115,116;1,191,116,117;1,193,117,118;1,195,118,119;1,197,119,120;1,199,120,121;
        1,200,121,110;1,180,110,99;1,160,99,88;1,140,88,77;1,120,77,66;1,100,66,55;1,80,55,44;1,60,44,33;1,40,33,22;1,20,22,11;
        1,20,11,10;1,18,10,9;1,16,9,8;1,14,8,7;1,12,7,6;1,10,6,5;1,8,5,4;1,6,4,3;1,4,3,2;1,2,2,1];
%%构造基函数
K=zeros((n+1)*(m+1),(n+1)*(m+1));
F=zeros((n+1)*(m+1),1);
for  element=1:2*n*m
    i=T(1,element);
    j=T(2,element);
    m=T(3,element);  
    xi=p(1,i);
    yi=p(2,i);
    xj=p(1,j);
    yj=p(2,j);
    xm=p(1,m);
    ym=p(2,m);
    %%单元面积
    delta=[xi,yi,1;xj,yj,1;xm,ym,1];
    delta_s=det(delta);
    s=0.5; 
    DELTA_S=(xi-xm)*(yj-ym);
%%    
syms x y

fun_c =@(x,y) (xi-xm)*x+(xj-xm)*y+xm+(yi-ym)*x+(yj-ym)*y+ym;

fun_ni= @(x, y)  6*((xi-xm)*x+(xj-xm)*y+xm+(yi-ym)*x+(yj-ym)*y+ym).*(x)/(2*s);  

fun_nj = @(x, y) 6*((xi-xm)*x+(xj-xm)*y+xm+(yi-ym)*x+(yj-ym)*y+ym).*(y)/(2*s); 

fun_nm = @(x, y) 6*((xi-xm)*x+(xj-xm)*y+xm+(yi-ym)*x+(yj-ym)*y+ym).*(-x - y + 1)/(2*s);

% 定义形状函数 (线性形状函数)
N1 = 1 - x - y;
N2 = x;
N3 = y;
% 计算形状函数的梯度
B = [diff(N1, x), diff(N2, x), diff(N3, x); 
     diff(N1, y), diff(N2, y), diff(N3, y)];

% 定义雅可比矩阵及其行列式
J = [xj - xi, xm - xi; yj - yi, ym - yi];
%J = [xi - xm, xj - xm; yi - ym, yj - ym];
detJ = det(J);
invJ = inv(J);

ymax=@(x) 1-x;
% 计算Ke (使用局部坐标系)
Ke = detJ * integral2(fun_c, 0, 1, 0,ymax)* B'*invJ'*invJ*B;

% 计算Fe
f1 = detJ * integral2(fun_ni, 0, 1, 0, ymax);
f2 = detJ * integral2(fun_nj, 0, 1, 0, ymax);
f3 = detJ * integral2(fun_nm, 0, 1, 0,ymax);
    
    %%
    Fe = [f1, f2, f3];
    K(T(:,element),T(:,element))= K(T(:,element),T(:,element))+Ke;
    F(T(:,element)) = F(T(:,element))+ Fe';

end
%%
for i=1:40
    if boundaryedge(i,1)==1
        j=boundaryedge(i,3);
        k=boundaryedge(i,4);
        K(j,:)=0;
        K(k,:)=0;
        K(j,j)=1;
        K(k,k)=1;
        F(j)=-p(1,j)^2-p(2,j)^2;
        F(k)=-p(1,k)^2-p(2,k)^2;
    end
end
u=K\F;
%%
matrix = reshape(u, 11, 11);
clear x,y
x=0:0.1:1;
y=0:0.1:1;
subplot(1,2,1);
[x,y]=meshgrid(x,y);
mesh(x,y,matrix);
hold on 
yy=-x.^2-y.^2;
subplot(1,2,2);
mesh(x,y,yy);

不知道是什么原因，计算结果有些误差

下面是数据第一个是计算结果，第二张是精确解

欢迎讨论。