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()
ug1.SetOrigin(minCorner1)
ug1.SetDimensions(dims1)
ug1.SetSpacing(dx1)
# 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
break
if blank:
ug1.BlankCell(cellId)
# Finer uniform grid
minCorner2 = globalOrigin
dims2 = [11]*3
dx2 = [0.5]*3
ug2 = vtk.vtkUniformGrid()
ug2.SetOrigin(minCorner2)
ug2.SetDimensions(dims2)
ug2.SetSpacing(dx2)
# 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
break
if blank:
ug2.BlankCell(cellId)
# Put the two grids in an AMR dataset
amr = vtk.vtkOverlappingAMR()
amr.Initialize(2, [1, 1])
amr.SetOrigin(globalOrigin)
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()
writer.SetFileName('example_amr.vthb')
writer.SetInputDataObject(amr)
writer.Write()
# 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()
writer.SetFileName('example_mg.vtm')
writer.SetInputDataObject(mg)
writer.Write()