We formulate a numerical framework to model the structural patterns emerged from the long-term highly viscous tectonic flow for both two and three spatial dimensions by coupling the discontinuous Galerkin level set method with a finite element Stokes-like flow solver. Our formulation, implemented with adaptive mesh refinement near the material interface, allows for accurate interface capturing and automatic handling of topological splitting and merging. Compared to particle-in-cell family of methods, the level set formulation has the advantage of retaining information on the interface geometry, less memory requirement and the savings of computational expense on the two-way particle-mesh information transfer. Furthermore, our formulation discretizes the level set in the same finite element framework as the flow solver, thus enabling us to fully exploit the advantages of the finite element method such as the flexibility of mesh geometry and the ease of handling anisotropic materials. In order to track the finite deformation in the modelling domain, passive tracer particles are generated at and around locations of interest, whose deformation can be accumulated through arbitrary time interval within the total modelled time span, thus offering a fully dynamical approach for modelling non-steady and inhomogeneous structural patterns. The material distribution and the finite deformation pattern generated from the numerical model can be directly compared with the geological map patterns and the field structural analyses, thus offering the possibility of ground-truthing the modelling results by field evidence.