Per-vertex normals (if the mesh has any) are scaled as well, using the inverse transpose of the scaling matrix: a non-uniform scale would otherwise leave normals pointing in a direction that no longer matches the surface, which shows up as wrong shading.
Examples
mesh <- generate_cuboid(c(0, 0, 0), c(1, 1, 1))
big <- scale_mesh(mesh, 3)
flat <- scale_mesh(mesh, c(2, 0.5, 1))