----- Original Message ----
From: Jonathan Guyer <[EMAIL PROTECTED]>
To: Multiple recipients of list <[email protected]>
Sent: Tuesday, 17 October, 2006 6:09:55 PM
Subject: Re: How to define such cell variables in FiPy?
>Yes, unstructured meshes. "Eighty-two" is a bit of an exaggeration,
>but FiPy is designed throughout to deal with unstructured meshes and
>so it does not know a priori how many neighbors there are. There are
>some efficiency shortcuts for structured meshes, but you should not
>generally assume that you know what the layout of the mesh is.
For structured meshes in the example as attached, I have done a program to
calculate s. Please run the program and then focus on the results displayed on
your screen for "curpos=656 to 671", where the s values should have been
changed to "1"s after line 73 in the program. However, in line 79, when I
check the s value, it is still zero. What is wrong with the program?
Thanks and regards,
hzhatlboro
Send instant messages to your online friends http://uk.messenger.yahoo.com
from fipy.variables.cellVariable import CellVariable
H=0.25
nx=16
ny=90
dx=H
dy=dx
from fipy.meshes.grid2D import Grid2D
mesh=Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phases=[
CellVariable(mesh=mesh, name='phi0', value=0),
CellVariable(mesh=mesh, name='phi1', value=0),
CellVariable(mesh=mesh, name='phi2', value=0),
CellVariable(mesh=mesh, name='phi3', value=0)
]
x, y=mesh.getCellCenters()[...,0], mesh.getCellCenters()[...,1]
phases[0].setValue(1, where=(y<41*H))
phases[1].setValue(1, where=(41*H<y) & (y<45*H) & (x<8*H))
phases[2].setValue(1, where=(41*H<y) & (y<45*H) & (x>8*H))
phases[3].setValue(1, where=(y>45*H))
def scalc(phases):
s=phases
for i in range(1):
s[i]=(phases[i]!=0)
for j in range(2,ny):
for k in range(1,nx+1):
curpos=(j-1)*nx+k-1
curposX=0.5*dx+(k-1)*dx
curposY=0.5*dy+(j-1)*dy
northpos=curpos+nx
northposX=curposX
northposY=curposY+dy
southpos=curpos-nx
southposX=curposX
southposY=curposY-dy
if k==0:
westpos=curpos+nx-1
westposX=curposX+(nx-1)*dx
westposY=curposY
eastpos=curpos+1
eastposX=curposX+dx
eastposY=curposY
elif k==nx-1:
westpos=curpos-1
westposX=curposX-dx
westposY=curposY
eastpos=curpos-nx+1
eastposX=curposX-(nx-1)*dx
eastposY=curposY
else:
westpos=curpos-1
westposX=curposX-dx
westposY=curposY
eastpos=curpos+1
eastposX=curposX+dx
eastposY=curposY
print "S ",i," curpos= ",curpos," S[i][curpos]= ",s[i][curpos]," S[i][north]= ",s[i][northpos]," S[i][south]= ",s[i][southpos]," S[i][west]= ",s[i][westpos]," S[i][east]= ",s[i][eastpos]
if ((s[i][curpos]==0) & (s[i][northpos]==0) & (s[i][southpos]==0) & (s[i][westpos]==0) & (s[i][eastpos]==0)):
#s[i].setValue(0, where=(x==curposX) & (y==curposY))
print "S ",i," value at ",curpos," does not need to change \n"
else:
if s[i][curpos]==0:
s[i].setValue(1, where=((x==curposX) & (y==curposY)))
print "\nS ",i," value at ",curpos,"is changed to 1"
print "x coordiate at ",curpos," is ",x[curpos]
print "curposX = ",curposX
print "y coordiate at ",curpos," is ",y[curpos]
print "curposY = ",curposY
print "Now S ",i," value at ",curpos,"is ",s[i][curpos],"\n"
return s
s=scalc(phases)
