常用网格单元的重心

有限单元的重心可以用来构造更加稳定的有限元函数空间. 例如在间断有限元方法中, 单元上的 \(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), \quad\text{(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), \quad\text{(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|}.\quad \text{(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, \quad\text{(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, $$ 因此, $$ \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), $$

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\):

$$ |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. $$

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

def polygon(coords):
    """coords - (i, n by 2) numpy array, coordinates of vertices
    
    Example:
        area, centroid = polygon(np.array([[0.,0], [0,-1],[1,0],[0,1],[-1,0]]))
    """
    n = coords.shape[0]
    new_coords = np.vstack((coords, coords[0:1,:]))
    x = new_coords[:,0]
    y = new_coords[:,1]
    vol = .5 * ( np.sum(x[0:n] * y[1:]) - \
                np.sum(y[0:n] * x[1:]) )
    centroid=np.zeros(2, dtype=np.float64)
    centroid[0] = np.sum((x[0:n]+x[1:])*(x[0:n]*y[1:]-x[1:]*y[0:n]))/(6*vol)
    centroid[1] = np.sum((y[0:n]+y[1:])*(x[0:n]*y[1:]-x[1:]*y[0:n]))/(6*vol)

    return vol, centroid

You may also like...

发表评论

电子邮件地址不会被公开。