python之不规则多面体的体积

我正在尝试将一种方法放在一起来计算一个不规则但在其他方面是凸的多面体的体积: 它使用三角剖分将多面体分割成多个子四面体(单形)并独立计算体积，然后将所有子体积值相加。

但是，在我的测试中，下面的单元 - 立方体得到了奇怪的结果。有人知道错误所在吗？

class Simplex(object): 
    def __init__(self,coordinates):   
        if not len(coordinates) == 4: 
            raise RuntimeError('You must provide only 4 coordinates!') 
        self.coordinates = coordinates 
 
    def volume(self): 
        ''' 
        volume: Return volume of simplex. Formula from http://de.wikipedia.org/wiki/Tetraeder 
        ''' 
        import numpy 
 
        vA = numpy.array(self.coordinates[1]) - numpy.array(self.coordinates[0]) 
        vB = numpy.array(self.coordinates[2]) - numpy.array(self.coordinates[0]) 
        vC = numpy.array(self.coordinates[3]) - numpy.array(self.coordinates[0]) 
 
        return numpy.abs(numpy.dot(numpy.cross(vA,vB),vC)) / 6.0   
''' 
Old code that did not work 
class Polyeder(object): 
    def __init__(self,coordinates): 
        if len(coordinates) < 4: 
            raise RuntimeError('You must provide at least 4 coordinates!') 
        self.coordinates = coordinates 
 
    def volume(self): 
        pivotCoordinate = self.coordinates[0] 
        volumeSum = 0 
 
        for i in xrange(1,len(self.coordinates)-3): 
            newCoordinates = [pivotCoordinate] 
            for j in xrange(i,i+3): 
                newCoordinates.append(self.coordinates[j]) 
            simplex = Simplex(newCoordinates) 
            volumeSum += simplex.volume() 
 
        return volumeSum 
''' 
 
class Polyeder(object): 
 
def __init__(self,coordinates): 
    ''' 
    Constructor 
    ''' 
 
    if len(coordinates) < 4: 
        raise RuntimeError('You must provide at least 4 coordinates!') 
 
    self.coordinates = coordinates 
 
 
def volume(self): 
    from pyhull.delaunay import DelaunayTri 
 
    delaunay = DelaunayTri(self.coordinates,joggle=True) 
    volume = 0 
    for vertices in delaunay.vertices: 
 
        coords = [self.coordinates[i] for i in vertices] 
        simplex = Simplex(coords) 
        volume += simplex.volume() 
 
 
    return volume 
 
coords = [] 
 
coords.append([0,0,0]) 
coords.append([1,0,0]) 
coords.append([0,1,0]) 
coords.append([0,0,1]) 
 
s = Simplex(coords) 
print s.volume() 
 
coords.append([0,1,1]) 
coords.append([1,0,1]) 
coords.append([1,1,0]) 
coords.append([1,1,1]) 
 
p = Polyeder(coords) 
print p.volume()  

旧结果打印输出为:

0.166666666667 
0.666666666667 

四面体的值应该是 1/6(正确)，但单位立方体的值应该是 1

新的结果是: 0.166666666667 1.0

请您参考如下方法:

我建议在数值积分公式上使用高斯求积法。这就是通常使用有限元方法完成的方式。您将从参数空间中的单位形状开始并将其转换为全局坐标。

您还可以考虑使用格林定理将体积积分转换为表面积分。以这种方式离散化复杂形状的表面会更容易。它特别适合带孔的复杂形状。