Generating 3D graphic by PyOpenGL

Question

I'm using PyOpenGL to generate the 3D sea surface according to "2D wave equation". The main purpose is to show the dynamic graphic of "2D wave equation".But it keeps telling me this error:

E:\WinPython-64bit-3.4.3.5\python-3.4.3.amd64\mytest>python seawave_2d_opengl.py
Traceback (most recent call last):
  File "E:\WinPython-64bit-3.4.3.5\python-3.4.3.amd64\lib\site-packages\OpenGL\GLUT\
special.py", line 130, in safeCall
    return function( *args, **named )
  File "seawave_2d_opengl.py", line 106, in Draw
    glVertex3f(x,y,z)
  File "E:\WinPython-64bit-3.4.3.5\python-3.4.3.amd64\lib\site-packages\OpenGL\p
latform\baseplatform.py", line 402, in __call__
    return self( *args, **named )
ctypes.ArgumentError: argument 3: <class 'TypeError'>: wrong type
GLUT Display callback <function Draw at 0x0000000004086950> with (),{} failed: r
eturning None argument 3: <class 'TypeError'>: wrong type

E:\WinPython-64bit-3.4.3.5\python-3.4.3.amd64\mytest>

Here is the code:

from numpy import linspace,zeros,sin,pi,exp,sqrt
from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *
import sys

def solver0(I, f, c, bc, Lx, Ly, nx, ny, dt, tstop, user_action=None):
    dx = Lx/float(nx)
    dy = Ly/float(ny)
    x = linspace(0, Lx, nx+1)  #grid points in x dir
    y = linspace(0, Ly, ny+1)  #grid points in y dir
    if dt <= 0:                #max time step?
        dt = (1/float(c))*(1/sqrt(1/dx**2 + 1/dy**2))
    Cx2 = (c*dt/dx)**2
    Cy2 = (c*dt/dy)**2  #help variables
    dt2 = dt**2

    up = zeros((nx+1,ny+1))  #solution array
    u = up.copy()            #solution at t-dt
    um = up.copy()           #solution at t-2*dt

    #set initial condition:
    t =0.0
    for i in range(0,nx):
        for j in range(0,ny):
            u[i,j] = I(x[i], y[j])
    for i in range(1,nx-1):
        for j in range(1,ny-1):
            um[i,j] = u[i,j] + \
                      0.5*Cx2*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + \
                      0.5*Cy2*(u[i,j-1] - 2*u[i,j] + u[i,j+1]) + \
                      dt2*f(x[i], y[j], t)
    #boundary values of um (equals t=dt when du/dt=0)
    i = 0
    for j in range(0,ny): um[i,j] = bc(x[i], y[j], t+dt)
    j = 0
    for i in range(0,nx): um[i,j] = bc(x[i], y[j], t+dt)
    i = nx
    for j in range(0,ny): um[i,j] = bc(x[i], y[j], t+dt)
    j = ny
    for i in range(0,nx): um[i,j] = bc(x[i], y[j], t+dt)

    if user_action is not None:
        user_action(u, x, y, t)   #allow user to plot etc.

    while t <= tstop:
        t_old = t
        t += dt

        #update all inner points:
        for i in range(1,nx-1):
            for j in range(1,ny-1):
                up[i,j] = -um[i,j] + 2*u[i,j] + \
                          Cx2*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + \
                          Cy2*(u[i,j-1] - 2*u[i,j] + u[i,j+1]) + \
                          dt2*f(x[i], y[j], t_old)

        #insert boundary conditions:
        i = 0
        for j in range(0,ny): up[i,j] = bc(x[i], y[j], t)
        j = 0
        for i in range(0,nx): up[i,j] = bc(x[i], y[j], t)
        i = nx
        for j in range(0,ny): up[i,j] = bc(x[i], y[j], t)
        j = ny
        for i in range(0,nx): up[i,j] = bc(x[i], y[j], t)

        if user_action is not None:
            user_action(up, x, y, t)

        um, u, up = u, up, um  #update data structures
    return u  #dt might be computed in this function
#Actually,the book wrote `return dt`,but I changed `dt` to `u`
def I(x, y):
    return exp(-(x-Lx/2.0)**2/2.0 -(y-Ly/2.0)**2/2.0)
def f(x, y, t):
    return sin(2*x*y*pi*t/Lx)  #defined by myself
def bc(x, y, t):
    return 0.0
#These three functions are some basic functions related to the first function "solver0"

Lx = 10
Ly = 10
c = 1.0
dt = 0
nx = 40
ny = 40
tstop = 20

#The following part is to generate 3D graphics,where I must make mistakes:
def init():
    glClearColor(1.0,1.0,1.0,0.0)  

def Draw():
    glClear(GL_COLOR_BUFFER_BIT)
    glColor3f(0,0,1.0)
    glBegin(GL_LINES)   
    for t in range(0,20,1):
        z = solver0(I, f, c, bc, Lx, Ly, nx, ny, dt, t, user_action=None)
        glVertex3f(x,y,z) 
 #x and y cannot be used here because they are not defined as global variables.
    glEnd()
    glFlush  

def Update():
    global t
    t += 0.1
    glutPostRedisplay() 

def main():    
    glutInit(sys.argv)
    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGBA)
    glutInitWindowSize(800,600)
    glutInitWindowPosition(100,50)
    glutCreateWindow("2D Wave Equations".encode("cp932"))
    init()
    glutDisplayFunc(Draw)
    glutIdleFunc(Update)
    glutMainLoop()   

main()

What mistake did I made? Can someone help me of this? :(

Answer 1

As you can guess from the stacktrace, one of the arguments (the third?!?) in glVertex3f(x,y,z) has a wrong type. The discussion in the comments made clear that z is a two dimensional ndarray while glVertex3f() expects scalars. It looks like solver0() computes an array of z values instead of one z-value per call.

EDIT I now sort of understand what solver0() does. The function should be documented in the book where it was printed in. Although Stackoverflow is not meant to give interpretations of copy and paste code, I'll give a little overview of what I think the function does:

1. Lx and Ly give the range of all x and y
2. nx and ny give the number of x and y value between 0 and Lx, Ly that are used.
3. The function computes an array of x and y values at wich the z-value (up) is computed.
4. It computes up for several time values from 0 to tstop with a step width dt.
5. If a user function user_action is given, it is called after up was computed. The userfunction is called with up, x, y, t as arguments.

To sum things up: One call of solver0 computes all x, y, and z values for a given range of x and y value and a given time span with a given resolution.