Ask Your Question
0

Forming a solid from a set of intersecting planes

asked 2012-06-12 11:26:42 +0200

ebs gravatar image

updated 2012-06-16 01:10:48 +0200

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 flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted
1

answered 2012-06-12 14:54:00 +0200

ndomes gravatar image

updated 2012-06-12 16:21:05 +0200

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  

"draw the solid  - start with vertices"        
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)
edit flag offensive delete link more

Comments

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.

ebs gravatar imageebs ( 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/

ndomes gravatar imagendomes ( 2012-06-14 10:22:04 +0200 )edit
1

answered 2012-06-16 06:51:27 +0200

ndomes gravatar image

updated 2012-06-16 09:11:58 +0200

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

edit flag offensive delete link more

Comments

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

ebs gravatar imageebs ( 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.

ndomes gravatar imagendomes ( 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

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

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2012-06-12 11:26:42 +0200

Seen: 689 times

Last updated: Jun 16 '12