Phase field models are a widely used approach to describe physical processes that are characterized by thin interfacial regions between large almost homogeneous domains. Important application areas of phase field models are transition processes of the state of matter and the separation of alloys. A fundamental property of these models is, that the transition and separation of phases is driven by a double-well potential with distinct minima for the different phases. Already the pioneering work of Cahn and Hilliard used a temperature dependent logarithmic potential that is differentiable with singular derivatives. If the temperature tends to zero it degenerates to the non-differentiable obstacle potential. The goal of this thesis is to develop methods for the efficient numerical solution of such equations that are also robust for nonsmooth potentials and anisotropic surface energies. These methods are derived for the Cahn-Hilliard equation that are prototypic for a multitude of such models. The main result of the thesis is the development of a fast iterative solver for nonlinear saddle point problems like the ones that arise from finite element discretization of Cahn-Hilliard equations. The solver relies on a reformulation of the problem as dual minimization problem whose energy functional is differentiable. The gradient of this functional turns out to be the nonlinear Schur complement of the saddle point problem. Generalized linearizations for the Schur complement are derived and used for a nonsmooth Newton method. Global convergence for this 'Schur Nonsmooth Newton' method and inexact versions is proved using the fact the equivalence to a descent method for the dual minimization problem. Each step of this method requires the solution of a nonlinear convex minimization problem. To tackle this problem the 'Truncated Nonsmooth Newton Multigrid' (TNNMG) method is developed. In contrast to other nonlinear multigrid methods the TNNMG method is significantly easier to implement and can also be applied to anisotropic problems while its convergence speed is in general comparable or sometimes even faster. Numerical examples show that the derived methods exhibit mesh independent convergence. Furthermore they turn out to be robust with respect to the temperature including the limiting case zero. The reason for this robustness is, that all methods do not rely on smoothness but on the inherent convex structure of the problems.