三角形单元上的二元函数的四点高斯积分公式

\(\int_{▲ijk} f(x,y) \,dxdy = S_▲(w_1f(x_1,y_1)+w_2f(x_2,y_2)+w_3f(x_3,y_3)+w_4f(x_4,y_4))\)
式子中\(f(x,y)\)是定义在\(▲ijk\)上的二元函数，\((x_1,y_1),(x_2,y_2),(x_3,y_3),(x_4,y_4)\)是\(▲ijk\)上的4个高斯积分点，\(S_▲\)是三角形\(▲ijk\)的面积。
\((x_1,y_1),(x_2,y_2),(x_3,y_3),(x_4,y_4)\)的计算公式为
\(x_1 = 0.6x_i+0.2x_j+0.2x_k\)
\(y_1 = 0.6y_i+0.2y_j+0.2y_k\)
\(x_2 = 0.2x_i+0.6x_j+0.2x_k\)
\(y_2 = 0.2y_i+0.6y_j+0.2y_k\)
\(x_3 = 0.2x_i+0.2x_j+0.6x_k\)
\(y_3 = 0.2y_i+0.2y_j+0.6y_k\)
\(x_4 = (1/3)x_i+(1/3)x_j+(1/3)x_k\)
\(y_4 =(1/3)y_i+(1/3)y_j+(1/3)y_k\)
\(w_1,w_2,w_3=25/48,w_4=-9/16\)
python代码如下

def two_dem_gauss_int(vertices,f):
    # vertices是2*3的矩阵 里面存储（x_i,y_i），（x_j,y_j）（x_k,y_k）
    x_i = vertices[0][0]
    x_j = vertices[0][1]
    x_k = vertices[0][2]
    y_i = vertices[1][0]
    y_j = vertices[1][1]
    y_k = vertices[1][2]
    x_1 = x_i*0.6+x_j*0.2+x_k*0.2
    y_1 = y_i*0.6+y_j*0.2+y_k*0.2
    x_2 = x_i*0.2+x_j*0.6+x_k*0.2
    y_2 = y_i*0.2+y_j*0.6+y_k*0.2
    x_3 = x_i*0.2+x_j*0.2+x_k*0.6
    y_3 = y_i*0.2+y_j*0.2+y_k*0.6
    x_4 = x_i*(1/3)+x_j*(1/3)+x_k*(1/3)
    y_4 =y_i*(1/3)+y_j*(1/3)+y_k*(1/3)
    tri_S = math.fabs(1/2*np.linalg.det(np.array([[1,1,1],[x_i,x_j,x_k],[y_i,y_j,y_k]])))
    gauss_int = tri_S*((25/48)*(f(x_1,y_1))+(25/48)*(f(x_2,y_2))+(25/48)*(f(x_3,y_3))+(-9/16)*(f(x_4,y_4)))
    return gauss_int

下面测试2个算例

算例一

以\((0,0),(2,0),(0,2)\)为顶点的三角形上定义的二元函数\(f(x,y)=2-x-y\)的积分
此积分的精确值为
\(\int_{▲(0,0),(2,0),(0,2)} 2-x-y \,dxdy=4/3\)
用如下的代码调用函数

import numpy as np
def f(x,y):
    return 2-x-y
print(two_dem_gauss_int(np.array([[0,2,0],[0,0,2]]),f))

程序输出为
1.3333333333333337

算例二
以\((1,1),(6,1),(5,3)\)为顶点的三角形上定义的二元函数\(f(x,y)=7\)的积分
此积分的精确值为
\(\int_{▲(1,1),(6,1),(5,3)} 7 \,dxdy=35\)
用如下的代码调用函数

import numpy as np
def f(x,y):
    return 7
print(two_dem_gauss_int(np.array([[1,6,5],[1,1,3]]),f))

程序输出为
34.99999999999999