We present a fundamental solution‐based finite‐element (FE) method to homogenize heterogeneous elastic medium, that is, fault zone, under static, and dynamic loading. This method incorporates Eshelby’s strain perturbation into FE weak forms. The resulting numerical model implicitly considers the existence of inhomogeneity bodies within each element, without introducing additional degrees of freedom. The new method is implemented within an open‐source FE package that is applicable to alternating seismic and aseismic cycles. To demonstrate this method, we modify a dynamic fault‐slip problem, hosted at Southern California Earthquake Center (SCEC), by introducing a fault zone that contains different microstructures than the host matrix. The preliminary results suggest that the fault‐zone microstructure orientation has effects on fault slip, seismic arrivals and waveform frequency contents.