Blender 3.x merge objects (non-destructive)

import bpy, bmesh

def merge_objects(objects, name=None, obj_matrix=None):
    '''
    Non-destructive merge - Objects will be merged IN PLACE wherever they are in the world.
    setting obj_matrix simply defines the new object's origin matrix,
    otherwise the first object in the list will be used as the new object's origin (visual
    geometry will not be affected).
    '''

    if name is None:
        name = "JoinedMesh"

    # create a temp mesh to copy everything into
    if name in bpy.data.meshes:
        bpy.data.meshes.remove(bpy.data.meshes[name])

    final_mesh = bpy.data.meshes.new(name)

    def add_data_to_mesh(o, me, extra_transform=None):
        '''
        Creates a bmesh for a given mesh, then adds it to the existing output mesh data.
        extra_transform is used for instances and nested instances.
        '''

        # 'bake' visual normals to custom split normals, regardless of what type of smoothing is being used. This ensures our final
        # mesh is visually consistent with the source objects- if something looked faceted and sharp, it should look faceted and sharp
        # in the merged object even if it was merged into something with weighted normals, or something that was had autosmoothing set
        # to 180 degrees.
        me.calc_normals_split()
        me.normals_split_custom_set([l.normal for l in me.loops])

        # add any new materials to the temp mesh and create a re-indexing map
        mat_index_lookup = {}
        if not me.materials:
            # adding a null material here to avoid a situation where a mesh that does not use a material gets automatically assigned the
            # first mat that was found (which is what happens if you use the built-in join operator)
            me.materials.append(None)
        for i, m in enumerate(me.materials):
            if m not in final_mesh.materials.values():
                mat_index_lookup[i] = len(final_mesh.materials)
                final_mesh.materials.append(m)
            else:
                new_idx = final_mesh.materials.values().index(m)
                if i != new_idx:
                    mat_index_lookup[i] = final_mesh.materials.values().index(m)

        bm = bmesh.new()
        bm.from_mesh(me)

        mat_face_reindex = [f for f in bm.faces if f.material_index in mat_index_lookup]
        for f in mat_face_reindex:
            f.material_index = mat_index_lookup[f.material_index]

        # transform the bmesh to match the object transform
        if extra_transform is not None:
            bmesh.ops.transform(bm, matrix=extra_transform @ o.evaluated_get(bpy.context.view_layer.depsgraph).matrix_world, verts=bm.verts)
        else:
            bmesh.ops.transform(bm, matrix=o.evaluated_get(bpy.context.view_layer.depsgraph).matrix_world, verts=bm.verts)

        # copy in the existing mesh data before writing everything back out
        bm.from_mesh(final_mesh)
        bm.to_mesh(final_mesh)
        return bm

    def add_mesh_object(co):
        return add_data_to_mesh(co.evaluated_get(bpy.context.view_layer.depsgraph), co.evaluated_get(bpy.context.view_layer.depsgraph).data, extra_transform=this_matrix_offset)

    def handle_object_recursive(o, extra_transform=None):
        bm = None
        if o.type == 'EMPTY' and o.instance_collection is not None:
            # loop through each object in the collection and add the data to the mesh
            for co in o.instance_collection.all_objects:
                this_matrix_offset = None
                if extra_transform is not None:
                    # if extra_transform isn't None, o is already a nested instance. need to offset each object in the source
                    # collection by the incoming extra_transform AND the local offset of the instance empty itself. since
                    # this function is recursive, any additional nested instances will inherit this transformation and will
                    # continuously update the matrix offset to match the visual transform of that particular geometry.
                    this_matrix_offset = extra_transform @ o.matrix_world
                else:
                    # if extra_transform is None, it means we're working with a first-level object inside an instance.
                    # We only need to offset the overall transformation by the instance empty, the contents of the instance will
                    # be handled recursively, which is all kicked off by setting the initial offset here.
                    this_matrix_offset = o.matrix_world
                if co.type == 'MESH':
                    bm = add_data_to_mesh(co.evaluated_get(bpy.context.view_layer.depsgraph), co.evaluated_get(bpy.context.view_layer.depsgraph).data, extra_transform=this_matrix_offset)
                elif co.type == 'EMPTY' and co.instance_collection is not None:
                    bm = handle_object_recursive(co, extra_transform=this_matrix_offset)
                elif co.type == 'CURVE':
                    bm = add_data_to_mesh(co, co.evaluated_get(bpy.context.view_layer.depsgraph).to_mesh(), extra_transform=this_matrix_offset)
                    o.evaluated_get(bpy.context.view_layer.depsgraph).to_mesh_clear()
        elif o.type == 'MESH':
            bm = add_data_to_mesh(o, o.evaluated_get(bpy.context.view_layer.depsgraph).data)
        elif o.type == 'CURVE':
            bm = add_data_to_mesh(o, o.evaluated_get(bpy.context.view_layer.depsgraph).to_mesh())
            o.evaluated_get(bpy.context.view_layer.depsgraph).to_mesh_clear()
        else:
            # TODO: handle other mesh-like types, such as metaballs.
            print(f"Skipping {o.name}, unsupported type ({o.type})")
        return bm



    # kick off the recursive bmesh builder
    bm = None
    for i, o in enumerate(objects):
        bm = handle_object_recursive(o)

    if bm is None:
        bpy.data.meshes.remove(final_mesh)
        return None

    # need autosmoothing turned on to see custom split normals
    final_mesh.use_auto_smooth = True

    # create the object itself and handle the matrix offset to get the origin into the requested location.
    ob = bpy.data.objects.new(name, final_mesh)
    bpy.context.scene.collection.objects.link(ob)

    if obj_matrix is not None:
        ob.matrix_world = obj_matrix
        bmesh.ops.transform(bm, matrix=obj_matrix.inverted(), verts=bm.verts)
    else:
        ob.matrix_world = objects[0].matrix_world
        bmesh.ops.transform(bm, matrix=objects[0].matrix_world.inverted(), verts=bm.verts)

    # done
    bm.to_mesh(final_mesh)
    bm.free()
    bpy.context.view_layer.update()
    return ob

context = bpy.context
selected_objects = context.selected_objects
result = merge_objects(selected_objects, name="MergedResult", obj_matrix=context.active_object.matrix_world if context.active_object else None)

if result:
    bpy.ops.object.select_all(action="DESELECT")
    result.select_set(True)
    context.view_layer.objects.active = result