def METHOD_NAME(geometry, npts_irregular):
    filename = os.path.join(mesh_dir, geometry)
    domain = Domain.from_file(filename=filename)
    domainh = discretize(domain, filename=filename, comm=comm)
    mapping = list(domainh.mappings.values())[0]
    space = mapping.space
    # Irregular tensor grid
    grid_i = [np.linspace(space.breaks[i][0], space.breaks[i][-1], npts_irregular, True) for i in range(space.ldim)]
    # This calls eval_fields
    mesh_part = mapping.build_mesh(grid_i)
    cell_indexes = [cell_index(space.breaks[i], grid_i[i]) for i in range(space.ldim)]
    starts, ends = space.local_domain
    actual_starts = tuple(np.searchsorted(cell_indexes[i], starts[i], side='left')
                          for i in range(space.ldim))
    actual_ends = tuple(np.searchsorted(cell_indexes[i], ends[i], side='right')
                         for i in range(space.ldim))
    index = tuple(slice(s, e, 1) for s,e in zip(actual_starts, actual_ends))
    shape = tuple(len(grid_i[i]) for i in range(space.ldim))
    # Saving in an hdf5 file to compare on root
    fh5 = h5.File(f'result_parallel_i.h5', mode='w', driver='mpio', comm=comm)
    for i, dset_name in enumerate(['x_mesh', 'y_mesh', 'z_mesh'][:space.ldim]):
        fh5.create_dataset(dset_name, shape=shape, dtype=float)
        fh5[dset_name][index] = mesh_part[i]
    fh5.close()
    # Check
    if rank == 0:
        domainh = discretize(domain, filename=filename, comm=None)
        mapping = list(domainh.mappings.values())[0]
        space = mapping.space
        mesh = mapping.build_mesh(grid_i)
        fh5 = h5.File(f'result_parallel_i.h5', mode='r')
        for i, dset_name in enumerate(['x_mesh', 'y_mesh', 'z_mesh'][:space.ldim]):
            dset = fh5[dset_name][...]
            assert np.allclose(dset, mesh[i], atol=ATOL, rtol=RTOL)
        fh5.close()
        os.remove('result_parallel_i.h5')