Как извлечь плоскость из 3D-переменной в FiPy (из 3D в 2D)

У меня есть переменная на 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.]

Я думаю, что где-то должна быть глупая ошибка, но я понятия не имею. Любая идея?


person Community    schedule 25.04.2015    source источник
comment
Для всех, кто интересуется, что такое fipy, это может быть следующее: github.com/usnistgov/fipy. Поскольку это кажется довольно нишевым пакетом, и ваш вопрос, похоже, касается того, как этот пакет работает, вам может повезти больше в списке рассылки: ctcms.nist.gov/fipy/documentation/MAIL.html. Вы можете включить ссылку на этот вопрос, если напишете в список рассылки, и вы можете опубликовать решение здесь, если получите полезный совет в другом месте.   -  person YXD    schedule 25.04.2015
comment
Забавно, на StackOverflow уже было 17 вопросов с пометкой FiPy, но это первый вопрос, в котором кто-то добавляет ссылку на проект. Собственно (если кому-то интересно) главная страница проекта - ctcms.nist.gov/ fipy / index.html со всей документацией и инструкциями по установке, а на GitHub есть только исходный код. Сказал, что, вероятно, это хорошая идея попробовать со списком рассылки, но я все еще надеюсь, что другой случайный пользователь FiPy может прочитать этот вопрос и дать мне объяснение этого странного поведения или хорошее обходное решение.   -  person    schedule 25.04.2015
comment
Мы отслеживаем как список рассылки, так и StackOverflow.   -  person jeguyer    schedule 02.05.2015


Ответы (1)


Центры ячеек расположены в точках z = 0,5, 1,5, 2,5, ... поэтому FiPy старается найти ближайшие ячейки к z = 0, 1, 2, ...

Пытаться

var2D.setValue(var3D((x2D, y2D, 
                      zcut * ones(var2D.mesh.numberOfCells) + dx/2.), order=0))
person jeguyer    schedule 02.05.2015
comment
Это работает, но не могли бы вы объяснить мне, почему в половине случаев используется левое значение для вычисления значения, а в половине случаев используется правильное значение? Расстояние между левым центром или правым центром всегда одинаково для каждой ячейки (по крайней мере, в этом простом случае). Я не мог понять причину. - person ; 02.05.2015
comment
Подпрограмма FiPy, которая вычисляет ближайший к точке идентификатор ячейки, отображается по адресу github.com/usnistgov/fipy/blob/develop/fipy/meshes/, использует rint (). Возможно, более интуитивный результат был бы получен с помощью fix (), но в любом случае это просто соглашение. - person jeguyer; 05.05.2015
comment
Спасибо за объяснение. - person ; 05.05.2015
comment
На самом деле, fix () будет хуже, потому что округляется до 2,9. - person jeguyer; 06.05.2015