Improving interpolation with MultiBlock and

[alexkaiser]

Description

I am combining blocks of RectilinearGrid data into a single data structure for further processing. The fields to plot are stored as cell data. If I convert to point data, then at boundaries between the blocks in the combined mesh, there is a visible seam when plotting in paraview. When plotting cell data directly, the cells are adjacent in space and no seams are visible. Code and an image of the seam in a plot follow.

Is there a way to add cells across the blocks so that paraview will apply its standard interpolation there or otherwise removed the seam? Perhaps this is because more cells are needed that go across boundaries of the blocks? Or some other thing with ghost zones?

dir_name = 'eulerian_vars1365'
files = natsorted(glob.glob(dir_name + "/*.vtr"))
blocks = pyvista.MultiBlock([pyvista.RectilinearGrid(f) for f in files])
mesh = blocks.combine()
mesh_points = mesh.cell_data_to_point_data()

Example Data

[akaszynski]

It's possible you might be able to merge the adjacent blocks with clean after extracting the surface.

[alexkaiser]

Calling extract_surface followed by clean is giving a mesh only for the boundaries of the blocks, and discarding the volumetric data. As in, I am only seeing the data at the location of the seams. The combined mesh, which has type UnstructuredGrid, does not appear to have a clean function.

Is there some way to re-assign these cleaned boundaries to the original mesh, or clean without extracting the surface?

Thank you,

[akaszynski]

What about combine(merge_points=True)?

[alexkaiser]

This reduces the number of points by ~85k, which is a plausible number to be on the internal boundaries, but the seam remains in the plots.

[alexkaiser]

There may be a problem with lack of data perfectly coinciding on the boundaries of the vtr files, so I end up with some redundant points. The total number of points should be (NX+1)(NY+1)(NZ+1) = 3164625 on the combined mesh, but instead there are 3266378.

If a point on the boundary of the block is then only connected to cells on its side of the boundary, its interpolated value is determined by only half the cells it should be. Then the redundant point only gets its neighbors on the other side, and the values do not agree, and this creates the seam in plotting. Does this make sense?

[alexkaiser]

I tried calling mesh.cell_centers, which gives point data of the fields at the cell centers.

Is there a way to get cells centered at these newly constructed points? (The dual of the original cells.) The data looks right but only plots in paraview as literal points, not a surface, because of the absence of appropriate cell data.

[alexkaiser]

I have found a solution. Using the cell centers, sort the points by x,y,z, then manually compute new cells connecting the centers. This is possible since it is a regular Cartesian mesh and the original numbers of cells, NX,NY,NZ are known. This works fine but felt relatively complicated and manual for something simple.

Is there a simpler way with filters to accomplish this? For example is there a tolerance parameter for combine(merge_points=True) that will combine points that differ only slightly? Reviewing the source, such a tolerance was not obvious, but perhaps easy to miss.

[akaszynski]

Is there a simpler way with filters to accomplish this? For example is there a tolerance parameter for combine(merge_points=True) that will combine points that differ only slightly? Reviewing the source, such a tolerance was not obvious, but perhaps easy to miss.

I've added this #1595. Please checkout the branch and let me know if that solved your issue.

[alexkaiser]

This works great. I set the tolerance to dx/10 and get exactly the right number of points. Thank you for introducing this prompt and effective fix!