code: def _energy_uniform_per_node(X, cells):
      
      dim = 2
      mesh = MeshTri(X, cells)
  
      star_integrals = numpy.zeros(mesh.node_coords.shape[0])
      # Python loop over the cells... slow!
      for cell, cell_volume in zip(mesh.cells["nodes"], mesh.cell_volumes):
          for idx in cell:
              xi = mesh.node_coords[idx]
              tri = mesh.node_coords[cell]
              val = quadpy.triangle.integrate(
                  lambda x: numpy.einsum("ij,ij->i", x.T - xi, x.T - xi),
                  tri,
                  # Take any scheme with order 2
                  quadpy.triangle.Dunavant(2),
              )
              star_integrals[idx] += val
  
      return star_integrals / (dim + 1)
docs: The CPT mesh energy is defined as
  
          sum_i E_i,
          E_i = 1/(d+1) * sum int_{omega_i} ||x - x_i||^2 rho(x) dx,
  
      see Chen-Holst. This method gives the E_i and  assumes uniform density, rho(x) = 1.