2

我正在尝试将一种方法放在一起来计算不规则但凸面体的体积：它使用三角剖分将多面体分成多个子四面体（单面体）并独立计算体积，然后对所有子体积值求和。

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

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

4

1 回答 1

1

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

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

于 2013-10-29T09:22:50.893 回答