def METHOD_NAME(self, uh, bc, index=np.s_[:]):
    gdof = self.number_of_global_dofs()
    phi = self.basis(bc, index=index) # (NQ, NC, ldof)
    cell2dof = self.dof.cell_to_dof(index=index)
    dim = len(uh.shape) - 1
    s0 = 'abcdefg'
    if self.doforder == 'nodes':
        # phi.shape == (NQ, NC, ldof)
        # uh.shape == (..., gdof)
        # uh[..., cell2dof].shape == (..., NC, ldof)
        # val.shape == (NQ, ..., NC)
        s1 = '...ci, {}ci->...{}c'.format(s0[:dim], s0[:dim])
        val = np.einsum(s1, phi, uh[..., cell2dof])
    elif self.doforder == 'vdims':
        # phi.shape == (NQ, NC, ldof)
        # uh.shape == (gdof, ...)
        # uh[cell2dof, ...].shape == (NC, ldof, ...)
        # val.shape == (NQ, NC, ...)
        s1 = '...ci, ci{}->...c{}'.format(s0[:dim], s0[:dim])
        val = np.einsum(s1, phi, uh[cell2dof, ...])
        return val