Loading [MathJax]/jax/output/HTML-CSS/jax.js
Ask Your Question
0

Forming a solid from a set of intersecting planes

asked 12 years ago

ebs gravatar image

updated 12 years ago

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

Aix+Biy+Ciz=Di 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!

Preview: (hide)

2 Answers

Sort by » oldest newest most voted
1

answered 12 years ago

ndomes gravatar image

updated 12 years ago

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)
Preview: (hide)
link

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 ( 12 years ago )

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

ndomes gravatar imagendomes ( 12 years ago )
1

answered 12 years ago

ndomes gravatar image

updated 12 years ago

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

Preview: (hide)
link

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 ( 12 years ago )

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 ( 12 years ago )

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 ( 12 years ago )

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: 12 years ago

Seen: 828 times

Last updated: Jun 16 '12