Abstract: This article describes a Lagrangian finite element method that simulates the behavior of liquids and solids in a unified framework. Local mesh improvement operations maintain a high-quality tetrahedral discretization even as the mesh is advected by fluid flow. We conserve volume and momentum, locally and globally, by assigning to each element an independent rest volume and adjusting it to correct for deviations during remeshing and collisions. Incompressibility is enforced with per-node pressure values, and extra degrees of freedom are selectively inserted to prevent pressure locking. Topological changes in the domain are explicitly treated with local mesh splitting and merging. Our method models surface tension with an implicit formulation based on surface energies computed on the boundary of the volume mesh. With this method we can model elastic, plastic, and liquid materials in a single mesh, with no need for explicit coupling. We also model heat diffusion and thermoelastic effects, which allow us to simulate phase changes. We demonstrate these capabilities in several fluid simulations at scales from millimeters to meters, including simulations of melting caused by external or thermoelastic heating.