Benchmark of Polygon Quality Metrics for Polytopal Element Methods
M. Attene, S. Biasotti, S. Bertoluzza, D. Cabiddu, M. Livesu, G. Patanè, M. Pennacchio, D. Prada, M. Spagnuolo
BBenchmark of Polygon Quality Metricsfor Polytopal Element Methods
M. Attene, S. Biasotti, S. Bertoluzza, D. Cabiddu, M. Livesu,G. Patanè, M. Pennacchio, D. Prada, M. Spagnuolo
Consiglio Nazionale delle RicercheIstituto di Matematica Applicata e Tecnologie InformaticheItaly
Abstract
Polytopal Element Methods (PEM) allow to solve differential equations on generalpolygonal and polyhedral grids, potentially offering great flexibility to mesh genera-tion algorithms. Differently from classical finite element methods, where the relationbetween the geometric properties of the mesh and the performances of the solver arewell known, the characterization of a good polytopal element is still subject to ongoingresearch. Current shape regularity criteria are quite restrictive, and greatly limit the setof valid meshes. Nevertheless, numerical experiments revealed that PEM solvers canperform well on meshes that are far outside the strict boundaries imposed by the currenttheory, suggesting that the real capabilities of these methods are much higher. In thiswork, we propose a benchmark to study the correlation between general 2D polygonalmeshes and PEM solvers. The benchmark aims to explore the space of 2D polygonalmeshes and polygonal quality metrics, in order to identify weaker shape-regularity cri-teria under which the considered methods can reliably work. The proposed tool is quitegeneral, and can be potentially used to study any PEM solver. Besides discussing thebasics of the benchmark, in the second part of the paper we demonstrate its applicationon a representative member of the PEM family, namely the Virtual Element Method,also discussing our findings.
Preprint submitted to Elsevier June 5, 2019 a r X i v : . [ c s . G R ] M a y . Introduction Solving a PDE on geometrically complex domains is a fundamental task of sci-entific computing. In real applications, the generation of a good discretization of thedomain is key to obtain quality results, but it can also be an extremely complex prob-lem, often even harder than the numerical solution of the discretized equations [15].To alleviate meshing issues, recent literature widely explored Polytopal Element Meth-ods (PEM), that is, methods for the numerical solution of PDEs based on polygonaland polyhedral grids. PEM approaches allow to: (i) achieve high flexibility in thetreatment of complex geometries; (ii) incorporate complex features at different scaleswithout triggering mesh refinement; (iii) automatically include hanging nodes (i.e., T-junctions); (iv) simplify refinement and coarsening operators.Similarly to standard Finite Elements, the performance of PEM depends on thequality of the underlying mesh, in terms of accuracy, stability, and effectiveness ofpreconditioning techniques. However, while the concept of shape regularity of classicalfinite elements is well understood [28, 6, 3], the characterization of a good polytopalelement is still subject to ongoing research.Recent works on shape regularity for PEM [18, 21] require each element in themesh to be star-shaped. Despite these tight constraints, numerical tests revealed thatPEM allows to reliably solve PDEs even on meshes made of bizarre non starred poly-gons, such as the ones shown in Figure 1. This suggests that there is a gap between thecurrent theory and the real capabilities of the solver, which promise to be much higher(thus allowing for much weaker requirements on the tessellation). However, unlikeclassical finite elements, the enormous freedom in terms of different shapes given bythe polytopal framework, makes it difficult to pinpoint exactly what geometric featuresdo or do not have a negative effect on the performance of the method.In this paper, we introduce a tool that allows to systematically explore the relationbetween the geometric properties of the mesh and the performances of PEM solvers.Our goal is twofold: first, we aim to identify more permissive shape-regularity criteriaunder which the considered method is effective; second, we wish to single out specificissues that negatively affect the results, with the final aim of designing better meth-2 igure 1: A PDE solved on a mesh made of snowflake-like elements (and their dual, that is elements withsnowflake holes). Despite the accuracy of the solution, these elements clearly violate any known shaperegularity criterion for polygonal elements. ods. More precisely, we propose a benchmark that studies the correlation between theperformance of a given PEM solver, and a wide set of polygon quality metrics.We considered a large set of geometric properties of polygons, spanning from areas,angles and edge lengths, to kernels, inscribed and circumscribed circles. The wholeset is reported in Table 1, and comprises 12 per polygon metrics, which become 36per mesh metrics considering minima, maxima, and averages. Concerning polygonalmeshes, available Voronoi based meshing tools (e.g. [11]) are not suited to aid ourstudy, because they produce convex elements which are not challenging enough tostress PEM solvers. Instead, we opted for a family of parametric elements explicitlydesigned to progressively stress one or more of the metrics listed in Table 1, enrichedwith random polygons to avoid the risk to bias the study with artificially constructedelements. Finally, for the solver we considered both the accuracy of the solution andthe conditioning of the associated linear system. The whole framework can be appliedto any PEM method, and is highly modular, thus favoring further extensions with newpolygons or metrics. 3o the best of our knowledge, this work is the first systematic study on the correla-tion between polygon quality metrics and numerical methods to solve PDEs on discretepolygonal grids. As such, we believe it will have great practical impact on two commu-nities. For the geometry processing community, metrics with a good correlation withthe solver could be incorporated into new meshing tools, producing polygonal gridsand refinement/coarsening operators that leverage all the flexibility granted by PEM,thus opening the door to new PDE-aware mesh processing tools. For the math com-munity, correlations will suggest directions to further improve the theory, developingmore permissive shape-regularity criteria for the current methods, and new methods toovercome the limitations of the current solvers.The outline of the paper is as follows. In Section 2, we recall classical resultsconcerning shape regularity criteria for FEM, and discuss how these concepts are ad-dressed in the PEM framework. In Section 3 we provide details on the structure andorganization of the benchmark. In Section 4 we use the benchmark to study a PEMsolver based on the Virtual Element Method, and discuss our findings. In Section 5,we suggest some directions of future research. Finally, in Appendix Appendix A, webriefly recall the definition and the main properties of the specific PEM solver we usedfor our numerical experiments.
2. Related works
As already mentioned in the introduction, the influence of the shape of the ele-ments on the classical finite element method is quite well understood [6]. In particularthe interested reader can refer to [3], where the equivalence between different shaperegularity criteria for triangular and tetrahedral elements is studied.The PEM family already counts quite a few methods, such as Mimetic Finite Dif-ferences [36, 4], Discontinuous Galerkin-Finite Element Method (DG-FEM) [1, 5],Hybridizable and Hybrid High-Order Methods [7, 10], Weak Galerkin Method [37],BEM-based FEM [24], Poly-Spline FEM [25], and Polygonal FEM [31], to name afew. 4 d Figure 2: Minimal polygon shape regularity is currently expressed in terms of the ratio between the maximalball inscribed in the kernel of the polygon ( d ), and the maximal ball inscribing the element ( D ). Shape-regularity for VEM..
As a representative example of the PEM family of ap-proaches, one may consider the Virtual Element Method (VEM) [33], which can beseen as an extension to FEM for handling general polytopal meshes. Current a-prioriestimates for VEM admit only star-shaped polytopes, and connect the approximationerror with the ratio between the radius of the circumscribed circle and the radius of thecircle inscribed in the kernel of the element (Figure 2). Moreover, the optimality of themethod is only achieved in an even stricter framework (see Assumption Appendix A.2in Appendix Appendix A).
PDEs in Computer Graphics..
Solutions of differential equations have been exten-sively used in computer graphics for a variety of applications, such as mesh param-eterization [12], computation of diffusion distances [22], Voronoi diagrams [14], andsmoothing [8, 19]. However, the problem of assessing the capabilities of a numericalmethod when paired with a particular mesh has been considered only in very recentyears. Restricting to the CG community, the work that is closely related to us is proba-bly [13], where a study on the FEM performances on hexhaderal meshes according tovarious geometric metrics is proposed. For this work, the authors could leverage both arich database of hexmeshes produced with various approaches (subsequently released5
NALYSIS SOLUTIONS(RAW DATA)
Ground Truth............... Numerical...............
PEM SOLVERMETRICS
PEM Solver
L2LinfK
Geometric
IC, CC, CR,AR, KE, KAR,PAR, MA, SE,ER, MPD, NPD
MESHESCODE/SCRIPTS
C++............ Matlab............ Bash............
POLYGONS
ParametricRandom
Figure 3: Our benchmark (gray shaded area) and how it relates with the PEM solver. in [2]), and well established quality metrics for hexahedral cells [30]. The benchmarkwe propose can be seen as its extension to the PEM case. However, this extension isfar from trivial, mostly because for the polygonal case the are no available resources interms of mesh databases, and no consensus has yet been reached on what are the rightmetrics to consider to evaluate the mesh quality with respect to the performances of thesolver. Schneider et al. [26] identify in mesh resolution, element quality, and basis or-der the three factors that affect the accuracy of standard FEM, and advocate the use ofhigher order basis to compensate badly shaped elements. PEM approaches promise tobe robuster than FEM against poor elements (Figure 1), and the benchmark we proposeaims to study how much geometric freedom one can gain, for a fixed basis order. Re-berol and Levy [23] propose an efficient voxel-based method to compare two solutionscomputed on alternative meshes of the same object. In our benchmark ground truthand numerical solution are always computed in the same mesh, therefore sophisticatedtechniques to transfer a solution from a mesh to another are not necessary.
3. Benchmark
We present here the core of our contribution, that is a full benchmark that allowsto investigate the correlation of a set of per polygon quality measures ( Γ geom ) withthe performances of a PEM method ( Γ pde ). The relationship between these entities is6omputed on a carefully crafted set of discrete polygonal meshes ( Γ mesh ). The maincomponents of the benchmark are therefore the three sets Γ geom , Γ pde , Γ mesh , the detailsof which are described in Sections 3.1, 3.2 and 3.3, respectively.Note that the tool we offer is aimed to be quite general, and can be used to evaluatepotentially any PEM solver. To this end, we point out that the solver is not part of thebenchmark, but is rather seen as a black box that takes in input a mesh, and returns thesolution of a PDE, as depicted in Figure 3.The major output of the benchmark is a set of measures, both related to the meshitself and to the PDE solver, computed on the meshes Γ mesh . Output data is providedto the user both as raw ASCII files and in statistical format. Statistical data amounts toa correlation analysis between the polygon quality measures in Γ geom and the perfor-mances of the solver, measured according to the metrics in Γ pde . Technical details onhow the statistical analysis is performed are given in Section 3.4. Such analysis comesin the form of a standard color-coded square matrix (Matlab scripts to load the data andcompute the analysis are included in the benchmark).As an additional result, the benchmark also outputs correlation analyses betweenthe geometric metrics in Γ geom , and between the performance estimators in Γ pde . Es-pecially for the geometric side, these correlations emphasize any possible redundancyin the set of metrics considered in the study, and therefore offer multiple optimizationchoices for the design of new meshing tools. To better understand this concept, let usassume that two geometric metrics, m and m , have a strong correlation to each other,and also assume that metric m has a strong correlation with the solver performances(e.g., accuracy, stability). Then, meshing algorithms may optimize either for m or for m , obtaining similar (positive) results. Considering that some geometric metrics aremuch more challenging to optimize than others, having such redundancy is of muchpractical utility.For a concrete example of correlation analysis produced with our benchmark, thereader can refer to Section 4, where the Virtual Element Method (VEM) is used as ablack box PEM solver. 7 C CCKEAR MASEMPD { Figure 4: Polygon measures we consider in our study: inscribed circle (IC); circumscribed circle (CC); poly-gon area (AR); kernel area (KE); minimum angle (MA); shortest edge length (SE); and minimum point topoint distance (MPD). Relative metrics obtained computing ratios between these quantities are also consid-ered. The full list is available at Table 1.
For the geometric part, our study is based on the measures depicted in Figure 4,which are considered both alone and combined to one another in order to obtain scaleinvariant metrics. The resulting full list of 12 per polygon metrics is presented in Ta-ble 1. Considering minimum, average, and maximum values of each metric, we obtaina total of 36 per mesh metrics. The benchmark provides C++ code to evaluate thesemetrics on polygonal meshes, as well as scripts to conveniently run these programs inbatch mode on entire datasets of domains.Given a polygon P , the list of per polygon metrics is the following: • Inscribed Circle (IC): it measures the radius of the biggest circle fully containedin P . Computing the maximum inscribed circle of a general polygon is a thoughgeometric problem. We start from a Voronoi diagram of the edges of P , and8 igure 5: Left: a polygon that does not admit a circumscribed circle. To be able to scale on general polygons,we define the radius of the circumscribed circle (CC) as the radius of the smallest circle containing thepolygon itself. Right: for a skinny triangle, CC is therefore the radius of the smallest circle that passesthrough the endpoints of its longest edge (green), and not of the circle passing through all its three vertices(red). select as center of the circle the corner in the diagram that is furthest from alledges. The minimum distance between such point and any of the edges of P isthe radius of the IC. Note that, differently from the point case, the diagram of aset of segments has curved boundaries between cells. For its computation, werely on the Boost Polygon Library [29]; • Circumscribed Circle (CC): it measures the radius of the smallest circle fullycontaining P . We compute it by treating the vertices of P as a point cloud, andthe running the Welzl’s algorithm to solve the minimum covering circle prob-lem [38]. Note that this does not correspond to the classical definition of circum-scribed circle, which does not necessarily exist for all polygons (Figure 5); • Circle Ratio (CR): it is the ratio between IC and CC. Differently from the previ-ous two, this measure does not depend on the scale of the polygon, and is alwaysdefined in the range ( , ] ; 9 etric Abbr. Range Trend Scaleinvariant Inscribed circle (radius) IC ( , ∞ ) – NoCircumscribed circle (radius) CC ( , ∞ ) – NoCircle ratio (IC/CC) CR [ , ] ↑ YesArea AR [ , ∞ ) – NoKernel area KE [ , ∞ ) – NoKernal area ratio (KE/AR) KAR [ , ] ↑ YesPerimeter area ratio PAR ( , ∞ ) ↓ YesMinimum angle MA ( , π ) ↑ YesShortest edge SE ( , ∞ ) – NoEdge ratio ER ( , ] ↑ YesMin point to point distance MPD ( , ∞ ) – NoNormalized point distance NPD ( , ] ↑ Yes
Table 1: List of 12 per polygon metrics used in our study. The metrics become 36 for a full polygonal mesh(minimum, maximum and averages of each metric are considered). For scale invariant measures, the fourthcolumn indicates whether optimal values are at the top ( ↑ ) or bottom ( ↓ ) of the definition range. • Area (AR): simply the area of the polygon P ; • Kernel Area (KE): it is the area of the kernel of the polygon. The kernel of P isdefined as the set of points p ∈ P from which the whole polygon is visible. If thepolygon is convex, the area of the polygon and the area of the kernel are equal.If the polygon is star-shaped, the area of the kernel is a positive number. If thekernel is non star-shaped, the polygon has no kernel and KE will be zero; • Kernel Area ratio (KAR): it is the ratio between the area of the kernel of P and itswhole area. For convex polygons, this ratio is always 1. For concave star-shapedpolygons, KAR is strictly defined in between 0 and 1. For non star-shaped poly-gons, KAR is always zero; 10 Area Perimeter Ratio (APR): often referred to as compactness of P , it is definedas 2 π ∗ area ( P ) perimeter ( P ) . This measure reaches its minimum for the most compact 2D shape (the circle),and grows for less compact polygons; • Minimum angle (MA): it is defined as the minimum inner angle of the polygon P ; • Shortest Edge (SE): it is the length of the shortest edge of P ; • Edge Ratio (ER): it is the ratio between the shortest and the longest edge of P ; • Minimum point to point distance (MPD): it is the minimum distance betweentwo (possibly non consecutive) points in P . Note that in case the two closestpoints in P are also consecutive, then MPD and SE are equivalent; • Normalized Point distance (NPD): it is the normalized version of MPD. To makethe shortest point to point distance independent from the scale of the polygon, weuse the diameter of the minimum circumscribed circle as normalization factor.This bounds NPD within the definition range ( , ] . To measure the performance of a PDE solver we compute the following quantities,where u and u h are the ground-truth solution and the solution computed with a PEMsolver, respectively: • Relative L ∞ -error : ε ∞ : = (cid:107) u − u h (cid:107) ∞ / (cid:107) u (cid:107) ∞ , (cid:107) v (cid:107) ∞ = max x ∈ Ω v ( x ) ; • Relative L -error : ε : = (cid:107) u − u h (cid:107) / (cid:107) u (cid:107) , where (cid:107) v (cid:107) = (cid:0) (cid:82) Ω v ( x ) d x (cid:1) / ; • Relative error in the discrete energy norm : ε S : = ( (cid:107) u − u h (cid:107) S ) / (cid:107) u (cid:107) S , with (cid:107) v (cid:107) S = v T S v , where S is the global stiffness matrix arising from a PEMdiscretization, and v is the vector collecting all the degrees of freedom of v ; • Multiplicative constant, C, and estimated convergence rate, k , for the relative L -error: for most PEM, it can be proved that ε = Ch p , (1)where C is a constant depending on the shape of the elements of the tessellation, h is the maximum point-to-point distance across all the elements, and p = ε ≈ log C + p log h ;therefore, given a sequence of M mesh refinements, C and p can be estimatedusing linear regression on the pairs ( log h , log ε ) i , i = , . . . , M ; • L condition number of the stiffness matrix : κ ( S ) = (cid:107) S (cid:107) (cid:107) S − (cid:107) , which represents an indicator of the numerical stability of the PEM solver. Thehigher κ , the less stable the PEM solver. Condition number heavily impacts12he performances of iterative solvers, and is therefore important especially whensolving big problems that cannot be plugged in a direct solver. In order to evaluate the dependence of the performance of a PEM solver on thegeometrical features of the underlying mesh, we propose a set of polygonal meshesexplicitly designed to progressively stress the metrics listed in Table 1. We achievedthis goal by generating a family of parametric polygons, P ( t ) , with t ∈ [ , ] . Eachpolygon has a baseline configuration that does not present critical geometric features( P ( ) ), and is progressively made worse by a deformer, controlled by the parameter t . We designed 8 different parametric polygons, depicted in Figure 6. For each poly-gon, we report both its evolution for growing values of t , and the geometric metricsit affects (Table 2). Note that since the metrics are not independent of one another,each polygon family controls multiple metrics. To minimize any bias due to artificiallyconstructed polygons, we enriched the family of elements with randomly generatedpolygons, created with CGAL [32]. A few examples are shown in Figure 7.Given both the parametric and random polygons, the 2D polygonal domains arecreated by placing an element at the center of a squared canvas, and filling the restof the domain with triangles, using [27]. Since we are interested in the response ofthe PEM solver with respect to the parametric and random polygons, we aim to fillthe rest of the domain with well shaped elements. We obtain the desired effect byimposing 20 ◦ as minimum inner angle for each triangle. Note that quality constraintsare satisfied only away from the central polygon, whereas triangles directly incident toit can be arbitrarily badly shaped, depending on its shape.The benchmark contains C++ code to generate single polygonal elements as wellas random elements. The program inputs the name of the polygon family, the value ofparameter t , and a binary flag to enable/disable the generation of the triangular canvassurrounding the polygon. Scripts to conveniently sample the parameter space and cre-ate a whole dataset with a single instruction are also included in the software package.To complete the mesh generation tools, the benchmark offers the possibility tocreate a full multi-resolution hierarchy, obtained by mirroring the squared domains13 =0 (rest configuration) t=0.25 t=0.5 t=0.75 t=1.0 C o m b C o n v e x i t y I s o t r o p y M a ze N s i d e dS t a r U L i k e Z e t a Figure 6: The eight families of parametric polygons used to generate the meshes in the dataset. Each polygonstarts from a rest configuration ( t =
0, left column), and progressively worsen for growing values of t ,stressing either one or multiple quality measures listed in Table 1. igure 7: A few examples of randomly generated polygons considered in our benchmark. As can be noticedthere is a big variety in terms of complexity of the shapes, spanning from simple convex polygons (top right)to extremely challenging ones (top left). multiple times. The users can choose how many levels in the hierarchy they want tocreate (e.g. to control average edge length and test the convergence of the solver). Afew example of polygonal meshes are shown in Figure 8. Figure 8: An example of 2D polygonal mesh used in the benchmark (left). A mesh hierarchy is also obtainedby mirroring the meshes, thus reducing the average edge length and enabling convergence checks (right).
To measure the correlation among the geometric metrics and/or performance met-rics, the benchmark offers MATLAB R (cid:13) scripts to load data, measure and visualize the15 o m b C o n v e x i t y Is o t r o p y M a ze N S i d e d S t a r U L i k e Z e t a IC (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) CCCR (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) AR (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) KE (cid:88) (cid:88) (cid:88) (cid:88) KAR (cid:88) (cid:88) (cid:88) (cid:88) (cid:88)
PAR (cid:88) (cid:88) (cid:88) (cid:88) MA (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) SE (cid:88) (cid:88) (cid:88) ER (cid:88) (cid:88) MPD (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88)
NPD (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88)
Table 2: Correlation between metrics per polygon and the parametric polygons used to generate the meshesin the benchmark. well-known Pearson correlation coefficient (see, for instance, [9]), defined as: ρ a , b = ∑ ni = ( a i − ¯ a )( b i − ¯ b ) (cid:112) ∑ ni = ( a i − ¯ a ) ∑ ni = ( b i − ¯ b ) where a i and b i are two observations for the mesh i (measures of either geometric orsolver quality in our case), n is the number of meshes considered, and ¯ a = ∑ n x i is themean of the observations done (analogous for ¯ b ). The correlation coefficient ρ ab ∈ [ − , ] measures the linear correlation between a and b . A value of 1 implies perfectlinear correlation, that is, that there exists a linear equation that perfectly describes therelationship between a and b , with all data points lying on a line for which the secondobservation increases as the first one increases. A value of minus 1 implies that all datapoints lie on a line for which the second observation decreases as the first one increases.A value of 0 implies that there is no linear correlation between these observations.16 .5. Extensions The whole benchmark is extremely modular, and its components can be easily ex-tended according to the user needs. Specifically, all the parametric polygons that popu-late the set Γ mesh implement the same virtual C++ class. Therefore, extending the fam-ily with a new polygon amounts to simply define its base shape, and then implementthe function that deforms it according to the value of parameter t ∈ [ , ] . Similarly,new geometric metrics can be added to the system by minimally extending the code-base. For the correlation analysis, MATLAB R (cid:13) routines to load and process the dataare included in the software package, and more complex analysis can be performedleveraging the capabilities of statistical toolboxes.
4. Benchmarking the Virtual Element Method
In this section we present the results obtained by coupling the benchmark describedin the previous section with a specific PEM solver. Specifically, we used a MATLAB R (cid:13) implementation of the lowest order Virtual Element Method (VEM) [20]. For a briefintroduction on VEM, the reader may refer to Appendix Appendix A and referencestherein.As a preliminary step, we investigate the existence of correlations between differentgeometric metrics (Section 4.4). Subsequently, we consider the correlation between ge-ometric and solver performance metrics (Section 4.5). For both these studies, we haverestricted our investigation to scale-invariant geometric metrics, namely, to CR, KAR,PAR, MA, ER, NPD (see Section 3.1 for details). We do not consider scale-dependentmetrics because we are interested in the effect of the shape of the elements only (theinfluence of element size on the performance of the PEM solver is well understood). In our numerical tests, we let the computational domain Ω be the unit square [ , ] × [ , ] , and consider the following model problem − ∆ u = f , f taken so that u ( x , y ) : = sin ( π x ) sin ( π y ) π is the exact solution. The MATLAB R (cid:13) VEM black-box prototype returns in outputan estimate of the condition number of the stiffness matrix, and of the relative errorbetween the exact and the numerical solution under various norms, as detailed in Sec-tion 3.2. This raw data is then processed with the scripts provided by the benchmark,feeding the correlation analysis that aims to study both the accuracy and numericalstability of the solver with respect to the various polygonal meshes generated by thebenchmark (Sections 4.4 and 4.5), identifying what geometric parameters influence thesolver performance.
For each of the eight base polygons, we picked 20 equally spaced samples in itsparameter domain [0, 1]. This produced 160 different polygons, classified in eightfamilies. Each of these polygons was used to generate a mesh for the unit square do-main, as detailed in Section 3.3. The 160 resulting meshes constitute our basic dataset D . Furthermore, the domain mirroring approach described in Section 3.3 was used togenerate four refined versions of D , that we call D , D , D and D , respectively. Notethat each D i has four times more polygons than D i − , and same ratio between numberof triangles and general polygons. Following a similar construction, we have also gen-erated a random dataset R composed of 100 meshes, as well as its associated level ofdetail mirrored hierarchy R , R , R , R . All the aforementioned meshes were gener-ated by using the scripts offered by the benchmark. Bigger (or smaller) datasets can beeasily obtained by asking such scripts to sample the space of parametric polygons andthe space of random polygons with a different frequency. The 6 metrics introduced in Section 3.2 are used to evaluate different aspects of thePEM solver. This analysis aims at studying the accuracy and numerical stability of theVEM solver defined on the different polygonal meshes of Ω , in order to identify which18 o n v e x i t y C o m b I s o t r o p y M a z e N s i d e d S t a r U li k e Z e t a R a nd o m Figure 9: ( y -axis) Solver metrics evaluated on the polygonal meshes ( x -axis) of the benchmark. geometric parameters of the meshes can influence the VEM solver. The results of ouranalysis are shown in Figure 9.We observe that the condition number explodes as the quality degenerates for thepolygons in the classes Convexity, Isotropy, Maze, U-like and
Zeta , which are the poly-gons with vanishing area for t =
1. While this behavior for these classes is to beexpected, it is somewhat surprising that the condition number does not explode for theclasses
Comb , Star and
N sided , which lead us to formulate the hypothesis that themetrics ER and NPD have little influence on this performance parameter.As far as the energy norm of the relative error is concerned, all the classes are af-fected to different degrees.
Isotropy , N sided and
Star are affected only minimally,while the most affected are
Maze and
U like , where a steep degradation of the errorcan be observed for very degenerate elements (surprisingly, the error for the
Zeta classdoes not seem to present the same behaviour, which deserves further investigation).Particularly interesting are the results for the
Convexity class, which displays the worstbehaviour as far as the condition number is concerned, but where the blowup-like be-haviour that we can observe in the two previously mentioned classes, does not occur.This leads us to think that the bad behaviour of these latter classes is indeed resultingfrom a bad performance of the VEM method itself, and it is not only a consequence ofthe bad condition number of the stiffness matrix.Finally we used the mirrored data set to estimate the convergence rate p and mul-19a) C o n v e x i t y C o m b I s o t r o p y M a z e N s i d e d S t a r U li k e Z e t a R ando m (b) C o n v e x i t y C o m b I s o t r o p y M a z e N s i d e d S t a r U li k e Z e t a R ando m Figure 10: Behaviour of the multiplicative constant of the (a) relative L -error and (b) estimated convergencerate of ε . tiplicative constant C of ε . The results, shown in Figure 10, are in accordance withthe theoretical estimates ( p =
2) also for the very degenerate cases and confirm therobustness of VEM already observed in the literature.
Before analyzing the correlation of the different mesh quality metrics with the per-formance metrics for the VEM method, we start by analyzing the correlations amongthe different mesh quality metrics. We focus on the non-triangular elements, thus wemeasure the metrics on these polygonal elements only. Figure 11 reports the correlationrelationship for both D and R , on the left and on the right respectively. As discussedin Section 3.4, two measurements positively correlate if ρ =
1, negatively correlate if ρ = − ρ =
0. When dealing with real data, a threshold is oftenconsidered to assess whether two measurements strongly correlate or not. In our set-ting, we assume { . , − . } the thresholds for strong positive and negative correlation,respectively, and { . , − . } the thresholds for weak positive and negative correlation,respectively. Considering the parametric polygons (Figure 11, left) we see that there isnot strong correlation among the quality metrics; only a few of them exhibit a moderatepositive correlation, while no negative correlation is present. Stated differently, the se-lected six metrics are rather independent measures. Figure 11 (right) reveals that, whenthe random polygons are measured, the correlation grows slightly (both in positive andnegative sense) with respect to the previous test. This behaviour is not unexpected be-cause random polygons are, on average, farther from the extreme configurations whencompared to the parameterized polygons. As a consequence, the measures cover a20 igure 11: Correlation among the six geometric metrics over the polygon elements of the dataset D (left)and over the polygon elements of the random dataset R (right). Legend: strong correlation, ρ ∈ ( . , ] ;weak correlation, ρ ∈ ( . , . ] ; strong inverse correlation, ρ ∈ [ − , − . ) ; weak inverse correlation, ρ ∈ [ − . , . ) . smaller range, which induces a higher correlation. Nonetheless, the two matrices inFigure 11 do not contradict each other. Finally, we have repeated the same experimentson the refined domains { D i } i = , and their random counterparts { R i } i = . We do notshow the results of these experiments as they are essentially identical to those on theunrefined domains. Even in this case, we have measured the six scale-independent metrics for all thepolygonal meshes in D . However, here we have measured the metrics both on thepolygonal element and over the triangular elements in its complement, and we haveselected the worst over all these values. Also, for each mesh we have run our VEM im-plementation and have compared the solution with the ground-truth using six differenterror metrics. By analyzing the correlation between geometric metrics and accuracymetrics on the set D we have found that (see Figure 12, left side): • KAR and MA are strong positively correlated with κ ( S ) , ε ∞ , ε , and ε S ; • CR is strongly positively correlated with κ ( S ) , ε , and ε S . And, though at aminor extent, with ε ∞ ; • ER is strongly positively correlated with ε and ε S ;21 C p K errinf errL2 errS C p K errinf errL2 errS [t]
Figure 12: Correlation among the geometric and accuracy metrics for the dataset D (left) and R (right).Legend: strong correlation, ρ ∈ ( . , ] ; weak correlation, ρ ∈ ( . , . ] ; strong inverse correlation, ρ ∈ [ − , − . ) ; weak inverse correlation, ρ ∈ [ − . , . ) . • APR and NPD are strongly negatively correlated with C .The same experiment was run on our random dataset R (see Figure 12, right side)and the results, though not contradicting the previous correlations, are less evident dueto smaller coverage of the geometric metrics. As for the geometry-geometry correlationexperiment, we have repeated the same tests on the refined domains { D i } i = , and ontheir random counterparts { R i } i = . The results of these experiments are essentiallyidentical to those on the unrefined domains.
5. Conclusions and future works
We have presented a benchmark to study the correlation between the geometricproperties of polygonal meshes, and the performances of Polytopal Element Methods(PEM) in 2D. The benchmark comprises meshes, metrics, code, and various scripts tosynthesize, load and process data, and can be coupled with any PEM solver.To the best of our knowledge this is the first tool to systematically study the cor-relation between polygonal meshes and PEM solvers. Considering the numerous ap-proaches of this kind recently introduced in literature [1, 5, 10, 37, 24, 25], we expectwide adoption of our tool.In Section 4 we used the benchmark to evaluate a specific PEM technique: theVirtual Element Method. The results of this study confirmed the robustness of the22ethod in terms of the L error accuracy, whereas the accuracy in the energy norm isaffected by some of the mesh quality metrics. In particular, some specific geometricmetrics (e.g. ER, NPD) seem to poorly affect the performances of VEM, whereasa specific class of polygons (i.e. Convexity ) revealed some weak spot to be furtherinvestigated, suggesting the the geometric metrics involved in the evolution of thatparticular shape, may be a good shortlist of polygon measures to deeply study withinthe VEM framework, also from a theoretical standpoint.
As mentioned in Section 3.5 the benchmark is modular, and can be extended in anumber of ways. At PEM level, we are obviously interested in testing the performanceof new solvers. For example, we are currently investigating the mutual influence ofshape regularity criteria and increasing polynomial degree. Moreover, the results ofour analysis, provide us with several interesting research directions relative to the the-oretical analysis of the method, such as the study of the effect, or lack thereof, of thepresence of small edges on the condition number.At mesh level, an interesting direction for future improvements consists in study-ing the space of 2D tilings [17, 16], and try to match it with the concept of parametricpolygons used in the benchmark. This would allow to create meshes fully made ofpolygons subject to study, without having to fill the canvas with triangles. Regardingthe correlation analysis, at the moment our study is limited to simple linear correla-tion. For future versions of the benchmark we plan to extend this part, and introducemore sophisticated (e.g. non-linear) correlation analysis. Finally, we are deeply inter-ested in extending the benchmark to the volumetric case, to study the performancesof PEM solvers also on general polyhedral meshes in 3D. Such extension is far fromtrivial though, as the number of geometric measures and associated metrics to be con-sidered is likely to be considerably higher. Same goes for the mesh generation phase,where densely sampling the space of parametric 3D polytopes may potentially lead toa combinatorial explosion. 23 ppendix A. The Virtual Element Method (VEM)
We briefly recall the definition and the main properties of the lowest order VirtualElement Method [33]. To fix ideas, we focus on the following elliptic model problem: − ∆ u = f in Ω , u = ∂Ω . (A.1)where Ω ⊂ R is a polygonal domain and f ∈ L ( Ω ) .We consider a family {T h } h of tessellations of Ω into a finite number of simplepolygons K , and let E h be the set of edges e of T h .For each polygon K ∈ T h , let the local space V K be defined as: V K = (cid:8) v ∈ H ( K ) : v | ∂ K ∈ C ( ∂ K ) , v | e ∈ P ( e ) ∀ e ∈ E K , ∆ v = K } , and the global virtual element space V h as V h = { v ∈ V : w | K ∈ V K ∀ K ∈ T h } . As degrees of freedom, uniquely identifying a function v h ∈ V h , we consider the valuesof v h at the vertices of the tessellation.Let ( · , · ) be the scalar product in L , a ( u , v ) = ( ∇ u , ∇ v ) then, using a Galerkinapproach, we would look for u h ∈ V h such that for all v h ∈ V h a ( u h , v h ) = (cid:90) Ω f v h dx . (A.2)However, both terms at the right and at the left hand side cannot be computed exactlywith the knowledge alone of the values of the degrees of freedom of u h and v h . On theother hand, setting, for each K ∈ T h a K ( u h , v h ) = (cid:90) K ∇ u h · ∇ v h dx we observe that, by using Green’s formula, given any v h ∈ V K and any p ∈ P ( K ) ⊆ V K a K ( p , v h ) = (cid:90) ∂ K v h ∂ p ∂ n . K v h is a known linear and ∂ p / ∂ n is a known constant, theright hand side can be computed exactly and directly from the value of degrees of free-dom of v h and p . This allows to define the “element by element” exactly computableprojection operator Π ∇ : V K −→ P ( K ) a K ( Π ∇ u h , q ) = a K ( u h , q ) ∀ q ∈ P ( K ) , and we clearly have a K ( u h , v h ) = a K ( Π ∇ u h , Π ∇ v h ) + a K ( u h − Π ∇ u h , v h − Π ∇ v h ) . The Virtual Element method stems from replacing the second term of the sum on theright hand, that cannot be computed exactly, with any computable symmetric bilinearform S satisfying for all v h with Π ∇ v h = c a K ( v h , v h ) ≤ S ( v h , v h ) ≤ c a K ( v h , v h ) , for two positive constants c and c , resulting in defining a Kh ( u , v ) = a K ( Π ∇ u , Π ∇ v ) + S ( u − Π ∇ u , v − Π ∇ v ) . The virtual element discretization of (A.2) yields the following discrete problem:
Problem Appendix A.1.
Find u h ∈ V h such thata h ( u h , v h ) = f h ( v h ) ∀ v h ∈ V h with a h ( u h , v h ) = ∑ K a Kh ( u h , v h ) .Problem Appendix A.1 is usually analyzed under the following assumption on thepolygons of the tessellation. Assumption Appendix A.2.
There exist constants γ , γ > such that (i) each element K ∈ T h is star-shaped with respect to a ball of radius ≥ γ h K , whereh K is the diameter of K; (ii) for each element K in T h the distance between any two vertices of K is ≥ γ h K . S [35]), including the simplest one (which weused in the numerical tests) of defining S in terms of the vectors of local degrees offreedom as the properly scaled euclidean scalar product.For details on the implementation as well as for the study of the convergence, sta-bility and robustness properties of the method we refer to [34, 33].Care is needed when computing the relative errors ε ∞ and ε for VEM. Since theanalytic expression of a VEM function is not known, we can only use its degrees offreedom and all the information we can deduce from them. In particular, for the lowestorder VEM, the degrees of freedom are the pointwise values at the mesh vertices P = { p i } i , i = , . . . , N vertices. Thus, to compute ε ∞ and ε , we define (cid:107) u − u h (cid:107) ∞ : = max i | u ( p i ) − u h ( p i ) | , and replace (cid:107) u − u h (cid:107) with (cid:107) u − Π ∇ u h (cid:107) . Acknowledgements
This work has been supported by the EU ERC Advanced Grant CHANGE, grantagreement No. 694515.
References [1] Antonietti, P.F., Cangiani, A., Collis, J., Dong, Z., Georgoulis, E.H., Giani, S.,Houston, P., 2016. Review of discontinuous Galerkin finite element methodsfor partial differential equations on complicated domains, in: Connections andChallenges in Modern Approaches to Numerical Partial Differential Equations.Springer, pp. 279–308.[2] Bracci, M., Tarini, M., Pietroni, N., Livesu, M., Cignoni, P., 2019. Hexalab.net:an online viewer for hexahedral meshes. Computer-Aided Design 110.[3] Brandts, J., Korotov, S., Kˇrížek, M., 2008. On the equivalence of regularity cri-teria for triangular and tetrahedral finite element partitions. Computers & Math-ematics with Applications 55, 2227 – 2233.264] Brezzi, F., Lipnikov, K., Shashkov, M., 2005. Convergence of mimetic finite dif-ference methods for diffusion problems on polyhedral meshes. SIAM J. Numer.Anal. 43, 1872–1896.[5] Cangiani, A., Georgoulis, E.H., Houston, P., 2014. hp-version discontinuousGalerkin methods on polygonal and polyhedral meshes. Math. Mod. Meth. Appl.Sci. 24, 2009–2041.[6] Ciarlet, P., 2002. The Finite Element Method for Elliptic Problems. Society forIndustrial and Applied Mathematics.[7] Cockburn, B., Dong, B., Guzmán, J., 2008. A superconvergent LDG-hybridizableGalerkin method for second-order elliptic problems. Math. Comp. 77, 1887–1916.[8] Desbrun, M., Meyer, M., Schröder, P., Barr, A.H., 1999. Implicit fairing of ir-regular meshes using diffusion and curvature flow, in: Proceedings of the 26thannual conference on Computer graphics and interactive techniques, Citeseer. pp.317–324.[9] Deza, M.M., Deza, E., 2009. Encyclopedia of Distances. Springer Berlin Heidel-berg.[10] Di Pietro, D.A., Ern, A., 2015. Hybrid high-order methods for variable-diffusionproblems on general meshes. Comptes Rendus Mathèmatique 353, 31–34.[11] Du, Q., Faber, V., Gunzburger, M., 1999. Centroidal voronoi tessellations: Ap-plications and algorithms. SIAM review 41, 637–676.[12] Floater, M.S., Hormann, K., 2005. Surface parameterization: a tutorial and sur-vey, in: Advances in multiresolution for geometric modelling. Springer, pp. 157–186.[13] Gao, X., Huang, J., Xu, K., Pan, Z., Deng, Z., Chen, G., 2017. Evaluating hex-mesh quality metrics via correlation analysis, in: Computer Graphics Forum, Wi-ley Online Library. pp. 105–116. 2714] Herholz, P., Haase, F., Alexa, M., 2017. Diffusion diagrams: Voronoi cells andcentroids from diffusion, in: Computer Graphics Forum, Wiley Online Library.pp. 163–175.[15] Hughes, T.J., Cottrell, J.A., Bazilevs, Y., 2005. Isogeometric analysis: Cad, finiteelements, nurbs, exact geometry and mesh refinement. Computer methods inapplied mechanics and engineering 194, 4135–4195.[16] Kaplan, C.S., Salesin, D.H., . Dihedral escherization, in: Proc. of Graphics Inter-face, pp. 255–262.[17] Kaplan, C.S., Salesin, D.H., 2000. Escherization, in: Proceedings of the27th annual conference on Computer graphics and interactive techniques, ACMPress/Addison-Wesley Publishing Co.. pp. 499–510.[18] Lipnikov, K., 2013. On shape-regularity of polyhedral meshes for solving pdes,in: Research Notes, 22nd International Meshing Roundtable.[19] Livesu, M., 2018. A heat flow relaxation scheme for n dimensional discrete hypersurfaces. Computers & Graphics 71, 124 – 131. doi: .[20] Mascotto, L., 2018. Ill-conditioning in the virtual element method: Stabilizationsand bases. Numerical Methods for Partial Differential Equations 34, 1258–1281.[21] Mu, L., Wang, X., Wang, Y., 2015. Shape regularity conditions for polyg-onal/polyhedral meshes, exemplified in a discontinuous galerkin discretization.Numerical Methods for Partial Differential Equations 31, 308–325.[22] Patané, G., 2016. Star-laplacian spectral kernels and distances for geometry pro-cessing and shape analysis, in: Computer Graphics Forum, Wiley Online Library.pp. 599–624.[23] Reberol, M., Lévy, B., 2018. Computing the distance between two finite elementsolutions defined on different 3d meshes on a gpu. SIAM Journal on ScientificComputing 40, C131–C155. 2824] Rjasanow, D.S., Weißer, S., 2013. Higher order BEM-based FEM on polygonalmeshes. SIAM J. Numer. Anal. 241, 103–115.[25] Schneider, T., Dumas, J., Gao, X., Botsch, M., Panozzo, D., Zorin, D., 2019.Poly-spline finite element method. ACM Transactions on Graphics .[26] Schneider, T., Hu, Y., Dumas, J., Gao, X., Panozzo, D., Zorin, D., 2018. Decou-pling simulation accuracy from mesh quality, in: SIGGRAPH Asia 2018 Techni-cal Papers, ACM. p. 280.[27] Shewchuk, J.R., 1996. Triangle: Engineering a 2d quality mesh generator anddelaunay triangulator, in: Applied computational geometry towards geometricengineering. Springer, pp. 203–222.[28] Shewchuk, J.R., 2002. What is a good linear finite element? interpolation, con-ditioning, anisotropy, and quality measures (preprint). University of California atBerkeley 73, 137.[29] Simonson, L., Suto, G., 2009. Geometry template library for stl-like 2d opera-tions. Colorado: GTL Boostcon .[30] Stimpson, C., Ernst, C., Knupp, P., Pébay, P., Thompson, D., 2007. The verdictlibrary reference manual. Sandia National Laboratories Technical Report 9.[31] Sukumar, N., Tabarraei, A., 2004. Conforming polygonal finite elements". Inter-national Journal for Numerical Methods in Engineering 61, 2045–2066.[32] The CGAL Project, 2018. CGAL User and Reference Manual. 4.13 ed.,CGAL Editorial Board. URL: https://doc.cgal.org/4.13/Manual/packages.htmlhttps://doc.cgal.org/4.13/Manual/packages.html