We develop a unified and easy to use framework to study robust fully discrete numerical methods for nonlinear degenerate diffusion equations partial derivative(t) u - L-sigma,L- mu[phi(u)] = f in R-N x (0, T), where L-sigma,L- mu is a general symmetric diffusion operator of Levy type and partial derivative is merely continuous and nondecreasing. We then use this theory to prove convergence for many different numerical schemes. In the nonlocal case most of the results are completely new. Our theory covers strongly degenerate Stefan problems, the full range of porous medium equations, and, for the first time for nonlocal problems, also fast diffusion equations. Examples of diffusion operators L-sigma,L- mu are the (fractional) Laplacians Delta and -(-Delta)alpha/2 for alpha is an element of (0, 2), discrete operators, and combinations. The observation that monotone finite difference operators are nonlocal Levy operators allows us to give a unified and compact nonlocal theory for both local and nonlocal linear and nonlinear diffusion equations. The theory includes stability, compactness, and convergence of the methods under minimal assumptions, including assumptions that lead to very irregular solutions. As a byproduct, we prove the new and general existence result announced in [F. del Teso, J. Endal, and E. R. Jakobsen, C. R. Math. Acad. Sci. Paris, 355 (2017), pp. 1154-1160]. We also present some numerical tests, but extensive testing is deferred to the companion paper [F. del Teso, J. Endal, and E. R. Jakobsen, SIAM T. Numer. Anal., 56 (2018), pp. 3611-3647] along with a more detailed discussion of the numerical methods included in our theory.