使用 scipy 数值积分计算 3d 形状的体积

我已经编写了一个用于计算立方体和半空间相交体积的函数，现在我正在为其编写测试。

我尝试像这样以数字方式计算体积：

integral = scipy.integrate.tplquad(lambda z, y, x: int(Vector(x, y, z).dot(normal) < distance),
                                   -0.5, 0.5,
                                   lambda x: -0.5, lambda x: 0.5,
                                   lambda x, y: -0.5, lambda x, y: 0.5,
                                   epsabs=1e-5,
                                   epsrel=1e-5)

...基本上我对整个立方体进行积分，每个点根据它是否位于半空间内而获得值 1 或 0。 这变得非常慢（每次调用超过几秒）并且不断向我发出警告，例如

scipy.integrate.quadpack.IntegrationWarning: The integral is probably divergent, or slowly convergent

有没有更好的方法来计算这个体积？

特征函数的积分在数学上是正确的，但不实用。这是因为大多数集成方案旨在集成多项式在某种程度上确实如此，因此所有“相对平稳”的功能都相当好。然而，特征函数一点也不平滑。多项式式的积分不会有任何结果。

更适合的方法是首先构建域的离散版本，然后简单地总结小四面体的体积。

3D 离散化可以通过以下方式完成皮加尔网格 https://github.com/nschloe/pygalmesh（矿井接口项目CGAL https://www.cgal.org/）。下面的代码将截止立方体离散化为

您可以通过减少来提高精度max_cell_circumradius and/or max_edge_size_at_feature_edges，但是网格划分将需要更长的时间。此外，您可以指定“特征边”来准确解析相交边。即使像元尺寸最粗，这也会给您带来完全正确的结果。

import pygalmesh
import numpy


c = pygalmesh.Cuboid([0, 0, 0], [1, 1, 1])
h = pygalmesh.HalfSpace([1.0, 2.0, 3.0], 4.0, 10.0)
u = pygalmesh.Intersection([c, h])
mesh = pygalmesh.generate_mesh(
    u, max_cell_circumradius=3.0e-2, max_edge_size_at_feature_edges=1.0e-2
)


def compute_tet_volumes(vertices, tets):
    cell_coords = vertices[tets]
    a = cell_coords[:, 1, :] - cell_coords[:, 0, :]
    b = cell_coords[:, 2, :] - cell_coords[:, 0, :]
    c = cell_coords[:, 3, :] - cell_coords[:, 0, :]
    # omega = <a, b x c>
    omega = numpy.einsum("ij,ij->i", a, numpy.cross(b, c))
    # https://en.wikipedia.org/wiki/Tetrahedron#Volume
    return abs(omega) / 6.0

vol = numpy.sum(compute_tet_volumes(mesh.points, mesh.get_cells_type("tetra")))
print(f"{vol:.8e}")

8.04956436e-01