def _energy_uniform_per_node(X, cells):
    """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.
    """
    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)
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.