Development of tectonomagnetic modeling on the basis of the piezomagnetic effect is reviewed for the period since the early 1990's. First, the basic theory is briefly summarized, in which the representation theorem or the surface integral representation for the piezomagnetic potential and the Green's function method are presented. In the 1990's, several field observations in earthquakes and volcanoes were interpreted with the aid of analytic solutions based on the Green's function method. A general formula was developed for an inclined rectangular fault with strike-slip, dip-slip and tensile faulting. The surface integral method has been applied to 2D and 3D models, as well as to fault models in the inhomogeneously magnetized crust. When the magnetic field is measured within a bore hole, the effect of magnetic poles around the hole should be taken into account. As a result, tectonomagnetic signals are much enhanced in a bore hole compared with on the ground surface. Finally, piezomagnetic field changes associated with the Parkfield fault model are introduced and the new aspect of the model is discussed.