Dear all,
thank you for providing FiPy -- while we have only just started to use
it and I am not ever sure we understand the basics sufficiently, it is
clearly a powerful tool.
Here is a question to which I couldn't find the answer in the manual,
examples, or by inspecting the variable and mesh classes:
Given a CellVariable that contains the solution of a (scalar) field
phi (through an instance of the CellVariable class), I would like to
know the value of phi at a particular position.
Let's further assume we are looking at a 2d simulation domain (so that
a position would be one point r=(x,y) in the x-y plane).
Conceptually, this would require to find the cell that contains the
position r, and then to retrieve the corresponding value from phi.
Custom code can be written to do this: for a regular mesh (as
fipy.Grid2D would return) this is fairly straightforward, but for an
irregular mesh this is a more expensive operation.
The question is: does FiPy provide any tools that help with this, or
is this something the user needs to build themselves?
I understand that this operation is rarely required and may not be
provided by fipy (in which case this is the anwser to my question),
but I wanted to check first to avoid duplication of effort.
Many thanks,
Hans
PS I attach an example that might make it clearer what I am trying to
do. The question is repeated in the comment block at the of the file.
--
Hans Fangohr
School of Engineering Sciences
University of Southampton
Phone: +44 (0) 238059 8345
Email: [email protected]
http://www.soton.ac.uk/~fangohr
#!/usr/bin/env python
import fipy
Ly = 20
Lx = Ly
nx = 20
ny = nx
dx = 1.
dy = dx
mesh = fipy.Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phi = fipy.CellVariable(name = "solution variable",
mesh = mesh,
value = 0.0)
#Basic equation
D = 1
eq = fipy.TransientTerm() == fipy.DiffusionTerm(coeff=D)
#boundary conditions
phi_outside = 0.1
x, y = mesh.getFaceCenters()
facesTopLeft = ((mesh.getFacesLeft() & (y > Ly / 2))
| (mesh.getFacesTop() & (x < Lx / 2)))
facesBottomRight = ((mesh.getFacesRight() & (y < Ly / 2))
| (mesh.getFacesBottom() & (x > Lx / 2)))
BCs = (fipy.FixedValue(faces=facesTopLeft, value=phi_outside),
fipy.FixedValue(faces=facesBottomRight, value=-phi_outside))
if __name__ == '__main__':
viewer = fipy.Viewer(vars=phi, datamin=-phi_outside*1., datamax=phi_outside*1,colorbar='vertical')
viewer.plot()
#and solve the equation by repeatedly looping in time:
explicit_max_timeStepDuration = 0.9 * dx*dy / (2 * D)
timeStepDuration = 10 * explicit_max_timeStepDuration
steps = 2
for step in range(steps):
eq.solve(var=phi,
boundaryConditions=BCs,
dt=timeStepDuration)
#at this point, we would like to know the value of phi at some position r=(x,y).
#is there a generic function for this, for example
#phi_at_xy = phi.probe(pos=(x,y))
#or maybe
#index = mesh.get_index_from_pos(pos=(x,y))
#phi_at_xy = phi.getValue(index)
if __name__ == '__main__':
viewer.plot()
