Сферы с периодическими границами на трехмерной сетке

Я пытаюсь понять, как определить периодические границы на сетке 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

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

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


python numpy mesh periodicity
person Oier Arcelus    schedule 25.05.2020    source источник

Ответы (1)


arrow_upward
1
arrow_downward

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

Ваш код:

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

введите здесь описание изображения

person imochoa    schedule 25.05.2020
comment
Большое спасибо, я придумал свое решение, где я проверяю, выйдет ли сфера из сетки или нет, и если да, я вычислю значения противоположного сайта, пересчитав переменную I, со второй вспомогательной сферой который был перемещен в заданном направлении. Это гораздо более сжато и элегантно. - person Oier Arcelus; 25.05.2020
comment
Вы также можете сделать это здесь, чтобы получить ускорение. Вам просто нужно отфильтровать список centers, чтобы он содержал только центры сфер, которые фактически пересекают основную сетку. Рад, что могу помочь :) Если вы довольны ответом, пометьте его как принятый - person imochoa; 25.05.2020