[FIXED] Spheres with periodic boundaries on a 3D mesh

Issue

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

Solution

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

Answered By - Simon