We present a fluid simulation method whereincompressibility is enforced through a holonomic constrainton the mass density. The method starts in aLagrangian particle formulation where the mass densityand other field quantities are represented by SmoothedParticle Hydrodynamics (SPH) kernel approximations.The density constraint is formulated as a regularizedmanybody constraint and is equivalent to very highsound speed. The system is integrated using a variationaldiscrete-time scheme, SPOOK, that includesconstraint regularization and stabilization. This constraintformulation of SPH enables systematic multiphysicsintegration, between rigid multibody physicsand fluids, where buoyancy falls out naturally. The fluidmodel results in a linear system of equations, whilemore general multiphysics systems result in a mixedlinear complementarity problem (MLCP) and we solvethese using iterative methods. The results demonstratenear perfect incompressibility, vastly improved stability,allowing for large time steps, and two orders of magnitudeimproved computational performance. Proof ofconcept is given for computer graphics applications andinteractive simulations.