The problem is that in the newest meshio, the input msh file split the cells in CellBlocks. See for example: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio where I combine all the blocks of the triangular elements. You need to do the same for tetrahedron elements (and triangular elements).
To make sure that you read in your correct mesh (before doing calculations) compute volume and surface areas (using assemble(Constant(1)*dx(domain=mesh, subdomain_data=... )) and similarly for dS and ds.