I/O from XDMF/HDF5 files in dolfin-x

I want to share what is now working much faster for me than my old approach via the coordinates.
In preprocessing (serial), read in your mesh with the standard XDMF routine, set values to a function, and write a text file using
igi = msh.geometry.input_global_indices

with the following structure

igi1 fx1 fy1 fz1
igi2 fx2 fy2 fz2
...
ign fxn fyn fzn

so the input node index followed by the x,y,z values of your vector field (or the one value of a scalar field).

Then, in parallel, use

def readfunction(f, V, datafile, msh):

    # block size of vector
    bs = f.vector.getBlockSize()

    # load data and input node indices
    data = np.loadtxt(datafile,usecols=(np.arange(1,bs+1)),ndmin=2)
    ind_file = np.loadtxt(datafile,usecols=(0),dtype=int)

    # index map and input indices
    im = np.asarray(V.dofmap.index_map.local_to_global(np.arange(V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts, dtype=np.int32)), dtype=PETSc.IntType)
    igi = msh.geometry.input_global_indices

    # since in parallel, the ordering of the dof ids might change, so we have to find the
    # mapping between original and new id via the coordinates
    ci = 0
    for i in im:

        ind = np.where(ind_file == igi[ci])[0]

        # only write if we've found the index
        if len(ind):
            for j in range(bs):
                f.vector[bs*i+j] = data[ind[0],j]

        ci+=1

    f.vector.assemble()
    f.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
1 Like