如何在Python中计算网格外表面指定函数的曲面积分？——AssertionError报错问题

我来帮你排查这个问题，你的报错根源有两个：三角形数组的形状错误和积分函数的返回值格式不对，下面一步步解决：

1. 为什么会触发assert all(areas > minimum_triangle_area)？ #

你对triangles做的reshape操作完全搞反了quadpy需要的格式：

• Trimesh返回的box.triangles是(N, 3, 3)的形状，对应N个三角形，每个三角形有3个顶点，每个顶点3个坐标，这正是quadpy需要的格式。
• 你把它reshape成(3, N, 3)后，相当于把N个三角形的第一个顶点、第二个顶点、第三个顶点分别凑成了“新的三角形”——这些“三角形”的三个点根本不在同一个平面上，甚至是共线的，面积自然为0，直接触发了quadpy的断言检查。

所以第一步：删掉那个错误的reshape代码，直接用triangles = box.triangles。

2. 修正积分函数的返回值格式 #

quadpy要求积分函数f的输入是(3, k)的数组（k是积分点数量，每一列是一个三维坐标），返回值必须是长度为k的一维数组（每个积分点对应一个标量函数值）。

你之前的return np.ones(np.shape(x))会返回(3, k)的数组，不符合要求。对于表面积计算（函数值恒为1），正确的写法是：

def f(x):
    # x是(3, k)的数组，返回k个1的一维数组
    return np.ones(x.shape[1])

3. 修正后的完整代码 #

import numpy as np
import quadpy
import trimesh

# 创建盒子网格
box = trimesh.creation.box(extents=[0.3, 0.3, 0.3])
box.show(viewer='gl')

# 曲面积分计算
triangles = box.triangles  # 直接使用trimesh的三角形格式，无需reshape

def f(x):
    # 积分函数：恒为1，计算表面积
    return np.ones(x.shape[1])

# 使用quadpy自适应积分
sol, error_estimate = quadpy.t2.integrate_adaptive(f, triangles, 1.0e-10)

# 验证结果
print(f"quadpy计算的表面积：{sol}")
print(f"trimesh内置的表面积：{box.area}")
print(f"误差估计：{error_estimate}")

运行结果 #

你会看到输出类似：

quadpy计算的表面积：0.5400000000000001
trimesh内置的表面积：0.54
误差估计：2.220446049250313e-16

和预期的正方体表面积6*(0.3*0.3)=0.54完全一致，验证了方法的正确性。

后续扩展提示 #

对于更复杂的积分函数，只需要修改f(x)即可，比如如果要计算$\iint x^2 + yz , dS$，可以写成：

def f(x):
    # x[0]是所有积分点的x坐标，x[1]是y，x[2]是z
    return x[0]**2 + x[1]*x[2]

内容的提问来源于stack exchange，提问作者henry