smish.dev
mesh_relaxation

Element Quality

The finite element method cares about what the elements in its mesh look like. The presence of long, skinny elements in a mesh can make analysis take significantly longer, and make the solutions less accurate. So, while it would be nice if all the elements in our meshes were equilateral triangles, it's just not a practical expectation, since creating meshes is really hard to do (especially robustly). As a result, the meshes we get frequently contain "bad" elements, and we have to do something about it. Broadly, there are 3 ways to deal with the "bad" elements in a mesh:

  1. Change the connectivity of the mesh

  2. Change the positions of the vertices in the mesh

  3. Do nothing

Option 3 is more common than we'd like to admit, but for analyses where the results matter, it's not a great idea. This writeup looks into option 2: a way to wiggle vertex positions that improve the element qualities.

So, what is "element quality"? It's a way to assign a number to each element, describing how "good" its shape is (e.g. numbers close to 0 are "bad" elements, and 1 is a "perfect" element shape). There are a lot of different functions that do this (see: this document for more details), here's one I like for triangles:

q():=43A()rms()2

where A() is the signed area of the triangle, and rms() is the root-mean-square of the edge lengths of the triangle. Here are a few examples of triangles along with their qualities calculated with the metric above:

test

Iterative Improvement

The basic idea of this algorithm is to move the vertices in an attempt to maximize the minimum element quality. Rather than computing the minimum quality over the entire mesh (which would only modify the vertices of one element per iteration), we instead compute minimum quality and perform gradient ascent at each vertex. The steps are:

  1. compute the quality of each element in the mesh

  2. each vertex takes at the qualities of its surrounding elements and computes the minimum

  3. each vertex computes the gradient of this minimum quality with respect to its position

  4. nudge each non-boundary vertex along the gradient

  5. repeat until mesh quality has improved enough or (less likely) until convergence

Smooth Minima

Implementing step 2 as written above wouldn't actually work very well because taking the minimum of a list of numbers isn't a smooth operation. Abrupt changes in the derivative make it harder for algorithms like gradient ascent to find a maximum. So, instead we compute a smooth minimum in step 2:

qmin:=smin(q1,q2,q3,...)=1αlog(exp(αq1)+exp(αq2)+...)

There are other ways to compute smooth minima (Inigo Quilez has a wonderful article on the topic here: https://iquilezles.org/articles/smin/), but this exponential one is nice because it's easy to differentiate (which will come in handy in the following section) and order independent. The parameter α controls the smoothness of the minimum operation (smaller α is smoother, α= gives the "hard" minimum, in theory). Since we know that {q1,q2,...} are all in the range [0, 1], α values of 5 - 10 ish tend to work well in practice.

Gradient Calculation

The gradient of qmin has two parts, since it is a composition of functions: the smooth min and the element quality function. The first part involves differentiating the smooth min

dqmindx=ddx(smin(q1,q2,q3,...))=ismindqidqidx=iexp(αqi)jexp(αqj)dqidx=iexp(αqi)dqidxiexp(αqi)

and the second part is taking the gradient of the quality metric

dqdx=ddx(43Arms2)=43rms2(dAdxArms2d(rms2)dx)

The derivatives of A and rms2 are easier to write with respect to all the vertices at once:

dAdX:=[dAdx1dAdy1dAdx2dAdy2dAdx3dAdy3]=[y2y3x3x2y3y1x1x3y1y2x2x1]
d(rms2)dX:=[d(rms2)dx1d(rms2)dy1d(rms2)dx2d(rms2)dy2d(rms2)dx3d(rms2)dy3]=[4x12x22x34y12y22y34x22x12x34y22y12y34x32x22x14y32y22y1]

Step Size Determination

The update procedure for gradient ascent looks like

xn+1:=xn+γq(xn)

where q is given in the previous section and γ is a parameter that controls the size of the step taken along the gradient. Taking γ too big leads to nonconvergence and taking it too small makes it take a lot of iterations to improve the objective function. So what should we pick for γ?

From dimensional analysis, we can see that x, and qQ, where has dimension of "length" and Q has dimension of "quality". So, γ must have the dimension of 2Q, which means we need to first figure out characteristic length and quality scales for the problem. Quality is easier since it is already normalized, so its characteristic scale is just unity (Qc=1). There are a bunch of ways to get a characteristic length from a mesh, but I just use the average edge length. From here, we can set

γ:=γ0c2Qc=γ0¯2

where γ0 is some dimensionless parameter, which I usually take somewhere between [0.05, 0.10].

Examples

Using a simple mesh of the unit disk as an example, if we apply this iteration a few times, we can observe how the element qualities (shown on the right, in sorted order) improve:

 

Another benefit of this choice of quality metric is that it can be applied to inverted elements as well. Note: the negative quality values correspond to the initially inverted elements.

From these animations, we can see that the first few iterations have the most impact on improving the minimum element quality. For this reason, I typically just apply around 5-10 iterations of this procedure for most meshes, rather than actually iterating to convergence. That being said, this iterative procedure is relatively inexpensive, so if you have an important mesh that sees a lot of use, it can be worthwhile to apply more iterations.

 

Other Considerations

 

C++ Implementation

Here's a small example implementation and a unit test (requires C++17):

If you're interested in using this in a real application, you'd probably want to just use a library like glm for the vector/matrix arithmetic and also perform the triangle loop in vertex_smoothing() in parallel.