Unstructured Geometric Multigrid in Two and Three Dimensions on Complex and Graded Meshes
UUnstructured Geometric Multigrid in Two and ThreeDimensions on Complex and Graded Meshes
Peter R. Brune † Matthew G. Knepley ‡ L. Ridgway Scott † October 29, 2018
Abstract
The use of multigrid and related preconditioners with the finite element method is oftenlimited by the difficulty of applying the algorithm effectively to a problem, especially when thedomain has a complex shape or the mesh has adaptive refinement. We introduce a simplificationof a general topologically-motivated mesh coarsening algorithm for use in creating hierarchies ofmeshes for geometric unstructured multigrid methods. The connections between the guaranteesof this technique and the quality criteria necessary for multigrid methods for non-quasi-uniformproblems are noted. The implementation details, in particular those related to coarsening,remeshing, and interpolation, are discussed. Computational tests on pathological test casesfrom adaptive finite element methods show the performance of the technique.
Multigrid methods enable the asymptotically optimal approximate numerical solution to a widearray of partial differential equations. Their major advantage is that they have, for a wide class ofproblems, O ( n ) complexity for solving a linear system of n degrees of freedom (DoFs) to a specifiedaccuracy. However, there are many stumbling blocks to their full implementation in conjunctionwith adaptive unstructured finite element methods.Algebraic multigrid (AMG) methods and widely-available libraries implementing them have madegreat strides in allowing for the widespread use of multilevel methods in computational science [37].This is due to their relatively black-box nature and ability to be used for a wide variety of commonproblems without much expert tuning. Interesting problems often must use meshes that have beenadaptively refined in order to resolve numerical error or scale of interest without regard to multileveldiscretization. In order to develop a fast multilevel method for such a problem, either an algebraicor geometric coarsening procedure must be employed. Practitioners faced with such problems are † Department of Computer Science, 1100 S. 58th St., University of Chicago, Chicago, IL 60637([email protected], [email protected]). Partial support from NSF grant DMS-0920960. ‡ Computation Institute, 5735 S. Ellis Ave., University of Chicago, Chicago, IL 60637, ([email protected]),Supported by the Office of Advanced Scientific Computing Research, Office of Science, U.S. Department of Energy,under Contract DE-AC02-06CH11357. a r X i v : . [ c s . NA ] A p r RADED UNSTRUCTURED MG
A fairly typical scheme for constructing a series of meshes for geometric multigrid methods involvesuniform refinement from some coarsest starting mesh. However, the goals of refinement for thepurpose of resolving discretization or modelling error and refinement for the purposes of creatinga multigrid hierarchy can easily be at odds with each other if care is not taken. If a mesh hasbeen refined to satisfy some analytical or adaptively determined error control, then it may havevastly different length scales existing in different parts of the mesh. It is because of this that thecombination of refinement in order to resolve physics and refinement for the purpose of multigridis fraught with peril.If one refines for the purpose of creating a multigrid hierarchy from an initial coarse mesh withsome grading to handle the numerical properties of the problem, the refinement of the finest meshwill not reflect the error properly. Conversely if some stages of the refinement for the physics areused for the multigrid hierarchy, one may be unable to guarantee that the meshes would satisfy thequality and size constraints required by multigrid. Because of this it is very difficult to make thetwo concepts, adaptive refinement and geometric multigrid, go hand in hand as a single adaptivemeshing strategy that satisfies the needs of both methods in the general case.In typical applications in computational science, one is often given a mesh that has already beenadapted to the physics of some problem and asked to solve some some series of linear equationsrepeatedly on this static mesh. Two examples where this happens quite frequently are in optimiza-tion problems, and in the solution of nonlinear equations by methods requiring the equations tobe linearized and solved repeatedly. In this case, the only available geometric multigrid methodsare based upon coarsening, and huge advantages may be reaped from precomputation of a series ofcoarser spaces that effectively precondition the problem.The need for a mesh refined to resolve error can be demonstrated in some fairly straightforwardcases. The need for a priori grading around a reentrant corner to resolve pollution effects [5] isa well-studied classical phenomenon [3, 10, 22] in adaptive finite element computation. Multigrid
RADED UNSTRUCTURED MG θ r ≥ π we can define a factor µ ≥ πθ r Then, given some constants C a > C b < ∞ related to the maximum length scale in the mesh,a mesh that will optimally resolve the error induced by the reentrant corner, for h as the lengthscale of the cell containing any point a distance r away from the corner, will satisfy: C a r − µ ≤ h ≤ C b r − µ This problem will be considered further in Section 4.Figure 1: Quasi-uniform (right) and a priori (left) refined meshes: 500 vertices
When the creation of a coarsened hierarchy is considered in the geometric sense, we mean that theproblem is reduced to finding some series of meshes { M ...M n } where M is the initial mesh andand M n is the coarsest mesh. This hierarchy is constructed starting from M . The differentialoperator is either rediscretized on each mesh in order to create the series of coarse problems, orinterpolators between function spaces on the meshes are used to create approximate coarse versionsof the original fine linear system, a process known as Galerkin multigrid [36]. The experiments inthis paper are all done by rediscretization.One restriction that we will maintain throughout this paper is that the meshes are simplicial andnode-nested. This means that the vertices of mesh M i +1 are some subset of the vertices of mesh M i . This is an assumption used by many mesh refinement and coarsening methods, as it is areasonable and common assumption that the vertices of the initial fine mesh adequately capturethe domain’s geometry. This is in contrast to fully nested multigrid hierarchies, where each coarsecell is the union of several fine cells. A series of non-nested unstructured coarse meshes do notnecessarily need to be node-nested either. For instance, the mesh hierarchy may be created byrepeatedly remeshing some CAD or exact description of the domain at different resolutions.There is some justification for extending the classical multigrid convergence results to the case ofnon-nested meshes [8, 9, 24]. Results of optimal multigrid convergence have been extended to non- RADED UNSTRUCTURED MG
For any mesh M and any cell τ in M let ρ τ be τ ’s incircle diameter and h τ its longest edge, anddefine the aspect ratio as AR ( τ ) = h τ ρ τ . (1)The quality of approximation [34] and matrix condition number [4] in FEM calculations dependstrongly on AR ( τ ). For some constant C AR , we must require that the aspect ratio of any given cellin { M ...M n } satisfy AR ( τ ) ≤ C AR . (2)The other half of the criteria are the local comparability conditions. Assume that if we have twomeshes, M k and M k − , we define, for some cell τ in M k S kτ = { T : T ∈ M k − , T ∩ τ (cid:54) = ∅} (3)which defines the set of cells in a mesh M k − overlapping a single cell τ in the next coarser mesh M k . We can state the local comparability conditions assup τ ∈ M k {| S kτ |} ≤ C o for some C o ≥ τ ∈ M k { sup T ∈ S kτ { h τ h T }} ≤ C l for some C l ≥ | M | = M Then, for some C m ≥ | M k | > C m | M k +1 | . (5)The implication is that at each coarser mesh there must be a sufficient decrease in the number ofcells. A geometric decrease in the number of cells over the course of coarsening is necessary in orderto have an O ( n ) method. RADED UNSTRUCTURED MG For node-nested coarsening of an unstructured mesh, we separate the process of mesh coarseninginto two stages: coarse vertex selection and remeshing. The problem of coarse vertex selection on amesh M involves creating some subset V c of the vertices V of M that is, by some measure, coarse.Several techniques have been described for coarse vertex selection on unstructured simplicial meshes.The most widely used of these methods are Maximum Independent Set (MIS) algorithms [2,11,17].These choose a coarse subset of the vertices of the mesh such that no two vertices in the subset areconnected by a mesh edge. The resulting set of vertices is then remeshed.Given some mesh M , we define the graph G M = ( V, E ), which consists of the edges E and vertices V of M . The MIS algorithms may then be described as finding some set of coarse vertices V MISc such that ( v a , v b ) (cid:54)∈ E for all pairs of vertices ( v a , v b ) ∈ V MISc . This may be implemented by a simple greedy algorithm where at each stage a vertex is selectedfor inclusion in V MISc and its neighbors are removed from potential inclusion. There are a coupleof issues with this method. The first being that there is no way to determine the size of the meshresulting from coarsening. The other is that there are no real guarantees on the spatial distributionof vertices in the resulting point set. It has been shown that MIS methods may, for very simpleexample meshes, degrade the mesh quality quite rapidly [25]. Other methods for choosing the coarsevertex selection have been proposed to mitigate this shortcoming, often based upon some notion oflocal increase in length scale [30].
The coarse vertex selection method we will focus on for this was developed by Miller et al. [25] andis referred to as function-based coarsening. The method begins by defining some spacing function Sp ( v ) over the vertices v ∈ V of M . Sp ( v ) is some approximation of the local feature size. Inpractice, the approximation of the local feature size based upon the nearest-mesh-neighbor distanceat each vertex is a reasonable and efficiently computable spacing function.Define β as the multiplicative coarsening parameter. β can be considered the the minimum factorby which Sp ( v ) is increased at v by the coarsening process. Here we choose β > β where β = √ β = √ β = β + (cid:15) for small (cid:15) would reproduce therepeated structured coarsening of an isotropic, structured, shape-regular mesh where the lengthscale is increased by two in every direction. β may also be tuned, changing the size of the setof resulting coarse vertices in a problem-specific fashion to account for mesh and function-spaceproperties, such as polynomial order.We say that a coarse subset of the vertices of M , V F BCc satisfies the spacing condition if β ( Sp ( v i ) + Sp ( v j )) < dist ( v i , v j ) for all pairs ( v i , v j ) ∈ V F BCc . (6)After determining V F BCc , M c may be created by some remeshing process. A hierarchy of meshesmay then be created by reapplying the coarsening procedure to each newly-minted coarse meshwith constant β and some recalculated Sp in turn to create a yet coarser mesh. This may be done RADED UNSTRUCTURED MG h i ( x ) being the length-scaleof the cell overlapping the point x for all x ∈ Ω:1 I h i +1 ( x ) ≤ h i ( x ) ≤ I h i +1 ( x ) for some I > M n and some small constant b , | M n | ≤ b. This implies that the coarsening procedure will be able to create a constant-sized coarse mesh.Because of these conditions and their parallels with the conditions on the multigrid hierarchy, thismethod for coarsening is particularly appealing.Parts of the function-based coarsening idea have been incorporated into other methods for meshcoarsening by Ollivier-Gooch [28]. Some similarities between the method we propose here and thisprevious work are that both use traversal of the mesh and local remeshing in order to coarsenmeshes that exhibit complex features. This method constructs the conflict graph G C = ( V, E C )where E C = { ( v i , v j ) : β ( Sp ( v i ) + Sp ( v j )) > dist ( v i , v j ) } . This graph is then coarsened with an MIS approach as shown above. Note that in the limit ofextreme grading or large β the G C can grow to be of size O ( | V | ). Our method avoids this bylocalizing and simplifying the notion of the spacing function and choice of vertex comparisons. Wemodify function based coarsening in a way that reliably guarantees optimal complexity irregardlessof mesh non-quasi-uniformity as discussed in Section 1.1 without dependence on the mesh and theparameter β . We describe here a greedy algorithm for determining an approximation V GCAc to the set of coarsevertices V F BCc satisfying a weakened notion of the spacing condition based upon G M . Talmor[38] proposed using shortest-weighted-path distance determined by traversing along mesh edges toaccomplish this. Instead of doing using the shortest-weighted-path approach, we choose to usethe edge connectivity to progressively transform G M into some final G ˜ M c , which approximates theconnectivity of the coarse mesh, M c . It should be noted that, beyond the initial G M formed at thestart of the algorithm, G M does not necessarily satisfy the condition that its edges represent the RADED UNSTRUCTURED MG (a) G M with Sp ( v ) (b) G M with βSp ( v )(c) Subalgorithm on v (d) Subalgorithm on v Figure 2: This shows coarsening of a patch of G M . Sp is defined, multiplied by β , and thesubalgorithm is invoked for two vertices, v and v shown in Figs. 2(c) and 2(d) in dark grey withvertices violating the graph spacing condition shown in light gray. RADED UNSTRUCTURED MG M and V GCA as it is created in this section isdiscussed in Section 2.3.We begin this process by modifying condition 6 to be β ( Sp ( v ) + Sp ( v )) < dist ( v , v ) for all ( v , v ) ∈ E of G M . (7)This restricts the spacing condition to take into account only distance between vertices connectedby graph edges rather than all pairs of vertices. This restriction makes sense considering thecalculation of the Sp ( v ) as nearest-neighbor distance is based upon adjacency of the mesh. This isan important simplification for complex domains, where the mesh may represent complex featuresof the domain such as holes, surfaces, and interfaces with different local characteristics on each sidethat may be hard to encode when the vertices are considered as a point cloud. In our model, theedges of the mesh are seen as encoding the topological connectivity of the domain. Spacing onlyapplies to topologically connected areas of the mesh as encoded in G M .Define F ( v ) = unknown , included , or excluded for all v ∈ V. The algorithm starts with F ( v ) = unknown for all v ∈ V . Vertices that must be included for somereasons such as the boundary issues described in Section 2.4, are set to be included automaticallyand visited first. The remaining v ∈ V are then visited in some arbitrary order. If F ( v ) = excluded ,the iteration skips the vertex. If F ( v ) = unknown , F ( v ) is changed to included .When F ( v ) is set to included , a subalgorithm is invoked that, given G M and some particular vertex v , transforms G M to some intermediate G vM . This subalgorithm corresponds to coarsening the areaaround v until the spacing function is satisfied for all N G vM ( v ). Each w ∈ N G M ( v ) are tested to seeif the edge ( v, w ) violates Eq. 7. In the case that the condition is violated and F ( w ) = unknown , w is removed from G M and F ( w ) is changed to excluded . A removed w ’s neighbors in G M areadded to N G M ( v ) by adding edges in G M between v and all u ∈ N G M ( w ) if u (cid:54)∈ N G M ( v ) already.This may also be considered as an edge contraction from w to v .The outcome of this subalgorithm, G vM , has that all w ∈ N G vM ( v ) such that F ( w ) = unknown obeyEq. 7. There is the possibility that for some v ∈ N G vM ( v ), F ( v ) = included and Eq. 7 is notsatisfied. This may arise due to the necessary inclusion of boundary vertices due to conditionsdescribed in Section 2.4. After the subalgorithm terminates, one is left with G vM . The algorithmthen visits some v (Figure 2(d)), which, if included , will create some G v M by modification of G vM .Once every vertex is visited, the whole algorithm terminates. At this point we may define V GCAc = { v : v ∈ V, F ( v ) = included } . Despite the fact that we are considering coarsening and remeshing as a preprocessing step, we stillhave the goal that the algorithm should be O ( n ) with n as the number of DoFs in the discretizedsystem. For piecewise linear finite elements, the number of DoFs in the system is exactly thenumber of vertices, | V | , in M . For all reasonably well-behaved meshes, we can additionally assumethat | E | = O ( | V | ) and | M | = O ( | V | ). This implies that n = O ( | V | ). The complexity of a singleinvocation of the subalgorithm may be O ( | V | + | E | ) if some large fraction of the mesh vertices arecontracted to some v in one step. However, the aggregate number of operations taken to reach G ˜ M c must be the order of the size of the original graph G M . This aggregate complexity bound isindependent of the order in which v are visited. While here we keep the assumed order in which RADED UNSTRUCTURED MG
Theorem 2.1.
Given a graph derived from a mesh G M , the graph coarsening algorithm will create V GCAc in O ( | V | + | E | ) time.Proof. The fact that the complexity of the algorithm depends on | V | at least linearly is apparentfrom the fact that every vertex is visited once. In order to show that the entire algorithm is O ( | V | + | E | ) we must establish that each edge e is visited at most a constant number of timesindependent of the size of the graph.As vertex v is visited and F ( v ) is set to included , the subalgorithm inspects each e = ( v, n ) for n ∈ N G M ( v ) to see if they satisfy the spacing condition with v . e is either deleted from G M if n and v do not satisfy the spacing condition, or left in place if the spacing condition is satisfied or if F ( n ) = included .Edges that are deleted are not ever considered again. Therefore, we must focus on edges that survivecomparison. Suppose that an edge e = ( v , v ) is considered a second time by the subalgorithm atthe visit to vertex v . As this is the second visit to consider e , F ( v ) must be included as in thefirst consideration of e necessarily happened during the visit to v . As both endpoints of e havenow been visited, there is no way e may be considered a third time. As each vertex is visited onceand the distance across each edge in G M is considered no more than twice, the algorithm runs in O ( | V | + | E | ) time. In the node-nested case, the problem of remeshing may be treated as a procedure taking a mesh M and some coarse subset of M ’s vertices V c and returning a new mesh M c that approximates thedomain spanned by M but only includes vertices V c . Traditional constrained remeshing methodsusing standard meshing technology require specialized information about the domain, includingthe location and topology of holes in the domain, that is not typically provided with a pre-gradedfine mesh. This inhibits automation since a great deal of expert input is necessary to preserve theboundary during coarsening. Instead, we choose to locally remove the set of vertices V \ V c from M one at a time to create M c .In the 2D case, we remove a vertices from M in-place by a simple series of geometric predicatesand edge flips. This procedure is completely localized to the neighborhood around the vertex beingremoved. This is done in a way that retains the Delaunay condition [27] if the initial mesh wasDelaunay. Similar but much more complex techniques exist in 3D [21], but are difficult to makecomputationally robust enough for the requirements of multigrid methods, where the creation of asingle bad cell may severely impact the strength of the resulting multigrid preconditioner. In thiswork, we relax the Delaunay condition and instead focus on mesh quality of the sort required bythe multigrid method as stated in Section 1.2.An extremely simple method of removing a vertex v from the mesh M is by quality-conserving edgecontraction. Similar approaches have been proposed for this problem previously [28]. This may bestated as choosing n ∗ ∈ N ( v ) such that for the set of cells C n that would be created during the RADED UNSTRUCTURED MG n ∗ = arg min n ∈ N ( v ) (max c ∈ C n AR ( c ))subject to the constraint max c ∈ C n AR ( c ) < C AR . Each resulting c must also be properly oriented in space.The computational work done to remove a single vertex is O ( | C v || N ( v ) | ) as each potential configu-ration C n has that, when compared to the set of cells adjacent to v , C v , | C n | < | C v | because at leasttwo cells adjacent to both n and v must be removed by each contraction. For each n ∈ N ( v ) onemust check the aspect ratio and orientation of each cell in the associated C n . Assuming boundedmesh degree, the removal of O ( | V | ) vertices would take O ( | V | ) time.Note that we are directly enforcing the maximum aspect ratio that a given mesh modification maycreate, which makes this akin to best effort remeshing. If it is impossible to satisfy the qualityconditions required of the mesh while removing a vertex, the vertex is left in place, at least forthe time being. Such vertices are revisited for removal if their link is modified by the removalof one of their neighbors, a process that may only happen O ( N ( v )) times, bounding the numberof attempts that may be made to remove a given vertex. Due to this, we cannot guarantee thatthis 3D remeshing finds a tetrahedralization including just the coarse vertices. However, we notethat our experience with this method is that the number of vertices that must be left in place issignificantly less than the typical number of Steiner points added through direct Delaunay remeshingusing Tetgen [35], and that this method produces tetrahedralizations entirely suitable for use withmultigrid methods. The results of a study of the quality of meshes created by this method is shownin Section 4.1.
Preservation of the general shape of the domain is important to the performance of the multigridmethod as the coarse problem should be an approximation of the fine problem. In the worstcase a coarser mesh could become topologically different from the finer mesh leading to very slowconvergence or even divergence of the iterative solution. If this is to be an automatic procedure,then some basic criteria for the shapes of the sequence of meshes must be enforced. Therefore, thevertex selection and remeshing algorithms must be slightly modified in order to take into accountthe mesh boundaries.First, we must define the features, in 2D and 3D, which require preservation. We choose thesefeatures by use of appropriate-dimensional notions of the curvature of the boundary. Techniqueslike these are widely used in computational geometry and computer graphics [15, 31]. In 2D, thefeatures we choose to explicitly preserve are the corners of the mesh. We consider any boundaryvertex with an angle differing from a straight line more than C K to be a corner. For the sake ofthis work we assume C K = π .In 3D, the discrete curvature at a boundary vertex is computed as the difference between 2 π andthe sum of the boundary facet angles tangent to that vertex. Vertices where the absolute valueof the discrete curvature is greater than C K are forced to be included. High-dihedral-angle edges, RADED UNSTRUCTURED MG (a) interior coarsening (b) boundary coarsening (c) final graph Figure 3: Interior and boundary coarsening stages. Dark gray vertices are included automatically;Light gray circles violate the spacing condition.corresponding to ridges on the surface of the domain, must also be preserved. We consider anyboundary edge with the absolute value of the dihedral angle greater than C K to be a boundaryridge.Our approach to protecting the boundary during the coarse vertex selection procedure is to separatevertex selection into two stages, one for the interior, and one for the boundary. In the interior stage,the boundary vertices are marked as included automatically, and therefore any interior verticesviolating the spacing condition with respect to a boundary vertex are removed from G M (Figure3(a)). All boundary vertices now respect the spacing condition with respect to all remaining interiorvertices, making the second stage, boundary vertex selection, entirely independent of the previouscoarsening of the interior. The boundary vertex selection procedure then operates independently ofthis, producing the completely coarsened graph G ˜ M c (Figure 3(c)). In 3D this process is repeatedonce again. During the boundary coarsening procedure, vertices lying on edges that have beenidentified as ridges are automatically included . The ridge vertices are then coarsened in turn.Corner vertices are automatically included during all stages of vertex selection, as shown in Figure3(b). For remeshing, the only modification necessary for boundary preservation is that the removalof a boundary or ridge vertex may only be done by contraction along an edge lying on the boundaryor ridge. This simple modification works especially well for meshes that have a lot of flat planes ontheir surfaces with sharp ridges separating them. Without this procedure, ridges may appear tornor flipped in the coarser meshes, changing the shape of the domain significantly. As the hierarchycoarse meshes are created, new edges or vertices may become marked as corners or as belonging toridges. We have shown that local traversal is a powerful tool for the construction of node-nested coarsemeshes. Traversal of the node-nested meshes may also be used to efficiently construct the interpo-lation operators I l for l ∈ [0 , n − I l is defined as the interpolation operator between the finiteelement function spaces V hc and V hf , which are defined over the coarse mesh M c and fine mesh M f respectively for the adjacent f = l and c = l + 1. We also assume that the domain Ω spanned byall M is connected. For the purposes of this paper, we assume that for each DoF i , associated withbasis function φ fi in V fh in the system has some nodal point x i ∈ R d and that a function f may be RADED UNSTRUCTURED MG V fh by constructing the coefficients f i = f ( x i ) · ψ fi ( x i ) . Given this assumption, we may construct the prolongation operator from functions in V ch to V fh byassociating each i with a cell τ ck in M c and constructing, for ψ fi and basis functions ψ cj supportedon τ ck : I ij = ψ cj ( x i ) · ψ fi ( x i ) . One can, for the sake of simplicity, associate every ψ fj with a single cell τ fi . The choice of cell isarbitrary between the cells the unknown is supported on. This makes the problem equivalent tofinding, for some set of x ∈ Ω, T c ( x ), defined as the cell in M c that x is located in. We also define˜ T c ( x ), which is some cell in M c that is, in some sense, nearby x . In addition to the points x i associated with all ψ f , we have the midpoint x τ fi of every cell τ fi ∈ M f .We define a two-level nested-loop traversal algorithm. The outer loop locates T c ( x τ f ) for each τ f ∈ M f , and the inner loop determines T c ( x i ) for all x i associated with ψ fi supported on τ fi .Once T c ( x i ) is resolved for all i , one may construct the nodal interpolation operator.The outer loop consists of breadth-first searches (BFSes) on the graph of the neighbor relations ofthe cells of M c in order to determine T c ( x τ k ) for each τ k ∈ M f . This is implemented by a simpleFIFO queue in which the neighbors of a given cell are pushed on the back of the queue when it isvisited. Enqueued and visited cells are marked as such and not enqueued again. We say that twocells τ c and τ c are neighbors if they share a vertex, edge, or facet. The BFS to locate τ k starts withthe cell ˜ T c ( x τ k ).As T c ( x τ fk ) is established for each cell τ fk , the inner loop is invoked. This loop consists of a BFS foreach ψ fi associated with τ fk and determines T c ( x i ) by BFS starting from ˜ T c ( x i ) = T c ( x τ fk ).The last ingredient in the algorithm is how to determine ˜ T c ( x τ fk ). One simple way of doing this issetting ˜ T c ( x τ fk ) = T c ( x τ fm ) for any cell τ fm that is neighboring τ fk . This notion of locality may beexploited by setting the values of ˜ T c ( τ n ) for all neighbors τ n of a cell τ upon determining T c ( τ ).When T c ( τ f ) for any τ f ∈ M f is determined, the connectivity of M c and M f may be effectivelytraversed in tandem in order to extend the search over the whole meshed domain.This process may be both started and made more efficient by exploiting the node-nested propertiesof the meshes. We have that meshes M f and M c will share some set of vertices V ( f,c ) due to ournode-nested coarsening procedure. Define τ fv to be a cell in mesh M f that is adjacent to some v ∈ V ( f,c ) and τ cv to be an arbitrary cell in M c adjacent to that same v . Then, one may initialize˜ T c ( x τ fv ) to be τ cv .An obervation about the complexity of this algorithm may be established assuming that the condi-tions in Section 1.2 are satisfied by the hierarchy of meshes. It should be obvious that the complexityis going to be bounded by the total number of BFSes done during the course of the algorithm mul-tiplied by the worst-case complexity of an invocation of the BFS. It’s easy to see that for | M l | cellsin the fine mesh and n unknowns, the geometric search must be done | M l | + n times.In order to bound the complexity of a given BFS to a constant, we must show how many stepsof traversal one must search outwards in order to locate some T c ( τ f ) given ˜ T c ( τ f ). This may be RADED UNSTRUCTURED MG dist ( x τ f , x τ f ) is going to be lessthan some maximum search radius r τ f = h τ f + h τ f given that they are adjacent. We also may puta lower limit on the minimum length scale of cells in M c that overlap τ f and τ f by the constant inEq. 4b and by the aspect ratio limit in Eq. 2 as h c min = C h τf C AR .This gives us that the number of cells that may fit between ˜ T c ( τ f ) and T c ( τ f ) is at most (cid:108) r fτ h c min (cid:109) ,which is independent of overall mesh size. Note that the distance between the center of a cell τ and some nodal points x i of ψ i supported on τ is always less than or equal to h τ . This extends thisconstant-size traversal bound to the inner loop also, making the entire algorithm run in O ( | M f | + n )time. In practice on isotropic meshes ˜ T ( τ f ) and T ( τ f ) are almost always in the same cell or onecell over, meaning that the topological search is an efficient method for building the interpolationoperators.An extension of the algorithm that has proven necessary on complex meshes is to allow for inter-polation to x i not located within a cell of the mesh. For example, the Pacman mesh (Fig. 4), whencoarsened, will have vertices that lie outside the coarser mesh on the round section. In order to dothis, the outer loop BFSes are replaced by a more complex arrangement where the BFSes searchfor the nearest cell rather than an exact location. This nearest cell then becomes T ( x τ ). The innerloop is replaced by a similar procedure, which then projects any x i that could not be exactly locatedto some ˜ x i on the surface of the nearest cell and modifies the basis function projection to use φ cj (˜ x i )instead of φ cj ( x i ).This procedure is easily extended to discretizations consisting of non-nodal DoFs by locating setsof quadrature points associated with the fine mesh cells instead of DoF nodes. This also allows forthe construction of pseudointerpolants satisfying smoothness assumptions, such as the Scott-Zhangtype interpolants [32]. We test the multigrid method using standard isotropic non-quasi-uniform examples from the liter-ature. The domains used are the Pacman in 2D and the Fichera corner in 3D. The problems andboundary conditions computed are designed to induce a corner singularity that makes refinementsuch as that described in Section 1.1 necessary.The standard model problem for the Pacman domain is [3] designed to have a pollution effect atthe corner singularity. The Pacman domain is a unit circle with a tenth removed. The reentrantangle is therefore θ r = π radians. The associated differential equation is − ∆ u = 0 in Ω ,u ( r, θ ) = r sin (cid:18) θ (cid:19) on ∂ Ω . For the Fichera corner model problem, we separate the boundary of the domain into the reentrantsurface ∂ Ω i and the outer surface ∂ Ω o . The Fichera corner mesh consists of a twice-unit cubecentered at the origin with a single octant removed. The reentrant angle is θ r = π radians RADED UNSTRUCTURED MG − ∆ u = 0 in Ω ,∂u∂n = g on ∂ Ω o ,u = 0 on ∂ Ω i . Where g ( x, y, z ) = g xy + g yz + g zy ,g x i x j ( r, θ ) = r cos( θ ) . where r and θ consist of the radius and angle when restricted to the x i , x j plane. Unfortunately,this problem does not have an analytical exact solution. Our goal is to show that we can effectivelycoarsen the a priori grading required to resolve the problem, so we will not be concerned with thesolution accuracy as a measure of performance. In order to motivate the use of the method with the multigrid solver, one must look if our imple-mentation of the method actually lives up to the mesh quality metrics stated as requirements forguaranteed multigrid performance. This study was done using large initial Pacman and Fichera cor-ner meshes coarsened using the same algorithm used for the multigrid studies in Section 4.2.The measurements of mesh quality we have taken correspond to the bounds on the meshes in themultigrid requirements. These are the number of vertices and cells in the mesh, the worst aspectratio in each mesh, the maximum number of cells in a mesh overlapping any given cell in the nextcoarser mesh, and the maximum change in local length-scale between a mesh and the next coarsermesh at any point in the domain. The meshes are created by repeated application of the coarseningalgorithm with β = 1 . β = 1 . RADED UNSTRUCTURED MG β = 1 . AR ( τ )) sup τ ∈ M k {| S kτ |} max x ∈ Ω (cid:18) h k − ( x ) h k ( x ) (cid:19) β = 1 . C AR = 60level cells vertices max( AR ( τ )) sup τ ∈ M k {| S kτ |} max x ∈ Ω (cid:18) h k − ( x ) h k ( x ) (cid:19) C AR for most of the levels. Further work on our incrediblysimple remeshing scheme should be able to improve this. However, we do not see successive degra-dation of the quality of the series of meshes as the coarsening progresses. We can assume that thequality constraints are reasonably met in both 2D and 3D. RADED UNSTRUCTURED MG L Error is computed based upon difference from the exact solution for the Pacman testproblem. (a) Pacman Problem Performance
DoFs levels MG ILU L Error762 3 6 76 9.897e-61534 4 6 129 3.949e-63874 5 6 215 9.807e-79348 5 7 427 2.796e-723851 6 7 750 1.032e-735660 7 7 1127 8.858e-878906 7 7 1230 3.802e-8139157 8 9 3262 1.535e-8175602 8 9 3508 3.440e-9 (b) Fichera Problem Performance
DoFs levels MG ILU2018 2 7 463670 2 8 609817 3 8 11227813 4 9 20758107 5 9 385107858 6 9 543153516 6 9 440206497 7 9 616274867 8 9 652
The performance of multigrid based upon the mesh series created by this algorithm using the twotest examples was carried out using the DOLFIN [23] finite element software modified to use thePCMG multigrid preconditioner framework from PETSc [6]. The operators were discretized usingpiecewise linear triangular and tetrahedral finite elements. The resulting linear systems were thensolved to a relative tolerance of 10 − . The solvers used were ILU-preconditioned GMRES for thestandard iterative case, shown in the ILU columns as a control. In the multigrid case we choseto use GMRES preconditioned with V-cycle multigrid. Three pre and post-smooths using ILU asthe smoother were performed on all but the coarsest level, for which a direct LU solve was used.For this to make sense, the coarsest mesh must have a small, nearly constant number of vertices;a condition easily satisfied by the mesh creation routine. We coarsen until the coarsest mesh isaround 200 vertices in 2D or 300 vertices in 3D.These experiments show that the convergence of the standard iterative methods becomes more andmore arduous as the meshes become more and more singularly refined. The singularity in 2D ismuch greater, due to the sharper reentrant angle, than it is in 3D, so the more severe difference inperformance between multigrid and the control method is to be expected. We see that the numberof multigrid cycles levels out to a constant for both 2D (Table 3(a)) and 3D (Table 3(b)). We alsosee that a steadily increasing number of multigrid levels are automatically generated and that themethod continues to behave as expected as the coarsening procedure is continually applied. We have proposed and implemented a method for mesh coarsening that is well suited to the construc-tion of geometric unstructured multigrid on 2D and 3D meshes. The technique used to constructthe mesh hierarchy is novel in that it combines the advantages of graph-based mesh coarsening
RADED UNSTRUCTURED MG
References [1]
Mark Adams , Evaluation of three unstructured multigrid methods on 3D finite element prob-lems in solid mechanics , Int. J. Numer. Meth. Engng, 55 (2000), p. 519534.[2]
Mark Adams and James W. Demmel , Parallel multigrid solver for 3D unstructured finiteelement problems , in Supercomputing ’99: Proceedings of the 1999 ACM/IEEE conference onSupercomputing (CDROM), New York, NY, USA, 1999, ACM, pp. 27+.[3]
Thomas Apel and Frank Milde , Comparison of several mesh refinement strategies nearedges , Communications In Numerical Methods In Engineering, 12 (1996), pp. 373–381.[4]
I. Babuˇska and A. K. Aziz , On the angle condition in the finite element method , SIAMJournal on Numerical Analysis, 13 (1976), pp. 214–226.
RADED UNSTRUCTURED MG
Ivo Babuˇska, R. B. Kellogg, and J. Pitk¨aranta , Direct and inverse error estimates forfinite elements with mesh refinements , Numerische Mathematik, 33 (1979), pp. 447–471–471.[6]
Satish Balay, Kris Buschelman, William D. Gropp, Dinesh Kaushik, Matthew G.Knepley, Lois Curfman McInnes, Barry F. Smith, and Hong Zhang , PETSc Webpage
Randolph Bank , Some variants of the Bank–Holst parallel adaptive meshing paradigm , Com-puting and Visualization in Science, 9 (2006), pp. 133–144.[8]
M. L. Bittencourt, C. C. Douglas, and R. A. Feij´oo , Non-nested and non-structuredmultigrid methods applied to elastic problems , Part II: The three-dimensional case. Submittedfor publication, (1998).[9] ,
Adaptive non-nested multigrid methods , Engineering Computations, 19 (2002), pp. 158–176.[10]
H. Blum and R. Rannacher , Extrapolation techniques for reducing the pollution effect ofreentrant corners in the finite element method , Numerische Mathematik, 52 (1988), pp. 539–564.[11]
T. Chan, B. Smith, and J. Zou , Multigrid and domain decomposition methods for unstruc-tured meshes , 1994.[12]
L. Chen, M. Holst, J. Xu, and Y. Zhu , Local Multilevel Preconditioners for EllipticEquations with Jump Coefficients on Bisection Grids , ArXiv e-prints, (2010).[13]
L. Chen and C. Zhang , A coarsening algorithm on adaptive grids by newest vertex bisectionand its applications , Journal of Computational Mathematics, 28 (2010), pp. 767–789.[14]
Edmond Chow , An unstructured multigrid method based on geometric smoothness , AmericanJournal of Mathematics, 65 (2001), pp. 197–215.[15]
N. Dyn, K. Hormann, S.-J. Kim, and D. Levin , Optimizing 3D Triangulations UsingDiscrete Curvature Analysis , Vanderbilt University, Nashville, TN, USA, 2001, pp. 135–146.[16]
D. Feuchter, I. Heppner, S. Sauter, and G. Wittum , Bridging the gap between ge-ometric and algebraic multigrid methods , Computing and visualization in science, 6 (2003),pp. 1–13.[17]
Herv´e Guillard , Node-nested multi-grid method with Delaunay coarsening , Research ReportRR-1898, INRIA, 1993.[18]
Jim E. Jones and Panayot S. Vassilevski , AMGE Based on Element Agglomeration ,SIAM Journal on Scientific Computing, 23 (2001), pp. 109–133.[19]
M. Jung , Some multilevel methods on graded meshes , Journal of Computational and AppliedMathematics, 138 (2002), pp. 151–171.[20]
Michael K¨oster and Stefan Turek , The influence of higher order FEM discretisations onmultigrid convergence , Computational Methods in Applied Mechanics, 6 (2005), pp. 221–232.[21]
Hugo Ledoux, Christopher M. Gold, and George Baciu , Flipping to robustly delete avertex in a delaunay tetrahedralization. , in ICCSA (1), 2005, pp. 737–747.
RADED UNSTRUCTURED MG
Xiaohai Liao and Ricardo H. Nochetto , Local a posteriori error estimates and adaptivecontrol of pollution effects , Numerical Methods for Partial Differential Equations, 19 (2003),pp. 421–442.[23]
A. Logg and G. N. Wells , DOLFIN: Automated finite element computing , ACM Transac-tions on Mathematical Software, 37 (2010), pp. 20:1–20:28.[24]
R. L¨ohner and K. Morgan , An unstructured multigrid method for elliptic problems , Inter-national Journal for Numerical Methods in Engineering, 24 (1987), pp. 101–115.[25]
Gary L. Miller, Dafna Talmor, and Shang-Hua Teng , Optimal coarsening of unstruc-tured meshes , J. Algorithms, 31 (1999), pp. 29–65.[26]
E. Morano, D. J. Mavriplis, and V. Venkatakrishnan , Coarsening strategies for un-structured multigrid techniques with application to anisotropic problems , SIAM J. Sci. Comput.,20 (1998), pp. 393–415.[27]
M. Mostafavi , Delete and insert operations in Voronoi/Delaunay methods and applications ,Computers & Geosciences, 29 (2003), pp. 523–530.[28]
Carl Ollivier-Gooch , Coarsening unstructured meshes by edge contraction , Int. J. Numer.Meth. Engng., 57 (2003), pp. 391–414.[29]
W. Rachowicz, D. Pardo, and L. Demkowicz , Fully automatic hp-adaptivity in threedimensions , Computer Methods in Applied Mechanics and Engineering, 195 (2006), pp. 4816–4842.[30]
B. Rey, K. Mocellin, and L. Fourment , A node-nested Galerkin multigrid method formetal forging simulation , Computing and Visualization in Science, 11 (2008), pp. 17–25.[31]
Szymon Rusinkiewicz , Estimating curvatures and their derivatives on triangle meshes ,3D Data Processing Visualization and Transmission, International Symposium on, 0 (2004),pp. 486–493.[32]
L. Ridgway Scott and Shangyou Zhang , Finite element interpolation of nonsmooth func-tions satisfying boundary conditions , Mathematics of Computation, 54 (1990), pp. 483–493.[33] ,
Higher-dimensional nonnested multigrid methods , Mathematics of Computation, 58(1992), pp. 457 – 466.[34]
Jonathan Richard Shewchuk , What is a good linear element? - interpolation, conditioning,and quality measures , in In 11th International Meshing Roundtable, 2002, pp. 115–126.[35]
Hang Si , TetGen website , 2007. http://tetgen.berlios.de/index.html.[36]
Barry Smith, Petter Bjørstad, and William Gropp , Domain Decomposition: ParallelMultilevel Methods for Elliptic Partial Differential Equations , Cambridge University Press,March 2004.[37]
K. St¨uben , A review of algebraic multigrid , J. Comput. Appl. Math., 128 (2001), pp. 281–309.[38]
Dafna Talmor , Well-Spaced Points for Numerical Methods , PhD thesis, Carnegie MellonUniversity, Pittsburgh, Pennsylvania, 1997.
RADED UNSTRUCTURED MG
Shangyou Zhang , Optimal-order nonnested multigrid methods for solving finite element equa-tions I: On quasi-uniform meshes , Mathematics of Computation, 55 (1990), pp. 23–36.[40] ,
Optimal-order nonnested multigrid methods for solving finite element equations II: Onnon-quasi-uniform meshes , Mathematics of Computation, 55 (1990), pp. 439–450.[41] ,