Partial differential equations with inequality constraints in the form of variational inequalities arise in variety of applications. In particular, we mention biofilm and phase transition models, in which the variational inequality appears in one of the equations of a larger system. In this talk, we consider the scalar model equation, and focus on the theoretical background as well as implementation of finite element and finite difference algorithms. In particular, we demonstrate optimal convergence. In addition, we consider and compare two solvers for discrete elliptic and parabolic variational inequalities with pointwise inequality constraints. The first solver is the classical relaxation method, and the second modern one is from the class of semismooth Newton methods.