## 1 Introduction
To more accurately characterize the anisotropic properties of subsurface media, in 1996, [PERSON] et al. analyzed the response of a vertical electric dipole transmitter in layered anisotropic formations [13]. In 2005, [PERSON] et al. conducted detailed studies on the interpretation of deep, complex geological anomalies, warning that the continued use of isotropic models for electromagnetic data interpretation could undermine the reliability of the results [14]. In 2011, [PERSON] et al. proposed a parallel, goal-oriented adaptive finite element method for rapidly solving 2.5D controlled-transmitter audio-frequency magnetotellurics (CSAMT) and 2D magnetotellurics (MT) modeling problems, which yielded high-precision solutions and is particularly suited for investigating the response characteristics of CSAMT in media with resistivity anisotropy [15]. Furthermore, in 2017, [PERSON] et al. derived and implemented a three-dimensional staggered-mesh finite difference numerical simulation method for CSAMT with axial resistivity anisotropy [16]. Given that subsurface media often exhibit varying resistivity in different directions, posing challenges for the interpretation of electromagnetic exploration data [17], the development of forward modeling algorithms capable of effectively simulating three-dimensional anisotropic resistivity features is of significant importance. In the numerical solution of electromagnetic response problems, the finite element method (FEM), the finite difference method (FDM), and the finite volume method (FVM) are the primary numerical techniques [18]. However, the reliance of the FDM on structured meshes limits its application in complex geological models, potentially leading to decreased computational accuracy and erroneous interpretations of predictive parameters [19]. Therefore, the FEM has been increasingly adopted in electromagnetic numerical simulations due to its ability to accurately model complex geological structures using unstructured meshes and enhance the precision of numerical solutions through high-order basis functions [18]. The FEM originated in the field of structural analysis. Although [PERSON] et al. provided early mathematical treatment of the method in 1943 [20], it was not until 1968 that the FEM was applied to solve electromagnetic problems [21]. The FEM can be categorized into nodal finite elements and edge-based finite elements. In the nodal finite element method, the degrees of freedom are associated with the nodes of the elements, whereas in the edge-based finite element method, the degrees of freedom are assigned to the edge-based parts within the elements, with the edge-based elements constructed and based on Whitney spaces [22]. The nodal elements enforce both the tangential and normal continuity of edge-based or scalar fields at interface boundaries [23]. However, in edge-based elements, only tangential continuity is enforced at interfaces, aligning with the characteristics of electromagnetic field problems [24]. In electromagnetic analysis, the required continuity is directly achieved using edge-based elements, while in nodal elements, it is obtained through potential formulations [25]. Therefore, when using edge-based elements, field variables can be directly computed, whereas for nodal elements, field variables must be derived from potentials during post-processing [26]. In the finite element method, the most commonly used mesh types are hexahedral [27,28] and tetrahedral meshes [29, 30, 31]. Among these, tetrahedral meshes are the most prevalent in three-dimensional unstructured meshes, as tetrahedra are the simplest shapes in 3D space and can fill any spatial region [32]. Traditional regular-mesh methods have resolution limitations constrained by mesh density, making it challenging to accurately depict the complex structures of subsurface media [33,34]. Finite element methods based on unstructured meshes can flexibly adapt to the intricate structures of subsurface media, effectively reducing computational complexity and time [35].