Open oarcelus opened 2 years ago
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?
Thanks! Turns out I arrived to a very similar solution myself.
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.