Export custom normals I am writing an importer and exporter for a custom binary format.
In this format I have triangulated mesh. Each vertex has position, normal, and UV.

The importer code sets them like this and it is displayed properly:

mesh = MyMesh(file_name)
bl_mesh = bpy.data.meshes.new(object_name)
bl_mesh.from_pydata(mesh.positions, [], mesh.indices)
bl_mesh.create_normals_split()
for l in bl_mesh.loops:
    l.normal[:] = mesh.normals[l.vertex_index]

bl_mesh.normals_split_custom_set_from_vertices(mesh.normals)
bl_mesh.use_auto_smooth = True


But the exporter code does not give me the same results.

mesh = MyMesh()
bl_mesh.calc_normals_split()
for l in bl_mesh.loops:
    mesh.indices.append(l.vertex_index)
    mesh.normals.insert(l.vertex_index, (*l.normal[:]))




I have tried accessing the obj.data.vertices[].normal but the normals there are wrong as well.

How do I get custom normals from the mesh?

UPD.
The full source code of the extension is in this repo. The model used is nvx/model.nvx from the repo. All interactions with blender objects happen in the io_nebula/bl_nvx.py
After hours of debugging I found out that list.insert function does not do what I expected (I was not alone in this assumption https://stackoverflow.com/q/25840177/3344612). And the order of vertices was messed because I inserted at a wrong index.

Answering my question. To get custom normals first call calc_normals_split() on the mesh and then access through loops[].normal with corresponding vertex indices in loops[].vertex_index.

To fix the mentioned code first resize the list and then modify with [] operator.

mesh.normals = [None] * len(bl_mesh.vertices)
for l in bl_mesh.loops:
    mesh.indices.append(l.vertex_index)
    mesh.normals[l.vertex_index] = l.normal[:]