flipped axes for MaterialGrid when using Block geometric objects

[oskooi]

It seems there is a bug in the MaterialGrid object in which the coordinate axes are flipped. This can be demonstrated when trying to stretch a circle with radius 0.3 into an ellipse with major/minor axis lengths of 0.6/0.3 as shown in the two figures below. The code used to generate these results is also provided. The correct dimensions of the Block object containing the MaterialGrid is given in the figure titles but the actual geometry shows that the coordinates are in fact flipped.

This could be a possible explanation for the unexpectedly large relative errors observed in #1568 (see comment). Unfortunately, the unit test in test_material_grid.py does not involve this particular use case which is probably why this bug has thus far gone unnoticed.

import meep as mp
import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

res = 100

sxy = 2.0
cell_size = mp.Vector3(sxy,sxy,0)

## radius of circle                                                                                               
rad = 0.3

design_region_res = 2*res
design_shape = mp.Vector3(1.0,1.0,0)
Nx = int(design_region_res*design_shape.x)
Ny = int(design_region_res*design_shape.y)
x = np.linspace(-0.5*design_shape.x,0.5*design_shape.x,Nx)
y = np.linspace(-0.5*design_shape.y,0.5*design_shape.y,Ny)
xv, yv = np.meshgrid(x,y)
design_params = np.sqrt(np.square(xv) + np.square(yv)) < rad
filtered_design_params = gaussian_filter(design_params, sigma=3.0, output=np.double)
beta = 1000
eta = 0.5
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                          mp.air,
                          mp.Medium(index=3.5),
                          weights=filtered_design_params,
                          do_averaging=True,
                          beta=beta,
                          eta=eta)

## try stretching the circle into an ellipse with major/minor axis lengths of 0.3/0.6 in x/y directions
stretch = mp.Vector3(0.3/rad,0.6/rad,0)
geometry = [mp.Block(center=mp.Vector3(),
                     size=stretch,
                     material=matgrid)]

sim = mp.Simulation(resolution=res,
                    cell_size=cell_size,
                    geometry=geometry)

## output entire epsilon grid                                       
x = np.linspace(-0.5*sxy,0.5*sxy,sxy*design_region_res)
y = np.linspace(-0.5*sxy,0.5*sxy,sxy*design_region_res)
z = np.array([0])
eps_grid = sim.get_epsilon_grid(x,y,z)
plt.figure()
plt.pcolor(x,y,eps_grid,cmap='Blues')
plt.gca().axis('equal')
plt.gca().set_aspect('equal','box')
plt.title('Block: {} x {}, MaterialGrid: {} x {}'.format(stretch.x,stretch.y,Nx,Ny))
plt.savefig('matgrid_stretch_bug.png',dpi=160,bbox_inches='tight')
plt.close()

[stevengj]

I think this is just a confusion about plotting?

When you get a 2d array from Meep it gives you X×Y. However, when you plot with pcolor, I think it shows the first dimension (the columns of the matrix) as "vertical", i.e. it expects Y×X. (Moreover, it plots the columns assuming that the first row corresponds to the "top" of the plot, whereas Meep's data in each dimension starts at the smallest coordinate, which would correspond to the "bottom" of the plot).

[oskooi]

Yes, good catch: the axes are indeed flipped by pcolor and there is no bug in MaterialGrid. To verify this, I used a rectangular cell (rather than a square cell in the example above) so that the axes can be clearly identified. With this change, stretching the circle along the y direction (which is the longer side of the rectangular cell) matches the horizontal axis in the plot.

import meep as mp
import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

res = 100

sx = 2.0
sy = 3.0
cell_size = mp.Vector3(sx,sy,0)

## radius of circle 
rad = 0.3

design_region_res = 2*res
design_shape = mp.Vector3(1.0,1.0,0)
Nx = int(design_region_res*design_shape.x)
Ny = int(design_region_res*design_shape.y)
x = np.linspace(-0.5*design_shape.x,0.5*design_shape.x,Nx)
y = np.linspace(-0.5*design_shape.y,0.5*design_shape.y,Ny)
xv, yv = np.meshgrid(x,y)
design_params = np.sqrt(np.square(xv) + np.square(yv)) < rad
filtered_design_params = gaussian_filter(design_params, sigma=3.0, output=np.double)
beta = 1000
eta = 0.5
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                          mp.air,
                          mp.Medium(index=3.5),
                          weights=filtered_design_params,
                          do_averaging=True,
                          beta=beta,
                          eta=eta)

stretch = mp.Vector3(0.3/rad,0.6/rad,0)
geometry = [mp.Block(center=mp.Vector3(),
                     size=stretch,
                     material=matgrid)]

sim = mp.Simulation(resolution=res,
                    cell_size=cell_size,
                    geometry=geometry)

## output epsilon_grid sampled at the same resolution as input grid
x = np.linspace(-0.5*sx,0.5*sx,sx*design_region_res)
y = np.linspace(-0.5*sy,0.5*sy,sy*design_region_res)
z = np.array([0])
eps_grid = sim.get_epsilon_grid(x,y,z)
plt.figure()
plt.pcolor(y,x,eps_grid,cmap='Blues')
plt.gca().axis('equal')
plt.gca().set_aspect('equal','box')
plt.title('Block: {} x {}, MaterialGrid: {} x {}'.format(stretch.x,stretch.y,Nx,Ny))
plt.savefig('matgrid_stretch_bug.png',dpi=160,bbox_inches='tight')
plt.close()