Hybrid mesh for magnetotelluric forward modeling based on the finite element method

Unstructured tetrahedral grids have been applied in magnetotelluric (MT) forward modeling using the finite element (FE) method because of their adaptability to complex anomalies. However, high-quality results require an extreme refinement of the near-surface area, which leads to excessive meshes and an increased degree of freedom (DoF) of the governing equation of the finite element system. To reduce the computational cost, we have developed a hybrid mesh based on triangular prisms and tetrahedrons. The required elements in the near-surface area are reduced because the quality of the triangular prism is not limited by the element aspect ratio. The deep area is discretized by tetrahedral elements to ensure the flexibility of the unstructured grids. The superiority of this hybrid mesh has been tested on a layered model, the DTM1 model and terrain relief models. The results show that the modeling efficiency has been improved, especially for high-frequency data. The accuracy of the modeling using the hybrid mesh is significantly higher than that of the tetrahedral mesh with a similar DoF. Usage of the hybrid mesh can be easily adapted to complex geoelectric models with strong terrain fluctuations, which requires less computational cost than using conventional unstructured elements.

Unstructured tetrahedral grids have been applied in magnetotelluric (MT) forward modeling using the finite element (FE) method because of their adaptability to complex anomalies. However, high-quality results require an extreme refinement of the near-surface area, which leads to excessive meshes and an increased degree of freedom (DoF) of the governing equation of the finite element system. To reduce the computational cost, we have developed a hybrid mesh based on triangular prisms and tetrahedrons. The required elements in the near-surface area are reduced because the quality of the triangular prism is not limited by the element aspect ratio. The deep area is discretized by tetrahedral elements to ensure the flexibility of the unstructured grids. The superiority of this hybrid mesh has been tested on a layered model, the DTM1 model and terrain relief models. The results show that the modeling efficiency has been improved, especially for high-frequency data. The accuracy of the modeling using the hybrid mesh is significantly higher than that of the tetrahedral mesh with a similar DoF. Usage of the hybrid mesh can be easily adapted to complex geoelectric models with strong terrain fluctuations, which requires less computational cost than using conventional unstructured elements.
Unstructured grids with tetrahedral elements are suitable for dividing complex underground anomalous bodies in terms of topography and bathymetry because they can fit arbitrary shapes of geological bodies well [1][2][3][4][5] . This type of grid has become a useful tool for discretizing the geoelectric model in the numerical simulation of three-dimensional geophysical electromagnetic (EM) fields using the finite element method (FEM) [6][7][8] and has been widely used in magnetotelluric (MT) forward modeling [9][10][11] . However, unstructured grids often need to be refined into extremely small elements to avoid ill-conditioning 12,13 . The refinement often results in too many redundant elements and reduces the computational efficiency of the FEM forward modeling.
To reduce the number of required elements, a scheme of nonconforming grids has been used to discretize the geoelectric model 14,15 . This kind of regular grid, whose element quality is not restricted by the aspect ratio of the element, has a variable scale and stable local refinement capabilities 16,17 . It helps to reduce the degrees of freedom (DoF) and is suitable for large-scale 3-D geophysical EM forward problems, such as MT forward modeling. A multiresolution approach suggested by Cherevatova et al. 13 allows grid refinement only in the horizontal directions and keeps the degree of refinement in the vertical direction constant. The stretching grid obtained from this multiresolution approach can help to improve calculation accuracy with fewer elements in the near-surface area . In addition, the longer the element is stretched, the more elements are saved, and the decrease in numerical accuracy caused by element stretching can be compensated by a higher-order element processing scheme 18,19 . However, some elements of the stretched grid with hanging nodes cannot share just one whole edge or face 13,14 , and the convenience and flexibility of mesh generation are not as good as traditional tetrahedral elements.
Triangular prism elements can be used to build the stretched grid in the near surface area 20 . However, the accuracy of the numerical solution will be reduced when the stretched grid extension is too long. In this research, second-order shape functions are used to discretize near-surface prismatic elements. The deep region is still freely discretized by the tetrahedral elements. Because triangular prism elements and tetrahedral elements can share a triangular surface, we developed a hybrid mesh with second-order shape functions containing two kinds of elements.
This paper is structured as follows. We first introduce our governing equations for forward modeling, the A -ϕ system with gauge fixing. Then, we describe how to create a hybrid mesh with triangular prism elements and tetrahedral elements and how we apply it to FEM forward modeling. We use a layered model to preliminarily Governing equation. The electric field can be decomposed into the magnetic vector potential A and electric scalar potential ϕ . This can reduce the weakness of the curl-curl electric equation [21][22][23][24] , accelerate the speed of solving the FE system with an iterative solver and obtain higher numerical accuracy. The relation of the potentials and fields can be represented as: Substituting Eq. (4) into Eq. (3) results in the following equation 25 : In conjunction with J = σE, by substituting Eq. (4) into the Gaussian divergence equation, the current divergence is zero and can be expressed as follows: When both the curl and divergence of A are given, A is uniquely defined up to a constant. Therefore, to obtain a unique solution, we add the Coulomb gauge condition 26,27 : The Coulomb gauge is imposed to strengthen the divergence limit and enhance the uniqueness of the FE equations. On the basis of Eq. (6) derived from Ampere's rule, we add a Lagrange multiplier term to eliminate the external current density divergence 21,23 . Equations (6) and (7) can be rewritten as follows: In Eq. (10), −iωσ is constant, so the divergence of the magnetic vector potential is specified. For the geophysical EM field propagation problem, we use homogeneous Dirichlet boundary conditions in the truncated model boundary.

Solving a linear system of equations.
We use the Galerkin method to discretize Eq. (9). The discretized expression is as follows: where N denotes the vector shape function. Ã , φ , and ˜ are the approximate solutions of A , ϕ , and , respectively. is all domains of the calculation model. Ŵ and γ are the outer and inner boundaries of the mesh, respectively. Similarly, the discrete forms of Eqs. (7) and (8) are obtained as follows 23 : where N denotes the scalar shape function. We use natural electromagnetic field sources to form the right side of the system. When the polarization direction of the magnetic field source is the y direction, the air top interface loads the magnetic field source H 0 = (0, 1, 0), the boundary condition of the outer boundary parallel to the The calculation area is subdivided into units of smaller subdomains. The vector and scalar potential of each element in the grid are represented by piecewise polynomial basis functions, and the expressions are listed as follows: where N j denotes the vector shape function on the jth edge in the element, N j denotes the scalar shape function of the jth node in the element, and n edge and n node represent the number of edges and nodes in the element, respectively. Combining Eqs. (11) to (13) and Eqs. (16) to (18), the matrix form of the system equations after finite element discretization is obtained as follows: where Ŵ air isthe air top interface. i, j = 1, ..., n edge and l, k = 1, ..., n node . We use the MUMPS direct solver to solve the A-ϕ system equations.

FEM with a hybrid mesh
Mesh generation. The geoelectric model can be divided into near-surface areas and deep areas. The nearsurface area is discretized by triangular prism elements, which can be seen as a boundary layer near the surface (Fig. 1a). The deep area is discretized by tetrahedrons, which lie below the boundary layer. The triangular prism and the tetrahedron can share a triangular face, which perfectly couples the triangular prism mesh and the tetrahedral mesh (Fig. 1b). We use the mesh generator in COMSOL to build the mesh, which is based on the Delaunay algorithm [28][29][30][31][32][33] . Other software capable of generating tetrahedrons, prisms, hexahedrons and other meshes can also be used to generate hybrid meshes, such as Gmsh 34 .
For the near-surface area of the geoelectric model, we use Delaunay triangulation algorithms to freely divide the Earth's surface into triangular meshes and then extend the triangular meshes in the vertical direction to be a top layer comprised of triangular prisms. The extending distance is approximately 1 to 2 times the skin depth. The near-surface layer is approximately vertical. The change in the underground electromagnetic field can be measured by the skin depth. After passing our experiment, the vertical size of each subdivision layer is set to be less than 1/3 of the skin depth.
The areas adjacent to the triangular prism layers are the deep area below it and the air layer area above it. These adjacent areas are discretized into free tetrahedrons, which couple the upper and lower interfaces of the triangular prism layer through the triangular plane. The sizes of the tetrahedrons are based on the sizes of the triangles on the coupling surface and gradually grow away from the surface. Element analysis of triangular prism and pyramid. The shape functions of the triangular prism elements ( Fig. 2) consist of vector edge-shaped functions and scalar node-shaped functions 36 . A triangular prism www.nature.com/scientificreports/  Shape functions of second-order triangular prism element. The quadratic triangular prism element is shown in Fig. 4. It has a total of 15 nodes, eighteen variables along edges and ten on faces in the element. V i (i = 1, 2, ..., 28) is the variable in the element. The element illustrated in Fig. 4 has fifteen scalar node-based shape functions 37 .
The scalar node-based shape functions of the triangular prism element of corner nodes are: The scalar node-based shape functions of the mid-edge of the triangle are: The scalar node-based shape functions of the mid-edge of the rectangle are: Ahagon et al. presented that the vector edge-shaped functions are derived from those scalar node-based shape functions 38,39 . Their relationship is: where the edge i-j is associated with the node j shape function. According to Eqs. 21-24, the 2nd vector edgeshaped functions of triangular prisms can be derived easily. Taking node 1, node 7 and node 11 as examples:  www.nature.com/scientificreports/ The other 2nd vector edge-shaped functions can be deduced in the same way as above. We can divide the vector edge-shaped functions into several categories as follows: The vector edge-shaped functions for the bottom and top triangle edges can be expressed as: The vector edge-shaped functions for the rectangle edges (volume vector functions) can be expressed as: The vector edge-shaped functions for the bottom and top triangle faces can be expressed as: The vector edge-shaped functions for the rectangle faces can be expressed as:

Advantages of hybrid mesh
We use the numerical example of a layered model to illustrate the superiority of triangular prism elements in MT forward modeling with FEM. The apparent resistivity and phase are the responses, which help us to analyze the accuracy of the calculation.
To simplify the analysis, we use the highest frequency to calculate the smallest skin depth, which is used to stretch the triangular prism layer in the near-surface area. Although this size is sufficient for the calculation of an EM field with lower frequencies, it does not affect the analysis of the superiority of the triangular prism to the tetrahedron.
Meshing the three-dimensional layered model. The layered model is divided into three layers (Fig. 5).
The resistivities of the first, second and third layers are 200 · m, 1000 · m, and 200 · m, respectively. The thicknesses of the first layer and the second layer are both 0.5 km. The measuring points are located from −1200 m to 1200 m along the X direction at Y = 0 m and Z = 0 m, with a distance of 300 m. The size of our entire calculation model is 20 km × 20 km × 70 km.
The sizes of mesh elements can be divided into horizontal sizes and vertical sizes and are controlled by the maximum unit size (MES) and growth rate. The control of these sizes mainly lies in the discretization of prismatic (27) ∇N 11 = 2ξ 1 ξ 2 ∇ζ + 2(1 + ζ )ξ 1 ∇ξ 2  For the near-surface area, the horizontal sizes are controlled by the triangles used for meshing the ground. We use the triangular MES (TriMES) and triangular growth rate (TriGR) to generate these triangles. We use the prism growth rate (PriGR) to control the growth multiple of the vertical size of the triangular prism and control the total thickness of the prismatic layer.
For the deep area, the sizes of the top layer of the tetrahedrons, which are coplanar with the coupling plane, are controlled by the size of the coupling triangles. The sizes of the other tetrahedrons are controlled by the tetrahedral growth rate (TetGR).
We use four types of meshes (Table 1) to discretize the three-dimensional layered model. Mesh A is a hybrid mesh that is divided into two parts: the upper part is a triangular prism mesh, and the lower part is a tetrahedral mesh. The horizontal size is determined by the size of the triangular element at the top. The TriMES is 1500 m, and the TriGR is 1.3. The vertical size of the hybrid mesh is determined by the vertical size of the triangular prism and the tetrahedron at the same time. The vertical size of the first prismatic layer is 8 m, and the PriGR is equal to the TetGR. The total number of prismatic layers is 10. Below the prismatic layers, the vertical size of the tetrahedral element is controlled by the TriMES and TetGR. The TetGR in hybrid mesh A is 1.5.
Meshes B1, B2, and B3 are tetrahedral meshes, and their granularities are different. Mesh B1 is the coarsestgrained mesh. The horizontal size of the mesh B1 element is the same as that of the triangular prism element, and the TriMES is 1500 m. The vertical size of the mesh B1 element in the first layer is approximately 500 m under  Table 1. Parameters of different meshes for the three-dimensional layered model.l layered model. MES means maximum unit size. The horizontal size is controlled by the triangular MES (TriMES) at the Earth's surface and triangular growth rate (TriGR). The TriGR represents the maximum element growth rate of the triangular element size. Vertical size is controlled by the thickness of the prismatic layer and tetrahedral growth rate (TetGR). The TetGR represents the maximum element growth rate of the tetrahedral element size. The prismatic growth rate (PriGR) represents the vertical prismatic element growth rate of the prism layer thickness The TriMES and the first layer thickness of the B2 mesh are 100 m and approximately 65 m, respectively. We compared the mesh refinement at measuring points (Fig. 6) and the grid diagram of the cross section (Fig. 7) of the four grids.
Higher-order shape function. MT responses are calculated in the frequency range of 1E5 Hz to 1E−1 Hz with hybrid mesh 20 of different orders. In total, 26 frequency points were calculated, of which 10 frequency points were equally spaced between 1E5 Hz and 1E4 Hz, and the other frequency points were equally spaced from 1E4 Hz to 1E−1 Hz. The apparent resistivity and phase change with frequency are shown in Fig. 8. When the shape function is first order, the maximum absolute errors of the apparent resistivity and phase are 27.8510% and 6.8170%, respectively. When the shape function is second order, the maximum absolute errors of the apparent resistivity and phase are 2.0968% and 0.4063%, respectively. All the computations are implemented in MAT-LAB 2019 with an Intel Core i5-10400 CPU @ 2.90 GHz and 32 GB RAM.
Mesh quality. MT responses are calculated in the frequency range of 1E4 Hz to 1E−1 Hz. The apparent resistivity and phase change with frequency are shown in Fig. 9. The accuracy of the apparent resistivity and phase calculated by different meshes is quite different in this frequency range. We use the root mean square error (RMSE) to measure the overall error of the apparent resistivity and phase, and the reference solution uses the semianalytical solution of the 1D layered model. For a mesh containing only tetrahedrons, the smaller the mesh granularity is, the smaller the calculated response misfit. For the layered model, the vertical size of the element is a key factor affecting the calculation accuracy. Considering the element quality, the horizontal size of the element needs to be consistent with its A hybrid mesh with a triangular prism can overcome this limitation. Even if the aspect ratio of the triangular prism element is greater than 100 times, high accuracy can be obtained. When the accuracy is not limited by the aspect ratio, the number of required elements is greatly reduced, and the required DoF number is also reduced. The DoF, mean absolute error (MAE) misfit and NE of the response results are all shown in Table 2. Table 2 shows the comparison of the FE system element number, DoF, calculation time and numerical solution accuracy information of different grids. A higher order is helpful to improve the accuracy of the numerical solution for the prismatic element. A comparison of the calculation results of mesh A and mesh B1 shows that mesh A has a greater DoF and more accurate numerical solutions. A comparison of the calculation results of meshes A and B2 shows that even if the DoF number of the tetrahedral mesh is eight times that of the hybrid mesh, the accuracy of the hybrid mesh is much higher than that of the tetrahedral mesh. A comparison of the results of meshes A and B3 shows that when the DoF number of the tetrahedral mesh is 56 times that of the hybrid mesh, the accuracy of the tetrahedral mesh can be at the same level as that of the hybrid mesh. This shows that the mesh type not only affects the size and DoF number of the system matrix but may also affect the solution properties of the system matrix.

System matrix quality.
The sparsity distribution of the system matrix of different meshes can reflect that the system matrix of the hybrid grid has better solution properties. For tetrahedral grids, as the number of DoFs increases, the bandwidth of the system matrix gradually weakens. The DoF of the hybrid grid is between the number of DoFs of meshes B1 and B2, but the bandwidth of its system matrix is notably better than that of the other two tetrahedral meshes.

Numerical example: DTM1
The hybrid mesh mainly uses triangular prism elements to process near-surface areas, which are often related to the calculation accuracy of high-frequency data. To verify the role of the hybrid mesh in improving the accuracy of the high-frequency data solution, we increase the complexity of the model on the basis of the three-dimensional layered model and select the DTM1 model [41][42][43][44] to show the advantages of the prism elements. A survey line is selected to record the responses of apparent resistivity and phase (Fig. 10). In the public DTM1 test data  We set up four different meshes in total, namely, a hybrid mesh and three tetrahedral meshes. The parameters of different meshes are shown in Table 3. In the hybrid mesh, the TriMES, which controls the size of the horizontal element, is set to 2000 m, and the thickness of the first layer of the triangular prism element is 318 m. In the tetrahedral mesh, the TriMES gradually changes from 2000 m to 500 m. The size of the tetrahedral element on the first layer is approximately equal to the TriMES on the surface. In the four meshes, the growth rate of both the triangular element and the tetrahedral element, TriGR and TetGR, is set to 1.5. The system equations established by TriMES's 1000 m tetrahedral mesh and TriMES's 2000 m hybrid mesh have the same DoF. The tetrahedral mesh of 500 m TriMES can be regarded as a dense tetrahedral mesh.
The calculation results of the apparent resistivity and phase at frequencies of 0.1 Hz, 1 Hz, and 10 Hz are shown in Fig. 11. When the frequency is 10 Hz (Fig. 11a, d), the denser the tetrahedral mesh is, the smaller the oscillation of the apparent resistivity curve and the phase curve, but neither is as smooth as the response curve of the hybrid mesh. When the frequency is 1 Hz (Fig. 11b, e), the oscillation of the response curves of the tetrahedral mesh is reduced, but there is still a gap between the smoothness of the response curve of the hybrid mesh. Up to the frequency of 0.1 Hz (Fig. 11c, f), the smoothness of the response curves of the tetrahedral mesh is equivalent to that of the hybrid mesh. When calculating the MT high-frequency electromagnetic field, the discrete equation using the hybrid mesh is smoother, and the calculation result has higher accuracy.
When the horizontal size of the hybrid mesh and the tetrahedral mesh are both approximately 2000 m, because the triangular prism element in the hybrid mesh can reduce the vertical size without reducing the mesh quality, the hybrid mesh has more DoFs and smoother responses. When we refine the tetrahedral mesh, we reduce its horizontal size to approximately 1000 m and set the vertical size to approximately 1000 m to ensure that the number of DoFs of the tetrahedral mesh is similar to that of the hybrid mesh. The apparent resistivity and phase curves calculated by this dense tetrahedral mesh still oscillate significantly. The tetrahedral mesh www.nature.com/scientificreports/ is further refined, and the horizontal and vertical sizes are both set to approximately 500 m. At this time, the number of DoFs is nearly twice that of the hybrid mesh, but the calculation result is still not as stable. This shows that the hybrid mesh is indeed more efficient than the pure tetrahedral mesh in calculating high-frequency data.

Numerical example: sloped terrain models
Due to the volume effect of the electromagnetic method, the calculation of the high-frequency electromagnetic response is more difficult when the terrain is undulating. We adjust the change in the terrain undulation angle to analyze the improvement in the calculation efficiency of the hybrid mesh under different terrain conditions when the frequency is 1E4 Hz (Fig. 12).  www.nature.com/scientificreports/ We set up three kinds of meshes, namely, a hybrid mesh, general tetrahedral mesh and dense tetrahedral mesh. The parameters of different meshes are shown in Table 4. The hybrid mesh and general tetrahedral mesh have the same horizontal element size, and their TriMESs are both 120 m. The horizontal element size of the dense tetrahedral mesh is smaller, and its TriMES is 50 m. The thickness of the first layer of triangular prism elements in the hybrid grid is 15 m. The vertical size of the first tetrahedron in the other two tetrahedral meshes is approximately equal to its horizontal size. Among these meshes, the mesh growth rates of triangles and tetrahedrons are both 1.5. The system matrix produced by the hybrid mesh and the tetrahedral mesh with the same horizontal element size has a similar number of DoFs, while the system matrix produced by the tetrahedral grid with the smallest horizontal element size has nearly twice the number of DoFs.
When the tilt angle changed, we calculated the responses of the xy component (Fig. 13) and the responses of the yx component (Fig. 14) with various types of meshes. For the xy component (Fig. 13), the accuracy of the resolution of the apparent resistivity accuracy is similar, and the calculation results of the general tetrahedral  www.nature.com/scientificreports/  The inconsistent point between the yx component and the xy component is that the response of the yx component has more dramatic oscillation (by comparing Figs. 13 and 14). The elevation on the slope changes along the X direction, which brings greater misfit to the yx component. It can be suggested that there may be a large misfit in the direction of large topographic relief.

Numerical example: terrain relief model
In the high-frequency range, the hybrid mesh can obtain a high-precision EM response regardless of whether it is in flat terrain or undulating terrain with slopes. To further verify the adaptability of the hybrid mesh to undulating terrains, we design a three-dimensional geoelectric model with complex undulating terrain, which is referred to as the real terrain undulations (Fig. 15), and compare the accuracy of the MT response calculated by the hybrid mesh and the tetrahedral mesh. Here, the frequency is 1E4 Hz.
We set up three kinds of meshes, namely, a hybrid mesh, general tetrahedral mesh and dense tetrahedral mesh. The parameters of different meshes are shown in Table 5. The hybrid mesh and general tetrahedral mesh have the same horizontal element size, and their TriMESs are both 100 m. The horizontal element size of the dense tetrahedral mesh is smaller, and its TriMES is 20 m. The TriGRs in the three kinds of meshes are all 1.5. The thicknesses of the first layer in the three kinds of meshes are 15 m, 80 m and 20 m. To make the system matrix generated by the hybrid mesh and general tetrahedral mesh have similar DoF numbers, we set the growth rates of the hybrid mesh and general tetrahedral mesh to 1.5 and 1.34, respectively. The TetGR in the dense tetrahedral mesh is 1.5.
The MT response results of the xy and yx components on the two survey lines are shown in Figs. 16 and 17. On survey line 1 (Fig. 16), the calculation result of the general tetrahedral mesh is the most inaccurate. The result of the hybrid mesh is closer to the result of the dense tetrahedral mesh. For the phase responses, the results of the dense tetrahedral mesh have obvious oscillations. This shows that the hybrid mesh can obtain stable calculation results while building a system matrix with a lower number of DoFs.
The response results on survey line 2 (Fig. 17) also show that when the DoF numbers are similar, the calculation results of the hybrid grid have higher accuracy than the calculation results of the general tetrahedral mesh. The dense tetrahedral grid can obtain more accurate calculation results, but the system matrix DoF number is nearly 2.6 times that of the hybrid grid. Additionally, compared with the phase results of the hybrid grid, the phase results of the dense tetrahedral grid also show stronger oscillations.

Conclusions
We have proposed a hybrid mesh strategy for MT forward modeling with the FEM, which contains triangular prism elements and tetrahedral elements. The use of triangular prism elements to discretize the near-surface area can significantly reduce the amount of required NE, which is caused by mesh quality constraints for tetrahedrons. It can improve the computational efficiency of high-frequency data.
Different geoelectric models are used to analyze the accuracy and efficiency of the hybrid mesh in calculating the MT responses. The three-dimensional layered model and the DTM1 model verified their accuracy and  www.nature.com/scientificreports/ applicability to general types of geoelectric models. A hybrid mesh can build an FE system matrix with fewer DoFs and obtain high accuracy in the solutions. Its superiority has also been verified in the example of terrain relief models. The triangular prism elements have the flexibility of horizontal refinement, which is controlled by the triangular elements. It also has stability in terms of the element quality because the rectangular elements are not affected by the aspect ratio. We believe that application of this type of hybrid mesh is not limited to MT forward modeling but can also be useful to other geophysical electromagnetic calculation problems. www.nature.com/scientificreports/

Data availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Appendix A
Shape functions of triangular prism element. The vector edge-shaped functions of the triangles are chosen as the Whitney-1 form basis function 46,47 and can be represented as W . The corresponding edge is determined by the node number indicated by the subscript W . Two nodes at both ends determine the edge represented by W . W can be expressed as: where ξ i (i = 1, 2, 3) is a plane triangle nodal shape function established by Lagrangian interpolation. The Whitney-1 shape function represents the area coordinates of the triangle element. ∇ξ i (i = 1, 2, 3) represents the gradient vector of the triangular element. The calculation of these shape functions can directly refer to the element analysis of triangular elements in the two-dimensional case, but the three-dimensional shape function of the triangular prism cannot refer to the two-dimensional case. The local coordinate ξ (Fig. 2) in the height direction needs to be added to the shape function. The vector edge-shaped functions for the bottom triangle edges (bottom vector functions) can be expressed as: The vector edge-shaped functions for the top triangle edges (top vector functions) can be expressed as: The vector edge-shaped functions of the rectangle edges (volume vector functions) can be expressed as: