3

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で説明されている方法を使用しようとしました

CellVariableのcallメソッドは、 callメソッドを介して渡された一連の座標への補間を可能にします( callメソッドは、関数呼び出しのように括弧を使用してアクセスされます)。渡された各座標に対応する一連の値を返します。順序引数は、補間の順序を決定するだけです。

これが実際にどのように機能するかはわかりませんが、私が理解していることから、これは単一の平面を 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.]

どこかに愚かな間違いがあるに違いないと思いますが、手がかりがありません。何か案が?

4

1 に答える 1

0

セルの中心は 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))
于 2015-05-02T00:21:25.147 に答える