Spheres with periodic boundaries on a 3D mesh

Question

I am trying to figure out how to define periodic boundaries on a numpy mesh.

Let's say I define a box of size 1x1x1, and I put a sphere of radius 0.25 inside. This sphere is not in the center but close enough to a border such that part of the sphere has to come out from the opposite side of the box.

For instance if code of the following

  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

I will get a 3D array of booleans where True values are the meshpoints that fall inside the sphere and False values are the meshpoints that fall inside the sphere. However this does not guarantee the periodic boundaries that I would like.

Is there any elegant way for this, without having to loop over each gridpoint

Simon · Accepted Answer · 2022-07-08T09:30:11.727

Sorry, but the answer provided by imochoa is not fully correct. The offset does not include all 26 periodic neighbor boxes but only 9:

import itertools
grid_size = 1.0
offsets = itertools.combinations_with_replacement([grid_size,0,-grid_size],r=3)
[(x_offset, y_offset, z_offset) for x_offset,y_offset, z_offset in offsets]
>>>[(1.0, 1.0, 1.0),
 (1.0, 1.0, 0),
 (1.0, 1.0, -1.0),
 (1.0, 0, 0),
 (1.0, 0, -1.0),
 (1.0, -1.0, -1.0),
 (0, 0, 0),
 (0, 0, -1.0),
 (0, -1.0, -1.0),
 (-1.0, -1.0, -1.0)]

This is missing offsets. e.g: (1.0,0.0,1.0).

I suggest the following approach to create the offsets:

off_x,off_y,off_z= np.meshgrid([-1,0,1],[-1,0,1],[-1,0,1],indexing='ij')
offs=np.array([off_x.flatten(),off_y.flatten(),off_z.flatten()]).T;
centers_fix = [(particle['x']+x_offset, particle['y']+y_offset,particle['z']+z_offset) for x_offset,y_offset, z_offset in offs]
I_fix=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_fix])

Comparison to the previous answer:

import numpy as np
import itertools
grid_size = 1.0

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.9, 'y':0.9,'z':0.9}

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

off_x,off_y,off_z= np.meshgrid([-1,0,1],[-1,0,1],[-1,0,1],indexing='ij')
offs=np.array([off_x.flatten(),off_y.flatten(),off_z.flatten()]).T;
centers_fix = [(particle['x']+x_offset, particle['y']+y_offset,particle['z']+z_offset) for x_offset,y_offset, z_offset in offs]
I_fix=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_fix])

from matplotlib import pyplot as plt
fig, ax = plt.subplots(3, 2,figsize=[10,10])
ax[0, 0].imshow(I[int(particle['x']*100),:,:], origin='lower')
ax[0, 0].set_title('I x')
ax[1, 0].imshow(I[:,int(particle['y']*100),:], origin='lower')
ax[1, 0].set_title('I y')
ax[2, 0].imshow(I[:,:,int(particle['z']*100)], origin='lower')
ax[2, 0].set_title('I z')
ax[0, 1].imshow(I_fix[int(particle['x']*100),:,:], origin='lower')
ax[0, 1].set_title('I_fix x')
ax[1, 1].imshow(I_fix[:,int(particle['y']*100),:], origin='lower')
ax[1, 1].set_title('I_fix y')
ax[2, 1].imshow(I_fix[:,:,int(particle['z']*100)], origin='lower')
ax[2, 1].set_title('I_fix z')

Comparision image

3D Image of original answer

3D Image of the suggested fix

Hi Simon, thanks for the addition. I realized this, my approach is a bit different. I just iterate the currently accepted solution thrice, such that the periodic images reproduce in the corners and edges of the target box. — Oier Arcelus, Jul 11 '22 at 10:27

imochoa · Answer 2 · 2020-05-25T13:40:29.177

One simple way to do that would be to replicate your circle in the neighboring periodic grids and check the distances from the meshpoints in your current grid to the centers in the neighboring grids:

Your Code:

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 added some example circle parameters:

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

Since your grid has a length of 1x1x1, I guess the spacing between points is 0.01 so:

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

You can double check it by visualizing a slice:

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

Or the full 3D grid (pretty slow)

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

Thank you so much, I came up with my own solution where I check whether the sphere will come out of the grid or not, and if so I will compute the values of the opposite site by recalculating the variable I, with a second auxiliary sphere that has been moved in a given direction. This is much more compressed and elegant though. — Oier Arcelus, May 25 '20 at 14:21
You could also do that here to get a speed boost. You would just have to filter the `centers` list so that it only contains the centers of spheres that actually intersect the main grid Glad I could help :) If you're happy with the answer consider marking it as "accepted" — imochoa, May 25 '20 at 14:56

Spheres with periodic boundaries on a 3D mesh

2 Answers2