寻找用于在镶嵌域上进行数值积分的 Python 包

Posted 2023-03-12

【中文标题】寻找用于在镶嵌域上进行数值积分的 Python 包

【英文标题】：Looking for Python package for numerical integration over a tessellated domain

【发布时间】：2011-08-21 22:01:24

【问题描述】：

我想知道是否有人知道基于 numpy/scipy 的 python 包，可以在镶嵌域（在我的具体情况下，一个由 voronoi 单元限制的二维域）上对复杂的数值函数进行数值积分？过去，我使用了一些来自 matlab 文件交换的包，但如果可能的话，我希望留在我当前的 python 工作流程中。 matlab例程是

http://www.mathworks.com/matlabcentral/fileexchange/9435-n-dimensional-simplex-quadrature

对于正交和网格生成使用：

http://www.mathworks.com/matlabcentral/fileexchange/25555-mesh2d-automatic-mesh-generation

任何关于网格生成以及对该网格进行数值积分的建议都将不胜感激。

【问题讨论】：

您的输入是分散的数据 xi, yi, zi，还是（如您所写）一次只有一个 Voronoi 单元，大约有 6 个边？另外，有多少点或单元格——1000、1000000？

@Denis：它是 N 个点 xi,yi 的集合，标记了 N 个 voronoi 细胞的中心，其中 N 大约为 10^2 到 10^3。每个 voronoi 单元的边数不能保证是特定的数字。

乔希，标签网格 qhull (- python) 怎么样？

查看quadrature。您会在这里找到许多三角形和 tets 的方案。对于网格生成，有mshr、pygmsh、meshpy、frentos等。

【参考方案1】：

数值积分通常在三角形或矩形等简单元素上定义。也就是说，您可以将每个 polgon 分解为一系列三角形。运气好的话，您的多边形网格有一个三角形对偶，这将使集成更加容易。

quadpy（我的一个项目）在许多领域进行数值积分，例如三角形：

import numpy
import quadpy


sol, error_estimate = quadpy.t2.integrate_adaptive(
    lambda x: numpy.exp(x[0]),
    numpy.array(
        [
            [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]],
            [[1.0, 0.0], [1.0, 0.0], [1.0, 0.0], [1.0, 0.0], [1.0, 0.0]],
            [[0.0, 1.0], [0.0, 1.0], [0.0, 1.0], [0.0, 1.0], [0.0, 1.0]],
        ]
    ),
    1.0e-10,
)

print(sol)

3.5914091422921017

您还可以通过为三角形提供数百种方案之一来进行非自适应集成。

【讨论】：

quadpy 是一个很好的应用程序集成包，特别是有限元方法，谢谢！

【参考方案2】：

这直接在三角形上积分，而不是 Voronoi 区域，但应该是接近的。 （运行不同数量的点来查看？）它也适用于 2d、3d ...

#!/usr/bin/env python
from __future__ import division
import numpy as np

__date__ = "2011-06-15 jun denis"

#...............................................................................
def sumtriangles( xy, z, triangles ):
    """ integrate scattered data, given a triangulation
    zsum, areasum = sumtriangles( xy, z, triangles )
    In:
        xy: npt, dim data points in 2d, 3d ...
        z: npt data values at the points, scalars or vectors
        triangles: ntri, dim+1 indices of triangles or simplexes, as from
http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html
    Out:
        zsum: sum over all triangles of (area * z at midpoint).
            Thus z at a point where 5 triangles meet
            enters the sum 5 times, each weighted by that triangle's area / 3.
        areasum: the area or volume of the convex hull of the data points.
            For points over the unit square, zsum outside the hull is 0,
            so zsum / areasum would compensate for that.
            Or, make sure that the corners of the square or cube are in xy.
    """
        # z concave or convex => under or overestimates
    npt, dim = xy.shape
    ntri, dim1 = triangles.shape
    assert npt == len(z), "shape mismatch: xy %s z %s" % (xy.shape, z.shape)
    assert dim1 == dim+1, "triangles ? %s" % triangles.shape
    zsum = np.zeros( z[0].shape )
    areasum = 0
    dimfac = np.prod( np.arange( 1, dim+1 ))
    for tri in triangles:
        corners = xy[tri]
        t = corners[1:] - corners[0]
        if dim == 2:
            area = abs( t[0,0] * t[1,1] - t[0,1] * t[1,0] ) / 2
        else:
            area = abs( np.linalg.det( t )) / dimfac  # v slow
        zsum += area * z[tri].mean(axis=0)
        areasum += area
    return (zsum, areasum)

#...............................................................................
if __name__ == "__main__":
    import sys
    from time import time
    from scipy.spatial import Delaunay

    npt = 500
    dim = 2
    seed = 1

    exec( "\n".join( sys.argv[1:] ))  # run this.py npt= dim= ...
    np.set_printoptions( 2, threshold=100, edgeitems=5, suppress=True )
    np.random.seed(seed)

    points = np.random.uniform( size=(npt,dim) )
    z = points  # vec; zsum should be ~ constant
    # z = points[:,0]
    t0 = time()
    tessellation = Delaunay( points )
    t1 = time()
    triangles = tessellation.vertices  # ntri, dim+1
    zsum, areasum = sumtriangles( points, z, triangles )
    t2 = time()

    print "%s: %.0f msec Delaunay, %.0f msec sum %d triangles:  zsum %s  areasum %.3g" % (
        points.shape, (t1 - t0) * 1000, (t2 - t1) * 1000,
        len(triangles), zsum, areasum )
# mac ppc, numpy 1.5.1 15jun:
# (500, 2): 25 msec Delaunay, 279 msec sum 983 triangles:  zsum [ 0.48  0.48]  areasum 0.969
# (500, 3): 111 msec Delaunay, 3135 msec sum 3046 triangles:  zsum [ 0.45  0.45  0.44]  areasum 0.892

【讨论】：

感谢您的回复。我很快会对此进行测试，看看它是否适用于我的用例。

@Josh，用已知积分的测试函数来比较方法怎么样？

【参考方案3】：

scipy.integrate.dblquad 怎么样？它使用自适应正交规则，因此您放弃对集成网格的控制。不知道这对您的应用来说是好是坏。

【讨论】：

我不确定dblquad 能否轻松处理我需要的集成域。

听起来你有一个有点凌乱的离散网格，你正在整合它。有限元社区可能对您有所帮助，但这不是我的专业领域。