Generalised primal-dual grids for unstructured co-volume schemes

12/03/2017 ∙ by Darren Engwirda, et al. ∙ Columbia University 0

The generation of high-quality staggered unstructured grids for computational simulation is considered; leading to the development of a new optimisation-based strategy for the construction of weighted `Regular-Power' tessellations appropriate for co-volume type numerical techniques. This new framework aims to extend the conventional Delaunay-Voronoi primal-dual structure; seeking to assemble generalised orthogonal tessellations with enhanced geometric quality. The construction of these grids is motivated by the desire to improve the performance and accuracy of numerical methods based on unstructured co-volume type schemes, including various staggered grid techniques for the simulation of fluid dynamics and hyperbolic transport. In this study, a hybrid optimisation strategy is proposed; seeking to optimise the geometry, topology and weights associated with general, two-dimensional Regular-Power tessellations using a combination of gradient-ascent and energy-based techniques. The performance of the new method is tested experimentally, with a range of complex, multi-resolution primal-dual grids generated for various coastal and regional ocean modelling applications.



There are no comments yet.


page 11

page 14

page 16

page 18

Code Repositories


JIGSAW is a Delaunay-based unstructured mesh generator for two- and three-dimensional geometries.

view repo


JIGSAW is a Delaunay-based unstructured mesh generator for two- and three-dimensional geometries.

view repo


JIGSAW(GEO) is an unstructured mesh generator for geophysical modelling.

view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Co-volume techniques, in which systems of partial differential equations are discretised using a pair of staggered, orthogonal computational grids, represent an important class of numerical methods; leading to a range of schemes for the simulation of complex physical phenomena, including various problems in fluid dynamics and hyperbolic transport. Such schemes include the classical Marker And Cell (MAC) method

harlow1965numerical for the solution of the Navier-Stokes equations, and various Arakawa-type schemes arakawa1977computational for the solution of the primitive equations in geophysical fluid dynamics. Due to the underlying staggered discretisation strategy, such schemes are typically imbued with a range of desirable conservation properties, preserving, for example, various mass-, energy- and rotation-based physical invariants in the simulation of fluid phenomena. Such behaviour makes co-volume schemes natural choices for the simulation of geophysical flows, leading to the development of high-fidelity algorithms for numerical weather prediction, ocean modelling, and planetary climate dynamics.

In an unstructured setting, the development of co-volume schemes leads to a challenging grid generation problem, requiring the construction of a compatible primal-dual pair that tessellates a given domain into a mesh of discrete elements. The pair is required to be a dual structure; a tessellation of the discrete points into a primal simplicial triangulation and its dual polyhedral complex . Here, duality is expressed as a relationship between the faces of the triangulation and its dual, such that each -simplex is associated with a dual vertex (a -cell), each -simplex is associated with a dual edge (a -cell) spanning between the two dual vertices associated with the adjacent -simplexes, and so on, as per memari2011parametrization . Such structures are also typically required to satisfy a number of additional conditions: preserving orthogonality between adjacent faces in and , conforming to general, user-defined constraints on grid-spacing, and consisting of optimally-shaped grid-cells, designed to minimise the error associated with discrete numerical operators. See Figure 1 for additional detail.

Despite the maturity of approaches for the generation of high-quality unstructured triangulations shewchuk1996triangle ; si2015tetgen ; jamin2015cgalmesh ; schoberl1997netgen ; geuzaine2009gmsh ; engwirda2017jigsaw , the construction of primal-dual structures for co-volume techniques remains a challenging task, with the quality of the dual polyhedral grid typically lagging behind that of the primal triangulation. In this study, a new framework for the generation of high-quality primal-dual grids is developed, based on a class of weighted Voronoi-type tessellations known as power diagrams aurenhammer1987power ; edelsbrunner2012algorithms . While this new approach is applicable to co-volume type schemes in general, a particular motivation is the generation of optimal, multi-resolution grids for unstructured ocean modelling and climate dynamics, focusing on techniques appropriate for the mimetic finite-difference/volume schemes presented by, for example, Ringler et al ringler2010unified ; ringler2013multi , Thurburn et al thurburn2009 , and Korn et al KORN2017525 ; KORN2017156 . In the present study, a two-dimensional implementation is developed, leading to the construction of complex, multi-resolution grids for coastal and regional ocean modelling applications.

1.1 Constraints on co-volume grids

Figure 1: Anatomy of an unstructured co-volume scheme for geophysical fluid dynamics, illustrating: (a) locally-orthogonal primal-dual grid staggering, (b) the mimetic finite-difference/volume formulation of Ringler et al ringler2010unified ; ringler2013multi , and (c) a comparison of well-centred and poorly-staggered configurations. In (a), polygonal cells in the dual grid are formed by joining the orthocentres associated with adjacent primal triangles. In (b), the unstructured TRSK scheme employs conservative cell-centred tracer quantities , edge-centred normal velocity components , and auxiliary vertex-centred vorticity variables . In (c), the well-centred (lower) and poorly-staggered (upper) configurations differ in the relationship between dual vertices and primal triangles. In the well-centred configuration, all dual vertices are interior to their parent triangles. In the poorly-staggered configuration, a dual vertex is exterior to its associated triangle (shaded). Note that adjacent primal/dual edges do not intersect in the poorly-staggered configuration.

Due to the nature of the underlying dual discretisation strategy, co-volume schemes impose a unique set of requirements on the associated meshing problem; requiring high-quality primal-dual tessellations that satisfy a number of geometrical constraints. While spatial-staggering and pair-wise orthogonality are basic prerequisites for such schemes, additional conditions are motivated by a desire to maximise the accuracy, stability and efficiency of the associated co-volume formulation. Here, such considerations are explored in the context of the mimetic finite-difference/volume scheme presented by Ringler et al ringler2010unified ; ringler2013multi and Thurburn et al thurburn2009 for the simulation of geophysical flows, though similar arguments are applicable to the optimisation of co-volume formulations in general.

The ‘TRSK’ formulation of Ringler et al is a generalisation of the structured Arakawa C-grid scheme to unstructured primal-dual pairs. In this scheme, fluid pressure, depth, and density degrees of freedom are located within the dual polygonal control volumes, and a set of orthogonal velocity vectors are positioned along the primal grid edges. Additional vorticity degrees of freedom are placed at the vertices of the dual grid cells. Such an arrangement facilitates the construction of a standard conservative finite-volume type scheme for the transport of cell-centred fluid properties, and a mimetic finite-difference formulation for the evolution of velocity components. Vorticity is computed by considering the discrete circulation of velocity about primal triangles. See Figure 

1 for details. Overall, this scheme is known to posses a variety of desirable conservation properties, conserving mass, potential vorticity and enstrophy, and preserving exact geostrophic balance ringler2010unified . The TRSK scheme is currently employed in the Model for Prediction Across Scales (MPAS) for both atmospheric and oceanic modelling skamarock2012multiscale ; ringler2010unified ; ringler2013multi .

Despite its unstructured nature, such a formulation is not applicable to general unstructured tessellations; requiring that a number of auxiliary geometrical constraints also be satisfied. Specifically, such schemes necessitate the use of grids that are not only staggered and pair-wise orthogonal but also well-centred and centroidal. These additional conditions are constraints on the geometry of the underlying primal triangles and their dual polygons.

A well-centred grid vanderzee2008well ; vanderzee2010well is one in which all dual vertices lie within the interior of their associated primal triangles. Such a condition guarantees that adjacent primal and dual edges intersect; inducing a dual structure that is nicely staggered. In the context of the unstructured TRSK scheme described previously, this configuration guarantees that a consistent computational stencil exists for evaluation of the various discrete numerical operators present in the governing equations, allowing such schemes to preserve physical invariants associated with the flow. Centroidal tessellations are those in which primal and dual vertices are positioned at the centroids of the associated staggered grid cells, with the vertices of primal triangles located at the centres-of-mass of the associated dual polyhedra and visa versa. Ideally, centroidal staggering should also be exhibited by the lower dimensional faces of the pair , with adjacent faces intersecting at their local centroids. Such structures can be used to construct compact spatial discretisation schemes that exhibit optimal numerical performance. The unstructured TRSK scheme described previously is known to achieve second-order accurate convergence on uniform primal-dual grids that are well-centred, centroidal and pair-wise orthogonal ringler2010unified .

1.2 Related work

A number of approaches for the generation of high-quality primal-dual grids have previously been pursued, typically focusing on the optimisation of conventional Delaunay-Voronoi type pairs. Such methods include a range of well-known variational techniques, including methods based on Centroidal Voronoi Tessellation (CVT) du1999centroidal ; du2002grid and Optimal Delaunay Triangulation (ODT) alliez2005variational ; chen2004optimal ; chen2011efficient that seek to reposition primal vertices to minimise an objective function that depends on the geometry and topology of the grid. While sophisticated CVT/ODT schemes have been shown to generate high quality centroidal Voronoi-type tessellations, the construction of structures that are additionally well-centred has proven to be a difficult task; requiring the generation of so-called non-obtuse triangulations in which all angles in the primal mesh are bounded below . In vanderzee2008well ; vanderzee2010well , Vanderzee et al have shown that a non-linear optimisation-based strategy can be used to create such non-obtuse triangulations, though the resulting primal-dual tessellations do not appear to be centroidal or necessarily of optimal quality.

More recently, generalised orthogonal primal-dual pairs have been investigated memari2011parametrization , with Mullen et al developing the so-called Hodge-Optimized Triangulation (HOT) mullen2011hot ; goes2014weighted framework in the context of computer graphics applications. The HOT formulation is a generalised variational approach, in which the minimisation of a coupled primal-dual energy functional is undertaken. Such a formulation is designed to optimise the performance of an associated Discrete Exterior Calculus (DEC) type formulation; aiming to generate generalised primal-dual pairs that are typically centroidal, well-centred and pair-wise orthogonal. In this work, Mullen et al employ a weighted primal-dual structure based on the power diagram and its associated regular triangulation, and construct a coupled, optimisation-based approach designed to minimise the combined HOT functional.

The use of weighted primal-dual pairs has also been pursued by Walton et al walton2017advances ; walton2013a ; walton2013b in which a multi-paradigm optimisation scheme is used to improve the well-centredness of two- and three-dimensional unstructured grids for co-volume type aerospace modelling. In this work, a novel modified ‘cuckoo’ type optimisation strategy is employed to minimise a dual mesh quality metric that measures the offset between dual vertices and their associated primal element centroids. These operations are combined with topological modifications, in which clusters of primal elements associated with sufficiently small dual edges are merged into polyhedra. Walton et al report good performance for a variety of two- and three-dimensional problems, and demonstrate that optimised primal-dual meshes often improve the stability and efficiency of co-volume type solvers for problems in aerodynamics and electromagnetism.

In the context of unstructured ocean-modelling, the generation of optimal co-volume grids have typically focused on the application of CVT-based methods, with Jacobsen et al presenting the massively-parallel MPI-SCVT algorithm jacobsen2013parallel for the generation of graded, Voronoi-type meshes on the sphere. These grids have been utilised for various multi-resolution studies using the Model for Prediction Across Scales (MPAS-O) ringler2008multiresolution ; ringler2013multi . Despite producing high-quality centroidal Voronoi structures, the CVT-based approach has not been found to reliably generate well-centred grids; leading to issues with the preservation of rotational invariants in the simulation of geophysical flows, as discussed previously. An alternative, hybrid approach has recently been described by the author in engwirda2017jigsaw , in which a Frontal-Delaunay refinement procedure is used to provide a high-quality initial Delaunay-Voronoi pair that is subsequently improved through an optimisation phase. Such an approach has been shown to successfully generate a range of complex, centroidal and well-centred Voronoi-type grids on the sphere.

In the present study, a generalisation of the approach presented by the author in engwirda2017jigsaw is pursued, but in which a weighted primal-dual structure is optimised in place of the conventional Delaunay-Voronoi pair. Firstly, a new dual mesh quality metric is introduced, designed to minimise the discretisation errors associated with a staggered mimetic discretisation of the primitive equations on unstructured grids. Specifically, the numerical errors associated with the discrete , and operators are analysed, and are used to inform the construction of a new dual grid quality metric. Mesh generation is cast as a coupled optimisation problem, seeking to select a set of vertex positions, weights and an associated mesh topology that maximises a set of primal-dual quality metrics. While the approach presented here is related to both the HOT framework of Mullen et al mullen2011hot ; goes2014weighted and the optimisation-based strategies of Walton et al walton2017advances ; walton2013a ; walton2013b , a number of important differences exist. Here, a monotone optimisation scheme is developed, ensuring that mesh quality metrics are increased monotonically throughout the optimisation process, guaranteeing that no new low-quality elements are created. Secondly, through the construction of a new scale-invariant dual quality metric, the treatment of non-uniform meshes is realised. A coupled optimisation schedule is employed here, interleaving various edge refinement and collapse operations, topological transformations and updates to grid geometry and weights. The performance of the new scheme is analysed experimentally, demonstrating that significantly enhanced primal-dual complexes can be generated for a range of challenging problems involving complex boundary definitions and non-uniform mesh-spacing constraints.

2 Regular triangulations & power diagrams

Figure 2: A comparison of various unstructured primal-dual pairs, showing: (a) a standard Delaunay-Voronoi tessellation, (b) an optimised Regular-Power structure, and (c) the distribution of vertex weights employed in (b). The position of triangle centroids and poorly-staggered elements is highlighted in (a) and (b), demonstrating that the weighted Regular-Power grid is both more centroidal and better self-centred than the corresponding Delaunay-Voronoi mesh. Note that though the underlying Delaunay triangulation is of high quality and is approximately centroidal, its corresponding Voronoi dual incorporates a range of undesirable features, including short edge-segments and non-optimal staggering.

Before describing the new primal-dual grid generation algorithm in full, several concepts associated with the construction of regular triangulations and their associated dual power diagrams are presented. These structures can be viewed as a weighted generalisation of the conventional Delaunay tessellation and its dual Voronoi complex. Such structures are also sometimes referred to as Laguerre diagrams or Dirichlet cell complexes (see, for example aurenhammer1987power ).

Definition 1 (Weighted points & power distance).

A weighted point set is defined as a pair , where are a set of points embedded in -dimensional Euclidean space and are an associated set of scalar weights. The power distance aurenhammer1987power , herein denoted , between an unweighted point and a weighted point is defined as , where is the standard Euclidean distance operator.

Definition 2 (Regular-Power tessellation).

Given a weighted point set , the power complex aurenhammer1987power is the union of polyhedral cells , where each . The associated primal complex is a simplicial triangulation of the weighted points , consisting of the union of simplexes , where each contains the vertices iff mullen2011hot . The primal complex is known as a regular triangulation of the weighted points .

The power complex is a Voronoi-like subdivision of space, where each cell defines a convex region for which the weighted point is closer, in a weighted sense, than all other points in . Such structures, along with their dual regular triangulations, offer a generalised framework for the construction of meshes for co-volume schemes memari2011parametrization ; mullen2011hot ; goes2014weighted ; providing a family of compatible, pair-wise orthogonal primal-dual structures paramaterised by the set of scalar weights . In this work, the selection of an optimal set of weights is sought; tuning the geometry of the power cells to construct a primal-dual structure with optimal characteristics. See Figure 2 for details. Detailed analysis of the properties of power complexes and regular triangulations has been investigated previously, and the reader is referred to, for example, studies by Aurenhammer aurenhammer1987power and Edelsbrunner edelsbrunner2012algorithms for additional detail.

Given a weighted point set , the associated Regular-Power tessellation is fully specified by determination of the position of vertices associated with the polyhedral dual. These points are known as orthocentres, and can be viewed as a weighted generalisation of the circumcentres associated with the standard Delaunay-Voronoi structure. In the present work, orthocentres are constructed with respect to both the edges and faces of the primal triangulation, as detailed in the following sections.

2.1 Weighted ‘edge’ orthocentres

Given a triangle in the weighted primal mesh , the orthocentre associated with the edge is the point of equal power distance to the endpoints


Constraining the point to the edge vector, and solving for the associated line parameter leads to an expression for , such that


Note that in unweighted configurations with , the expression for reduces to a simple equation for the edge midpoint , with , as expected. Note also that is sensitive to the difference in the weights and rather than their magnitude, with equally weighted configurations again recovering .

2.2 Weighted ‘face’ orthocentres

Given a triangle in the weighted primal mesh , the orthocentre associated with the face is the point of equal power distance to the three corner vertices


Following suitable algebraic manipulation111While expressions for cell orthocentres can be written in various forms, note that the linear system shown here does not involve absolute coordinate/weight values, but is instead written in terms of relative differences. As per Shewchuk shewchuk1999lecture , such a formulation is designed to improve numerical robustness when computations are executed in finite-precision arithmetic., these expressions lead to the following system of linear equations for the point


where denotes the difference . Note that in unweighted configurations with , the point is equivalent to the centre of the circumscribing ball associated with . Noting again that the expression for is a function of the local weight differences , , equally-weighted triangles with reduce to

. Such considerations demonstrate that Regular-Power structures differ from conventional Delaunay-Voronoi pairs only in cases where there is a non-uniform distribution of weights throughout the grid.

3 Discretisation errors in a model co-volume scheme

To motivate discussions of optimality in various primal-dual grid configurations, an analysis of the numerical errors associated with a ‘model’ co-volume scheme is presented, based on a mimetic-type discretisation of a nonlinear momentum equation and a hyperbolic transport law


Here is a two-dimensional velocity field , is a passive scalar distribution advected with the flow and is an unspecified potential function. denotes a derivative with respect to time , and . Following, for example, Ringler et al ringler2010unified , the momentum equation (7) can be expressed in so-called vector-invariant form via suitable manipulation of the following vector identity


Substituting (9) into (7), and defining the absolute vorticity and the kinetic energy , the vector-invariant form of the model system can be written as


where is the perpendicular velocity field. As per, for example, Ringler et al ringler2010unified ; ringler2013multi , such a formulation allows for the construction of mimetic discretisations that explicitly conserve rotational and energetic invariants associated with the system, and is typically preferred for the discretisation of geophysical flows as a result.

Figure 3: Mimetic operators defined on a staggered unstructured grid, showing (a) the reconstruction of edge fluxes for the discrete divergence operator defined on the dual cells, (b) evaluation of a discrete normal gradient operator associated with primal edges, and (c) computation of discrete vorticity on the primal cells. The various edge- and face-centred primal-dual ‘defect’ terms are illustrated, consistent with the discretisation errors associated with each operator.

3.1 A staggered mimetic formulation

Following Ringler et al ringler2010unified , the coupled system (10)–(11) can be discretised via a mimetic strategy, in which a set of normal velocity components are positioned on and aligned with the edges of the primal triangulation, conserved cell-mean quantities are located within the cells of the dual grid and a discrete vorticity distribution is computed on the triangles in the primal complex. As per the discussions presented in Section 1.1, such a formulation can be interpreted as an unstructured variant of the well-known Arakawa-type C-grid scheme, with the nonlinear momentum equation (10) expressed in terms of mimetic finite-difference operators, and the continuity equation (11) discretised following a conservative finite-volume approach. See Figure 1 for additional details.

Here, the numerical errors associated with such a formulation are analysed, and such analysis is used to inform the design of optimised primal-dual grid configurations. The effect of various geometrical offsets, here termed ‘defects’, are considered:

  • — the offset between dual vertices and primal cell centres.

  • — the offset between dual edges and primal edge midpoints.

  • — the offset between primal vertices and dual cell centres.

  • — the offset between primal edges and dual edge midpoints.

‘Perfectly’ regular primal-dual pairs consisting of equilateral triangles and uniform hexagonal dual cells satisfy and . In subsequent sections, the impact of ‘imperfect’ grids where , , and are analysed in detail.

3.2 Errors associated with

A finite-volume approximation for the discrete divergence operator in the continuity equation (11) can be expressed as a summation of cell-edge fluxes. Integrating over a given cell in the dual grid


where the edge index iterates over the edges of the cell . Use of a suitable quadrature rule to evaluate the integrals over edges in (12) leads to a fully discrete method. In practice, (12) is typically discretised using a midpoint rule, with and evaluated at a single quadrature point on each edge . Given a perfectly regular tessellation, with , such a formulation achieves second-order accuracy, with the normal velocity components coincident with the quadrature points positioned at the midpoints of the dual edges. Numerical accuracy is degraded as grid distortion increases, with , leading to misalignment between and the associated quadrature points. See Figure 3a for details. Note also that the discretisation becomes non-interpolatory when the grid is not well-centred, with adjacent primal and dual edges failing to intersect in such cases.

3.3 Errors associated with

A finite-difference approximation for the discrete gradient operator in the momentum equation (10) can be formed by exploiting the orthogonality of the underlying primal-dual tessellation. Expressing the momentum equation in terms of normal velocity components


it is clear that an expression for the normal gradients, rather than full gradient vector, is required to be evaluated on each primal edge . Noting that the normal velocity components are positioned at the midpoints of primal edges, the two-point finite-difference


where is the length of the primal edge , and , are values of centred at the two adjacent primal vertices, achieves second-order accuracy. See Figure 3b for detail. Additional consideration is needed when , are not specified directly at primal vertices, but are instead reconstructed from cell-mean data stored on the dual cells. Such an arrangement is consistent with, for example, evaluation of the pressure gradient force in geophysical flows. In such cases, second-order accuracy is maintained only when the tessellation is centroidal, with primal vertices coincident with dual cell centres, such that .

3.4 Errors associated with , ,

A discrete approximation for the nonlinear terms in the momentum equation (10) requires a sequence of composite numerical operations; reconstructing the full velocity vector field from its normal components and evaluating a discrete vorticity distribution. Various vector field reconstruction techniques exist, including Perot’s method PEROT200058

, least-squares interpolation and mixed finite-element formulations

PEIXOTO2014185 , or the TRSK operators introduced in Thurburn et al thurburn2009 . As per the analysis of Peixoto peixoto2016accuracy , these methods typically require perfectly regular tessellations (i.e. ) to achieve higher-order accuracy. The accuracy of the reconstructed perpendicular velocity field and the kinetic energy distribution are constrained accordingly.

A finite-volume type vorticity distribution can be reconstructed over the cells in the primal mesh by exploiting the alignment of the discrete velocity components with primal edges


Here, is the integral-mean vorticity induced within a given primal cell and is the unit tangent vector associated with a primal edge , where, due to the orthogonality of the grid, . Noting that the tangential velocity components are positioned at the midpoints of the primal edges, the cell-mean vorticity can be evaluated to second-order accuracy on the primal cells by discretising (15) using a midpoint quadrature rule. Evaluation of the rotational contribution to the momentum equation is completed by construction of a suitable interpolation operator, expressing the vorticity at the primal edge midpoints in terms of the adjacent primal cell-mean values. Following Ringler et al ringler2010unified , a simple arithmetic average , with , the cell-mean values adjacent to a given primal edge , achieves second-order accuracy when the primal-dual pair is perfectly centred. Such configurations require . See Figure 3c for detail.

3.5 Implications for primal-dual grid staggering

Clearly, the magnitude of the discretisation error associated with various unstructured mimetic numerical schemes ringler2010unified ; ringler2013multi ; thurburn2009 ; KORN2017525 ; KORN2017156 is a strong function of the regularity of the underlying primal-dual grid, with higher-order accuracy achieved only for uniform tessellations with and . In general, imperfections in the grid staggering induce geometrical offsets between the edges and faces of adjacent primal and dual cells, leading to a degradation in the accuracy of the various discrete divergence, gradient and/or rotational terms defined by the model system (10)–(11), as detailed previously. The minimisation of such grid defects is highly desirable as a result.

By construction, conventional Delaunay-Voronoi pairs satisfy , guaranteeing that Voronoi edges cross primal edges at their midpoints. Additionally, so-called Centroidal Voronoi Tessellations (CVTs) ensure that primal vertices are positioned at the centre of their adjacent dual Voronoi cells, corresponding to . Such structures, though, do not explicitly bound the offsets between primal edges and Voronoi edge midpoints, , or the offsets between dual vertices and primal element centroids, . Despite consisting of very high quality primal triangulations, high-quality centroidal Delaunay-Voronoi grids do not necessarily contain optimal dual structure, and are often not well-centred as a result. In the following sections, a generalised primal-dual formulation based on the Regular-Power pair is pursued, in which the distribution of vertex weights is carefully selected to minimise the various primal-dual defect terms. Specifically, by relaxing the constraint associated with conventional Voronoi structures, higher-quality Power diagrams are sought; designed to more evenly minimise the full set of defects , , and . These optimised primal-dual tessellations are expected to lead to improved performance of various staggered unstructured numerical formulations as a result.

4 Weight selection

Before addressing the generation of optimal primal-dual tessellations in full, the task of static weight selection is explored; seeking to choose a set of scalar weights that optimise the geometry of the dual grid associated with a given set of points . Here, the two-dimensional problem is considered, where and the primal-dual tessellation is a pair of two-dimensional complexes, consisting of collections of triangles and polygons, respectively. In the present work, the weight selection task is cast as an optimisation problem; seeking to find a set of weights that maximise a given dual grid quality metric.

4.1 A dual grid quality metric

Figure 4: Anatomy of the dual cost metric (16), showing the contributions to for a single triangle . Here, the triangle and edge centroids and are drawn in red, and the associated weighted orthocentres and are drawn in blue. Non-uniform vertex weights are applied at the vertices of . The dual quality metric is a function of the local face- and edge-centred defects and ; representing the geometrical offset between adjacent faces in the generalised pair .

Given a triangle in the weighted primal mesh , a local dual grid quality metric is proposed; designed to measure the relative geometrical defect between polygonal cells in the orthogonal dual and an ideal centroidal configuration


Here, , are the weighted orthocentres associated with the face and edge segments of the triangle , the points , are the corresponding face and edge centroids, and is a local characteristic length, taken as the mean of the incident edge lengths . Such a scaling serves to non-dimensionalise the various defect terms in (16), and results in a metric that is scale-invariant. The linear coefficients and weight the face- and edge-centred contributions to . In this study, a simple average is employed, with , indicating that equal ‘importance’ is assigned to both face- and edge-centred defects. See Figure 4 for additional detail.

The metric is a smooth function of the geometry of the triangle and the weights centred at its vertices. A maximum value is attained when the dual cell is optimally staggered with respect to the triangle ; spanning between its centroid and edge midpoints, such that and . Conversely, as the relative face- and/or edge-centred defects increase as the local primal-dual staggering becomes less centroidal. The function aims to provide a combined measure of the quality of the staggering between primal and dual grid cells, with the first term in (16) accounting for the defect between the dual grid vertices and triangle centroids, , and the second term the mean defect between dual grid edges and triangle edge midpoints, . Such considerations are designed to minimise the numerical error of the model co-volume discretisation scheme described previously, where a mutually centroidal staggering between the primal and dual grid cells, vertices and edges is optimal.

Optimisation of is expected to lead to configurations with and . When combined with a CVT-like optimisation to also improve the quality of the primal tessellation (i.e. ), the combined scheme is expected to generate very high-quality staggered unstructured grids, exceeding the capabilities of conventional Delaunay-Voronoi pairs. In particular, direct minimisation of is expected to lead to the construction of higher quality dual structures and the generation of well-centred tessellations.

In the context of geophysical fluid dynamics, a detailed error analysis for the mimetic formulation of Ringler et al ringler2010unified ; ringler2013multi is presented by Peixoto in peixoto2016accuracy , and, in conjunction with the analysis presented here in Section 3, shows that the various defects between primal element centroids and dual vertices, and primal and dual edge midpoints are significant sources of overall numerical error; limiting the order of accuracy of the underlying numerical operators. The dual quality metric is intended as a ‘penalty’ operator for such primal-dual offsets, and is expected to improve the accuracy of such staggered unstructured discretisation schemes.

4.2 Selection of locally-optimal weights

Given a weighted primal triangulation , the task of constructing a high-quality orthogonal dual can be cast as an optimisation problem


In general, (18) is a global, non-convex optimisation problem, with the weight at each vertex contributing to the non-linear metric , computed for all adjacent triangles . Rather than seeking solutions to this global problem directly, a simpler, locally-optimal approach is pursued here, inspired by the constrained, gradient-ascent type methods introduced for primal mesh smoothing in, for example, work by Freitag and Ollivier-Gooch freitag1997tetrahedral and Klinger and Shewchuk klingner2008aggressive .

Specifically, the solution of a sequence of local, decoupled problems is considered, with a local optimisation of each given vertex weight computed incrementally. Noting that each weight has local influence, contributing to the metrics for the triangles adjacent to only, a steepest-ascent type update is employed, with




Here, the index is taken as a loop over the primal triangles incident to . The scalar step length is computed via a line search along the gradient ascent vector and, in this study, is taken as the first value that leads to an improvement in the worst-case incident-quality metric . Note that the ascent vector is taken as the gradient of the worst incident quality metric, rather than some mean measure, and that the scheme is biased toward providing updates in a worst-first manner as a result.

A local line search is performed via a simple bisection strategy, iteratively reducing the step-length such that , where is the local line search iterate. The initial guess is computed by considering a first-order Taylor expansion in local grid quality metrics


Here, is a taken as a mean quality value over the local incident set , with the corresponding

an estimate of the weight perturbation required to improve the worst metric until it is equal to the local mean measure. A limited line search is employed, testing iterations

until a successful step is found. If no such improvement is identified, the weight is left unchanged, with . Note that such a procedure is guaranteed to result in monotone improvements in the grid quality function ; ensuring that the worst-case grid quality metric in the incident set is never decreased. This behaviour is here referred to as a quality-preserving update strategy.

The local scheme described in (19)–(27) can be used to construct pseudo-optimal solutions to the global problem (18) via iteration. In this work, a weakly stochastic, Gauss-Seidel type update is employed, consisting of a series of global sweeps over the full set of vertex weights, with each weight updated incrementally using the local scheme described in (19)–(27). Within each outer iteration, vertices are visited in a randomised order; seeking to reduce the likelihood of the iteration becoming trapped in order-induced local optima. Consistent with a Gauss-Seidel type philosophy, updated weight values are used in subsequent metric evaluations as soon as they are available, such that contributes to the calculation of for all .

While such a procedure is not guaranteed to find a globally optimal solution to (18), it does inherit the quality-preserving nature of the local gradient-ascent scheme; ensuring that the worst-case grid quality metrics are improved monotonically. Noting that the accuracy and performance of a numerical discretisation scheme typically scales with the worst quality elements in a grid, rather than some mean measure, it is suggested that such behaviour is highly desirable; consistent with the development of so-called smart node smoothing schemes for primal meshes, as per freitag1997tetrahedral ; klingner2008aggressive .

5 Optimisation of primal-dual grids

Building on the weight-selection procedure introduced in Section 4, the full primal-dual grid generation algorithm is presented; seeking to find a set of vertex positions and associated vertex weights that induce a primal-dual pair with optimal geometrical and topological characteristics. Consistent with the weight-selection strategy presented previously, the grid generation task is cast as an optimisation problem; seeking to maximise a pair of primal-dual grid quality metrics. All primal and dual optimisation operations are designed to be quality-preserving; ensuring that improvements in the local, worst-case incident grid quality metrics and/or are achieved where possible. Such behaviour is designed to ensure the overall robustness of the scheme; resulting in a coupled optimisation strategy designed to improve the worst-case elements in a given primal-dual grid.

5.1 A primal grid quality metric

Given a triangle in the weighted primal mesh , the area-length ratio Parthasarathy94 is employed as a local primal grid quality metric ; designed to measure the relative distortion and irregularity of the primal grid cells.


Here, is the signed area of the triangle and is the associated root-mean-square edge length. The area-length ratio is a robust, scalar measure of triangle shape quality, and is normalised to achieve a score of for ideal elements. The area-length ratio decreases with increasing cell distortion, achieving a score of for degenerate elements and for ‘tangled’ elements with reversed orientation.

As per Shewchuk shewchuk2002good , the area-length metric is a measure of the condition number of various numerical operators and transformation matrices expressed on the primal grid. Compared to other primal mesh quality criteria, it is numerically well-behaved, is relatively cheap to evaluate and will detect all classes of low quality simplexes. The area-length metric and its generalisations have been successfully employed in a number of studies, including, for example, work by Klinger and Shewchuck klingner2008aggressive , Gosselin and Olliver-Gooch Gosselin2011 , as well as previous investigations by the author engwirda2017jigsaw .

5.2 Selection of locally-optimal vertex positions

Considering, firstly, the geometrical optimality of the primal grid , a mesh-smoothing procedure is undertaken; seeking to reposition vertices to improve both the primal grid quality metrics and conformance to a user-defined mesh-spacing function , defining the distribution of desired primal edge lengths throughout the domain. A hybrid, two-pass strategy is pursued; employing a variation of the Optimal Delaunay Triangulation (ODT) scheme of Chen et al chen2004optimal ; chen2011efficient as a first-pass update, and falling back to a constrained, gradient-ascent type approach when required.

5.2.1 ‘Optimal’ weighted triangulations

Given a vertex in the primal mesh , an initial update is attempted using a variation on the ODT strategy of Chen et al chen2004optimal ; chen2011efficient . In the ODT formulation, primal vertices are repositioned such that an element-wise energy-functional is minimised; leading to the optimal piecewise linear reconstruction of quadratic functionals over an associated Delaunay triangulation . The ODT scheme is known to be an effective Delaunay-based mesh optimisation strategy; leading to high quality Delaunay triangulations and Voronoi diagrams that can be adapted to general, spatially-dependent mesh-density functions chen2004optimal ; chen2011efficient .

In the conventional ODT formulation, primal vertices are repositioned to the weighted mean of adjacent triangle circumcentre coordinates; exploiting the orthogonal structure of the standard Delaunay-Voronoi tessellation. In the current work, such a strategy is modified to incorporate the generalised, weighted primal-dual structure pursued throughout; replacing the triangle circumcentres with the associated element orthocentres. This modification leads to a scheme in which primal vertex positions are updated as a weighted sum of the adjacent dual vertex coordinates; maintaining conceptual consistency with the conventional, unweighted ODT strategy. The resulting scheme is denoted Optimal Weighted Triangulation (OWT)222Though additional methods are not pursued in the present study, note that a generalisation of a range of existing variational mesh optimisation schemes, including the popular Centroidal Voronoi Tessellation (CVT) du1999centroidal ; du2002grid could be achieved through similar modifications to the definition of dual vertices; replacing element circumcentre coordinates with those of the associated weighted orthocentre. throughout. Adapting the ODT-style update strategy of Chen and Holst chen2011efficient , primal vertices are repositioned such that


Here, denotes the star of — the local set of elements adjacent to . denotes the spacing-weighted area of the triangle , and the summation of such terms over the set . The points are the orthocentres associated with the triangles and is a relaxation factor, computed via a local line search. In this study, the weighted-area factors are computed using a 1-point quadrature rule


Conventionally (see, for example chen2011efficient ), such terms are computed with respect to an imposed mesh-density function rather than a mesh-spacing function . Here, the relation is used to express density in terms of spacing, where is the dimensionality of the mesh.

A local relaxation procedure is undertaken in an attempt to determine a quasi-optimal update. A series of local iterates are considered, setting for ; terminating as soon an improvement in the worst-case incident-quality metric is achieved. This constrained update strategy is designed to achieve monotonic behaviour; guaranteed to only enact vertex updates that improve the worst-case incident grid quality metrics. In cases where no improvement is possible, an alternative gradient-ascent type strategy is pursued.

5.2.2 Direct gradient-ascent

While the OWT-type iteration is effective in improving the quality and mesh-spacing conformance of a grid on average, it is not guaranteed to improve the metrics in all cases. Therefore, an additional gradient-ascent type strategy is pursued freitag1997tetrahedral ; klingner2008aggressive ; engwirda2017jigsaw , in which vertex positions are adjusted based on the local gradients of incident element-quality functions , following a procedure similar to that developed for weight-selection in Section 4. Specifically, a given vertex is repositioned along a local steepest-ascent vector, chosen to improve the quality of the worst incident element




Consistent with the steepest-ascent strategy outlined in Section 4, the index is taken as a loop over the primal triangles incident to the vertex , is a scalar step-length, and the ascent vector is taken as the gradient of the worst incident quality metric. Such a strategy is designed to provide updates in a worst-first manner as a result.

A local line search is performed via a simple bisection strategy, iteratively reducing the step-length such that , where is the local line search iterate. The initial guess is computed by considering a first-order Taylor expansion in local grid quality metrics


Here, is a taken as a mean quality value over the local incident set , with the corresponding an estimate of the perturbation required to improve the worst metric until it is equal to this local mean measure. A limited line search is employed, testing iterations until a successful step is found. If no such improvement is identified, the vertex is left unchanged, with . Consistent with previous discussions, such a procedure is guaranteed to monotonically improve the worst-case grid quality metric in the incident set, constituting a quality-preserving strategy.

5.3 Incremental updates to mesh topology

Figure 5: Topological operations for primal grid optimisation, showing (left) an edge-flip, (middle) an edge-collapse, and (right) an edge-refinement operation. Grid configurations before and after each update are shown in the upper and lower panels, respectively.

Following perturbations to the geometry of the primal and/or dual grid, adjustments must be made to the underlying mesh topology, such that the triangulation remains a valid regular tessellation and the dual its orthogonal power diagram. While it is possible to simply re-compute the full weighted triangulation after each adjustment is made, such an approach would impose a significant computational burden, especially when considering that a majority of updates involve small perturbations to the vertex positions or weights; unlikely to require large-scale changes to the grid topology. In this work, an alternative strategy is pursued, based on local element-wise transformations, known as topological flips.

Given a pair of adjacent triangles a local re-triangulation can be achieved by flipping the local connectivity about the shared edge , forming a new edge between the opposing vertices . This operation replaces the existing triangle pair with a new set . See Figure 5a for details. In the present work, the iterative application of such edge-flipping operations is used to incrementally transform the topology of the primal triangulation; ensuring that the local weighted tessellation criterion is satisfied throughout the optimisation process. Specifically, given a general triangulation , possibly violating the local weighting criterion, a cascade of edge-flips are employed to transform into a valid regular triangulation . For each adjacent pair , a flip is enacted if a local violation of the weighted criterion is detected. New triangles created by successful edge flips are subsequently re-examined until no further modifications are required. This approach is a generalisation of the standard incremental edge-flipping algorithms for Delaunay triangulations, previously presented by various authors, including, for example, Lawson lawson77flips .

Given a pair of adjacent triangles , the local weighted criterion is violated when the power distance and between the staggered vertices and the opposing orthocentres is less than the power-radii of the orthogonal balls associated with and . Such a constraint is equivalent to requiring that adjacent cells in the dual power diagram be non-overlapping. To prevent issues with exact floating-point comparisons, a tolerance-based comparison strategy is employed. Specifically, both opposing vertices are required to penetrate their opposing orthoball by a squared distance before a flip is enacted, with and in the current double-precision implementation. See Figure 5a for details.

5.4 Refinement & edge collapse

In addition to updates to the bulk geometry and topology of the grid, mesh quality and mesh spacing conformance can often be further improved through the addition and/or removal of vertices. In the present work, such behaviour is achieved through a set of edge-refinement and edge-collapse operations, generalising the methodology presented by the author in engwirda2017jigsaw to handle general weighted primal-dual tessellations.

Given an edge in the primal tessellation , a collapse operation is achieved by merging the two vertices to a local mean position and re-triangulating the local cavity incident to the edge. In the present work, vertices are merged to an average of adjacent element orthocentre coordinates, such that for all local triangles . While such an approach is slightly more computationally intensive than simpler strategies based on collapsing vertices to primal edge midpoints, the use of averaged orthocentres has been found to be substantially more effective in practice. See Figure 5b for details.

Given an edge in the primal tessellation , a refinement operation is achieved by inserting a new vertex , positioned at the centre of the orthoball associated with the lower quality adjacent triangle . Insertion of the new vertex induces a re-triangulation of the local cavity ; constructed by expanding about in a local greedy fashion. Starting from the initial cavity , where are the triangles adjacent to the edge , additional elements are added to the cavity in a breadth-first manner, with a new, unvisited neighbour added if doing so will improve the worst-case grid-quality metric of the re-triangulated configuration. The final cavity is therefore a locally optimal stencil for the insertion of . In practice, this iterative deepening of typically convergences in one or two iterations, and can be implemented efficiently using local adjacency information. See Figure 5c for details.

Both the edge-collapse and edge-refinement operations are designed to achieve an improvement in worst-case primal grid quality metrics. A sweep over primal mesh edges is performed, with both a collapse and refinement operation is first simulated for each edge. Given an edge , a subsequent collapse or refinement operation is then accepted if the re-triangulated cavity is of improved quality, such that , where is a quality improvement threshold, and the indices and iterate over the local cavities, considering all and .

5.5 A coupled optimisation schedule

1:function OptimisePrimalDual()
2:     // Optimise a given primal-dual pair , employing
3:     // a combination of geometrical, topological and weight-
4:     // based operations.
5:     for  do .5 // coupled outer iterations
6:         for  do .5 // vertex + weight updates
8:              for  do .5 // vertex updates
9:                  Find update direction
10:                  Perform line search
11:                  Monotonic updates
12:              end for
13:              for  do .5 // weight updates
14:                  Find update direction
15:                  Perform line search
16:                  Monotonic updates
17:              end for
18:         end for
19:         // Update grid topology
21:         // Refine/collapse edge
24:     end for
25:     return optimised primal-dual complexes and
26:end function
Algorithm 1 A Coupled Primal-dual Optimisation Strategy

The full primal-dual grid optimisation procedure is realised as a combination of the various geometrical and topological operations described previously; organised into a particular iterative optimisation schedule. See Algorithm 1 for details. Each outer iteration consists of a fixed set of operations: eight sweeps to update vertex positions and weights, an iterative edge-flipping scan to restore the local ‘weighted’ triangulation criterion, and, finally, a single pass of edge refinement/collapse operations. It was found that combining multiple vertex and weight updates with a single stage of topological transformation led to enhanced results in practice, allowing substantial improvement to the geometry of the grid to be achieved before making discrete changes to its topology. In this study, a maximum of sixteen outer iterations were employed. In typical cases, it was found that additional outer iterations did not further improve the worst-case grid quality metrics. Each vertex- and weight-update pass is implemented as a composite operation, with the variational OWT-type technique supplemented with local gradient-ascent iterations as required. Given a vertex in the primal mesh, an OWT update is always attempted first, with a subsequent gradient-ascent step employed in cases where local grid-quality metrics are sufficiently improved y the initial OWT step. Updates to the vertices and weights

occur sequentially, with a single, linear sweep over the grid vertices followed by a single pass over the mesh weights. Each vertex- and weight-update pass is arranged to follow a weakly stochastic, Gauss-Seidel type philosophy; with vertices and weights visited in a randomised order, and updated position and weight values leveraged in subsequent calculations soon as they are available. The coupled optimisation schedule employed here is not based on a particular theoretical derivation, but is simply a set of heuristic choices that have proven to be effective in practice, building on previous primal mesh optimisation techniques

freitag1997tetrahedral ; klingner2008aggressive . A similar mesh improvement strategy, combining a series of geometrical and topological operations, was previously pursued by the author in engwirda2017jigsaw for the optimisation of conventional Delaunay-Voronoi structures.

6 Experimental results

The performance of the primal-dual optimisation algorithms presented in Sections 4 and 5 was investigated experimentally, with the methods employed to generate a set of optimised Regular-Power tessellations for a series of benchmark problems. A pair of test-cases appropriate for regional ocean modelling applications were investigated. All schemes have been implemented in C++. The full algorithm has been implemented as part of the JIGSAW meshing package, and is currently available online JIGSAW-GEO or by request from the author. All tests were completed on a Linux platform, using a single core of an Intel i7 processor. Visualisation and post-processing was completed using MATLAB.

6.1 Preliminaries

Primal-dual grids for the Box test-case, showing: (i) the initial grid, and (ii) a generalised dual grid generated using the weight-selection algorithm. Both the primal triangulation Primal-dual grids for the Box test-case, showing: (i) the initial grid, and (ii) a generalised dual grid generated using the weight-selection algorithm. Both the primal triangulation
Figure 6: Primal-dual grids for the Box test-case, showing: (i) the initial grid, and (ii) a generalised dual grid generated using the weight-selection algorithm. Both the primal triangulation and the dual power diagram are shown here. Dual cells are coloured by their relative power, with blue tones indicating negative values, orange tones positive values, and white tones associated with values approaching zero. Triangles highlighted in red violate the well-centred constraint, and do not contain their own dual vertices. Note that this is a dual-only example; no optimisation of the primal grid is conducted.
Primal-dual grids for the Box test-case, showing: (i) the initial grid, and (ii) a generalised dual grid generated using the weight-selection algorithm. Both the primal triangulation Primal-dual grids for the Box test-case, showing: (i) the initial grid, and (ii) a generalised dual grid generated using the weight-selection algorithm. Both the primal triangulation
Figure 7: Histograms of primal-dual quality metrics for the Box test-case, showing data for: (a) the initial grid, and (b) a generalised dual grid generated using the weight-selection algorithm. Note that this is a dual-only example; no optimisation of the primal grid is conducted..

Three benchmark problems were considered; the first designed to assess the performance of the weight selection strategy described in Section 4 and the second and third designed to test the coupled primal-dual strategy developed in Section 5. In the first example, a generalised dual grid is generated for a given static primal tessellation, using the locally-optimal weight-selection strategy to improve the quality of the dual grid geometry while keeping the underlying primal mesh unchanged. In the second and third examples, the fully coupled optimisation problem is considered; employing the full suite of optimisation techniques to select optimal vertex positions, vertex weights and mesh topology. In all cases, the initial primal mesh consists of a high-quality conforming Delaunay triangulation; generated using the JIGSAW library. All problems are initialised with , implying and . The primal and dual grids associated with each test-case are shown in Figures 6.1, 8 and 9, with associated distributions of grid-quality metrics presented in Figures 6.1 and 10. Grid quality statistics and computational performance is summarised in Table 1.

For all test problems, detailed statistics on primal and dual cell quality are presented, including histograms of the primal and dual quality metrics and defined in (16) and (22). The distribution of angles in the primal mesh as well as the relative-length metric are also presented, where represents conformance to the imposed mesh-spacing function sampled at each primal edge . Histograms highlight the minimum, maximum, and mean values of the relevant distributions as appropriate. Additional measures of spread in the distribution of angles and mesh-spacing conformance are also provided, defined as the Mean Absolute Deviation in and , where


where and denote the means of their respective distributions. Smaller values of and indicate distributions that are more tightly clustered about their mean values, consistent with meshes of higher average quality. A measure of the distribution of weights in each grid is also provided, represented by what is here called the relative power associated with each cell in the dual tessellation


Here, is the scalar weight assigned to the primal vertex associated with a given dual cell and is the area of . Noting that the weights can be interpreted as the squared-radii of circles positioned at the primal vertices aurenhammer1987power , the relative power is a measure of the relative magnitude of the weights applied throughout a mesh. is a non-dimensional metric; allowing for a comparison of weights applied throughout graded tessellations where cells are of varying sizes. Contours of are presented for each test-case considered in the following sections; illustrating the distribution of weights generated by the new primal-dual optimisation scheme.

6.2 Standalone dual grid optimisation

In the ‘Box’ test-case shown in Figure 6.1, the performance of the weight-selection algorithm developed in Section 4 is investigated; seeking to generate an optimised Power diagram given a static primal tessellation. This test problem was designed as a simple verification case and consists of a coarse Delaunay triangulation of a concentric square geometry. The Delaunay-refinement algorithm available in the JIGSAW package engwirda2017jigsaw ; JIGSAW-GEO was used to provide the initial configuration; defined as a quality Delaunay tessellation, such that all angles in the primal elements exceed . The associated unperturbed dual is the conventional Voronoi diagram. A set of primal-dual pairs are shown in Figure 6.1, illustrating the effect of weight optimisation on the geometry of the grid. Dual cells are coloured by their relative power , showing the distribution of weights throughout the grid. Poorly staggered elements are also highlighted, identifying primal triangles that do not contain their associated dual vertices. Overall mesh quality is quantified in Figure 6.1, where histograms of primal and dual grid quality metrics illustrate the distribution of , and before and after the application of the weight-selection scheme.

6.2.1 Discussions

A comparison of the primal-dual pairs shown in Figure 6.1 reveals the effectiveness of the proposed weight-selection procedure, demonstrating that a generalised dual structure can be generated that is both more centroidal and more well-centred than the standard Delaunay-Voronoi configuration. In the initial state, 11 triangles are poorly staggered, with their associated dual vertices lying exterior to the elements themselves. The optimised primal-dual tessellation contains no poorly staggered triangles, with dual vertices shifted toward their associated triangle centroids through the selection of an appropriate set of vertex weights. Improvements to the structure of the dual are confirmed through an analysis of the grid quality metrics presented in Figure 6.1, confirming that both the minimum and mean values of are improved in the generalised weighted configuration. In this case, and , respectively. Note that a slight change in the primal metric is also recorded; due to differences in the topology of the associated Delaunay and Regular triangulations. Though a direct optimisation of the primal grid is not employed in this example, the topology of the triangulation is modified as required; ensuring that the primal-dual pair remains orthogonal. In this case, optimisation of the vertex weights leads to a single edge flip in the Regular triangulation.

6.3 Coupled primal-dual optimisation



& & cpu (sec) & & & & &

australia/coral sea & 54,017 & – & 0.525 & 0.986 & 0.621 & 0.970 & 2073
australia/coral sea (-p) & – & 5.25 & 0.564 & 0.993 & 0.633 & 0.980 & 180
australia/coral sea (-pd) & – & 8.54 & 0.857 & 0.997 & 0.633 & 0.983 & 8

antarctic peninsula & 145,311 & – & 0.450 & 0.983 & 0.594 & 0.965 & 6401
antarctic peninsula (-p) & – & 20.7 & 0.450 & 0.992 & 0.594 & 0.979 & 333
antarctic peninsula (-pd) & – & 27.4 & 0.862 & 0.997 & 0.604 & 0.981 & 11

Table 1: Summary of results for the East Australia and Antarctic Peninsula test-cases, enumerating: the size of the primal triangulation , algorithm runtime CPU, minimum and mean primal/dual grid quality metrics and , and the number of poorly staggered triangles . Data is provided for the initial configuration, primal-only optimisation (-p), and coupled primal-dual optimisation (-pd).
Figure 8: An optimised primal-dual grid for the Australian test-case, showing (clockwise from top-left): (a) the coastal geometry and open-boundary definition, (b) the weighted triangulation (c) contours of the mesh spacing function , and (d) the dual power diagram . Here, dual cells are coloured by their relative power, with blue tones indicating negative values, orange tones positive values, and white tones associated with values approaching zero.

In the ‘Australian’ and ‘Antarctic’ test-cases shown in Figures 8 and 9, the performance of the coupled primal-dual optimisation algorithm developed in Section 5 is investigated; seeking to generate optimal primal-dual pairs appropriate for a co-volume type discretisation scheme. These test problems were designed to explore the effectiveness of the proposed methods for the generation of very high-quality staggered unstructured grids for regional ocean modelling applications, focusing on the East Australian Shelf & Coral Sea region in the ‘Australian’ example, and the Antarctic Peninsula & Weddell Sea region in the ‘Antarctic’ problem. In both cases, a two-dimensional, coastally conforming geometry was generated by employing a locally-centred stereographic projection; a conformal Spherical-to-Cartesian mapping designed to minimise distortion about a region of interest. Additional bathymetric data was retrieved from the well-known ETOPO1 dataset amante2009etopo1 ; providing a high-resolution representation of the local sea-floor. These test-cases will form a basis for future modelling experiments using the Model for Predication Across Scales (MPAS-O) ringler2008multiresolution ; ringler2013multi , in which an unstructured ocean circulation model will be employed to study ocean dynamics in a multi-resolution framework; transitioning from global to coastal scales.

In both problems, the Frontal-Delaunay refinement algorithm provided by the JIGSAW package engwirda2017jigsaw ; JIGSAW-GEO was used to produce an initial configuration, consisting of a high-quality, graded Delaunay triangulation and its associated Voronoi dual. The Frontal-Delaunay algorithm is known to produce particularly high-quality Delaunay triangulations engwirda2016off ; guaranteeing that all angles in the primal triangles exceed , and that the distribution of element sizes is a tight approximation to the target grid-spacing function . Use of the high-quality Delaunay-Voronoi grid generated by the Frontal-Delaunay refinement procedure can be thought of as providing a ‘good’ initial condition for the subsequent primal-dual optimisation phase. In this study, grid-spacing was scaled based on the maximum expected barotropic wave speed


where is the mean ocean depth, is the acceleration due to gravity and is a user-defined parameter used to control effective mesh resolution. The scalars and are lower and upper bounds, here set to and . Such grid-spacing is designed to achieve a quasi-uniform CFL restriction on model time-step throughout the domain and to provide high mesh resolution in shallow regions on the shelf and in the coastal zone. Noting that bathymetry, and hence ocean depth is typically non-smooth, a additional gradient-limiting process is applied to the raw grid spacing function to generate a Lipschitz-continuous grid-spacing function , subject to


Following the work of Persson Persson06SizeFunc , the smooth function is taken as the steady-state solution of the associated Hamilton-Jacobi equation


In this study, a scalar smoothing parameter is used to globally, and isotropically limit spatial gradients. The gradient-limited size function becomes more uniform as . Here, is used throughout.

Figure 9: An optimised primal-dual grid for the Antarctic test-case, showing (clockwise from top-left): (a) the coastal geometry and open-boundary definition, (b) the weighted triangulation (c) contours of the mesh spacing function , and (d) the dual power diagram . Here, dual cells are coloured by their relative power, with blue tones indicating negative values, orange tones positive values, and white tones associated with values approaching zero.

6.3.1 Discussions

A set of optimised primal-dual pairs are presented in Figures 8 and 9, showing: (a) the geometry of each test problem, (b) the primal triangulation , (c) contours of the target grid-spacing function , and (d) detail of the associated dual structure . Results are shown after the application of the primal-dual optimisation procedure. Dual cells are coloured by their relative power , showing the distribution of weights throughout the grid. Mesh quality is quantified in Figure 10, where histograms of primal and dual grid quality metrics illustrate the distribution of , , and before and after the application of mesh optimisation. Results for three different configurations are presented: (i) the initial configuration, (ii) the result of primal-only optimisation, and (iii) the result of the new coupled primal-dual optimisation method presented in Section 5. The primal-only scheme employed in (ii) is based on an optimisation of the primal grid only, leading to the generation of standard Delaunay-Voronoi pairs with . This primal-only scheme contains the full set of operations descried in Section 5: the hybrid ODT-type vertex smoothing operator, topological flips, and edge refinement/collapse operations, but does not include the weight-selection strategy described in Section 4. This primal-only procedure is similar to the strategy presented by the author in engwirda2017jigsaw for the optimisation of Delaunay-Voronoi pairs, and is used here to quantify the added benefit derived from a coupled optimisation of the full primal-dual grid structure compared to that of the primal tessellation alone.

Figure 10: Histograms of primal-dual quality metrics for the Australian (left) and Antarctic (right) test problems, showing metrics for (from top to bottom): (a) the initial grid, (b) application of the primal-only optimisation scheme, and (c) application of the new coupled primal-dual approach.

The primal-dual pairs presented in Figures 8 and 9 show that the coupled optimisation scheme is effective in generating high-quality, graded primal-dual tessellations for problems involving complex boundary definitions, such as coastlines, and non-uniform mesh spacing constraints. Visually, grid quality appears to be very high, with the dual grids consisting of a majority of regular hexagonal polygons, smoothly transitioning between regions of high and low resolution. The effectiveness of the coupled scheme is confirmed through an analysis of the grid quality metrics presented in Figure 10, showing that optimisation improves the minimum and mean values for both the primal and dual metrics and , as well as sharpening the distribution of angles in the primal mesh about and mesh-spacing conformance such that . In the ‘Australian’ test problem, application of the primal-dual scheme improves grid quality such that , , , and . Similar improvement was reported for the ‘Antarctic’ test-case, with , , , and . The full set of results are summarised in Table 1.

Compared to the primal-only scheme, application of the coupled optimisation strategy was found to induce larger improvements in the primal and dual grid quality metrics for both test problems; demonstrating the effectiveness of such an approach. Interestingly, as well as significant improvements to the dual metrics , the coupled algorithm was able to improve the primal triangulations slightly more than the primal-only scheme, resulting in sharper distributions of and slightly improved . Such results present interesting questions for the development of ODT/CVT schemes, showing that, despite such schemes being derived using the standard Delaunay-Voronoi pair, they may benefit in practice from use of the generalised, weighted structures employed here. The coupled optimisation strategy was also found to be significantly more effective than the primal-only scheme in enhancing the well-centredness of the grids, reducing the number of poorly staggered triangles by a factor of 25-30 in both cases. The small number of poorly-staggered elements not removed via optimisation were observed to span highly non-convex segments of the boundary, with all triangle vertices fixed on the bounding contour. In the present implementation, such configurations offer few opportunities for optimisation, admitting perturbations to weight values only. It is expected that these cases may be remediated by employing more sophisticated updates for boundary constrained entities, for example, allowing mesh vertices to ‘slide’ along the bounding contour. In practice, the generation of complex, well-centred Delaunay-Voronoi grids is known to be extremely difficult vanderzee2008well ; vanderzee2010well , showcasing the utility of the new Regular-Power structures presented here. Overall, these results demonstrate that the coupled, primal-dual optimisation scheme can be used to generate very high-quality staggered, pair-wise orthogonal unstructured grids that are both more centroidal and better self-centred than standard Delaunay-Voronoi type approaches.

Turning, lastly, to the distribution of vertex weights generated by the coupled scheme, the relative power metric is used to colour the cells of the dual grids presented in Figures 8 and 9. These results show that the distribution of weights throughout each mesh is typically smooth, except in regions where the spacing of the grid is changed abruptly, or in the neighbourhood of topological inconsistencies. It appears that dual cells with a lower than optimal valence (i.e. 5-cells and smaller) are typically associated with larger, negative weights. Conversely, those with higher than optimal valence (i.e. 7-cells and larger) are typically associated with larger positive weights. Larger magnitude weights were also observed in the neighbourhood of boundary constraints, where it is often more difficult to optimise the geometry and/or topology of the grid via other means. Perfectly regular hexagonal configurations were found to induce more uniform distributions of weights; typically tending towards zero. Such behaviour is consistent with expectations; showing that local weight gradients are induced in patches of the grid where a conventional Delaunay-Voronoi pair would define sub-optimal configurations.

7 Conclusions & future work

A new optimisation-based algorithm designed to generate very high-quality staggered unstructured grids has been described, focusing on the development of techniques for the construction of weighted Regular-Power tessellations appropriate for unstructured co-volume type numerical formulations. Employing a combination of geometrical, topological and weight-based local optimisation operators, the new coupled primal-dual mesh optimisation schemes have been shown to generate generalised staggered unstructured grids that are simultaneously pair-wise orthogonal, centroidal and are typically well-centred. The performance of the new schemes have been tested experimentally, demonstrating that improved primal-dual complexes are generated compared to conventional optimised Delaunay-Voronoi type tessellations. Specifically, results confirm that the generalised Regular-Power pairs generated by the new optimisation algorithms are typically both more centroidal and are better self-centred.

A range of work is slated for future investigations, focusing on further improvements to the grid generation algorithm itself, as well as an application of the new methodology to various problems in ocean modelling and computational physics more broadly. Work on the primal-dual optimisation algorithm is proceeding in a number of directions, including: (i) an investigation of alternative, quasi-Newton based optimisation schemes chen2011efficient ; chen2014revisiting ; exploring whether further enhancements to grid quality can be achieved through the application of globally-coupled optimisation procedures, (ii) the generation of higher-dimensional primal-dual tessellations through an extension of the Regular-Power formulation to problems in , and (iii) improvements to the computational efficiency of the implementation, through the application of parallel programming patterns. Such work will be made available through future releases of the JIGSAW meshing library; currently under development by the author JIGSAW-GEO .

It is expected that the methods described here can be generalised to higher-dimensional problems, both for triangulations embedded in and for complexes of higher topological degree (i.e. tetrahedral meshes). In the case of general -simplex meshes, additional work is needed — firstly, requiring the definition of an appropriate dual quality metric and secondly, necessitating the development of methods to update general -simplex topologies. It is proposed that a straightforward extension of the dual metric (16) described here could be derived for higher degree elements, consisting of a weighted sum of the geometrical ‘defects’ associated with each -face per simplex. For example, a dual quality metric for tetrahedral complexes would consist of a weighted sum of the six edge-, four face- and single element-centred defect terms associated with each tetrahedron. It is proposed that optimisation of such a metric would reduce the geometrical offsets between the various -faces in a given primal-dual pair, leading to configurations that are more mutually ‘centred’.

It is expected that topological updates may be more challenging for higher-degree problems. Specifically, it is known ChengDeyShewchuk that algorithms based on local element flips are not guaranteed to recover the topologies of higher-dimensional Delaunay and/or Regular tessellations; potentially becoming trapped in local topological optima. It is proposed that a full re-triangulation would instead be computed at each iteration, based on, for example, efficient implementations of Bowyer-Watson type tessellation schemes Boissonnat02CGALTria . While such strategies may increase computational expense compared to the flipping-based scheme described here, they are not expected to fundamentally alter the composition of the resulting primal-dual optimisation scheme, which would consist of optimisation-driven vertex and weight updates, a global re-triangulaton step, and various edge- and face- collapse and refinement operations.

The use of generalised primal-dual structures as a basis for computational simulation is a relatively new area of interest, and a wide range of investigations are expected to be fruitful. It is expected that the application of optimal tessellations that are more centroidal and better self-centred than conventional Delaunay-Voronoi pairs may lead to improved numerical performance; enhancing the accuracy, stability and efficiency of the underlying formulations. Presently, the generation of optimal unstructured pairs for problems in geophysical fluid dynamics is a key focus; seeking to optimise the performance of the mimetic co-volume scheme employed in the Model for Prediction Across Scales (MPAS-O) ringler2008multiresolution ; ringler2010unified ; ringler2013multi for the simulation of ocean dynamics and climate processes. The application of generalised grids to other problems in computational physics is also of interest, including, for example, applicability to co-volume type formulations for the solution of Maxwell’s equations (e.g. yee1966numerical ; walton2017advances ), in addition to various methods formulated using Discrete Exterior Calculus.

8 Acknowledgements

This work was conducted at the Centre for Climate Systems Research (CCSR) at Columbia University and NASA’s Goddard Institute for Space Studies. The author wishes to thank Scott Mitchell for his thoughts on an initial draft, as well as the three anonymous reviewers for their helpful comments and feedback.

9 Code availability

The JIGSAW mesh generation package used in this study is available online:

10 References


  • (1) F. H. Harlow, J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, The physics of fluids 8 (12) (1965) 2182–2189.
  • (2) A. Arakawa, V. R. Lamb, Computational design of the basic dynamical processes of the UCLA general circulation model, Methods in computational physics 17 (1977) 173–265.
  • (3) P. Memari, P. Mullen, M. Desbrun, Parametrization of generalized primal-dual triangulations, in: Proceedings of the 20th International Meshing Roundtable, Springer, 2011, pp. 237–253.
  • (4) J. R. Shewchuk, Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator, in: Applied computational geometry towards geometric engineering, Springer, 1996, pp. 203–222.
  • (5) H. Si, TetGen, a Delaunay-based quality tetrahedral mesh generator, ACM Transactions on Mathematical Software (TOMS) 41 (2) (2015) 11.
  • (6) C. Jamin, P. Alliez, M. Yvinec, J.-D. Boissonnat, CGALmesh: a generic framework for Delaunay mesh generation, ACM Transactions on Mathematical Software (TOMS) 41 (4) (2015) 23.
  • (7) J. Schöberl, NETGEN: An advancing front 2D/3D-mesh generator based on abstract rules, Computing and visualization in science 1 (1) (1997) 41–52.
  • (8) C. Geuzaine, J.-F. Remacle, 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) (2009) 1309–1331.
  • (9) D. Engwirda, JIGSAW-GEO (1.0): locally orthogonal staggered unstructured grid generation for general circulation modelling on the sphere, Geoscientific Model Development 10 (6) (2017) 2117.
  • (10) F. Aurenhammer, Power diagrams: properties, algorithms and applications, SIAM Journal on Computing 16 (1) (1987) 78–96.
  • (11) H. Edelsbrunner, Algorithms in combinatorial geometry, Vol. 10, Springer Science & Business Media, 2012.
  • (12) T. Ringler, J. Thuburn, J. Klemp, W. Skamarock, A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured C-grids, Journal of Computational Physics 229 (9) (2010) 3065–3090.
  • (13) T. Ringler, M. Petersen, R. Higdon, D. Jacobsen, P. Jones, M. Maltrud, A multi-resolution approach to global ocean modeling, Ocean Modelling 69 (2013) 211–232.
  • (14) J. Thuburn, T. Ringler, W. Skamarock, J. Klemp, Numerical representation of geostrophic modes on arbitrarily structured C-grids, Journal of Computational Physics 228 (22) (2009) 8321 – 8335.
  • (15) P. Korn, Formulation of an unstructured grid model for global ocean dynamics, Journal of Computational Physics 339 (2017) 525 – 552.
  • (16) P. Korn, S. Danilov, Elementary dispersion analysis of some mimetic discretizations on triangular C-grids, Journal of Computational Physics 330 (2017) 156 – 172.
  • (17) W. C. Skamarock, J. B. Klemp, M. G. Duda, L. D. Fowler, S.-H. Park, T. D. Ringler, A multiscale nonhydrostatic atmospheric model using centroidal Voronoi tesselations and C-grid staggering, Monthly Weather Review 140 (9) (2012) 3090–3105.
  • (18) E. Vanderzee, A. N. Hirani, D. Guoy, E. Ramos, Well-centered planar triangulation–an iterative approach, in: Proceedings of the 16th International Meshing Roundtable, Springer, 2008, pp. 121–138.
  • (19) E. Vanderzee, A. N. Hirani, D. Guoy, E. A. Ramos, Well-centered triangulation, SIAM Journal on Scientific Computing 31 (6) (2010) 4497–4523.
  • (20) Q. Du, V. Faber, M. Gunzburger, Centroidal Voronoi tessellations: applications and algorithms, SIAM review 41 (4) (1999) 637–676.
  • (21) Q. Du, M. Gunzburger, Grid generation and optimization based on centroidal Voronoi tessellations, Applied Mathematics and Computation 133 (2) (2002) 591–607.
  • (22) P. Alliez, D. Cohen-Steiner, M. Yvinec, M. Desbrun, Variational tetrahedral meshing, in: ACM Transactions on Graphics (TOG), Vol. 24, ACM, 2005, pp. 617–625.
  • (23) L. Chen, J.-C. Xu, Optimal Delaunay triangulations, Journal of Computational Mathematics (2004) 299–308.
  • (24) L. Chen, M. Holst, Efficient mesh optimization schemes based on optimal Delaunay triangulations, Computer Methods in Applied Mechanics and Engineering 200 (9) (2011) 967–984.
  • (25) P. Mullen, P. Memari, F. de Goes, M. Desbrun, HOT: Hodge-optimized triangulations, ACM Transactions on Graphics (TOG) 30 (4) (2011) 103.
  • (26) F. d. Goes, P. Memari, P. Mullen, M. Desbrun, Weighted triangulations for geometry processing, ACM Transactions on Graphics (TOG) 33 (3) (2014) 28.
  • (27) S. Walton, O. Hassan, K. Morgan, Advances in co-volume mesh generation and mesh optimisation techniques, Computers & Structures 181 (2017) 70–88.
  • (28) S. Walton, O. Hassan, K. Morgan, Reduced order mesh optimisation using proper orthogonal decomposition and a modified cuckoo search, International Journal for Numerical Methods in Engineering 93 (5) (2013) 527–550.
  • (29) S. Walton, O. Hassan, K. Morgan, Selected engineering applications of gradient free optimisation using cuckoo search and proper orthogonal decomposition, Archives of Computational Methods in Engineering 20 (2) (2013) 123–154.
  • (30) D. Jacobsen, M. Gunzburger, T. Ringler, J. Burkardt, J. Peterson, Parallel algorithms for planar and spherical Delaunay construction with an application to centroidal Voronoi tessellations, Geoscientific Model Development 6 (4) (2013) 1353–1365.
  • (31) T. Ringler, L. Ju, M. Gunzburger, A multiresolution method for climate system modeling: application of spherical centroidal Voronoi tessellations, Ocean Dynamics 58 (5-6) (2008) 475–498.
  • (32) J. R. Shewchuk, Lecture notes on geometric robustness, Interpolation, Conditioning, and Quality Measures. In Eleventh International Meshing Roundtable.
  • (33) B. Perot, Conservation Properties of Unstructured Staggered Mesh Schemes, Journal of Computational Physics 159 (1) (2000) 58 – 89.
  • (34) P. S. Peixoto, S. R. Barros, On vector field reconstructions for semi-Lagrangian transport methods on geodesic staggered grids, Journal of Computational Physics 273 (2014) 185 – 211.
  • (35) P. S. Peixoto, Accuracy analysis of mimetic finite volume operators on geodesic grids and a consistent alternative, Journal of Computational Physics 310 (2016) 127–160.
  • (36) L. A. Freitag, C. Ollivier-Gooch, Tetrahedral mesh improvement using swapping and smoothing, International Journal for Numerical Methods in Engineering 40 (21) (1997) 3979–4002.
  • (37) B. M. Klingner, J. R. Shewchuk, Aggressive tetrahedral mesh improvement, in: Proceedings of the 16th international meshing roundtable, Springer, 2008, pp. 3–23.
  • (38) V. N. Parthasarathy, C. M. Graichen, A. F. Hathaway, A Comparison of Tetrahedron Quality Measures, Finite Elements in Analysis and Design 15 (3) (1994) 255–261.
  • (39) J. Shewchuk, What is a good linear finite element? interpolation, conditioning, anisotropy, and quality measures (preprint), University of California at Berkeley (2002) 70.
  • (40) G. S., O. C., Tetrahedral mesh generation using delaunay refinement with non‐standard quality measures, International Journal for Numerical Methods in Engineering 87 (8) (2011) 795–820.
  • (41) C. L. Lawson, Software for Surface Interpolation, in: J. R. Rice (Ed.), Mathematical Software III, Academic Press, New York, 1977, pp. 161–194.
  • (42) D. Engwirda, JIGSAW(GEO): Unstructured grid generation for geophysical modelling, (2017).
  • (43) C. Amante, B. W. Eakins, ETOPO1 1 arc-minute global relief model: procedures, data sources and analysis, US Department of Commerce, National Oceanic and Atmospheric Administration, National Environmental Satellite, Data, and Information Service, National Geophysical Data Center, Marine Geology and Geophysics Division Colorado, 2009.
  • (44) D. Engwirda, D. Ivers, Off-centre Steiner points for Delaunay-refinement on curved surfaces, Computer-Aided Design 72 (2016) 157–171.
  • (45) P.-O. Persson, Mesh Size Functions for Implicit Geometries and PDE-based Gradient Limiting, Engineering with Computers 22 (2) (2006) 95–109.
  • (46) Z. Chen, W. Wang, B. Lévy, L. Liu, F. Sun, Revisiting Optimal Delaunay Triangulation for 3D graded mesh generation, SIAM Journal on Scientific Computing 36 (3) (2014) A930–A954.
  • (47) S. W. Cheng, T. K. Dey, J. R. Shewchuk, Delauay Mesh Generation, Taylor & Francis, New York, 2013.
  • (48) J. D. Boissonnat, O. Devillers, S. Pion, M. Teillaud, M. Yvinec, Triangulations in CGAL, Comput. Geom. Theory Appl. 22 (1-3) (2002) 5–19.
  • (49) K. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Transactions on antennas and propagation 14 (3) (1966) 302–307.