卷不規則多面體

Question

我試圖把一個方法一起以計算不規則但否則凸pyhedron的體積： 它使用三角測量來分裂多面體成多個子四面體（單工）和獨立地計算體積，然後總結所有子卷值。

但是，我得到了奇怪的結果爲單位 - 在我的測試下面的立方體。任何人都知道錯誤所在的位置？

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

Answer 1

我會推薦在數值積分公式上使用高斯積分。這是通常使用有限元方法完成的方式。您將從參數空間中的單位形狀開始，並將其轉換爲全局座標。

您也可以考慮使用格林定理將體積積分轉換爲曲面積分。分離形狀複雜的表面將會更容易。它特別適用於帶孔的複雜形狀。