有限单元法的优点是对形状不规则的地质体和地形起伏模拟能力强，但它的计算程序编制相当复杂。要编制出一个效度好、适应能力强的高密度电阻率法正演有限元算法程序，需要做很多理论研究和实际工作才行。
2.1.2 变分问题的离散化
为了用有限单元法求解偏微分方程的边值问题，首先需要将其等价的变分问题离散化。也就是要将（2.21）式离散化，这包括求解区?的离散化（网格剖分）和泛函J(v)的离散化（导出线性方程）两个组成部分。 ① 网格剖分
由于受计算机的计算量、储存量及计算速度的限制，所有的数值计算方法都需要将无限大地中的电场分布限定在有限的求解区内进行讨论，并将连续的求解区域离散化，即网格剖分。
原则上讲，网格剖分可以根据所研究的地电条件灵活地将求解区限定为任何形状，并灵活采用任意适当的网格来剖分求解区。在二维情况下，最常用的办法是将求解区剖分为一系列互不重迭的三角单元，每个单元的顶点称为节点。剖分的规则是：1)整个剖分范围应该是越大越好，剖分范围越大，计算精度也就越高。但剖分得越大，计算量增大，计算速度降低。因此，在作网格剖分时，既要考虑精度，又要照顾速度；2)各单元的节点只能与相邻单元的节点相重，而不能成为相邻单元的边内点；3)要使每个单元内介质的电导率为常数，在不同电导率介质的内分界面上或整个求解区的外边界上，使三角形单元的一个边相互衔接，以这样的折线去逼近边界线；4)三角形单元最好是接近正三角形，不要使其中一个角很小或出现很大的钝角；5)在电场变化剧烈，电位参数的二阶或高阶导数大的地方（如：在异常体边沿、场源附近等），单元宜剖分的细一些，而在电场变化平缓的地方（如：远离异常体和场源），单元面积可取得大一些。当然，单元面积的变化应是渐变的；。
为使程序简化，L.R.Rijo采用了图2.3所示比较规则的矩形求解区和三角形剖分网格。试算结果表明，这种在矩形网格中布置交叉对称三角形剖分网格可以足够近似地模拟一般常见的不平地形和电性异常体，同时又能节省计算量，是经典的一种剖分方法。
在完成上述单元剖分的同时，于求解区内形成了若干个节点（包括内节点和边节点），待定的傅氏电位函数
v(?,x,z)在这些节点上的值，
称为节点函数值，自然它们是待定的未知值。这样，便用求解区有限个待定节点函数值来近似表征待定函数在连续空间中的分布。这些待定的节点函数值v1，v2，?，vM（M—总结点数），组成一组独立变量。所谓求待定函数的数值解，就是要确定这些节点的函数值。
②泛函的离散化 1）线性插值
为了计算二维变分问题（2.21）式中积分形式的泛函J(v)，需要知道整个求解区Ω内的v函数值。前面说到，求解区经过网格剖分后，可以用求解区内有限个待定节点函数值来近似表征待定函数在连续空间中的分布。这样，对于二维情况，就可以利用各节点的函数值在各单元内作线性插值来逐个单元计算v，也就是说，将v表示成该单元三个节点之函数值的函数，具体算法如下
设第e个单元的三个节点按逆时针方向的编号分别为i、j、m，其坐标记为
(xi,zi)、(xj,zj)、(xm,zm)，对应的节点函数值为vi、vj、vm（参见图2.2）。
图2.3 L.R.Rijo网格剖分示意图 假设各单元内，函数v(?,x,z)是线性变化的，即
v(x,z)?a?bx?cz （因为积分变量?是一个常数，不妨把 v(?,x,z)记为v(x,z)）
（2.24）
式中的的系数a、b、c可由单元上的三个节点的函数值和坐标算出。对于单元e，三个节点上都满足（2.24）式，故有
?vi?a?bxi?czi??vj?a?bxj?czj （2.25） ??vm?a?bxm?czm按克莱姆(Cramer)法则求解上述线性方程组，得
vixixjxmzizm1xizizj zm a?vjvmzj?1xj1xm?(aivi?ajvj?amvm)/2? （2.26a）
同样，有
b?(bivi?bjvj?bmvm)/2? （2.26b） c?(civi?cjvj?cmvm)/2? （2.26c）
式中
ai?xjzm?xmzjaj?xmzi?xizmam?xizj?xjzi单元面积
11xj21xm1xizibi?zj?zmbj?zm?zibm?zi?zjci?xm?xj??cj?xi?xm? （2.27） cm?xj?xi?? ??zj?zm1(bicj?bjci) （2.28） 2将（2.26a、2.26b、2.26c）式代入（2.24）式，便得到单元e内函数v线性插值的近似表示式
v(x,z)?1?ai?bix?ciz?vi??aj?bjx?cjz?vj??am?bmx?cmz?vm2???
（2.29） 或
v(x,z)?Ni(x,z)vi?Nj(x,z)vj?Nm(x,z)vm （2.30）
式中的Ni、Nj、Nm称为基函数，分别为
Ni(x,z)?(ai?bix?ciz)/2???Nj(x,z)?(aj?bjx?cjz)/2?? （2.31） Nm(x,z)?(am?bmx?cmz)/2???(x,z)?e
如果在每个单元上都作出了v(x,z)的这种近似，就能得到整个求解区内
v(x,z)的总体近似函数，这个函数在各个单元内是线性的，即它是一个分片为线
性的函数。对任意两个相邻单元来说，近似函数在公共边上的值，被两个端节点的函数值唯一决定，所以，总体近似函数在整个求解区自然是连续的。
2）单元分析
至此，二维变分问题（2.21）式中的泛函——整个求解区Ω上的积分J(v)，可以分解为各个单元e上的积分之和，即
J(v)??Je(v) （2.32）
e而Je(v)不难根据各单元e内函数v的线性插值近似表示式（2.29）或（2.30）算出。因此，首先应从分析一个三角单元e出发。对于当前问题的变分提法（2.21式），三个顶点分别为i，j，m的三角形单元e（参见图2.2）上的泛函
???v??v?Je??????()2?()2??2v2??2fv?ds????v2dl
?x?z??e???2?v???v=?e???()2?()2?ds??2?e??v2ds
?x?z?e?e-??v?Ik?(x?xk)?(z?zk)ds????v2dl (2.33)
ek?1?2N式中，?e表示单元e内的电导率值，由于它在每一个单元是常数，因此?e可以提到积分外。同样，积分变量?也可以提到积分外。
根据（2.24）式及（2.26b）、（2.26c），有
?v?b?(bivi?bjvj?bmvm)/2? （2.34a） ?x?v?c?(civi?cjvj?cmvm)/2? （2.34b） ?z可见，它们只与节点坐标及节点函数值有关，在一个单元内是常数，故可提到积分号外。于是，可将（2.33）式中第一项积分表示为
?v2???v2()?()?ds ????x?z?e??v???v=?()2?()2???dxdz
?z?e??x?1(bivi?bjvj?bmvm)2?(civi?cjvj?cmvm)2 (2.35) 4???同时，根据狄拉克?函数的性质，有
???(x?xA)?(z?zA)dxdz??e?1?0A?eA?e (2.36)
于是（2.33）式中积分的第三项为
IiviIjvjImvmv?Ik?(x?xk)?(z?zk)ds??? (2.37) ??nnnk?1ijmeN式中，Ii、Ij、Im分别表示在i、j、m点的供电电流强度。若场源（供电电极）不在某个节点上，则该点上的电流值为零；ni、nj、nm分别表示以节点i、j、（2.33）式的第二项积分 m为顶点的单元数。
2v??ds???Ni(x,y)vi?Nj(x,y)vj?Nm(x,y)vmdxdz ee??2 ?vi2222??N(x,z)dxdz?v(N(x,z))dxdz j??j??iee22?vm??(Nm(x,z))dxdz?2vivj??Ni(x,z)Nj(x,z)dxdzee?2vjvm??Nj(x,z).Nm(x,z)dxdze(2.
?2vmvi??Nm(x,z)Ni(x,z)dxdze38 )
为了计算上式中的积分，我们首先来看基函数Ni(x,z)、
Nj(x,z)、Nm(x,z)的几何意
义。对于三角单元e（i,j,m）中的任一点P(x,y)（见图2.4），由（2.27）式可算得
图2.4 基函数几何意义示图