[Numpy-discussion] First Numerical Python code...comments?

[Torgil Svensson]

> # Fill in SK with proper values
> for i in range(numel):
>   RL = xord[i+1] - xord[i]
>   RK = tens/RL
>   SK[i][0] = SK[i][0] + RK
>   SK[i][1] = SK[i][1] - RK
>   IP1 = i+1
>   SK[IP1][0] = SK[IP1][0] + RK

To get advantage of numpy you'll need to calculate with arrays instead of
individual numbers:

RL=xord[1:]-xord[:-1]
RK=tens/RL
SK[:-1,0]+=RK
SK[:-1,1]-=RK
SK[1:,0]+=RK

hope I got that right... but you'll get the point anyway

//Torgil

On 3/4/07, Tyler Hayes <hayes.tyler at gmail.com> wrote:
>
> Hello All:
>
> Well, I recently made the switch to Python for scientific computing
> (using NumPy and Matplotlib/PyLab) and got my first program to work.
> It is a simple 1D finite element code for a tightly stretched wire and
> is a textbook example.
>
> To solve it, I also created a symmetric, banded matrix solver module
> from a Fortran90 routine using f2py (thank you Pearu for your help the
> other day to get f2py working).
>
> Anyways, my background is in Fortran and I was got the program to run,
> but I can't help but wonder if I did not use the most efficient method
> for calcualting the matrix values, etc.
>
> I would appreciate your comments on my code that I have attached
> below. My concern is that I used some "for" loops as I would have
> using Fortran, but perhaps there is a more intelligent (i.e.,
> efficient) way using native NumPy matirx routines?
>
> I would love to hear what you think.
>
> Also, while I'm at it. I have a simple question about importing NumPy,
> SciPy, and PyLab. If I import PyLab, it says I also import NumPy. Do I
> ever need to call in SciPy as well? Or is this the same thing?
>
> Cheers,
>
> t.
>
> The "code"
>
> #!/usr/bin/env python
> #
> # This is the example code from Thompson for a tight wire problem
> #
> # Coded by Tyler Hayes March 1, 2007
> #
> #################################################################
>
> # Load the pylab module for plotting (it also imports NumPy)
> from pylab import *
> # Load the sparse matrix solver symgauss created with f2py
> import symgauss as G
> """
> This is the example program from Thompson which will plot the
> displacement of a tight wire with constant loading.
>
> DEPENDS:
> This code requires the symgauss module to solve the sparse
> symmetric matrix.
> """
>
> # Initialize data
> #----------------------------------------------------------------
>
> # Parameters
> tens  = 20000.0
> numel = 8
> numnp = numel + 1
> IB    = 2
> BLAST = 1.0e6
>
> # Boundary conditions
> npbc     = zeros(numnp,float)
> npbc[0]  = 1.0
> npbc[-1] = 1.0
>
> # Displacement
> U = [0.0]*numnp
>
> # Element X-coordinates
> xord     = zeros(numnp,float)
> xord[1]  = 2.0
> xord[2]  = 4.0
> xord[3]  = 5.0
> xord[4]  = 6.0
> xord[5]  = 7.0
> xord[6]  = 8.0
> xord[7]  = 10.0
> xord[8]  = 12.0
>
> # Force data vector
> F    = -1.0*ones(numnp,float)
> F[1] = -2.0
> F[2] = -1.5
> F[-2]= -2.0
> F[-3]= -1.5
>
> # Initialize SK --- the sparse symmetrix stiffness matrix "K"
> SK = zeros((numnp,IB),float)
>
>
> # Main routine
> #----------------------------------------------------------------
>
> # Fill in SK with proper values
> for i in range(numel):
>    RL = xord[i+1] - xord[i]
>    RK = tens/RL
>    SK[i][0] = SK[i][0] + RK
>    SK[i][1] = SK[i][1] - RK
>    IP1 = i+1
>    SK[IP1][0] = SK[IP1][0] + RK
>
> # Set boundary conditions
> for i in range(numnp):
>    if npbc[i] == 1:
>        SK[i][0] = SK[i][0]*BLAST
>        F[i] = U[i]*SK[i][0]
>
> # Call Symmetric Sparse matrix solver
> U = G.symgauss(SK,F,numnp,IB)
>
>
> # Visualize the data
> #----------------------------------------------------------------
>
> # Plot data using matplotlib
> plot(xord,U)
> title('Deflection of a tighly stretched wire')
> xlabel('Distance along wire (normalized)')
> ylabel('Deflection of wire (normalized)')
> grid(True)
> show()
>
>
> --
> Tyler Joseph Hayes
> 600 Talbot St. -- Apt. 812
> London, Ontario
> N6A 5L9
>
> Tel   : 519.435.0967
> Fax  : 519.661.3198
> Cell  : 416.655.7897
> email: thayes at uwo.ca
>
> GPG Key ID# 0x3AE05130
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20070304/1584cea0/attachment.html>