Due to the mathematical guarantees of mesh quality, there have been extensive efforts to develop Delaunay-based algorithms that can efficiently handle arbitrarily complicated geometry for 2D domains [1]}, [2]}. However, extending these algorithms to 3D is challenging. In the past three decades, many studies have been conducted regarding the construction of tetrahedral partitions given complex domains [3]}, [4]}, [5]}, [6]}. Different in the ways of handling domain boundaries and improving mesh quality, the tetrahedralization methodologies can be roughly characterized into three groups: boundary-constrained, Delaunay refinement, and variational methods; see [7]} for more detailed descriptions.
