# Mesh Density and Geodesic Tortuosity in Planar Triangular Tesselations Devoted to Fracture Mechanics

Joffrey Lhonneur, Nawfal Blal, Yann Monerie

► To cite this version:

Joffrey Lhonneur, Nawfal Blal, Yann Monerie. Mesh Density and Geodesic Tortuosity in Planar Triangular Tesselations Devoted to Fracture Mechanics. 2024. hal-03703549v4

HAL Id: hal-03703549

https://hal.science/hal-03703549v4

Preprint submitted on 20 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mesh Density and Geodesic Tortuosity in Planar Triangular Tesselations Devoted to Fracture Mechanics

Joffrey Lhonneur, Nawfal Blal, Yann Monerie

► To cite this version:

Joffrey Lhonneur, Nawfal Blal, Yann Monerie. Mesh Density and Geodesic Tortuosity in Planar Triangular Tesselations Devoted to Fracture Mechanics. Journal of Theoretical, Computational and Applied Mechanics, 2024, pp.1-25. 10.46298/jtcam.9768. hal-03703549v5

HAL Id: hal-03703549

https://hal.science/hal-03703549v5

Submitted on 25 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mesh Density and Geodesic Tortuosity in Planar Triangular Tessellations Devoted to Fracture Mechanics

Joffrey LHONNEUR<fnref n="1" />,<fnref n="2" />, Nawfal BLAL<fnref n="2" />,<fnref n="4" />, and Yann MONERIE<fnref n="2" />,<fnref n="3" />

In fracture mechanics, the mesh sensitivity is a key issue. It is particularly true concerning cohesive-volumetric finite element methods in which the crack path and the overall behavior are respectively influenced by the mesh topology and the mesh density. Poisson-Delaunay tessellations parameters, including the edge length distributions, were widely studied in the literature but very few works concern the mesh density and topology in Delaunay type meshes suitable for finite element simulations, which is of crucial interest for practical use.

Starting from previous results concerning Poisson-Delaunay tessellations and studying in detail the Lloyd relaxation algorithm, we propose estimates for the probability density functions of the edge length and triangle top angles sets. These estimates depend both on the intensity of the underlying point process and on an efficiency index associated to the global quality of the mesh. The global and local accuracies of these estimates are checked for various standard mesh generators. Finally the mesh density and geodesic tortuosity are estimated for standard random or structured triangular meshes typically used in finite element simulations.

These results provide practical formulas to estimate bias introduced by the mesh density and topology on the results of cohesive-volumetric finite element simulations.

Keywords triangular tessellations, edge length distribution, angle distribution, mesh density, mesh tortuosity, Delaunay mesh

## Short nomenclature

### Nomenclature

$\tau$ denotes the intensity of a point process, defined in ??.

$N_E$ denotes the number of edges of a mesh.

$\mathcal{E}$ denotes the efficiency index of a mesh, introduced in ??.

$f_G(\tau; x)$ is the probability density function (pdf) of a quantity $x$ for a graph/tessellation $G$ with intensity $\tau$; this pdf has mean $\mu_G(\tau; x)$ and variance $\sigma_G^2(\tau; x)$.

$Z_G$ is the mesh density of a graph/tessellation $G$ defined as the total length of an edge set divided by the area $A$ of the discretized planar domain.

$T_G(n)$ denotes the tortuosity of a particular continuous path between two nodes of a tessellation $G$ composed by $n$ connected edges.

PD and St stand for Poisson-Delaunay and Standard type meshes respectively.

## 1 Introduction

The geometrical features statistics of Poisson-Delaunay planar triangular tessellations have been studied by numerous authors among which Meijering 1953; Miles 1970; Sibson 1980; Chang et al. 1984; Watanabe 2008; Liu et al. 2009. Their works were mostly devoted to network/graph theories or computational geometry and the context of finite element simulations has been less studied.

<footnote n="1">¹ IRSN, Fontenay-aux-Roses, France</footnote>

<footnote n="2">² MIST Laboratory, IRSN-CNRS-Univ. Montpellier, France</footnote>

<footnote n="3">³ LMGC, Univ. Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="4">⁴ Univ. Lyon, INSA-Lyon, CNRS UMR5259, LaMCoS, F-69621, France</footnote>
**Identifiers**

DOI 10.46298/jtcam.9768

HAL hal-03703549v5

**History**

Received Jul 7, 2022

Accepted May 6, 2024

Published Sept 25, 2024

**Associate Editor**

Julien RÉTHORÉ

**Reviewers**

Jérémie GIRARDOT

Mahmoud SHAQFA

**Open Review**

HAL hal-04653166

**Supplementary Material**

See addendum

**Licence**

CC BY 4.0

©The Authors

# Mesh density and geodesic tortuosity in planar triangular tessellations devoted to fracture mechanics

Joffrey LHONNEUR<fnref n="1" />,<fnref n="2" />, Nawfal BLAL<fnref n="2" />,<fnref n="4" />, and Yann MONERIE<fnref n="2" />,<fnref n="3" />

In fracture mechanics, the mesh sensitivity is a key issue. It is particularly true concerning cohesive-volumetric finite element methods in which the crack path and the overall behavior are respectively influenced by the mesh topology and the mesh density. Poisson-Delaunay tessellations parameters, including the edge length distributions, were widely studied in the literature but very few works concern the mesh density and topology in Delaunay type meshes suitable for finite element simulations, which is of crucial interest for practical use. Starting from previous results concerning Poisson-Delaunay tessellations and studying in detail the Lloyd relaxation algorithm, we propose estimates for the probability density functions of the edge length and triangle top angles sets. These estimates depend both on the intensity of the underlying point process and on an efficiency index associated to the global quality of the mesh. The global and local accuracies of these estimates are checked for various standard mesh generators. Finally the mesh density and geodesic tortuosity are estimated for standard random or structured triangular meshes typically used in finite element simulations. These results provide practical formulas to estimate bias introduced by the mesh density and topology on the results of cohesive-volumetric finite element simulations.

**Keywords:** triangular tessellations, edge length distribution, angle distribution, mesh density, mesh tortuosity, Delaunay mesh

## 1 Introduction

The geometrical features statistics of Poisson-Delaunay planar triangular tessellations have been studied by numerous authors among which (Meijering 1953; Miles 1970; Sibson 1980; Chang et al. 1984; Watanabe 2008; Liu et al. 2009). Their works were mostly devoted to network/graph theories or computational geometry and the context of finite element simulations has been less studied. Tessellations suitable for finite element simulations are built to avoid badly shaped elements in order to preserve the quality of the finite element approximation, see e.g. (Shewchuk 2002) for a review. Roughly speaking, the standard mesh generators lead to less “disordered” tessellations than the Poisson-Delaunay ones. To the best of our knowledge, there exists no estimates of the geometrical features statistics of this type of tessellations, although insights on such estimates could prove useful, for example, for fracture mechanics simulations.

The mesh dependency is indeed a crucial problem in fracture mechanics. In continuum damage mechanics as well as in surface fracture methods (e.g. in DEM simulations (Nguyen et al. 2017; Liu et al. 2021; Manorosoa et al. 2023)), softening behaviors lead to ill-posed problems and the field equations which describe the motion of the body lose hyperbolicity as soon as softening occurs (Lasry et al. 1988; Bourdin et al. 2008). Many regularization methods have been proposed, most of them introducing a characteristic length scale (nonlocal formulations, higher gradient theories, extended or enriched finite element methods, etc.) or a characteristic time scale (volume or surface viscosity), see e.g. (Chaboche et al. 2001; Lu et al. 2009) for a review.

<footnote n="1">¹IRSN, Fontenay-aux-Roses, France</footnote>

<footnote n="2">²MIST Laboratory, IRSN-CNRS-Université de Montpellier, France</footnote>

<footnote n="3">³LMGC, Université de Montpellier, CNRS, France</footnote>

<footnote n="4">⁴Université de Lyon, INSA-Lyon, CNRS UMR5259, LaMCoS, France</footnote>
Tessellations suitable for finite element simulations are built to avoid badly shaped elements in order to preserve the quality of the finite element approximation, see e.g. Shewchuk (2002) for a review. Roughly speaking, the standard mesh generators lead to less “disordered” tessellations than the Poisson-Delaunay ones. To the best of our knowledge, there exists no estimates of the geometrical features statistics of this type of tessellations, although insights on such estimates could prove useful, for example, for fracture mechanics simulations.

The mesh dependency is indeed a crucial problem in fracture mechanics. In continuum damage mechanics as well as in surface fracture methods (e.g. in DEM simulations (Nguyen et al. 2017; Liu et al. 2021; Manorosoa et al. 2023)), softening behaviors lead to ill-posed problems and the field equations which describe the motion of the body lose hyperbolicity as soon as softening occurs (Lasry et al. 1988; Bourdin et al. 2008). Many regularization methods have been proposed, most of them introducing a characteristic length scale (nonlocal formulations, higher gradient theories, extended or enriched finite element methods, etc.) or a characteristic time scale (volume or surface viscosity), see e.g. Chaboche et al. (2001); Lu et al. (2009) for a review.

Two particular tessellations geometrical features are studied in the present work: the *mesh density* and the *mesh geodesic tortuosity*. These quantities play an important role in surface damage approaches and in particular in the Cohesive-Volumetric Finite Element (CVFE) method (Dugdale 1960; Barenblatt 1962; Tijssens et al. 2000; Perales et al. 2008). This method consists in inserting surface cohesive elements (a traction-separation law) between two adjacent elements of a considered tessellation. The interface between the adjacent elements are governed by cohesive zone models which can either be described by *intrinsic* or *extrinsic* laws.

Within intrinsic model, the CZM curve comprises two parts: an artificial elastic part with a finite initial slope and a softening part describing the failure process. They are pre-inserted in the underlying mesh of the simulation leading to inherent crack activations without any need to add additional criteria. On the other hand, the extrinsic models have infinite initial slope, i.e. no initial elastic part, and the cohesive elements are activated, or incorporated, on the fly once the crack criterion is reached.

These cohesive elements can be inserted on a predefined crack path or between each element of a given mesh. The later is often used when the crack path is not known a priori and has to emerge naturally from the loading path. Despite the success of these models, their use suffers from a twofold mesh dependency.

Firstly, in the CVFE method, cracks are forced to propagate along the tessellation elements boundaries. Thus, a dependence on the mesh topology appears (Xu et al. 1994; Tijssens et al. 2000): the computed crack path length is always greater or equal to the theoretical one as a straight path between two nodes of a planar triangular mesh is not always accessible (Pele et al. 2023). As a consequence, the computed overall damage energy, proportional to the crack length (resp. surface) in 2D (resp. 3D) may be overestimated. The *geodesic tortuosity* is introduced in order to quantify this overestimate. It is defined as the ratio between the length of the shortest path between two nodes of the mesh composed of connected edges and the euclidian distance between these two nodes. It is a deterministic variable for structured meshes and a stochastic variable characterized by a mean and a variance for random meshes. Various solutions were proposed in literature to avoid this dependency on the mesh topology. Among others, the extended finite element method (Dolbow et al. 2000), the pin-wheel based meshes (Papoulia et al. 2006) and K-means meshes (Rimoli et al. 2012) are very promising.

Secondly, the CVFE approach exhibits a *mesh size dependency* when considering an *intrinsic* cohesive law characterised by a finite initial cohesive stiffness. The latter is related to the characteristic length scale introduced by cohesive zone models, since they are defined with the help of a surface energy and a peak stress. A rigorous micromechanical-based approach has been proposed by Blal et al. (2011); Blal et al. (2012); Blal et al. (2012) in order to limit this mesh size sensitivity. In this approach, the cohesive zone models are considered as inclusions distributed in a matrix (the set of the bulk elements). The spatial distribution of these inclusions strictly corresponds to the one of the edges (resp. surfaces) in a 2-D (resp. 3-D) mesh. For any random planar meshes with an equiprobable distribution of the orientation of the edges, the spatial distribution of the inclusions is isotropic in space and in orientation. Since these inclusions have zero thickness, they are characterized by a surface density, referred in this paper as *mesh density*.
<table><thead><tr><th colspan="2">Acronyms</th><th></th><th></th></tr></thead><tbody><tr><td>PD</td><td>Poisson-Delaunay mesh</td><td>f<sub>G</sub>(τ; x)</td><td>probability density function of a quantity x for a graph/tessellation G with intensity τ</td></tr><tr><td>PDF</td><td>probability density function</td><td>N<sub>E</sub></td><td>number of edges of a mesh</td></tr><tr><td>St</td><td>Standard mesh</td><td>T<sub>G</sub>(n)</td><td>tortuosity of a particular continuous path between two nodes of a tessellation G composed by n connected edges</td></tr><tr><td colspan="2"><strong>Symbols</strong></td><td>Z<sub>G</sub></td><td>mesh density of a graph/tessellation G defined as the total length of an edge set divided by the area A of the discretized planar domain</td></tr><tr><td>ε</td><td>efficiency index of a mesh, see Section 3.3</td><td></td><td></td></tr><tr><td>μ<sub>G</sub>(τ; x)</td><td>mean of f<sub>G</sub>(τ; x)</td><td></td><td></td></tr><tr><td>σ<sup>2</sup><sub>G</sub>(τ; x)</td><td>variance of f<sub>G</sub>(τ; x)</td><td></td><td></td></tr><tr><td>τ</td><td>intensity of a point process, see Section 2</td><td></td><td></td></tr></tbody></table>

Two particular tessellations geometrical features are studied in the present work: the *mesh density* and the *mesh geodesic tortuosity*. These quantities play an important role in surface damage approaches and in particular in the Cohesive-Volumetric Finite Element (CVFE) method (Dugdale 1960; Barenblatt 1962; Tijssens et al. 2000; Perales et al. 2008). This method consists in inserting surface cohesive elements (a traction-separation law) between two adjacent elements of a considered tessellation. The interface between the adjacent elements are governed by cohesive zone models which can either be described by *intrinsic* or *extrinsic* laws.

Within intrinsic model, the CZM curve comprises two parts: an artificial elastic part with a finite initial slope and a softening part describing the failure process. They are pre-inserted in the underlying mesh of the simulation leading to inherent crack activations without any need to add additional criteria. On the other hand, the extrinsic models have infinite initial slope, i.e. no initial elastic part, and the cohesive elements are activated, or incorporated, on the fly once the crack criterion is reached.

These cohesive elements can be inserted on a predefined crack path or between each element of a given mesh. The later is often used when the crack path is not known a priori and has to emerge naturally from the loading path. Despite the success of these models, their use suffers from a twofold mesh dependency.

Firstly, in the CVFE method, cracks are forced to propagate along the tessellation elements boundaries. Thus, a dependence on the mesh topology appears (Xu et al. 1994; Tijssens et al. 2000): the computed crack path length is always greater or equal to the theoretical one as a straight path between two nodes of a planar triangular mesh is not always accessible (Pele et al. 2023). As a consequence, the computed overall damage energy, proportional to the crack length (resp. surface) in 2D (resp. 3D) may be overestimated. The *geodesic tortuosity* is introduced in order to quantify this overestimate. It is defined as the ratio between the length of the shortest path between two nodes of the mesh composed of connected edges and the Euclidean distance between these two nodes. It is a deterministic variable for structured meshes and a stochastic variable characterized by a mean and a variance for random meshes. Various solutions were proposed in literature to avoid this dependency on the mesh topology. Among others, the extended finite element method (Dolbow et al. 2000), the pin-wheel based meshes (Papoulia et al. 2006) and K-means meshes (Rimoli et al. 2012) are very promising.

Secondly, the CVFE approach exhibits a *mesh size* dependency when considering an *intrinsic* cohesive law characterised by a finite initial cohesive stiffness. The latter is related to the characteristic length scale introduced by cohesive zone models, since they are defined with the help of a surface energy and a peak stress. A rigorous micromechanical-based approach has been proposed by Blal et al. (2011); Blal et al. (2012a); Blal et al. (2012b) in order to limit this mesh size sensitivity. In this approach, the cohesive zone models are considered as inclusions distributed in a matrix (the set of the bulk elements). The spatial distribution of these inclusions strictly corresponds to the one of the edges (resp. surfaces) in a 2-D (resp. 3-D) mesh. For any random planar meshes with an equiprobable distribution of the orientation of the edges, the spatial distribution of the inclusions is isotropic in space and in orientation. Since these inclusions have zero thickness, they are characterized by a surface density, referred in this paper as *mesh density*. It is defined as the total length of the edge set divided by the area of the discretized domain and is directly related to the intensity of the underlying point process. As it is the case for the geodesic tortuosity, it is a deterministic variable for structured meshes and a stochastic variable characterized by a mean and a variance for random meshes.
It is defined as the total length of the edge set divided by the area of the discretized domain and is directly related to the intensity of the underlying point process. As it is the case for the geodesic tortuosity, it is a deterministic variable for structured meshes and a stochastic variable characterized by a mean and a variance for random meshes.

We want to emphasize that the main objective of this paper is to estimate this mesh density and the geodesic tortuosity for standard triangular tessellations suitable for finite element simulations. In turn, most of the content of the present paper concerns geometrical statistics matters. Fracture mechanics is only addressed through practical applications of the mesh density and geodesic tortuosity estimates.

Introducing some notations (??) and starting from standard results concerning Poisson-Delaunay tessellations, we show that the Lloyd relaxation algorithm leads to random triangular tessellations with edge length and triangle top angles distributions having the same means as Poisson-Delaunay tessellations but smaller variances (??). The “quality” of the resulting meshes is also studied using a global efficiency index. Estimates are proposed for both edge length and triangle top angles distributions of any random triangular tessellation with arbitrary efficiency index and intensity of the underlying point process (??). The local robustness and the validity range of these estimates are analyzed in ?? and ?? . Moreover, we show that the edge length and the triangle top angles distributions of a Lloyd relaxed random triangular tessellation, are suitable for standard mesh generators. Finally, the mesh density and geodesic tortuosity are estimated for random tessellations (??). Additional results concerning the mesh density of usual structured tessellations are available in ?? . Practical estimates are exhibited.

2 Notations

A Delaunay triangle is obtained by joining three points of a set of points P if and only if the resulting triangle does not contain any other points of P in the interior of its circumcircle. For a N points generating Poisson point process (intensity τ = N/A , where A is the area covered by the set of points), ?? holds. For a bounded Delaunay tessellation, the number of faces in a triangulation is not exactly twice the number of points and the edge set intensity and the face set intensity given in ?? are only approximations in that case. We assume in the sequel that these approximations are accurate for high intensity point processes even for bounded tessellations.

<table><caption>Table 1: Densities of the point set, of the edge set and of the face set of a Delaunay Tessellation based on a generating Poisson point process with intensity $\tau$. According to Meijering (1953) and Miles (1970).</caption>
  <thead>
    <tr>
      <th>Intensities</th>
      <th>point set</th>
      <th>edge set</th>
      <th>face set</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Delaunay - G<sub>D</sub></td>
      <td>&tau;</td>
      <td>3&tau;</td>
      <td>2&tau;</td>
    </tr>
  </tbody>
</table>

3 Poisson-Delaunay (PD) tessellation and Lloyd relaxed random
triangular tessellation

3.1 Poisson-Delaunay tessellation: probability density functions of edge lengths and triangle top angles

If the set P corresponds to a Poisson point process with intensity $\tau$, Collins (1968), Miles (1970) and Sibson (1980) have shown that the probability density functions of both the length $l$ of the Delaunay edge and the top angle $\alpha$ of a Delaunay triangle are:

$$
f_{\text{PD}}(\tau; l)(x) = \frac{1}{3}\pi x \tau \left( \text{erfc}\left(\frac{1}{2}\sqrt{\pi x \sqrt{\tau}}\right) + x\sqrt{\tau}e^{-\frac{1}{4}\pi x^2 \tau} \right), \quad \text{for } x \ge 0, \qquad (1)
$$

and

$$
f_{\text{PD}}(\alpha)(x) = 4 \left[ (\pi - x) \cos(x) + \sin(x) \right] \frac{\sin(x)}{3\pi}, \quad \text{for } x \in [0, \pi], \qquad (2)
$$

where erfc is the complementary error function.
We want to emphasize that the main objective of this paper is to estimate this mesh density and the geodesic tortuosity for standard triangular tessellations suitable for finite element simulations. In turn, most of the content of the present paper concerns geometrical statistics matters. Fracture mechanics is only addressed through practical applications of the mesh density and geodesic tortuosity estimates.

Introducing some notations, see Section 2, and starting from standard results concerning Poisson-Delaunay tessellations, we show that the Lloyd relaxation algorithm leads to random triangular tessellations with edge length and triangle top angles distributions having the same means as Poisson-Delaunay tessellations but smaller variances, see Section 3. The “quality” of the resulting meshes is also studied using a global efficiency index. Estimates are proposed for both edge length and triangle top angles distributions of any random triangular tessellation with arbitrary efficiency index and intensity of the underlying point process, in Section 4. The local robustness and the validity range of these estimates are analyzed in Appendices A and B. Moreover, we show that the edge length and the triangle top angles distributions of a Lloyd relaxed random triangular tessellation, are suitable for standard mesh generators. Finally, the mesh density and geodesic tortuosity are estimated for random tessellations in Section 5. Additional results concerning the mesh density of usual structured tessellations are available in Appendix C. Practical estimates are exhibited.

## 2 Notations

A Delaunay triangle is obtained by joining three points of a set of points P if and only if the resulting triangle does not contain any other points of P in the interior of its circumcircle. For a N points generating Poisson point process (intensity $\tau = N/A$, where A is the area covered by the set of points), Table 1 holds. For a bounded Delaunay tessellation, the number of faces in a triangulation is not exactly twice the number of points and the edge set intensity and the face set intensity given in Table 1 are only approximations in that case. We assume in the sequel that these approximations are accurate for high intensity point processes even for bounded tessellations.

<table><caption>**Table 1** Densities of the point set, of the edge set and of the face set of a Delaunay Tessellation based on a generating Poisson point process with intensity $\tau$. According to (Meijering 1953; Miles 1970).</caption><thead><tr><th>Intensities</th><th>point set</th><th>edge set</th><th>face set</th></tr></thead><tbody><tr><td>Delaunay - G<sub>D</sub></td><td>τ</td><td>3τ</td><td>2τ</td></tr></tbody></table>

## 3 Poisson-Delaunay tessellation and Lloyd relaxed random triangular tessellation

### 3.1 Poisson-Delaunay tessellation: probability density functions of edge lengths and triangle top angles

If the set P corresponds to a Poisson point process with intensity $\tau$, Collins (1968); Miles (1970); Sibson (1980) showed that the probability density functions of both the length l of the Delaunay edge and the top angle $\alpha$ of a Delaunay triangle are

$$f_{\text{PD}}(\tau; l)(x) = \frac{1}{3}\pi x \tau (\text{erfc}(\frac{1}{2}\sqrt{\pi x}\sqrt{\tau}) + x\sqrt{\tau}e^{-\frac{1}{4}\pi x^2\tau}), \quad \text{for } x \ge 0, \qquad (1)$$

and

$$f_{\text{PD}}(\alpha)(x) = \frac{4}{3\pi}[(\pi - x)\cos x + \sin x]\sin x, \quad \text{for } x \in [0, \pi], \qquad (2)$$

where erfc is the complementary error function. The corresponding means and variances read

$$\mu_{\text{PD}}(\tau; l) = \frac{32}{9\pi\sqrt{\tau}} \quad (\text{length mean}), \quad \sigma_{\text{PD}}^2(\tau; l) = \frac{405\pi - 32^2}{81\pi^2\tau} \quad (\text{length variance}), \qquad (3)$$
The corresponding means and variances read:

$$ \mu_{\text{PD}}(\tau; l) = \frac{32}{9\pi\sqrt{\tau}} \quad (\text{length mean}), \quad \sigma_{\text{PD}}^2(\tau; l) = \frac{405\pi - 32^2}{81\pi^2\tau} \quad (\text{length variance}), \quad (3) $$

and

$$ \mu_{\text{PD}}(\alpha) = \frac{\pi}{3} \quad (\text{top angle mean}), \quad \sigma_{\text{PD}}^2(\alpha) = \frac{\pi^2}{9} - \frac{5}{6} \quad (\text{top angle variance}). \quad (4) $$

These pdf are represented in ?? and compared to histograms of a particular PD tessellation.

<figure><img src="image_2.png" /><figcaption>Figure 1: Left: Delaunay tessellation of a Poisson point process ($\tau = 1000$). Middle: Density function of the length $l$ of a typical Delaunay edge. Right: Density function of a top angle $\alpha$ of a typical Delaunay triangle; black lines correspond to ?? and ??.</figcaption></figure>

We underline that both the mean and the variance of the edge lengths depend on $\tau$. However, the triangle top angle distribution does not depend on the intensity of the point process. With the same number of points generated by the point process, the larger the domain, the greater the variance of the edge length distribution.

## 3.2 Lloyd relaxation

The Lloyd algorithm (Lloyd 1982) is an iterative scheme for computing Centroidal Voronoi tessellations. Centroidal Voronoi tessellations are Voronoi tessellations of a bounded geometric domain $\Omega$ such that the generating points of the tessellations are also the centroids (mass centers) of the corresponding Voronoi regions. The Lloyd algorithm consists in: i) selecting an initial set of points (seeds) according to a point process, ii) constructing the Voronoi tessellation associated to these seeds, iii) computing the mass centroids of the constructed Voronoi zones and iv) consider these centroids as the seeds of the new Voronoi tessellation. The process is iterated, see ???. At convergence, the obtained Voronoi cells compose the so-called Centroidal Voronoi tessellation. The dual of this centroidal Voronoi tessellation is as close as possible to a tessellation of equilateral triangles.

Since Lloyd iterations lead to triangles as close as possible to equilateral ones and knowing that the mean of triangle top angles in a PD tessellation corresponds to the top angle of equilateral triangles, see ??, we consider that the triangle top angles mean remains unchanged through the Lloyd iterations. The same goes with the edge length mean following Chang et al. (1984): denoting by the subscript $L(i)$ parameters of the tessellation obtained by the Lloyd algorithm at iteration $i$:

$$ \mu_{L(i)}(\tau; l) = \mu_{\text{PD}}(\tau; l) \quad \text{and} \quad \mu_{L(i)}(\alpha) = \mu_{\text{PD}}(\alpha). \quad (5) $$

Despite its wide application in many problems, the convergence of the Lloyd algorithm still remains an open issue. Numerous (local and global) convergence properties are presented in Du et al. (2006) and in the references cited therein. The Lloyd algorithm may be viewed as a fixed point iteration of a mapping from a set of distinct seeds to the mass centers. The global convergence is only guaranteed in one-dimensional space for any positive and smooth density
and

$$ \mu_{\text{PD}}(\alpha) = \frac{\pi}{3} \quad (\text{top angle mean}), \quad \sigma_{\text{PD}}^2(\alpha) = \frac{\pi^2}{9} - \frac{5}{6} \quad (\text{top angle variance}). \qquad (4) $$

These pdf are represented in Figure 1 and compared to histograms of a particular PD tessellation.

<figure><img src="image_2.png" /><figcaption>**Figure 1** Geometrical statistics of a PD tessellation. Left: Delaunay tessellation of a Poisson point process ($\tau = 1000$). Middle: Density function of the length $l$ of a typical Delaunay edge. Right: Density function of a top angle $\alpha$ of a typical Delaunay triangle; black lines correspond to Equations (1) and (2).</figcaption></figure>

We underline that both the mean and the variance of the edge lengths depend on $\tau$. However, the triangle top angle distribution does not depend on the intensity of the point process. With the same number of points generated by the point process, the larger the domain, the greater the variance of the edge length distribution.

## 3.2 Lloyd relaxation

The Lloyd algorithm (Lloyd 1982) is an iterative scheme for computing Centroidal Voronoi tessellations. Centroidal Voronoi tessellations are Voronoi tessellations of a bounded geometric domain $\Omega$ such that the generating points of the tessellations are also the centroids (mass centers) of the corresponding Voronoi regions. The Lloyd algorithm consists in: i) selecting an initial set of points (seeds) according to a point process, ii) constructing the Voronoi tessellation associated to these seeds, iii) computing the mass centroids of the constructed Voronoi zones and iv) consider these centroids as the seeds of the new Voronoi tessellation. The process is iterated, see Figure 2. At convergence, the obtained Voronoi cells compose the so-called *Centroidal Voronoi tessellation*. The dual of this centroidal Voronoi tessellation is as close as possible to a tessellation of equilateral triangles.

Since Lloyd iterations lead to triangles as close as possible to equilateral ones and knowing that the mean of triangle top angles in a PD tessellation corresponds to the top angle of equilateral triangles, see Equation (4), we consider that the triangle top angles mean remains unchanged through the Lloyd iterations. Following (Chang et al. 1984), the same goes with the edge length mean

$$ \mu_{L(i)}(\tau; l) = \mu_{\text{PD}}(\tau; l) \quad \text{and} \quad \mu_{L(i)}(\alpha) = \mu_{\text{PD}}(\alpha). \qquad (5) $$

where the subscript $L(i)$ denotes parameters of the tessellation obtained by the Lloyd algorithm at iteration $i$.

Despite its wide application in many problems, the convergence of the Lloyd algorithm still remains an open issue. Numerous (local and global) convergence properties are presented in (Du et al. 2006) and in the references cited therein. The Lloyd algorithm may be viewed as a fixed point iteration of a mapping from a set of distinct seeds to the mass centers. The global convergence is only guaranteed in one-dimensional space for any positive and smooth density function. For multidimensional spaces, no rigorous convergence can be pronounced. Moreover, when the convergence holds, the speed of convergence is very low and can be significantly improved using dedicated algorithms (Liu et al. 2009).
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Lloyd relaxation algorithm. Top row: successive Lloyd relaxations of a Voronoi tessellation (intensity $\tau = 500$). Second row: dual graph. Third row: pdf of the edge length in the dual graph (bar charts). Bottom row: pdf of the triangle top angle in the dual graph (bar charts). Black lines over bar charts are estimates (??) and (??). From left to right various iterations of the Lloyd algorithm: initial tessellation ($\bar{\mathcal{E}} = 0.005$), 2d iteration ($\bar{\mathcal{E}} = 0.41$), 30th iteration ($\bar{\mathcal{E}} = 0.70$). Retained boundary elements may introduce a small bias. Thus edge lengths and summit angles of triangles containing a boundary node were not considered for computing efficiency indexes and establishing bar charts.</figcaption></figure>

function. For multidimensional spaces, no rigorous convergence can be pronounced. Moreover, when the convergence holds, the speed of convergence is very low and can be significantly improved using dedicated algorithms (Liu et al. 2009).

However, we show in the sequel that Lloyd relaxed tessellations lead to high quality triangles. For a sufficient number of iterations, this quality is similar to the quality of triangular mesh generators devoted to finite element computations. This quality is estimated with a global efficiency index. The associated variances of both edge length and triangle top angle sets are also estimated.
<figure><img src="image_3.png" /><figcaption>**Figure 2** Lloyd relaxation algorithm. Top row: successive Lloyd relaxations of a Voronoi tessellation (intensity $\tau = 500$). Second row: dual graph. Third row: pdf of the edge length in the dual graph (bar charts). Bottom row: pdf of the triangle top angle in the dual graph (bar charts). Black lines over bar charts are estimates (14) and (15). From left to right various iterations of the Lloyd algorithm: initial tessellation ($\tilde{\mathcal{E}} = 0.005$), 2d iteration ($\tilde{\mathcal{E}} = 0.41$), 3dth iteration ($\tilde{\mathcal{E}} = 0.70$). Retained boundary elements may introduce a small bias. Thus edge lengths and summit angles of triangles containing a boundary node were not considered for computing efficiency indexes and establishing bar charts.</figcaption></figure>

However, we show in the sequel that Lloyd relaxed tessellations lead to high quality triangles. For a sufficient number of iterations, this quality is similar to the quality of triangular mesh generators devoted to finite element computations. This quality is estimated with a global efficiency index. The associated variances of both edge length and triangle top angle sets are also estimated.

### 3.3 Efficiency index and variance of edge length and triangle top angle sets in PD tessellations improved by the Lloyd relaxation algorithm

The quality of a planar triangular mesh suitable for finite element calculations can be evaluated by two main ways. The first one concerns a topological measure of this quality: spatial homogeneity of the area of each triangle or of the edge length set, respect of minimal or maximal angles, etc.
### 3.3 Efficiency index and variance of edge length and triangle top angle sets in PD tessellations improved by the Lloyd relaxation algorithm

The quality of a planar triangular mesh suitable for finite element calculations can be evaluated by two main ways. The first one concerns a topological measure of this quality: spatial homogeneity of the area of each triangle or of the edge length set, respect of minimal or maximal angles, etc. The second one refers directly to the quality of the finite element approximation. Various authors have thus proposed some interpolation quality measures, see Shewchuk (2002) for example. Since this later type of measure depends on the degree of interpolation in the finite element strategy, we focussed our attention on topological indexes.

Following George et al. (2000), the efficiency index of a triangular tessellation reads:

$$ \mathcal{E} = e^{\langle \eta \rangle}, \quad \text{with} \quad \left\{ \begin{array}{ll} \eta = \bar{l} - 1 & \text{if } \bar{l} < 1 \\ \eta = 1/\bar{l} - 1 & \text{if } \bar{l} \ge 1 \end{array} \right. \qquad (6) $$

where $\bar{l}$ is the adimensional length of an arbitrary edge $e_i$ ($\rho(\mathbf{x})$ is a prescribed size field):

$$ \bar{l} = \int_{e_i} \frac{1}{\rho(\mathbf{x})} dl, \qquad (7) $$

and $\langle \gamma \rangle$ denotes the average of any quantity $\gamma = \{\gamma_1, \dots, \gamma_{N_E}\}$ defined on the set of edges (according to ??, $N_E = 3\tau A$ if $\tau$ is the intensity of the underlying point process over the area $A = |\Omega|$):

$$ \langle \gamma \rangle = \frac{1}{N_E} \sum_{i=1}^{N_E} \gamma_i. \qquad (8) $$

The efficiency index $\mathcal{E}$ belongs to $[0, 1]$ and should be as close as possible to 1 in order to have a mesh suitable for finite element simulations.

In the sequel, the size field $\rho(\mathbf{x})$ is chosen uniform: $\rho(\mathbf{x}) = \mu_{\text{PD},l}(\tau)$, $\forall \mathbf{x} \in \Omega$. Propagating the probability density function of the edge length ?? through the efficiency index ??, one obtains the PD tessellation approximated efficiency index ($\forall \tau$):

$$ \mathcal{E}_{\text{PD}} \approx 0.725. \qquad (9) $$

In the sequel, efficiency indexes, edge length and triangle top angle variances are computed. In order to avoid bias induced by boundary elements, only triangles whose summit are inner tessellation nodes have been considered for computing these values.

The Lloyd algorithm leads to a triangular tessellation whose triangle elements are as close as possible to a tiling of equilateral triangles. As a consequence, both the edge length and triangle top angle variances diminish with Lloyd iterations (see ?? on the left). Thus Lloyd relaxations improves the efficiency index, as the efficiency index increases when the edge length variance decreases (see ?? on the right).

A semi-empirical relationship can be exhibited between the variances and the efficiency for the Lloyd algorithm. Starting from a PD tessellation and remarking that for a structured tessellation based on equilateral triangles the efficiency index is optimal $\mathcal{E}^+ = 1$ and both the edge length and triangle top angle variances vanish $\sigma(\tau, l)^2 = 0$ and $\sigma(\alpha)^2 = 0$, we can derive semi-empirical relationships (see ??):

$$ \frac{\sigma_{L(i)}(\tau; l)^2}{\sigma_{\text{PD}}(\tau; l)^2} \approx \frac{(1 - \bar{\mathcal{E}}_{L(i)})^2}{1 + \bar{\mathcal{E}}_{L(i)}} \quad \text{and} \quad \frac{\sigma_{L(i)}(\alpha)^2}{\sigma_{\text{PD}}(\alpha)^2} \approx \frac{(1 - \bar{\mathcal{E}}_{L(i)})^{3/2}}{(1 + \bar{\mathcal{E}}_{L(i)})^{1/2}}, \qquad (10) $$

where $\bar{\mathcal{E}} \in [0, 1]$ denotes the reduced efficiency index of any index $\mathcal{E}$ greater than or equal to $\mathcal{E}_{\text{PD}}$:

$$ \bar{\mathcal{E}} = \frac{\mathcal{E} - \mathcal{E}_{\text{PD}}}{\mathcal{E}^+ - \mathcal{E}_{\text{PD}}} = \frac{\mathcal{E} - \mathcal{E}_{\text{PD}}}{1 - \mathcal{E}_{\text{PD}}}. \qquad (11) $$
The second one refers directly to the quality of the finite element approximation. Various authors have thus proposed some interpolation quality measures, see (Shewchuk 2002) for example. Since this later type of measure depends on the degree of interpolation in the finite element strategy, we focussed our attention on topological indexes.

Following Frey et al. (2000), the efficiency index of a triangular tessellation reads

$$ \mathcal{E} = e^{\langle \eta \rangle} \quad \text{with} \quad \begin{cases} \eta = \bar{l} - 1 & \text{if } \bar{l} < 1 \\ \eta = 1/\bar{l} - 1 & \text{if } \bar{l} \ge 1 \end{cases} \qquad (6) $$

where

$$ \bar{l} = \int_{e_i} \frac{1}{\rho(\mathbf{x})} d\bar{l}, \qquad (7) $$

is the adimensional length of an arbitrary edge $e_i$ ($\rho(\mathbf{x})$ is a prescribed size field) and

$$ \langle \gamma \rangle = \frac{1}{N_E} \sum_{i=1}^{N_E} \gamma_i \qquad (8) $$

denotes the average of any quantity $\gamma = \{\gamma_1, \dots, \gamma_{N_E}\}$ defined on the set of edges (according to Table 1, $N_E = 3\tau A$ if $\tau$ is the intensity of the underlying point process over the area $A = |\Omega|$).

The efficiency index $\mathcal{E}$ belongs to $[0, 1]$ and should be as close as possible to 1 in order to have a mesh suitable for finite element simulations.

In the sequel, the size field $\rho(\mathbf{x})$ is chosen uniform: $\rho(\mathbf{x}) = \mu_{\text{PD},l}(\tau)$, $\forall \mathbf{x} \in \Omega$. Propagating the probability density function of the edge length Equation (1) through the efficiency index Equation (6), one obtains the PD tessellation approximated efficiency index ($\forall \tau$):

$$ \mathcal{E}_{\text{PD}} \approx 0.725. \qquad (9) $$

In the sequel, efficiency indexes, edge length and triangle top angle variances are computed. In order to avoid bias induced by boundary elements, only triangles whose summit are inner tessellation nodes have been considered for computing these values.

The Lloyd algorithm leads to a triangular tessellation whose triangle elements are as close as possible to a tiling of equilateral triangles. As a consequence, both the edge length and triangle top angle variances diminish with Lloyd iterations (see Figure 3, left). Thus Lloyd relaxations improves the efficiency index, as the efficiency index increases when the edge length variance decreases (see Figure 3, right).

<figure><img src="image_4.png" /><figcaption>**Figure 3** The progressive structuration of a mesh during a Lloyd relaxation. Left: edge length (black marks) and triangle top angle (gray marks) variances for the Lloyd algorithm (normalized by the corresponding variances of the Poisson-Delaunay tessellation respectively $\sigma_{\text{PD}}^2(\tau, l)$ and $\sigma_{\text{PD}}^2(\alpha)$) as a function of the successive iterations. Right: efficiency index $\mathcal{E}$ of the Lloyd algorithm as a function of the successive iterations ($\tau = 300$).</figcaption></figure>

A semi-empirical relationship can be exhibited between the variances and the efficiency for the Lloyd algorithm. Starting from a PD tessellation and remarking that for a structured tessellation based on equilateral triangles the efficiency index is optimal $\mathcal{E}^+ = 1$ and both the
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Left: edge length (black marks) and triangle top angle (gray marks) variances for the Lloyd algorithm (normalized by the corresponding variances of the Poisson-Delaunay tessellation respectively $\sigma_{PD}^2(\tau, l)$ and $\sigma_{PD}^2(\alpha)$) as a function of the successive iterations. Right: efficiency index $\mathcal{E}$ of the Lloyd algorithm as a function of the successive iterations ($\tau = 300$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Edge length and triangle top angle variances for the Lloyd algorithm (normalized by the corresponding variances of the Poisson-Delaunay tessellation respectively $\sigma_{PD}^2(\tau, l)$ and $\sigma_{PD}^2(\alpha)$) as a function of the reduced efficiency index $\bar{\mathcal{E}}$. Points: numerical results starting from 10 different PD tessellations ($\tau = 1000$). Lines: semi-empirical relationships (??).</figcaption></figure>

The corresponding means and variances of edge length and triangle top angle sets of a Lloyd tessellation read (combining ??, ??, ?? and ??):

$$ \mu_{L(i)}(\tau; l) = \frac{32}{9\pi\sqrt{\tau}} \quad \text{and} \quad \sigma_{L(i)}^2(\tau; l) \approx \frac{(1 - \bar{\mathcal{E}}_{L(i)})^2}{1 + \bar{\mathcal{E}}_{L(i)}} \left( \frac{405\pi - 32^2}{81\pi^2\tau} \right), \quad (12) $$

and

$$ \mu_{L(i)}(\alpha) = \frac{\pi}{3} \quad \text{and} \quad \sigma_{L(i)}^2(\alpha) \approx \frac{(1 - \bar{\mathcal{E}}_{L(i)})^{3/2}}{(1 + \bar{\mathcal{E}}_{L(i)})^{1/2}} \left( \frac{\pi^2}{9} - \frac{5}{6} \right). \quad (13) $$

The next section proposes a probability density functions for both edge length and triangle top angle sets for standard triangular random meshes. These meshes try to avoid as much as possible distorted elements to ensure the accuracy of the finite element method.

# 4 Edge length and triangle top angle distributions in standard triangular random meshes

## 4.1 Standard meshing strategies

In the huge literature corresponding to mesh generators, three main approaches can be exhibited for unstructured (or random triangular) meshes (octree methods are not taking into account):

* *Delaunay tessellations.* Starting from a Poisson point process or using iterative George-Frey-like algorithms, this situation leads to $f_{PD}(\tau; l)$.

* *Frontal algorithms.* Using advancing-front techniques, meshes are generated starting from boundaries and adding progressively some “ideal” points as Steiner-like points (Rebay 1993; Lohner et al. 1988). The added points are chosen in order to have a new sequence
edge length and triangle top angle variances vanishe $\sigma(\tau, l)^2 = 0$ and $\sigma(\alpha)^2 = 0$, we can derive semi-empirical relationships (see Figure 4):

$$ \frac{\sigma_{L(i)}(\tau; l)^2}{\sigma_{PD}(\tau; l)^2} \approx \frac{(1 - \bar{\mathcal{E}}_{L(i)})^2}{1 + \bar{\mathcal{E}}_{L(i)}} \quad \text{and} \quad \frac{\sigma_{L(i)}(\alpha)^2}{\sigma_{PD}(\alpha)^2} \approx \frac{(1 - \bar{\mathcal{E}}_{L(i)})^{3/2}}{(1 + \bar{\mathcal{E}}_{L(i)})^{1/2}}, \qquad (10) $$

where $\bar{\mathcal{E}} \in [0, 1]$ denotes the reduced efficiency index of any index $\mathcal{E}$ greater than or equal to $\mathcal{E}_{PD}$:

$$ \bar{\mathcal{E}} = \frac{\mathcal{E} - \mathcal{E}_{PD}}{\mathcal{E}^{+} - \mathcal{E}_{PD}} = \frac{\mathcal{E} - \mathcal{E}_{PD}}{1 - \mathcal{E}_{PD}}. \qquad (11) $$

<figure><img src="image_5.png" /><figcaption>**Figure 4** Edge length and triangle top angle variances for the Lloyd algorithm (normalized by the corresponding variances of the Poisson-Delaunay tessellation respectively $\sigma_{PD}^2(\tau, l)$ and $\sigma_{PD}^2(\alpha)$) as a function of the reduced efficiency index $\bar{\mathcal{E}}$. Points: numerical results starting from 10 different PD tessellations ($\tau = 1000$). Lines: semi-empirical relationships (10).</figcaption></figure>

The corresponding means and variances of edge length and triangle top angle sets of a Lloyd tessellation read (combining Equations (3) to (5) and (10)):

$$ \mu_{L(i)}(\tau; l) = \frac{32}{9\pi\sqrt{\tau}} \quad \text{and} \quad \sigma_{L(i)}^2(\tau; l) \approx \frac{(1 - \bar{\mathcal{E}}_{L(i)})^2}{1 + \bar{\mathcal{E}}_{L(i)}} \left( \frac{405\pi - 32^2}{81\pi^2\tau} \right), \qquad (12) $$

and

$$ \mu_{L(i)}(\alpha) = \frac{\pi}{3} \quad \text{and} \quad \sigma_{L(i)}^2(\alpha) \approx \frac{(1 - \bar{\mathcal{E}}_{L(i)})^{3/2}}{(1 + \bar{\mathcal{E}}_{L(i)})^{1/2}} \left( \frac{\pi^2}{9} - \frac{5}{6} \right). \qquad (13) $$

The next section proposes a probability density functions for both edge length and triangle top angle sets for standard triangular random meshes. These meshes try to avoid as much as possible distorted elements to ensure the accuracy of the finite element method.

# 4 Edge length and triangle top angle distributions in standard triangular random meshes

## 4.1 Standard meshing strategies

In the huge literature corresponding to mesh generators, three main approaches can be exhibited for unstructured (or random triangular) meshes (octree methods are not taking into account):

* *Delaunay tessellations.* Starting from a Poisson point process or using iterative George-Frey-like algorithms, this situation leads to $f_{PD}(\tau; l)$.

* *Frontal algorithms.* Using advancing-front techniques, meshes are generated starting from boundaries and adding progressively some “ideal” points as Steiner-like points (Rebay 1993; Lohner et al. 1988). The added points are chosen in order to have a new sequence of elements as close as possible as those of the previous sequence corresponding to the advancing-front (see e.g. Oliveira (1988) for a review). Far away from the boundaries of the domain $\partial\Omega$, these frontal algorithms lead to semi-structured meshes with hexagonal shape.

* *Ad hoc optimized meshes.* Each mesh generator has its own strategy to construct elements having a good shape. Mathematical indicators of this good shape can be local (minimal angle, minimal aspect ratio, etc.) or global as Equation (6). As for Delaunay tessellations, these mesh generators lead to random meshes.
of elements as close as possible as those of the previous sequence corresponding to the advancing-front (see e.g. De Oliveira (1988) for a review). Far away from the boundaries of the domain $\partial\Omega$, these frontal algorithms lead to semi-structured meshes with hexagonal shape.

* Ad hoc optimized meshes. Each mesh generator has its own strategy to construct elements having a good shape. Mathematical indicators of this good shape can be local (minimal angle, minimal aspect ratio, etc.) or global as ?? . As for Delaunay tessellations, these mesh generators lead to random meshes.

Since, we have here interest for random meshes, we consider in the sequel the “ad hoc optimized meshes” and in particular the so-called “MeshAdapt” strategy proposed by Gmsh (Geuzaine et al. 2009) as a representative strategy of triangular random meshes generators. Starting from a PD tessellation, the “MeshAdapt” strategy consists in four local mesh modifications:

* Edge Splitting: long edges ($\bar{l}_{e_i} > 1.4$) are splitted if the two new obtained edges will have a sufficient adimensional size (min($\bar{l}_{e_j}$, $\bar{l}_{e_k}$) > 0.7).

* Edge Collapsing: short edges ($\bar{l}_{e_i} < 0.7$) are collapsed.

* Edge Swapping: an edge between two triangles is swapped if the efficiency index ?? will improved.

* Triangle re-positioning: each triangle (vertex) is moved inside the cavity made of all its surrounding triangles in order to maximize the worst element quality (Freitag et al. 1997).

## 4.2 Edge length and triangle top angle distributions in random triangular tessellations with any efficiency index and intensity

The Lloyd relaxation algorithm depicted in ?? and ?? shows that when the efficiency index of a random triangular mesh increases, the probability density functions of the edge length and triangle top angle sets seem to be contracted and shifted while globaly keeping their shapes. Moreover edge lengths and triangle top angles means are preserved. We assume that this observation is valid at any intensity $\tau$ and any efficiency index $\mathcal{E}$ and we propose to define the pdf of any random triangular tessellation edge length and triangle top angle sets (referred with the subscript St for Standard mesh generators) as follows:

$$f_{\text{St}}(\bar{\mathcal{E}}, \tau; l)(x) = k(\bar{\mathcal{E}}; l)f_{\text{PD}}(\tau; l)(\tilde{x}(\bar{\mathcal{E}}, x; l)), \quad (14)$$

and

$$f_{\text{St}}(\bar{\mathcal{E}}; \alpha)(x) = k(\bar{\mathcal{E}}; \alpha)f_{\text{PD}}(\alpha)(\tilde{x}(\bar{\mathcal{E}}, x; \alpha)), \quad (15)$$

where $f_{\text{PD}}(\tau; l)$ and $f_{\text{PD}}(\alpha)$ are prolonged to be defined over the whole set of real numbers:

$$\forall x < 0, f_{\text{PD}}(\tau; l)(x) = 0 \quad \text{and} \quad \forall x \in ]-\infty; 0[ \cup ]\pi; +\infty[, f_{\text{PD}}(\alpha)(x) = 0, \quad (16)$$

with

$$\begin{cases} \tilde{x}(\bar{\mathcal{E}}, x; l) = k(\bar{\mathcal{E}}; l)x + p(\bar{\mathcal{E}}; l) \\ \tilde{x}(\bar{\mathcal{E}}, x; \alpha) = k(\bar{\mathcal{E}}; \alpha)x + p(\bar{\mathcal{E}}; \alpha) \end{cases}, \quad (17)$$

the contraction parameters being defined as follows:

$$k(\bar{\mathcal{E}}; l) = \frac{\sigma_{\text{PD}}(\tau; l)}{\sigma_{\text{St}}(\bar{\mathcal{E}}, \tau; l)} = \frac{(1 + \bar{\mathcal{E}})^{1/2}}{1 - \bar{\mathcal{E}}} \quad (18)$$

and

$$k(\bar{\mathcal{E}}, \alpha) = \frac{\sigma_{\text{PD}}(\alpha)}{\sigma_{\text{St}}(\bar{\mathcal{E}}, \alpha)} = \frac{(1 + \bar{\mathcal{E}})^{1/4}}{(1 - \bar{\mathcal{E}})^{3/4}}, \quad (19)$$

and the shift parameters by:

$$p(\bar{\mathcal{E}}, \tau; l) = (1 - k(\bar{\mathcal{E}}; l))\mu_{\text{PD}}(\tau; l) \quad (20)$$
Since, we have here interest for random meshes, we consider in the sequel the “ad hoc optimized meshes” and in particular the so-called “MeshAdapt” strategy proposed by Gmsh (Geuzaine et al. 2009) as a representative strategy of triangular random meshes generators. Starting from a PD tessellation, the “MeshAdapt” strategy consists in four local mesh modifications:

* Edge Splitting: long edges ($\tilde{l}_{e_i} > 1.4$) are split if the two new obtained edges will have a sufficient adimensional size (min($\tilde{l}_{e_j}, \tilde{l}_{e_k}$) > 0.7).

* Edge Collapsing: short edges ($\tilde{l}_{e_i} < 0.7$) are collapsed.

* Edge Swapping: an edge between two triangles is swapped if the efficiency index Equation (6) will improved.

* Triangle re-positioning: each triangle (vertex) is moved inside the cavity made of all its surrounding triangles in order to maximize the worst element quality (Freitag et al. 1997).

## 4.2 Edge length and triangle top angle distributions in random triangular tessellations with any efficiency index and intensity

The Lloyd relaxation algorithm depicted in Figures 2 and 3 shows that when the efficiency index of a random triangular mesh increases, the probability density functions of the edge length and triangle top angle sets seem to be contracted and shifted while globally keeping their shapes. Moreover edge lengths and triangle top angles means are preserved. We assume that this observation is valid at any intensity $\tau$ and any efficiency index $\mathcal{E}$ and we propose to define the pdf of any random triangular tessellation edge length and triangle top angle sets (referred with the subscript St for Standard mesh generators) as

$$f_{\text{St}}(\bar{\mathcal{E}}, \tau; l)(x) = k(\bar{\mathcal{E}}; l)f_{\text{PD}}(\tau; l)(\tilde{x}(\bar{\mathcal{E}}, x; l)) \quad (14)$$

and

$$f_{\text{St}}(\bar{\mathcal{E}}; \alpha)(x) = k(\bar{\mathcal{E}}; \alpha)f_{\text{PD}}(\alpha)(\tilde{x}(\bar{\mathcal{E}}, x; \alpha)) \quad (15)$$

where $f_{\text{PD}}(\tau; l)$ and $f_{\text{PD}}(\alpha)$ are prolonged to be defined over the whole set of real numbers:

$$\forall x < 0, \quad f_{\text{PD}}(\tau; l)(x) = 0 \quad \text{and} \quad \forall x \in ]-\infty; 0[ \cup ]\pi; +\infty[, \quad f_{\text{PD}}(\alpha)(x) = 0 \quad (16)$$

with

$$\tilde{x}(\bar{\mathcal{E}}, x; l) = k(\bar{\mathcal{E}}; l)x + p(\bar{\mathcal{E}}; l) \\ \tilde{x}(\bar{\mathcal{E}}, x; \alpha) = k(\bar{\mathcal{E}}; \alpha)x + p(\bar{\mathcal{E}}; \alpha), \quad (17)$$

the contraction parameters being defined as

$$k(\bar{\mathcal{E}}; l) = \frac{\sigma_{\text{PD}}(\tau; l)}{\sigma_{\text{St}}(\bar{\mathcal{E}}, \tau; l)} = \frac{(1 + \bar{\mathcal{E}})^{1/2}}{1 - \bar{\mathcal{E}}} \quad (18)$$

and

$$k(\bar{\mathcal{E}}, \alpha) = \frac{\sigma_{\text{PD}}(\alpha)}{\sigma_{\text{St}}(\bar{\mathcal{E}}, \alpha)} = \frac{(1 + \bar{\mathcal{E}})^{1/4}}{(1 - \bar{\mathcal{E}})^{3/4}}, \quad (19)$$

and the shift parameters by $p(\bar{\mathcal{E}}, \tau; l) = (1-k(\bar{\mathcal{E}}; l))\mu_{\text{PD}}(\tau; l)$ and $p(\bar{\mathcal{E}}; \alpha) = (1-k(\bar{\mathcal{E}}; \alpha))\mu_{\text{PD}}(\alpha) = \frac{\pi}{3}(1-k(\bar{\mathcal{E}}; \alpha))$ where $\bar{\mathcal{E}}$ stands again for the reduced efficiency index given in Equation (11).

The accuracy of estimates (14) and (15) is shown in Figures 2 and 5 to 7. Since we have interest in standard mesh generators, four new cases are considered in Figures 6 and 7 covering various standard efficiency indexes and intensities (efficiency indexes associated to each mesh generator have been chosen arbitrarily): 1/ a mesh generated with the frontal algorithm of Gmsh with low intensity but high efficiency index, 2/ a mesh generated with Abaqus (2008) with very high intensity and high efficiency index, 3/ a mesh obtained with the pre-processor of the Cast3m (2011) finite element code with high intensity and intermediate efficiency index, 4/ a mesh drawn with the mesh generator Triangle (Shewchuk 1996) with high intensity and low efficiency index.

Figures 5 to 7 show that estimates (14) and (15) lead to accurate estimates of the edge length and triangle top angle means. However variances estimates are accurate only for tessellations whose efficiency index is not too high (typically equal or at least less than 0.92, see Figures 6 and 7). As a conclusion, estimates (14) and (15) are accurate for standard random meshes with efficiency index between 0.725 and 0.92 (i.e. with reduced efficiency index at least less than 0.71).
and

$$p(\bar{\mathcal{E}}; \alpha) = (1 - k(\bar{\mathcal{E}}; \alpha))\mu_{\text{PD}}(\alpha) = (1 - k(\bar{\mathcal{E}}; \alpha))\frac{\pi}{3}, \qquad (21)$$

where $\bar{\mathcal{E}}$ states again for the reduced efficiency index given in ??.

Accuracy of estimates (??) and (??) is shown in ??, ??, ?? and ??.

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Histograms of the edge length and triangle top angle sets obtained with the Gmsh mesh generator (Geuzaine et al. 2009), with “MeshAdapt” strategy. First line: edge lengths distributions. Second line: triangle top angle distributions. From left to right the meshes with increasing intensity $\tau$, (coarser mesh $\tau = 525$, intermediate mesh $\tau = 831$, finest mesh $\tau = 1611$) black lines refer to the standard mesh estimators for the density functions of both the edge length ?? and the triangle top angle ?? sets. The efficiency index of the finest mesh is about $\mathcal{E} \approx 0.86$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Edge length probability density: accuracy of estimate (??) (black line). Meshes obtained with: 1/ the “frontal” algorithm of Gmsh ($\tau = 750, \mathcal{E} = 0.976$), 2/ the Abaqus software ($\tau = 1386, \mathcal{E} = 0.981$), 3/ the Cast3m software ($\tau = 1102, \mathcal{E} = 0.917$), 4/ the Triangle mesh generator ($\tau = 1008, \mathcal{E} = 0.863$).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5** Histograms of the edge length and triangle top angle sets obtained with the Gmsh mesh generator (Geuzaine et al. 2009), with “MeshAdapt” strategy. First line: edge lengths distributions. Second line: triangle top angle distributions. From left to right the meshes with increasing intensity $\tau$, (coarser mesh $\tau = 525$, intermediate mesh $\tau = 831$, finest mesh $\tau = 1611$) black lines refer to the standard mesh estimators for the density functions of both the edge length (14) and the triangle top angle (15) sets. The efficiency index of the finest mesh is about $\mathcal{E} \approx 0.86$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6** Edge length probability density: accuracy of estimate (14) (black line). Meshes obtained with: 1/ the “frontal” algorithm of Gmsh ($\tau = 750$, $\mathcal{E} = 0.976$), 2/ the Abaqus software ($\tau = 1386$, $\mathcal{E} = 0.981$), 3/ the Cast3m software ($\tau = 1102$, $\mathcal{E} = 0.917$), 4/ the Triangle mesh generator ($\tau = 1008$, $\mathcal{E} = 0.863$).</figcaption></figure>
Since we have interest for standard mesh generators, four new cases are considered in ?? and ?? covering various standard efficiency indexes and intensities (efficiency indexes associated to each mesh generator have been chosen arbitrarily):

1/ a mesh generated with the *frontal* algorithm of Gmsh with low intensity but high efficiency index, 2/ a mesh generated with Abaqus (2008) with very high intensity and high efficiency index, 3/ a mesh obtained with the pre-processor of the Cast3m (2011) finite element code with high intensity and intermediate efficiency index, 4/ a mesh drawn with the mesh generator Triangle (Shewchuk 2005) with high intensity and low efficiency index.

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Top angle probability density: accuracy of estimate (??) (black line). Meshes obtained with: 1/ the “frontal” algorithm of Gmsh ($\tau = 750, \mathcal{E} = 0.976$), 2/ the Abaqus software ($\tau = 1386, \mathcal{E} = 0.981$), 3/ the Cast3m software ($\tau = 1102, \mathcal{E} = 0.917$), 4/ the Triangle mesh generator ($\tau = 1008, \mathcal{E} = 0.863$).</figcaption></figure>

??, ?? and ?? show that estimates (??) and (??) lead to accurate estimates of the edge length and triangle top angle means. However variances estimates are accurate only for tessellations whose efficiency index is not too high (typically equal or at least less than 0.92, see ?? and ??). As a conclusion, estimates (??) and (??) are accurate for standard random meshes with efficiency index between 0.725 and 0.92 (i.e. with reduced efficiency index at least less than 0.71).

# 5 Density of random triangular meshes

The mesh density Z is defined as the sum of the edge lengths $\sum_{i=1}^{N_E} l_i$ divided by the area A of the meshes domain $\Omega$:

$$Z = \frac{1}{A} \sum_{i=1}^{N_E} l_i$$

where $N_E$ is the number of edges of the mesh.

For practical purpose, we showed that the pdf of the length l of the typical edge is given by ???. Due to the central limit theorem, the sum $\sum_{i=1}^{N_E} l_i$ thus follows a normal distribution with mean $N_E \times \mu_{\text{St}}(\tau; l)$ and variance $N_E \times \sigma_{\text{St}}^2(\tilde{\mathcal{E}}, \tau; l)$ when the number of edges $N_E$ is large enough (Billingsley 1995). Then the mesh density Z follows a normal distribution with mean $(N_E/A) \times \mu_{\text{St}}(\tau; l)$ and variance $(N_E/A^2) \times \sigma_{\text{St}}^2(\tilde{\mathcal{E}}, \tau; l)$.

Recalling that $N_E = 3\tau A$, one obtains the probability density function $f_{\text{St}}(\tilde{\mathcal{E}}, \tau; Z)$ of the mesh density as a normal distribution with mean $\mu_{\text{St}}(\tau; Z)$ equal to $3\tau\mu_{\text{St}}(\tau; l)$ and variance $\sigma_{\text{St}}^2(\tilde{\mathcal{E}}, A; Z)$ equal to $(3\tau/A)\sigma_{\text{St}}^2(\tilde{\mathcal{E}}, \tau; l)$:

$$f_{\text{St}}(\tilde{\mathcal{E}}, \tau; Z) = \mathcal{P}\mathcal{D}\mathcal{F} \left[ N \left( \mu_{\text{St}}(\tau; Z), \sigma_{\text{St}}^2(\tilde{\mathcal{E}}, A; Z) \right) \right] \quad (22)$$
<figure><img src="image_8.png" /><figcaption>Figure 7 Top angle probability density: accuracy of estimate (15) (black line). Meshes obtained with: 1/ the “frontal” algorithm of Gmsh (τ = 750, E = 0.976), 2/ the Abaqus software (τ = 1386, E = 0.981), 3/ the Cast3m software (τ = 1102, E = 0.917), 4/ the Triangle mesh generator (τ = 1008, E = 0.863).</figcaption></figure>

5 Density of random triangular meshes

The mesh density

$$Z = \frac{1}{A} \sum_{i=1}^{N_E} l_i \qquad (20)$$

is defined as the sum of the edge lengths divided by the area of the meshes domain Ω, where $N_E$ is the number of edges of the mesh.

For practical purpose, we showed that the pdf of length $l$ of the typical edge is given by Equation (14). Due to the central limit theorem, the sum $\sum_{i=1}^{N_E} l_i$ follows a normal distribution with mean $N_E \times \mu_{\text{St}}(\tau; l)$ and variance $N_E \times \sigma_{\text{St}}^2(\bar{\mathcal{E}}, \tau; l)$ when the number of edges $N_E$ is large enough (Billingsley 1995). Then the mesh density $Z$ follows a normal distribution with mean $(N_E/A) \times \mu_{\text{St}}(\tau; l)$ and variance $(N_E/A^2) \times \sigma_{\text{St}}^2(\bar{\mathcal{E}}, \tau; l)$.

Recalling that $N_E = 3\tau A$, one obtains the probability density function $f_{\text{St}}(\bar{\mathcal{E}}, \tau; Z)$ of the mesh density as a normal distribution with mean $\mu_{\text{St}}(\tau; Z)$ equal to $3\tau\mu_{\text{St}}(\tau; l)$ and variance $\sigma_{\text{St}}^2(\bar{\mathcal{E}}, A; Z)$ equal to $(3\tau/A)\sigma_{\text{St}}^2(\bar{\mathcal{E}}, \tau; l)$:

$$f_{\text{St}}(\bar{\mathcal{E}}, \tau; Z) = \mathcal{P}\mathcal{D}\mathcal{F}[N(\mu_{\text{St}}(\tau; Z), \sigma_{\text{St}}^2(\bar{\mathcal{E}}, A; Z))] \qquad (21)$$

with

$$\mu_{\text{St}}(\tau; Z) = \frac{32\sqrt{\tau}}{3\pi} \quad (\text{mean}) \quad \text{and} \quad \sigma_{\text{St}}^2(\bar{\mathcal{E}}, A; Z) \approx \frac{(1 - \bar{\mathcal{E}})^2}{1 + \bar{\mathcal{E}}} \left( \frac{405\pi - 32^2}{27\pi^2 A} \right) \quad (\text{variance}). \qquad (22)$$

It is worth noting that the variance of this mesh density no more depends on the intensity $\tau$, but now depends on the area $A$ of the meshed domain. In particular, $\sigma_Z^2$ depends on the unit length used to define the boundary of the domain $\partial\Omega$.

The ratio in percentage between the standard deviation and the mean of a standard mesh density is given by

$$\frac{\sigma_{\text{St}}(\bar{\mathcal{E}}, \tau; Z)}{\mu_{\text{St}}(\tau; Z)} \times 100 \approx \left(0.49 \frac{1 - \bar{\mathcal{E}}}{\sqrt{1 + \bar{\mathcal{E}}}}\right) \frac{1}{\sqrt{N_E}} \times 100. \qquad (23)$$
with

$$ \mu_{\text{St}}(\tau; Z) = \frac{32\sqrt{\tau}}{3\pi} \text{ (mean)} \quad \text{and} \quad \sigma_{\text{St}}^2(\bar{\mathcal{E}}, A; Z) \approx \frac{(1 - \bar{\mathcal{E}})^2}{1 + \bar{\mathcal{E}}} \left( \frac{405\pi - 32^2}{27\pi^2 A} \right) \text{ (variance).} \quad (23) $$

It is worth noting that the variance of this mesh density no more depends on the intensity $\tau$, but now depends on the area A of the meshed domain. In particular, $\sigma_Z^2$ depends on the unit length used to define the boundary of the domain $\partial\Omega$.

The ratio in percentage between the standard deviation and the mean of a standard mesh density is given by:

$$ \frac{\sigma_{\text{St}}(\bar{\mathcal{E}}, \tau; Z)}{\mu_{\text{St}}(\tau; Z)} \times 100 \approx \left( 0.49 \frac{1 - \bar{\mathcal{E}}}{\sqrt{1 + \bar{\mathcal{E}}}} \right) \frac{1}{\sqrt{N_E}} \times 100. \quad (24) $$

This relation leads to the following result:

$$ \forall \mathcal{E} \ge 0.86, \quad \frac{\sigma_{\text{St}}(\bar{\mathcal{E}}, \tau; Z)}{\mu_{\text{St}}(\tau; Z)} \times 100 \le \frac{20}{\sqrt{N_E}}. \quad (25) $$

In fracture mechanics, $N_E$ is sufficiently large to consider that the variance of Z is negligible. We then consider for practical purpose that Z can be understand as a scalar value depending exclusively on the intensity of the underlying point-process or equivalently on the edge length mean:

$$ Z_{\text{St}}(\tau) = \frac{32\sqrt{\tau}}{3\pi} = \frac{1024}{27\pi^2\mu_{\text{St}}(\tau; l)}. \quad (26) $$

The accuracy of estimate (??) is illustrated in ??.

<figure><img src="image_9.png" /><figcaption>**Figure 8:** The mesh density Z as a function of the intensity of the point-process $\tau$ (line: $\mu_{\text{St}}(\tau; Z)$ in ??, unitary area). Gmsh results with “MeshAdapt” strategy (points) and with “Frontal” strategy (open circles); Cast3m result (cross); Triangle result (filled square); Abaqus result (open square). Corresponding meshes are given in ?? and ??.</figcaption></figure>

The dependence of the mesh density on the square root of the intensity $\tau$ is in accordance with the Vaschy-Buckingham theorem. Indeed, the intensity of the point process and the total area of the discretized domain fully describe the tessellation geometrical features statistics. Thus, an intensive feature such as the mesh density Z would only depend on the point process intensity $\tau$. Using the Vaschy-Buckingham theorem, it is then expected that the density Z - homogeneous to the inverse of a length - would be proportional to the square root of the point process intensity.

# 6 Geodesic tortuosity of random Delaunay meshes

## 6.1 Shortest path approximation

Let A and B be two nodes of a planar triangular mesh. Among all the paths which are composed of a collection of adjacent edges and whose extremities are nodes A and B, there exists at least one path with a minimum length. Let an edge inclination be the angle between this edge and direction AB, whose value has to be taken between $-\frac{\pi}{2}$ and $\frac{\pi}{2}$.
This relation leads to the result

$$ \forall \varepsilon \ge 0.86, \quad \frac{\sigma_{\text{St}}(\tilde{\varepsilon}, \tau; Z)}{\mu_{\text{St}}(\tau; Z)} \times 100 \le \frac{20}{\sqrt{N_E}}. \qquad (24) $$

In fracture mechanics, $N_E$ is sufficiently large to consider that the variance of Z is negligible. We then consider for practical purpose that Z can be understand as a scalar value depending exclusively on the intensity of the underlying point-process or equivalently on the edge length mean:

$$ Z_{\text{St}}(\tau) = \frac{32\sqrt{\tau}}{3\pi} = \frac{1024}{27\pi^2\mu_{\text{St}}(\tau; l)}. \qquad (25) $$

The accuracy of estimate (25) is illustrated in Figure 8. The dependence of the mesh density on

<figure><img src="image_9.png" /><figcaption>**Figure 8** Mesh density Z as a function of the intensity of the point-process $\tau$ (line: $\mu_{\text{St}}(\tau; Z)$ in Equation (21), unitary area). Gmsh results with the “MeshAdapt” strategy (points) and with the “Frontal” strategy (open circles); Cast3m result (cross); Triangle result (filled square); Abaqus result (open square). Corresponding meshes in Figures 6 and 7.</figcaption></figure>

the square root of the intensity $\tau$ is in accordance with the Vaschy-Buckingham theorem. Indeed, the intensity of the point process and the total area of the discretized domain fully describe the tessellation geometrical features statistics. Thus, an intensive feature such as the mesh density Z would only depend on the point process intensity $\tau$. Using the Vaschy-Buckingham theorem, it is then expected that the density Z—homogeneous to the inverse of a length—would be proportional to the square root of the point process intensity.

# 6 Geodesic tortuosity of random Delaunay meshes

## 6.1 Shortest path approximation

Let A and B be two nodes of a planar triangular mesh. Among all the paths which are composed of a collection of adjacent edges and whose extremities are nodes A and B, there exists at least one path with a minimum length. Let an edge inclination be the angle between this edge and direction AB, whose value has to be taken between $-\frac{\pi}{2}$ and $\frac{\pi}{2}$.

The orientation with respect to a given direction AB of the edges of a PD mesh is uniformly distributed over $[-\frac{\pi}{2}; \frac{\pi}{2}]$. It is reasonably supposed that this observation remains valid for meshes obtained from a PD tessellation after few Lloyd relaxations (as long as the mesh remains sufficiently unstructured). In the sequel, such meshes will be referred as “unstructured Delaunay-like meshes”. As a consequence, in such meshes only the distance $d(A, B)$ has to be considered for estimating the shortest path length between A and B. Considering a fixed value for the distance $d(A, B)$ is somehow the same as fixing the number of edges n composing the shortest path. For sake of simplicity, the choice is therefore made to characterise a shortest path only by the number n of edges it contains.

At that point, the edge length and triangle top angle distributions are the only geometrical data on which we can rely on for estimating the shortest path length between two nodes of a given planar triangular mesh. However, to the best of our knowledge, there is no explicit mathematical expression relying these quantities. The difficulty lies in the fact that a shortest path is a “global” quantity in the way that there is no explicit rule followed individually by each edge. Thus all the edges of the shortest path have to be found at the same time and cannot be chosen independently from each other.
The orientation with respect to a given direction AB of the edges of a PD mesh is uniformly distributed over [-π/2; π/2]. It is reasonably supposed that this observation remains valid for meshes obtained from a PD tessellation after few Lloyd relaxations (as long as the mesh remains suffi- ciently unstructured). In the sequel, such meshes will be referred as “unstructured Delaunay-like meshes”. As a consequence, in such meshes only the distance d(A, B) has to be considered for estimating the shortest path length between A and B. Considering a fixed value for the distance d(A, B) is somehow the same as fixing the number of edges n composing the shortest path. For sake of simplicity, the choice is therefore made to characterise a shortest path only by the number n of edges it contains.

At that point, the edge length and triangle top angle distributions are the only geometrical data on which we can rely on for estimating the shortest path length between two nodes of a given planar triangular mesh. However, to the best of our knowledge, there is no explicit mathematical expression relying these quantities. The difficulty lies in the fact that a shortest path is a “global” quantity in the way that there is no explicit rule followed individually by each edge. Thus all the edges of the shortest path have to be found at the same time and cannot be chosen independently from each other.

Facing this difficult problem, it has been considered to focus on *n*-edges paths verifying the following property, hereafter refered as “inclination rule”:

$$ Each edge inclination of the path is as close as possible to zero. \qquad (27) $$

The idea is to consider paths as close as possible to the shortest one and whose description is made easier by the inclination rule they fulfill. Thanks to this rule the probability density function of the tortuosity of such a path - i.e. the ratio between the path length and the length of its projection along the direction AB - can be estimated (see ??).

?? represents on its first column examples of a shortest path between two nodes A and B of a PD tessellation and a Lloyd relaxed PD tessellation. A distinction is made between edges following the inclination rule (blue ones) and those who are not (red ones). On the same figure, on the second column, corresponding histograms of the percentage of the path edges which are not following the inclination rule are displayed.

Although all the edges are not following the inclination rule, a major part of them does. To consider that the shortest path is following the rule (??) is then not optimal but seems acceptable for estimating the shortest path tortuosity.

## 6.2 Definition of the tortuosity of a path following the inclination rule

Among all edges connected to a given node, at least one is associated to a minimum inclination (along direction AB) absolute value. The inclination θ of this particular edge is referred in the sequel as “best inclination” (see ??).

Edges best inclinations of a path following the inclination rule (??) are particular events of the random variable θ. For sake of simplicity, we consider in the sequel that these events are statistically independent. The tortuosity of a *n*-edges path of a given tessellation G following the inclination rule is given by:

$$ T_G(n) = \frac{\sum_{i=1}^{n} l_i}{\sum_{i=1}^{n} l_i \cos(\theta_i)}, \qquad (28) $$

where $l_i$ are the edge lengths of the path and $\theta_i$ its edge inclinations.

In the particular case where the path is a straight line, $\forall i \in [1; n]$, $\theta_i = 0$, the tortuosity is equal to 1. For other case the tortuosity is greater than 1.
Facing this difficult problem, it has been considered to focus on *n*-edges paths verifying the following property, hereafter referred to as “inclination rule”:

$$ \textit{Each edge inclination of the path is as close as possible to zero.} \qquad (26) $$

The idea is to consider paths as close as possible to the shortest one and whose description is made easier by the inclination rule they fulfill. Thanks to this rule the probability density function of the tortuosity of such a path—i.e. the ratio between the path length and the length of its projection along the direction *AB*—can be estimated (see Section 6.4).

Figure 9 represents on its first column examples of a shortest path between two nodes *A* and *B* of a PD tessellation and a Lloyd relaxed PD tessellation. A distinction is made between edges following the inclination rule (blue ones) and those who are not (red ones). On the same figure, on the second column, corresponding histograms of the percentage of the path edges which are not following the inclination rule are displayed.

<figure><img src="image_10.png" /><figcaption>**Figure 9** Relevance of the use of the inclination rule. First column: Representation of the shortest path between two nodes *A* and *B* in a PD tessellation (top) and a Lloyd relaxed tessellation (bottom) of a Poisson point process ($\tau = 300$). Blue edges follow the inclination rule, red ones do not. On the second column: over 1000 PD tessellations (top) and over 1000 Lloyd relaxed tessellations (bottom) of a Poisson point process ($\tau = 300$), histograms of the percentage of edges of the shortest path which are not following the inclination rule.</figcaption></figure>

Although all the edges are not following the inclination rule, a major part of them does. To consider that the shortest path is following the rule (26) is then not optimal but seems acceptable for estimating the shortest path tortuosity.

## 6.2 Definition of the tortuosity of a path following the inclination rule

Among all edges connected to a given node, at least one is associated to a minimum inclination (along direction *AB*) absolute value. The inclination $\theta$ of this particular edge is referred to in the sequel as “best inclination”, see Figure 10.

Edges best inclinations of a path following the inclination rule (26) are particular events of the random variable $\theta$. For sake of simplicity, we consider in the sequel that these events are statistically independent. The tortuosity of a *n*-edges path of a given tessellation *G* following the inclination rule is given by

$$ T_G(n) = \left( \sum_{i=1}^{n} l_i \cos \theta_i \right)^{-1} \sum_{i=1}^{n} l_i, \qquad (27) $$

where $l_i$ are the edge lengths of the path and $\theta_i$ its edge inclinations.
<figure><img src="image_10.png" /><figcaption>**Figure 9:** First column: Representation of the shortest path between two nodes A and B in a PD tessellation (top) and a Lloyd relaxed tessellation (bottom) of a Poisson point process ($\tau = 300$). Blue edges follow the inclination rule, red ones do not. On the second column: over 1000 PD tessellations (top) and over 1000 Lloyd relaxed tessellations (bottom) of a Poisson point process ($\tau = 300$), histograms of the percentage of edges of the shortest path which are not following the inclination rule.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10:** Schematic definition of the incident angle $\gamma$, bisector orientation $\beta$ and best inclination $\theta$.</figcaption></figure>

## 6.3 Probability density function of the best inclination in Delaunay-type meshes

Except for trivial cases, the half straight line from a node $P$ with direction $AB$ crosses exactly one of the triangle elements surrounding $P$. We respectively denote as $\beta$ and $\gamma$ the bisector orientation with respect to the direction $AB$ and the angle of this triangle at top $P$, hereafter called the incident angle (see ??).

Node $P$ is surrounded by triangles having their own angle at summit $P$. The orientation with respect to a given direction $AB$ of the edges of an unstructured Delaunay-like mesh may reasonably supposed uniformly distributed over $[-\frac{\pi}{2}; \frac{\pi}{2}]$. As a consequence, the probability that a summit angle corresponds to a particular value of $\gamma$ is proportional to the value of this angle. In other words, it is more likely that among the surrounding triangles at node $P$, the triangle crossed by the half straight line from node $P$ with direction $AB$ be the one having the highest angle at summit $P$. Considering then that the angles at summit $P$ are independent events of the
In the particular case where the path is a straight line, $\forall i \in [1; n]$, $\theta_i = 0$, the tortuosity is equal to 1. For the other cases, the tortuosity is greater than 1.

<figure><img src="image_11.png" /><figcaption>**Figure 10** Definition of the incident angle $\gamma$, bisector orientation $\beta$ and best inclination $\theta$.</figcaption></figure>

## 6.3 Probability density function of the best inclination in Delaunay-type meshes

Except for trivial cases, the half straight line from a node P with direction AB crosses exactly one of the triangle elements surrounding P. We respectively denote as $\beta$ and $\gamma$ the bisector orientation with respect to the direction AB and the angle of this triangle at top P, hereafter called the incident angle (see Figure 10).

Node P is surrounded by triangles having their own angle at summit P. The orientation with respect to a given direction AB of the edges of an unstructured Delaunay-like mesh may reasonably supposed uniformly distributed over $[-\frac{\pi}{2}; \frac{\pi}{2}]$. As a consequence, the probability that a summit angle corresponds to a particular value of $\gamma$ is proportional to the value of this angle. In other words, it is more likely that among the surrounding triangles at node P, the triangle crossed by the half straight line from node P with direction AB be the one having the highest angle at summit P. Considering then that the angles at summit P are independent events of the random variable $\alpha$ (as defined in Section 3), the following expression for the pdf of $\gamma$ is obtained:

$$f_{\text{St}}(\bar{\mathcal{E}}; \gamma)(x) = \frac{k(\bar{\mathcal{E}}; \alpha)}{\frac{\pi}{3} - p(\bar{\mathcal{E}}; \alpha)} \times x f_{\text{St}}(\bar{\mathcal{E}}; \alpha)(x), \quad (28)$$

where the prefactor ensures that $\int_{0}^{\pi} f_{\text{St}}(\bar{\mathcal{E}}; \gamma)(x) dx = 1$. By construction, the bisector orientation verifies

$$\beta \in [-\gamma/2; \gamma/2]. \quad (29)$$

Since the edges inclinations follow a uniform distribution, the conditional density of $\beta$ knowing $\gamma = \gamma_0$ is derived from Equation (29):

$$g_{\text{St}}(\bar{\mathcal{E}}; \beta|\gamma = \gamma_0)(x) = \begin{cases} \frac{1}{\gamma_0} & \text{if } x \in [-\frac{\gamma_0}{2}, \frac{\gamma_0}{2}] \\ 0 & \text{if not} \end{cases} = \begin{cases} \frac{1}{\gamma_0} & \text{if } |x| \le \frac{\pi}{2} \text{ and } \gamma_0 \in [2|x|, \pi] \\ 0 & \text{if not.} \end{cases} \quad (30)$$

Using then the law of total probability, the pdf of $\beta$ reads

$$f_{\text{St}}(\bar{\mathcal{E}}; \beta)(x) = \int_{0}^{\pi} g_{\text{St}}(\bar{\mathcal{E}}; \beta|\gamma = z)(x) f_{\text{St}}(\bar{\mathcal{E}}; \gamma)(z) dz \quad (31)$$

which, according to Equations (28) and (30) could be written as

$$\begin{aligned} f_{\text{St}}(\bar{\mathcal{E}}; \beta)(x) &= \int_{2|x|}^{\pi} \frac{f_{\text{St}}(\bar{\mathcal{E}}; \gamma)(z)}{z} dz = \left(\frac{k(\bar{\mathcal{E}}; \alpha)}{\frac{\pi}{3} - p(\bar{\mathcal{E}}; \alpha)}\right) \int_{2|x|}^{\pi} f_{\text{St}}(\bar{\mathcal{E}}; \alpha)(z) dz \\ &= \left(\frac{k(\bar{\mathcal{E}}; \alpha)}{\frac{\pi}{3} - p(\bar{\mathcal{E}}; \alpha)}\right) (1 - F_{\text{St}}(\bar{\mathcal{E}}; \alpha)(2|x|)), \end{aligned} \quad (32)$$
random variable $\alpha$ (as defined in ??), the following expression for the pdf of $\gamma$ is obtained:

$$f_{\text{St}}(\bar{\mathcal{E}}; \gamma)(x) = \frac{k(\bar{\mathcal{E}}; \alpha)}{\frac{\pi}{3} - p(\bar{\mathcal{E}}; \alpha)} \times x f_{\text{St}}(\bar{\mathcal{E}}; \alpha)(x), \quad (29)$$

where the prefactor ensures that $\int_0^\pi f_{\text{St}}(\bar{\mathcal{E}}; \gamma)(x) \, dx = 1$.

By construction, the bisector orientation $\beta$ verifies:

$$\beta \in [-\frac{\gamma}{2}, \frac{\gamma}{2}]. \quad (30)$$

Since the edges inclinations follow a uniform distribution, the conditional density of $\beta$ knowing $\gamma = \gamma_0$ is derived from ??:

$$\begin{aligned} g_{\text{St}}(\bar{\mathcal{E}}; \beta|\gamma = \gamma_0)(x) &= \begin{cases} \frac{1}{\gamma_0} & \text{if } x \in [-\frac{\gamma_0}{2}, \frac{\gamma_0}{2}] \\ 0 & \text{if not} \end{cases} \\ &= \begin{cases} \frac{1}{\gamma_0} & \text{if } |x| \le \frac{\pi}{2} \text{ and } \gamma_0 \in [2|x|, \pi], \\ 0 & \text{if not} \end{cases} \end{aligned} \quad (31)$$

Using then the law of total probability, the pdf of $\beta$ reads:

$$f_{\text{St}}(\bar{\mathcal{E}}; \beta)(x) = \int_{0}^{\pi} g_{\text{St}}(\bar{\mathcal{E}}; \beta|\gamma = z)(x) f_{\text{St}}(\bar{\mathcal{E}}; \gamma)(z) dz, \quad (32)$$

which, according to ?? and ??, could be written as follow:

$$\begin{aligned} f_{\text{St}}(\bar{\mathcal{E}}; \beta)(x) &= \int_{2|x|}^{\pi} \frac{f_{\text{St}}(\bar{\mathcal{E}}; \gamma)(z)}{z} dz = \left( \frac{k(\bar{\mathcal{E}}; \alpha)}{\frac{\pi}{3} - p(\bar{\mathcal{E}}; \alpha)} \right) \int_{2|x|}^{\pi} f_{\text{St}}(\bar{\mathcal{E}}; \alpha)(z) dz \\ &= \left( \frac{k(\bar{\mathcal{E}}; \alpha)}{\frac{\pi}{3} - p(\bar{\mathcal{E}}; \alpha)} \right) (1 - F_{\text{St}}(\bar{\mathcal{E}}; \alpha)(2|x|)), \end{aligned} \quad (33)$$

where $F_{\text{St}}(\bar{\mathcal{E}}; \alpha)$ is the cumulative density function of $\alpha$ given by:

$$F_{\text{St}}(\bar{\mathcal{E}}; \alpha)(x) = \frac{1}{3\pi} \left[ \pi + 2\tilde{x}(\bar{\mathcal{E}}, x; \alpha) + (\tilde{x}(\bar{\mathcal{E}}, x; \alpha) - \pi) \cos(2\tilde{x}(\bar{\mathcal{E}}, x; \alpha)) - \frac{3}{2} \sin(2\tilde{x}(\bar{\mathcal{E}}, x; \alpha)) \right], \quad \text{if } \tilde{x}(\bar{\mathcal{E}}, x; \alpha) \in [0, \pi], \quad (34)$$

if $\tilde{x}(\bar{\mathcal{E}}, x; \alpha) \in [0, \pi]$, and by

$$F_{\text{St}}(\bar{\mathcal{E}}; \alpha)(x) = 0, \quad \text{if } \tilde{x}(\bar{\mathcal{E}}, x; \alpha) \le 0 \quad (35)$$

and

$$F_{\text{St}}(\bar{\mathcal{E}}; \alpha)(x) = 1, \quad \text{if } \tilde{x}(\bar{\mathcal{E}}, x; \alpha) \ge \pi, \quad (36)$$

where $\tilde{x}(\bar{\mathcal{E}}, x; \alpha)$ is defined in ??.

Considering ?? and denoting by $i$ the $i^{th}$ edge of a path following the inclination rule (??), a relation can be exhibited between $\theta_i$, $\beta_i$ and $\gamma_i$:

$$\text{if } \beta_i \ge 0, \quad \theta_i = \beta_i - \frac{\gamma_i}{2} \le 0, \quad \text{else, } \theta_i = \beta_i + \frac{\gamma_i}{2} \ge 0, \quad (37)$$

which can be rewritten as follows:

$$\text{if } \theta_i \ge 0, \quad \beta_i = \theta_i - \frac{\gamma_i}{2} \le 0, \quad \text{else, } \beta_i = \theta_i + \frac{\gamma_i}{2} \ge 0. \quad (38)$$
where $F_{\text{St}}(\bar{\mathcal{E}}; \alpha)$ is the cumulative density function of $\alpha$ given by

$$F_{\text{St}}(\bar{\mathcal{E}}; \alpha)(x) = \begin{cases} 0 & \text{if } \tilde{x}(\bar{\mathcal{E}}, x; \alpha) \le 0 \\ \frac{1}{3\pi} \left[ \pi + 2\tilde{x}(\bar{\mathcal{E}}, x; \alpha) + (\tilde{x}(\bar{\mathcal{E}}, x; \alpha) - \pi) \cos(2\tilde{x}(\bar{\mathcal{E}}, x; \alpha)) - \frac{3}{2} \sin(2\tilde{x}(\bar{\mathcal{E}}, x; \alpha)) \right] & \text{if } \tilde{x}(\bar{\mathcal{E}}, x; \alpha) \in [0, \pi] \\ 1 & \text{if } \tilde{x}(\bar{\mathcal{E}}, x; \alpha) \ge \pi \end{cases} \quad (33)$$

where $\tilde{x}(\bar{\mathcal{E}}, x; \alpha)$ is defined in Equation (17).

Considering Figure 10 and denoting by $i$ the $i^{th}$ edge of a path following the inclination rule (26), a relation can be exhibited between $\theta_i$, $\beta_i$ and $\gamma_i$:

$$\text{if } \beta_i \ge 0, \quad \theta_i = \beta_i - \gamma_i/2 \le 0, \quad \text{else } \theta_i = \beta_i + \gamma_i/2 \ge 0, \quad (34)$$

which can be rewritten as

$$\text{if } \theta_i \ge 0, \quad \beta_i = \theta_i - \gamma_i/2 \le 0, \quad \text{else } \beta_i = \theta_i + \gamma_i/2 \ge 0. \quad (35)$$

The symmetry of relations (34) and (35) implies that $\beta$ and $\theta$ have the same density

$$f_{\text{St}}(\bar{\mathcal{E}}; \theta)(x) = \left( \frac{k(\bar{\mathcal{E}}; \alpha)}{\frac{\pi}{3} - p(\bar{\mathcal{E}}; \alpha)} \right) (1 - F_{\text{St}}(\bar{\mathcal{E}}; \alpha)(2|x|)). \quad (36)$$

Figure 11 compares the distributions of $\gamma$, $\beta$ and $\theta$ for a PD tessellation and a Lloyd relaxed PD tessellation. A very good agreement is obtained between the actual distributions and the different pdf expressions defined in Equations (28), (32) and (36).

<figure><img src="image_12.png" /><figcaption>Figure 11 Accuracy of the estimated probability density functions of angles $\gamma$, $\beta$ and $\theta$. First line: angles distributions of a PD tessellation ($\tau = 500$). Second line: angles distributions of a Lloyd relaxed PD tessellation ($\mathcal{E} = 0.89$, $\tau = 500$). From left to right: $\gamma$, $\beta$ and $\theta$ distributions. Black lines: density functions expressed in Equations (28), (32) and (36).</figcaption></figure>

## 6.4 Tortuosity probability density function of a path following the inclination rule

Equation (27) can be rewritten as

$$T_{\text{St}}(n) = \left( \sum_{i=1}^{n} w_i \cos \theta_i \right)^{-1} \quad (37)$$
The symmetry of relations (??) and (??) implies that $\beta$ and $\theta$ have the same density:

$$f_{\text{St}}(\bar{\mathcal{E}}; \theta)(x) = \left( \frac{k(\bar{\mathcal{E}}; \alpha)}{\frac{\pi}{3} - p(\bar{\mathcal{E}}; \alpha)} \right) (1 - F_{\text{St}}(\bar{\mathcal{E}}; \alpha)(2|x|)). \quad (39)$$

?? compares the distributions of $\gamma$, $\beta$ and $\theta$ for a PD tessellation and a Lloyd relaxed PD tessellation. A very good agreement is obtained between the actual distributions and the different pdf expressions defined in ??, ?? and ??.

<figure><img src="image_12.png" /><figcaption>Figure 11: First line: angles distributions of a PD tessellation ($\tau = 500$). Second line: angles distributions of a Lloyd relaxed PD tessellation ($\mathcal{E} = 0.89$, $\tau = 500$). From left to right: $\gamma$, $\beta$ and $\theta$ distributions. Black lines: density functions ??, ?? and ??.</figcaption></figure>

## 6.4 Probability density function of the tortuosity of a path following the inclination rule

?? can be rewritten as follows:

$$T_{\text{St}}(n) = \frac{1}{\sum_{i=1}^{n} w_i \cos(\theta_i)}, \quad (40)$$

with weigth $w_i$ defined by the proportion of the $i^{th}$ edge length in the total length of the path:

$$l_i = w_i \sum_{k=1}^{n} l_k. \quad (41)$$

It is worth noting that with such a definition, summation of $w_i$ values is equal to 1.

?? represents $w_i$ values versus $\cos(\theta_i)$ values of a PD tessellation (left) and a Lloyd relaxed PD one (right). A correlation can be exhibited between $w_i$ and $\cos(\theta_i)$ since low values of $\cos(\theta_i)$ seem to be associated to low values of $w_i$. The following hypothesis is therefore considered in the sequel for purpose of simplification:

$$\textit{The expected value of } w_i \textit{ knowing that } \cos(\theta_i) = c_0 \textit{ is proportional to } c_0. \quad (42)$$
with weight $w_i$ defined by the proportion of the $i^{th}$ edge length in the total length of the path:

$$l_i = w_i \sum_{k=1}^{n} l_k. \qquad (38)$$

It is worth noting that with such a definition, summation of $w_i$ values is equal to 1.

Figure 12 represents $w_i$ values versus $\cos \theta_i$ values of a PD tessellation (left) and a Lloyd relaxed PD one (right). A correlation can be exhibited between $w_i$ and $\cos \theta_i$ since low values of $\cos \theta_i$ seem to be associated to low values of $w_i$. The following hypothesis is therefore considered in the sequel for purpose of simplification:

$$\text{The expected value of } w_i \text{ knowing that } \cos \theta_i = c_0 \text{ is proportional to } c_0. \qquad (39)$$

As a consequence, it is necessary to compute the correlation coefficient (which would differs

<figure><img src="image_13.png" /><figcaption>**Figure 12** Weights $w_i$ versus $\cos \theta_i$ of a PD tessellation (left) and a Lloyd relaxed PD one (right).</figcaption></figure>

from zero) between weights $w_i$ and $\cos \theta_i$ sets for the computation of the probabilistic density function of the tortuosity.

Splitting $I$, the interval of accessible values for $\cos \theta_i$, into $m$ intervals $I_k$, $k \in [1; m]$ of equal lengths, the sum at denominator of Equation (37) reads

$$\sum_{i=1}^{n} w_i \cos \theta_i = \sum_{k=1}^{m} \left( \sum_{i=1}^{n_k} w_i \cos \theta_i \right) \qquad (40)$$

where $n_k$ corresponds to the number of $\cos \theta$ realizations which are in $I_k$. Therefore, if $n_k$ is large enough and according to the central limit theorem, the sum can be approximated by

$$\sum_{i=1}^{n} w_i \cos \theta_i \approx \sum_{k=1}^{m} n_k \langle w \cos \theta \rangle_k = \sum_{k=1}^{m} \left( \sum_{i=1}^{n_k} \langle w \cos \theta \rangle_k \right) \qquad (41)$$

where notation $\langle X \rangle_k$ stands for the expected value of the random variable $X$ inside the interval $I_k$ and $\theta$ and $w$ are respectively the random variables associated to the values $\theta_i$ and $w_i$. If the interval $I_k$ is narrow enough, values $\cos \theta_i$ in $I_k$ can be considered as equal to a constant value corresponding to the center $c_k$ of interval $I_k$. Then $\langle w \cos \theta \rangle_k \approx \langle w \rangle_k c_k$ and estimate (41) becomes

$$\sum_{i=1}^{n} w_i \cos \theta_i \approx \sum_{k=1}^{m} \left( \sum_{i=1}^{n_k} \langle w \rangle_k c_k \right). \qquad (42)$$

Referring to the hypothesis (39), and remembering that the summation of $w_i$ values is equal to 1, values of $\langle w \rangle_k$ can be expressed as

$$\langle w \rangle_k \approx \left( \sum_{i=1}^{m} n_i c_i \right)^{-1} c_k \qquad (43)$$
<figure><img src="image_13.png" /><figcaption>**Figure 12:** Weights values $w_i$ versus $\cos(\theta_i)$ values of a PD tessellation (left) and a Lloyd relaxed PD one (right)</figcaption></figure>

As a consequence, it is necessary to compute the correlation coefficient (which would differs from zero) between weigths $w_i$ and $\cos(\theta_i)$ sets for the computation of the probabilistic density function of the tortuosity.

Splitting $I$, the interval of accessible values for $\cos(\theta_i)$, into $m$ intervals $I_k$, $k \in [1; m]$ of equal lengths, the sum at denominator of ?? reads:

$$ \sum_{i=1}^{n} w_i \cos(\theta_i) = \sum_{k=1}^{m} \left( \sum_{i=1}^{n_k} w_i \cos(\theta_i) \right), \qquad (43) $$

where $n_k$ corresponds to the number of $\cos(\theta)$ realizations which are in $I_k$.

Therefore if $n_k$ is large enough and according to the central limit theorem, the sum can be approximated by:

$$ \sum_{i=1}^{n} w_i \cos(\theta_i) \approx \sum_{k=1}^{m} n_k \langle w \cos(\theta) \rangle_k = \sum_{k=1}^{m} \left( \sum_{i=1}^{n_k} \langle w \cos(\theta) \rangle_k \right), \qquad (44) $$

where notation $\langle X \rangle_k$ stands for the expected value of the random variable $X$ inside the interval $I_k$ and $\theta$ and $w$ are respectively the random variables associated to the values $\theta_i$ and $w_i$.

If the interval $I_k$ is narrow enough, values $\cos(\theta_i)$ in $I_k$ can be considered as equal to a constant value corresponding to the center $c_k$ of interval $I_k$. Then $\langle w \cos(\theta) \rangle_k \approx \langle w \rangle_k c_k$ and estimate (??) becomes:

$$ \sum_{i=1}^{n} w_i \cos(\theta_i) \approx \sum_{k=1}^{m} \left( \sum_{i=1}^{n_k} \langle w \rangle_k c_k \right). \qquad (45) $$

Referring to the hypothesis (??), and remembering that the summation of $w_i$ values is equal to 1, values of $\langle w \rangle_k$ can be expressed as follows:

$$ \langle w \rangle_k \approx \frac{c_k}{\sum_{i=1}^{m} n_i c_i}, \qquad (46) $$

where $c_i$ and $n_i$ are respectively the center and the number of elements of the interval $I_i$.

Therefore the following estimate is accurate when the number of edges $n$ is large enough and considering that hypothesis (??) is true:

$$ \sum_{i=1}^{n} w_i \cos(\theta_i) \approx \sum_{i=1}^{n} \frac{\cos(\theta_i)^2}{\sum_{k=1}^{n} \cos(\theta_k)}. \qquad (47) $$
where $c_i$ and $n_i$ are respectively the center and the number of elements of the interval $I_i$. Therefore, the following estimate is accurate when the number of edges $n$ is large enough and considering that hypothesis (39) is true:

$$ \sum_{i=1}^{n} w_i \cos \theta_i \approx \left( \sum_{k=1}^{n} \cos \theta_k \right)^{-1} \sum_{i=1}^{n} \cos^2 \theta_i. \qquad (44) $$

This estimate leads to the approximation

$$ T_{\text{St}}(n) \approx \left( \sum_{i=1}^{n} \cos^2 \theta_i \right)^{-1} \sum_{i=1}^{n} \cos \theta_i \qquad (45) $$

for the tortuosity of an n-edges path following the inclination rule (26). In the sequel the numerator and the denominator in Equation (45) are respectively denoted by X and Y. Using the academic result $\sigma_{X+Y}^2 = \sigma_X^2 + \sigma_Y^2 + 2\text{Cov}(X, Y)$ on the variance of the sum of two correlated variables, the correlation coefficient between variables X and Y can be expressed as

$$ \rho_{X,Y} = \frac{\text{Cov}(X, Y)}{\sigma_X \sigma_Y} = \frac{\sigma_{X+Y}^2 - \sigma_X^2 - \sigma_Y^2}{2\sigma_X \sigma_Y} \qquad (46) $$

where the variances are estimated by the central limit theorem as

$$ \sigma_{X+Y}^2 \approx n\sigma_{\cos \theta+\cos^2 \theta}^2, \quad \sigma_X^2 \approx n\sigma_{\cos \theta}^2 \quad \text{and} \quad \sigma_Y^2 \approx n\sigma_{\cos^2 \theta}^2. \qquad (47) $$

Relations (46) and (47) imply that $\rho_{X,Y}$ does not depend on $n$. Equation (36) allows then to directly estimate $\rho_{X,Y}$ in function of $\bar{\mathcal{E}}$.

Having estimated $\rho_{X,Y}$ and given the fact that $T_{\text{St}}(n)$ is a ratio between two correlated normal random variables, Hinkley (1969) gives the estimate of the probabilistic density function of $T_{\text{St}}(n)$

$$ f_{\text{St}}(\bar{\mathcal{E}}, n; T)(x) \approx \frac{\sqrt{1 - \rho_{X,Y}^2}}{\pi \sigma_X \sigma_Y a(x)^2} \exp\left(-\frac{c}{2(1 - \rho_{X,Y}^2)}\right) + \frac{b(x)d(x)}{\sqrt{2\pi}\sigma_X \sigma_Y a(x)^3} \left[ \phi\left(\frac{b(x)}{a(x)\sqrt{1 - \rho_{X,Y}^2}}\right) - \phi\left(-\frac{b(x)}{a(x)\sqrt{1 - \rho_{X,Y}^2}}\right) \right], \qquad (48) $$

where

$$ a(x) = \sqrt{\frac{x^2}{\sigma_X^2} - \frac{2\rho_{X,Y}x}{\sigma_X \sigma_Y} + \frac{1}{\sigma_Y^2}}, \qquad b(x) = \frac{\mu_X x}{\sigma_X^2} - \frac{\rho_{X,Y}(\mu_X + \mu_Y x)}{\sigma_X \sigma_Y} + \frac{\mu_Y}{\sigma_Y^2} \qquad (49) $$

$$ c = \frac{\mu_X^2}{\sigma_X^2} - \frac{2\rho_{X,Y}\mu_X\mu_Y}{\sigma_X \sigma_Y} + \frac{\mu_Y^2}{\sigma_Y^2}, \qquad d(x) = \exp\left(\frac{b(x)^2 - ca(x)^2}{2(1 - \rho_{X,Y}^2)a(x)^2}\right) \qquad (50) $$

$$ \mu_X = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos x f_{\text{St}}(\bar{\mathcal{E}}; \theta)(x) dx, \qquad \mu_Y = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos^2 x f_{\text{St}}(\bar{\mathcal{E}}; \theta)(x) dx, \qquad (51) $$

with the normal cumulative distribution function

$$ \phi(y) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{y} \exp(-\frac{1}{2}u^2) du. \qquad (52) $$

Figure 13 shows that Equation (48) is accurate for PD tessellations while considering n-edges path following the inclination rule (26) with $n \ge 20$. For Lloyd relaxed PD tessellation, Equation (48) slightly underestimates the actual tortuosity of paths following the inclination rule.

Figure 14 shows that estimate (48) leads to an overestimate of the actual geodesic tortuosity. This overestimate is explained by the fact that the geodesic path is, by definition, always associated to the minimum tortuosity. In addition, the overestimate tends to decrease as the quality index increases. This observation is in accordance with Figure 9 which shows that the inclination rule
This estimate leads to the following approximation for the tortuosity of an *n*-edges path following the inclination rule (??):

$$ T_{\text{St}}(n) \approx \frac{\sum_{i=1}^{n} \cos(\theta_i)}{\sum_{i=1}^{n} \cos(\theta_i)^2}. \qquad (48) $$

In the sequel the numerator and the denominator in ?? are respectively denoted by X and Y. Using an academic result on the variance of the sum of two correlated variables:

$$ \sigma_{X+Y}^2 = \sigma_X^2 + \sigma_Y^2 + 2\text{Cov}(X, Y), \qquad (49) $$

the correlation coefficient between variables X and Y can be expressed as follows:

$$ \rho_{X,Y} = \frac{\text{Cov}(X, Y)}{\sigma_X \sigma_Y} = \frac{\sigma_{X+Y}^2 - \sigma_X^2 - \sigma_Y^2}{2\sigma_X \sigma_Y}, \qquad (50) $$

where variances are estimated by the central limit theorem as follows:

$$ \sigma_{X+Y}^2 \approx n\sigma_{\cos(\theta)+\cos(\theta)^2}^2 \quad \sigma_X^2 \approx n\sigma_{\cos(\theta)}^2 \quad \text{and} \quad \sigma_Y^2 \approx n\sigma_{\cos(\theta)^2}. \qquad (51) $$

Relations (??) and (??) imply that $\rho_{X,Y}$ does not depend of $n$. ?? allows then to directly estimate $\rho_{X,Y}$ in function of $\bar{\mathcal{E}}$.

Having estimated $\rho_{X,Y}$ and given the fact that $T_{\text{St}}(n)$ is a ratio between two correlated normal random variables, Hinkley (1969) gives the estimate of the probabilistic density function of $T_{\text{St}}(n)$:

$$ f_{\text{St}}(\bar{\mathcal{E}}, n; T)(x) \approx \frac{\sqrt{1 - \rho_{X,Y}^2}}{\pi \sigma_X \sigma_Y a(x)^2} \exp \left( -\frac{c}{2(1 - \rho_{X,Y}^2)} \right) + \frac{b(x)d(x)}{\sqrt{2\pi}\sigma_X \sigma_Y a(x)^3} \left[ \phi \left( \frac{b(x)}{a(x)\sqrt{1 - \rho_{X,Y}^2}} \right) - \phi \left( \frac{b(x)}{a(x)\sqrt{1 - \rho_{X,Y}^2}} \right) \right], \qquad (52) $$

where

$$ a(x) = \sqrt{\frac{x^2}{\sigma_X^2} - \frac{2\rho_{X,Y}x}{\sigma_X \sigma_Y} + \frac{1}{\sigma_Y^2}}, \quad b(x) = \frac{\mu_X x}{\sigma_X^2} - \frac{\rho_{X,Y}(\mu_X + \mu_Y x)}{\sigma_X \sigma_Y} + \frac{\mu_Y}{\sigma_Y^2}, \qquad (53) $$

$$ c = \frac{\mu_X^2}{\sigma_X^2} - \frac{2\rho_{X,Y}\mu_X\mu_Y}{\sigma_X \sigma_Y} + \frac{\mu_Y^2}{\sigma_Y^2}, \quad d(x) = \exp\left(\frac{b(x)^2 - ca(x)^2}{2(1 - \rho_{X,Y}^2)a(x)^2}\right), \qquad (54) $$

$$ \mu_X = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos(x) f_{\text{St}}(\bar{\mathcal{E}}; \theta)(x) dx, \quad \mu_Y = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos(x)^2 f_{\text{St}}(\bar{\mathcal{E}}; \theta)(x) dx, \qquad (55) $$

and $\phi$ is the normal cumulative distribution function:

$$ \phi(y) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{y} \exp\left(-\frac{1}{2}u^2\right) du. \qquad (56) $$

?? shows that ?? is accurate for PD tessellations while considering *n*-edges path following the inclination rule (??) with $n \ge 20$. For Lloyd relaxed PD tessellation, ?? slightly underestimates the actual tortuosity of paths following the inclination rule.
<figure><img src="image_14.png" /><figcaption>**Figure 13** Distributions of the tortuosity of *n*-edges paths following the inclination rule. Top: distributions obtained from a PD tessellation (τ = 500). Bottom: distributions obtained from a Lloyd relaxed PD tessellation (E = 0.86, τ = 500). Left: 20-edges path, middle: 50-edges path and right: 100-edges path. Black line: tortuosity pdf estimate (48).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14** Accuracy of the tortuosity probability density function. Left: geodesic tortuosity distribution of a PD tessellation (τ = 500). Right: geodesic tortuosity distribution of a Lloyd relaxed PD tessellation (E = 0.89, τ = 500). Black lines: geodesic tortuosity pdf estimates (48).</figcaption></figure>

is more often respected by meshes having a high quality index. As a result, estimate (48) appears to be accurate for unstructured Delaunay-like meshes suitable to finite element simulations (see Figure 14, right). This statement as been verified by a complementary study carried out on meshes generated by the Gmsh Delaunay algorithm.

Figure 15 shows comparisons between the geodesic tortuosity pdf estimate (48) and distributions obtained by a Monte Carlo simulation using the Dijkstra algorithm on 1000 Gmsh Delaunay meshes for several path lengths. The pdf estimates correctly mimic the geodesic tortuosity distributions of Gmsh Delaunay meshes with a slight overestimate: the average value obtained for Gmsh meshes is about 1.043 while estimate (48) gives an average value of 1.051.

It could have been possible to use regressions on the results of the Monte Carlo simulation to obtain estimates of the tortuosity value. However, this method comes with an additional difficulty: finding the underlying tortuosity probability density distribution.

# 7 Practical use of the results in computational fracture mechanics

## 7.1 Application of the mesh density estimate

When using finite initial cohesive stiffnesses (intrinsic cohesive laws) in a CVFE simulation, an artificial compliance is introduced between each finite elements. In turn, an overall stiffness loss is observed. The loss may be reduced by considering sufficiently high initial cohesive stiffnesses.
?? shows that estimate (??) leads to an overestimate of the actual geodesic tortuosity. This overestimate is explained by the fact that the geodesic path is, by definition, always associated to the minimum tortuosity.

In addition, the overestimate tends to decrease as the quality index increases. This observation is in accordance with ?? which shows that the inclination rule is more often respected by meshes having a high quality index.

As a result, estimate (??) appears to be accurate for unstructured Delaunay-like meshes suitable to finite element simulations (see ??, right). This statement as been verified by a complementary study carried out on meshes generated by the Gmsh Delaunay algorithm. ?? shows comparisons between the geodesic tortuosity pdf estimate (??) and distributions obtained by a Monte Carlo simulation using the Dijkstra algorithm on 1000 Gmsh Delaunay meshes for several path lengths.

<figure><img src="image_14.png" /><figcaption>**Figure 13:** Distributions of the tortuosity of *n*-edges paths following the inclination rule. Top: distributions obtained from a PD tessellation (τ = 500). Bottom: distributions obtained from a Lloyd relaxed PD tessellation (E = 0.86, τ = 500). Left: 20-edges path, middle: 50-edges path and right: 100-edges path. Black line: tortuosity pdf estimate (??).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14:** Left: geodesic tortuosity distribution of a PD tessellation (τ = 500). Right: geodesic tortuosity distribution of a Lloyd relaxed PD tessellation (E = 0.89, τ = 500). Black lines: geodesic tortuosity pdf estimates (??).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 15** Tortuosity pdf estimate compared to tortuosity histograms associated with Gmsh meshes. top: histograms of geodesic tortuosities of cracks containing distinct number of edges *n* computed via a Monte Carlo simulation using Dijkstra's algorithm. Black lines: geodesic tortuosity pdf estimates (48); bottom: zoom on particular shortest paths obtained on meshes generated by the Gmsh Delaunay algorithm (from left to right: paths containing respectively 10, 20 and 30 edges). Results obtained on Gmsh meshes ($\mathcal{E}$ = 0.89).</figcaption></figure>

Blal et al. 2012b have proposed practical formulas for choosing initial cohesive stiffnesses (hereafter denoted by $C_N < \infty$ and $C_T < \infty$ respectively for the normal and tangential decohesion modes) ensuring to maintain as enough rigidity as wanted in a homogeneous linear elastic isotropic material (Young modulus *E*, Poisson ratio *v*).

Denoting by *R* the “tolerated” ratio between the apparent overall stiffness obtained by using the CVFE approach and the theoretical one, the practical formulas read

$$ \frac{C_N}{E} \geq Z\left(\frac{R}{1-R}\right)\left(\frac{1}{3-6v}\right) \qquad (53) $$

and

$$ \frac{C_T}{C_N} = 2\left(\frac{1-2v}{1+3v}\right). \qquad (54) $$

Equation (53) shows that the initial cohesive stiffness in normal mode has to be chosen proportionally to the mesh density *Z*. Equation (54) permits to maintain the overall Poisson effect. For random planar triangular meshes suitable for finite element simulations, the mesh density is given by Equation (25) as function of the underlying point process intensity $\tau$. Additional mesh densities associated to several structured meshes are given in Appendix C.

## 7.2 Application of the geodesic tortuosity estimate

In CVFE approach, crack paths are restrained to propagate along elements boundaries (edges). If the localization and/or the orientation of the crack path is not known a priori, an isotropic tessellation is used, typically a Delaunay tessellation. A Delaunay tessellation seldom includes the theoretical path followed by the crack. In turn, a more tortuous crack path is obtained by the CVFE approach leading to an overestimate of the crack length (or crack surface in 3D).

As the theoretical overall fracture energy $W_{th}$ is proportional to the crack length times the specimen width (in 2D), the computed overall fracture energy $W_{comp}$ is also overestimated. It is reasonably expected that the fracture energy and equivalently the crack length overestimates are
<figure><img src="image_16.png" /><figcaption>**Figure 15:** Obtained on Gmsh meshes ($\mathcal{E} = 0.89$): top: histograms of geodesic tortuosities of cracks containing distinct number of edges $n$ computed via a Monte Carlo simulation using Dijkstra's algorithm. Black lines: geodesic tortuosity pdf estimates (??); bottom: zoom on particular shortest paths obtained on meshes generated by the Gmsh Delaunay algorithm (from left to right: paths containing respectively 10, 20 and 30 edges).</figcaption></figure>

The pdf estimates correctly mimic the geodesic tortuosity distributions of Gmsh Delaunay meshes with a slight overestimate: the average value obtained for Gmsh meshes is about 1.043 while estimate (??) gives an average value of 1.051.

It could have been possible to use regressions on the results of the Monte Carlo simulation to obtain estimates of the tortuosity value. However, this method comes with an additional difficulty: finding the underlying tortuosity probability density distribution.

# 7 Practical use of the results in computational fracture mechanics

## 7.1 Application of the mesh density estimate

When using finite initial cohesive stiffnesses (intrinsic cohesive laws) in a CVFE simulation, an artificial compliance is introduced between each finite elements. In turn, an overall stiffness loss is observed. The loss may be reduced by considering sufficiently high initial cohesive stiffnesses.

Blal et al. 2012b have proposed practical formulas for choosing initial cohesive stiffnesses (hereafter denoted by $C_N < \infty$ and $C_T < \infty$ respectively for the normal and tangential decohesion modes) ensuring to maintain as enough rigidity as wanted in a homogeneous linear elastic isotropic material (Young modulus $E$, Poisson ratio $\nu$).

Denoting by $R$ the “tolerated” ratio between the apparent overall stiffness obtained by using the CVFE approach and the theoretical one, the practical formulas read:

$$ \frac{C_N}{E} \geq Z \left( \frac{R}{1-R} \right) \left( \frac{1}{3-6\nu} \right) \qquad (57) $$

and

$$ \frac{C_T}{C_N} = 2 \left( \frac{1-2\nu}{1+3\nu} \right) \qquad (58) $$

?? shows that the initial cohesive stiffness in normal mode has to be chosen proportionally to the mesh density $Z$. ?? permits maintaining the overall Poisson effect.
of the same order of magnitude as the tortuosity estimate $T_{St}(n)$. Thus, denoting by $L_{th}$ the theoretical length of the crack path and by $L_{comp}$ the computed one, the relations

$$ \frac{W_{\text{comp}}}{W_{\text{th}}} (n) = \frac{L_{\text{comp}}}{L_{\text{th}}} (n) \approx T_{\text{St}}(n) \qquad (55) $$

are obtained, where *n* is the number of edges composing the crack path.

From Equation (48), confidence intervals of geodesic tortuosity estimate $T_{St}(n)$ as well as its expected value can be derived and estimated for distinct values of *n*, see Figure 16. Table 2

<figure><img src="image_17.png" /><figcaption>**Figure 16** Estimate of the geodesic tortuosity expected value (continuous line) and 95 % (dashed line) and 99 % (dotted line) confidence intervals in function of the efficiency index. Left: for 100-edges paths. Right: for 1000-edges paths.</figcaption></figure>

quantifies the dependency of the geodesic tortuosity estimate $T_{St}$ distribution to the number of edges *n* composing the crack path. Considering Equation (55) and table 2, we expect that an

<table><caption>**Table 2** Estimates of the geodesic tortuosity expected value and 95 % and 99 % confidence intervals in function of the number of edges ($\mathcal{E} = 0.89$).</caption><thead><tr><th>Number of edges</th><th>Mean tortuosity</th><th>95 % certainty interval</th><th>99 % certainty interval</th></tr></thead><tbody><tr><td>20</td><td>1.0512</td><td>[1.0279; 1.0760]</td><td>[1.0212; 1.0844]</td></tr><tr><td>50</td><td>1.0510</td><td>[1.0362; 1.0665]</td><td>[1.0317; 1.0716]</td></tr><tr><td>100</td><td>1.0510</td><td>[1.0404; 1.0617]</td><td>[1.0371; 1.0653]</td></tr><tr><td>1000</td><td>1.0509</td><td>[1.0474; 1.0542]</td><td>[1.0464; 1.0553]</td></tr></tbody></table>

overestimate of the fracture energy of about 5 % may be obtained while using an unstructured Delaunay-like mesh in a CVFE simulation.

# 8 Conclusion

The main objective of this paper was to derive accurate estimates for both the mesh density and geodesic tortuosity of planar triangular tessellations suitable for finite element simulations in the context of cohesive-based fracture mechanics. The following results were obtained.

* New estimates for the edge length and the triangle top angle distributions were proposed for Poisson-Delaunay tessellations and Lloyd relaxed ones for arbitrary efficiency indexes and intensities, Equations (14) and (15).

* The edge length and triangle top angle distributions of random triangular meshes obtained by standard mesh generators (Abaqus, Gmsh, Cast3M and Triangle) have been compared to these estimates depending on the intensity of the underlying point process and its efficiency index. These estimates are accurate for meshes with efficiency index at least less than 0.92. In other words, estimates are no longer accurate for semi-structured meshes (obtained for example by the GMSH “Frontal” algorithm and Abaqus). Moreover, these estimates are robust enough to be used at a local scale, see Appendix A.

* Based on these practical estimates, the distribution of the mesh density was derived as a normal law in Equation (21). Its mean grows as the square root of the intensity of the point process, while
For random planar triangular meshes suitable for finite element simulations, the mesh density is given by ?? as function of the underlying point process intensity $\tau$. Additional mesh densities associated to several structured meshes are given in ??.

## 7.2 Application of the geodesic tortuosity estimate

In CVFE approach, crack paths are restrained to propagate along elements boundaries (edges). If the localization and/or the orientation of the crack path is not known a priori, an isotropic tessellation is used, typically a Delaunay tessellation.

A Delaunay tessellation seldom includes the theoretical path followed by the crack. In turn, a more tortuous crack path is obtained by the CVFE approach leading to an overestimate of the crack length (or crack surface in 3D).

As the theoretical overall fracture energy $W_{th}$ is proportional to the crack length times the specimen width (in 2D), the computed overall fracture energy $W_{comp}$ is also overestimated.

It is reasonably expected that the fracture energy and equivalently the crack length overestimates are of the same order of magnitude as the tortuosity estimate $T_{St}(n)$. Thus, denoting by $L_{th}$ the theoretical length of the crack path and by $L_{comp}$ the computed one, the following relations are obtained:

$$ \frac{W_{comp}}{W_{th}}(n) = \frac{L_{comp}}{L_{th}}(n) \approx T_{St}(n), \qquad (59) $$

where $n$ is the number of edges composing the crack path.

From ??, confidence intervals of geodesic tortuosity estimate $T_{St}(n)$ as well as its expected value can be derived and estimated for distinct values of $n$ (see ??).

<figure><img src="image_17.png" /><figcaption>Figure 16: Estimate of the geodesic tortuosity expected value (continuous line) and 95% (dashed line) and 99% (dotted line) confidence intervals in function of the efficiency index. Left: for 100-edges paths. Right: for 1000-edges paths.</figcaption></figure>

?? quantifies the dependency of the geodesic tortuosity estimate $T_{St}$ distribution to the number of edges $n$ composing the crack path.

Considering ??, and ??, we expect that an overestimate of the fracture energy of about 5% may be obtained while using an unstructured Delaunay-like mesh in a CVFE simulation.

## 8 Conclusion

The main objective of this paper was to derive accurate estimates for both the mesh density and geodesic tortuosity of planar triangular tessellations suitable for finite element simulations in the context of cohesive-based fracture mechanics. The following results were obtained.
its variance does not depend on this intensity but decreases when the area of the discretized
domain increases. The expected value of the mesh density does not depend on the global
efficiency index. A dependence to the efficiency index is observed for the variance of the mesh
density. However, for standard meshes used in finite elements simulations, the variance of the
mesh density appears to be negligible in comparison to its expected value. It is then stated
that the mesh density can be viewed as a singular function of the intensity of the underlying
point-process (equal to the mesh density expected value). The mesh densities of usual structured
triangular meshes are also given in Appendix C.

* For unstructured Delaunay-like meshes (whose edge orientations may be considered uniformly distributed), a tortuosity estimate has been derived from an hypothetical “orientation rule” followed by each edge of a shortest path between two nodes. Since the orientation rule is not exactly fulfilled by shortest paths, the tortuosity is slightly overestimated. However, for application we have in mind, the overestimate is small enough for considering that the geodesic tortuosity estimate (48) is accurate.

* Overestimate of both the overall fracture energy and cracks tortuosities computed by a CVFE simulation are expected while using an unstructured Delaunay-like mesh. Using the geodesic tortuosity estimate $T_{\text{St}}$, we deduced that this overestimate is approximately equal to 5 %. We emphasize the fact that the results of this work are valid only for “sufficiently disordered” homogeneous planar triangular meshes. They may be improved by considering the case of inhomogeneous Poisson point processes or three-dimensional Poisson-Delaunay tessellations.

In three-dimensional tessellations, the notion of geodesic plane might be introduced as
well as an associated tortuosity. To the best of our knowledge, there are no known results for
estimating the three-dimensional tortuosity distribution. To provide such an estimate seems to be
far more challenging than providing the estimate of the two-dimensional tortuosity.

A Local robustness of the edge length and triangle top angle distribu-
tions Equation (14) and Equation (15)

Since we have in mind fracture mechanics of complex microstructures, involving local hetero-
geneities, gradients, cavities, etc., we have to check at which scale or for which local intensity $\tau$
estimates (14) and (15) are accurate enough. The accuracy is quantified by the standard deviations
of the means $\mu_{\text{St}}(\tau; l)$ and $\mu_{\text{St}}(\alpha)$, respectively denoted by $\sigma_{\mu_{\text{St}}}(\tau; l)$ and $\sigma_{\mu_{\text{St}}}(\alpha)$, and by the
standard deviations of the variances $\sigma_{\text{St}}^2(\tau; l)$ and $\sigma_{\text{St}}^2(\alpha)$, respectively denoted by $\sigma_{\sigma_{\text{St}}^2}(\tau; l)$ and
$\sigma_{\sigma_{\text{St}}^2}(\alpha)$, for 100 samples. Each sample corresponds to a subset of a Gmsh random triangular mesh
("MeshAdapt" strategy, efficiency index $\mathcal{E} \approx 0.86$) consisting of a disk with random center and
radius $R$, see Figure A.1. The mesh of Figure A.1 has an efficiency index of $\mathcal{E} \approx 0.86$.

<figure><img src="image_18.png" /><figcaption>**Figure A.1** Six circular subset samples with random center and radius $R \approx 3.3\mu_{\text{St}}(\tau; l)$.</figcaption></figure>

Figure A.2 shows that, when the radius *R* is larger than the expected mean of the edge length, *R* ≥ *μ*<sub>St</sub>(*τ*; *l*), the standard deviation of the mean *σ*<sub>μ<sub>St</sub></sub>(*τ*; *l*) is lower than about 5 % of the actual edge length mean and the standard deviations of the variances *σ*<sub>St</sub><sup>2</sup>(*l*) and *σ*<sub>St</sub><sup>2</sup>(*α*) are respectively lower than 4 % and 12 % of the corresponding actual variances. In addition, the standard deviation of the mean *μ*<sub>St</sub>(*α*), which is not represented in Figure A.2, is lower than the machine epsilon
<table><caption>**Table 2:** Estimates of the geodesic tortuosity expected value and 95% and 99% confidence intervals in function of the number of edges ($\mathcal{E} = 0.89$).</caption><thead><tr><td>Number of edges</td><td>Mean tortuosity</td><td>95% certainty interval</td><td>99% certainty interval</td></tr></thead><tbody><tr><td>20</td><td>1.0512</td><td>[1.0279; 1.0760]</td><td>[1.0212; 1.0844]</td></tr><tr><td>50</td><td>1.0510</td><td>[1.0362; 1.0665]</td><td>[1.0317; 1.0716]</td></tr><tr><td>100</td><td>1.0510</td><td>[1.0404; 1.0617]</td><td>[1.0371; 1.0653]</td></tr><tr><td>1000</td><td>1.0509</td><td>[1.0474; 1.0542]</td><td>[1.0464; 1.0553]</td></tr></tbody></table>

* New estimates for the edge length and the triangle top angle distributions were proposed for Poisson-Delaunay tessellations and Lloyd relaxed ones for arbitrary efficiency indexes and intensities (?? and ??).

* The edge length and triangle top angle distributions of random triangular meshes obtained by standard mesh generators (Abaqus, Gmsh, Cast3M and Triangle) have been compared to these estimates depending on the intensity of the underlying point process and its efficiency index. These estimates are accurate for meshes with efficiency index at least less than 0.92. In other words, estimates are no longer accurate for semi-structured meshes (obtained for example by the GMSH “Frontal” algorithm and Abaqus). Moreover, these estimates are robust enough to be used at a local scale (see ??).

* Based on these practical estimates, the distribution of the mesh density was derived as a normal law (??). Its mean grows as the square root of the intensity of the point process, while its variance does not depend on this intensity but decreases when the area of the discretized domain increases. The expected value of the mesh density does not depend on the global efficiency index. A dependence to the efficiency index is observed for the variance of the mesh density. However, for standard meshes used in finite elements simulations, the variance of the mesh density appears to be negligible in comparison to its expected value. It is then stated that the mesh density can be viewed as a singular function of the intensity of the underlying point-process (equal to the mesh density expected value). The mesh densities of usual structured triangular meshes are also given in ??.

* For unstructured Delaunay-like meshes (whose edge orientations may be considered uniformly distributed), a tortuosity estimate has been derived from an hypothetic “orientation rule” followed by each edge of a shortest path between two nodes. Since the orientation rule is not exactly fulfilled by shortest paths, the tortuosity is slightly overestimated. However, for application we have in mind, the overestimate is small enough for considering that the geodesic tortuosity estimate (??) is accurate.

* Overestimate of both the overall fracture energy and cracks tortuosities computed by a CVFE simulation are expected while using an unstructured Delaunay-like mesh. Using the geodesic tortuosity estimate $T_{st}$, we deduced that this overestimate is approximatively equal to 5%.

We emphasize the fact that the results of this work are valid only for “sufficiently disordered” homogeneous planar triangular meshes. They may be improved by considering the case of inhomogeneous Poisson point processes or three-dimensional Poisson-Delaunay tessellations.

In three-dimensional tessellations, the notion of geodesic plane might be introduced as well as an associated tortuosity. To the best of our knowledge, there are no known results for estimating the three-dimensional tortuosity distribution. To provide such an estimate seems to be far more challenging than providing the estimate of the two-dimensional tortuosity.

## Acknowledgement

The authors are indebted to P. Querre for fruitful conversations and useful comments.
whatever $R \ge \mu_{\text{St}}(\tau; l)$. Finally, the local robustness of Equations (14) and (15) is such that they
can be applied at the length scale of a few triangles in an unstructured Delaunay-like mesh.

<figure><img src="image_19.png" /><figcaption>**Figure A.2** Local robustness of the estimates of the edge length and triangle top angle probability density functions. Left: standard deviation of the mean of the edge length divided by its actual mean with respect to $R/\mu_{\text{St}}(\tau; l)$. Right: standard deviation of the variance of the edge length (resp. the triangle top angle) divided by its actual variance with respect to $R/\mu_{\text{St}}(\tau; l)$ (black dots/resp. gray dots). The standard deviations are computed over 100 samples of circular random subsets belonging to a Gmsh triangular mesh ("MeshAdapt" strategy, efficiency index $\mathcal{E} \approx 0.86$).</figcaption></figure>

B Some hints on the range of validity for estimates (14) and (15)

The proposed estimates (14) and (15) are strictly devoted to unstructured Delaunay-like meshes.
These estimates are not able to deal with:

* *Badly shaped elements.* Efficiency indexes lower than the minimal value $\mathcal{E}_{\text{PD}}$ corresponding to PD tessellations (see Equation (11)) can not be taken into account by these estimates. They are restricted to “good” enough meshes (i.e. with an efficiency index greater than 0.725).

* *Structured or quasi-structured meshes.* Structured meshes lead to discrete values of the edge length and triangle top angle sets. Their distributions thus correspond to successive Dirac delta functions. Estimates (14) and (15) are restricted to “sufficiently unstructured” meshes (i.e. with an efficiency index at least less than 0.92).

To underline these two limitations, planar hand-generated meshes are considered as counterex-
amples. These meshes consist in cross-triangular square patterns stretched in one direction
(degradation of the shape) with random perturbation of the inner points (in order to control the
regularity), tessellations remaining valid. Three cases are considered.

With highly random perturbation, the mesh is no more structured and estimates (14) and
(15) exhibit a rough validity. In addition, these estimates do not accurately predict the tail of
the edge length and tirangle angle distributions when badly shaped elements are numerous;
see Figure B.3, left. With small random perturbation, the edge length distribution consists in
multi-modal normal distribution. Estimates (14) and (15) are able to predict the mode with higher
intensity, but are not able to predict the other modes; see Figure B.3, middle. With even smaller
random perturbation, a quasi-structured mesh is obtained and estimates (14) and (15) fail to
predict the quasi-Dirac combs; see Figure B.3, right.

Thus, the use of estimates (14) and (15) is restricted to standard random triangular tessellations
whose efficiency index is greater than or equal to the one of a PD tessellation and at least less
than 0.92 (see Section 4.2).

C Complementary results for structured triangular meshes

The estimate of the total length of the edges of a tessellation containing $N_P$ points regularly
distributed in an area A has been presented in Watanabe (2008) for various measures. Their results
concern some triangular meshes, square meshes and hexagonal meshes. We focus on the mesh
## References

Abaqus (2008). User's Manual Version 6.7

Barenblatt, G. I. (1962). "The Mathematical Theory of Equilibrium of Crack in Brittle Fracture". *Advances in Applied Mechanics* 7, pp. 55-129. DOI: 10.1016/S0065-2156(08)70121-2

Billingsley, P. (1995). *Probability and Measure*. Wiley-Interscience

Blal, N., L. Daridon, Y. Monerie, and S. Pagano (2011). "Criteria on the Artificial Compliance Inherent to the Intrinsic Cohesive Zone". *Comptes Rendus - Mecanique* 339.12, pp. 789-795. DOI: 10.1016/j.crme.2011.10.001

Blal, N., L. Daridon, Y. Monerie, and S. Pagano (2012a). "Artificial Compliance Inherent to the Intrinsic Cohesive Zone Models: Criteria and Application to Planar Meshes". *International Journal of Fracture* 178, pp. 71-83. DOI: 10.1007/s10704-012-9734-y

Blal, N., L. Daridon, Y. Monerie, and S. Pagano (2012b). "Micromechanical-Based Criteria for the Calibration of Cohesive Zone Parameters". *Journal of Computational and Applied Mathematics* 246, pp. 206-214. DOI: 10.1016/j.cam.2012.10.031

Bourdin, B., G. A. Francfort, and J.-J. Marigo (2008). "The Variational Approach to Fracture". *Journal of Elasticity* 91, pp. 5-148. DOI: 10.1007/s10659-007-9107-3

Cast3m (2011). *Finite Element Software*

Chaboche, J.-L., F. Feyel, and Y. Monerie (2001). "Interface Debonding Model: A Viscous Regularization with a Limited Rate Dependency". *International Journal of Solids and Structures* 38.18, pp. 3127-3160. DOI: 10.1016/S0020-7683(00)00053-6

Chang, R. C. and R. C. T. Lee (1984). "On the Average Length of Delaunay Triangulations". *Bit Numerical Mathematics* 24, pp. 269-273

Collins, R. (1968). "A Geometrical Sum Rule for Two-Dimensional Fluid Correlation Functions". *Journal of Physics C: Solid State Physics* 1.6, pp. 1461-1471. DOI: 10.1088/0022-3719/1/6/301

De Oliveira, S. G. (1988). "A Review on Delaunay Refinement Techniques". *ICCSA'12 Proceedings of the 12th International Conference on Computational Science and Its Applications* (2009). Berlin, Germany, pp. 172-187. DOI: 10.1007/978-3-642-31125-3_14

Dolbow, J., N. Moës, and T. Belytschko (2000). "Discontinuous Enrichment in Finite Elements with a Partition of Unity Method". *Finite Elements in Analysis and Design* 36.3-4, pp. 235-260. DOI: 10.1016/S0168-874X(00)00035-4

Du, Q., M. Emelianenko, and L. Ju (2006). "Convergence of the Lloyd Algorithm for Computing Centrical Voronoi Tessellations". *SIAM Journal on Numerical Analysis* 44.1, pp. 102-119. DOI: 10.1137/040617364

Dugdale, D. S. (1960). "Yielding of Steel Sheets Containing Slits". *Journal of the Mechanics and Physics and Solids* 8.2, pp. 100-104. DOI: 10.1016/0022-5096(60)90013-2

Freitag, L. and C. Ollivier-Gooch (1997). "Tetrahedral Mesh Improvement Using Swapping and Smoothing". *International Journal for Numerical Methods in Engineering* 40.21, pp. 3979-4002. DOI: 10.1002/(SICI)1097-0207(19971115)40:21<3979::AID-NME251>3.0.CO;2-9

George, P.-L. and P. Frey (2000). *Mesh Generation*. Hermes

Geuzaine, C. and J.-F. Remacle (2009). "Gmsh: A 3-d Finite Element Mesh Generator with Built-in Pre- and Post-Processing Facilities". *International Journal for Numerical Methods in Engineering* 79.11, pp. 1309-1331. DOI: 10.1002/nme.2579

Hinkley, D. V. (1969). "On the Ratio of Two Correlated Normal Random Variables". *Biometrika* 56.3, pp. 635-639. DOI: 10.1093/biomet/56.3.635

Lasry, D. and T. Belytschko (1988). "Localization Limiters in Transient Problems". *International Journal of Solids and Structures* 24.6, pp. 581-597. DOI: 10.1016/0020-7683(88)90059-5

Liu, G.-Y., W.-J. Xu, N. Govender, and D. N. Wilke (2021). "Simulation of Rock Fracture Process Based on GPU-accelerated Discrete Element Method". *Powder Technology* 377, pp. 640-656. DOI: 10.1016/j.powtec.2020.09.009

Liu, Y., W. Wang, B. Lévy, F. Sun, D.-M. Yan, L. Lu, and C. Yang (2009). "On Centroidal Voronoi Tessellation - Energy Smoothness and Fast Computation". *ACM Transactions on Graphics* 28.4, pp. 1-17. DOI: 10.1145/1559755.1559758

Lloyd, L. (1982). "Least Square Quantization in PCM". *IEEE Transactions on Information Theory* 28.2, pp. 129-137. DOI: 10.1109/TIT.1982.1056489
<figure><img src="image_20.png" /><figcaption>**Figure B.3** Limit of validity for estimates (14) and (15) in case of quasi-structured meshes. Pdf of the edge length (first row) and triangle top angle (second row) sets for hand-generated cross-triangular square meshes, stretched in the horizontal direction and randomly perturbed of the inner points. Left: stretching/elongation ratio of 3 with high distortion, efficiency index $\mathcal{E} \approx 0.78$, intensity $\tau = 466$. Middle: stretching ratio of 2 with intermediate distortion, efficiency index $\mathcal{E} \approx 0.84$, intensity $\tau = 472$. Right: stretching ratio of 1 with small distortion, efficiency index $\mathcal{E} \approx 0.87$, intensity $\tau = 458$.</figcaption></figure>

density Z of structured triangular meshes usually used in finite element discretizations (Tijssens et al. 2000).

The total length of edges in a structured triangular tessellation based on a hexagonal (H) point process reads (see Table C.1)

$$ \sum_{i=1}^{N_E} l_i = 3 \sqrt{\frac{2\tau}{\sqrt{3}}} A. \qquad (C.1) $$

The corresponding mesh density is thus

$$ Z_H = 3 \sqrt{\frac{2}{\sqrt{3}}} \sqrt{\tau}. \qquad (C.2) $$

Regular points distributed along a square lattice (see Table C.1) lead to a structured “Triangle Quadrilateral” (TQ) pattern whose the total length of the edges reads

$$ \sum_{i=1}^{N_E} l_i = (2 + \sqrt{2}) \sqrt{\tau} A, \qquad (C.3) $$

and the mesh density,

$$ Z_{TQ} = (2 + \sqrt{2}) \sqrt{\tau}. \qquad (C.4) $$

In fracture mechanics, “Cross-Triangle Quadrilateral” (CTQ) patterns are often used (see Table C.1). For that mesh topology, the mesh density is deduced from the previous one:

$$ Z_{CTQ} = 2(1 + \sqrt{2}) \sqrt{\tau}. \qquad (C.5) $$

This CTQ pattern can be stretched in one direction and the underlying quadrilateral is not a square but a rectangle (aspect ratio $s \ge 1$). In this case, the mesh density becomes

$$ Z_{CTQs} = \left( 2 \sqrt{1 + \frac{1}{s^2}} + 1 + \frac{1}{s} \right) \sqrt{\tau}. \qquad (C.6) $$
Lohner, R. and P. Parikh (1988). "Generation of the Three-Dimensional Unstructured Grids by Advancing-Front Method". *26th Aerospace Sciences Meeting (Jan. 11-14, 1988)*. Reno, NV, U.S.A., pp. 1135-1149. DOI: 10.2514/6.1988-515

Lu, X., J.-P. Bardet, and M. Huan (2009). "Numerical Solutions of Strain Localization with Nonlocal Softening Plasticity". *Computer Methods Applied Mechanics Engineering* 198.47-48, pp. 3702-3711. DOI: 10.1016/j.cma.2009.08.002

Manorosoa, Z. A., A. Chrysochoos, A. Jelea, Y. Monerie, and F. Perales (2023). "Atomistic study of the properties and behavior under uniaxial tensile stress of some representative UO₂ grain boundaries". *Computational Materials Science* 228, p. 112314. DOI: 10.1016/j.commatsci.2023.112314

Meijering, J. L. (1953). "Interface Area, Edge Length, and Number of Vertices in Crystal Components in Boolean Models". *Philips Research Report* 8, pp. 270-290

Miles, R. E. (1970). "On the Homogeneous Planar Poisson Point Process". *Mathematical Biosciences* 6, pp. 85-127. DOI: 10.1016/0025-5564(70)90061-1

Nguyen, N. H. T., H. H. Bui, G. D. Nguyen, and J. Kodikara (2017). "A Cohesive Damage-Plasticity Model for DEM and Its Application for Numerical Investigation of Soft Rock Fracture Properties". *International Journal of Plasticity* 98, pp. 175-196. DOI: 10.1016/j.ijplas.2017.07.008

Papoulia, K. D., S. A. Vavasis, and P. Ganguly (2006). "Spatial Convergence of Crack Nucleation Using a Cohesive Finite-Element Model on a Pinwheel-Based Mesh". *International Journal for Numerical Methods in Engineering* 67.1, pp. 1-16. DOI: 10.1002/nme.1598

Pele, K., J. Baccou, L. Daridon, J. Liandrat, T. Le Gouic, Y. Monerie, and F. Perales (2023). "A probabilistic model for fast-to-evaluate 2D crack path prediction in heterogeneous materials". *Journal of Theoretical, Computational and Applied Mechanics* 15, pp. 1-21. DOI: 10.46298/jtcam.8322

Perales, F., S. Bourgeois, A. Chrysochoos, and Y. Monerie (2008). "Two Field Multibody Method for Periodic Homogenization in Fracture Mechanics of Non Linear Composites". *Engineering Fracture Mechanics* 75.11, pp. 3378-3398. DOI: 10.1016/j.engfracmech.2007.07.017

Rebay, S. (1993). "Efficient Unstructured Mesh Generation by Means of Delaunay Triangulation and Bowyer-Watson Algorithm". *Journal of Computational Physics* 106.1, pp. 125-138. DOI: 10.1006/jcph.1993.1097

Rimoli, J., J. Rojas, and F. Khemani (2012). "On the Mesh Dependency of Cohesive Zone Models for Crack Propagation Analysis". *53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference*. Honolulu, Hawaii. DOI: 10.2514/6.2012-1536

Shewchuk, J. (2002). *What Is a Good Linear Finite Element? - Interpolation, Conditioning, Anisotropy, and Quality Measures*. In Proc. of the 11th International Meshing Roundtable

Shewchuk, J. (2005). *Triangle - a Two-Dimensional Quality Mesh Generator and Delaunay Triangulator*

Sibson, R. (1980). "The Dirichlet Tessellation as an Aid in Data Analysis". *Scandinavian Journal of Statistics* 7.1, pp. 14-20

Tijssens, M. G. A., L. Sluys, and E. van der Gissen (2000). "Numerical Simulation of Quasi-Brittle Fracture Using Damaging Cohesive Surfaces". *European Journal of Mechanics A/Solids* 19.5, pp. 761-779. DOI: 10.1016/S0997-7538(00)00190-X

Watanabe, D. (2008). "Evaluating the Configuration and the Travel Efficiency on Proximity Graphs as Transportation Networks". *Forma* 23.2, pp. 81-87

Xu, X.-P. and A. Needleman (1994). "Numerical Simulations of Fast Crack Growth in Brittle Solids". *Journal of the Mechanics and Physics of Solids* 42.9, pp. 1397-1434. DOI: 10.1016/0022-5096(94)90003-5
<figure><img src="image_21.png" /></figure>

$$Z_H = 3 \sqrt{\frac{2}{\sqrt{3}}} \sqrt{\tau}$$

<figure><img src="image_22.png" /></figure>

$$Z_{CTQ} = 2(1 + \sqrt{2})\sqrt{\tau}$$

<figure><img src="image_23.png" /></figure>

$$Z_{TQ} = (2 + \sqrt{2})\sqrt{\tau}$$

<figure><img src="image_24.png" /></figure>

$$Z_{CTQs} = (2\sqrt{1 + \frac{1}{s^2}} + 1 + \frac{1}{s})\sqrt{\tau}$$

<caption>**Table C.1** Mesh density of the structured triangular meshes usually used in finite element discretizations.</caption>

Table C.1 summarizes the mesh densities for structured meshes usually used in finite element discretizations. It is worth noting that these mesh densities for structured triangular meshes evolve as $\sqrt{\tau}$, exactly as the means of the mesh densities for random triangular meshes do. In particular, the mean $\mu_{St}(\tau; Z)$ of standard random triangular meshes given in Equation (21) is very close to the mesh density of TQ tessellations: the prefactor of $\mu_{St}(\tau; Z)$ is about 3.40, while the prefactor of $Z_{TQ}(\tau)$ is about 3.41:

$$Z_{TQ}(\tau) \approx \mu_{St}(\tau; Z). \qquad (C.7)$$

From the point of view of the mesh density, a Triangle Quadrilateral pattern accurately mimics a standard random triangular mesh (keeping in mind that the distribution of the edge length set is bimodal in TQ tessellations). Moreover, the mesh density Z of stretched CTQ (CTQs) decreases as the stretch factor s increases, see Figure C.4.

<figure><img src="image_25.png" /><figcaption>**Figure C.4** Mesh density Z as a function of the point process intensity $\tau$ for different types of structured triangular meshes. Gray line: Cross-Triangle Quadrilateral mesh (C.5). Gray dashed line: stretched Cross-Triangle Quadrilateral mesh (C.6) with stretch factor s = 2. Black dashed line: Triangle Quadrilateral mesh (C.4). Black line: hexagonal mesh (C.2).</figcaption></figure>

## References

* Abaqus (2008). User's Manual Version 6.7.

* Barenblatt, G. I. (1962). The mathematical theory of equilibrium of crack in brittle fracture. *Advances in Applied Mechanics* 7:55–129. [DOI], [HAL].

* Billingsley, P. (1995). *Probability and Measure*. Wiley-Interscience.

* Blal, N., L. Daridon, Y. Monerie, and S. Pagano (2011). Criteria on the artificial compliance inherent to the intrinsic cohesive zone. *Comptes Rendus Mécanique* 339(12):789–795. [DOI], [OA].

* Blal, N., L. Daridon, Y. Monerie, and S. Pagano (2012a). Artificial compliance inherent to the intrinsic cohesive zone models: criteria and application to planar meshes. *International Journal of Fracture* 178(1-2):71–83. [DOI], [HAL].

* Blal, N., L. Daridon, Y. Monerie, and S. Pagano (2012b). Micromechanical-based criteria for the calibration of cohesive zone parameters. *Journal of Computational and Applied Mathematics* 246:206–214. [DOI], [HAL].

* Bourdin, B., G. Francfort, and J.-J. Marigo (2008). The variational approach to fracture. *Journal of Elasticity* 91(1-3):5–148. [DOI], [OA].

* Cast3m (2011). *Finite Element Software*.

* Chaboche, J., F. Feyel, and Y. Monerie (2001). Interface debonding model: A viscous regularization with a limited rate dependency. *International Journal of Solids and Structures* 38(18):3127–3160. [DOI], [HAL].
# A Local robustness of the edge length and triangle top angle distributions ?? and ??

Since we have in mind fracture mechanics of complex microstructures, involving local heterogeneities, gradients, cavities, etc., we have to check at which scale or for which local intensity $\tau$ estimates (??) and (??) are accurate enough. The accuracy is quantified by the standard deviations of the means $\mu_{St}(\tau; l)$ and $\mu_{St}(\alpha)$, respectively denoted by $\sigma_{\mu_{St}}(\tau; l)$ and $\sigma_{\mu_{St}}(\alpha)$, and by the standard deviations of the variances $\sigma_{St}^2(\tau; l)$ and $\sigma_{St}^2(\alpha)$, respectively denoted by $\sigma_{\sigma_{St}^2}(\tau; l)$ and $\sigma_{\sigma_{St}^2}(\alpha)$, for 100 samples. Each sample corresponds to a subset of a Gmsh random triangular mesh (“MeshAdapt” strategy, efficiency index $\mathcal{E} \approx 0.86$) consisting of a disk with random center and radius $R$, see ?? . The mesh of ?? has an efficiency index of $\mathcal{E} \approx 0.86$.

<figure><img src="image_18.png" /><figcaption>**Figure 17:** Six circular subset samples with random center and radius $R \approx 3.3\mu_{St}(\tau; l)$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 18:** Left: standard deviation of the mean of the edge length divided by its actual mean with respect to $R/\mu_{St}(\tau; l)$. Right: standard deviation of the variance of the edge length (resp. the triangle top angle) divided by its actual variance with respect to $R/\mu_{St}(\tau; l)$ (black dots/resp. gray dots). The standard deviations are computed over 100 samples of circular random subsets belonging to a Gmsh triangular mesh (“MeshAdapt” strategy, efficiency index $\mathcal{E} \approx 0.86$).</figcaption></figure>

?? shows that, when the radius $R$ is larger than the expected mean of the edge length, $R \ge \mu_{St}(\tau; l)$, the standard deviation of the mean $\sigma_{\mu_{St}}(\tau; l)$ is lower than about 5% of the actual edge length mean and the standard deviations of the variances $\sigma_{\sigma_{St}^2}(l)$ and $\sigma_{\sigma_{St}^2}(\alpha)$ are respectively lower than 4% and 12% of the corresponding actual variances. In addition, the standard deviation of the mean $\mu_{St}(\alpha)$, which is not represented in ??, is lower than the machine epsilon whatever $R \ge \mu_{St}(\tau; l)$. Finally, the local robustness of ?? and ?? is such that they can be applied at the length scale of a few triangles in an unstructured Delaunay-like mesh.
Chang, R. C. and R. C. T. Lee (1984). On the average length of Delaunay triangulations. *BIT* 24(3):269-273. [DOI].

Collins, R (1968). A geometrical sum rule for two-dimensional fluid correlation functions. *Journal of Physics C: Solid State Physics* 1(6):1461-1471. [DOI].

Dolbow, J., N. Moës, and T. Belytschko (2000). Discontinuous enrichment in Finite Elements with a Partition of Unity Method. *Finite Elements in Analysis and Design* 36(3-4):235-260. [DOI], [HAL].

Du, Q., M. Emelianenko, and L. Ju (2006). Convergence of the Lloyd algorithm for computing centrical Voronoi tessellations. *SIAM Journal on Numerical Analysis* 44(1):102-119. [DOI].

Dugdale, D. S. (1960). Yielding of steel sheets containing slits. *Journal of the Mechanics and Physics of Solids* 8(2):100-104. [DOI].

Freitag, L. and C. Ollivier-Gooch (1997). Tetrahedral mesh improvement using swapping and smoothing. *International Journal for Numerical Methods in Engineering* 40(21):3979-4002. [DOI].

Frey, P. J. and P.-L. George (2000). *Mesh Generation*. Wiley. [DOI].

Geuzaine, C. and J.-F. Remacle (2009). Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. *International Journal for Numerical Methods in Engineering* 79(11):1309-1331. [DOI], [HAL].

Hinkley, D. V. (1969). On the ratio of two correlated normal random variables. *Biometrika* 56(3):635-639. [DOI].

Lasry, D. and T. Belytschko (1988). Localization limiters in transient problems. *International Journal of Solids and Structures* 24(6):581-597. [DOI].

Liu, G.-Y., W.-J. Xu, N. Govender, and D. N. Wilke (2021). Simulation of rock fracture process based on GPU-accelerated Discrete Element Method. *Powder Technology* 377:640-656. [DOI].

Liu, Y., W. Wang, B. Lévy, F. Sun, D.-M. Yan, L. Lu, and C. Yang (2009). On centroidal Voronoi tessellation - Energy smoothness and fast computation. *ACM Transactions on Graphics* 28(4):1-17. [DOI], [HAL].

Lloyd, S. (1982). Least square quantization in PCM. *IEEE Transactions on Information Theory* 28(2):129-137. [DOI], [HAL].

Lohner, R. and P. Parikh (1988). Generation of the three-dimensional unstructured grids by advancing-front method. *26th Aerospace Sciences Meeting (Reno, Nevada, USA, Jan. 11-14, 1988)*. American Institute of Aeronautics and Astronautics, pp 1135-1149. [DOI].

Lu, X., J.-P. Bardet, and M. Huang (2009). Numerical solutions of strain localization with nonlocal softening plasticity. *Computer Methods in Applied Mechanics and Engineering* 198(47-48):3702-3711. [DOI].

Manorosoa, Z. A., A. Chrysochoos, A. Jelea, Y. Monerie, and F. Perales (2023). Atomistic study of the properties and behavior under uniaxial tensile stress of some representative UO₂ grain boundaries. *Computational Materials Science* 228:112314. [DOI].

Meijering, J. L. (1953). Interface area, edge length, and number of vertices in crystal components in boolean models. *Philips Research Report* 8:270-290.

Miles, R. (1970). On the homogeneous planar Poisson point process. *Mathematical Biosciences* 6:85-127. [DOI].

Nguyen, N., H. Bui, G. Nguyen, and J. Kodikara (2017). A cohesive damage-plasticity model for DEM and its application for numerical investigation of soft rock fracture properties. *International Journal of Plasticity* 98:175-196. [DOI].

Oliveira, S. Gonzaga de (1988). A review on Delaunay refinement techniques. *Lecture Notes in Computer Science* (2009). Springer, pp 172-187. [DOI].

Papoulia, K. D., S. A. Vavasis, and P. Ganguly (2006). Spatial convergence of crack nucleation using a cohesive finite-element model on a pinwheel-based mesh. *International Journal for Numerical Methods in Engineering* 67(1):1-16. [DOI].

Pele, K., J. Baccou, L. Daridon, J. Liandrat, T. Le Gouic, Y. Monerie, and F. Péralès (2023). A probabilistic model for fast-to-evaluate 2D crack path prediction in heterogeneous materials. *Journal of Theoretical, Computational and Applied Mechanics* 15:1-21. [DOI], [OA].

Perales, F., S. Bourgeois, A. Chrysochoos, and Y. Monerie (2008). Two field multibody method for periodic homogenization in fracture mechanics of non linear composites. *Engineering*
# B Some hints on the range of validity for estimates (??) and (??)

The proposed estimates (??) and (??) are strictly devoted to unstructured Delaunay-like meshes. These estimates are not able to deal with:

* *Badly shaped elements.* Efficiency indexes lower than the minimal value $E_{PD}$ corresponding to PD tessellations (see ??) can not be taken into account by these estimates. They are restricted to “good” enough meshes (i.e. with an efficiency index greater than 0.725).

* *Structured or quasi-structured meshes.* Structured meshes lead to discrete values of the edge length and triangle top angle sets. Their distributions thus correspond to successive Dirac delta functions. Estimates (??) and (??) are restricted to “sufficiently unstructured” meshes (i.e. with an efficiency index at least less than 0.92).

To underline these two limitations, planar hand-generated meshes are considered as counterexamples. These meshes consist in cross-triangular square patterns stretched in one direction (degradation of the shape) with random perturbation of the inner points (in order to control the regularity), tessellations remaining valid. Three cases are considered.

<figure><img src="image_20.png" /><figcaption>**Figure 19:** Limit of validity for estimates (??) and (??) in case of quasi-structured meshes. Pdf of the edge length (first row) and triangle top angle (second row) sets for hand-generated cross-triangular square meshes, stretched in the horizontal direction and randomly perturbed of the inner points. Left: stretching/elongation ratio of 3 with high distortion, efficiency index $E \approx 0.78$, intensity $\tau = 466$. Middle: stretching ratio of 2 with intermediate distortion, efficiency index $E \approx 0.84$, intensity $\tau = 472$. Right: stretching ratio of 1 with small distortion, efficiency index $E \approx 0.87$, intensity $\tau = 458$.</figcaption></figure>

With highly random perturbation, the mesh is no more structured and estimates (??) and (??) exhibit a rough validity. In addition, these estimates do not accurately predict the tail of the edge length and tirangle angle distributions when badly shaped elements are numerous; see ??, left. With small random perturbation, the edge length distribution consists in multi-modal normal distribution. Estimates (??) and (??) are able to predict the mode with higher intensity, but are not able to predict the other modes; ??, middle. With even smaller random perturbation, a quasi-structured mesh is obtained and estimates (??) and (??) fail to predict the quasi-Dirac combs; ??, right.

Thus, the use of estimates (??) and (??) is restricted to standard random triangular tessellations whose efficiency index is greater than or equal to the one of a PD tessellation and at least less than 0.92 (see ??).
*Fracture Mechanics* 75(11):3378-3398. [DOI], [HAL].

Rebay, S. (1993). Efficient unstructured mesh generation by means of Delaunay triangulation and Bowyer-Watson algorithm. *Journal of Computational Physics* 106(1):125-138. [DOI].

Rimoli, J., J. Rojas, and F. Khemani (2012). On the mesh dependency of cohesive zone models for crack propagation analysis. *53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference & 20th AIAA/ASME/AHS Adaptive Structures Conference & 14th AIAA*. Honolulu, Hawaii: American Institute of Aeronautics and Astronautics. [DOI], [HAL].

Shewchuk, J. R. (1996). Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator. *Applied Computational Geometry Towards Geometric Engineering*. Ed. by M. C. Lin and D. Manocha. Lecture Notes in Computer Science. Springer, pp 203-222. [DOI].

Shewchuk, J. R. (2002). What is a Good Linear Element? Interpolation, Conditioning, and Quality Measures. *Proceedings of the 11th International Meshing Roundtable, IMR 2002* (Ithaca, New York, USA, Sept. 15-18, 2002). Ed. by N. Chrisochoides, pp 115-126. [HAL].

Sibson, R. (1980). The Dirichlet tessellation as an aid in data analysis. *Scandinavian Journal of Statistics* 7(1):14-20.

Tijssens, M. G., B. L. Sluys, and E. van der Giessen (2000). Numerical simulation of quasi-brittle fracture using damaging cohesive surfaces. *European Journal of Mechanics - A/Solids* 19(5):761-779. [DOI], [HDL].

Watanabe, D. (2008). Evaluating the configuration and the travel efficiency on proximity graphs as transportation networks. *Forma* 23(2):81-87. [OA].

Xu, X.-P. and A. Needleman (1994). Numerical simulations of fast crack growth in brittle solids. *Journal of the Mechanics and Physics of Solids* 42(9):1397-1434. [DOI].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long

as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

<figure><img src="image_26.png" /></figure>

**Authors' contributions** Author 1 is the major contributor to the paper writing. Author 1 has created most of the Figures in the paper and has provided all the results concerning the geodesic tortuosity of planar triangular meshes. Authors 2 and 3 have provided the results concerning the density of planar triangular meshes.

**Supplementary Material** All the results and figures presented in the paper are reproducible from a python script (workflows) available at the permalink 10.5281/zenodo.13833960. The dataset also contains the figures and tortuosity values (data).

**Acknowledgements** The authors are indebted to P. QUERRE for fruitful conversations and useful comments.

**Funding** This work was funded by the French Institute for Radiation Protection and Nuclear Safety (IRSN) and the University of Montpellier, France.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.
## C Complementary results for structured triangular meshes

The estimate of the total length of the edges of a tessellation containing $N_P$ points regularly distributed in an area A has been presented in Watanabe (2008) for various measures. Their results concern some triangular meshes, square meshes and hexagonal meshes. We focus on the mesh density Z of structured triangular meshes usually used in finite element discretizations (Tijssens et al. 2000).

The total length of edges in a structured triangular tessellation based on a hexagonal (H) point process reads (see ??):

$$ \sum_{i=1}^{N_E} l_i = 3 \sqrt{\frac{2\tau}{\sqrt{3}}} A. \qquad (60) $$

The corresponding mesh density is thus

$$ Z_H = 3 \sqrt{\frac{2}{\sqrt{3}}} \sqrt{\tau}. \qquad (61) $$

Regular points distributed along a square lattice (see ??) lead to a structured “Triangle Quadrilateral” (TQ) pattern whose the total length of the edges reads:

$$ \sum_{i=1}^{N_E} l_i = (2 + \sqrt{2}) \sqrt{\tau} A, \qquad (62) $$

and the mesh density

$$ Z_{TQ} = (2 + \sqrt{2}) \sqrt{\tau}. \qquad (63) $$

In fracture mechanics, “Cross-Triangle Quadrilateral” (CTQ) patterns are often used (see ??). For that mesh topology, the mesh density is deduced from the previous one:

$$ Z_{CTQ} = 2 (1 + \sqrt{2}) \sqrt{\tau}. \qquad (64) $$

This CTQ pattern can be streched in one direction and the underlying quadrilateral is not a square but a rectangle (aspect ratio $s \ge 1$). In this case, the mesh density becomes:

$$ Z_{CTQ_s} = \left( 2 \sqrt{1 + \frac{1}{s^2}} + 1 + \frac{1}{s} \right) \sqrt{\tau}. \qquad (65) $$

<table><caption>**Table 3:** Mesh density of the structured triangular meshes usually used in finite element discretizations</caption><tr><td><img src="image.png" alt="A triangular mesh with a square hole in its center."></td><td>$Z_H = 3 \sqrt{\frac{2}{\sqrt{3}}} \sqrt{\tau}$</td><td><img src="image.png" alt="A triangular mesh with a square hole in its center."></td><td>$Z_{TQ} = (2 + \sqrt{2}) \sqrt{\tau}$</td></tr><tr><td><img src="image.png" alt="A triangular mesh with a square hole in its center."></td><td>$Z_{CTQ} = 2 (1 + \sqrt{2}) \sqrt{\tau}$</td><td><img src="image.png" alt="A triangular mesh with a square hole in its center."></td><td>$Z_{CTQ_s} = (2 \sqrt{1 + \frac{1}{s^2}} + 1 + \frac{1}{s}) \sqrt{\tau}$</td></tr></table>

?? summarizes the mesh densities for structured meshes usually used in finite element discretizations. It is worth noting that these mesh densities for structured triangular meshes evolve as $\sqrt{\tau}$, exactly as the means of the mesh densities for random triangular meshes do. In particular, the mean $\mu_{St}(\tau; Z)$ of standard random triangular meshes given in ?? is very close to
<figure><img src="image_21.png" /><figcaption>**Figure 20:** The mesh density Z as a function of the point process intensity τ for different types of structured triangular meshes. Gray line: Cross-Triangle Quadrilateral mesh ?? . Gray dashed line: stretched Cross-Triangle Quadrilateral mesh ?? with stretch factor s = 2. Black dashed line: Triangle Quadrilateral mesh ?? . Black line: hexagonal mesh ?? .</figcaption></figure>

the mesh density of TQ tessellations: the prefactor of µ<sub>St</sub>(τ; Z) is about 3.40, while the prefactor of Z<sub>TQ</sub>(τ) is about 3.41:

$$
Z_{TQ}(\tau) \approx \mu_{St}(\tau; Z). \qquad (66)
$$

From the point of view of the mesh density, a Triangle Quadrilateral pattern accurately mimics a
standard random triangular mesh (keeping in mind that the distribution of the edge length set is
bimodal in TQ tessellations). Moreover, the mesh density Z of stretched CTQ (CTQs) decreases
as the stretch factor s increases, see ??.