# Forming a solid from a set of intersecting planes

suppose I have equations of planes in 3D space in the form:

$A_ix + B_iy + C_iz = D_i$ for some values of i

how can use sage to find the solid enclosed by them? I need a plot of the surface of this solid. Actually I have many such sets of planes each forming a solid like this.(all concentric)

For example:

planes at x = -1, 1 in yz plane; planes at y = -1, 1 in xz plane and planes at z = -1, 1 in xy plane intersect and form a unit cube centred at the origin.

ty ndomes! But, its taking very long. I have like 12 or more equations (seen a maximum of 30 but can be more). I am now thinking of a different approach. We can easily make a function:

def inside(point , planes):
if the point is inside all the planes:
return 1
else:
return 0


so, now is there a good function in sagemath that can show all the region in space for which f(point) gives 1. (yea, its okay if the solids are filled as long as they're transparent enough. . Naively Taking lots of points and finding this solid(Like monte carlo) is very inefficient

Actually the program is taking very long time for these equations.

[-x - y - 1 == 0, -x - z - 1 == 0, -x + z - 1 == 0, -x + y - 1 == 0, -y- z - 1== 0, -y + z - 1 == 0, y - z - 1 == 0, y + z - 1 == 0, x - y - 1== 0, x - z - 1 == 0, x + z - 1 == 0, x + y - 1 == 0]


so, check it once. Thanks!

edit retag close merge delete

Sort by ยป oldest newest most voted
var("x y z")

P = vector((x,y,z))

eqns  = [ x==-1, x==1, y==-1, y==1, z==-1, z==1 ]
#eqns  = [ z==0, x+y+z ==1, -x+y+z == 1, -x-y+z==1, x-y+z==1]

indizes = range(len(eqns))

vertices = []
for C in combinations(indizes,3):
lse = [eqns[i] for i in C]
sol = solve(lse,x,y,z,solution_dict=True)
if sol:
vertices.append(P.substitute(sol[0]))

#for V in vertices:    print V

G = point3d(vertices,size=10,color='red')

" faces "
for eq in eqns:
poly = []
for V in vertices:
if(eq.substitute({x:V[0],y:V[1],z:V[2]})):
if not V in poly:
poly.append(V)

"""TODO    poly should be convex"""
shuffle(poly)
# print poly
try:
G += polygon(poly,alpha=0.5,color=(random(),random(),random()))
except:
pass

G.show(aspect_ratio=1)

more

I think we can find all vertices on each face and sort them according to polar angle about the centroid of the face. Then we can plot each face separately.

( 2012-06-13 00:44:36 +0200 )edit

I found another way to reorder the polygons. See http://www.hr.shuttle.de:9000/home/pub/104/

( 2012-06-14 10:22:04 +0200 )edit

If you write an equation as list [A,B,C,D] instead of Ax + By + C*z == D, you can easily use matrix operations to find the vertices. This seems to be much faster.

equations =[[1,0,0,1],[1,0,0,-1],[0,1,0,1],[0,1,0,-1],[0,0,1,1],[0,0,1,-1]]

def get_matrix(L,K):
M = matrix(QQ,3,3)
M.set_row(0,L[K[0]][:3])
M.set_row(1,L[K[1]][:3])
M.set_row(2,L[K[2]][:3])
return  M

def get_vector(L,K):
return vector(QQ,(L[K[0]][3],L[K[1]][3],L[K[2]][3]))

def find_vertices_2(eqns):
vertices = []
for C in combinations(range(len(eqns)),3):
M = get_matrix(eqns,C)
b = get_vector(eqns,C)
if M.is_invertible():
P = M.solve_right(b)
if not P in vertices:
vertices.append(P)
return vertices


I compared this function find_vertices_2 to a function find_vertices_1 based on solving systems of linear equations:

timeit('find_vertices_1(eqns)',number=20)


20 loops, best of 3: 1.3 s per loop

timeit('find_vertices_2(equations)',number=20)


20 loops, best of 3: 15.8 ms per loop

more

Does your code go into an infinite loop if the polygon is concave? I think my planes can give rise to concave polygons. (Actually I am trying to get the brillouin zones of a 3d lattice). See this for example: http://www.wolfram.com/products/publicon/samples/default/HTMLFiles/DefaultSample_90.gif

( 2012-06-16 08:13:51 +0200 )edit

I took for granted that all faces are convex. May be I misunderstood you. There is still your idea to order the vertices of each face by the respective rotation angle.

( 2012-06-16 09:10:10 +0200 )edit

ok I"m marking this as solved because I think I can complete the code for sorting the vertices and joining them in order

( 2012-06-18 07:16:12 +0200 )edit