`bmesh.ops.subdivide_edges` corrupts vertex references depending on subdivision order

I'm working on a mesh processing tool using BMesh in Blender, where I subdivide a set of low-resolution edges to match the resolution of a high-resolution mesh.

The process uses bmesh.ops.subdivide_edges with a variable number of cuts per edge, calculated proportionally based on edge length. To apply this, I group edges by the number of cuts, then run:

for num_cuts, edges in edges_by_cuts.items():
    geometry = bmesh.ops.subdivide_edges(self.bm, edges=edges, cuts=num_cuts)
    ...

After each subdivision, I retrieve the new BMVert instances and track them via a static mapping system.

The problem

The vertices returned from bmesh.ops.subdivide_edges can become unreliable (invalid or mismatched), depending on the order in which cuts are applied. The default order for my script was the following: [7, 9, 2, 3, 4, 11, 14, 13, 1, 5, 17, 8, 12, 6, 10]. Running my script would correctly subdivide edges based on the number of cuts defined in the list above, but when trying to select the new geometry, i.e. the new vertices created by the subdivision, the vertices of the edges with 7 cuts (the first subdivision) would be wrongly selected. I decided to run some tests to try understanding what happened.

Tests

Test 1

Running only the first subdivision with 7 cuts would allow me to correctly select the newly created vertices. Running only the two first subdivisions would cause the vertices of the edges with 7 cuts to be wrongly selected.

Conclusion: The bug must occur between the first and second subdivisions only.

Test 2

Running _build_static_mapping only one time after all subdivisions allow me to select only the vertices from the very last subdivision.

Conclusion: I should be calling _build_static_mapping for each subdivision.

Test 3

Running all subdivisions except the first one would cause all the newly created vertices to be correctly selected.

Conclusion: The bug must be caused by the specific number of 7 cuts.

Test 4

When running small cuts first by sorting the order, vertices of edges with 1 or 2 cuts are wrongly selected. The rest is correctly selected.

Conclusion: The bug does seem to occur during the first subdivisions.

Test 5

Replacing the number of cuts 7 with 1 or 6 cuts, vertices of edges with originally 7 cuts gets wrongly selected.

Conclusion: The number of cuts 7 does not seem to be causing the bug.

Test 6

Running the subdivision with 7 cuts at the very end allows me to successfully select all newly created vertices.

Conclusion: I need your help to figure out what's going on.

TL;DR

When using bmesh.ops.subdivide_edges on multiple groups of edges with different cuts values, the order in which I apply the subdivisions affects whether the returned vertices are correct. Some combinations (e.g., processing cuts=7 early) cause invalid or mismatched vertex selections. Splitting the process and running cuts=7 last consistently solves the issue. I'm looking to understand why this happens and whether there's a more robust solution.

After running the tests, I'm left with these questions:

• Why does the selection of subdivision vertices become inconsistent depending on the cut order?
• Is this a known limitation or bug in BMesh?
• Is there a more reliable way to get the new vertices created by subdivision?

The current fix is to replace the subdivision loop with:

for num_cuts, edges in sorted(edges_by_cuts.items(), key=lambda item: item[0] == 7)

I'd like to understand and find a more generic solution. Any help or clarification would be greatly appreciated.

You can find the detailed code below.

from typing import Set, Tuple
from collections import defaultdict
import bmesh
import bpy
import math

class MeshCleaner:
    def __init__(self, obj):
        self.obj = obj
        self.bm = None
        self.static_to_bmvert = {}  # static index -> bmvert
        self.bmvert_to_static = {}  # bmvert -> static index

        self.bm = bmesh.new()
        self.bm.from_mesh(obj.data)
        self.bm.verts.ensure_lookup_table()

        mesh_verts = obj.data.vertices

        # Build mapping between static indices and bmesh vertices
        for bmvert, mesh_vert in zip(self.bm.verts, mesh_verts):
            static_index = mesh_vert.index
            self.static_to_bmvert[static_index] = bmvert
            self.bmvert_to_static[bmvert] = static_index

    def clean_mesh(self, data):
        low_res_edges: Set[Tuple[int, int]] = data.get('lre')
        high_res_edges: Set[Tuple[int, int]] = data.get('hre')

        original_low_static_verts = set()
        for v1, v2 in low_res_edges:
            original_low_static_verts.add(v1)
            original_low_static_verts.add(v2)

        original_high_static_verts = set()
        for v1, v2 in high_res_edges:
            original_high_static_verts.add(v1)
            original_high_static_verts.add(v2)

        new_static_vertices = self._subdivide_low_res_edges_bmesh(
            low_res_edges=low_res_edges,
            high_res_edges=high_res_edges
        )

        vertex_indices = []
        vertex_indices.extend(new_static_vertices)
        # vertex_indices.extend(list(original_low_static_verts))
        # vertex_indices.extend(list(original_high_static_verts))

        self._finish_bmesh()
        self._select_vertices(vertex_indices=vertex_indices)


    def _subdivide_low_res_edges_bmesh(self, low_res_edges, high_res_edges):
        current_points = len(low_res_edges)
        target_points = len(high_res_edges)
        new_points = target_points - current_points

        if new_points <= 0:
            raise ValueError("No subdivisions needed.")

        # Measure lengths of low-res edges
        edge_lengths = []
        total_length = 0.0
        edge_map = {}

        static_set = {frozenset(edge) for edge in low_res_edges}
        
        for edge in self.bm.edges:
            static_indices = [self.bmvert_to_static.get(v) for v in edge.verts]
            if frozenset(static_indices) in static_set:
                v1 = self.static_to_bmvert[static_indices[0]].co
                v2 = self.static_to_bmvert[static_indices[1]].co
                length = (v1 - v2).length
                edge_lengths.append((edge, length))
                total_length += length
                edge_map[edge] = static_indices
                logger.info(f"{len(edge_map)} , {edge} , {edge_map[edge]}")

        # Proportional subdivision assignment
        ideal_subdivs = [(edge, length / total_length * new_points) for edge, length in edge_lengths]

        floored = [(edge, int(math.floor(val))) for edge, val in ideal_subdivs]
        remainder = new_points - sum(val for _, val in floored)

        fractions = sorted(
            [(edge, val - math.floor(val)) for edge, val in ideal_subdivs],
            key=lambda x: -x[1]
        )

        subdivs = {edge: val for edge, val in floored}
        for i in range(remainder):
            edge = fractions[i][0]
            subdivs[edge] += 1

        # Group edges by number of cuts
        edges_by_cuts = defaultdict(list)
        for edge, cuts in subdivs.items():
            cuts = max(1, cuts)
            edges_by_cuts[cuts].append(edge)

        # Subdivide edges with bmesh.ops.subdivide_edges
        all_new_verts = []
        static_indices = []

        for num_cuts, edges in edges_by_cuts.items():
            geometry = bmesh.ops.subdivide_edges(self.bm, edges=edges, cuts=num_cuts)

            bmverts = [v for v in geometry['geom'] if isinstance(v, bmesh.types.BMVert)]
            valid_bmverts = [v for v in bmverts if v.is_valid]
            all_new_verts.extend(valid_bmverts)
            self._build_static_mapping()
            static_indices.extend([self.bmvert_to_static[v] for v in valid_bmverts if v in self.bmvert_to_static])            

        return static_indices

    def _finish_bmesh(self):
        self.bm.to_mesh(self.obj.data)
        self.bm.free()
        self.bm = None

    def _select_vertices(self, vertex_indices=None):
        if vertex_indices:
            if bpy.context.mode != 'OBJECT':
                bpy.ops.object.mode_set(mode='OBJECT')
            for idx in vertex_indices:
                if 0 <= idx < len(self.obj.data.vertices):
                    self.obj.data.vertices[idx].select = True

    def _prepare_selection(self, mode):
        if bpy.context.mode != 'EDIT':
            bpy.ops.object.mode_set(mode='EDIT')
        bpy.ops.mesh.select_mode(type=mode)
        bpy.ops.mesh.select_all(action='DESELECT')
        bpy.ops.object.mode_set(mode='OBJECT')

    def _build_static_mapping(self):
        self.bmvert_to_static = {v: i for i, v in enumerate(self.bm.verts)}