首页 > 解决方案 > 如何四面体每个面的体积偏移

问题描述

四面体平面的体积如何偏移?

计算三角锥每个面的平行面,计算每个面在四个点的交点。体积是根据交点计算的坐标计算的。

①我该怎么办?

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

标签: pythonsympy

解决方案


你会如何用微积分做这样的事情?在考虑 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 维空间。在您的情况下,您的四面体是被 2d 形状切割的 3d 形状。使用绝对值将使您的形状保持在 x 轴上方(因此为正)。


推荐阅读