use 'norbs' returned by 'hamiltonian_submatrix' in 'kpm.conductivity'

This does not require assuming that site families have 'norbs' set.

Closes #313.

kwant/kpm.py

@@ -1051,10 +1051,11 @@ def _velocity(hamiltonian, params, op_type, positions):
     elif isinstance(op_type, str):
         direction = directions[op_type]
         if isinstance(hamiltonian, system.System):
-            operator = hamiltonian.hamiltonian_submatrix(params=params,
-                                                         sparse=True)
-            positions = np.array([site.pos for site in hamiltonian.sites
-                                  for iorb in range(site.family.norbs)])
+            operator, norbs, norbs = hamiltonian.hamiltonian_submatrix(
+                params=params, sparse=True, return_norb=True
+            )
+            positions = np.vstack([[hamiltonian.pos(i)] * norb
+                                   for i, norb in enumerate(norbs)])
         elif positions is not None:
             operator = coo_matrix(hamiltonian, copy=True)
         displacements = positions[operator.col] - positions[operator.row]