def METHOD_NAME(self, bc, cellidx=None):
    mesh = self.mesh
    n, t = mesh.edge_frame()
    NE = mesh.number_of_edges()
    idx = np.array([(0, 0), (0, 1), (1, 1)])
    TT = np.zeros((NE, 3, 3), dtype=np.float)
    TT[:, 0, ...] = np.prod(t[:, idx], axis=-1) 
    TT[:, 1, ...] = np.sum(t[:, idx]*n[:, idx[:, [1, 0]]], axis=-1)/2
    TT[:, 2, ...] = np.prod(n[:, idx], axis=-1) 
    NC = mesh.number_of_cells()
    if len(phi0.shape) == 1:
        shape0 = (NC,)+ phi0.shape + (3, 2, 2)
        shape1 = (NC, -1, 2, 2)
    elif len(phi0.shape) == 2:
        shape0 = (phi0.shape[0], NC, phi0.shape[1], 3, 2, 2)
        shape1 = (phi0.shape[0], NC, -1, 2, 2)
    else:
        shape0 = phi0.shape[:-1] + (NC, phi0.shape[-1], 3, 2, 2)
        shape1 = phi0.shape[:-1] + (NC, -1, 2, 2)
    phi = np.zeros(shape0, dtype=np.float)
    dof = self.scalarspace.dof
    multiIndex = dof.multiIndex 
    isEdgeDof = (multiIndex == 0) & (multiIndex != p)
    cell2edge = mesh.ds.cell_to_edge()
    phi0 = self.scalarspace.METHOD_NAME(bc)
    for i in range(3):
        VAL = np.einsum('ijk, kmn->ijmn', TT[cell2edge[:, i]], self.T)
        phi[..., isEdgeDof[:, i], :, :, :] = np.einsum('...k, ijmn->...ikjmn', phi0[..., isEdgeDof[:, i]], VAL) 
    isOtherDof = isEdgeDof[:, 0] | isEdgeDof[:, 1] | isEdgeDof[:, 2]
    phi[..., isOtherDof, :, :, :] = np.einsum('...k, jmn, i->...ikjmn', phi0[..., isOtherDof], self.T, np.ones(NC))
    return phi.reshape(shape1)