三维高阶单元的数值积分优化方法 | 有限元基础编程
本期分享的是有限元编程过程中，对于三维高阶单元的数值积分优化方法。
看上去是枯燥的理论分析，实则添加了数值实现和可视化部分，力争有趣活泼，吸引更多的小伙伴了解有限元与数学的梦幻联动。比如，下面这样的：
 |  |
|---|---|
 |  |
 |  |
 |  |
背景描述
在常规三维高阶等参单元计算单元刚度时，常常使用的是高斯积分方案即： 、 、 。实际编程时，只需要三重循环即可完成，但是随即而来的就是每个单元需要计算较多的积分点，在整个系统刚度方程求解时，就会显得计算量较大。
这时，本篇推文应运而生，翻阅上古大佬的有限元理论文献，IRONS 教授对于三维高阶单元的数值积分方案进行数学上面的改进，整篇文献就一页半！（大佬就是大佬）文献 title:QUADRATURE RULES FOR BRICK BASED FINITE ELEMENTS，缩略图如下：
积分策略
对于一个三重积分的表达式，可拆解为如下形式进行求值：
文献中分别规定了六种积分方法，通过更改积分点的位置和权重来实现不同方法的求积。
感兴趣的小伙伴可直接翻阅原文献对于每个规则的定义以及误差分析，我直接将该篇文献的数值思想用 Python 实现，方便移植到自己的有限元程序中。
数值实现
现假设被积函数 ，使用文献中的方法求解：
对于每条规则的定义，代码注释都标注的十分清晰！
# -*- coding:UTF-8 -*-
# 定义被积函数 f(x, y, z)
def f(x, y, z):
return x**2 + y**2 + z**2
def rule_6(f):
"""
规则6：6点规则
----------------------------------
规则6基于六个面中点的数值积分法。积分点为六个坐标轴上的正负单位点，
每个点的权重相等，总积分由六个点的函数值的平均值乘以一个系数得到。
公式：
∫∫∫ f(ξ,η,ζ) dξ dη dζ ≈ (8/6) * (f(-1,0,0) + f(1,0,0) + f(0,-1,0) +
f(0,1,0) + f(0,0,-1) + f(0,0,1))
"""
B = 1.3333333333333333 # 8/6
a = 1.0
points_a = [
(-a, 0, 0), (a, 0, 0),
(0, -a, 0), (0, a, 0),
(0, 0, -a), (0, 0, a)
]
result_a = sum(f(x, y, z) for x, y, z in points_a)
result = B * result_a
return result
def rule_14(f):
"""
规则14：14点规则
----------------------------------
规则14使用14个积分点：6个面中点和8个立方体对角线上的点。
该规则用于提高积分精度，权重为B6和C8，积分点的坐标由b和c表示。
公式：
∫∫∫ f(ξ,η,ζ) dξ dη dζ ≈ B6 * (f(-b,0,0) + f(b,0,0) + f(0,-b,0) +
f(0,b,0) + f(0,0,-b) + f(0,0,b))
+ C8 * (f(-c,-c,-c) + f(c,-c,-c) + f(-c,c,-c) +
f(c,c,-c) + f(-c,-c,c) + f(c,-c,c) +
f(-c,c,c) + f(c,c,c))
其中，B6 = 0.886426593，C8 = 0.335180055，b = 0.795822426，c = 0.758786911
"""
B6 = 0.886426593
C8 = 0.335180055
b = 0.795822426
c = 0.758786911
points_b = [
(-b, 0, 0), (b, 0, 0),
(0, -b, 0), (0, b, 0),
(0, 0, -b), (0, 0, b)
]
points_c = [
(-c, -c, -c), (c, -c, -c),
(-c, c, -c), (c, c, -c),
(-c, -c, c), (c, -c, c),
(-c, c, c), (c, c, c)
]
result_b = sum(f(x, y, z) for x, y, z in points_b)
result_c = sum(f(x, y, z) for x, y, z in points_c)
result = B6 * result_b + C8 * result_c
return result
def rule_15a(f):
"""
规则15a：15点规则
----------------------------------
规则15a采用15个积分点，包括中心点、六个面中点和八个立方体对角线上的点。
该规则的积分精度较高，权重为A1、B6和C8。
公式：
∫∫∫ f(ξ,η,ζ) dξ dη dζ ≈ A1 * f(0,0,0)
+ B6 * (f(-1,0,0) + f(1,0,0) + f(0,-1,0) +
f(0,1,0) + f(0,0,-1) + f(0,0,1))
+ C8 * (f(-c,-c,-c) + f(c,-c,-c) + f(-c,c,-c) +
f(c,c,-c) + f(-c,-c,c) + f(c,-c,c) +
f(-c,c,c) + f(c,c,c))
其中，A1 = 1.564444444，B6 = 0.355555556，C8 = 0.537777778，c = 0.674199862
"""
A1 = 1.564444444
B6 = 0.355555556
C8 = 0.537777778
a = 1.0
c = 0.674199862
point_a = [(0, 0, 0)]
points_b = [
(-a, 0, 0), (a, 0, 0),
(0, -a, 0), (0, a, 0),
(0, 0, -a), (0, 0, a)
]
points_c = [
(-c, -c, -c), (c, -c, -c),
(-c, c, -c), (c, c, -c),
(-c, -c, c), (c, -c, c),
(-c, c, c), (c, c, c)
]
result_a = sum(f(x, y, z) for x, y, z in point_a)
result_b = sum(f(x, y, z) for x, y, z in points_b)
result_c = sum(f(x, y, z) for x, y, z in points_c)
result = A1 * result_a + B6 * result_b + C8 * result_c
return result
def rule_15b(f):
"""
规则15b：15点规则（另一种形式）
----------------------------------
规则15b与规则15a类似，但积分点和权重有所不同。
该规则的积分点包括中心点、六个面中点和八个立方体对角线上的点，
但权重和坐标均与规则15a不同。
公式：
∫∫∫ f(ξ,η,ζ) dξ dη dζ ≈ A1 * f(0,0,0)
+ B6 * (f(-b,0,0) + f(b,0,0) + f(0,-b,0) +
f(0,b,0) + f(0,0,-b) + f(0,0,b))
+ C8 * (f(-c,-c,-c) + f(c,-c,-c) + f(-c,c,-c) +
f(c,c,-c) + f(-c,-c,c) + f(c,-c,c) +
f(-c,c,c) + f(c,c,c))
其中，A1 = 0.712137436，B6 = 0.686227234，C8 = 0.396312395，
b = 0.848418011，c = 0.727662441
"""
A1 = 0.712137436
B6 = 0.686227234
C8 = 0.396312395
b = 0.848418011
c = 0.727662441
point_a = [(0, 0, 0)]
points_b = [
(-b, 0, 0), (b, 0, 0),
(0, -b, 0), (0, b, 0),
(0, 0, -b), (0, 0, b)
]
points_c = [
(-c, -c, -c), (c, -c, -c),
(-c, c, -c), (c, c, -c),
(-c, -c, c), (c, -c, c),
(-c, c, c), (c, c, c)
]
result_a = sum(f(x, y, z) for x, y, z in point_a)
result_b = sum(f(x, y, z) for x, y, z in points_b)
result_c = sum(f(x, y, z) for x, y, z in points_c)
result = A1 * result_a + B6 * result_b + C8 * result_c
return result
def rule_19(f):
"""
规则19：19点规则
----------------------------------
规则19使用了19个积分点：一个中心点、六个面中点和十二个边中点。
该规则提高了积分的精度，并且适用于复杂的多项式函数。
公式：
∫∫∫ f(ξ,η,ζ) dξ dη dζ ≈ A1 * f(0,0,0)
+ B6 * (f(-b,0,0) + f(b,0,0) + f(0,-b,0) +
f(0,b,0) + f(0,0,-b) + f(0,0,b))
+ D12 * (f(-d,-d,0) + f(d,-d,0) + f(-d,d,0) +
f(d,d,0) + f(-d,0,-d) + f(d,0,-d) +
f(-d,0,d) + f(d,0,d) + f(0,-d,-d) +
f(0,d,-d) + f(0,-d,d) + f(0,d,d))
其中，A1 = 2.074074074，B6 = -0.24691358，D12 = 0.617283951，
b = 0.774596669，d = 0.774596669
"""
A1 = 2.074074074
B6 = -0.24691358
D12 = 0.617283951
b = 0.774596669
d = 0.774596669
point_a = [(0, 0, 0)]
points_b = [
(-b, 0, 0), (b, 0, 0),
(0, -b, 0), (0, b, 0),
(0, 0, -b), (0, 0, b)
]
points_d = [
(-d, -d, 0), (d, -d, 0),
(-d, d, 0), (d, d, 0),
(-d, 0, -d), (d, 0, -d),
(-d, 0, d), (d, 0, d),
(0, -d, -d), (0, d, -d),
(0, -d, d), (0, d, d)
]
result_a = sum(f(x, y, z) for x, y, z in point_a)
result_b = sum(f(x, y, z) for x, y, z in points_b)
result_d = sum(f(x, y, z) for x, y, z in points_d)
result = A1 * result_a + B6 * result_b + D12 * result_d
return result
def rule_27a(f):
"""
规则27a：27点规则
----------------------------------
规则27a使用了27个积分点：一个中心点、六个面中点、八个立方体对角线点、
和十二个边中点。这是一个非常精确的积分规则，适用于高精度要求的计算。
公式：
∫∫∫ f(ξ,η,ζ) dξ dη dζ ≈ A1 * f(0,0,0)
+ B6 * (f(-b,0,0) + f(b,0,0) + f(0,-b,0) +
f(0,b,0) + f(0,0,-b) + f(0,0,b))
+ C8 * (f(-c,-c,-c) + f(c,-c,-c) + f(-c,c,-c) +
f(c,c,-c) + f(-c,-c,c) + f(c,-c,c) +
f(-c,c,c) + f(c,c,c))
+ D12 * (f(-d,-d,0) + f(d,-d,0) + f(-d,d,0) +
f(d,d,0) + f(-d,0,-d) + f(d,0,-d) +
f(-d,0,d) + f(d,0,d) + f(0,-d,-d) +
f(0,d,-d) + f(0,-d,d) + f(0,d,d))
其中，A1 = 0.788073483，B6 = 0.499369002，C8 = 0.478508449，
D12 = 0.032303742，b = 0.848418011，c = 0.652816472，d = 1.106412899
"""
A1 = 0.788073483
B6 = 0.499369002
C8 = 0.478508449
D12 = 0.032303742
b = 0.848418011
c = 0.652816472
d = 1.106412899
point_a = [(0, 0, 0)]
points_b = [
(-b, 0, 0), (b, 0, 0),
(0, -b, 0), (0, b, 0),
(0, 0, -b), (0, 0, b)
]
points_c = [
(-c, -c, -c), (c, -c, -c),
(-c, c, -c), (c, c, -c),
(-c, -c, c), (c, -c, c),
(-c, c, c), (c, c, c)
]
points_d = [
(-d, -d, 0), (d, -d, 0),
(-d, d, 0), (d, d, 0),
(-d, 0, -d), (d, 0, -d),
(-d, 0, d), (d, 0, d),
(0, -d, -d), (0, d, -d),
(0, -d, d), (0, d, d)
]
result_a = sum(f(x, y, z) for x, y, z in point_a)
result_b = sum(f(x, y, z) for x, y, z in points_b)
result_c = sum(f(x, y, z) for x, y, z in points_c)
result_d = sum(f(x, y, z) for x, y, z in points_d)
result = A1 * result_a + B6 * result_b + C8 * result_c + D12 * result_d
return result
def rule_27G(f):
"""
规则27G：3x3x3积-高斯规则
----------------------------------
规则27G是标准的三重高斯积分，使用了3x3x3个积分点，
适用于较高精度的积分计算，积分点权重分别为5/9, 8/9, 5/9。
公式：
∫∫∫ f(ξ,η,ζ) dξ dη dζ ≈ ∑∑∑ f(x,y,z) * wx * wy * wz
其中，积分点为：±√(3/5) 和 0，权重为：5/9, 8/9, 5/9
"""
points = [-0.7745966692414834, 0.0, 0.7745966692414834] # sqrt(3/5)
weights = [0.5555555555555556, 0.8888888888888888,
0.5555555555555556] # 5/9, 8/9, 5/9
result = 0.0
for x in points:
for y in points:
for z in points:
result += f(x, y, z) * weights[points.index(x)] * \
weights[points.index(y)] * weights[points.index(z)]
return result
def rule_64G(f):
"""
规则64G：4x4x4积-高斯规则
----------------------------------
规则64G是一个更高精度的高斯积分规则，使用了4x4x4个积分点，
适用于需要更高精度的积分计算，积分点权重分别为特定的四个值。
公式：
∫∫∫ f(ξ,η,ζ) dξ dη dζ ≈ ∑∑∑ f(x,y,z) * wx * wy * wz
其中，积分点为：±0.861136311594053, ±0.339981043584856，权重为：0.347854845137454, 0.652145154862546
"""
points = [-0.861136311594053, -0.339981043584856,
0.339981043584856, 0.861136311594053]
weights = [0.347854845137454, 0.652145154862546,
0.652145154862546, 0.347854845137454]
result = 0.0
for x in points:
for y in points:
for z in points:
result += f(x, y, z) * weights[points.index(x)] * \
weights[points.index(y)] * weights[points.index(z)]
return result
# 计算并输出所有规则的积分结果
print("规则6积分结果:", rule_6(f))
print("规则14积分结果:", rule_14(f))
print("规则15a积分结果:", rule_15a(f))
print("规则15b积分结果:", rule_15b(f))
print("规则19积分结果:", rule_19(f))
print("规则27a积分结果:", rule_27a(f))
print("规则27G积分结果:", rule_27G(f))
print("规则64G积分结果:", rule_64G(f))
根据大学数学知识可知：，程序运行结果如下：
规则6积分结果: 8.0
规则14积分结果: 8.000000001701789
规则15a积分结果: 7.999999997028948
规则15b积分结果: 8.000000001027692
规则19积分结果: 8.000000001411939
规则27a积分结果: 7.999999981223534
规则27G积分结果: 8.000000000000004
规则64G积分结果: 8.000000000000014
绘图
可借助 Pyvista 绘制每个规则下的积分点位置：
| |
|---|---|
| |
| |
| |
绘制代码如下：
import pyvista as pv
import numpy as np
def plot_integration_points(points, title):
plotter = pv.Plotter()
# 创建一个立方体
cube = pv.Cube(center=(0, 0, 0), x_length=2, y_length=2, z_length=2)
plotter.add_mesh(cube, style='wireframe', color='blue')
# 绘制积分点
for point in points:
sphere = pv.Sphere(radius=0.05, center=point)
plotter.add_mesh(sphere, color='red')
# 绘制坐标平面
x_plane = pv.Plane(center=(0, 0, 0), direction=(1, 0, 0), i_size=2, j_size=2)
y_plane = pv.Plane(center=(0, 0, 0), direction=(0, 1, 0), i_size=2, j_size=2)
z_plane = pv.Plane(center=(0, 0, 0), direction=(0, 0, 1), i_size=2, j_size=2)
plotter.add_mesh(x_plane, color='red', opacity=0.3)
plotter.add_mesh(y_plane, color='green', opacity=0.3)
plotter.add_mesh(z_plane, color='blue', opacity=0.3)
# 设置标题
plotter.add_text(title, position='upper_edge', color='black', font_size=10)
# 显示图形
plotter.show()
# 定义积分点
# 规则6的积分点
rule_6_points = np.array([[-1, 0, 0], [1, 0, 0], [0, -1, 0], [0, 1, 0], [0, 0, -1], [0, 0, 1]])
plot_integration_points(rule_6_points, "Rule 6 Integration Points")
# 规则14的积分点
b = 0.795822426
c = 0.758786911
rule_14_points = np.array([[-b, 0, 0], [b, 0, 0], [0, -b, 0], [0, b, 0], [0, 0, -b], [0, 0, b],
[-c, -c, -c], [c, -c, -c], [-c, c, -c], [c, c, -c],
[-c, -c, c], [c, -c, c], [-c, c, c], [c, c, c]])
plot_integration_points(rule_14_points, "Rule 14 Integration Points")
# 规则15a的积分点
c_15a = 0.674199862
rule_15a_points = np.array([[-1, 0, 0], [1, 0, 0], [0, -1, 0], [0, 1, 0], [0, 0, -1], [0, 0, 1],
[-c_15a, -c_15a, -c_15a], [c_15a, -c_15a, -c_15a],
[-c_15a, c_15a, -c_15a], [c_15a, c_15a, -c_15a],
[-c_15a, -c_15a, c_15a], [c_15a, -c_15a, c_15a],
[-c_15a, c_15a, c_15a], [c_15a, c_15a, c_15a],
[0, 0, 0]])
plot_integration_points(rule_15a_points, "Rule 15a Integration Points")
# 规则15b的积分点
b_15b = 0.848418011
c_15b = 0.727662441
rule_15b_points = np.array([[-b_15b, 0, 0], [b_15b, 0, 0], [0, -b_15b, 0], [0, b_15b, 0],
[0, 0, -b_15b], [0, 0, b_15b],
[-c_15b, -c_15b, -c_15b], [c_15b, -c_15b, -c_15b],
[-c_15b, c_15b, -c_15b], [c_15b, c_15b, -c_15b],
[-c_15b, -c_15b, c_15b], [c_15b, -c_15b, c_15b],
[-c_15b, c_15b, c_15b], [c_15b, c_15b, c_15b],
[0, 0, 0]])
plot_integration_points(rule_15b_points, "Rule 15b Integration Points")
# 规则19的积分点
b_19 = 0.774596669
rule_19_points = np.array([[-b_19, 0, 0], [b_19, 0, 0], [0, -b_19, 0], [0, b_19, 0],
[0, 0, -b_19], [0, 0, b_19],
[-b_19, -b_19, 0], [b_19, -b_19, 0],
[-b_19, b_19, 0], [b_19, b_19, 0],
[-b_19, 0, -b_19], [b_19, 0, -b_19],
[-b_19, 0, b_19], [b_19, 0, b_19],
[0, -b_19, -b_19], [0, b_19, -b_19],
[0, -b_19, b_19], [0, b_19, b_19],
[0, 0, 0]])
plot_integration_points(rule_19_points, "Rule 19 Integration Points")
# 规则27a的积分点
b_27a = 0.848418011
c_27a = 0.652816472
d_27a = 1.106412899
rule_27a_points = np.array([[-b_27a, 0, 0], [b_27a, 0, 0], [0, -b_27a, 0], [0, b_27a, 0],
[0, 0, -b_27a], [0, 0, b_27a],
[-c_27a, -c_27a, -c_27a], [c_27a, -c_27a, -c_27a],
[-c_27a, c_27a, -c_27a], [c_27a, c_27a, -c_27a],
[-c_27a, -c_27a, c_27a], [c_27a, -c_27a, c_27a],
[-c_27a, c_27a, c_27a], [c_27a, c_27a, c_27a],
[-d_27a, -d_27a, 0], [d_27a, -d_27a, 0],
[-d_27a, d_27a, 0], [d_27a, d_27a, 0],
[-d_27a, 0, -d_27a], [d_27a, 0, -d_27a],
[-d_27a, 0, d_27a], [d_27a, 0, d_27a],
[0, -d_27a, -d_27a], [0, d_27a, -d_27a],
[0, -d_27a, d_27a], [0, d_27a, d_27a],
[0, 0, 0]])
plot_integration_points(rule_27a_points, "Rule 27a Integration Points")
# 规则27G的积分点
g_points = [-np.sqrt(3/5), 0, np.sqrt(3/5)]
rule_27G_points = np.array([[x, y, z] for x in g_points for y in g_points for z in g_points])
plot_integration_points(rule_27G_points, "Rule 27G Integration Points")
# 规则64G的积分点
g_points_64 = [-0.861136311594053, -0.339981043584856, 0.339981043584856, 0.861136311594053]
rule_64G_points = np.array([[x, y, z] for x in g_points_64 for y in g_points_64 for z in g_points_64])
plot_integration_points(rule_64G_points, "Rule 64G Integration Points")
最后，希望通过这个小案例，可以带给大家有限元方面的“小惊喜”，在关注电脑高性能计算的同时，不要忘记有限元理论&数学分析方面的精华。
嗨~ 我是木木
很高兴认识你
以文字为媒介，分享有限元数值思想的点点滴滴
“积跬步，至千里”
登录后免费查看全文
著作权归作者所有，欢迎分享，未经许可，不得转载
首次发布时间：2024-09-15
最近编辑：2月前
硕士
有限元爱好者
相关推荐
热门文章
