二维FEM建模及求解（入门级：数理方程到代码编写）

前言：好久没有更新了，最近一直在看关于脑电正问题建模的问题，重点是关于FEM头模型搭建的问题。网上好多资料要不就是纯理论讲解，要不就是直接上打断代码的。所以趁着现在的理解还是在由少到多的阶段，记录一下，面向群体是跟我一样基础比较薄弱，想从理论到代码编写捋一遍FEM方法的同志们。

先把要解决的问题摆出来：

1）相关数理方程化简步骤；

2）离散化，用什么原理？合理性？误差来源；

3）代码实现。

注：1）数理方程就是将物理现象搭建为数学模型，很有趣的一门课，感兴趣的可以了解；2）关于基函数的理解还不够到位，如果有人想了解基函数相关内容请移步。

1. 数学模型搭建和相关处理

数学模型为：

将公式弱化（试函数为v）：

 函数弱化思想见相应参考文章。

采用分部积分有：

注：公式弱化还能很好反映原微分公式的前提是积分的面积足够小，可以参考高数积分定理的证明思想，这也是需要进行离散化的原因。

2. 离散化（三角形剖分）

 需要求解的u用uh表示，j 表示第j个element，每个element包含3个顶点。为插值基函数，为3*1的关于坐标x、y的线性函数矩阵，为相应的权值（常数），为1*3的矩阵。

注：1）线性函数：是因为我们采用的是线性插值；

        2）3*1的函数矩阵是因为每个element由3个节点组成。第 j 个element可由如下公式表示：

 其中，N的推理和计算见下个文档。

对离散化的公式应用弱形式有：（取的转置为试函数）

 

3. 代码实现

代码实际上就是对和进行迭代积分。

具体见如下：（已经进行了详细注释）

function u = FEM_solver(nds,els,bcs)

nnd = size(nds,1);u = zeros(nnd,1);K =zeros (nnd,nnd);f = zeros(nnd,1);
for j = 1:size(els,1)
    %sti = stima(nds(els(j,:),:))
    K(els(j,:),els(j,:)) = K(els(j,:),els(j,:)) + stima(nds(els(j,:),:));    % 仅对角线上有值,stima为第j个K
    f(els(j,:)) = f(els(j,:)) + ones(3,1) * det([1,1,1;nds(els(j,:),:)'])/6;
end
freends = setdiff(1:nnd,bcs);     % 去除边界上的值，即边界条件u=0的描述。
u(freends) = K(freends,freends)\f(freends);

    function sti = stima(vertices)   
        % 计算第j个K
        ndim = size(vertices,2);
        J = [ones(1,ndim+1);vertices'];   % J可逆 ，构建目的：1）用于求偏导；2）计算element的面积
        B = J\[zeros(1,ndim);eye(ndim)];   % 方程组求解：求试函数的偏导
        sti = det(J).*B*B'/prod(1:ndim);   % 试函数的转置与试函数的乘积再乘以element的面积就是K
    end       
end

调用主程序为：

 其中网格划分采用的是distmesh工具包进行的。

相应链接为：DistMesh - A Simple Mesh Generator in MATLABhttp://persson.berkeley.edu/distmesh/

% 采用distmesh工具包进行网格生成
clear
figure
fd=@(p) sqrt(sum(p.^2,2))-1;
[p,t]=distmesh2d(fd,@huniform,0.2,[-1,-1;1,1],[]);
e=boundedges(p,t);
% 采用FEM_solver求解器求解
figure
nodes = p;
element = t;
boundaries = e;
x = FEM_solver(nodes,element,boundaries);
pp = {'FaceColor','Interp'};
trisurf(element,nodes(:,1),nodes(:,2),x,pp{:});
axis square;
colorbar

结果为：

1）网格划分结果：

2）FEM结果：

 

 最后，回答一下一开始提出的问题：

1）数学推理主要用到采用积分的形式对微分函数进行弱化，简化计算；然后就是采用分部积分将二阶导数转化为一阶。

2）离散化是采用基函数和相应系数的方式进行的，离散化是为了配合函数弱化进行的。误差来源是采用线性插值使得小范围内的u出现了平面。合理性在于当网格足够小时，误差将越来越小。

最后的最后，我也是刚入门，若是有疑问希望不吝赐教。

参考的文章：（这三篇文章给了我很大的启示，在此谢谢他们的分享）

你的第一个有限元求解器：仅十行 MATLAB - 知乎

什么是弱形式？ - 知乎

有限元分析简介及伽辽金法_feyily的博客-CSDN博客_伽辽金有限元方法