Modeling and simulation of large deformation problems and material flow has become an important area of research in geomechanics and geotechnical engineering. Typical applications include debris flows, post-failure of slopes, and pile driving. These problems are often characterized by the interaction of multiple, physically distinct materials on a hierarchy of spatial scales. Material interfaces (free surfaces, contact discontinuities) dynamically evolve, be generated, or coalesce. Moreover, due to its multiphase nature and internal structure, the mechanical behavior of soil can be either solid-like or fluid-like, depending on the underlying flow regime. In such situations the classical Lagrangian finite element method (FEM) suffers from serious losses in accuracy or even simulation break-off because material interfaces are tried to be maintained by the spatial discretization. Several alternative methods have been proposed during the past decades to meet the challenges, but most of them retain the Lagrangian formulation. In this work, a multi-material Eulerian formulation is presented in which the materials and their interfaces may move relative to the computational mesh. The implementation subdivides each solution step into a Lagrangian step, rezone step, and remap step. The Lagrangian step employs almost standard FEM technology with implicit time integration. Mesh elements intersecting with material interfaces contain a mixture of materials, each of which has its own set of constitutive equations. The interfaces are reconstructed and tracked by using the volume of fluid (VOF) method to achieve a reasonable accuracy of the flux-based, conservative remap step. This enables the simulation of very large deformations, shear and vorticity in a natural manner. Numerical examples concerned with geotechnical benchmark problems are presented which highlight the capabilities of the method and potential advantages over other approaches.