Code from Poisson subdomain problem tutorial not working

10

The problem from the tutorial was that the way of separating out the CellBlock file into line and triangular elements (looks like it’s creating a sub-dictionary?) did not work. I replaced

mesh_of_triangles = meshio.Mesh(points=points[:, :2],
                                cells={'triangle': cells['triangle']},
                                cell_data={'triangle': 
                                    {'subdomain': 
                                        cell_data['triangle']
                                        ['gmsh:physical']}},
                                field_data=field_data) 

meshio.write("poisson_subdomain_triangle.xdmf", mesh_of_triangles )

mesh_of_lines = meshio.Mesh(points=points[:, :2],
                                cells={'line': cells['line']},
                                cell_data={'line': {'subdomain':
                                    cell_data['triangle']
                                    ['gmsh:physical']}},
                                field_data=field_data)

with (from the post here)

line_cells = []
for cell in mesh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif  cell.type == "line":
        if len(line_cells) == 0:
            line_cells = cell.data
        else:
            line_cells = np.vstack([line_cells, cell.data])

line_data = []
for key in mesh.cell_data_dict["gmsh:physical"].keys():
    if key == "line":
        if len(line_data) == 0:
            line_data = mesh.cell_data_dict["gmsh:physical"][key]
        else:
            line_data = np.vstack([line_data, mesh.cell_data_dict["gmsh:physical"][key]])
    elif key == "triangle":
        triangle_data = mesh.cell_data_dict["gmsh:physical"][key]

triangle_mesh = meshio.Mesh(points=mesh.points, cells={"triangle": triangle_cells})
line_mesh =meshio.Mesh(points=mesh.points,
                           cells=[("line", line_cells)],
                           cell_data={"name_to_read":[line_data]})

which seems to overcome the problem I originally encountered. However, now I get a different error, which I think warrants its own thread.