I am using h5py to dump meshes with vertices, cells, facets, domain markers, facet markers etc into on single hdf5 file and manually write xdmf compliant xml files for visualization in ParaView without the overhead of redundantly storing parts of the mesh in multiple files.
But, reading back meshes, which are just the Mesh arrays dumped to hdf5 datasets, is tediously slow using MeshEditor.
Is there a faster way to generate Meshes from vertex and cell arrays by passing the whole arrays?
Until now I was following
https://fenicsproject.org/qa/3886/initializing-with-user-defined-coordinates-connectivity
for the mesh creation.
A reduced example for some timings:
from dolfin import *
import h5py, numpy, time
big_number = int(1e3)
mesh = UnitSquareMesh(big_number,big_number)
basedim = mesh.geometry().dim()
num_vertices = mesh.num_vertices()
num_cells = mesh.num_cells()
tac = time.process_time()
hdf5file = HDF5File(mpi_comm_self(),'test1.xdmf','w')
hdf5file.write(mesh,'mesh')
hdf5file.close()
del hdf5file
tec = time.process_time()
print('HDF5File writing {:.2e}s'.format(tec-tac))
hdf5file = HDF5File(mpi_comm_self(),'test1.xdmf','r')
mesh2 = Mesh()
hdf5file.read(mesh2,'mesh',False)
del hdf5file
tic = time.process_time()
print('\nHDF5File reading {:.2e}s'.format(tic-tec))
print('\nHDF5File total {:.2e}s'.format(tic-tac))
tac = time.process_time()
h5file=h5py.File('test.h5','w')
h5file.create_dataset('/vertices',(num_vertices,basedim), data=mesh.coordinates())
h5file.create_dataset('/cells', (num_cells,basedim+1), data=numpy.array(mesh.cells(), dtype=numpy.uintp))
h5file.close()
del h5file
tec = time.process_time()
print('\nh5py writing {:.2e}s'.format(tec-tac))
h5file = h5py.File('test.h5','r')
h5vertices = numpy.array(h5file['/vertices'])
shape_vertices = h5vertices.shape
h5cells = numpy.array(h5file['/cells'])
shape_cells = h5cells.shape
h5file.close()
del h5file
tic = time.process_time()
print('\nh5py raw reading {:.2e}s'.format(tic-tec))
mesh1=Mesh()
edit=MeshEditor()
edit.open(mesh1, basedim, basedim)
edit.init_vertices(shape_vertices[0])
edit.init_cells(shape_cells[0])
for ii in range(shape_vertices[0]):
edit.add_vertex(ii, h5vertices[ii])
for ii in range(shape_cells[0]):
edit.add_cell(ii, h5cells[ii])
edit.close()
toc = time.process_time()
print('h5py mesh editor {:.2e}s'.format(toc-tic))
print('h5py total reading {:.2e}s'.format(toc-tec))
print('\nh5py total {:.2e}s'.format(toc-tac))
Essentially I don't know a good way to create a mesh from given arrays.
Is there something like a reasonable fast
mesh = Mesh(h5vertices, h5cells)
Also, curiosly HDF5File.write(mesh,'mesh') apparently does a lot more than just writing vertices and cells as datasets.
