# Variable resolution Poisson-disk sampling for meshing discrete fracture
networks

Johannes Krotz 11footnotemark: 1 Matthew R. Sweeney22footnotemark: 2 Carl W.
Gable33footnotemark: 3 Jeffrey D. Hyman44footnotemark: 4 Juan M.
Restrepo55footnotemark: 5

###### Abstract

We present the near-Maximal Algorithm for Poisson-disk Sampling (nMAPS) to
generate point distributions for variable resolution Delaunay triangular and
tetrahedral meshes in two and three-dimensions, respectively. nMAPS consists
of two principal stages. In the first stage, an initial point distribution is
produced using a cell-based rejection algorithm. In the second stage, holes in
the sample are detected using an efficient background grid and filled in to
obtain a near-maximal covering. Extensive testing shows that nMAPS generates a
variable resolution mesh in linear run time with the number of accepted
points. We demonstrate nMAPS capabilities by meshing three-dimensional
discrete fracture networks (DFN) and the surrounding volume. The discretized
boundaries of the fractures, which are represented as planar polygons, are
used as the seed of 2D-nMAPS to produce a conforming Delaunay triangulation.
The combined mesh of the DFN is used as the seed for 3D-nMAPS, which produces
conforming Delaunay tetrahedra surrounding the network. Under a set of
conditions that naturally arise in maximal Poisson-disk samples and are
satisfied by nMAPS, the two-dimensional Delaunay triangulations are guaranteed
to only have well-behaved triangular faces. While nMAPS does not provide
triangulation quality bounds in more than two dimensions, we found that low-
quality tetrahedra in 3D are infrequent, can be readily detected and removed,
and a high quality balanced mesh is produced.

###### keywords:

Discrete Fracture Network, Maximal Poisson-disk Sampling , Mesh Generation ,
Conforming Delauany Triangulation

††journal: Journal of Computational and Applied Mathematics

[label1]organization=Department of Mathematics, University of
Tennessee,city=Knoxville, postcode=37919, state=Tennessee, country=USA

[label2]organization=Computational Earth Science (EES-16), Earth and
Environmental Sciences, Los Alamos National Laboratory,city=Los Alamos,
postcode=87545, state=New Mexico, country=USA

[label3]organization=Oak Ridge National Laboratory,city=Oak Ridge,
postcode=37830, state=Tennessee, country=USA

##  1 Introduction

There are a number of methods used to model flow and the associated transport
of chemical species in low-permeability fractured rock, such as shale and
granite. The most common are continuum models, which use effective medium
parameters [21, 40, 54, 55, 66, 68] and discrete fracture network/matrix (DFN)
models, where fractures and the networks they form are explicitly represented
[10, 42, 56]. In the DFN methodology, individual fractures are represented as
planar $N-1$ dimensional objects embedded within an $N$ dimensional space.
Each fracture in the network is resolved with a computational mesh and the
governing equations for flow and transport are solved thereon. While the
explicit representation of fractures allows DFN models to represent a wider
range of transport phenomena and makes them a preferred choice when linking
network attributes to flow properties compared to continuum methods [25, 31,
29], it also leads to unique and complex issues associated with mesh
generation. Both conforming methods, where the mesh conforms to intersections
[30, 50, 51], and non-conforming methods, which use more complex
discretization schemes so the mesh does not need to be conforming [5, 18, 57,
58], have been developed. If the volume representing the rock matrix
surrounding the fracture network also needs to be meshed, then the
complications associated with mesh generation are compounded for both
techniques [4].

We present the near-Maximal Algorithm for Poisson-disk Sampling (nMAPS) that
generates a conforming variable resolution triangulation mesh on and around
three-dimensional discrete fracture networks. Although the theory of Poisson
disk sampling to produce a point distribution is straightforward, details of
the implementation, algorithm, termination criteria, and maximizing efficiency
are nuanced and often problem-dependent. nMAPS addresses these issues in the
context of mesh generation using a two-stage approach. The first stage is
based on the framework presented by Dwork et al., [14] and uses a rejection
algorithm to generate an initial Poisson-disk sample point set. The second
phase is based on the framework presented by Mitchell et al., [47] and adds
additional points to fill gaps in the covering, which maximizes point density
without violating the restrictions of a Poisson-disk sampling. nMAPS combines
the two methods to achieve near-maximal coverage with a run time that scales
linearly in the number of accepted points. nMAPS efficiency results from a
novel rejection technique that uses a quadrilateral or hexahedral background
mesh search based on already sampled points. This method significantly
shortens the time needed to reach a sufficient, i.e., near-maximal, point
density due to a large reduction in computational operations. Once the point
set is generated, the conforming Delaunay algorithm presented by Murphy et
al., [49] is used to generate a Delaunay triangulation in two or three
dimensions. Although nMAPS does not guarantee that the point distribution’s
density is maximal, i.e., no additional points can be added without violating
the restrictions on distances between points, the sample sets are sufficiently
maximal to produce high-quality meshes. Comparison with previously implemented
meshing techniques for variable resolution conforming triangulations showed
that nMAPS produces meshes comparable in quality but requires considerably
less computational time. Moreover, the nMAPS framework is presented in a
general manner that can be easily extended beyond mesh generation for DFNs.

In section 2, we describe the challenges associated with DFN mesh generation
and the general properties of maximal Poisson-disk sampling. In section 3, we
provide a detailed explanation of nMAPS, for both 2D fracture networks and 3D
volume meshing. In section 4, we assess the quality of the mesh and
performance, i.e., run times, of 2D and 3D examples using nMAPS. In section 5,
we provide a few concluding remarks.

##  2 Background

###  2.1 Discrete Fracture Networks & Mesh Generation

Due to the epistemic uncertainty associated with hydraulic and structural
properties of subsurface fractured media, fracture network models are
typically constructed stochastically [52, 53, 55]. In the DFN methodology,
individual fractures are placed into the computational domain with locations,
sizes, and orientations that are sampled from appropriate distributions based
on field site characterizations. The fractures form an interconnected network
embedded within the porous medium. Each fracture must be meshed for
computation so that the governing equations for flow and transport can be
numerically integrated to simulate physical phenomena of interest.

Formally, each fracture in a DFN can be represented as a planar straight-line
graph (PSLG) composed of a set of line segments that represent the boundary of
the fracture and a set of line segments that represent where other fractures
intersect it. In this manner, each fracture can be described by a set of
boundary points on the PSLG, denoted $\\{p\\}$, and a set of intersection
lines $\\{\ell_{i,j}\\}$, where the subscripts $i$ and $j$ indicate that this
line corresponds to the intersection between the $i$th and $j$th fractures.
Once $\\{p\\}$ and $\\{\ell_{i,j}\\}$ are obtained for every fracture in the
network, a point distribution covering each fracture must be generated. If a
conforming numerical scheme is used, then all cells of $\\{\ell_{i,j}\\}$ are
discretized lines in the mesh that must coincide between intersecting
fractures. So long as minimum feature size constraints are met, a conforming
triangulation method, such as that presented in Murphy et al., [49], can be
implemented to connect the vertices such that all lines of intersection form a
set of connected edges in a triangulation.

In general, one wants to properly resolve all relevant flow and transport
properties of interest while minimizing the number of nodes in the mesh, and
these two goals compete. The first of these conditions depends partially on
the mesh quality, which is controlled partially by the second. A starting
point for the notion of mesh quality is the minimum angle condition, i.e.,
that the smallest angle should be bounded away from zero, as poorly shaped
elements can affect the condition number of the linearized system of equations
[1, 63, 69]. While a uniform resolution mesh is straightforward to generate
and the minimum angle condition easily met, it is computationally more
expensive than a variable resolution mesh refined in areas of interest, which
can be tailored to reduce the number of nodes in the mesh. A variable
resolution mesh can be appropriate for single-phase flow simulations or in
particle tracking simulations where the spatially variable resolution does not
adversely affect transport properties. Specifically, Eulerian formulations of
transport will have spatially variable numerical diffusion on a variable
resolution mesh, which means that the mesh needs to have a sufficiently small
refinement to accurately capture solute fronts [2]. Variable resolution mesh
generation is more complex than its uniform counterpart. One of the principal
complications is creating a smooth transition of cell sizes. Jumps in the
computed fields of interest and other numerical artifacts can occur if the
transition is not sufficiently smooth.

Most methods to generate a conforming DFN mesh use a uniform point
distribution and modify connectivity locally to conform to intersections [50,
51]. When using a conforming mesh, the numerical methods for resolving flow
and transport in the network are typically simpler and have fewer degrees of
freedom compared to non-conforming mesh methods [20]. Similarly, almost all
non-conforming numerical methods use a uniform resolution, but some create
variable resolutions across fractures, still uniform within a single plane, in
an attempt to reduce the number of total nodes in the mesh [6]. Using a
variable mesh resolution in non-conforming schemes could drastically reduce
the number of nodes in the mesh while retaining the the ability to retain a
high order of accuracy. However, it is rarely implemented due to the
associated meshing complications [8].

The generation of a variable resolution, unstructured conforming mesh is quite
rare, even with the advantages noted above. One technique in use is the
Features Rejection Algorithm for Meshing (FRAM) that addressed the issues
associated with conforming DFN mesh creation by coupling it with network
generation [30]. Through this technique, FRAM allows for the creation of a
variable resolution mesh that smoothly coarsens away from intersections where
pressure gradients in flow simulations are typically the highest. FRAM has
been implemented in the computational suite dfnWorks [33], which has been used
to probe fundamental aspects of geophysical flows and transport in fractured
media [27, 28, 35, 36, 44, 62] as well as practical applications including
hydraulic fracturing operations [26, 37, 43], inversion of micro-seismicity
data for characterization of fracture properties [48], the long term storage
of spent civilian nuclear fuel [25], and geo-sequestration of carbon dioxide
into depleted reservoirs [32]. However, the current method for mesh generation
uses an inefficient iterative refinement technique to produce the point
distribution used to generate the mesh. Sweeps of a refinement algorithm are
applied to an initially coarse triangulation based on the boundary set
$\\{p\\}$. If an edge in the mesh is larger than the current maximum edge
length, then a new point is added to the mesh at the midpoint of that edge to
split it into two new edges. In practice, the edge splitting is done using
Rivara refinement [60, 61]. The resulting field is then smoothed using
Laplacian smoothing in combination with Lawson flipping [38]. This process is
repeated until all edges meet the assigned target edge length, which could be
a spatially variable field based on the distance to $\\{\ell_{i,j}\\}$, for
example. While the resulting mesh quality is quite good, the process is rather
slow and cumbersome.

If these implementation complexities can be addressed, then the superior
modeling qualities of variable resolutions can be made practical. To
accomplish this, we present nMAPS, where the final vertex distribution is
directly created rather than iteratively derived. While the method was
initially designed to specifically improve FRAM, we provide the details in a
general format such that it can be implemented for mesh generation in general.
Specifically, it can be used as the basis for any DFN meshing methodology,
including both conforming and non-conforming techniques. Details are given for
Delaunay triangulations, which are of importance in many two-point flux finite
volume solvers as they are used to generate the Voronoi control volumes on
which these solvers compute. In the next section, we recount the properties of
maximal Poisson-disk sampling that we used to design and implement nMAPS,
including theoretical bounds on mesh gradation that ensure high-quality
variable mesh resolutions.

###  2.2 Maximal Poisson-disk Sampling

Meshes produced from a point distribution that is dense yet cluster free have
provable high quality bounds under specific conditions [7, 12, 16]. Maximal
Poisson-disk samples fulfill these conditions [22]. Similar quality bounds can
be established for sphere-packings, whose radii are Lipschitz continuous with
respect to their location [45, 46, 64]. Traditionally, Poisson-disk sampling
is performed using an expensive dart-throwing algorithm [13]. These algorithms
struggle to achieve maximality as the probability to select a free spot
becomes decreasingly small as the number of points sampled ($n$) increases.
The algorithm in [47], which is based on these dart-throwing algorithms, was
the first to guarantee maximality. It achieves maximality with run times of
$O(n\log(n))$ by using a regular grid for acceleration and sampling from
polygonal regions in its second phase. In practice, performance close to
$O(n)$ has been reported [15, 17, 47]. Another algorithm that is not based on
dart-throwing was proposed in [9]. While the method does not guarantee
maximality, it did show linear performance in the number of points sampled.
This algorithm was extended to variable radii by [14]. Other authors have
provided algorithms that produce variable Poisson-disk samplings on
3D-surfaces [23, 24]. Their triangulation-based algorithm runs in linear time,
and while theoretically not guaranteeing maximality, their numerical
experiments indicate that near-maximality can certainly be achieved. A summary
of recent developments in this area, along with a comparison of different
methods, can be found in [67].

Maximal Poisson-disk samplings $X$ on a domain $\Omega\subseteq\mathbb{R}^{d}$
are random selections of points $X=\\{{\mathbf{x}}_{i}\\}_{i=1}^{n}$, that
fulfill the following properties:

  1. 1.

empty disk property:

$\forall i\neq j\in\\{1,...,n\\}:|{\mathbf{x}}_{i}-{\mathbf{x}}_{j}|>r.$

We will call $r$ the inhibition radius,

  2. 2.

maximality:

$\Omega=\bigcup_{i=1}^{n}B_{R}({\mathbf{x}}_{i}),$

where
$B_{\varepsilon}({\mathbf{x}})=\\{{\mathbf{y}}\in\Omega:|{\mathbf{x}}-{\mathbf{y}}|<{\mathbf{y}}\\}$
is the open ball of radius $\varepsilon$ around ${\mathbf{x}}$. $R$ will be
called the coverage radius. [47]

Intuitively, the empty disk property says that every sample point is at the
center of a $d$-dimensional ball or disk that does not contain any other
points of the sampling. Maximality implies that these balls cover the whole
domain, i.e., there is no point $y\in\Omega$ that is not already contained in
one of the balls around a point in the sample.

It is useful to generalize these definitions, such that both the inhibition
and the coverage radius depend on the sampling points, i.e.,
$r=r({\mathbf{x}}_{i},{\mathbf{x}}_{j})$ and
$R=R({\mathbf{x}}_{i},{\mathbf{x}}_{j})$ for all
${\mathbf{x}}_{i},{\mathbf{x}}_{j}\in X$. We hereafter refer to this construct
as a variable radius maximal Poisson-disk sampling, and we refer to a Poisson-
disk sampling with constant radii as a fixed-radii maximal Poisson-disk
sampling [47].

A common approach is to assign each point ${\mathbf{x}}\in\Omega$ a positive
radius $\rho({\mathbf{x}})$ and have $r({\mathbf{x}}_{i},{\mathbf{x}}_{j})$ be
a function of $\rho({\mathbf{x}}_{i})$ and $\rho({\mathbf{x}}_{j})$. Natural
choices for $r({\mathbf{x}}_{i},{\mathbf{x}}_{j})$ are, for example,
$\rho({\mathbf{x}}_{i})$ or $\rho({\mathbf{x}}_{j})$ for $i<j$, thereby
determining the inhibition radius depending on the ordering on $X$. Order
independent options include
$\min(\rho({\mathbf{x}}_{i}),\rho({\mathbf{x}}_{j})),\max(\rho({\mathbf{x}}_{i}),\rho({\mathbf{x}}_{j}))$
or $\rho({\mathbf{x}}_{i})+\rho({\mathbf{x}}_{j})$. The last of these options
corresponds to a sphere packing [47]. The coverage radius can but does not
have to be different from $\rho$.

Delaunay triangulations maximize the smallest angle of all triangulations
based on a generating point distribution [41]. Since numerical errors in many
applications tend to increase if these angles are small [69], Delaunay
triangulations often are a triangulation of choice. Moreover, the dual of the
Delaunay triangulation is a Voronoi tessellation, which in a certain sense is
optimal for two-point flux finite volume solvers [19] that are commonly used
in subsurface flow and transport simulators such as fehm [70], tough2 [59],
and pflotran [39]. In the case of maximal Poisson-disk samplings, we can go
one step further and give a lower bound on these angles. We provide a brief
summary of the proofs found in [47] and highlight the most important results
we use. We then proceed with the new bounds.

######  Lemma 1.

The smallest angle $\alpha$ in any triangle is greater than
$\arcsin\left(\frac{r}{2R}\right)$, where $r$ is the length of the shortest
edge and $R$ the radius of the circumcircle or

$\displaystyle\sin(\alpha)\geq\frac{r}{2R}$ (1)

###### Proof.

This is a direct corollary of the central angle theorem. ∎

This Lemma allows us to give explicit bounds for maximal Poisson-disk
samplings. While we will focus entirely on inhibition radii given by
$r({\mathbf{x}}_{i},{\mathbf{x}}_{j})=\min(\rho({\mathbf{x}}_{i}),\rho({\mathbf{x}}_{j}))$,
where $\rho({\mathbf{x}})$ is some positive function, comparable results can
be found for different $r({\mathbf{x}}_{i},{\mathbf{x}}_{j})$ in a similar
fashion.

######  Lemma 2.

Let $\varepsilon\geq 0$ and $\rho:\mathbb{R}^{d}\rightarrow\mathbb{R}$ ($d\geq
2$) be a positive Lipschitz continuous function with Lipschitz constant $L$
with $L\varepsilon<1$. Let $X$ be a variable maximal Poisson-disk sampling on
the domain $\Omega\subset\mathbb{R}^{d}$ with inhibition radius
$r({\mathbf{x}},{\mathbf{y}})=\min(\rho({\mathbf{x}}),\rho({\mathbf{y}}))$ and
coverage radius
$R({\mathbf{x}},{\mathbf{y}})\leq(1+\varepsilon)r({\mathbf{x}},{\mathbf{y}})$
with $\varepsilon>0$. Let the triangle $\Delta$ be an arbitrary element of the
Delaunay triangulation of $X$ ($d=2$) or an arbitrary 2-dimensional face of a
cell of the Delaunay triangulation of $X$ ($d>2$). If the circumcenter of
$\Delta$ is contained in $\Omega$, then each angle $\alpha$ of $\Delta$ is
greater or equal to $\arcsin\left(\frac{1-L-\varepsilon
L}{2+2\varepsilon}\right)$ or

$\displaystyle\sin(\alpha)\geq\frac{1-L-\varepsilon L}{2+2\varepsilon}.$

###### Proof.

Let $\alpha$ be the smallest angle of $\Delta$ and
${\mathbf{x}},{\mathbf{y}}\in X$ be the vertices of the shortest edge of
$\Delta$, i.e., the vertices opposite to $\alpha$. Without loss of generality,
assume $\rho({\mathbf{x}})\leq\rho({\mathbf{y}})$. Since $X$ is a Poisson-disk
sampling
$|{\mathbf{x}}-{\mathbf{y}}|\geq\min\left(\rho({\mathbf{x}}),\rho({\mathbf{y}})\right)=\rho({\mathbf{x}})$.
Now, let ${\mathbf{z}}\in\Omega$ be the circumcenter of $\Delta$. Since $X$ is
maximal, there exists ${\mathbf{v}}\in X$ with
$|{\mathbf{z}}-{\mathbf{v}}|\leq
R({\mathbf{z}},{\mathbf{v}})\leq(1+\varepsilon)\rho({\mathbf{z}})$. Because
$\Delta$ was retrieved from a Delaunay triangulation ${\mathbf{v}}$ cannot be
contained in the interior of $\Delta$’s circumcircle. Hence,

$\displaystyle|{\mathbf{z}}-{\mathbf{x}}|\leq|{\mathbf{z}}-{\mathbf{v}}|\leq(1+\varepsilon)\rho({\mathbf{z}})\leq(1+\varepsilon)\left(\rho({\mathbf{x}})+L|{\mathbf{z}}-{\mathbf{x}}|\right).$

Rearranging this inequality yields

$\displaystyle|{\mathbf{z}}-{\mathbf{x}}|\leq\rho({\mathbf{x}})\frac{1+\varepsilon}{1-L-\varepsilon
L}.$

The result follows by applying Lemma 1 after noticing that
$|{\mathbf{x}}-{\mathbf{y}}|$ is the length of the shortest edge and that
$|{\mathbf{z}}-{\mathbf{x}}|$ is the radius of the circumcirle. ∎

######  Remark.

Note that for $d>2$ the same result is true, if we assume the circumcenter of
the $d$-simplex, of which $\Delta$ is a face, is contained in $\Omega$ instead
of the circumcenter of $\Delta$ itself. The proof is identical.

######  Remark.

While this result allows us to control the quality of 2D-triangulations of
maximal Poisson-disk samplings, it can also be used to gauge how close a given
Poisson-disk sampling is to being maximal.

The previous Lemma only gives us bounds on all triangles, if their
circumcenters are contained in $\Omega$. The next two Lemmas will give
sufficient conditions to guarantee exactly this as long as $\Omega$ is a
polytope.

######  Lemma 3.

Let $\Omega\subset\mathbb{R}^{2}$ be a polygonal region and $X$ a maximal
Poisson-disk sampling containing all vertices of $\Omega$. Let the inhibition
radius $r({\mathbf{x}},{\mathbf{y}})$ be defined like in the previous lemma.
Furthermore, let the coverage radius of $X\cap\delta\Omega$ fulfill
$R^{\delta}({\mathbf{x}},{\mathbf{y}})<\frac{r({\mathbf{x}},{\mathbf{y}})}{\sqrt{2}(1+L)}$,
i.e.,
$|{\mathbf{x}}-{\mathbf{y}}|<{\frac{\sqrt{2}}{1+L}}r({\mathbf{x}},{\mathbf{y}})$
for all ${\mathbf{x}},{\mathbf{y}}\in\delta\Omega\cap X$. Then the
circumcenter of all triangles in the Delaunay triangulation of $X$ are
contained in $\bar{\Omega}$.

###### Proof.

Suppose this claim is wrong. Let $\Delta$ be a triangle in the Delaunay
triangulation with circumcenter ${\mathbf{z}}\notin\Omega$. For this to be
possible the circumcircle needs to be cut into at least two pieces by
$\delta\Omega$, separating ${\mathbf{z}}$ and the vertices of $\Delta$. Since
$\Delta$ is part of a Delaunay triangulation and all vertices of $\Omega$ are
part of the sampling, this is done by (at least) one segment of a straight
line, i.e., $\delta\Omega$ contains a secant of the circumcircle.

Next, let ${\mathbf{b}}_{1},{\mathbf{b}}_{2}\in\delta\Omega$ be the two
boundary points closest to the circumcircle on either side of that line
segment and let $B$ be the disk bounded by the circumcircle. Note that
$\bar{B}\cap\Omega$ contains $\Delta$ and is itself entirely contained in the
disk of radius
$\frac{1}{2}|{\mathbf{b}}_{1}-{\mathbf{b}}_{2}|<\frac{1}{2}{\frac{\sqrt{2}}{1+L}}r({\mathbf{b}}_{1},{\mathbf{b}}_{2})$
around $\frac{1}{2}({\mathbf{b}}_{1}+{\mathbf{b}}_{2})$.

Now, let ${\mathbf{x}}\notin\\{{\mathbf{b}}_{1},{\mathbf{b}}_{2}\\}$ be a
vertex of $\Delta$ and let
${\mathbf{b}}\in\\{{\mathbf{b}}_{1},{\mathbf{b}}_{2}\\}$ be the point of the
two, that is closer to ${\mathbf{x}}$. We already established that
${\mathbf{x}}$ lies within the just mentioned ball around
$\frac{1}{2}({\mathbf{b}}_{1}+{\mathbf{b}}_{2})$. Let ${\mathbf{x}}_{p}$ be
the projection of ${\mathbf{x}}$ onto the line segment connecting
${\mathbf{b}}_{1}$ and ${\mathbf{b}}_{2}$. Then
$|{\mathbf{x}}-{\mathbf{x}}_{p}|<\frac{1}{2}{\frac{\sqrt{2}}{1+L}}r({\mathbf{b}}_{1},{\mathbf{b}}_{2})$,
because ${\mathbf{x}}$ lies within the circle of that radius,
$|{\mathbf{x}}_{p}-{\mathbf{b}}|<\frac{1}{2}{\frac{\sqrt{2}}{1+L}}r({\mathbf{b}}_{1},{\mathbf{b}}_{2})$,
because ${\mathbf{b}}$ is the closer of the two points
${{\mathbf{b}}_{1},{\mathbf{b}}_{2}}$ and therefore

$\displaystyle|{\mathbf{x}}-{\mathbf{b}}|=\sqrt{|{\mathbf{x}}-{\mathbf{x}}_{p}|^{2}+|{\mathbf{b}}-{\mathbf{x}}_{p}|^{2}}<\frac{r({\mathbf{b}}_{1},{\mathbf{b}}_{2})}{1+L}\leq\frac{\rho({\mathbf{b}})}{1+L}\leq\rho({\mathbf{b}}).$
(2)

Since
$|{\mathbf{x}}-{\mathbf{b}}|\geq\min(\rho({\mathbf{x}}),\rho({\mathbf{b}}))$
this implies $|{\mathbf{x}}-{\mathbf{b}}|\geq\rho({\mathbf{x}})$. However
assuming this and applying the Lipschitz condition on (2) gives us

$\displaystyle|{\mathbf{x}}-{\mathbf{b}}|<\frac{\rho({\mathbf{b}})}{1+L}\leq\frac{1}{1+L}\left(\rho({\mathbf{x}})+L|{\mathbf{x}}-{\mathbf{b}}|\right)\leq\frac{1+L}{1+L}|{\mathbf{x}}-{\mathbf{b}}|,$

which is a contradiction. ∎

######  Remark.

Although Lemma 3 generalizes to higher dimensions, it is not very practical
because it is difficult to guarantee the bounds on $R^{\delta}$ if the
boundary is more than 1-dimensional. However, it is still possible to get some
bounds on the radii of the circumcircles and then, using Lemma 1, on the
angles, if the distance of non-boundary points is greater than some lower
bound $a>0$. In fact, using notation from the previous proof, let $\Delta$
again be an $d$-simplex with circumcenter outside of $\Omega$ and
${\mathbf{x}}\notin\delta\Omega$ on of its points. Since the circumsphere of
any simplex in a Delaunay triangulation does not contain any other points the
radius of the intersection with $\delta\Omega$ is bounded by $R^{\delta}$.
Using simple geometric arguments, one can show that this feature forces the
radius of $\Delta$’s circumsphere $R$ to fulfill the following inequality

$\displaystyle R^{2}\leq(R-a)^{2}+\left(R^{\delta}\right)^{2}\Rightarrow
R\leq\frac{a^{2}+\left(R^{\delta}\right)^{2}}{2a}\leq\frac{\left(R^{\delta}\right)^{2}}{a}.$
(3)

If $R^{\delta}({\mathbf{x}},{\mathbf{y}})<r({\mathbf{x}},{\mathbf{y}})$ there
is a lower bound on $a$, continuously depending on $R^{\delta}$, solely due to
the fact, that we have a Poisson-disk sampling. If
$R^{\delta}=R^{\delta}({\mathbf{x}}_{p},{\mathbf{b}})$ is any bigger, $a$
needs to be bounded artificially. This implies that the angle bounds change
continuously, if the conditions for Lemma 3 cannot be met they still can be
relatively controlled by the choice of the artificial bound on $a$.

Under the conditions of the previous Lemmas, the simplices of the two-
dimensional Delaunay triangulation resulting from the maximal Poisson-disk
sampling are guaranteed to only have well-behaved triangular faces. In three
or more dimensions, however, this does not imply that the simplices themselves
are well-behaved. It is possible for a 3D Delaunay triangulation to contain
slivers, which are tetrahedra whose four points are all positioned
approximately on the equator of their circumsphere. In [11] slivers are
characterized as tetrahedra whose points are all close to a plane and whose
orthogonal projection onto that plane is a quadrilateral. In [3] slivers are
equivalently classified as tetrahedra with a dihedral angle close to
$180^{\circ}$ containing their own circumcenter. All faces of a sliver can be
equilateral triangles, yet they can still have arbitrarily small dihedral
angles. This latter feature can cause large numerical errors in the
computational physics simulations performed on the mesh. Therefore, slivers
should be eschewed as much as possible. While slivers cannot be entirely
avoided, one can show that if the points
${\mathbf{x}},{\mathbf{y}},{\mathbf{z}},{\mathbf{w}}$ of a maximal Poisson-
disk sampling form a sliver, the distance between ${\mathbf{w}}$ and the plane
spanned by ${\mathbf{x}},{\mathbf{y}},{\mathbf{z}}$ is very small [11]. We use
this property to minimize slivers around the DFN and the faces of the
surrounding volume by using the 2D sampling of the DFN as the seed of the 3D
sampling and enforce a minimum distance between the DFN and subsequent points
in the 3D sample. This procedure leads to slivers being scarce in the initial
3D triangulation. Moreover, given any three points in the mesh, the
probability a fourth point will produce is a sliver quite small in practice.
Therefore, if the points producing slivers are removed and resampled, the
resulting new triangulation will most likely have fewer slivers than the
previous one. We tested this rejection/resampling algorithm and found it
improves the overall quality of the mesh.

##  3 Method: The near-Maximal Algorithm for Poisson-disk Sampling (nMAPS)

###  3.1 nMAPS Workflow overview

The nMAPS workflow contains the following high-level steps:

  1. 1.

Generation of a DFN using dfnWorks [34].

  2. 2.

Deconstruct the DFN into individual fractures/polygons.

  3. 3.

Generation of a 2D-variable-radii Poisson-disk sampling on each polygon.

  4. 4.

Construct a conforming Delaunay triangulation of each fracture based on the
Poisson-disk sample point set.

  5. 5.

Merge individual fracture meshes into a single DFN.

  6. 6.

Generation of a 3D variable radii Poisson-disk sampling of the volume
surrounding the meshed DFN.

  7. 7.

Mesh generation of the resulting 3D sample set, identify low-quality
tetrahedra and remove two of their nodes that are not located in the DFN mesh.

  8. 8.

Repeat steps 6-7 with the remaining nodes as the seed set until all slivers
are removed.

A graphical version of the workflow is presented in Figure 1. DFN generation
was discussed in section 2.1, so we begin with a description of the 2D
sampling method.

1\. Deconstruct the DFN into individual fractures 2\. Sample points each
fracture independently (A-D) 3\. Re-assemble the DFN and eliminate
conflicts4\. Sample points in the surrounding volume 5\. Recursively remove
slivers

A. Identify boundary and intersections B. Sample points along boundaryC. Use
boundary points as the seed for nMAPS D. When the algorithm terminates, find
the undersampled regions and sample new points therein. Repeat from step C.

Figure 1:  Overview of the nMAPS workflow from the creation of DFN to the
final mesh (left 1-5) and during 2D-sampling (right A-D).

###  3.2 Two-Dimensional Sampling Method

We generate a 2D Poisson-disk sampling in a successive manner using a
rejection method. First, we provide a general overview of the method. This
method can be performed on every fracture in the network independent of the
other fractures. In each step, a new point candidate is generated. It is
accepted if it does not break the empty disk property with any of the already
accepted points. For the sampling in two dimensions, we use a variable
inhibition radius that increases linearly with distance from intersections on
the fracture. Specifically, we reject a candidate point ${\mathbf{y}}$, if
there is an already excepted point ${\mathbf{x}}$ such that the condition

$\displaystyle|{\mathbf{x}}-{\mathbf{y}}|\geq
r({\mathbf{x}},{\mathbf{y}})=\min(\rho({\mathbf{x}}),\rho({\mathbf{y}}))$ (4)

is violated. In this equation, $\rho({\mathbf{x}})$ as a piecewise linear
function given by

$\displaystyle\rho({\mathbf{x}})=\rho(D({\mathbf{x}}))=\left\\{\begin{matrix}\frac{H}{2}&\text{for}&D({\mathbf{x}})\leq
FH\\\ A(D({\mathbf{x}})-FH)+\frac{H}{2}&\text{for}&FH\leq
D({\mathbf{x}})\leq(R+F)H\\\
(AR+\frac{1}{2})H&&\text{otherwise}\end{matrix}\right.$ (5)

Here $D({\mathbf{x}})$ is the Euclidean distance between ${\mathbf{x}}$ and
the closest intersection. $H,A,R$, and $F$ are parameters that determine the
global minimum distance between two points ($H/2$), the distance around an
intersection where the local inhibition radius remains at its minimum value
($FH$), the global maximum inhibition radius ($ARH+H/2$), and the slope at
which the inhibition radius grow with $D({\mathbf{x}})$ ($A$). Since $\rho(D)$
is piecewise linear, it is a Lipschitz-function with Lipschitz-constant $A$.

If the sampling has a coverage radius
$R({\mathbf{x}},{\mathbf{y}})\leq(1+\varepsilon)r({\mathbf{x}},{\mathbf{y}})$
for some $\varepsilon>0$, then the conditions of (2) hold. To satisfy the
conditions of (3) as well, and thereby ensure angle bounds on all triangles in
a Delaunay triangulation, we begin by sampling points along the boundary and
enforce a maximum distance of
$\frac{r({\mathbf{x}},{\mathbf{y}})}{\sqrt{2}(1+L)}$ between points. Next, we
generate new candidates by randomly sampling within an annulus around an
already accepted point; this is illustrated in Figure 2. The minimum distance
another point could be from the center point and still preserve the empty disk
property determines the inner radius. The maximum distance a point be from to
the center in a maximal sampling determines the outer radius. For our choice
of inhibition radius, and assuming it has the same radius as the coverage
radius, these distances are $\frac{\rho({\mathbf{x}})}{1+A}$ for the inner
radius and $\frac{2\rho({\mathbf{x}})}{1-A}$ for the outer radius.

Figure 2: Visualisation of a single sampling step. The current point is in the
center, the new candidates in the annulus ($k=4$). The inner circle is bounded
by the inhibition radius of the current point. The outer circle is bounded by
the maximum distance a point can be away from the current point if the
Poisson-disk sampling was maximal.

We will now go over the individual steps of the 2D algorithm in detail. These
steps are also presented in the pseudocode in Algorithm 1 (Section 3.4) and
are illustrated in A-D of Fig. 1. The notation used in the pseudocode is found
in the table at the start of the same section. First, a 1D Poisson-disk
sampling along the boundary of the polygon is generated as a seed set (line
2). Next, $k$ candidate points at a time (line 12) around each already
accepted point are sampled and determined whether they are accepted or not
(lines 13 through 21). $k$ is a positive integer and a user-defined parameter.
We use cell lists to find points around a candidate that could potentially
cause a violation of the empty disk. A visual depiction of these cells is
provided in Fig. 3(a). The size of these cells is chosen to contain at most
one point, which allows us to skip distance calculations with points beyond a
certain cutoff (line 16). Unlike the previously mentioned algorithms, we label
cells containing points as occupied and those cells that are too close to an
accepted point to contain the candidate point. Specifically, if a candidate
${\mathbf{x}}$ lies in a cell $C$ and any other cell $D$ with $diam(C\cup
D)\leq r_{in}({\mathbf{x}})$ is occupied, then ${\mathbf{x}}$ is immediately
rejected as it conflicts with the point in $D$ (line 13). On the other hand,
if $dist(C,D)>\rho({\mathbf{x}})$, then there is no need to calculate the
distance between ${\mathbf{x}}$ and any potential element of $D$ because they
cannot violate the empty disk property. An example of this property is shown
in Fig. 3(b). We use this technique to our advantage in two ways. First, it
allows us to reject many candidates without calculating any distances to
nearby points, which provides substantial speedup compared to previous methods
and especially for large values of $k$. Second, unmarked cells that contain
space for another point are easily identified, which allows us to find under-
sampled regions after the initial sweep terminates (line 32). If a point is
accepted, then it is added to the end of the sample set (line 22). Note that
this newly accepted point will be used as a seed for sampling another $k$
point when its turn comes up in the queue. Thus, the method grows the sample
set within the primary while loop. If all $k$ candidates around a point are
rejected, then we move on to the next already accepted point (line 29). The
algorithm terminates once every accepted point has been used as a sampling
center (line 30). Next, we detect unmarked or under-sampled cells and randomly
place points within them (line 33). Then, the main algorithm is restarted with
the seed set, including these newly added points, and ends once the
termination criterion is met once again (line 46). While this process of
resampling sweeps in under-sampled cells can be performed multiple times, we
found that a single resampling was sufficient to increase the quality of the
sampling to acceptable levels, and additional sweeps were unnecessary.

Once the point distribution is created, the conforming Delaunay triangulation
method of Murphy et al., [49] is used to generate the mesh. We recount the
general idea of the method here for completeness. To create a conforming
Delaunay triangulation that preserves the lines of fracture intersections as a
set of triangle edges to be created, it is sufficient that the circumscribed
circle of each segment of the discretized line of intersection be empty of any
other point in the distribution prior to connecting the mesh. To achieve this
sufficient condition, any point within the circumscribed circle of each
segment of the discretized lines of intersection is removed from the point
distribution. Next, a two-dimensional unconstrained Delaunay triangulation
algorithm is used to connect the resulting point set. Because of the
construction method, i.e., empty regions around the lines of intersection, the
line segments that represent lines of fracture intersection naturally emerge
in the triangulation. In turn, the Delaunay triangulation will conform to all
of the fracture intersection line segments. Once every fracture polygon is
triangulated, they are all joined together to form the mesh of the entire DFN.

(a)  (b)

Figure 3: Visualization of how the grid is used to find possibly conflicting
points. The new candidate is red, already accepted points are green, and cells
containing conflicting points are shaded grey. The red circle shows the
inhibition radius of the candidate, and blue circles show the furthest cells a
point in the center cell could conflict with.

###  3.3 Three-dimensional Sampling Method

nMAPS uses a similar method to that presented in the two-dimensional section
to generate a point distribution in three-dimensional space. The primary
difference is that candidates are generated on a spherical shell around
accepted nodes instead of an annulus. The 3D variant of $\rho({\mathbf{x}})$
given by

$\displaystyle\rho({\mathbf{x}})=\rho(D({\mathbf{x}}))=\left\\{\begin{matrix}\rho_{2}({\mathbf{x}}_{p})&\text{
for }&D({\mathbf{x}})\leq F\rho_{2}({\mathbf{x}}_{p})\\\ \\\
A(D({\mathbf{x}})-Fr_{2D}({\mathbf{x}}_{p}))+\frac{H}{2}&\text{ for
}&F\rho_{2D}({\mathbf{x}}_{p})\leq
D({\mathbf{x}})\leq\frac{\rho_{max}-r_{2}({\mathbf{x}}_{p})}{A}\\\ \\\
\rho_{max}&&\text{otherwise}\end{matrix}\right.$ (6)

where ${\mathbf{x}}_{p}$ the fracture point closest to ${\mathbf{x}}$,
$\rho_{2}({\mathbf{x}}_{p})$ is its 2D inhibition radius on the fracture, and
$D({\mathbf{x}})$ is the distance between ${\mathbf{x}}$ and
${\mathbf{x}}_{p}$. Just as in the 2D case, equation (6) is a piecewise linear
function in $D({\mathbf{x}})$. It is constant within a distance of
$\rho_{2}({\mathbf{x}}_{p})F$ of the fracture network and increases linearly
with a slope of $A$ until a given maximal inhibition radius of $\rho_{max}$ is
reached. In practice these parameters can be chosen to be the same as their 2D
counterparts. One difference from the 2D case is that in addition to rejecting
all candidates that violate the empty-disk property (4), we also reject a
candidate ${\mathbf{x}}$ if it is within a distance of $\rho({\mathbf{x}})/2$
to a boundary or fracture. This latter piece prevents slivers from having
three nodes located on a single fracture or the domain’s boundary. In turn,
this limits the circumradius of tetrahedra with circumcenter outside of the
domain; cf. lemma 3 and subsequent remark for additional details.

The pseudocode of the 3D-sampling method provided is in Algorithm 2 in Section
3.4. The initial sampling process is identical to the 2D version, therefore,
we will focus on the initialization and resampling, which is where the methods
differ. The seed set is made up of points sampled along the boundary of the 3D
volume and those of the DFN generated by the 2D algorithm (line 2). Neighbor
cells can still be used in the same way as in 2D to speed up the rejection of
candidates. Unlike in 2D, a maximal Poisson-disk sampling in 3D does not
guarantee sliver-free triangulation, which is why we do not use the cell lists
to find undersampled cells in 3D. Instead, once the algorithm terminates, the
resulting sampling is triangulated (line 10), slivers identified (line 11),
and two nodes of every sliver, with a preference for nodes that are neither on
a boundary or a fracture, are removed (line 12). While the definition of a
sliver given earlier in Section 2.2 allows for a bit of leeway in what is
considered a small or large dihedral angle, we successfully replaced
tetrahedra with dihedral angles outside of $[8^{\circ},170^{\circ}]$ and
aspect ratios bigger than $0.2$. It is worth noting that more traditional ways
of sliver-removal like perturbation [65] or exudation [11] can break the empty
disk property and are not used in nMAPS. Then the algorithm is restarted with
the remaining nodes as the seed set (line 15). This process is repeated until
all slivers are removed (line 16). With this approach, we have obtained
triangulations with no elements of dihedral angles of less than $8^{\circ}$,
examples provided in the next section. The method for generating the
conforming mesh is similar to that for the 2D case, but spheres around
triangle cells of the fracture planes are excavated.

###  3.4 Pseudocode for the 2D and 3D sampling algorithms

Notation for Pseudocodes:
---
Input:
•$D^{3}\subset\mathbb{R}^{3}$: |  Cubical Domain ($\ast$)
•$DFN\subset D^{3}$: |  Generated by dfnWorks ($\ast$)
•$F_{l}\subset\mathbb{R}^{3}$ : |  $l$-th fracture of the DFN
•$q_{l,m}^{(1)}$ and $q_{l,m}^{(1)}$: |  Endpoints of intersection between fractures $F_{l}$ and $F_{m}$
User defined parameters:
•$H/2$: |  Minimum distance between points
•$F$ : |  $HF$ distance round intersections with constant density
•$R$: |  $ARH+H/2$ is maximum distance between points
•$A$: |  Maximum slope of inhibition radius
•$k$: |  Number of concurrently sampled candidates
Additional notation:
•$G$: |  Square cells covering $F_{l}$ with $diam(g)\leq H/2$ for all $g\in G$.
•$\rho({\mathbf{x}})$: |  Piecewise linear function defined in (5)-2D or (6)-3D
•$r({\mathbf{x}},{\mathbf{y}})$: |  Inhibition radius $\min(\rho({\mathbf{x}}),\rho({\mathbf{y}}))$
•$R({\mathbf{x}},{\mathbf{y}})$: |  Coverage radius
•$C({\mathbf{x}})\in G$: |  Grid cell containing the point ${\mathbf{x}}$.
•$N^{+}({\mathbf{x}})$ : |  $\\{g\in G:dist(C({\mathbf{x}}),g)\leq\rho({\mathbf{x}})\\}$: Cells that can contain points ${\mathbf{y}}$ with $|{\mathbf{x}}-{\mathbf{y}}|\leq r({\mathbf{x}},{\mathbf{y}})$
•$N^{-}({\mathbf{x}})$: |  $\\{g\in G:diam(g\cup C({\mathbf{x}}))\leq\frac{\rho({\mathbf{x}})}{1+A}\\}$: Cells, where for all their points $y$ $|{\mathbf{x}}-{\mathbf{y}}|\leq r({\mathbf{x}},{\mathbf{y}})$
•$G_{occ}$ |  $\bigcup\limits_{{\mathbf{x}}\in X}N^{-}({\mathbf{x}})$: Cells on which $X$ is already maximal.
•$\mathcal{T}(X)$: |  Delaunay triangulation of $X$ ($\ast$)
Output: |
•$X$ : |  Poisson-disk sampling on the $l$-th fracture
$(\ast)$: Only used in 3D sampling |

Algorithm 1

2D Poisson-disk sampling

1: Initializing:

2: $X\subset F_{l}$ $\triangleright$ Generate a 1D Poisson-disk sampling with
$R({\mathbf{x}},y)\leq\frac{r({\mathbf{x}},{\mathbf{y}})}{\sqrt{2}(1+L)}$along
boundary $\delta F_{l}$ as seed.

3: for ${\mathbf{x}}\in X$ do

4: $G_{occ}\leftarrow G_{occ}\cup N^{-}({\mathbf{x}})$ $\triangleright$
Initialize occupied cells

5: end for

6: Sampling:

7: $i\leftarrow 1$ $\triangleright$ Start sampling at first accepted point

8: $n\leftarrow|X|$ $\triangleright$ Will increase as more points are accepted

9: while $i\leq n$ do

10: repeat

11: for $j\in\\{1,...,k\\}$ do

12: ${\mathbf{p}}_{j}\in F_{l}$ $\triangleright$ Generate $k$ new candidate
points on the annulus around ${\mathbf{x}}_{i}$

13: if $C({\mathbf{p}}_{j})\in G_{occ}$ then

14: reject ${\mathbf{p}}_{j}$ $\triangleright$ Cell already blocked by
existing point’s inhibition radius

15: else

16: for ${\mathbf{y}}\in N^{+}({\mathbf{p}}_{j})$ do

17: if $|{\mathbf{p}}_{j}-{\mathbf{y}}|<r({\mathbf{p}}_{j},{\mathbf{y}})$ then

18: reject ${\mathbf{p}}_{j}$ $\triangleright$ Empty disk property violated

19: end if

20: end for

21: end if

22: if $p_{j}$ was not rejected then

23: $X\leftarrow X\cup\\{{\mathbf{p}}_{j}\\}$$\triangleright$ Accept
${\mathbf{p}}_{j}$ and add it to the sampling set

24: $G_{occ}\leftarrow G_{occ}\cup N^{-}({\mathbf{p}}_{j})$ $\triangleright$
Update occupied cells

25: $n\leftarrow n+1$ $\triangleright$ Ensures sampling around newly accepted
points

26: end if

27: end for

28: until All $k$ of the ${\mathbf{p}}_{j}$ are rejected

29: $i\leftarrow i+1$ $\triangleright$ Start sampling around next accepted
point

30: end while$\triangleright$ Terminate here or start resampling

Algorithm

Continuation of Algorithm 1 (2D Resampling)

31: Resampling: (Optional)

32: for $C\in G\setminus G_{occ}$ do

33: ${\mathbf{p}}\in C$ $\triangleright$ Generate a random candidate on each
cell

34: for ${\mathbf{y}}\in N^{+}({\mathbf{p}})$ do

35: if $|{\mathbf{p}}-{\mathbf{y}}|<r({\mathbf{p}},{\mathbf{y}})$ then

36: reject ${\mathbf{p}}$ $\triangleright$ Empty disk property violated

37: end if

38: end for

39: if ${\mathbf{p}}$ was not rejected then

40: $X\leftarrow X\cup\\{{\mathbf{p}}\\}$$\triangleright$ Accept
${\mathbf{p}}$ and add it to the sampling set

41: $G_{occ}\leftarrow G_{occ}\cup N^{-}({\mathbf{p}})$ $\triangleright$
Update occupied cells

42: $n\leftarrow n+1$

43: end if

44: end for

45:

46: Rerun algorithm from line 9 (Note: $i$ is not reset)

Algorithm 2

3D Poisson-disk sampling + resampling

1: Initializing:

2: $X\subset D^{3}$ $\triangleright$ Use Algorithm 1 to generate a Poisson-
disk sampling on $\delta D^{3}$ and the DFN

3: for ${\mathbf{x}}\in X$ do

4: $G_{occ}\leftarrow G_{occ}\cup N^{-}({\mathbf{x}})$ $\triangleright$
Initialize occupied cells

5: end for

6: Sampling:

7: The sampling process in 3D is the same as the 2D method (Algorithm 1)
except for the following:

  * A.

New candidates are generated on a spherical shell instead of an annulus

  * B.

A candidate $p\notin\delta D^{3}$ is rejected if $dist(p,\delta
D^{3})<\rho(p)/2$

8: Resampling: (Optional)

9: repeat

10: for $T\in\mathcal{T}(X)$ do

11: if $T$ is a sliver then

12: $X\leftarrow X\setminus\\{{\mathbf{x}},{\mathbf{y}}\\}$, where
${\mathbf{x}},{\mathbf{y}}\in T$ are two random points not contained in the
boundary or the DFN $\triangleright$ Minimum distance of points to DFN and
boundary (7B above) ensures that this is possible.

13: end if

14: end for

15: Rerun algorithm from line 6

16: until $\mathcal{T}(X)$ contains no more slivers.

##  4 Results

In this section, we provide examples of nMAPS and compare its performance with
previous methods.

###  4.1 Two-dimensional Examples

We begin with a network composed of four disc-shaped fractures. Figure 4 shows
the triangulation produced using a variable radius sampling on a single
fracture that contains three intersections. Triangles are colored by their
maximum edge length to demonstrate how their size increases with distance from
the intersections. Figure 5, shows the triangulation of a constant-radius
sampling (uniform resolution mesh) on the fracture but re-assembled into the
network.

Figure 4: Triangulation of variable radii Poisson-disk sampling on fracture
with three intersections. Parameters $H=0.01$, $R=40$, $A=0.1$, and $F=1$.
Triangles are colored according to their maximum edge length. The lines of
intersection are shown as spheres.  Figure 5: Triangulation of a uniform
resolution Poisson-disk sampling based mesh reassembled into the whole DFN.

Figure 6 shows a meshed network that contains 25 fractures whose radii are
generated from an exponential distribution with a decay exponent of 0.3. The
largest number of intersections on a fracture is eight within the network;
this is not a constraint of generation nor the sampling technique. The
inhibition radius parameters are $H=0.1$, $A=0.1$, $F=1$, and $R=40$. The mesh
contains 23,195 nodes and 47,367 triangles. The quality of the triangulation
is presented in the histograms shown in Figure 7: distribution of (a) minimum
angle, (b) maximum angle, and (c) aspect ratio. With the exception of two
elements, all of the triangles have a minimum angle greater than $27^{\circ}$.
The theoretical minimum angle in a maximal Poisson-disk sampling with
Lipschitz constant $A=0.1$ is $27.04^{\circ}$. The two exceptions are
$25^{\circ}$ and $26^{\circ}$. In terms of the maximum angle, there are very
few triangles with angles larger than $110^{\circ}$ and none larger than
$120^{\circ}$. The largest maximum angle theoretically possible in a maximal
Poisson-disk sampling with this Lipschitz-constant is $125.92^{\circ}$. The
vast majority of aspect ratios are greater than $0.8$ with only a marginal
number of triangles having a value less than $0.6$ and none less than $0.47$.

Figure 6: Triangulation of a variable radii Poisson-disk sampling reassembled
into original DFN. The network contains 25 fractures whose radii are generated
from an exponential distribution with a decay exponent of 0.3. Parameters used
in the sampled: $H=0.1$,$R=40$,$A=0.1$,$F=1$. The mesh contains 23195 nodes
and 47367 triangles. The minimum angle is $\geq 25^{\circ}$, the maximum angle
$\leq 120^{\circ}$, and all aspect ratios are $\geq 0.47$.
$0$$200$$400$$600$$800$$1000$$1200$$1400$$25$$30$$35$$40$$45$$50$$55$$60$

no. of triangles

min. angle
[∘]$0$$200$$400$$600$$800$$1000$$1200$$1400$$60$$70$$80$$90$$100$$110$$120$

no. of triangles

max. angle
[∘]$0$$200$$400$$600$$800$$1000$$1200$$1400$$1600$$1800$$2000$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$

no. of triangles

aspect ratio(a)(b)(c) Figure 7: Histograms of selected quality measures of the
triangulation of variable radii Poisson-disk sampling on a fracture with three
intersections. Parameters: $H=0.01$, $R=40$, $A=0.1$, and $F=1$. (a): Minimum
angle ($\geq 25^{\circ}$), (b): Maximum angle ($\leq 120^{\circ}$), (c):
Aspect ratio ($\geq 0.47$)

###  4.2 Run Time Analysis

Next, we present an analysis of the run time and quality of the sampling on a
DFN for varying sample sizes, variations of the parameter $k$, and different
numbers of resampling attempts. All data is generated using the DFN shown in
Fig. 6. All samplings are performed on Fujitsu Laptop with 4 2.5 GHz intel
Core i5 processors and 16GB of RAM Different numbers of nodes were achieved by
changing the parameter ${H}$, the minimum distance between points. All data is
from independent samplings. The plot in Figure 8 shows the run time prior to
the resampling process against the number of points sampled. The color
corresponds to the value of the parameter $k$, which controls the number of
concurrent samples. We see an increase in run time with increasing values of
$k$. The run times for samples with the same $k$ are positioned along straight
lines of slope one, indicating a linear dependence of the total run time and
the number of points sampled. The red lines in the plot have a slope of one
for reference.

Figure 9 shows the relation between the parameter $k$ and the run time. Colors
correspond to different numbers of points. As already established, the run
time increases linearly with the number of points sampled. The run time in
terms of $k$ exhibits a slightly sublinear behavior. The linear fit (black) of
the data on the log-log plot has a slope of $0.7(9)\pm 0.00(7)$. While this
fitting error of $\approx 9\%$ is not insignificant, comparing the data to the
two lines of slope $1$ (red) in the plot indicates that the run time does not
increase more than linearly with $k$.

$0.1$$1$$10$$10^{2}$$10^{3}$$10^{2.5}$$10^{3}$$10^{3.5}$$10^{4}$$10^{4.5}$$10^{5}$

run time [s]

no. of nodes$5$$10$$20$$40$$80$$160$$320$

k

Figure 8: Run time of nMAPS plotted as a function of points sampled prior to
the resampling process. Data points are generated using the same DFN.
Different point densities are generated by changing the minimum inhibition
radius $\frac{H}{2}$ between every pair of points. Data points are colored by
the value of $k$. Other parameters are set to $A=0.1,R=40,F=1$. Comparison to
lines of slope $1$ (red) indicates that the run time increases at an
approximately linear rate.  $10^{-4}$$10^{-3}$$10^{-2}$$10$$10^{2}$$10^{3}$

run time per node [s]

k$0$$10000$$20000$$30000$$40000$$50000$$60000$$70000$

no. of nodes

Figure 9: Run time per node of nMAPS plotted as a function of concurrently
sampled points $k$ prior to the resampling process. Data points generated
using the same DFN. Different point densities generated by changing the
minimum inhibition radius $\frac{H}{2}$ between every pair of points. Data
points are colored by on the total number of points sampled. Other parameters
are set to $A=0.1,R=40,F=1$. Linear fit (black) with slope $0.7(9)\pm
0.00(7)$. Comparison to lines of slope $1$(red) indicate sublinear behavior.

Figure 10 shows a comparison of runtime for nMAPS variable-radii sampling with
and without direct rejection of candidates using the background grid to find
nearby points. Recall that the original method presented by [9] and [14] did
not use direct rejection implemented in this manner. Thus, nMAPS without
direct rejection of candidates is equivalent to the implementation of those
algorithms. Data points generated by nMAPS with direct rejection are
represented by a filled circle, and data points generated without direct
rejection are empty squares. All data points are colored by their $k$ value.
nMAPS with direct rejection has a shorter run time for every pair of data
points. The difference between the methods increases with larger values of $k$
because more candidates are rejected with larger values of $k$, but the direct
rejection method does so without a distance calculation. For $k=5$, the speed
difference between the algorithms is slightly less than a factor of 2,
whereas, for $k=160$, the advantage has grown to about an order of magnitude.
These plots highlight the performance gained by using the additional features
found in nMAPS compared to previous methods.

$1$$10$$100$$10^{2.5}$$10^{3}$$10^{3.5}$$10^{4}$$10^{4.5}$

run time [s]

no. of nodes$5$$10$$20$$40$$80$$160$

k

Figure 10: Comparison of run time for an implementation of [9, 14] (squares)
and nMAPS algorithm (circles) plotted as a function of number of nodes and $k$
value. Data points generated using the same DFN and different point densities
using $H$. Data points are colored depending on the value of $k$. Other
parameters are set to $A=0.1,R=40,F=1$.

A comparison of nMAPS with the method used in dfnWorks [33] to create a
variable resolution point distribution is shown in Table 1. The method
implemented in dfnWorks uses an iterative Rivara refinement algorithm to
generate the nodes of the mesh as described in Section 2.1. We consider three
different DFN to characterize the difference between the methods. The first is
the deterministic network of four ellipses shown in Fig. 5. The second is the
network with fractures sampled from an exponential distribution containing 25
fractures shown in Fig. 6. The final network is composed of a single family of
disc-shaped fractures whose fracture lengths are sampled from a truncated
power-law with exponent 1.8, minimum length 1 m, maximum length 25 m within a
cubic domain with sides of length 100 m. There are 8,417 fractures in this
network. In the first two examples, the algorithms were run on a MacBook Pro
laptop with 8 2.9 GHz Intel Core i9 processors and 32GB of RAM. The third
example was run on a Linux server with 64 AMD Opteron(TM) Processor 6272
(1469.697 MHz) and 252GB of RAM. Note that nMAPS is intrinsically
parallelizable, and is implemented to mesh each fracture independently on a
separate processor. The original meshing technique in dfnWorks is parallelized
in the same manner. The mesh resolution and setup were consistent between the
methods so that the number of nodes in the final mesh is roughly the same. In
all cases, nMAPS was faster than the iterative method, and the speedup
improved with the number of fractures. However, differences in network
properties also likely play a role in the speedup, which is a feature that we
do not explore in this study.

DFN Description | $\begin{matrix}\text{Deterministic}\\\ \text{Ellipses}\end{matrix}$ | $\begin{matrix}\text{Exponential}\\\ \text{Distribution}\end{matrix}$ | $\begin{matrix}\text{Truncated}\\\ \text{Power-law}\end{matrix}$
---|---|---|---
Number of Fractures | 4 | 25 | 8417
Mesh Resolution | Uniform | Variable | Variable
Number of Nodes | $\approx 1\cdot 10^{4}$ | $\approx 1\cdot 10^{5}$ | $\approx 3\cdot 10^{6}$
Number of Processors | 4 | 8 | 32
Iterative - Run Time | 10.64 s | 57.23 s | 47.47 m
nMAPS - Run Time | 4.46 s | 9.91 s | 6.47 m
Speed-up | 2.4 | 5.8 | 7.8
Table 1: Comparison of run time between previous iterative method with nMAPS
for mesh generation. Number of nodes between methods differs slightly. In all
cases, nMAPS outperforms the iterative method.

###  4.3 Quality and resampling

$18000$$20000$$22000$$24000$$26000$$28000$$5$$10$$20$$40$$80$$160$$320$$640$

no. of nodes

k$0$$1$$2$$3$$4$$5$

repetitions

Figure 11: The total number of points accepted after resampling is plotted as
a function of $k$ and colored by the number of resamplings. Other parameters
are set to $A=0.1,R=40,F=1$.

The maximality and density of nMAPS depends on the choice of $k$ and the
number of times resampling is performed. Figure 11 shows the total number of
points accepted after a different number of resampling sweeps for various
values of $k$. Foremost, notice that the density of points increases with $k$.
This rate of growth is largest for small values of $k$ and is lower at higher
values. On the lower end of the $k$ scale, resampling increases the point
density significantly, whereas there is barely a difference when $k>100$. The
first resampling is particularly effective in adding additional nodes, whereas
the difference between each additional resampling decreases thereafter. Given
that resampling adds negligible run time to the original sampling process due
to the efficient background mesh look-up in nMAPS, there is a trade-off
between higher $k$ and more resampling sweeps. Adopting the latter method,
i.e., a larger number of sweeps, can yield better performance in terms of
obtaining a higher density with shorter run times. For example, a run at $k=5$
with a few resampling sweeps results in a density comparable to a run with
more than $10$ times higher $k$ performed without resampling, but the former
is significantly faster. Recall that run time increased linear with $k$.
Similar conclusions can be reached when looking at the quality of resulting
triangulations rather than just the density of the Poisson-disk sampling.

Figure 12 shows the smallest minimum angle in a triangulation using nMAPS for
variable $k$ with different numbers of resampling sweeps. We can see for
$k\gtrapprox 80$ this angle appears to be at around $25^{\circ}$ regardless of
the number of repetitions. The theoretical bound for a maximal Poisson-disk
sampling (with $r({\mathbf{x}},{\mathbf{y}})=R({\mathbf{x}},{\mathbf{y}})$)
for the settings used to generate these data points is $27.04^{\circ}$.
Solving the the angle bounds from Lemma 2 for $\varepsilon$ shows us that in
this sampling
$R({\mathbf{x}},{\mathbf{y}})\lesssim(1+0.1)r({\mathbf{x}},{\mathbf{y}})$.
Given that nMAPS is a stochastic method and identical inhibition/coverage
radii are not exactly guaranteed, these results can be considered very good.
While the quality of triangulations for smaller $k$ values without resampling
is significantly lower, it is worth emphasizing that a single repetition
resolves this issue and produces triangulations with qualities on par with
those produced using significantly higher $k$ values. Thus, nMAPS can be run
using a single or low double-digit value of $k$, perform a single resampling
sweep, and will produce a triangulation just as good as a higher value of $k$
would have produced but in a fraction of the time.

$18$$19$$20$$21$$22$$23$$24$$25$$26$$27$$5$$10$$20$$40$$80$$160$$320$$640$

min. angle [∘]

k$0$$1$$2$$3$$4$$5$

repetitions

Figure 12: Smallest minimum angle of the triangulation produced using nMAPS.
Points are colored by the number of resamplings. Increasing $k$ or the number
of sweeps increases the minimum angle in the mesh. The latter of these
requires much less time.

###  4.4 Three-Dimensional Example

While nMAPS is primarily designed to optimize 2D sampling, we conclude with an
example where these 2D samplings are combined with a 3D sampling of the
surrounding volume to showcase how nMAPS can be used to produce high-quality
3D triangulations as well. We use a simple network of seven square fractures
for clarity in the visualization. The mesh produced by the 3D algorithm is
shown in Fig. 14. The tetrahedra are colored according to their maximum edge
length to highlight how the point density depends on the distance from the
DFN. Note that both the mesh of the DFN and the volume are variable resolution
depending on the distance from the fracture intersections.

$0$$500$$1000$$1500$$2000$$2500$$3000$$3500$$0$$10$$20$$30$$40$$50$$60$$70$

no. of tetrahedra

min. angle
[∘]$0$$500$$1000$$1500$$2000$$2500$$3000$$3500$$70$$80$$90$$100$$110$$120$$130$$140$$150$$160$$170$

no. of tetrahedra

max. angle
[∘]$0$$500$$1000$$1500$$2000$$2500$$3000$$3500$$0$$0.2$$0.4$$0.6$$0.8$$1$

no. of tetrahedra

aspect ratio(a)(b)(c) Figure 13: Histograms of quality measures of the 3D mesh
produced using nMAPS around a seven fracture DFN. Parameters: $H=0.01$,
$R=40$, $A=0.1$, and $F=1$. (a): Minimum angle (all values $\geq 8^{\circ}$),
(b): Maximum angle (all values $\leq 165^{\circ}$), (c): Aspect Ratio (all
values $\geq 0.2$)

Figure 14: (Left) 3D Mesh produced by nMAPS of and around a seven fracture
DFN. Parameters: $H=0.25$, $R=100$, $A=0.125$, and $F=1$. (Right) Close up of
the conforming mesh. Tetrahedra are colored according to their maximum edge
length.

Histograms presented in Fig. 13 show the distribution of mesh quality measures
of the tetrahedra in the 3D triangulation. Tetrahedra with either a dihedral
angle of less than $8^{\circ}$ or an aspect ratio of less than $0.2$ are
discarded before the sampling algorithm was restarted. Histogram (a) shows the
distribution of the minimum dihedral angle of each tetrahedron. As expected,
no dihedral angle below $8^{\circ}$ remains, while the vast majority exceeds
values of $30^{\circ}$. Histogram (b) shows that despite not optimizing with
respect to the maximum dihedral angle, none of these angles exceed
$165^{\circ}$. Histogram (c) shows a sharp cut-off at $0.2$ in the
distribution of aspect ratios, indicating that the aspect ratio is likely to
have been the driving factor for a majority of the resamplings. This example
was run through the sliver-removal resampling process $17$ times to obtain its
triangulation quality. In each of these steps, only $200$ or fewer of $\approx
50000$ nodes were removed before the resampling. This low value of removed
points indicates both the scarcity of slivers in samples generated through
Algorithm 2 and that the vertices of these slivers can successfully be removed
and replaced in a way that does not give rise to new slivers.

##  5 Conclusions

We have presented the near-Maximal Algorithm for Poisson-disk Sampling (nMAPS)
to produce a point distribution designed to generate high-quality variable
resolution Delaunay triangulations in two and three dimensions. We provided
the theoretical basis on which nMAPS is built as well as details of its
implementation. nMAPS uses efficient rejection and resampling techniques to
achieve near maximality and linear run time scaling in the number of points
accepted. We demonstrated that nMAPS could successfully generate variable-
radii Poisson-disk samples on polygonal regions, networks of polygons, and the
surrounding volume they are embedded in. Meshes generated using the point
distributions produced by nMAPS show a quality nearly matching theoretical
bounds for maximal Poisson-disk samplings, in which coverage and inhibition
radii coincide. It is worth noting that near maximality is reached for a
coverage radius just slightly larger than the inhibition radius. We determined
that mesh quality produced by nMAPS is comparable to the method previously
used in dfnWorks but runs in a significantly shorter time. Thus, nMAPS is
significantly faster than the previous conforming variable mesh strategies due
to our efficient rejection techniques that omit costly distance calculations.
It achieves mesh quality only marginally worse than what is theoretically
possible. Moreover, it provides an iterative method where slivers in 3D volume
meshes can be removed entirely from the domain within certain bounds.

It is worth mentioning that nMAPS is fast and simple to run in a parallel
fashion in the context of mesh generation for DFNs, further improving the
overall run time performance. nMAPS is intrinsically parallelizable by working
on each fracture independently on a different processor, as was done in the
performance comparison in Table 1. Based on the grid structure used to accept
and reject candidates, 2D and 3D can also be further parallelized by dividing
their domain into several pieces, which could be sampled individually on
different processors while needing to communicate only cell information along
the boundaries of the split domains. However, once these point distributions
are produced, they all must reside on a single processor to connect them into
a Delaunay mesh.

As a final comment, it’s worth noting that nMAPS is not restricted to variable
resolution DFN mesh generation for conforming numerical schemes. As shown in
Fig. 5, it can be easily used to generate a uniform resolution mesh if
desired. Moreover, nMAPS can be readily used to create meshes of fractures for
non-conforming methods. One merely needs to skip the step in the algorithm
where any point within the circumscribed circle of each segment of the
discretized lines of intersection is removed from the point distribution.
Omitting this step and performing a triangulation, not necessarily a Delaunay
triangulation, will produce a mesh suitable for most non-conforming numerical
schemes. Finally, the general framework of nMAPS can be applied to efficiently
produce arbitrary triangulations in two and three dimensions so long as the
sampling can be constrained to within the domain.

##  6 Acknowledgments

J.K. gratefully acknowledges support from the 2020 National Science Foundation
Mathematical Sciences Graduate Internship to conduct this research at Los
Alamos National Laboratory. J.D.H. and M.R.S. gratefully acknowledges support
from the LANL LDRD program office Grant Numbers #20180621ECR and #20220019DR,
the Department of Energy Basic Energy Sciences program (LANLE3W1), and the
Spent Fuel and Waste Science and Technology Campaign, Office of Nuclear
Energy, of the U.S. Department of Energy. M.R.S. would also like to thank
support from the Center for Nonlinear Studies. J.M.R. received support from
DOE, Contract No. DE-AC05-00OR22725. Los Alamos National Laboratory is
operated by Triad National Security, LLC, for the National Nuclear Security
Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).
This work was prepared as an account of work sponsored by an agency of the
United States Government. The views and opinions of authors expressed herein
do not necessarily state or reflect those of the United States Government or
any agency thereof, its contractors or subcontractors. LAUR # LA-UR-21-24804.

## References

  * [1] I. Babuska and A. K. Aziz.  On the angle condition in the finite element method.  SIAM Journal on Numerical Analysis, 13(2):214–226, 1976.
  * [2] David A Benson, Tomás Aquino, Diogo Bolster, Nicholas Engdahl, Christopher V Henri, and Daniel Fernandez-Garcia.  A comparison of Eulerian and Lagrangian transport and non-linear reaction algorithms.  Advances in Water Resources, 99:15–37, 2017.
  * [3] Marshall Bern, Paul Chew, David Eppstein, and Jim Ruppert.  Dihedral bounds for mesh generation in high dimensions.  In Proceedings of the Sixth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’95, pages 189–196, USA, 1995. Society for Industrial and Applied Mathematics.
  * [4] Inga Berre, Florian Doster, and Eirik Keilegavlen.  Flow in fractured porous media: A review of conceptual models and discretization approaches.  Transport in Porous Media, Oct 2018.
  * [5] Stefano Berrone, Sandra Pieraccini, and Stefano Scialo.  A PDE-constrained optimization formulation for discrete fracture network flows.  SIAM J. Sci. Comput., 35(2):B487–B510, 2013.
  * [6] Stefano Berrone, Stefano Scialò, and Fabio Vicini.  Parallel meshing, discretization, and computation of flow in massive discrete fracture networks.  SIAM Journal on Scientific Computing, 41(4):C317–C338, 2019.
  * [7] Charles Boivin and Carl Ollivier-gooch.  Guaranteed-quality triangular mesh generation for domains with curved boundaries, 2001.
  * [8] Andrea Borio, Alessio Fumagalli, and Stefano Scialò.  Comparison of the response to geometrical complexity of methods for unstationary simulations in discrete fracture networks with conforming, polygonal, and non-matching grids.  Computational Geosciences, 25(1):143–162, 2021.
  * [9] Robert Bridson.  Fast Poisson disk sampling in arbitrary dimensions.  In SIGGRAPH ’07, 2007.
  * [10] M. C. Cacas, E. Ledoux, G. De Marsily, A. Barbreau, P. Calmels, B. Gaillard, and R. Margritta.  Modeling fracture flow with a stochastic discrete fracture network: Calibration and validation: 2. The transport model.  Water Resour. Res., 26(3):491–500, 1990.
  * [11] Siu-Wing Cheng, Tamal K. Dey, Herbert Edelsbrunner, Michael A. Facello, and Shang-Hua Teng.  Sliver exudation.  In Proceedings of the Fifteenth Annual Symposium on Computational Geometry, SCG ’99, pages 1–13, New York, NY, USA, 1999. Association for Computing Machinery.
  * [12] L. Paul Chew.  Guaranteed-quality Delaunay meshing in 3d (short version).  In Proceedings of the Thirteenth Annual Symposium on Computational Geometry, SCG ’97, pages 391–393, New York, NY, USA, 1997. Association for Computing Machinery.
  * [13] Robert L. Cook.  Stochastic sampling in computer graphics.  ACM Trans. Graph., 5(1):51–72, January 1986.
  * [14] Nicholas Dwork, Corey A. Baron, Ethan M.I. Johnson, Daniel O’Connor, John M. Pauly, and Peder E.Z. Larson.  Fast variable density poisson-disc sample generation with directional variation for compressed sensing in mri.  Magnetic Resonance Imaging, 77:186–193, 2021.
  * [15] Mohamed S. Ebeida, Andrew A. Davidson, Anjul Patney, Patrick M. Knupp, Scott A. Mitchell, and John D. Owens.  Efficient maximal Poisson-disk sampling.  ACM Trans. Graph., 30(4), July 2011.
  * [16] Mohamed S. Ebeida, Scott A. Mitchell, Andrew A. Davidson, Anjul Patney, Patrick M. Knupp, and John D. Owens.  Efficient and good Delaunay meshes from random points.  Computer-Aided Design, 43(11):1506–1515, 2011.  Solid and Physical Modeling 2011.
  * [17] Mohamed S. Ebeida, Scott A. Mitchell, Anjul Patney, Andrew A. Davidson, and John D. Owens.  A simple algorithm for maximal Poisson-disk sampling in high dimensions.  Computer Graphics Forum, 2012.
  * [18] J Erhel, J-R de Dreuzy, and B Poirriez.  Flow simulation in three-dimensional discrete fracture networks.  SIAM J. Sci. Comput., 31(4):2688–2705, 2009.
  * [19] Robert Eymard, Thierry Gallouët, and Raphaèle Herbin.  Finite volume methods.  Handbook of numerical analysis, 7:713–1018, 2000.
  * [20] Alessio Fumagalli, Eirik Keilegavlen, and Stefano Scialò.  Conforming, non-conforming and non-matching discretization couplings in discrete fracture network simulations.  Journal of Computational Physics, 376:694–712, 2019.
  * [21] HH Gerke and M Th Van Genuchten.  A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media.  Water Rescour. Res., 29(2):305–319, 1993.
  * [22] Jianwei Guo, Dong-Ming Yan, Li Chen, Xiaopeng Zhang, Oliver Deussen, and Peter Wonka.  Tetrahedral meshing via maximal Poisson-disk sampling.  Computer Aided Geometric Design, 43:186–199, 2016.  Geometric Modeling and Processing 2016.
  * [23] Jianwei Guo, Dong-Ming Yan, Xiaohong Jia, and Xiaopeng Zhang.  Efficient maximal poisson-disk sampling and remeshing on surfaces.  Computers & Graphics, 46:72–79, 2015.  Shape Modeling International 2014.
  * [24] Jianwei Guo, Dongming Yan, Guanbo Bao, Weiming Dong, Xiaopeng Zhang, and Peter Wonka.  Efficient triangulation of poisson-disk sampled point sets.  The Visual Computer, 30(6-8):773–785, may 2014.  KAUST Repository Item: Exported on 2020-10-01 Acknowledgements: This research was partially funded by National Natural Science Foundation of China (Nos. 61372168, 61172104, 61331018, and 61271431), the KAUST Visual Computing Center, and the National Science Foundation.
  * [25] Teklu Hadgu, Satish Karra, Elena Kalinina, Nataliia Makedonska, Jeffrey D. Hyman, Katherine Klise, Hari S. Viswanathan, and Yifeng Wang.  A comparative study of discrete fracture network and equivalent continuum models for simulating flow and transport in the far field of a hypothetical nuclear waste repository in crystalline host rock.  Journal of Hydrology, 553:59 – 70, 2017.
  * [26] J. D. Hyman, , S. Karra, J. W. Carey, C. W. Gable, H. S. Viswanathan, E. Rougier, and Z. Lei.  Discontinuities in effective permeability due to fracture percolation.  Mech. Mater., 119:25 – 33, 2018.
  * [27] J. D. Hyman.  Flow channeling in fracture networks: Characterizing the effect of density on preferential flow path formation.  Water Resources Research, 2020.
  * [28] J. D. Hyman, M. Dentz, A. Hagberg, and P. Kang.  Emergence of stable laws for first passage times in three-dimensional random fracture networks.  Phys. Rev. Lett., 123(24):248501, 2019.
  * [29] J. D. Hyman, M. Dentz, A. Hagberg, and P. Kang.  Linking structural and transport properties in three-dimensional fracture networks.  J. Geophys. Res. Sol. Ea., 2019.
  * [30] J. D. Hyman, C. W. Gable, S. L. Painter, and N. Makedonska.  Conforming Delaunay triangulation of stochastically generated three dimensional discrete fracture networks: A feature rejection algorithm for meshing strategy.  SIAM J. Sci. Comput., 36(4):A1871–A1894, 2014.
  * [31] J. D. Hyman and J. Jiménez-Martínez.  Dispersion and mixing in three-dimensional discrete fracture networks: Nonlinear interplay between structural and hydraulic heterogeneity.  Water Resources Research, 54(5):3243–3258, 2018.
  * [32] J. D. Hyman, Joaquin Jimenez-Martinez, Carl W Gable, Philip H Stauffer, and Rajesh J Pawar.  Characterizing the impact of fractured caprock heterogeneity on supercritical CO2 injection.  Transp. Porous Media, 131(3):935–955, 2020.
  * [33] Jeffrey D Hyman, Satish Karra, Nataliia Makedonska, Carl W Gable, Scott L Painter, and Hari S Viswanathan.  dfnworks: A discrete fracture network framework for modeling subsurface flow and transport.  Computers & Geosciences, 84:10–19, 2015.
  * [34] Jeffrey D Hyman, Satish Karra, Nataliia Makedonska, Carl W Gable, Scott L Painter, and Hari S Viswanathan.  dfnworks: A discrete fracture network framework for modeling subsurface flow and transport.  Computers & Geosciences, 84:10–19, 2015.
  * [35] Jeffrey D. Hyman, Harihar Rajaram, Shriram Srinivasan, Nataliia Makedonska, Satish Karra, Hari Viswanathan, and Gowri Srinivasan.  Matrix diffusion in fractured media: New insights into power law scaling of breakthrough curves.  Geophys. Res. Lett., 46(23):13785–13795, 2019.
  * [36] P. Kang, J. D. Hyman, W. S. Han, and M. Dentz.  Anomalous transport in three-dimensional discrete fracture networks: Interplay between aperture heterogeneity and particle injection modes.  Water Resour. Res., 2020.
  * [37] S Karra, N Makedonska, HS Viswanathan, SL Painter, and JD Hyman.  Effect of advective flow in fractures and matrix diffusion on natural gas production.  Water Resour. Res., 51(10):8646–8657, 2015.
  * [38] Ahmed Khamayseh and Andrew Kuprat.  Anisotropic smoothing and solution adaption for unstructured grids.  International Journal for Numerical Methods in Engineering, 39(18):3163–3174, 1996.
  * [39] P.C. Lichtner, G.E. Hammond, C. Lu, S. Karra, G. Bisht, B. Andre, R.T. Mills, and J. Kumar.  PFLOTRAN user manual: A massively parallel reactive flow and transport model for describing surface and subsurface processes.  Technical report, (Report No.: LA-UR-15-20403) Los Alamos National Laboratory, 2015.
  * [40] Peter Lichtner and Satish Karra.  Modeling multiscale-multiphase-multicomponent reactive flows in porous media: Application to co2 sequestration and enhanced geothermal energy using PFLOTRAN.  In Al-Khoury, R., Bundschuh, J. (eds.) Computational Models for CO2 Geo-sequestration & Compressed Air Energy Storage (http://www.crcnetbase.com/doi/pdfplus/10), pages 81–136. CRC Press, 2014.
  * [41] Yehong Liu and Guosheng Yin.  The Delaunay triangulation learner and its ensembles.  Computational Statistics & Data Analysis, page 107030, 2020.
  * [42] JCS Long, JS Remer, CR Wilson, and PA Witherspoon.  Porous media equivalents for networks of discontinuous fractures.  Water Resour. Res, 18(3):645–658, 1982.
  * [43] A. E. Lovell, S. Srinivasan, S. Karra, D. O’Malley, N. Makedonska, H. S. . Viswanathan, G. Srinivasan, J. W. Carey, and L. P. Frash.  Extracting hydrocarbon from shale: An investigation of the factors that influence the decline and the tail of the production curve.  Water Resour. Res., 2018.
  * [44] N. Makedonska, J. D. D Hyman, S. Karra, S. L Painter, C. W. W Gable, and H. S Viswanathan.  Evaluating the effect of internal aperture variability on transport in kilometer scale discrete fracture networks.  Adv. Water Resour., 94:486–497, 2016.
  * [45] Gary L. Miller, Dafna Talmor, Shang-Hua Teng, and Noel Walkington.  A Delaunay based numerical method for three dimensions: Generation, formulation, and partition.  In Proceedings of the Twenty-Seventh Annual ACM Symposium on Theory of Computing, STOC ’95, pages 683–692, New York, NY, USA, 1995. Association for Computing Machinery.
  * [46] Gary L. Miller, Dafna Talmor, Shang-Hua Teng, Noel Walkington, and Han Wang.  Control volume meshes using sphere packing: Generation, refinement and coarsening.  In In Fifth International Meshing Roundtable, pages 47–61, 1996\.
  * [47] Scott Mitchell, Alexander Rand, Mohamed Ebeida, and Chandrajit Bajaj.  Variable radii Poisson-disk sampling.  Proceedings of the 24th Canadian Conference on Computational Geometry, CCCG 2012, 01 2012.
  * [48] M. K. Mudunuru, S. Karra, N. Makedonska, and T. Chen.  Sequential geophysical and flow inversion to characterize fracture networks in subsurface systems.  Stat. Anal. Data. Min., 10(5):326–342, 2017.
  * [49] Michael Murphy, David M Mount, and Carl W Gable.  A point-placement strategy for conforming Delaunay tetrahedralization.  International Journal of Computational Geometry & Applications, 11(06):669–682, 2001.
  * [50] H. Mustapha and K. Mustapha.  A new approach to simulating flow in discrete fracture networks with an optimized mesh.  SIAM J. Sci. Comput., 29:1439, 2007.
  * [51] Hussein Mustapha, Roussos Dimitrakopoulos, Thomas Graf, and Abbas Firoozabadi.  An efficient method for discretizing 3d fractured media for subsurface flow and transport simulations.  International Journal for Numerical Methods in Fluids, 67(5):651–670, 2011.
  * [52] Engineering National Academies of Sciences, Medicine, et al.  Characterization, modeling, monitoring, and remediation of fractured rock.  National Academies Press, 2020.
  * [53] National Research Council.  Rock fractures and fluid flow: contemporary understanding and applications.  National Academy Press, 1996.
  * [54] Shlomo P Neuman and Joseph S Depner.  Use of variable-scale pressure test data to estimate the log hydraulic conductivity covariance and dispersivity of fractured granites near oracle, arizona.  J. Hydrol., 102(1-4):475–501, 1988.
  * [55] S.P. Neuman.  Trends, prospects and challenges in quantifying flow and transport through fractured rocks.  Hydrogeol. J., 13(1):124–147, 2005.
  * [56] A. Wille Nordqvist, Y. W. Tsang, C. F. Tsang, Björn Dverstorp, and Johan Andersson.  A variable aperture fracture network model for flow and transport in fractured rocks.  Water Resources Research, 28(6):1703–1713, 1992.
  * [57] G Pichot, J Erhel, and J-R de Dreuzy.  A mixed hybrid mortar method for solving flow in discrete fracture networks.  Appl. Anal., 89(10):1629–1643, 2010.
  * [58] G Pichot, J Erhel, and J-R de Dreuzy.  A generalized mixed hybrid mortar method for solving flow in stochastic discrete fracture networks.  SIAM J. Sci. Comput., 34(1):B86–B105, 2012.
  * [59] Karsten Pruess, Curtis M Oldenburg, and GJ Moridis.  Tough2 user’s guide version 2.  1999\.
  * [60] M Cecilia Rivara.  Algorithms for refining triangular grids suitable for adaptive and multigrid techniques.  International journal for numerical methods in Engineering, 20(4):745–756, 1984.
  * [61] Maria-Cecilia Rivara.  Mesh refinement processes based on the generalized bisection of simplices.  SIAM J. Numer. Anal., 21(3):604–613, 1984.
  * [62] T. Sherman, J. D. Hyman, M. Dentz, and D. Bolster.  Characterizing the influence of fracture density on network scale transport.  J. Geophys. Res. Sol. Ea., 125(1):e2019JB018547, 2020.  e2019JB018547 10.1029/2019JB018547.
  * [63] Gilbert Strang, George J. Fix, and D. S. Griffin.  An Analysis of the Finite-Element Method.  Journal of Applied Mechanics, 41(1):62–62, 03 1974.
  * [64] Dafna Talmor, Guy Blelloch, Alan M. Frieze, Noel J. Walkington, and Shang hua Teng.  Well-spaced points for numerical methods.  Technical report, 1997.
  * [65] Jane Tournois, Camille Wormser, Pierre Alliez, and Mathieu Desbrun.  Interleaving Delaunay refinement and optimization for practical isotropic tetrahedron mesh generation.  In ACM SIGGRAPH 2009 Papers, SIGGRAPH ’09, New York, NY, USA, 2009\. Association for Computing Machinery.
  * [66] YW Tsang, CF Tsang, FV Hale, and B Dverstorp.  Tracer transport in a stochastic continuum model of fractured media.  Water Resour. Res, 32(10):3077–3092, 1996.
  * [67] Dong Ming Yan, Jian Wei Guo, Bin Wang, Xiao Peng Zhang, and Peter Wonka.  A survey of blue-noise sampling and its applications.  Journal of Computer Science and Technology, 30(3):439–452, may 2015\.  Publisher Copyright: © 2015, Springer Science+Business Media New York.
  * [68] Robert W Zimmerman, Gang Chen, Teklu Hadgu, and Gudmundur S Bodvarsson.  A numerical dual-porosity model with semianalytical treatment of fracture/matrix flow.  Water Resour. Res, 29(7):2127–2137, 1993.
  * [69] M. ZLAMAL.  On the finite element method.  Numerische Mathematik, 12:394–409, 1968.
  * [70] G Zyvoloski.  FEHM: A control volume finite element code for simulating subsurface multi-phase multi-fluid heat and mass transfer.  Los Alamos Unclassified Report LA-UR-07-3359, 2007.