python - 3D 网格上具有周期性边界的球体

Question

我试图弄清楚如何在 numpy 网格上定义周期性边界。

假设我定义了一个大小为 1x1x1 的盒子，我在里面放了一个半径为 0.25 的球体。这个球体不在中心，但离边界足够近，以至于球体的一部分必须从盒子的另一侧出来。

例如，如果以下代码

  import numpy as np

  x_ = np.linspace(0,1,100)
  y_ = np.linspace(0,1,100)
  z_ = np.linspace(0,1,100)

  X,Y,Z = np.meshgrid(x_,y_,z_,indexing='ij')

  I = (X-particle['x'])**2 + (Y-particle['y'])**2 + (Z-particle['z'])**2 < particle['r']**2

我将得到一个 3D 布尔数组，其中 True 值是落在球体内的网格点，而 False 值是落在球体内的网格点。然而，这并不能保证我想要的周期性边界。

有什么优雅的方法，而不必遍历每个网格点

Answer 1

一种简单的方法是在相邻的周期性网格中复制你的圆圈，并检查从当前网格中的网格点到相邻网格中的中心的距离：

您的代码：

import numpy as np

x_ = np.linspace(0,1,100)
y_ = np.linspace(0,1,100)
z_ = np.linspace(0,1,100)

X,Y,Z = np.meshgrid(x_,y_,z_,indexing='ij')

我添加了一些示例圆参数：

particle = {'r':0.25, 'x':0.3, 'y':0.5,'z':0.8}

由于您的网格长度为 1x1x1，我猜点之间的间距为 0.01，因此：

import itertools
grid_size = 1.0
offsets = itertools.combinations_with_replacement([grid_size,0,-grid_size],r=3)
centers = [(particle['x']+x_offset, particle['y']+y_offset,particle['z']+z_offset) for x_offset,y_offset, z_offset in offsets]
I=np.logical_or.reduce([(X-c_x)**2 + (Y-c_y)**2 + (Z-c_z)**2 < particle['r']**2 for c_x, c_y, c_z in centers])

您可以通过可视化切片来仔细检查它：

from matplotlib import pyplot as plt
plt.imshow(I[:,50,:])

或完整的 3D 网格（很慢）

%matplotlib notebook
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.voxels(I)
plt.show()