A thread-parallel algorithm for anisotropic mesh adaptation

Solution discretisation errors are closely related to the size and the shape of the elements. However, in general meshes are generated using a priori information about the problem under consideration when the solution error estimates are not yet available. This may be problematic because:

• 1.

  Solution errors may be unacceptably high.

• 2.

  Parts of the solution may be over-resolved, thereby incurring unnecessary computational expense.

A solution to this is to compute appropriate local error estimates, and use this information to compute a field on the mesh which specifies the local mesh resolution requirement. In the most general case this is a metric tensor field so that the resolution requirements can be specified anisotropically; for a review of the procedure see [frey2005anisotropic]. Size gradation control can be applied to this metric tensor field to ensure that there are not abrupt changes in element size [alauzet2010size].

As discretisation errors are dependent upon element shape as well as size, a number of different measures of element quality have been proposed, e.g. [buscaglia1997anisotropic, vasilevskii1999adaptive, agouzal1999adaptive, tam2000anisotropic, pain_tetrahedral_2001]. For the work described here, we use the element quality measure for triangles proposed by Vasilevskii et al. [vasilevskii1999adaptive]:

where $|\triangle|_{M}$ is the area of element $\triangle$ and $|\partial\triangle|_{M}$ is its perimeter as measured with respect to the metric tensor $M$ as evaluated at the element’s centre. The first factor ($I$) is used to control the shape of element $\triangle$. For an equilateral triangle with sides of length $l$ in metric space, $|\triangle|=l^{2}\sqrt{3}/4$ and $|\partial\triangle|=3l$; and so $I=1$. For non-equilateral triangles, $I<1$. The second factor ($II$) controls the size of element $\triangle$. The function $F$ is smooth and defined as:

which has a single maximum of unity with $x=1$ and decreases smoothly away from this with $F(0)=F(\infty)=0$. Therefore, $II=1$ when the sum of the lengths of the edges of $\triangle$ is equal to $3$, e.g. an equilateral triangle with sides of unit length, and $II<1$ otherwise. Hence, taken together, the two factors in (1) yield a maximum value of unity for an equilateral triangle with edges of unit length, and decreases smoothly to zero as the element becomes less ideal.

The mesh is adapted through a series of local operations: edge collapse (§2.4.1); edge refinement (§2.4.2); element-edge swaps (§2.4.3); and vertex smoothing (§2.4.4). While the first two of these operations control the local resolution, the latter two operations are used to improve the element quality.

Algorithm 1 gives a high level view of the anisotropic mesh adaptation procedure as described by [li20053d]. The inputs are $\mathcal{M}$, the piecewise linear mesh from the modelling software, and $\mathcal{S}$, the node-wise metric tensor field which specifies anisotropically the local mesh resolution requirements. Coarsening is initially applied to reduce the subsequent computational and communication overheads. The second stage involves the iterative application of refinement, coarsening and mesh swapping to optimise the resolution and the quality of the mesh. The algorithm terminates once the mesh optimisation algorithm converges or after a maximum number of iterations has been reached. Finally, mesh quality is fine-tuned using some vertex smoothing algorithm (e.g. quality-constrained Laplacian smoothing [freitag1996comparisonw], optimisation-based smoothing [freitag1995cient]), which aims primarily at improving worst-element quality.

Before presenting the adaptive algorithms, it is necessary to give a brief description of the data structures used to store mesh-related information. Following that, we present a set of techniques which help us avoid hazards and data races and guarantee fast and safe concurrent read/write access to mesh data.

Every edge can be processed and refined without being affected by what happens to adjacent edges. Being free from structural hazards, the only issue we are concerned with is thread safety when updating mesh data structures. Refining an edge involves the addition of a new vertex to the mesh. This means that new coordinates and metric tensor values have to be appended to coords and metric and adjacency information in NNList has to be updated. The subsequent element split leads to the removal of parent elements from ENList and the addition of new ones, which, in turn, means that NEList has to be updated as well. Appending new coordinates to coords, metric tensors to metric and elements to ENList is done using the thread worklist strategy described in Section 3.1.4, while updates to NNList and NEList can be handled efficiently using the deferred operations mechanism.

The two stages, namely edge refinement and element refinement, of our threaded implementation are described in Algorithm 5 and Algorithm 6, respectively. The procedure begins with the traversal of all mesh edges. Edges are accessed using NNList, i.e. for each mesh vertex $V_{i}$ the algorithm visits $V_{i}$’s neighbours. This means that edge refinement is a directed operation, as edge $\overline{V_{i}V_{j}}$ is considered to be different from edge $\overline{V_{j}V_{i}}$. Processing the same physical edge twice is avoided by imposing the restriction that we only consider edges for which $V_{i}$’s ID is less than $V_{j}$’s ID. If an edge is found to be longer than desired, then it is split in the middle (in metric space) and a new vertex $V_{n}$ is created. $V_{n}$ is associated with a pair of coordinates and a metric tensor. It also needs an ID. At this stage, $V_{n}$’s ID cannot be determined. Once an OpenMP thread exits the edge refinement phase, it can proceed (without synchronisation with the other threads) to fix vertex IDs and append the new data it created to the mesh. The thread captures the number of mesh vertices $index=NNodes$ and increments it atomically by the number of new vertices it created. After capturing the index, the thread can assign IDs to the vertices it created and also copy the new coordinates and metric tensors into coords and metric, respectively.

Before proceeding to element refinement, all split edges are accumulated into a global worklist. For each split edge $\overline{V_{i}V_{j}}$, the original vertices $V_{i}$ and $V_{j}$ have to be connected to the newly created vertex $V_{n}$. Updating NNList for these vertices cannot be deferred. Most edges are shared between two elements, so if the update was deferred until the corresponding element were processed, we would run the risk of committing these updates twice, once for each element sharing the edge. Updates can be committed immediately, as there are no race conditions when accessing NNList at this point. Besides, for each split edge we find the (usually two) elements sharing it. For each element, we record that this edge has been split. Doing so makes element refinement much easier, because as soon as we visit an element we will know immediately how many and which of its edges have been split. An array of length NElements stores this type of information.

During mesh refinement, elements are visited in parallel and refined independently. It should be noted that all updates to NNList and NEList are deferred operations. After finishing the loop, each thread uses the worklist method to append the new elements it created to ENList. Once again, no thread synchronisation is needed.

This parallel refinement algorithm has the advantage of not requiring any mesh colouring and having low synchronisation overhead as compared with Freitag’s task graph approach.

Because any decision on whether to collapse an edge is strongly dependent upon what other edges are collapsing in the immediate neighbourhood of elements, an operation task graph for coarsening has to be constructed. Edge collapse is based on the removal of vertices, i.e. the elemental operation for edge collapse is the removal of a vertex. Therefore, the operation task graph $\mathcal{G}$ is the mesh itself.

Figure 6 demonstrates what needs to be taken into account in order to perform parallel coarsening safely. It is clear that adjacent vertices cannot collapse concurrently, so a distance-1 colouring of the mesh is sufficient in order to avoid structural hazards. This colouring also enforces processing of vertices topologically at least every other one which prevents skewed elements forming during significant coarsening[de1999parallel, li20053d].

An additional consideration is that vertices which are two edges away from each other share some common vertex $V_{common}$. Removing both vertices at once means that $V_{common}$’s adjacency list will have to be modified concurrently by two different threads, leading to data races. These races can be avoided using the deferred operations mechanism.

Algorithm 7 illustrates a thread parallel version of mesh edge collapse. Coarsening is divided into two phases: the first sweep through the mesh identifies what edges are to be removed, see Algorithm 8; and the second phase actually applies the coarsening operation, see Algorithm 9. Function coarsen_identify($V_{i}$) takes as argument the ID of a vertex $V_{i}$, decides whether any of the adjacent edges can collapse and returns the ID of the target vertex $V_{t}$ onto which $V_{i}$ should collapse (or a negative value if no adjacent edge can be removed). coarsen_kernel($V_{i}$) performs the actual collapse, i.e. removes $V_{i}$ from the mesh, updates vertex adjacency information and removes the two deleted elements from the element list.

Parallel coarsening begins with the initialisation of array $dynamic\_vertex$ which is defined as:

At the beginning, the whole array is initialised to -2, so that all mesh vertices will be considered for collapse.

In each iteration of the outer coarsening loop, $coarsen\_identify\_kernel$ is called for all vertices which have been marked for (re-)evaluation. Every vertex for which $dynamic\_vertex[V_{i}]\geq 0$ is said to be dynamic or active. At this point, a reduction in the total number of active vertices is necessary to determine whether there is anything left for coarsening or the algorithm should exit the loop.

Next up, we find the maximal independent set of active vertices $\mathcal{I}_{m}$. Working with independent sets not only ensures safe parallel execution, but also enforces the every other vertex rule. For every active vertex $V_{r}\in\mathcal{I}_{m}$ which is about to collapse, the local neighbourhood of all vertices $V_{a}$ formerly adjacent to $V_{r}$ changed and target vertices $dynamic\_vertex[V_{a}]$ may not be suitable choices any more. Therefore, when $V_{r}$ is erased, all its neighbours are marked for re-evaluation. This is how propagation of coarsening is implemented.

Algorithm 9 describes how the actual coarsening takes place in terms of modifications to mesh data structures. Updates which can lead to race conditions have been pointed out. These updates are deferred until the end of processing of the independent set. Before moving to the next iteration, all deferred operations are committed and colouring is repaired because edge collapse may have introduced inconsistencies.

The data dependencies in edge swapping are virtually identical to those of edge coarsening. Therefore, it is possible to reuse the same thread parallel algorithm as for coarsening in the previous section with slight modifications

In order to avoid maintaining edge-related data structures (e.g. edge-node list, edge-edge adjacency lists etc.), an edge can be expressed in terms of a pair of vertices. Just like in refinement, we define an edge $E_{ij}$ as a pair of vertices $(V_{i},V_{j})$, with $ID(V_{i})<ID(V_{j})$. We say that $E_{ij}$ is outbound from $V_{i}$ and inbound to $V_{j}$. Consequently, the edge $E_{ij}$ can be marked for swapping by adding $V_{j}$ to $marked\_edges[V_{i}]$. Obviously, a vertex $V_{i}$ can have more than one outbound edge, so unlike $dynamic\_vertex$ in coarsening, $marked\_edges$ in swapping needs to be a vector of sets (std::vector< std::set<int> >).

The algorithm begins by marking all edges. It then enters a loop which is terminated when no marked edges remain. The maximal independent set $\mathcal{I}_{m}$ of active vertices is calculated. A vertex is considered active if at least one of its outbound edges is marked. Following that, threads process all active vertices of $\mathcal{I}_{m}$ in parallel. The thread processing vertex $V_{i}$ visits all edges in $marked\_edges[V_{i}]$ one after the other and examines whether they can be swapped, i.e. whether the operation will improve the quality of the two elements sharing that edge. It is easy to see that swapping two edges in parallel which are outbound from two independent vertices involves no structural hazards.

Propagation of swapping is similar to that of coarsening. Consider the local patch in Figure 3 and assume that a thread is processing vertex $V_{0}$. If edge $\overline{V_{0}V_{1}}$ is flipped, the two elements sharing that edge change in shape and quality, so all four edges surrounding those elements (forming the rhombus in bold) have to be marked for processing. This is how propagation is implemented in swapping.

One last difference between swapping and coarsening is that $\mathcal{I}_{m}$ needs to be traversed more than once before proceeding to the next one. In the same example as above, assume that all edges adjacent to $V_{0}$ are outbound and marked. If edge $\overline{V_{0}V_{1}}$ is flipped, adjacency information for $V_{1}$, $V_{2}$ and $V_{3}$ has to be updated. These updates have to be deferred because another thread might try to update the same lists at the same time (e.g. the thread processing edge $\overline{V_{C}V_{1}}$). However, not committing the changes immediately means that the thread processing $V_{0}$ has a stale view of the local patch. More precisely, NEList[$V_{2}$] and NEList[$V_{3}$] are invalid and cannot be used to find what elements edges $\overline{V_{0}V_{2}}$ and $\overline{V_{0}V_{3}}$ are part of. Therefore, these two edges cannot be processed until the deferred operations have been committed. On the other hand, the rest of $V_{0}$’s outbound edges are free to be processed. Once all threads have processed whichever edges they can for all vertices of the independent set, deferred operations are committed and threads traverse the independent set again (up to two more times in 2D) to process what had been skipped before.

Algorithm 10 illustrates the colouring based algorithm for mesh smoothing. In this algorithm the graph $\mathcal{G}(\mathcal{V},\mathcal{E})$ consists of sets of vertices $\mathcal{V}$ and edges $\mathcal{E}$ that are defined by the vertices and edges of the computational mesh. By computing a vertex colouring of $\mathcal{G}$ we can define independent sets of vertices, $\mathcal{V}^{c}$, where $c$ is a computed colour. Thus, all vertices in $\mathcal{V}^{c}$, for any $c$, can be updated concurrently without any race conditions on dependent data. This is clear from the definition of the smoothing kernel in Section 2.4.4. Hence, within a node, thread-safety is ensured by assigning a different independent set $\mathcal{V}^{c}$ to each thread.