We developed a novel method to efficiently model complex fault (zone) across aseismic and seismic cycles, where Intersecting faults are dealt with cross-link constraint method, and off-fault inhomogeneity is dealt with Eshely’s inclusion method.