Using functions from the math or numpy modules

Dear Fipy users,

After considering many possible solutions for my problem, I am now back 
to FiPy. I could find almost every tool I need in FiPy, but there is at 
least one problem remaining.

The manual of FiPy says that functions such as "exp", etc. are available 
from "fipy.tools.numerix", I've tried that and as the attached script 
shows, it works perfectly. However, if "exp" is taken from the "math" or 
the "numpy" module, FiPy gives an error.

That's a problem for me because I want to use functions such as 
"scipy.integrate.quadrature" and "scipy.special.gamma" in my PDEs.

I noticed, for example, that the output of "math.exp" is a "float" 
object and the output of "fipy.tools.numerix.exp" is a "numpy.float64" 
object. I tried converting so that my function returns this type of 
object, but I haven't had much success with that method.

Is there anything I missed?

Many thanks,

Etienne Rivard

Btw., fans of matlab's pdepe solver should recognize the example in the 
script.

import pdb
#import customio
#reload(customio)
from fipy.tools import numerix
from fipy import *
import math
import numpy

#Filenames
u1 = "u1.dat"
u2 = "u2.dat"

# Reset
if vars().has_key('u1'):
	del u1, u2
if vars().has_key('eq1'):
	del eq1, eq2

# Tolerance
tol = 0.01

# Mesh
length = 1.
nx = 20
dx = 1./nx
mesh = Grid1D(dx=dx, nx=nx)
dt = 0.0021
ti = 0.
tf = 2.
steps = (tf-ti)/dt

# Initial Conditions
# u1(x,0) = 1
# u2(x,0) = 0t
u1 = CellVariable(mesh = mesh, value = 1., hasOld = 1, name = 'u1')
u2 = CellVariable(mesh = mesh, value = 0., hasOld = 1, name = 'u2')

# Boundary Conditions
# du1_dx(0,t) = 0
# u1(1,t) = 1
# u2(0,t) = 0
# du2_dx(1,t) = 0
BC1 = (FixedFlux(faces=mesh.getFacesLeft(), value=0.),
       FixedValue(faces=mesh.getFacesRight(), value=1.))

BC2 = (FixedValue(faces=mesh.getFacesLeft(), value=0.),
       FixedFlux(faces=mesh.getFacesRight(), value=0.))


# Equations
# du1_dt = 0.024 * d2u1_dx2 - F(u1-u2)
# du2_dt = 0.170 * d2u2_dx2 + F(u1-u2)
# F(y) = exp(5.73*y) - exp(-11.46*y)
def F(y):
	return numerix.exp(5.73*y) - numerix.exp(-11.46*y)

eq1 = TransientTerm() == ImplicitDiffusionTerm(coeff=0.024) - F(u1-u2)
eq2 = TransientTerm() == ImplicitDiffusionTerm(coeff=0.170) + F(u1-u2)

t = ti
#customio.vectorscat("u1.dat", t, mesh.getCellCenters()[:,0], u1.value, append=False)
#customio.vectorscat("u2.dat", t, mesh.getCellCenters()[:,0], u2.value, append=False)

while t < tf:
	# Update time
	t = t + dt
	# Calculate
	converged = 0
	while converged == 0:
		res1 = eq1.sweep(var=u1,
	                     boundaryConditions=BC1,
	                     dt=dt)
		res2 = eq2.sweep(var=u2,
	                     boundaryConditions=BC2,
	                     dt=dt)
		resmax = max(res1, res2)
		#print resmax
		if  resmax < tol:
			converged = 1
	
	# Write values to file 
	#customio.vectorscat("u1.dat", t, mesh.getCellCenters()[:,0], u1.value)
	#customio.vectorscat("u2.dat", t, mesh.getCellCenters()[:,0], u2.value)
	timestr = t.__str__() + "/" + tf.__str__()
	print timestr
	u1.updateOld()
	u2.updateOld()