我试图将一种方法放在一起以计算不规则但凸面的凸面体的体积:
它使用三角剖分将多面体拆分为多个子四面体(简单)并独立计算体积,然后对所有子体积值求和。
但是,在测试中,我得到以下单元的奇怪结果-多维数据集。任何人都知道错误在哪里吗?
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
最佳答案
我建议对公式进行数值积分时使用高斯求积。这就是通常使用有限元方法完成的方式。您将从参数空间中的单位形状开始,然后将其转换为全局坐标。
您可能还考虑使用格林定理将体积积分转换为表面积分。这样可以使复杂形状的表面离散化。它特别适合带孔的复杂形状。
关于python - 不规则多面体的体积,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/19647345/