Geomechanical restoration methods to date rely on considering the rock properties as fully elastic and the faults as frictionless surfaces. However, the linear elastic hypothesis is valid only for very small deformations and assumes reversibility and conservation of energy, which is not valid for rocks at geological time scales. Moreover, considering faults as frictionless surfaces ignores the mechanical processes associated with fault development and evolution. In order to achieve geomechanical restoration of large deformations at time-scales of millions of years, we use non-elastic laws, and we model fault displacement not as a localised surface but as a shear zone with specific material parameters. At this time-scale, rocks can behave as highly viscous fluids, so we use creeping flow to simulate their behaviour, similarly to what is done for Earth mantle simulations. The Stokes equations describing this behaviour generally involve significant deformations, which are difficult to handle numerically with Lagrangian methods such as those used in geomechanical restoration. Therefore, we use adaptative meshes such as Arbitrary Lagrangian Eulerian (ALE), which have shown their value for simulating large deformations, removing the necessity for regular remeshing when solving elastic equations with Finite Elements on tetrahedral meshes. We implement these methods in a new restoration scheme based on different physical considerations. To test and benchmark our method, we use it to perform forward and backward simulations on an analogical model. This model was previously restored using elastic laws and free-slip faults. It allows us to compare the two methodologies and their advantages and downfalls.