Tag PolyData faces about pyvista HOT 4 OPEN

from pyvista import Plane
m0 = Plane(i_resolution=1, j_resolution=2)
m1 = Plane(i_resolution=1, j_resolution=3)
m0.cell_data['tag'] = 0
m1.cell_data['tag'] = 1
m = m0.merge(m1)
m.cell_data['tag']

m0.cell_data['tag'] = ['front' for _ in range(m0.n_cells)]
m1.cell_data['tag'] = ['back' for _ in range(m1.n_cells)]
print(m0.merge(m1).cell_data['tag'])

def mesh_between_profiles(start_profile_xyz_array, end_profile_xyz_array, no_nodes, tag):
    # Initial part of the function to find mesh_xyz_array and faces
    # ...
    # Final part of the function
    mesh_polydata = pv.PolyData()
    mesh_polydata.points = mesh_xyz_array
    mesh_polydata.point_data['tag'] = [tag]*mesh_polydata.n_points
    mesh_polydata.faces = faces
    mesh_polydata.cell_data['tag'] = [tag]*mesh_polydata.n_cells
    return mesh_polydata

from pyvista import Plane
from numpy import zeros, array

m0 = Plane(i_resolution=2, j_resolution=2, center=(0, 0, 0))
m1 = Plane(i_resolution=2, j_resolution=2, center=(1, 0, 0))
m2 = Plane(i_resolution=2, j_resolution=2, center=(2, 0, 0))

# Map tag names to numbers for convenience and speed
tag_ids = dict(front=1, side=2, back=3)
m0.cell_data['tag'] = tag_ids['front']
m1.cell_data['tag'] = tag_ids['side']
m2.cell_data['tag'] = tag_ids['back']
merged = m0.merge([m1, m2]).clean()

# Convert from categorical tags to indicator variables
# Make a (n_cells, n_tags) array which is 1 if that face has that tag, otherwise 0
tag_data = zeros((merged.n_cells, len(tag_ids)))
for i, tag_id in enumerate(tag_ids.values()):
    tag_data[:, i] = merged.cell_data['tag'] == tag_id

merged.cell_data['tag_data'] = tag_data

# Use the 'cell_data_to_point_data' filter
# Notably this filter seems to ignore string-valued arrays
ctp = merged.cell_data_to_point_data()
points_tag_data = ctp.point_data['tag_data']  # (n_points, n_tags)

# Convert back to names
# Because of the averaging, the values in `points_tag_data` will be 0.0 if that point
# never belonged to a face with that tag, 1.0 if all of the faces had that tag,
# and some number inbetween for mixed face tags
names = array(list(tag_ids.keys()))
merged.point_data['tag_names'] = [
    ', '.join(names[avg_indicators > 0])
    for avg_indicators in points_tag_data
]

# Initial part of the script generating the PolyData objects stored as in the list meshes
# ...
# Merge all the mesh segments together
merged_box_beam_mesh = meshes[0].merge(meshes[1:])
# Clean the obtained mesh by merging points closer than a specified tolerance
cleaned_box_beam_mesh = merged_box_beam_mesh.clean(tolerance=element_length / 100)

# Tag the points based on the faces they belong to
# Step 1: gather the indices of points that form each face
cells = cleaned_box_beam_mesh.faces.reshape(-1, 5)[:, 1:]  # assume quad cells
point_indices = cells.flatten()  # flatten the cells array to get a list of point indices, repeated per cell occurrence
cell_tags_repeated = np.repeat(cleaned_box_beam_mesh.cell_data['tag'], 4)  # array of the same shape as point_indices, where each cell tag is repeated for each of its points
# Step 2: Map cell tags to point tags using an indirect sorting approach
sort_order = np.argsort(point_indices)  # get the sort order to rearrange the point indices
sorted_point_indices = point_indices[sort_order]  # sort the point indices
sorted_tags = cell_tags_repeated[sort_order]  # sort the cell tags in the same order
_, boundaries_indices = np.unique(sorted_point_indices, return_index=True)  # find the boundaries between different points in the sorted point indices
# Step 3: Split the sorted tags array at these boundaries to get lists of tags for each point
tags_split = np.split(sorted_tags, boundaries_indices[1:])  # split the sorted tags array at the boundaries to get lists of tags for each point
point_tags_list = np.array([', '.join(tags) for tags in tags_split])  # convert each list of tags into a comma-separated string
cleaned_box_beam_mesh.point_data['tags'] = point_tags_list  # apply the tags to the point_data of the mesh

nodes_xyz_array = cleaned_box_beam_mesh.points  # get xyz coordinates of all nodes
top_skin_points_indices = np.flatnonzero(np.char.find(cleaned_box_beam_mesh.point_data['tags'], "top skin")!=-1)  # find indices of top skin points by looking for "top skin" tag