Keywords: intrinsic mesh, discrete isogeometric triangulation, geodesic, shortest path, A* search, coarsening, direct high-order meshing

1 Introduction↩︎

The numerical representation of surfaces plays a crucial role across a wide range of domains, including geometric modeling and design, rendering, surface analysis, geometric processing, and computational physics. Broadly speaking, two main categories of surface representations can be distinguished: continuous representations and discrete representations.

The first category consists of continuous surfaces, typically originating from design blueprints or analytical descriptions. In computer-aided design (CAD), surfaces are most commonly defined by explicit parametric mappings \[(x,y,z) = \mathbf{x}(u,v),\] where \((u,v)\) belong to a parametric domain, usually rectangular. In practice, these parametric patches are often trimmed, meaning that only a subset of the parameter domain is retained. The trimming curves are themselves typically defined as parametric curves in the \((u,v)\) space. As a result, CAD models are generally composed of collections of trimmed parametric patches that must be assembled to form a complete surface.

Other types of continuous representations also exist. For instance, implicit surfaces are defined as the zero level set of a scalar function \[f(x,y,z) = 0,\] as commonly encountered in level-set methods, signed distance functions, or algebraic surface descriptions.

The second major category corresponds to discrete surface representations. These arise naturally from numerical approximation, geometric processing, or measurement procedures such as scanning and reconstruction. Typical examples include polygonal meshes (most commonly triangle meshes), triangle soups, and point clouds. Such representations are sometimes referred to as organic surfaces, emphasizing that they do not rely on an underlying analytical description but rather on sampled geometric data.

From the perspective of computational physics, our primary objective is to construct meshes suitable for numerical simulation starting from a wide range of surface representations, including CAD models, implicit descriptions, or discrete datasets. These meshes should satisfy two essential requirements: first, they must exhibit high quality to ensure the accuracy, robustness, and efficiency of the underlying numerical methods; second, they must remain as faithful as possible to the original geometry, whether it is defined in a continuous or discrete manner.

In this work, our objective is to construct isogeometric meshes, i.e. meshes that represent exactly the underlying geometry. When the surface is given by a parametric description, this objective can in principle be achieved in a straightforward manner. Consider a surface defined by an injective mapping \[\mathbf{x}(u,v) : \Omega \subset \mathbb{R}^2 \rightarrow \mathbb{R}^3 ,\] where \(\Omega\) is, for instance, a rectangular domain in the parametric space. A non-overlapping triangulation of \(\Omega\) can then be mapped through \(\mathbf{x}(u,v)\) to obtain a triangulation of the surface that represents the geometry exactly. In other words, an isogeometric triangulation can be obtained by triangulating the parametric domain and transporting this triangulation to physical space.

In practice, however, most mesh generation procedures rely on linear triangles constructed directly in physical space, where only the mesh vertices lie on the surface. The resulting meshes therefore approximate the geometry rather than representing it exactly. Nevertheless, when a parametric representation is available, the construction of isogeometric triangulations remains conceptually simple.

The situation is significantly more challenging when the surface is defined in a discrete manner. In this case, no global parametric mapping is available, and the notion of exact geometric representation must be reconsidered. The purpose of this paper is precisely to construct isogeometric triangulations for discrete surfaces, following an approach that is conceptually similar to standard finite element mesh generation, but adapted to preserve the intrinsic geometry of the discrete model.

Starting from a watertight triangulation of a closed surface, we propose to build an intrinsic triangulation of that surface, in which the edges correspond to discrete geodesics and the faces are composed of sets of small triangles inherited from the initial discrete geometry, as illustrated in Fig. 2. In this sense, the new triangulation is defined intrinsically on the surface while remaining fully consistent with the underlying discrete representation.

Our approach does not rely on modifying the original triangulation. Instead, it closely follows the philosophy used for meshing CAD surfaces, where a new mesh is constructed on top of the geometry through local mesh operations. The intrinsic triangulation is progressively obtained using classical local transformations such as edge swaps, edge splits, edge collapses, and Delaunay refinement. These operations are performed intrinsically on the surface, allowing the triangulation to evolve while preserving the underlying discrete geometry and maintaining a watertight representation.

We show in this paper that the proposed approach is both robust and sufficiently fast to be used in a mesh generation context, where users typically expect algorithms capable of handling large datasets with short turnaround times. The construction relies only on local operations and discrete geodesic computations, making it well suited for practical large-scale applications.

As a by-product of this intrinsic construction, we also propose a procedure to derive from the intrinsic mesh a representation more directly usable in the context of finite element analysis, namely a high-order mesh. The key feature of this approach is that the high-order geometry is obtained directly from the intrinsic triangulation, without relying on the classical paradigm in which a piecewise linear mesh is first generated and subsequently curved.

In that sense, and without much exaggeration, the method can be viewed as a first attempt to generate high-order meshes without ever passing through the curvature of a classical piecewise linear mesh. The geometry is instead defined directly from the intrinsic structure of the surface, with high-order elements naturally following discrete geodesic edges and intrinsically defined faces.

2 Related Work↩︎

Simple meshing algorithms have been extensively developed for planar surfaces, often based on Delaunay triangulations [1]–[4]. These methods aim to represent surfaces using the fewest possible elements while remaining as close as possible to the exact geometry and satisfying element quality and size criteria, which may vary spatially. Such planar methods can be naturally extended to curved surfaces by computing distances and edges in three-dimensional Euclidean space, as demonstrated in [2]. In this context, element size criteria can be defined to ensure that the mesh accurately captures the surface’s underlying features, including curvature and local geometric details [5].

However, although these methods can be used for meshing, they typically rely on a parametrization of the original surface, which must be available and may introduce artifacts. In [6], [7], remeshing methods are proposed that rely exclusively on local mesh operations without a parametrization. Nevertheless, inserting new points involves a local approximation of the surface, which can introduce inaccuracies even when the input surface is exactly known. Regardless of the method, the resulting mesh represents curved geometries using straight-sided elements, meaning that accurately capturing such geometries often requires a large number of elements.

2.1 Intrinsic Meshing↩︎

To improve surface remeshing methods, geodesic distances computed directly on the surface can be incorporated. This allows meshing operations to account for the actual geometric distances between points. Several methods [8]–[13] use intrinsic distances to perform more accurate operations such as vertex insertion/repositioning or surface parameterization. However, in most of these approaches, the mesh is still represented with straight-sided edges, which limits the advantages of using geodesic information. In contrast, other methods have been proposed to exploit geodesic distances more extensively. In these approaches, the connection between two points is no longer a straight line in 3D space but a piecewise linear curve computed directly on the surface. As a result, the surface is no longer approximated but is exactly represented by the generated elements.

The first intrinsic methods can be viewed as direct extensions of Delaunay triangulations to piecewise linear surfaces, defining intrinsic Delaunay triangulations [14], [15] and refinement schemes [16]. In these cases, intrinsic triangulations and all their properties are straightforward to compute because the intrinsic triangles can always be flattened onto the plane. As a result, angles and edge lengths can be computed using standard planar formulas. Moreover, the use of such developable triangles in computation is simple, since they can be treated as standard straight-sided elements without any additional effort. However, the developable nature of these elements imposes significant limitations, as intrinsic triangles alone cannot represent regions with nonzero Gaussian curvature. A direct consequence is that these methods do not allow vertex removal in the mesh and, therefore, can only be used for refinement, not for coarsening. This is particularly unfortunate, as one of the primary goals of meshing is often to reduce the number of elements compared to the input discrete surface. It should be noted that intrinsic coarsening was addressed in [17]. However, rather than explicitly computing geodesics, the authors locally modify the geometry by flattening vertices, thereby reducing the problem to a planar configuration in regions where vertex removal is required. This approach results in a mesh in which edges are no longer geodesic.

By adopting a more general framework for intrinsic triangles, it becomes possible to work with non-developable intrinsic triangles. For instance, [18] constructs a geodesic Delaunay triangulation on sampled surface points to achieve isotropic mesh simplification. Nevertheless, it should be noted that intrinsic Delaunay triangulations do not necessarily exist for an arbitrary set of sampled points on a surface [13], [19]. In this work, we present an approach for constructing an intrinsic triangulation of a surface that enforces constraints on both element size and quality. This defines a complete intrinsic meshing framework capable of performing both refinement and coarsening. An illustration of these capabilities is provided in Fig. 3, where the different elements are distinguished by various colors. This mesher is based on local optimizations leading to a parametrization-independent and robust method. We also demonstrate how the resulting piecewise-linear elements can be used not only for straight-sided mesh simplification but also for high-order mesh generation.

2.2 Geodesic Distances↩︎

A central requirement of the proposed intrinsic meshing approach is the ability to compute geometric quantities directly on the surface. Local operations require geodesic distances between pairs of points and intrinsic angles associated with triples of points. These angles correspond to the angle on the surface between two geodesics meeting at a point For Delaunay-type refinement, another key operation is the computation of geodesic circumcenters. Given three vertices of an intrinsic triangle, a geodesic circumcenter is a point on the surface that is equidistant, in geodesic distance, from the three vertices. It therefore plays the role of the classical Euclidean circumcenter, but in the intrinsic geometry of the surface. For developable intrinsic triangles [16], lengths and angles can be directly inferred from previously computed information, which significantly simplifies the computation of these geometric features. In contrast, since our focus here is on non-developable intrinsic triangles, the approach requires a mechanism to evaluate geodesic distances on the surface. Once these distances are obtained, other geometric constructs, such as geodesic lines, angles and triangle circumcircles, can be derived.

Numerous approximate methods for computing geodesic distances have been developed in the computer graphics community, many of which are based on variations of the Fast Marching Method [20], [21]. For our purposes, particularly when dealing with complex surface geometries, it is crucial to reduce ambiguities and prevent geodesics from intersecting or becoming indistinguishable. Moreover, each geodesic must be unique. To achieve this, we employ an exact geodesic distance computation rather than an approximate one. A comprehensive survey of both exact and approximate techniques for computing geodesic distances and paths on surfaces can be found in [22].

2.2.1 Exact geodesic distance computation↩︎

The problem of computing shortest paths from a single source to all destinations on a polyhedral surface was first studied in [23], where an algorithm with a prohibitive asymptotic time complexity of \(O(n^5)\) was proposed, with \(n\) denoting the number of edges of the surface. A significant advancement was later proposed in [24] with the Mitchell–Mount–Papadimitriou (MMP) algorithm, which attains a worst-case time complexity of \(O(n^2 \log n)\) and can be interpreted as a continuous analogue of Dijkstra’s shortest path algorithm [25]. An alternative approach, the Chen–Han (CH) algorithm, was proposed in [26], achieving \(O(n^2)\) time complexity. However, empirical comparisons in [27] show that the MMP algorithm significantly outperforms the CH algorithm in practice, both in terms of computational time and memory usage. In the same study, an improved version of the CH algorithm, referred to as Improved Chen-Han (ICH), is introduced. This algorithm sacrifices the theoretical time complexity advantage of CH by incorporating a priority queue propagation strategy similar to that of MMP, resulting in time complexity of \(O(n^2 \log n)\). The use of a priority queue in ICH allows for effective filtering of unnecessary paths, ultimately producing an algorithm similar to MMP in both time and space complexity, although ICH shows slightly improved performance in some cases [27]. Furthermore, [28] shows that all these algorithms can achieve improvements in both time and space by computing distances using an edge-based structure rather than the classical half-edge structure. Additionally, CH and ICH can be optimized to reduce memory usage by restricting the computation of shortest paths to mesh vertices only. However, this optimization is not applicable here, since shortest paths to points on edges or faces are also required. Consequently, within the scope of this work, the ICH and MMP algorithms exhibit no significant differences.

2.2.2 Optimization based on the application↩︎

Significant research has been conducted to improve the performance of algorithms for the “single-source, all-destinations” problem, notably through techniques such as parallelization strategies [29] and grouped information propagation [30], [31]. While these approaches can accelerate computation, more effective optimizations are possible when the objective is to compute the distance to a specific point on the surface, rather than to all points, as considered in this work. In [11], this is achieved by restricting propagation to an elliptical region around the source and destination, using an approximation of the geodesic distance computed via fast marching methods. This approach was further refined in [32], where a Dijkstra search limited to edges was combined with their own approximate method employing merged windows. The authors also drew inspiration from the A* search algorithm [33], even using the term “continuous A* search.” Specifically, they used Euclidean distances as a lower bound on the distance from a point to the destination. Importantly, these approximations are still only used to filter which windows are propagated.

In this work, the Euclidean distance will be used as a lower bound on the remaining distance from an intermediate point to the destinations. However, its primary role will be to determine the order in which windows are propagated, aligning the approach more closely with the original A* algorithm. This also leads to a considerably simpler method. A comparable idea is explored in [34] with the ICH algorithm. This method will also be extended to compute a point that is equidistant from three given points. This will enable the computation of circumcenters, as in [9], but without requiring recursive mesh refinement. Since a public implementation of the MMP algorithm is available² from one of the authors of [32], and because performance differences between the MMP and ICH algorithms are minor, the experimental results will be based on this implementation of MMP. Adapting the approach to the ICH algorithm would require only slight modifications.

3 Intrinsic Meshing↩︎

During the generation of a mesh, the distribution of vertices is governed by a subtle balance: geometry defines where vertices can be placed, the characteristic length determines their spacing, and quality constraints adjust their relative placement. Consequently, resolving all vertex positions in a single step is fundamentally challenging. As a result, most algorithms adopt an iterative approach to progressively modify the mesh toward the desired refinement and quality. In conventional unstructured meshing methods, the fundamental operations are point insertion, point removal and modification of the connectivity of elements. However, in the context considered here, these operations are constrained by the feasibility and computational cost of evaluating the required geometric elements for each operation. As a result, the local optimization operations considered in this work are limited to edge swapping, edge collapse, edge split, and, finally, triangle split.

3.1 Intrinsic Edge Swap↩︎

Given a set of points sampled on a surface, there exist multiple valid triangulations that can be constructed. Here, “valid” refers to triangulations composed of non-overlapping, non-intersecting triangles that preserve the topology of the underlying surface. However, not all triangulations exhibit the same properties in terms of field interpolation and numerical computation. In particular, nearly degenerate triangles tend to extrapolate values far from their vertices, leading to increased approximation errors and reduced numerical stability.

3.1.1 Delaunay Triangulation↩︎

In dimension \(d\), the Delaunay triangulation of a set of points is a triangulation of their convex hull such that the circumsphere of each \(d\)-simplex contains no other point in its interior. In two dimensions, this condition reduces to the empty circumcircle property for triangles. In low dimensions, Delaunay triangulations can be constructed efficiently using well-established algorithms, which explains their central role in practical meshing methods. The planar Delaunay triangulation is particularly important for mesh generation because it satisfies a max-min angle property: among all triangulations of the same point set, it maximizes the smallest angle. This tends to avoid skinny triangles, which is beneficial for numerical interpolation and simulation. In two dimensions, another important property is that the Delaunay triangulation can be reached from any triangulation of the same point set by a finite sequence of local edge swaps. More precisely, if \(n\) denotes the number of vertices, at most \(O(n^2)\) edge swaps are required in the worst case.

3.1.2 Intrinsic Delaunay Triangulation↩︎

An intrinsic Delaunay triangulation is a triangulation whose vertices lie on a surface and whose edges are represented by geodesic paths on that surface, such that each triangle satisfies an empty geodesic circumcircle condition. More precisely, for every intrinsic triangle, there exists a point on the surface that is equidistant, in geodesic distance, from its three vertices, and the corresponding geodesic disk contains no other vertex of the triangulation in its interior. The classical two-dimensional edge-flipping algorithm suggests a natural extension to surfaces. Instead of flipping straight edges in the plane, one can perform intrinsic edge swaps, where the new edge is recomputed as a geodesic path on the surface. Each swap is driven by a local intrinsic Delaunay criterion: if an edge is not locally Delaunay, the algorithm attempts to replace it by the opposite geodesic edge. In this way, the mesh is iteratively modified to enforce local Delaunayness, while remaining entirely defined on the surface. Fig. 4 illustrates the difference between a classical (extrinsic) edge swap and an intrinsic one.

In this work, mesh edges are defined as the shortest geodesic paths between vertices. This choice eliminates the ambiguity arising from the existence of multiple geodesics connecting the same pair of vertices and provides a robust framework for computing these geodesics consistently. However, restricting edges to shortest paths also limits the set of representable geodesics. Consequently, constructing intrinsic meshes based solely on shortest paths imposes constraints on the admissible triangulations, as discussed in Sec. 3.1.3 and 3.1.5.

3.1.3 Edge Swapping Feasibility↩︎

In the planar setting, it is always possible to swap an edge that is not locally Delaunay, thus verifying that an edge should be swapped is sufficient to guarantee that the operation is feasible. On curved surfaces, this property no longer holds. Indeed, the resulting edge after a swap may fail to remain within the cavity formed by the two adjacent triangles, as it is the case in Fig. 6. It is therefore essential to verify that the candidate edge lies entirely within this cavity. Moreover, to preserve the validity of the mesh, edges must remain non-intersecting except at shared vertices. As described in Alg. 5, both conditions can be enforced through an edge-intersection query: first, by confirming that the current and candidate edges do intersect, and second, by verifying that the candidate edge does not intersect the boundary edges of the cavity.

3.1.4 Edge Swapping Heuristic Criterion↩︎

The Delaunay criterion can also be directly generalized to non-planar surfaces by replacing circumcircles with geodesic distance-based circumcircles. For each edge, one could compute the circumcenters of both adjacent triangles and check whether the Delaunay criterion is locally satisfied. However, the geodesic circumcenter of a triangle on a curved surface is neither guaranteed to be unique nor to even exist, see Fig. 7. Moreover, even if it exists, its computation might be heavy and also might not be possible due to numerical precision limitations as it will be discussed in Sec. 4.4. As a consequence, the geodesic Delaunay criterion is not employed as it is. Instead, a heuristic based on readily available information will be adopted.

The proposed heuristic is required to satisfy several key properties. First, it must be consistent with the Delaunay criterion for developable triangles, which corresponds to the planar scenario. Second, it must ensure coherence, preventing newly introduced edges from being immediately swapped back. Third, it should not presuppose that the candidate edge resulting from a swap lies within the cavity, as this condition does not always hold. Finally, the heuristic must depend exclusively on intrinsic quantities simple to compute, specifically the edge lengths and angles defined on the surface.

In the planar case, a simple implementation of the Delaunay criterion consists in evaluating the two angles opposite to the edge under consideration and checking whether their sum exceeds \(\pi\). If this condition is satisfied, the edge is flipped. Equivalently, this procedure can be interpreted as minimizing the sum of the angles opposite to edges, since the total sum of angles in a quadrilateral is equal to \(2\pi\). Motivated by this observation, the heuristic proposed in this work selects the configuration that minimizes the sum of the angles opposite to the edge. This approach does not rely on any assumptions regarding the candidate edges and depends only on intrinsic angles.

3.1.5 Unreachable Delaunay Configurations↩︎

As it was discussed in Sec. 3.1.3, some configurations prohibit an edge flip when the resulting geodesic would lie outside the cavity or intersect other edges, thereby altering the mesh validity. In such cases, reaching a Delaunay configuration may become impossible if the swap cannot be avoided. An illustration of this situation is provided in Fig. 10 (a).

A possible strategy to recover a Delaunay configuration consists in splitting edges that would otherwise be swapped but cannot be flipped due to topological constraints. Although this operation introduces an additional vertex, it locally modifies the mesh connectivity, thereby enabling the construction of an intrinsic mesh in situations where the surface topology would otherwise prevent it. As illustrated in Fig. 10 (b), inserting a vertex at the midpoint of such problematic edges is sufficient to resolve the issue and reach a Delaunay configuration.

3.1.6 Existence and Uniqueness of Intrinsic Delaunay Triangulation↩︎

Beyond the limitations in representing certain Delaunay configurations, there exist situations in which a valid Delaunay configuration may fail to exist due to the topology of the underlying surface. This issue is independent of the representation of intrinsic triangles and is inherently topological, as discussed in [35]. In that work, the existence of a valid Delaunay triangulation is established under a condition on the dual Voronoi diagram, namely the closed ball property: (1) each Voronoi cell is homeomorphic to a planar disk, (2) the intersection of any two Voronoi cells contains at most one Voronoi edge, and (3) the intersection of any three Voronoi cells contains at most one Voronoi vertex.

Moreover, a constructive procedure is proposed to enforce the existence of a Delaunay triangulation by iteratively inserting additional points along Voronoi edges and pseudo-bisectors, i.e., loci of points admitting multiple shortest paths to a Voronoi vertex. Although the point insertion strategy described in Sec. 3.1.5 differs from that of [35], it is reasonable to expect that the existence of a Delaunay triangulation can still be enforced, provided that a sufficient number of points are inserted.

3.1.7 Edge Swapping Algorithm↩︎

In general, the order in which edges are swapped does not affect the final result, as the sequence of swaps converges to a unique Delaunay configuration. However, in the present context, additional points may be inserted to enforce the existence of such a configuration. Consequently, the order in which points are inserted can influence the final triangulation, potentially leading to different configurations.

Since there is no natural or canonical strategy to order the edge swaps and point insertions, we adopt a simple deterministic approach. Specifically, edge swaps are performed according to a predefined ordering (e.g., based on vertex numbering), ensuring reproducible behavior of the algorithm.

The complete edge-swapping procedure is detailed in Alg. 11. The function doWeSwap applies the heuristic Delaunay criterion to the neighboring triangles, while canWeSwap checks for intersections involving the candidate edge, as described in Alg. 5 and 8.

3.2 Intrinsic Edge Split and Collapse↩︎

A fundamental requirement for meshing algorithms is the capability to enforce spatially varying refinement across the geometry, commonly specified via a characteristic length field. The objective of the meshing process is then to produce edges whose dimensionless lengths are close to unity, so that regions requiring finer resolution correspond to smaller prescribed characteristic lengths. This is accomplished by locally adjusting element sizes.

Edge refinement can be achieved by splitting an edge, which replaces two triangles with four, thereby reducing local edge lengths. Conversely, edge coarsening can be performed by collapsing an edge into a single vertex, eliminating two triangles and updating the connectivity of adjacent elements. These operations are inherently local, as they only require the computation of the new vertex and its incident edges, and they naturally generalize to the intrinsic setting, see Fig. 12. Following each split or collapse, the restoration of a Delaunay configuration involves testing edge flips only within the modified neighborhood. These local edge operations are therefore employed to enforce the prescribed element size constraints.

3.2.1 Edge Split and Collapse Criterion↩︎

The decision to split or collapse an edge is based solely on its dimensionless length relative to a prescribed target length. Since this target value cannot, in general, be satisfied exactly, it is necessary to define an admissible range within which element sizes are considered acceptable. In this context, the minimum and maximum allowable characteristic lengths must be chosen to be at most \(1/\sqrt{2}\) and at least \(\sqrt{2}\) times the target length, respectively, in order to prevent immediate reversal of operations (i.e., collapsing an edge directly after splitting it, and vice versa). This condition, however, remains weak: even when it is satisfied, certain configurations may still lead to oscillatory sequences of splits and collapses. Consequently, in practice, wider tolerances are typically required to ensure robust behavior.

To provide greater flexibility, the characteristic length in this work is defined as a range of admissible element sizes, specified by a minimum and a maximum characteristic length. A single target value (e.g., \(l_c\)) may still be prescribed and systematically converted into an interval by introducing fixed margins (e.g., \([0.5 l_c, 2 l_c]\)). Additionally, the formulation allows for relaxed constraints by omitting bounds on the minimum and/or maximum element size. The primary limitation of this approach is that, if the admissible range is too narrow, the method may be unable to satisfy all constraints simultaneously and consequently fail to converge.

3.2.2 Edge Split Feasibility↩︎

When splitting an edge, a first verification consists in ensuring that the two newly created edges lie entirely within the cavity, i.e., that they do not intersect any boundary edges. In addition, it must be verified that these new edges do not intersect each other.

After the split operation, the mesh must be restored to a Delaunay configuration. To this end, adjacent edges are iteratively examined and flipped when necessary, following a procedure similar to Alg. 11, but restricted to the local neighborhood. If, during this process, an edge cannot be flipped and cannot be further subdivided, the operation is aborted. This limitation implies that, in certain configurations, it may not be possible to enforce the prescribed size constraints through edge splitting alone.

3.2.3 Edge Collapse Feasibility↩︎

The constraints on edge collapses are similar to those for edge splits but are generally stricter. Since the purpose of an edge collapse is to increase the local mesh size by removing a vertex, the operation is aborted if the subsequent edge-flipping phase introduces a new vertex, as this would contradict the intended coarsening.

Additionally, a topological check must be performed to ensure the mesh remains valid after the collapse. Specifically, if the cavity in which the edge is to be collapsed is not a topological disk, the collapse would alter the mesh topology and render it invalid as discussed in Sec. 3.1.6. This condition can be efficiently verified using the triangles and boundary edges of the cavity, as described in Alg. 15.

When collapsing an edge, the new vertex can be placed anywhere within the cavity. However, selecting one of the edge’s endpoints is preferable, as it enables reuse of previously computed geodesics. Also, permitting the edge to collapse to its midpoint can facilitate simplification in cases where collapsing to either of its endpoints is not feasible. Subsequently, the selected point is the one that maximizes a quality criterion that will be later defined in Sec. 3.3.1. After the edge collapse, adjacent edges must also be checked for possible flips in order to restore a Delaunay configuration. This operation may affect a region larger than the local cavity, since the midpoint can be inserted and the Delaunay criterion is enforced only heuristically.

3.2.4 Unstable Delaunay Configurations↩︎

When modifying an intrinsic mesh through vertex insertion or removal, the edge-flipping algorithm may fail to converge to a stable configuration. This limitation arises from the heuristic nature of the Delaunay criterion, which can permit repeated edge swaps and lead to infinite cycling, as illustrated in Fig. 16. In practice, such situations can be detected by imposing a limit on the number of swaps, allowing the operation to be safely rejected.

3.2.5 Edge Split and Collapse Algorithms↩︎

Unlike edge swapping, the order of edge splits and collapses is crucial. Arbitrary splitting can lead to oversampling and may affect convergence, whereas uncontrolled collapsing can prevent further operations, limiting the method’s coarsening capabilities. To address this, edge splitting should begin with the longest edges, while edge collapsing should start with the shortest.

These operations are implemented using a queue of edges sorted by length. After a split or collapse, both the affected edge and any edges no longer in the mesh must be removed from this queue, and any newly created edges eligible for further operations must be inserted. The full edge splitting and collapsing procedures are detailed in Alg. 17 and 18.

3.3 Intrinsic Triangle Split↩︎

Using only edge swap, split, and collapse operations, a mesher can achieve its primary goals. Splits and collapses adjust edges toward the target size defined by characteristic lengths, while swaps improve triangulation quality via a Delaunay-like criterion. An example intrinsic mesh is shown in Fig. 19. For some applications, this mesh may be sufficient, but additional control may be required to constrain other mesh properties.

Some triangles remain challenging to approximate in regions with concentrated geometric details. Even when edges meet the target lengths, triangle areas can remain large, as sizing depends solely on edge lengths. A quality criterion is therefore needed to identify such triangles, and a triangle split operation is introduced to refine them.

3.3.1 Triangle Quality↩︎

In conventional meshes, triangle angles are the primary measure of element quality, as extreme angles can significantly degrade computational performance. Other metrics, such as aspect ratio or circumradius-to-shortest-edge ratio, also provide valuable information but are ultimately determined by the angles in the plane. In intrinsic meshes, however, these relationships no longer hold: different quality measures can convey complementary information.

Nonetheless, angle-based quality remains a useful and straightforward indicator. Metrics involving the circumcircle or the area of an intrinsic triangle are generally impractical, as their computation is costly. For simplicity, we therefore define a quality \(q_0\) based on the intrinsic angle \(\theta\), analogous to the planar case: \[q_0(\theta) = \frac{3}{\pi} \min\left(\theta, \frac{\pi - \theta}{2}\right)\] where \(q_0 = 1\) for \(\theta=\pi/3\) and \(q_0\) becomes negative for angles that are negative or greater than \(\pi\). The quality of a triangle is then defined as the minimum of the quality measures associated with its interior angles.

3.3.2 Triangle Split Criterion↩︎

In addition to maximizing triangle quality, it is desirable to enforce a lower bound ensuring that all elements satisfy a minimum quality requirement. To this end, given prescribed minimum and maximum angles \(\theta_{\text{min}}\) and \(\theta_{\text{max}}\), we define a quality measure \(q\) such that negative \(q\) are considered invalid and therefore the corresponding triangle must be refined: \[q(\theta) = \min\left( \frac{\theta - \theta_{\text{min}}}{\pi/3 - \theta_{\text{min}}}, \frac{\theta_{\text{max}} - \theta}{\theta_{\text{max}} - \pi/3} \right)\]

When enforcing a mesh quality criterion, all operations that could violate it must be prevented. For a triangle-based criterion, collapse and split operations must include additional checks to avoid creating low-quality triangles, as detailed in Alg. 37 and 38. Similarly, for an edge-based criterion, edge swaps would also need to be monitored to prevent the introduction of poorly shaped edges.

3.3.3 Triangle Split Feasibility↩︎

When splitting a triangle, a new vertex is introduced to eliminate it. The placement of this vertex is critical: the removed triangle must not reappear after subsequent edge swaps, and the insertion should avoid the creation of poorly shaped elements. Classical Delaunay refinement methods [1], [2] place the vertex at the circumcenter of the triangle. This choice guarantees, by the Delaunay criterion, that the triangle will not reappear, while also maximizing the distance to existing vertices, thereby promoting good element quality. In this context, the computational cost associated with evaluating the intrinsic circumcircle is considered acceptable.

A practical difficulty arising from this choice is that the circumcenter may lie outside the triangle to be removed, requiring the identification of the triangle or edge that contains it. Additional edge swaps are then necessary to recover a Delaunay configuration in which the target triangle is effectively removed. Consequently, intersection checks are performed both after the insertion and during the sequence of edge swaps, similarly to the procedures used for edge split and collapse operations.

Although a Delaunay triangulation is guaranteed to exist, practical limitations may arise. Numerical issues can prevent the computation of circumcenters, and the heuristic Delaunay criterion may fail to eliminate the targeted triangle after edge swaps. Consequently, triangle splitting is not always successful.

3.3.4 Triangle Split Algorithm↩︎

In conventional Delaunay refinement, splitting triangles in order of increasing quality, starting with the worst, has been shown to reduce the total number of mesh elements [3]. A similar approach is adopted here, using an ordered triangle-splitting procedure as outlined in Alg. 22.

3.4 Ordering of the local intrinsic operations↩︎

When optimizing a mesh with respect to both element size and quality, an important factor is the ordering of the local operations. This ordering is chosen to minimize computational cost while maintaining satisfactory results.

The procedure begins by enforcing a Delaunay configuration via edge swaps, ensuring that such a configuration exists. Subsequently, it is beneficial to collapse as many edges as possible, thereby reducing the total number of triangles and edges. Edge splits are then applied, as they are relatively inexpensive, followed by triangle splits on the partially processed mesh. Because each operation is immediately followed by local edge swaps, a global edge swap procedure is only required at the initial stage.

This sequence of operations is then iterated until convergence, yielding a mesh that satisfies the prescribed constraints wherever feasible. The complete procedure is summarized in Alg. 23 and an example of intrinsic meshing is represented in Fig. 24.

4 Geodesic Computations↩︎

The intrinsic meshing approach relies on the ability to compute geodesic distances, shortest paths between vertices, and geodesic circumcenters, i.e., points equidistant from three given vertices on the surface. In this work, these computations are performed using continuous Dijkstra methods, specifically the Mitchell-Mount-Papadimitriou (MMP) or Improved Chen-Han (ICH) algorithms. The differences between these two algorithms are briefly discussed in Sec. 4.1.2. To further enhance the efficiency of the method, we propose several minor modifications to these algorithms, specifically designed to reduce the computational cost for the tasks considered here.

4.1 Geodesic Distance↩︎

As introduced, continuous Dijkstra algorithms compute exact shortest geodesic distances from given source points on vertices, edges, or faces to any point on a simplicial surface. They resemble the classic graph-based Dijkstra algorithm [25], propagating distance information in order of increasing distance. This ensures destinations are reached with minimal distance while avoiding unnecessary propagation. Unlike the graph case, where only scalar values are propagated between vertices, the continuous case requires information sufficient to determine distances to a continuous set of points.

4.1.1 Triangle Unfolding↩︎

None

Distance propagation on a triangulated surface occurs along edges connecting adjacent faces, as these edges constitute the only direct links between faces. Furthermore, knowledge of the distance from a source to an edge is sufficient to determine the distance to any point on the surface, since the remaining distance from a point to the edge can be computed within the corresponding face. Consequently, only distances to edges are considered in the following.

By rotating a face around a shared edge, adjacent faces can be aligned within a common plane on opposite sides of the edge, reducing inter-face distance computation to a planar problem. This procedure can be iteratively applied to unfold sequences of faces. However, in such an unfolding, distances from the source to an edge are valid only within a conic region defined by the extremities of the traversed edges, see Fig. 25.

Each unfolding therefore defines a continuous set of points on an edge, referred to as an interval or window, which is reachable via a geodesic corresponding to that unfolding. Computing the distance to any point within an interval requires only the relative position of the source in the unfolded plane. Also unfolding an adjacent face corresponds to projecting the interval onto the new edge and rotating the source coordinates in the unfolded plane with respect to this edge. Hence, each unfolding can be fully represented by its associated interval, and the process of computing subsequent unfoldings is equivalent to propagating intervals across neighboring faces.

4.1.2 Interval Intersections↩︎

Since the algorithm is applicable to any surface embedded in \(\mathbb{R}^3\), regions with nonplanar curvature can give rise to special cases. One such case occurs when two unfoldings generate overlapping intervals along an edge. There are two primary strategies for handling this situation, which distinguish the MMP and CH approaches.

In MMP, the overlap is resolved by computing the point on the edge that is equidistant to the sources of the intersecting intervals. Then the interval sizes are reduced so that they no longer overlap.

In contrast, CH retains overlapping intervals and defers the selection to the distance computation stage, ultimately choosing the interval that yields the shortest path. This strategy leads CH to generate more intervals, as no interval is pruned by another. The improved CH (ICH) approach introduces a filtering mechanism to discard intervals that will not contribute to the shortest path, thereby achieving performance comparable to MMP.

4.1.3 Geometric Shadow Propagation↩︎

A second special case arises when a window is propagated through a saddle or boundary point. In such situations, certain portions of edges cannot be reached directly, as they lie in a geometric shadow where no direct path exists from the source. Nevertheless, these regions can still be connected via the saddle or boundary point. Consequently, the resulting window is represented with a pseudo-source located at the corresponding position of the saddle or boundary point. Additionally, the window requires an auxiliary parameter \(d\) to track the distance from the pseudo-source to the original source, ensuring that the total distance to the source can still be computed efficiently.

4.1.4 Propagation Order↩︎

The performance difference between CH and ICH or MMP lies in the order in which intervals are propagated: in ICH and MMP, intervals are processed according to the smallest distance first. This ordering not only enables early termination upon reaching the destination but also facilitates effective pruning of intervals.

Therefore, to implement propagation in this distance-priority order, each window is associated with a parameter \(min\), which records the smallest distance from any point of the window to the source. Once all windows with smaller distances have been propagated, the distance to the destination is guaranteed to be minimal.

4.1.5 Performance Considerations↩︎

In the worst case, the number of windows is \(O(n^2)\), as each vertex may split a window into two. Since windows are processed in order of increasing distance, a priority queue is required, leading to an overall complexity of \(O(n^2 \log n)\). In practice, however, the number of windows is significantly smaller due to pruning and intersection effects, and propagation is typically restricted to the vicinity of the destination rather than the entire surface, resulting in much lower computational cost.

On top of that, the information associated with each window is relatively lightweight: its corresponding edge \(e\), the interval range (\(x_{\text{start}}\), \(x_{\text{end}}\)), the pseudo-source position (\(x_{\text{pseudo}}\), \(y_{\text{pseudo}}\)), the distance \(d\) from the pseudo-source to the original source, and the minimal distance value \(min\) within the window.

None

4.2 Geodesic Shortest Path↩︎

The shortest path is recovered by backtracking from the destination to the source, following windows across edges that minimize the distance. This process identifies a sequence of points on edges and vertices forming the path, and consists of iterating over candidate edges and selecting the shortest path through the corresponding windows.

4.2.1 Candidate Edges↩︎

The set of candidate windows depends on the type of surface point considered. For a face, distances are evaluated to windows on its edges. For an edge, they are evaluated on the opposite edges of adjacent faces. At a vertex, only opposite edges are considered, as adjacent edges would introduce ambiguity due to equal distances. The same procedure applies to MMP and ICH, with the difference that more windows are present per edge in ICH since windows are not intersected. However, since intervals remain unchanged in ICH, one can store additional information identifying their origin interval, similarly to Dijkstra’s algorithm on graphs. The backtrace then reduces to following the chain of ancestor intervals.

4.2.2 Distance through a window↩︎

The distance from a point to the source through a window is computed in the unfolded plane, where the pseudo-source is reflected across the edge. The minimal distance is then given either by the direct path through the window or via its endpoints.

None

4.2.3 Performance Considerations↩︎

The backtracing cost is proportional to the number of windows, bounded by \(O(n)\), and more importantly to the number of points \(k\) in the path. Its worst-case complexity is \(O(nk)\), which remains negligible compared to the cost of distance propagation. Consequently, computational efficiency must primarily be achieved during the propagation stage.

4.3 Geodesic A* Search↩︎

Introduced in [33], the A* search algorithm improves upon Dijkstra’s method on graphs by incorporating information about the goal’s location to guide the search. While Dijkstra computes shortest paths to all nodes, A* focuses on a single target, which aligns with our needs. The algorithm prioritizes nodes based on an estimate of the total path length given by the exact distance from the source and a heuristic lower bound on the distance to the destination. Consequently, A* applies to graphs where such admissible heuristics exist.

In our setting, an admissible heuristic is readily available, as the Euclidean distance between two points provides a lower bound on their geodesic distance. Accordingly, the window minimum distance can be modified to reflect the estimated total distance from source to destination through the window. This estimate is obtained by reusing the point-to-source distance computation introduced for backtracing, see Sec. 4.2.2. This adjustment is sufficient to guide the propagation, as the priority queue relies on these minimal values to select windows. Moreover, after the destination is first reached, propagation must continue for windows whose estimated total distance remains smaller than the current best solution, ensuring optimality as in the classical MMP and ICH algorithms.

Fig. 30 compares the A* search with a classical continuous Dijkstra algorithm. The A* method propagates distances selectively toward the destination, reducing the number of windows to only \(5\%\) of the original. This reduction occurs because edges far from the source generate many intervals due to intermediate vertices, most of which are avoided with the guided search. In addition to fewer windows, the priority queue size is smaller, as fewer windows are active at the propagation front. Overall, computing geodesics on a model with a thousand faces requires only \(3\%\) of the time compared to the classical MMP algorithm. This speedup is critical for enabling the intrinsic meshing approach presented in this work.

4.4 Geodesic Circumcenter↩︎

A fundamental component of Delaunay refinement is the circumcenter of a triangle. As discussed in Sec. 3.1.6, the existence and uniqueness of the circumcenter are not guaranteed for arbitrary configurations. However, in the following, such a circumcenter is assumed to exist and be unique.

For surfaces embedded in three-dimensional space, the notion of a circumcircle generalizes naturally as the locus of points equidistant from the triangle’s vertices. This formulation allows geodesic distances to be used to define and compute the circumcenter.

A method for computing circumcenters was proposed in [9], where triangles containing the circumcenter are recursively refined using intervals along the edges until the circumcenter is captured by an edge. However, since only circumcenters are required, rather than the full Voronoi diagram, a simplified strategy is adopted here.

Again, distances are propagated from the three vertices, and candidate triangles potentially containing the circumcenter are identified using edge intervals. Unlike the previous approach, all possible circumcenter positions arising from combinations of edge windows are evaluated. Each candidate is then verified to ensure it lies within the triangle and satisfies the equidistance condition. The only requirement for this method is the ability to compute a circumcenter from three distance windows.

In the planar case, circumcenters are obtained by intersecting edge bisectors, which are straight lines. On curved surfaces, bisectors are generally non-linear: points equidistant from two pseudo-sources satisfy a hyperbolic equation as represented in Fig. 31. Nevertheless, an analytical solution for the intersection of such bisectors can be derived, as detailed in Appendix 7.

To improve the propagation speed of distance intervals toward the circumcenter, a strategy similar to that presented in Section 4.3 can be employed. In this context, however, the final position of the circumcenter is unknown and constitutes the quantity to be determined. As an estimate, one may use the center of the smallest sphere passing through the three points as a target. This heuristic performs well in nearly planar settings, but it becomes unreliable in the presence of strong surface curvature, thereby limiting the efficiency of this optimization.

None

5 Experimental Results↩︎

5.1 Effect of constraints↩︎

Fig. 32 illustrates various test cases for intrinsic remeshing under varying angular and minimum characteristic lengths constraints. A primary observation is that these intrinsic elements, even if topologically defined as triangles, may exhibit complex geometries. Because they are defined as a subdivision of the underlying input manifold, a single intrinsic triangle often spans a collection of original mesh faces leading to discontinuous changes in its shape. However, as the sampling density increases, the geometric complexity of these intrinsic triangles decreases. In this high-density regime, the elements converge toward simpler shapes.

Furthermore, the results demonstrate that imposing strict angular bounds effectively drives adaptive refinement. In regions of high geometric curvature, typically associated with critical simulation phenomena, fitting large intrinsic triangles that satisfy prescribed angular limits is difficult. Consequently, the algorithm naturally produces a denser triangulation in these zones to satisfy the constraints.

Finally, while the prescribed angular bounds are not guaranteed due to the numerical challenges of locating intrinsic circumcenters and avoiding intersecting geodesics on highly curved surfaces, the majority of elements successfully adhere to the constraints. It should be noted, however, that strictly enforcing these bounds significantly increases the total element count.

5.2 Robustness↩︎

As previously discussed, all results were computed using an implementation of the MMP algorithm provided by one of the authors of [32]. The primary modification to this implementation consists in adapting the minimum distance computation to incorporate an A* search variant of the algorithm. In addition, the structure of the output points was revised: instead of returning three-dimensional coordinates, the algorithm now produces parametric coordinates on the surface, which offer a more robust representation. Finally, the numerical tolerances were adjusted to better reflect the limitations of floating-point arithmetic.

To evaluate the robustness of our method, we conducted a large-scale study using the Thingi10k dataset [36]. We employed the dataset’s filtering tools to extract all models that satisfy three strict criteria required by the method: the surfaces must be manifold, closed and consistently oriented. Indeed, the orientation of triangles is used to correctly compute intrinsic angles. The filtering process removes meshes with edges that are not adjacent to two triangles with the same orientation. However, the challenging ‘real-world’ artifacts such as high-aspect-ratio/degenerate triangles, varying sampling densities and complex topologies that are characteristic of the Thingi10k collection are not filtered out. The resulting test suite provides a rigorous benchmark, presented in Table 1 for the robustness of our implementation on topologically sound but geometrically complex surfaces.

The largest possible element size and angular bounds between \(20\) and \(140\) degrees are enforced as constraints throughout the process. It should be noted that the computation of circumcenters may occasionally fail due to numerical inaccuracies. However, such failures can be safely ignored without significantly affecting the overall algorithmic outcome. Consequently, these cases are not classified as failures. In contrast, a run is considered a failure only if a geodesic cannot be computed or if the limits of 2 hours of compute time or 10 GB of memory are exceeded.

Because circumcenters are not always available, the prescribed angle constraints cannot be systematically enforced. Moreover, certain geometries intrinsically prevent the satisfaction of these constraints due to their underlying shape characteristics. As a result, among all successful test cases, only \(27.5\%\) strictly satisfy the prescribed angular bounds, while the remaining cases exhibit one or more angles outside the specified range. Despite this, all successful cases converge within fewer than 25 iterations of the main loop.

5.3 Performance↩︎

All computations were carried out on a Genoa (3.25 GHz) CPU using a single-threaded implementation. The corresponding computation times are presented in Fig. 33. The two main computational bottlenecks of the method are the evaluation of geodesic paths and the computation of circumcenters. Even if these constructs are stored to prevent redundant calculations, they remain the dominant contributors to the overall computational cost, primarily due to the expense of evaluating geodesic distances on the surface.

The performance of the meshing method could thus be substantially enhanced through multithreading or GPU acceleration. A straightforward approach would be to use parallel implementations of continuous Dijkstra algorithms. However, due to the inherently sequential nature of Dijkstra’s ordered propagation, parallelization will only achieve limited speed-up.

A more ambitious strategy would involve performing intrinsic local optimization in parallel. This approach, however, necessitates a more sophisticated implementation to ensure that local optimizations do not interfere with other regions of the mesh being updated simultaneously. In [37], such an implementation was shown to be efficient for large-scale standard meshing methods, while also highlighting the inherent challenges of parallelizing such computations.

5.4 Direct High-Order Meshing↩︎

One important application of intrinsic meshes is the fitting of high-order polynomial elements. This is of significant interest because it provides a direct approach for constructing high-order elements by taking the true geometry and curvature into account from the beginning. This contrasts with classical high-order meshing methods, where linear simplices are first generated and then curved through an optimization procedure to obtain high-order elements [38]. An example is shown in Fig. 34, where the surface of a trabecular bone is first meshed using intrinsic meshing, and the resulting triangles are then fitted with high-order polynomial elements to produce a high-order mesh composed of classical elements suitable for simulation with standard methods.

It should be noted that these results are obtained solely by constraining edge lengths and angles. In Sec. 3.3.1, the notion of triangle quality was introduced to recover a measure of element quality suitable for numerical simulations. Indeed, triangle angles are commonly used as indicators of simulation accuracy in classical numerical methods, although this assumption may not necessarily hold for intrinsic triangulations. Consequently, depending on the intended application of the intrinsic mesh, it may be advantageous to introduce additional quality measures for both triangles and edges. More generally, the method described here is flexible and can accommodate alternative edge or triangle quality criteria.

6 Perspectives↩︎

In this work, we presented a complete method for intrinsic remeshing, from the technical computation of shortest-path geodesics and intrinsic circumcenters, to the development of a remeshing algorithm that leverages local optimization operations to enforce various constraints on the resulting mesh. More sophisticated operations could be considered, such as locally repositioning vertices within their cavities to further improve mesh quality. Such strategies could also enable more advanced remeshing approaches, including optimal Delaunay triangulation [16].

It was also demonstrated that, once an intrinsic triangulation is obtained, it can readily serve as a basis for fitting different types of elements, such as high-order polynomials. Of particular interest, however, is the potential to use these intrinsic meshes directly, without additional processing. Since the mesh preserves the exact geometry, computations performed on it can fully leverage this accuracy.

An important direction for future work is the extension of the proposed method to open surfaces, as its current restriction to closed surfaces constitutes a significant limitation. In addition, the method does not adequately capture sharp features, which may be interpreted as embedded edges on the surface. Consequently, enabling the representation of boundaries and embedded feature edges remains a key challenge. This difficulty stems from the fact that such edges do not necessarily coincide with shortest paths between vertices and therefore may not correspond to geodesics on the surface. As a result, intrinsic remeshing on open surfaces, as well as the treatment of embedded edges, will require specialized handling of these features.

Acknowledgments↩︎

This research was funded in part through the ARC grant for Concerted Research Actions (ARC Discrete-IGA 23/27-08), financed by the Wallonia-Brussels Federation of Belgium. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region. Human trabecular bone samples were provided by M.A. Hartmann and S. Blouin (Ludwig Boltzmann Institute of Osteology, Vienna). Micro-CT scanning was performed by E. Pedrinazzi (Mechanics of Biological and Bio-Inspired Materials, ULiège).

7 Circumcenter Computation With Pseudo-Sources↩︎

As detailed in Sec. 4.4, it is necessary to compute the position of the circumcenter, if it exists, defined by three pseudo-sources lying in the same plane. Each pseudo-source is characterized by its position in the plane and its distance to an associated source. The objective is to determine the coordinates of a point that is equidistant from the three associated sources. This problem reduces to solving the following system of equations for the unknowns \(x\) and \(y\), given three pseudo-sources \(A\), \(B\), and \(C\): \[\begin{align} &d_A + \sqrt{(x - x_A)^2 + (y - y_A)^2} \\ = &d_B + \sqrt{(x - x_B)^2 + (y - y_B)^2} \\ = &d_C + \sqrt{(x - x_C)^2 + (y - y_C)^2} \end{align}\]

A first simplification of the system consists in applying a translation of the coordinate system relative to the first pseudo-source. Additionally, if the pseudo-sources are ordered by increasing distance to their associated sources, the problem can be reformulated in terms of a shifted total distance: \[\begin{gather} \begin{cases} \sqrt{{x'}^2 + {y'}^2} = \, d_B + \sqrt{(x' - x'_B)^2 + (y' - y'_B)^2} \\ \sqrt{{x'}^2 + {y'}^2} = \, d_C + \sqrt{(x' - x'_C)^2 + (y' - y'_C)^2} \end{cases} \end{gather}\] with \(d' = d - d_A\), \(x' = x - x_A\) and \(y' = y - y_A\). A rotation of the coordinate system eliminates an additional parameter. The final step in this simplification is to scale the system by the inverse of \(x'_B\), thereby removing this parameter from the formulation: \[\begin{cases} X = (x'_B x' + y'_B y') \big/ (\sqrt{{x'_B}^2 + {y'_B}^2} x'_B) \\ Y = (- y'_B x' + x'_B y') \big/ (\sqrt{{x'_B}^2 + {y'_B}^2} x'_B) \end{cases}\] In summary, the problem can be expressed in the following simplified form: \[\begin{align} \begin{cases} \sqrt{X^2 + Y^2} = D_B + \sqrt{(X - 1)^2 + Y^2} \\ \sqrt{X^2 + Y^2} = D_C + \sqrt{(X - X_C)^2 + (Y - Y_C)^2} \end{cases} \\\text{with} \quad D_B = d'_B / x'_B \qquad D_C = d'_C / x'_B \\ X_C = (x'_B x'_C + y'_B y'_C) \big/ (\sqrt{{x'_B}^2 + {y'_B}^2} x'_B) \\ Y_C = (x'_B y'_C - y'_B x'_C) \big/ (\sqrt{{x'_B}^2 + {y'_B}^2} x'_B) \\\end{align}\]

Assuming that the parameters \(D_B\) and \(D_C\) are nonzero, the resulting system of two equations with two unknowns describes the intersection of two hyperbolas, each having a pseudo-source as a focus. Furthermore, since both \(D_B\) and \(D_C\) are positive, the solution is confined to the intersection of specific branches of these hyperbolas. As a result, two distinct solutions are expected. By isolating the last term in each equation and then squaring, a linear relationship between \(X\) and \(Y\) can be established. Crucially, this squaring does not introduce spurious solutions, as all terms involved are nonnegative: \[\begin{align} 2 D_B Y_C Y &= 2 (D_C - D_B X_C) X \\ &+ \big( D_C (D_B^2 - 1) + D_B (X_C^2 + Y_C^2 - D_C^2) \big) \end{align}\] Assuming that \(Y_C\) is nonzero, \(Y\) can be substituted from the derived relation into the first equation. This reduces the system to a quadratic equation in \(X\), which yields the two expected solutions. The problem can be solved using symbolic resolution, and the resulting solutions are presented in Table 2.

The above reasoning did not consider the cases where \(D_B\), \(D_C\), or \(Y_C\) are equal to zero. These scenarios correspond, respectively, to the intersection of one line and one hyperbola, the intersection of two lines and the intersection of two hyperbolas whose foci are collinear. Nonetheless, it remains possible to apply a similar approach to determine the solutions for these special cases. All valid solutions are presented in Table 2.

8 Additional pseudo-codes↩︎

References↩︎

---

1. The implementation code is openly available within the official Gmsh repository at https://gitlab.onelab.info/gmsh/gmsh.↩︎

2. Danil Kirsanov (2025). Exact geodesic for triangular meshes (https://www.mathworks.com/matlabcentral /fileexchange/18168-exact-geodesic-for-triangular-meshes), MATLAB Central File Exchange. Retrieved May 16, 2025.↩︎