Labels on generated mesh

[oarcelus]

I am using the 'generate_from_array' function to mesh a numpy array containing labels. Everything works correctly, but I am unable to extract the labels for each of the generated tetrahedra. The medit files that are generated from the 'c3t3.output_to_medit(medit_file);' call, indeed keeps the tetrahedra labels. But cannot manage to find them from the mesh objects generated with pygalmesh.

[uvilla]

I am using pygalmesh to generate meshes for FEniCS 2019.1. Here's my code snippet, where data is the 3D numpy array and
h is the list of voxels sizes in x, y, and z direction.

import numpy as np
import pygalmesh
import dolfin as dl

def generate_mesh(data, h):
    

    dd = data.astype(dtype=np.uint16)

    cell_sizes_map = {}
    cell_sizes_map['default'] = 0.5
        
    

    mesh = pygalmesh.generate_from_array(dd, h, max_cell_circumradius=cell_sizes_map,
                                         max_facet_distance=.5*h[0], verbose=True)

    
    dd_unique = np.unique(dd[:]).astype(np.uint32)

    cells = mesh.get_cells_type('tetra')
    old_labels = mesh.get_cell_data("medit:ref", 'tetra')
    labels = dd_unique[old_labels]
   
    print(labels.shape, cells.shape)

    mesh.cells = [meshio.CellBlock('tetra', cells)] 
    mesh.cell_data["c_labels"] = [labels]

    fname='tmp{0}'.format(np.random.randint(10000,100000-1))
    meshio.xdmf.write(fname+".xdmf", mesh)
    
    dl_mesh = dl.Mesh()
    with dl.XDMFFile(fname+".xdmf") as fid:
        fid.read(dl_mesh)
        geo_dim = dl_mesh.geometry().dim()
        c_labels = dl.MeshFunction('size_t', dl_mesh, geo_dim)
        c_labels.rename("c_labels", "c_labels")
        fid.read(c_labels, "c_labels")
        
    os.remove(fname+".xdmf")
    os.remove(fname+".h5")
    

    return dl_mesh, c_labels

Note that cgalmesh changes the labels of your subdomain and uses a consecutive list of integers.
For examples, if you subdomain labels in data are {10, 20, 30}, the labels returned by cgalmesh will be {0, 1, 2}.

The lines of code:

dd_unique = np.unique(dd[:]).astype(np.uint32)
old_labels = mesh.get_cell_data("medit:ref", 'tetra')
labels = dd_unique[old_labels]

Can be used to restore the original labels.

Also, in my code I am removing from the mesh information regarding all faces, because of the need of my finite element library.

    cells = mesh.get_cells_type('tetra')
    mesh.cells = [meshio.CellBlock('tetra', cells)] 
    mesh.cell_data["c_labels"] = [labels]

@nschloe : Do you think that the relabeling of element labels should be included in generate_from_array or should stay in the user code?

[oarcelus]

Thanks! Turns out I arrived to a very similar solution myself.