3

我有一个有限元程序的结果，该程序在三维空间中的规则间隔网格位置给出了各种感兴趣的测量值（例如，温度、密度、压力）。

值沿每个坐标等距分布，但不同坐标的间距可能不同。例如，

x1 = [0, 0.1, 0.2, ..., 1.0]      (a total of NX1 pts) 
x2 = [0, 0.5, 1.0, ..., 20]       (a total of NX2 pts) 
x3 = [0, 0.2, 0.4, ..., 15]       (a total of NX3 pts)

软件输出的结果格式如下：

x1_1, x2_1, x3_1, f_x, g_x, h_x
x1_1, x2_1, x3_2, f_x, g_x, h_x
x1_1, x2_1, x3_3, f_x, g_x, h_x
...
x1_1, x2_2, x3_1, f_x, g_x, h_x
x1_1, x2_2, x3_2, f_x, g_x, h_x
x1_1, x2_2, x3_3, f_x, g_x, h_x
...
x1_2, x2_1, x3_1, f_x, g_x, h_x
x1_2, x2_1, x3_2, f_x, g_x, h_x
x1_2, x2_1, x3_3, f_x, g_x, h_x
...

其中 f_x、g_x、h_x 是在特定网格点处感兴趣的度量。

我想转换上述数据格式并获得 f、g 和 h 的 (NX1 x NX2 x NX3) numpy 数组。

一些结果集相当大（80 x 120 x 100）。

有没有人有任何提示可以有效地进行这种转换？

4

2 回答 2

1

假设您将整个数组作为datashape数组读入内存(Nx1 * Nx2 * Nx3, 6)。

data = np.loadtxt('data.txt', dtype=float, delimiter=',')

如果如您的示例所示，这些点是按字典顺序生成的，则只需将列抓取到f，g然后重新调整h它们的形状：

f = data[:, 3].reshape(Nx1, Nx2, Nx3)
g = data[:, 4].reshape(Nx1, Nx2, Nx3)
h = data[:, 5].reshape(Nx1, Nx2, Nx3)

如果您需要弄清楚和是什么Nx1，您可以使用：Nx2Nx3np.unique

Nx1 = np.unique(data[:, 0]).shape[0]
Nx2 = np.unique(data[:, 1]).shape[0]
Nx3 = np.unique(data[:, 2]).shape[0]

如果不能保证点的顺序，一个更强大的解决方案是np.unique用于提取网格值的索引：

Nx1, idx1 = np.unique(data[:, 0], return_inverse=True)
Nx1 = Nx1.shape[0]
Nx2, idx2 = np.unique(data[:, 1], return_inverse=True)
Nx2 = Nx2.shape[1]
Nx3, idx3 = np.unique(data[:, 2], return_inverse=True)
Nx3 = Nx3.shape[0]

f = np.empty((Nx1, Nx2, Nx3))
f[idx1, idx2, idx3] = data[:, 3]
g = np.empty((Nx1, Nx2, Nx3))
g[idx1, idx2, idx3] = data[:, 4]
h = np.empty((Nx1, Nx2, Nx3))
h[idx1, idx2, idx3] = data[:, 5]

这将为f,g和h, 而不是原始数组的视图创建新数组data，因此它将使用更多内存。

当然，您应该使用循环或列表推导，而不是我上面的丑陋代码重复三遍！

于 2013-02-27T15:56:27.603 回答
0

无论如何，您都必须遍历整个文件，那么为什么不初始化数组并输入值呢？

棘手的部分是，如果你想要一个 shape 的数组(NX1,NX2,NX3)，但如果你的x1,x2,x3值是floats，那么你必须以某种方式索引你的数组。也许存在一个数据结构，但你可以使用类似的东西

def xyz_index((x,y,z),(n1,n2,n3)):
    """ return integer indices for x,y,z position
        given a constant step """
    return tuple(map(int,[x/n1,y/n2,z/n3]))

import numpy as np
NX1,NX2,NX3 =  (80, 120, 100)
ns = n1, n2, n3 =   (.1,.5,.2)
x1, x2, x3 = np.arange(0,1+n1,n1), np.arange(0,20+n2,n2), np.arange(0,15+n3,n3),

data = np.empty((NX1,NX2,NX3),dtype=[('f',float),('g',float),('h',float)])
with open(filename,'r') as f:
    for line in f:
        x,y,z,f,g,h = map(float,line.split(', '))
        data[xyz_index((x,y,z),ns)] = (f,g,h)

然后您可以按如下方式访问您的数据：

对于h点的 - 值x1,x2,x3 = .2, .5, 0.，使用

data[xyz_index((.2,.5,0),ns)]['h']

如果没有['h']，这将返回一个(f,g,h)包含上述内容的元组dtype。

如果没有索引，它将返回一个(NX1,NX2,NX3)包含所有值的数组h。

现在我看了一下，如果n1, n2, n3总是相同的，你可能想在你的xyz_index函数中定义它们，这样你就不必ns每次都传入：

data[xyz_index(.2,.5,0)]['h']
于 2013-02-27T15:33:41.290 回答