常用网格单元的重心

有限单元的重心可以用来构造更加稳定的有限元函数空间. 例如在间断有限元方法中, 单元上的 \(k\) 阶单项式基函数可以定义为 $$ \{ x^{i}y^j: \ \text{with }0\leq i+j\leq k \}. $$ 也可以定义为一种数值更加稳定的形式 $$ \{ (x-\bar{x})^{i}(y-\bar{y})^j: \ \text{with }0\leq i+j\leq k \}, $$ 这里 \((\bar{x},\bar{y})\) 为单元的重心. 我们知道, 如果质量是均匀分布的, 每一个过重心 (质心) 的直线都将几何体分为质量相同的两个部分. 这个结论对空间中若干个散点(质量分别集中在点上)系统同样成立. 散点系统的重心是各个点坐标的算术平均. 但是, 不是所有的几何体的重心也是各顶点的算术平均, 因为质量是均匀分布在几何体上而不是各个顶点上.

单纯形

三角形, 平行四边形, 四面体, 正六面体等几何体通常被称为单纯形, 其重心坐标即为各顶点坐标的算术平均:

• 二维空间. 设 \(A_i(x_i,y_i),i=0,1,2,…,n-1\) 为所有顶点, 则重心坐标为

$$ (\bar{x},\bar{y})=\left(\frac{x_0+x_1+\dots+x_{n-1}}{n},\frac{y_0+y_1+\dots+y_{n-1}}{n}\right). $$

• 三维空间. 设 \(A_i(x_i,y_i, z_i),i=0,1,2,…,n-1\) 为所有顶点, 则重心坐标为

$$ (\bar{x},\bar{y}, \bar{z}) = \left(\frac{x_0+x_1+\dots+x_{n-1}}{n},\frac{y_0+y_1+\dots+y_{n-1}}{n}, \frac{z_0+z_1+\dots+z_{n-1}}{n}\right). $$

多边形

不规则多边形的重心一般不是顶点坐标的算术平均, 遵循如下计算公式: 设 \(A_i(x_i,y_i),\ i=0,1,2,…,n-1\) 为按逆(顺)时针排列的所有顶点, 则重心坐标 \((\bar{x},\bar{x})\) 为 $$ \bar{x} = \frac{1}{6A} \sum_{i=0}^{n-1} (x_i+x_{i+1})(x_iy_{i+1}-x_{i+1}y_i), \tag{Eq. 1} $$ $$ \bar{y} = \frac{1}{6A}\sum_{i=0}^{n-1}(y_i+y_{i+1})(x_iy_{i+1}-x_{i+1}y_i), \tag{Eq. 2} $$ 这里 \(A\) 为多边形的面积, 定义为 $$ A = \frac{1}{2} \sum_{i=0}^{n-1} (x_iy_{i+1}-x_{i+1}y_i), $$ 且另 \(x_n=x_0, \ y_n=y_0\).

为什么不是算术平均?

下面我们从理论上看一下. 如果 \((\bar{x},\bar{y})\) 为多边形 \(D\) 的重心, 则其满足: $$ \iint_D (x-\bar{x})\text{d}x\text{d}y=0,\quad \iint_D (y-\bar{y})\text{d}x\text{d}y=0, $$ 那么重心坐标为 (\(|D|:=\iint_D \text{d}x\text{d}y\) 为 \(D\) 的面积) $$ \bar{x}=\frac{\iint_D x \text{d}x\text{d}y}{|D|}, \quad \bar{y}=\frac{\iint_D y \text{d}x\text{d}y}{|D|}. \tag{Eq. 3} $$

Step 1: 计算 \(\iint_D x \text{d}x\text{d}y\). 由 Green 公式可得

$$ \iint_D x \text{d}x\text{d}y = \iint_D \frac{\partial}{\partial x}\left(\frac{x^2}{2}\right) \text{d}x\text{d}y = \int_{\partial D}\frac{x^2}{2}n_x dS, \tag{Eq. 4} $$ 这里 \(n_x\) 为 \(D\) 的单位外法向量 \(n=(n_x,n_y)\) 的第一个分量.

Step 2: 计算线积分. 多边形的第 \(i\) 个边 \(e_i=A_iA_{i+1}, i=0,1,\ldots,n-1\) 长度为 $$ L_i = \sqrt{(x_{i+1}-x_i)^2+(y_{i+1}-y_i)^2}, $$ 且 \(e_i\) 上的单位外法向量为 $$(n_x,n_y)=(y_{i+1}-y_i, x_i-x_{i+1})/L_i .$$ 为计算线积分, 我们作如下变量变换

$$ x(t) = x_i+(x_{i+1}-x_i)t, \quad y(t) = y_i+(y_{i+1}-y_i)t, \quad 0\leq t\leq 1, $$ 因此, $$ \begin{split} \int_{e_i}\frac{x^2}{2}n_x dS &=\int_0^1 \frac{1}{2} (x_i+(x_{i+1}-x_i)t)^2 \frac{(y_{i+1}-y_i)}{L_i} L_i dt \\ &=\frac{1}{6} (x_i+x_{i+1})(x_iy_{i+1}-x_{i+1}y_i), \end{split} $$

Step 3: 将上式从 \(i=0,1,2,\ldots,n-1\) 连加, 可得

$$ \iint_D x \text{d}x\text{d}y = \frac{1}{6}\sum_{i=0}^{n-1} (x_i+x_{i+1})(x_iy_{i+1}-x_{i+1}y_i), $$ 上式除以面积 \(|D|\) (后面会证明 \(|D|=A\)) 就得到重心坐标 \(\bar{x}\), 这就证明了 (Eq. 1). 类似地可以证明 (Eq. 2).

Step 4: 同样地可以得到多边形的面积就是上述定义的 \(A\):

$$ \begin{split} |D|&=\iint_D \text{d}x\text{d}y = \iint_D \frac{\partial}{\partial x}(x/2)+ \frac{\partial}{\partial y}(y/2) \text{d}x\text{d}y \\ &= \int_{\partial D}\frac{x}{2}n_x+\frac{y}{2}n_y dS \\ &= \frac{1}{2}\sum_{i=0}^{n-1}(x_iy_{i+1}-x_{i+1}y_i)=A. \end{split} $$

附录 A: 多边形重心计算 python 程序