## 1 Introduction  
Computer-aided design (CAD) models have played crucial roles in technical fields, such as numerical simulation, additive manufacturing and visualization. Accordingly, the triangular mesh, generated from CAD models, has become the dominant type of functional element due to its advantages of efficiency, simplicity and flexibility. However, triangular meshes, which are always transformed directly from a CAD model, usually have the defects of excessively small or large angles, and redundant or irregular vertices. These defects tend to reduce the accuracy of the numerical simulation based on triangular meshes. Thus, it is necessary to establish effective remeshing algorithms to optimize the generated triangular meshes, which is important for the downstream applications. [PERSON] et al. [1] summarized many theories and applications of remeshing algorithms, which can be divided into the quantity-dominant simplification algorithms and quality-dominant optimization algorithms. In the aspects of simplification algorithms, [PERSON] et al. [2] presented a serious of basic mesh operations for mesh simplification. [PERSON] and [PERSON] [3] also proposed a kind of quadratic error metrics (QEM) simplification method, which is also regarded as the general model simplification algorithm. As for the optimization algorithms, [PERSON] and [PERSON] [4] first introduced the Laplace-Beltrami (LB) operator to improve the quality of the irregular triangular meshes. Based on the LB operator, they were able to construct area-equalizing triangulations with vertex smoothing operation and optimize the original mesh, which was regarded as the basic optimization algorithms. Furthermore, [PERSON] et al. [5] transformed the original judgment condition for edges into the angles, enabling the simplification after the mesh optimization. [PERSON] et al. [6] have realized the constraint of the shape of the triangle by varying the size of bubble based on the dynamics simulation method. [PERSON] et al. [7] also tried to plan a more reasonable distribution of the number of vertex-to-vertex connections and optimize the meshes by constraining the range of angles. Most of these algorithms were based on the method of area-equalizing triangulation, which were used to achieve optimization and simplification of meshes in forms of regular triangles. Although the regular triangles can improve the computational accuracy of simulation results, the overemphasis on the evenly constructed triangular elements may limit the reduction of data storage and cause destruction to the geometric feature information of the original model. It is considered that the accuracy of numerical simulations can be improved when the cotangent formula \(\cot\alpha+\cot\beta\) of the discrete LB operators has non-negative weight. According to this theory, a geodesic distance based on the discrete LB operator was used to simulate the thermal diffusion by [PERSON] et al. [8]. Moreover, [PERSON] and [PERSON] [9] proved that the Delaunay triangulation could always satisfy the non-negative weight for discrete LB operators, which was defined as the local property of the Delaunay triangulation. Correspondingly, the global property of the Delaunay triangulation was defined as the maximization of all the minimum angles in the triangulation, which can be used to indirectly constrain the shape of the triangle and avoid narrow triangles [10]. It is obvious that the Delaunay triangulation has the advantages in protecting the geometric features and simplifying the original meshes. However, it is difficult to achieve the Delaunay triangulation on the Riemannian manifold, considering that the geodesic paths are difficult to be directly represented with existing data structures. To solve this problem, [PERSON] et al. [11] defined a special concept of Delaunay mesh that, as long as the sampling criterion acquired a relatively lower bound among the distances of vertices in the Voronoi diagram, the Delaunay mesh can be obtained. Both [PERSON] et al. [12] and [PERSON] et al.