寻找Python套件以在细分领域进行数值积分 [英] Looking for Python package for numerical integration over a tessellated domain

问题描述

我想知道是否有人知道基于numpy/scipy的python软件包，以数字方式将复杂的数值函数集成在棋盘形区域(在我的特定情况下是由voronoi单元界定的2D区域)上?过去，我使用了matlab文件交换中的几个软件包，但如果可能的话，希望保留在当前的python工作流程中. Matlab例程是

I was wondering if anyone knew of a numpy/scipy based python package to numerically integrate a complicated numerical function over a tessellated domain (in my specific case, a 2D domain bounded by a voronoi cell)? In the past I used a couple of packages off of the matlab file exchange, but would like to stay within my current python workflow if possible. The matlab routines were

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

用于使用以下方法生成正交和网格:

for the quadrature and mesh generation using:

http://www.mathworks.com/matlabcentral/fileexchange/25555-mesh2d -自动生成网格

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

Any suggestions on mesh generation and then numerical integration over that mesh would be appreciated.

推荐答案

这直接在三角形上积分，而不是Voronoi区域，但应靠近三角形. (用不同数量的点运行才能看到?)它也可以在2d，3d ...

This integrates over triangles directly, not the Voronoi regions, but should be close. (Run with different numbers of points to see ?) Also it works in 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

这篇关于寻找Python套件以在细分领域进行数值积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！