У меня есть переменная на 3D-сетке, и я пытаюсь вырезать план. Я удивлен, что этот вопрос не задавали раньше, это выглядит простой и распространенной проблемой, но я не нашел хорошего решения. Буду признателен за любой совет.
Допустим, у меня есть параллелепипед 3x3x5 и я пытаюсь извлечь z-плоскость.
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
Вывод:
[ 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.]
3D-переменная выглядит правильно заполненной.
Сначала я попытался использовать способ, описанный по этой ссылке http://permalink.gmane.org/gmane.comp.python.fipy/1576
Метод call CellVariable включает интерполяцию к набору координат, переданному через метод call (доступ к методу call просто осуществляется с помощью круглых скобок, например в вызове функции). Он возвращает набор значений, соответствующих каждой из переданных координат. Аргумент order просто определяет порядок интерполяции.
Я не уверен, как это работает на самом деле, но, насколько я понимаю, это должно интерполировать одну плоскость с порядком 0, поэтому оно должно извлекать точное значение на определенной плоскости. Пожалуйста, поправьте меня, если я ошибаюсь.
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")
Как ни странно не работает. Четные индексы допустимы, но нечетные - это копии смежных плоскостей.
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.]
Я думаю, что где-то должна быть глупая ошибка, но я понятия не имею. Любая идея?
fipy
, это может быть следующее: github.com/usnistgov/fipy. Поскольку это кажется довольно нишевым пакетом, и ваш вопрос, похоже, касается того, как этот пакет работает, вам может повезти больше в списке рассылки: ctcms.nist.gov/fipy/documentation/MAIL.html. Вы можете включить ссылку на этот вопрос, если напишете в список рассылки, и вы можете опубликовать решение здесь, если получите полезный совет в другом месте. - person YXD   schedule 25.04.2015