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:

where $A(\mathrm{△})$$A(\triangle)$ is the signed area of the triangle, and ${\ell }_{rms}(\mathrm{△})$$\ell_{rms}(\triangle)$ 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:

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:

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 $\alpha$$\alpha$ controls the smoothness of the minimum operation (smaller $\alpha$$\alpha$ is smoother, $\alpha =\mathrm{\infty }$$\alpha = \infty$ gives the "hard" minimum, in theory). Since we know that $\{q_1,q_2,...\}$$\{q_1, q_2, ...\}$ are all in the range [0, 1], $\alpha$$\alpha$ values of 5 - 10 ish tend to work well in practice.

Gradient Calculation

The gradient of $q_{\text{min}}$$q_{\text{min}}$ 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

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

The derivatives of $A$$A$ and ${\ell }_{rms}^2$$\ell_{rms}^2$ are easier to write with respect to all the vertices at once:

Step Size Determination

The update procedure for gradient ascent looks like

where $\mathrm{\nabla }q$$\nabla q$ is given in the previous section and $\gamma$$\gamma$ is a parameter that controls the size of the step taken along the gradient. Taking $\gamma$$\gamma$ 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 $\gamma$$\gamma$?

From dimensional analysis, we can see that $x\propto \ell$$x \propto \ell$, and $\mathrm{\nabla }q\propto \frac{Q}{\ell }$$\nabla q \propto \frac{Q}{\ell}$, where $\ell$$\ell$ has dimension of "length" and $Q$$Q$ has dimension of "quality". So, $\gamma$$\gamma$ must have the dimension of $\frac{{\ell }^2}{Q}$$\frac{\ell^2}{Q}$, 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 ($Q_c=1$$Q_c = 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

where ${\gamma }_0$$\gamma_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

• This idea generalizes immediately to tetrahedron meshes as well, and possibly other element types for an appropriate choice of quality metric.

• Boundary nodes (except "corners") can also be allowed to move each iteration by either:

  • projecting back onto the original geometry

  • constraining the position changes along the local tangent plane

• It's also possible to mark/fix interior vertices to ensure the iterative updates preserve important mesh features (e.g. changes in material)

C++ Implementation

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

​x
#include <cmath>
#include <array>
#include <tuple>
#include <vector>
#include <iostream>

using int3 = std::array<int,3>;
using double2 = std::array<double,2>;
using triangle = std::array<double2,3>;

struct triangle_mesh {
  std::vector< double2 > vertices;
  std::vector< int3 > triangles;
};

inline double dot(double2 x, double2 y) {
  return x[0] * y[0] + x[1] * y[1];
}

inline double det(double2 x, double2 y) {
  return x[0] * y[1] - x[1] * y[0];
}

inline double2 operator-(double2 x, double2 y) {
  return {x[0] - y[0], x[1] - y[1]};
}

inline double2 operator*(double2 x, double scale) {
  return {x[0] * scale, x[1] * scale};
}

inline triangle operator*(triangle t, double scale) {
  return triangle{{{t[0][0] * scale, t[0][1] * scale},
                   {t[1][0] * scale, t[1][1] * scale},
                   {t[2][0] * scale, t[2][1] * scale}}};
}

inline triangle operator/(triangle t, double scale) {
  return triangle{{{t[0][0] / scale, t[0][1] / scale},
                   {t[1][0] / scale, t[1][1] / scale},
                   {t[2][0] / scale, t[2][1] / scale}}};
}

inline triangle operator-(triangle x, triangle y) {
  return triangle{{{x[0][0] - y[0][0], x[0][1] - y[0][1]},
                   {x[1][0] - y[1][0], x[1][1] - y[1][1]},
                   {x[2][0] - y[2][0], x[2][1] - y[2][1]}}};
}

inline void operator+=(double2 & x, double2 y) {
  x[0] += y[0];
  x[1] += y[1];
}

std::tuple<double, triangle> quality_and_gradient(const triangle & x) {

  double2 L01 = x[1] - x[0];
  double2 L12 = x[2] - x[1];
  double2 L20 = x[0] - x[2];

  double top = det(L20, L01);
  auto dtop_dx = triangle{{{ x[1][1] - x[2][1], -x[1][0] + x[2][0]},
                           {-x[0][1] + x[2][1],  x[0][0] - x[2][0]},
                           { x[0][1] - x[1][1], -x[0][0] + x[1][0]}}};

  double bot = dot(L01, L01) + dot(L12, L12) + dot(L20, L20);
  auto dbot_dx =
      triangle{{{ 4 * x[0][0] - 2 * x[1][0] - 2 * x[2][0],
                  4 * x[0][1] - 2 * x[1][1] - 2 * x[2][1]},
                {-2 * x[0][0] + 4 * x[1][0] - 2 * x[2][0],
                 -2 * x[0][1] + 4 * x[1][1] - 2 * x[2][1]},
                {-2 * x[0][0] - 2 * x[1][0] + 4 * x[2][0],
                 -2 * x[0][1] - 2 * x[1][1] + 4 * x[2][1]}}};

  constexpr double scale = 3.4641016151377543864;  // 2 * sqrt(3)

  return {
    (top / bot) * scale,
    (dtop_dx / bot - dbot_dx * (top / (bot * bot))) * scale
  };

}

void vertex_smoothing(triangle_mesh & mesh,              // mesh is modified in-place
                      const std::vector<int> & interior, // 1 = interior, 0 = boundary
                      double alpha,                      // smoothing parameter
                      double gamma,                      // gradient ascent step size
                      int num_iterations) {              

  int num_triangles = mesh.triangles.size();


  for (int k = 0; k < num_iterations; k++) {

    std::vector< double > scale(mesh.vertices.size(), 0.0);
    std::vector< double2 > grad(mesh.vertices.size(), double2{});

    for (int i = 0; i < num_triangles; i++) {

      int3 elem_ids = mesh.triangles[i];
      triangle tri = {
        mesh.vertices[elem_ids[0]],
        mesh.vertices[elem_ids[1]],
        mesh.vertices[elem_ids[2]]
      };

      auto [Q, dQdX] = quality_and_gradient(tri);

      double expQ = expf(-alpha * Q);

      dQdX = dQdX * expQ;
      for (int j = 0; j < 3; j++) {
        double2 g = dQdX[j];
        grad[elem_ids[j]] += g;
        scale[elem_ids[j]] += expQ;
      }
    };

    for (int i = 0; i < grad.size(); i++) {
      mesh.vertices[i] += grad[i] * (interior[i] * gamma / scale[i]);
    }

  }

}

int main() {

  double c = cos(M_PI / 3.0);
  double s = sin(M_PI / 3.0);

  // hexagon mesh with initially eccentric interior node
  //
  //       before                after
  //
  //       2 * * 1              2 * * 1
  //      *  *   **            * *   * *
  //     *     * * *          *   * *   *
  //    3--------6--0  --->  3-----6-----0
  //     *     * * *          *   * *   *
  //      *  *   **            * *   * *
  //       4 * * 5              4 * * 5
  //
  triangle_mesh mesh {
    {{1,0}, {c,s}, {-c,s}, {-1,0}, {-c,-s}, {c,-s}, {0.5,0}},
    {{0,1,6},{1,2,6},{2,3,6},{3,4,6},{4,5,6},{5,0,6}}
  };
  
  std::vector< int > interior = {0, 0, 0, 0, 0, 0, 1};

  double alpha = 10.0;
  double l_c = 1.0; // this specific mesh has characteristic length 1
  double gamma = 0.1f * l_c * l_c;

  for (int i = 0; i < 5; i++) {
    vertex_smoothing(mesh, interior, alpha, gamma, 1);
    std::cout << mesh.vertices[6][0] << " " << mesh.vertices[6][1] << std::endl;
  }
  
  // prints: 
  // 0.413154 7.82536e-19
  // 0.349373 3.77991e-18
  // 0.298876 9.43538e-18
  // 0.257198 8.07657e-18
  // 0.222032 8.77416e-18
  // ...
  //
  // so: interior vertex moves toward the optimal location (0,0) with each iteration

}

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.