Previous works have demonstrated that inverting magnetotelluric (MT) data jointly with gravity data can synergize the high lateral resolution of gravity and the vertical resolution of MT. However, the existing joint stabilizers usually work for structured grids instead of unstructured ones that are more powerful for characterizing complex geology. Here, we utilize the local Pearson correlation coefficient (LPCC) for the joint inversion of gravity and MT data based on unstructured grids. We first establish a background mesh by discretizing the research area into virtual rectangular grids and then enhance the structured similarity between density and resistivity via the LPCC. Compared to existing joint constraints, our method is more flexible in solving multiscale joint inversions thanks to the adjustable subdomain size. The synthetic experiments show that the joint stabilizer can recover the subsurface targets at a higher resolution, especially for gravity data, than the standalone inversions. This method is further applied to the joint inversion of gravity and MT data from the Yellowstone area and the inverted density and resistivity models are structurally consistent. Then, based on the inverted subsurface structures and incorporating the existing research, we infer that the two inverted zones with low density and low resistivity correspond to the partially molten rhyolitic and basaltic, respectively. The proposed joint constraint can be further extended to the inversion of other geophysical data.