如何计算四面体每个面的体积偏移量

Question

如何计算四面体平面的体积偏移量

计算三角棱锥体每个面的平行面，并计算每个面在四个点处的交点。根据交点计算的坐标计算体积

①我该怎么办

from sympy import *
def myIntersection_of_3_planes(a,b,c):
    ab = a.intersection(b)
    abc = ab[0].intersection(c)
    return abc[0]

def myUnitVector(C):
    D=Matrix([[0], [0], [0]])
    myL=sqrt(C[0]**2+C[1]**2+C[2]**2)
    D[0]=C[0]/myL
    D[1]=C[1]/myL
    D[2]=C[2]/myL
    return D

def my3Point_Offset(myT0,myT1,myT2,my03):
    mmy03=myCrossProduct(Matrix(myT1)-Matrix(myT0),Matrix(myT2)-Matrix(myT0))
    mmy03=myUnitVector(mmy03)
    mmy03[0]=float(mmy03[0])*my03
    mmy03[1]=float(mmy03[1])*my03
    mmy03[2]=float(mmy03[2])*my03
    myT0S=Matrix([[0], [0], [0]])
    myT1S=Matrix([[0], [0], [0]])
    myT2S=Matrix([[0], [0], [0]])
    myT0S[0]=myT0[0]-mmy03[0]
    myT0S[1]=myT0[1]-mmy03[1]
    myT0S[2]=myT0[2]-mmy03[2]
    myT1S[0]=myT1[0]-mmy03[0]
    myT1S[1]=myT1[1]-mmy03[1]
    myT1S[2]=myT1[2]-mmy03[2]
    myT2S[0]=myT2[0]-mmy03[0]
    myT2S[1]=myT2[1]-mmy03[1]
    myT2S[2]=myT2[2]-mmy03[2]
    return myT0S,myT1S,myT2S

def my4Point_Intersection(myT0,myT1,myT2,myT3,my01,my02,my03):
        aa=my3Point_Offset(myT0, myT1, myT2, my03[0])
        bb=my3Point_Offset(myT0, myT2, myT3, my01[0])
        cc=my3Point_Offset(myT0, myT3, myT1, my02[0])
        p11, p12, p13 = map(Point3D, [aa[0],aa[1],aa[2]])
        p21, p22, p23 = map(Point3D, [bb[0],bb[1],bb[2]])
        p31, p32, p33 = map(Point3D, [cc[0],cc[1],cc[2]])
        a = Plane(p11, p12, p13)
        b = Plane(p21, p22, p23)
        c = Plane(p31, p32, p33)
        ans=myIntersection_of_3_planes(a, b, c)
        return ans

def myCrossProduct(A, B):
    return expand(Matrix([
            [A[1,0]*B[2,0]-A[2,0]*B[1,0]],
            [A[2,0]*B[0,0]-A[0,0]*B[2,0]],
            [A[0,0]*B[1,0]-A[1,0]*B[0,0]]
                                ]))

def myIntersection_of_3_planes(a,b,c):
    ab = a.intersection(b)
    abc = ab[0].intersection(c)
    return abc[0]

def myVolumeTetrahedron(myTetrahedron):
    my_PQRS=Matrix(myTetrahedron)
    my_S   =Matrix([myTetrahedron[3],myTetrahedron[3],myTetrahedron[3]])
    return float(det(my_PQRS[:3, :3].T - my_S.T)/6)

myTetrahedron=[[0,0,0],[10,0,0],[0,10,0],[0,0,10]]
myOffset     =[[0]    ,[0]     ,[0]     ,[0]     ]
myT0=my4Point_Intersection(myTetrahedron[0],myTetrahedron[1],myTetrahedron[2],myTetrahedron[3], myOffset[1],myOffset[2],myOffset[3])
myT1=my4Point_Intersection(myTetrahedron[1],myTetrahedron[2],myTetrahedron[3],myTetrahedron[0], myOffset[2],myOffset[3],myOffset[0])
myT2=my4Point_Intersection(myTetrahedron[2],myTetrahedron[3],myTetrahedron[0],myTetrahedron[1], myOffset[3],myOffset[0],myOffset[1])
myT3=my4Point_Intersection(myTetrahedron[3],myTetrahedron[0],myTetrahedron[1],myTetrahedron[2], myOffset[0],myOffset[1],myOffset[2])
myTetrahedron=[myT0,myT1,myT2,myT3]
print("#",'{:.3f}'.format(myVolumeTetrahedron(myTetrahedron)))
#
myTetrahedron=[[0,0,0],[10,0,0],[0,10,0],[0,0,10]]
myOffset     =[[1]    ,[1]     ,[1]    ,[1]      ]
myT0=my4Point_Intersection(myTetrahedron[0],myTetrahedron[1],myTetrahedron[2],myTetrahedron[3], myOffset[1],myOffset[2],myOffset[3])
myT1=my4Point_Intersection(myTetrahedron[1],myTetrahedron[2],myTetrahedron[3],myTetrahedron[0], myOffset[2],myOffset[3],myOffset[0])
myT2=my4Point_Intersection(myTetrahedron[2],myTetrahedron[3],myTetrahedron[0],myTetrahedron[1], myOffset[3],myOffset[0],myOffset[1])
myT3=my4Point_Intersection(myTetrahedron[3],myTetrahedron[0],myTetrahedron[1],myTetrahedron[2], myOffset[0],myOffset[1],myOffset[2])
myTetrahedron=[myT0,myT1,myT2,myT3]
print("#",'{:.3f}'.format(myVolumeTetrahedron(myTetrahedron)))
# (10*10)/2*10/3=166.667
# -166.667
# -119.877

②我在找一个有官方公式的网页

吼叫

我试着用numpy

from sympy import *
import numpy as np
def myIntersection_of_3_planes(a,b,c):
    ab = a.intersection(b)
    abc = ab[0].intersection(c)
    return abc[0]

def myUnitVector(C):
    return np.array(C)/sqrt(C[0]**2+C[1]**2+C[2]**2)

def my3Point_Offset(myT0,myT1,myT2,my03):
    mmy03=my03*np.array(myUnitVector(myCrossProduct(Matrix(myT1)-Matrix(myT0),Matrix(myT2)-Matrix(myT0))))
    return Matrix((np.array(myT0) - np.array(mmy03))[0]), \
           Matrix((np.array(myT1) - np.array(mmy03))[0]), \
           Matrix((np.array(myT2) - np.array(mmy03))[0])

def my4Point_Intersection(myT0,myT1,myT2,myT3,my01,my02,my03):
        aa=my3Point_Offset(myT0, myT1, myT2, my03[0])
        bb=my3Point_Offset(myT0, myT2, myT3, my01[0])
        cc=my3Point_Offset(myT0, myT3, myT1, my02[0])
        p11, p12, p13 = map(Point3D, [aa[0],aa[1],aa[2]])
        p21, p22, p23 = map(Point3D, [bb[0],bb[1],bb[2]])
        p31, p32, p33 = map(Point3D, [cc[0],cc[1],cc[2]])
        a = Plane(p11, p12, p13)
        b = Plane(p21, p22, p23)
        c = Plane(p31, p32, p33)
        ans=myIntersection_of_3_planes(a, b, c)
        return ans

def myCrossProduct(A, B):
    return expand(Matrix([
            [A[1,0]*B[2,0]-A[2,0]*B[1,0]],
            [A[2,0]*B[0,0]-A[0,0]*B[2,0]],
            [A[0,0]*B[1,0]-A[1,0]*B[0,0]]
                                ]))

def myIntersection_of_3_planes(a,b,c):
    ab = a.intersection(b)
    abc = ab[0].intersection(c)
    return abc[0]

def myVolumeTetrahedron(myTetrahedron):
    my_PQRS=Matrix(myTetrahedron)
    my_S   =Matrix([myTetrahedron[3],myTetrahedron[3],myTetrahedron[3]])
    return float(det(my_PQRS[:3, :3].T - my_S.T)/6)

myTetrahedron=[[0,0,0],[10,0,0],[0,10,0],[0,0,10]]
myOffset     =[[0]    ,[0]     ,[0]     ,[0]     ]
myT0=my4Point_Intersection(myTetrahedron[0],myTetrahedron[1],myTetrahedron[2],myTetrahedron[3], myOffset[1],myOffset[2],myOffset[3])
myT1=my4Point_Intersection(myTetrahedron[1],myTetrahedron[2],myTetrahedron[3],myTetrahedron[0], myOffset[2],myOffset[3],myOffset[0])
myT2=my4Point_Intersection(myTetrahedron[2],myTetrahedron[3],myTetrahedron[0],myTetrahedron[1], myOffset[3],myOffset[0],myOffset[1])
myT3=my4Point_Intersection(myTetrahedron[3],myTetrahedron[0],myTetrahedron[1],myTetrahedron[2], myOffset[0],myOffset[1],myOffset[2])
myTetrahedron=[myT0,myT1,myT2,myT3]
print("#",'{:.3f}'.format(myVolumeTetrahedron(myTetrahedron)))
#
myTetrahedron=[[0,0,0],[10,0,0],[0,10,0],[0,0,10]]
myOffset     =[[1]    ,[1]     ,[1]     ,[1]     ]
myT0=my4Point_Intersection(myTetrahedron[0],myTetrahedron[1],myTetrahedron[2],myTetrahedron[3], myOffset[1],myOffset[2],myOffset[3])
myT1=my4Point_Intersection(myTetrahedron[1],myTetrahedron[2],myTetrahedron[3],myTetrahedron[0], myOffset[2],myOffset[3],myOffset[0])
myT2=my4Point_Intersection(myTetrahedron[2],myTetrahedron[3],myTetrahedron[0],myTetrahedron[1], myOffset[3],myOffset[0],myOffset[1])
myT3=my4Point_Intersection(myTetrahedron[3],myTetrahedron[0],myTetrahedron[1],myTetrahedron[2], myOffset[0],myOffset[1],myOffset[2])
myTetrahedron=[myT0,myT1,myT2,myT3]
print("#",'{:.3f}'.format(myVolumeTetrahedron(myTetrahedron)))
# (10*10)/2*10/3=166.667
# -166.667
# -165.517

Answer 1

你会如何用微积分做这样的事情？在考虑3维空间之前，考虑2D：抛物线下的线切割区域。

y1=(x-1)^2
y2=x

你可以很容易地想象出面积是直线积分和抛物线积分之间的差

integral(y1) = 1/3(x-1)^3
integral(y2) = 1/2(x^2) ...

integral(y1)|(bounds) – integral(y2)|(bounds) = 1/3(x-1)^3 |(bounds) – 1/2(x^2) |(bounds)

但是，我们会注意到，这个示例包括x轴下方的负区域。这是可以用绝对值治愈的

同样的原理也适用于n维空间。在您的情况下，您的四面体是一个三维形状，由二维形状切片。使用绝对值将使形状保持在x轴上方（因此为正值）