We present a newly developed algorithm for three-dimensional (3D) magnetotelluric model based on the vector finite element method. In this paper, unstructured grids which are composed of irregular tetrahedrons are used in our finite element model, which can be refined locally and adaptively according to the complex geometry of computational domain or subsurface structure. For obtaining more accurate solutions, secondary field rather than total field is numerically computed, which makes the errors limited to relatively small secondary field. Traditional node-based finite element method does not satisfy the condition that the normal electrical field is discontinuous at the interface of electric separatrix and the electrical current density is divergence-free throughout the regions without source, which obviously violates Maxwell equations. In order to overcome these drawbacks of node-based finite element, vector finite element method is employed to solve the secondary field-based partial differential equation. Moreover, in this study, the heterogeneous permeability is taken into consideration in our algorithm, as a consequence, which can deal with heterogeneous permeability model. The accuracy of our approach is verified by comparing with previously published numerical simulations of a COMMEMI-3D model. The advantages of our approach are also illustrated by the numerical simulations of model with arbitrary topography and complicated anomalous body. In addition, the simulation results of the irregular anomaly and the complex model of arbitrary terrain are both ideal results. It proves that the effect of 3D terrain is more serious and complex than that of two-dimensional (2D) terrain, and dealing with 3D models or data with 2D algorithm may bring in large errors. In the area where the magnetic permeability is abnormal, the magnetic permeability has an important influence on the numerical simulation results, and the magnetic permeability must be treated as an independent parameter in the magnetotellurics survey.