【问题标题】:Generate a mesh on a cylindrical surface and bin data to it在圆柱面上生成网格并将数据分箱
【发布时间】:2017-10-05 06:31:12
【问题描述】:

我不确定如何在 Python 中解决以下问题:
我有一个圆柱面,我想mesh 然后bin 数据到网格的每个块。我有撞击圆柱表面的光线的位置数据(X、Y、Z)。在圆柱体上生成网格后,我必须计算网格中每个块中的光线(数据点)数量。

气缸参数如下:
半径 = 0.1m,高度 = 0.15m,中点:[X=0,Y=0,Z=0]

要生成网格,可以使用来自How to generate regular points on cylindrical surface 的答案以及下面复制的代码: `

import numpy as np
def make_cylinder(radius, length, nlength, alpha, nalpha, center, orientation):
"""
radius = radius of cylinder,
length = cylinder height, nlength = number of lengthwise divisions.
alpha = total degrees of cylinder, nalpha = number of circumferential divisions.
center = [X,Y,Z] coordinates of cylinder's midpoint.

"""
#Create the length array
I = np.linspace(0, length, nlength)
#Create alpha array avoid duplication of endpoints
#Conditional should be changed to meet your requirements
if int(alpha) == 360:
    A = np.linspace(0, alpha, num=nalpha, endpoint=False)/180*np.pi
else:
    A = np.linspace(0, alpha, num=nalpha)/180*np.pi

#Calculate X and Y
X = radius * np.cos(A)
Y = radius * np.sin(A)

#Tile/repeat indices so all unique pairs are present
pz = np.tile(I, nalpha)
px = np.repeat(X, nlength)
py = np.repeat(Y, nlength)

points = np.vstack(( pz, px, py )).T

#Shift to center
shift = np.array(center) - np.mean(points, axis=0)
points += shift

#Orient tube to new vector
#Grabbed from an old unutbu answer
def rotation_matrix(axis,theta):
    a = np.cos(theta/2)
    b,c,d = -axis*np.sin(theta/2)
    return np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],
                     [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],
                     [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]])

ovec = orientation / np.linalg.norm(orientation)
cylvec = np.array([1,0,0])

if np.allclose(cylvec, ovec):
    return points

#Get orthogonal axis and rotation
oaxis = np.cross(ovec, cylvec)
rot = np.arccos(np.dot(ovec, cylvec))

R = rotation_matrix(oaxis, rot)
return points.dot(R)

` 现在剩下的就是循环遍历“光线数据”以找到击中每个网格块的光线数量。

我最初的思考过程如下:

data = np.genfromtxt('ray_data.csv', delimiter=',') 
num_data_rows, num_data_cols = np.shape(data)

for i in range (num_data_rows): #Loop through the data

这就是我卡住的地方。如前所述,“射线数据”是一个 CSV 文件,其中包含撞击圆柱表面的每条射线的位置 (X,Y,Z)。请参阅提供的链接:Ray data sample for cylindrical surface

我只需要弄清楚如何检查光线落在我的网格中的位置。 每个块中的光线数量将乘以一个常数(每条光线的功率)以获得每个块中的功率(以瓦特为单位)。然后将该值除以块的面积,得到热通量(W/m^2)。

我需要的最终输出是一个数组,其中包含网格的每个块的质心以及相应的热通量值。

任何想法如何解决这个问题? 我相信与pandas 合作也是一种选择。

【问题讨论】:

  • 您只对表面感兴趣,所以我建议将您的数据和圆柱表面转换为二维空间
  • 将数据转换为 2D 空间仍然会给我留下将光线分箱到网格中并确定每个箱的质心的问题。
  • 你的“光线”只是点。我是否更正了“射线”实际上是指某些射线(我认为是矢量)与圆柱体表面的交点?
  • 另外,如果圆柱体的直径远大于您在长轴和短轴上的网格间距,那么我会 1)找到每个箱的质心(或简单地重新解释您的箱角作为移动网格的质心)2)找到光线与表面(您可能已经拥有)的交点,3)通过找到每个交点的最近邻质心来确定相应的箱。最后一步可以使用 KD 树(例如scipy.spatial.cKDTree)有效地计算。

标签: python pandas numpy binning


【解决方案1】:

正如评论中已经提出的,我建议将您的表面转换为 2d 空间。这样,您就可以轻松地对数据进行分组。

import numpy as np
import pandas as pd

# generate some random rays (for which we will just assume they hit the surface)

rays = pd.DataFrame(
  np.random.uniform(-1,1,(8000,4)),
  columns=["x", "y", "z", "intensity"]
)

# transform x and y to polar coordinates while dropping the radius
# (again, just assuming they hit the surface)

rays["phi"] = rays.T.apply(lambda row: np.arctan2(row.y, row.x)).T

rays.head() 现在看起来如下所示。 zphi 是 2d 中光线的表示。 intensity 就是你所说的“射线的力量”。

x         y         z  intensity       phi
0 -0.237026 -0.634709 -0.889694   0.362156 -1.928199
1 -0.481137 -0.446912  0.687224   0.268080 -2.393056
2 -0.805538  0.068678  0.272009   0.990947  3.056541
3  0.549282 -0.330665  0.318683  -0.150776 -0.541886
4 -0.215676 -0.030922 -0.478929   0.408720 -2.999190

现在,只需创建 bin 并对数据进行分组。最后总结所有强度。

z_bins = np.arange(0, 1, .1)
phi_bins = np.arange(-np.pi, np.pi, np.pi/10)

result = rays.groupby([
  pd.cut(rays.phi, phi_bins), 
  pd.cut(rays.z, z_bins), 
]).intensity.sum()

result.head() 如下所示:

phi               z         
(-3.142, -2.827]  (0, 0.1]      0.719154
                  (0.1, 0.2]   -1.733479
                  (0.2, 0.3]    2.073013
                  (0.3, 0.4]    1.967453
                  (0.4, 0.5]    0.001312

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2011-02-10
    • 2021-05-14
    • 2020-09-06
    • 1970-01-01
    • 1970-01-01
    • 2014-12-12
    • 2021-07-05
    • 1970-01-01
    相关资源
    最近更新 更多