Сферы с границами периода c на 3D-сетке

Question

Я пытаюсь выяснить, как определить границы периода c на numpy мне sh.

Допустим, я определяю поле размером 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

Я получу трехмерный массив логических значений, где значения True - это точки сетки, попадающие внутрь сферы, а значения False - точки сетки, которые попадают внутрь сферы. Однако это не гарантирует границ периодаi c, которые я хотел бы.

Есть ли какой-нибудь элегантный способ для этого, без необходимости l oop для каждой точки сетки

Answer 1

Один простой способ сделать это - воспроизвести ваш круг в соседних сетках периодов c и проверить расстояния от узловых точек в вашей текущей сетке до центров в соседних сетках:

Ваш код :

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,:])

Или полную трехмерную сетку (довольно медленно)

%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()