We formulate a numerical framework, in both 2d and 3d, to model the structural patterns emerging from creeping viscous flow typically encountered in long‐term ductile lithospheric deformation by coupling the discontinuous Galerkin level set method with a finite element Stokes‐like flow solver. The level set formulation has the advantage of retaining information on the interface geometry, decreased memory requirement and improved computational efficiency from the two‐way particle‐mesh information transfer compared to particle‐in‐cell methods. Furthermore, our formulation fully exploits the advantages of the finite element method (e.g., the flexibility of mesh geometry and the ease of handling anisotropic materials) by using a unified finite element framework. The novelty of our formulation is the capability to offer a fully dynamic approach for modeling structural patterns resulting from a tectonic flow that is non‐steady and inhomogeneous. The material distribution and the finite deformation patterns predicted from the numerical model can be directly compared with geological map patterns (e.g., lithological distribution at specified depths and on cross‐sections) and field structural analyses (e.g., foliation, lineation and strain patterns), thus offering the possibility of ground‐truthing the modeling results by field evidence. As examples for potential applications of our method, we apply our method to the modeling of a competent inclusion in simple shear flow, as well as a Rayleigh‐Taylor type density overturn. Our models demonstrate good agreement with previous 2d benchmark results and produce 3d lithological and deformation patterns comparable to field observations.