We reproduce the general behavior of complicated bubble and droplet motions using the variational level set formulation introduced by the authors earlier. Our approach here ignores inertial effects; thus the motion is only correct as an approximation for very viscous problems. However, the steady states are true equilibrium solutions. Inertial forces will be added in future work. The problems include: soap bubbles colliding and merging, drops falling or remaining attached to a (generally irregular) ceiling, and liquid penetrating through a funnel in both two and three dimensions. Each phase is identified with a particular "level set" function. The zero level set of this function is that particular phase boundary. The level set functions all evolve in time through a constrained gradient descent procedure so as to minimize an energy functional. The functions are coupled through physical constraints and through the requirements that different phases do not overlap and vacuum regions do not develop. Both boundary conditions and inequality constraints are cast in terms of (either local or global) equality constraints. The gradient projection method leads to a system of perturbed (by curvature, if surface tension is involved) Hamilton-Jacobi equations coupled through a constraint. The coupling is enforced using the Lagrange multiplier associated with this constraint. The numerical implementation requires much of the modern level set technology; in particular, we achieve a significant speed up by using the fast localization algorithm of H.