Blanking in AMR datasets

I’m trying to understand how cell blanking works in AMR datasets. I’m using the script below to generate a dataset with two uniform grid blocks. I write those grids out as a vtkOverlappingAMR dataset and as a multi-block dataset. Then I load the two files into ParaView and make a slice at X=3. The blanking only works in the multi-block dataset, but not in the AMR dataset. I’ve attached an image showing this, created in ParaView 5.9.0.


Did I do something wrong in creating the AMR dataset? Is vtkOverlappingAMR the right type for this data? Is there any advantage to using an AMR type vs. multi-block? For large datasets, will ParaView be faster to create slices, iso-surfaces, etc. on the AMR data?

import numpy as np
import vtk
from vtk.util.numpy_support import numpy_to_vtk

globalOrigin = [0.0]*3

# Coarse uniform grid
minCorner1 = globalOrigin
dims1 = [11]*3
dx1 = [1.0]*3
ug1 = vtk.vtkUniformGrid()

# Set iblanking on coarse grid.  Blank where the coarse grid overlaps the
# fine grid, which we'll construct in the next step.
ib1 = np.ones(ug1.GetNumberOfPoints())
for i in range(6):
    for j in range(6):
        for k in range(6):
            pid = vtk.vtkStructuredData.ComputePointId(ug1.GetDimensions(), [i, j, k])
            ib1[pid] = 0.0
for cellId in range(ug1.GetNumberOfCells()):
    ptIds = vtk.vtkIdList()
    ug1.GetCellPoints(cellId, ptIds)
    blank = True
    for i in range(ptIds.GetNumberOfIds()):
        pid = ptIds.GetId(i)
        if ib1[pid] == 1.0:
            blank = False
    if blank:

# Finer uniform grid
minCorner2 = globalOrigin
dims2 = [11]*3
dx2 = [0.5]*3
ug2 = vtk.vtkUniformGrid()

# Set iblanking on fine grid.  Blank a hole in the middle of the box.
ib2 = np.ones(ug2.GetNumberOfPoints())
for i in range(3,7):
    for j in range(3,7):
        for k in range(3,7):
            pid = vtk.vtkStructuredData.ComputePointId(ug2.GetDimensions(), [i, j, k])
            ib2[pid] = 0.0
for cellId in range(ug2.GetNumberOfCells()):
    ptIds = vtk.vtkIdList()
    ug2.GetCellPoints(cellId, ptIds)
    blank = True
    for i in range(ptIds.GetNumberOfIds()):
        pid = ptIds.GetId(i)
        if ib2[pid] == 1.0:
            blank = False
    if blank:

# Put the two grids in an AMR dataset
amr = vtk.vtkOverlappingAMR()
amr.Initialize(2, [1, 1])

box1 = vtk.vtkAMRBox(minCorner1, dims1, dx1, globalOrigin)
amr.SetAMRBox(0, 0, box1)
amr.SetDataSet(0, 0, ug1)
amr.SetRefinementRatio(0, 2)
amr.SetSpacing(0, dx1)

box2 = vtk.vtkAMRBox(minCorner2, dims2, dx2, globalOrigin)
amr.SetAMRBox(1, 0, box2)
amr.SetDataSet(1, 0, ug2)
amr.SetRefinementRatio(1, 2)
amr.SetSpacing(1, dx2)

# Write out the AMR dataset
writer = vtk.vtkXMLUniformGridAMRWriter()

# Now write them out as a multi-block dataset
mg = vtk.vtkMultiBlockDataSet()
mg.SetBlock(0, ug1)
mg.GetMetaData(0).Set(vtk.vtkCompositeDataSet.NAME(), 'block1')
mg.SetBlock(1, ug2)
mg.GetMetaData(1).Set(vtk.vtkCompositeDataSet.NAME(), 'block2')

writer = vtk.vtkXMLMultiBlockDataWriter()

Blanking seems to work fine with my example dataset in 5.9.1-1863-g8d7fa82. In my original post, I was using 5.9.0-RC2.