3

I have a variable on a 3D mesh and I am trying to cut a plan. I am surprised this question hasn't been asked before, it looks an easy and common problem but I haven't found any good way. I would appreciate any advice.

Let's say that I have a parallelepiped 3x3x5 and that I am trying to extract a z-plane.

from fipy import *
from numpy import *
#Describes a 3x3x5 mesh
nx = 3
ny = 3
nz = 5

dx = 1
dy = 1
dz = 1

#Creates a 3D mesh and a 2D mesh to store the plane
mesh3D = Grid3D(dx, dy, dz, nx, ny, nz)
mesh2D = Grid2D(dx, dy, nx, ny)
#Defines the same variable, in 3D and in 2D
var2D = CellVariable(mesh = mesh2D)
var3D = CellVariable(mesh = mesh3D)

#Fills the 3D matrix
for i in xrange(nx*ny*nz*dx*dy*dz):
    var3D[i] = i

print var3D

Output:

[  0.   1.   2.   3.   4.   5.   6.   7.   8.   9.  10.  11.  12.  13.  14.  15.  16.  17.  18.  19.  20.  21.  22.  23.  24.  25.  26.  27.  28.  29.  30.  31.  32.  33.  34.  35.  36.  37.  38.  39.  40.  41.  42.  43.  44.]

The 3D variable looks correctly filled.

First, I tried to use the way described at this link http://permalink.gmane.org/gmane.comp.python.fipy/1576

The call method of a CellVariable enables interpolation to a set of coordinates passed in via the call method (the call method is just accessed using parentheses like in a function call). It returns a set of values for corresponding to each of the coordinates that were passed in. The order argument just determines the order of the interpolation.

I am not sure of how this actual works, but from what I understand this should interpolate a single plane with 0 order, so it should extract the exact value on a specific plane. Please correct me if I am wrong.

x3D, y3D, z3D = mesh3D.getCellCenters()
x2D, y2D =  mesh2D.getCellCenters()

for zcut in xrange(nz*dz):
    var2D.setValue(var3D((x2D, y2D, zcut * ones(var2D.getMesh().getNumberOfCells())), order=0))
    print "z-plane = %d" % zcut
    print var2D

raw_input("Press any key to close")

Strangely it doesn't work. The even indexes are ok but the odd ones are a copy of the contiguous planes.

z-plane = 0
[ 0.  1.  2.  3.  4.  5.  6.  7.  8.]
z-plane = 1
[ 0.  1.  2.  3.  4.  5.  6.  7.  8.]
z-plane = 2
[ 18.  19.  20.  21.  22.  23.  24.  25.  26.]
z-plane = 3
[ 18.  19.  20.  21.  22.  23.  24.  25.  26.]
z-plane = 4
[ 36.  37.  38.  39.  40.  41.  42.  43.  44.]

I think there must be a stupid mistake somewhere, but I have no clue. Any idea?

YXD
  • 31,741
  • 15
  • 75
  • 115
  • For anyone else wondering what `fipy` is, it seems to be this: https://github.com/usnistgov/fipy . Since this appears to be a rather niche package and your question seems to be about how that package works, you may have more luck on the mailing list: http://www.ctcms.nist.gov/fipy/documentation/MAIL.html . You can include a link to this question if you write to the mailing list and you can post a solution here if you get helpful advice elsewhere. – YXD Apr 25 '15 at 11:11
  • Funny, there were already 17 questions tagged as FiPy on StackOverflow but this is the first one in which someone adds reference to the project. Actually (if someone is interested) the main page of the project is http://www.ctcms.nist.gov/fipy/index.html with all the documentation and instructions for the installation, while on GitHub there is only the source code. Said that, probably it is a good idea to try with the mailing list, but I still hope that another "casual FiPy user" may read this question and give me an explanation of this strange behaviour or a good workaround. –  Apr 25 '15 at 13:54
  • We do monitor both the mailing list and StackOverflow. – jeguyer May 02 '15 at 00:22

1 Answers1

0

The cell centers are located at z=0.5, 1.5, 2.5, ... so FiPy is doing it's best to find the nearest cells to z=0, 1, 2, ...

Try

var2D.setValue(var3D((x2D, y2D, 
                      zcut * ones(var2D.mesh.numberOfCells) + dx/2.), order=0))
jeguyer
  • 2,379
  • 1
  • 11
  • 15
  • That works, but could you please explain me why half of the time it uses the left value to calculate the value and half of the time it uses the right value? The distance between the left center or the right center is always the same for every cell (at least in this simple case). I couldn't get the reason. –  May 02 '15 at 17:27
  • The FiPy routine that calculates the nearest cell ID to a point, seen at https://github.com/usnistgov/fipy/blob/develop/fipy/meshes/uniformGrid3D.py#L530, uses NumPy's [rint()](http://docs.scipy.org/doc/numpy/reference/generated/numpy.rint.html). An arguably more intuitive result would be obtained from [fix()](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fix.html), but it's just a convention, either way. – jeguyer May 05 '15 at 14:51
  • Actually, [fix()](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fix.html) would be worse, because it rounds 2.9 to 2. – jeguyer May 06 '15 at 17:10