# 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 <[EMAIL PROTECTED]> 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: [EMAIL PROTECTED]
GPG Key ID# 0x3AE05130
_______________________________________________
Numpy-discussion mailing list
[email protected]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________
Numpy-discussion mailing list
[email protected]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
