SUMMARY The geometrical complexity of non-planar faults is known to affect not only fault slip rates but also crustal deformation and seismicity in the surrounding region, but implementing non-planar faults in numerical models has been challenging. We have developed a 3-D viscoelastoplastic finite element model to simulate long-term, steady-state fault slip along conceptual restraining bends and crustal deformation around them. By incorporating plastic yielding in both fault zones and the surrounding crust, this model can predict the quasi-steady state stress field and strain partitioning while avoiding pathological stress build-up encountered in traditional viscoelastic models. Here we detail the formulation of this model and explore major model parameters. The model results show that restraining bends tend to impede fault slip, increase shear stress in the surrounding crust and localize strain in a pair of belts off the fault zone. These effects are enhanced by high viscosity in the lower crust and upper mantle, high aspect ratio of bend width over fault length, fast relative motion of fault blocks and high fault friction. These model results may help explain the diffuse deformation surrounding the Big Bend of the San Andreas Fault and the more confined deformation around the Lebanon Bend of the Dead Sea Fault.