docs: def Wm(self):
          """Return the smoothing regularization matrix Wm of the grid
  
          """
          centroids = self.get_element_centroids()
  
          Wm = scipy.sparse.csr_matrix(
              (self.nr_of_elements, self.nr_of_elements))
          # Wm = np.zeros((self.nr_of_elements, self.nr_of_elements))
          for i, nb in enumerate(self.element_neighbors):
              for j, edges in zip(nb, self.element_neighbors_edges[i]):
                  # side length
                  edge_coords = self.nodes['presort'][edges][:, 1:]
                  edge_length = np.linalg.norm(
                      edge_coords[1, :] - edge_coords[0, :]
                  )
                  distance = np.linalg.norm(centroids[i] - centroids[j])
  
                  # main diagonal
                  Wm[i, i] += edge_length / distance
                  # side diagonals
                  Wm[i, j] -= edge_length / distance
          return Wm