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:
Change the connectivity of the mesh
Change the positions of the vertices in the mesh
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
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:
compute the quality of each element in the mesh
each vertex takes at the qualities of its surrounding elements and computes the minimum
each vertex computes the gradient of this minimum quality with respect to its position
nudge each non-boundary vertex along the gradient
repeat until mesh quality has improved enough or (less likely) until convergence
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
The gradient of
and the second part is taking the gradient of the quality metric
The derivatives of
The update procedure for gradient ascent looks like
where
From dimensional analysis, we can see that
where
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.
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)
Here's a small example implementation and a unit test (requires C++17):
x
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.