An Implicit Adaptive FDTD Mesh Generation Techniquebased on Tetrahedrons

Weiran Zhang $^{1}$ , Zikun Xu $^{1}$ , Huaiyun Peng $^{2}$ , Juan Chen $^{1*}$ , and Chunhui Mou $^{1}$

¹School of Information and Communication Engineering
Xi’an Jiaotong University, Xi’an, 710049, China
1538162321@qq.com, *chen.juan.0201@mail.xjtu.edu.cn
²National Key Laboratory of Electromagnetic Environment
Qingdao, 266107, China

Submitted On: February 21, 2023; Accepted On: October 12, 2023

ABSTRACT

A novel implicit adaptive FDTD mesh generation method based on tetrahedrons is proposed in this paper. According to the vertex coordinates of tetrahedrons which make up an object, non-uniform grid lines are generated first. These grid lines are constrained by the structure of the object and follow three rules mentioned in the paper. The first rule is to find demarcation points of the object and drop grid lines on these points. The second rule is to make sure all mesh sizes are less than one-tenth of the wavelength by adding more grid lines. The last rule is to densify mesh at the fine structure of the object. Then by comparing the positional relationship between center points of Yee cells and tetrahedrons, the object can be discretized by Yee cells. Finally, numerical examples are given to verify the validity and accuracy of this novel method.

Index Terms: Adaptive mesh generation, FDTD, tetrahedron.

I. INTRODUCTION

Maxwell’s equations reveal the universal law of electromagnetic phenomena in nature. All electromagnetic problems can be attributed to the solution of Maxwell’s equations when given boundary conditions. The finite-difference time-domain (FDTD) method, as one of the three common electromagnetic field numerical calculation methods, starts from Maxwell’s curl equations in differential form and uses Yee cells to get results [1]. Electric field and magnetic field are distributed alternately in space and time, and they are solved gradually in time. The FDTD algorithm is widely used in time domain analysis of electromagnetic scattering and antenna design because it is easy to comprehend and implement [2].

When using the FDTD method to study the electromagnetic characteristics of a target, the discrete object model consisting of Yee cells is needed. For a simple model, Yee cells can be obtained by dividing the model manually. However, for a complex model, it is very complicated to generate Yee cells artificially and easy to produce errors. Therefore, using the theory of computer graphics and computer-aided computation to generate Yee cells is necessary [3].

Currently, there are not so many types of commercial software on the market that generate Yee cells. Some famous graphics software such as AutoCAD, SolidWorks, or Gmsh are used to construct three-dimensional (3D) models of targets and to provide triangular or tetrahedral meshes in different file formats [4]. Thus, some scholars have already proposed a part of theory and approaches to convert other mesh formats into Yee cells. In 1993, W. Sun proposed a ray tracing method of 3D FDTD surface mesh generation based on triangular mesh [5]. In 1996, W. Heinrich proposed an optimal FDTD mesh generation method, which laid the foundation for the emergence of adaptive mesh, but this method could only solve the simple shape model with relatively regular shape [6].

In recent years, the number of generation methods for Yee cells based on tetrahedral mesh has also increased. Y. Gong proposed a scheme using the theory of convex geometry to judge the position relationship between Yee cell and each tetrahedral element of the object [7]. However, each judgment needs to compare eight vertices of Yee cell and discuss six circumstances, which is time consuming and complex. Liang Hui proposed an improved method which only needs to judge the position relationship between the center point of Yee cells and the tetrahedron [8]. The target model was reconstructed by Yee cells effectively without multiple calculation, but only a sphere example was verified. Besides, this method was only suitable for convex bodies and cannot be used to handle more complicated structures. Z. Yu also adopted an algorithm to judge the positional relationship between points and tetrahedrons, and verified an aircraft model by Yee cells [9]. Nonetheless, only uniform mesh is used in the paper, while how to divide the object with non-uniform mesh is not mentioned. But non-uniform mesh is also needed in practice because it can save total calculation time while ensuring the accuracy of subdivision.

Therefore, a novel and effective method which can convert the tetrahedral data into Yee cells is proposed in this paper. According to the tetrahedral mesh data, we find the coordinates of demarcation points and drop grid lines on them. In addition, we guarantee all mesh sizes are less than one-tenth of the wavelength, and we densify mesh only at the precise structure of the object. Finally, the non-uniform Yee cells meshes are generated. Three examples are given to prove the accuracy of the subdivision method.

II. ADAPTIVE MESH GENERATION

There are two steps for the adaptive mesh generation based on tetrahedral data. The first step is to divide the target object with non-uniform grid lines. The second step is to discretize the target object by Yee cells.

Figure 1: The optimal bounding box.

A. Generate non-uniform grid lines

In this section, the method of how to place grid lines reasonably is presented. At present, many sorts of computer aided design (CAD) software can provide tetrahedral meshes, which can be used as front-end data for non-uniform grid lines placement. Gmsh is a commercial finite element mesh generator, which supports reading of various file formats and can generate tetrahedral meshes automatically. The tetrahedral data exported from Gmsh includes the coordinates of all vertex and vertex index numbers of each tetrahedron.

First, the 3D boundary of the target should be found to determine the setting range of non-uniform grid lines. According to the vertex coordinates, we can easily find the maximum and the minimum values of the object along x-, y-, z-axes. Depending on these values, an optimal bounding box can be constructed as shown in Fig. 1. All grid lines should be limited in this box to reduce the dividing time and internal storage.

The next step is to place adaptive grid lines in the box. There are three rules for the subdivision.

The first rule is to make sure grid lines fall on demarcation points. We find a novel way to seek out the demarcation points depending on the data of tetrahedrons. The coordinates of these points are all the places where objects have right-angle structures. Bringing in these points can make the partition more accurate, and this method is suitable for many different shapes of objects. For an object having slit structure, we can accurately locate the position of the gap, and get the width and length of the gap. For the object similar to the structure of a rectangular patch antenna, the position of the patch, as well as the width and length of the lines, can also be found precisely. Take a model for example. Two cuboid boxes are nested together. The outer cuboid box has a hollow structure in the middle, while the inner cuboid box is solid. Figure 2 shows the cross section of this model on the x-z plane. The yellow part represents the outer entity box, the blue part represents the inner one, and the orange dots represent the demarcation points.

Figure 2: The cross section to illustrate demarcation points.

It can be easily seen that these orange dots are the key points that distinguish the physical object from the air, where the structural changes happen. If the coordinates of these points can be obtained and the grid lines are guaranteed to fall in these points, the accuracy of subdivision can be improved. Although these points can be clearly seen from Fig. 2, we can only use the front-end data of tetrahedrons without seeing the model, so an algorithm is needed to implement it.

Take the x direction for example. A tetrahedron has four points and four surfaces. Three points in different positions can determine a plane. If three points of a tetrahedron have the same coordinates on the x-axis, the surface formed by these three points of this tetrahedron must be perpendicular to the x-axis. At this time, record the coordinate value of this surface in the x direction and store it in an array. As is shown in Fig. 3, record the x coordinate value of surface ABC.

Figure 3: The algorithm to find demarcation points.

After traversing all tetrahedrons and getting all qualified x coordinate values in the array, we sort this array and eliminate the duplicate values. Then we get the final demarcation points in the x direction. The demarcation points in the y and z directions can be obtained in the same way as those in the x direction. After finding the demarcation points and putting grid lines on these points, a sketch map is shown in Fig. 4.

Figure 4: The grid lines after the first rule.

The second rule is to find the upper limit for the size of the maximum mesh. Taking wavelength as the benchmark, it is to find the minimum number of grids in each coordinate direction per wavelength, while the wavelength refers to the highest frequency of simulation. According to the FDTD algorithm, this number means the spatial sampling rate for the input signal, so it has a strong influence on the quality of the results and on the calculation time. Increasing it leads to a higher accuracy, but unfortunately it also increases the total calculation time. A good compromise between calculation time and the achievable accuracy is the default value of ten. Therefore, the size of the maximum mesh should be less than one-tenth of the wavelength.

In order to satisfy the second rule, we have to add new grid lines. Because the intervals of the grid lines obtained based on the first rule may be greater or less than one-tenth of wavelength. For the grids less than one-tenth of wavelength, no processing is required. For the grids larger than one-tenth of wavelength, new grid lines have to be added evenly until all of them are less than one-tenth of wavelength. The process can be described from Figs. 5 (a) and 5 (b).

Figure 5: The grid lines after processing the second rule: (a) Original schematic diagram and (b) grid lines added.

The third rule is to adjust the mesh size depending on the tetrahedrons because the density of tetrahedrons corresponds to the fineness of the object structure. In addition, we believe that the tetrahedral mesh obtained by Gmsh is reasonable and accurate. It is wise to use small grids where the tetrahedrons are dense and large grids where the tetrahedrons are sparse. Therefore, we come up with a new algorithm.

Take the x direction for example. As shown in Fig. 6, if a tetrahedron ABCD falls completely between two adjacent grid lines in the x direction, it indicates that this grid is imprecise. A new grid line should be added in the middle of this grid which splits it into two new grids, and both are half of the original grid.

Figure 6: The algorithm description of the third rule.

Figure 7: The process for the third rule: (a) Original schematic diagram and (b) grid line added.

Continue judging whether the tetrahedron still falls completely in the new grid. If so, do that step again, and this time the grid size becomes a quarter of the original one. Until the tetrahedron does not fall completely into only one grid, the work is done.

The previous model is still used to illustrate this process. In Fig. 7 (a), supposing a tetrahedron ABCD (only three points A, B, C are shown in the cross section) falls completely in a grid, it means that this grid is imprecise. A new grid line should be added as shown in Fig. 7 (b).

In addition, we also consider an extreme case so as to make the algorithm more general and robust. If a tetrahedron falls in a grid and the size of the tetrahedron is too small compared with the size of the grid, the work of adding new grid lines will be repeated, resulting in excessive density. For example, in a microstrip antenna, the thinnest layer of material needs to be divided into two grids along the thickness direction, while a smaller mesh is unnecessary. Therefore, we add a new step that allows operators themselves to input the expected minimum grid size into the program to avoid too small a mesh. When the mesh size calculated by the algorithm is smaller than the expected value, the program will retain the previous grid lines without further processing.

To summarize, the adaptive mesh generation of the object needs to generate grid lines in the bounding box, and these lines should follow the three rules as mentioned above. The whole process can be described as: Find the demarcation points of the object according to the first rule, and drop the grid lines on these points. Then add grid lines evenly to ensure that all grids are less than one-tenth of the wavelength, satisfying the second rule. Finally, the sizes of tetrahedral meshes are compared with the sizes of grids, and more grid lines are added only at the fine structure of the object. All grid lines for the subdivision are shown in Fig. 8.

Figure 8: The final subdivision result.

B. Use Yee cells to reconstruct the target object

From part A, grid lines are generated. Supposing the shadow part represents the entity object in Fig. 9, the key for the object discretization by Yee cells is to determine which Yee cells are located in the object. That is, to compare the position of every center point of Yee cells with the actual location of the object.

Figure 9: The grid lines and Yee cells.

Before we import the algorithm of determining whether a point is inside an object in three dimensions, we first introduce a special situation judging whether a point is on a plane in two dimensions in order to make it easier to understand. And in this situation, it also has a practical use because in some structures the metal layer of the object can be seen as non-thickness, which can reduce the amount of mesh. We just need to mark the index numbers of grids where the metal is and set the tangential electric fields to zero on these locations.

According to the front-end data, the surface of the object is formed by several triangles. The key to judge the positional relationship between a point and a plane is to judge it between a point and a triangle. If a point is in one of the triangles forming the plane, the point must be in that plane.

Suppose the three vertices of a triangle are A, B, C, and the vertices are in anticlockwise direction. P is the point to be tested. According to the mathematical formula, if the results of three cross-product expressions have the same sign, the point P is in the triangle ABC.

{\begin{matrix} C P \times C A \\ A P \times A B \\ B P \times B C \end{matrix}} .

(1)

After we have finished the case in two dimensions, now let’s look at the circumstance in three dimensions. Based on the front-end data, the object is formed by several tetrahedrons. The key to judge the positional relationship between a point and an object is to judge it between a point and a tetrahedron. We use directed volume to finish the judgment in three dimensions, and this is an analogy to the two dimensions.

Suppose there are four vertices of a tetrahedron A, B, C, D, and the vertices are in anticlockwise direction. P is the point to be tested.

Assuming the coordinates of the five points A, B, C, D, P are, respectively,

A = (x 1, y 1, z 1),

B = (x 2, y 2, z 2),

C = (x 3, y 3, z 3),

D = (x 4, y 4, z 4),

P = (x, y, z),

and according to the cross products algorithm, the formula for calculating the directed volume of tetrahedron can be regarded as the calculation of determinant. For example, the formula for calculating the volume of tetrahedron ABCD is

V_{a b c d} = | \begin{matrix} x_{1} & y_{1} & z_{1} & 1 \\ x_{2} & y_{2} & z_{2} & 1 \\ x_{3} & y_{3} & z_{3} & 1 \\ x_{4} & y_{4} & z_{4} & 1 \end{matrix} | .

(2)

We can use the coordinates of point P to replace points A, B, C, and D, respectively and get the results of the other directed volume of tetrahedrons BCDP, ACDP, ABDP, ABCP in sequence. Judge and compare whether these five values have the same sign. If they have the same sign, it is proved that the point P is located in the tetrahedron ABCD, and vice versa [10].

Using the algorithm mentioned above, we can mark the Yee cells that are located inside the tetrahedron, and these Yee cells must also be located inside the object. By traversing all tetrahedrons, we can find all Yee cells which are inside the object, and according to the positions of them we can restore the object.

III. EXAMPLES AND RESULTS

In this section, we illustrate the efficiency and the robustness of the adaptive mesh generation method for FDTD simulation. The first model is an electrically large aircraft model, which is chosen to verify the stability of the method in processing a large number of meshes. The second model is a rectangular microstrip patch array, which is used to show the accuracy of the non-uniform grid generation technique. A schematic of a human skull serves as the third model, demonstrating the ability of our algorithm to handle complicated models. In all examples, we use CPML as absorbing conditions. For the meshes in absorbing boundary, we use uniform meshes, and the sizes are one-tenth of the wavelength in air. The meshes between object and absorbing boundary are also one-tenth of the wavelength in air. These meshes are all added directly in our FDTD computation process. In order to present partition results more clearly, we don’t add these meshes in subdivision pictures. The total number of meshes is only calculated in the object area.

A. The aircraft model

Figure 10: Aircraft model and subdivision results: (a) Original model, (b) subdivision results and (c) a part of discretized rear side of the aircraft.

We use an aircraft model to verify the performance of the program based on the meshing algorithm mentioned before. The fuselage of the plane is about 20 m. The height is 3.78 m, and the wingspan is 13.20 m. The simulation frequency range is from 0.1 GHz to 0.5 GHz. According to the second rule in the last part, the maximum size of the mesh should be less than one-tenth of the wavelength, we finally choose it as 0.05 m. The minimum size of the mesh is 0.015 m. For all meshes, 0.05 m refers to the largest side length, while 0.015 m refers to the smallest side length. The smallest side and the largest side, however, may not exist in one single mesh. On the basis of adaptive mesh generation algorithm, the number of grids in the modeling area is 643 $\times$ 393 $\times$ 79 = 19,963,221. The original model is shown in Fig. 10 (a), while the model discretized by Yee cells is shown in Fig. 10 (b). An enlarged image of the discretized rear side of the aircraft is shown in Fig. 10 (c).

The radar cross section (RCS) of this model is calculated by using the electromagnetic simulation software CST. We use the FIT engine in CST with hexahedral meshes. Meanwhile, the FDTD method is also used to calculate the RCS based on the subdivision mesh in Fig. 10. We also employ the message passing interface (MPI) method to save time. This method of parallel processing may distribute the entire computation work among several processor nodes, increasing calculation speed. After setting the incident plane wave from the belly of the aircraft, the RCS results over frequency are shown in Fig. 11. The RCS results on phi = 0 and phi = 90 plane at 250 MHz are shown in Fig. 12.

Figure 11: Radar cross section results over frequency: (a) Forward RCS and (b) backward RCS.

We can see from Figs. 11 and 12 that the simulation results obtained by CST software and the FDTD method are consistent with each other, which verifies the stability of the adaptive mesh generation method proposed in this paper.

Figure 12: Radar cross section results at 250 MHz: (a) phi = 0 and (b) phi = 90.

B. 2 $\times$ 2 rectangular microstrip patch antenna array

The second model is a patch antenna array. The antenna contains three parts: patch, dielectric substrate, and ground floor. The relative permittivity of the substrate is 2.2, and the thickness is 6 mm. The structure of the antenna is shown in Fig. 13, and the values of the parameters are given in Table 1.

Figure 13: Top view of the patch antenna array.

Table 1: Parameters of the patch antenna array

Parameter	Length(mm)	Parameter	Length(mm)
Metal thickness	0.12	Lp	36.8
Wm	8	Wp	48
W1	1.7	Sy	101
Lm	23	Sx	88
L1	45.3

Figure 14: Patch antenna array model and subdivision results: (a) Original model, (b) tetrahedral meshes, (c) subdivision results and (d) feed lines partially enlarged.

The simulation frequency range is from 2.2 GHz to 2.7 GHz, and the operating frequency is 2.4 GHz. The minimum mesh size is 0.12 mm. According to the adaptive mesh generation algorithm, the number of grids in modeling area is 44 $\times$ 62 $\times$ 5 = 13,640. The original model is shown in Fig. 14 (a). The tetrahedral meshes for the patch and substrate are separately shown in Fig. 14 (b). The model discretized by Yee cells is shown in Fig. 14 (c), and the thin feed lines partial of the model is enlarged in Fig. 14 (d).

Figure 15: S $_{11}$ parameter over frequency.

Figure 16: Radiation patterns at 2.4 GHz: (a) phi = 0 and (b) phi = 90.

We also use two ways to obtain results. One is by using the CST software, the other is by using FDTD method based on the subdivision cells in Fig. 14. The comparison results of S $_{11}$ and radiation patterns are shown in Figs. 15 and 16.

It can be seen from Figs. 15 and 16 that two groups of results are consistent with each other. Because the feed lines in the network are very narrow, a little change in the width of lines can lead to the electromagnetic performance variation of the antenna. Therefore, the high precision of the mesh generation algorithm proposed in this paper can be verified through this example.

C. The human skull model

We use this example to illustrate the performance of our method in dealing with complicated models. The whole model consists of many organs and skin tissues: skull, trachea, brain, eyes, tongue, teeth, jaw, etc. [11]. The simulation frequency range is from 1 GHz to 5 GHz. The maximum size of the mesh should be less than one-tenth of the wavelength; we finally chose 2 mm. The minimum size of the mesh is 0.5 mm. The total number of grids in the modeling area is approximately 7,400,000. The original model is shown in Fig. 17, while the model discretized by Yee cells is shown in Fig. 18.

Figure 17: Human skull model: (a) Whole model and (b) inside model without skin and fat.

Figure 18: Subdivision results: (a) Whole model subdivision and (b) inside model subdivision.

Using our approach, the discretized skull model of Yee cells is in good agreement with the original one, as shown in Figs. 17 and 18. The findings can demonstrate that our technique is capable of handling intricate models, since the skull model is made up of a variety of organ forms with different shapes.

IV. CONCLUSION

In this paper, a novel method for adaptive FDTD mesh generation based on front-end tetrahedral data is presented. The first step is to divide a target object with grid lines. These grid lines are constrained in the bounding box and are placed in a non-uniform way based on the structure of the object. The second step is to compare the positional relationship between center points of Yee cells and the object, which can realize the object discretization by Yee cells. Relevant principles are explained, and two simulation examples are verified. The first example shows the stability and reliability of this method in processing a large-scale example. The second one shows the accuracy of the non-uniform mesh technique. The correctness of the research content is proved.

ACKNOWLEDGMENT

This work was supported by National Natural Science Foundations of China (No. 61971340 and 62122061) and by National Key Research and Development Program of China (No. 2020YFA0709800).

REFERENCES

[1] K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equation in isotropic media,” IEEE Trans, vol. 14, 1966.

[2] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd edition, 2005.

[3] Y. Srisukh, J. Nehrbass, F. L. Teixeira, J.-F. Lee, and R. Lee, “An approach for automatic grid generation in three-dimensional FDTD simulations of complex geometries,” IEEE Antennas and Propagation Magazine, vol. 44, no. 4, pp. 75-80, Aug. 2002, doi: 10.1109/MAP.2002.1043151.

[4] H. Zhu, C. Gao, H.-L. Chen, Z.-H. Shi, and X.-H. Xu, “The research on FDTD mesh generation and visualization technology,” 2012 6th Asia-Pacific Conference on Environmental Electromagnetics (CEEM), Shanghai, China, pp. 282-284, 2012, doi: 10.1109/CEEM.2012.6410622.

[5] W. Sun, C. A. Balanis, M. P. Purchine, and G. Barber, “Three-dimensional automatic FDTD mesh generation on a PC,” Proceedings of IEEE Antennas and Propagation Society International Symposium, vol. 1, pp. 30-33, 1993, doi: 10.1109/APS.1993.385409.

[6] W. Heinrich, K. Beilenhoff, P. Mezzanotte, and L. Roselli, “Optimum mesh grading for finite-difference method,” IEEE Transactions on Microwave Theory and Techniques, vol. 44, no. 9, pp. 1569-1574, Sep. 1996, doi: 10.1109/22.536606.

[7] Y. Gong, Z. Wu, and Z. Dai, “Convex object Yee cell model building based on convex geometries,” Journal of Tsinghua University, vol. 47, no. 9, pp. 1521-1525, 2007.

[8] L. Hui and S. Zuxun, “FDTD mesh-generating and visual realization based on triangle-patch,” Computer Simulation, vol. 26, no. 11, pp. 106-109 (in Chinese), 2009.

[9] Z. Yu, W. Zhao, L. Xu, and X. Shi, “A novel mesh generation method for FDTD without ray-tracing,” 2017 Sixth Asia-Pacific Conference on Antennas and Propagation (APCAP), Xi’an, China, pp. 1-3, 2017, doi: 10.1109/APCAP.2017.8420508.

[10] J. W. Boardman, “Analysis, understanding, and visualization of hyperspectral data as convex sets in n space,” Proceedings of SPIE - The International Society for Optical Engineering, vol. 2480, pp. 14-22, 1995.

[11] S. N. Makarov, G. M. Noetscher, J. Yanamadala, M. W. Piazza, and S. Louie, “Virtual human models for electromagnetic studies and their applications,” IEEE Reviews in Biomedical Engineering, vol. 10, pp. 95-121 [Online], 2017. Available: https://www.nevaelectromagnetics.com/vhp-female-2-2.

BIOGRAPHIES

Weiran Zhang was born in Shanxi, China. He received the B.S. degree from Xi’an Jiaotong University, Xi’an, China, in 2017, in information engineering. He is currently working in Xi’an Jiaotong University as a postgraduate. His research interests include mesh generation and computational electromagnetics.

Zikun Xu was born in Jiujiang, China. He received the B.S. and M.S. degrees from Xi’an Jiaotong University, Xi’an, China, in 2023, all in electromagnetic field and microwave technology. His research interests are parallel FDTD methods based on MPI.

Huaiyun Peng was born in Hubei, China. He received the Ph.D. degree from Xidian University, Xi’an, China, in 2017, in radio physics.

He is currently working in National Key Laboratory of Eletromagnetic Environment, China Research Institute of Radiowave Propagation, Qingdao, China, as a researcher. Her research interests include radio wave propagation, numerical method for electromagnetic field, and electromagnetic scattering.

Juan Chen was born in Chongqing, China. She received the Ph.D. degree from Xi’an Jiaotong University, Xi’an, China, in 2008, in electromagnetic field and microwave technology. She is currently working in Xi’an Jiaotong University, Xi’an, China, as a professor. Her research interests include the computational electromagnetics, microwave device design, etc.

Chunhui Mou was born in Yantai, China. She received the B.S. and M.S. degrees from Xidian University, Xi’an, China, in 2012 and 2015, and the Ph.D. degree from Xi’an Jiaotong University, Xi’an, China, in 2023, all in electromagnetic field and microwave technology. She is currently working in Xi’an Jiaotong University, Xi’an, China, as a postdoctoral researcher. Her research interests include the fast FDTD method, FDTD mesh generation method, and multi-physical field calculation.