Clustering cells by similarity.

[bangerth]

@GrahamBenHarper gave a talk in our seminar today that made me think about how we can do our cell similarity things better.

What he is doing is this, in essence: Up front, all cells are classified as "similar" (to a certain tolerance) to other cells by clustering cells based on the difference in the Jacobian matrices of their mappings. So translated cells have a small distance and consequently end up in the same cluster. Then, instead of re-computing mappings, shape values, shape gradients, etc., on every cell, he only looks up the corresponding values for each cluster. In his examples, the number of clusters on many meshes is not very large, a small fraction of the total number of cells, so the fact that he doesn't have to recompute a bunch of things really pays off.

This whole thing made me think of our scheme in FEValues: We also consider cell similarity, but only to the last cell we visited. And we don't do it at all when we run in multithreaded mode because then the order in which we visit cells is no longer predictable. But if we clustered all cells up front, then we could switch the cell similarity test back on, and only compute the FEValues information once for each cluster, using a deterministically chosen "representative cell" for each cluster.

[bangerth]

In practice, many meshes will about as many clusters as there are cells, if they are totally unstructured. In cases like this, it may be useful to see whether we can identify -- say -- 20 or 50 clusters. If a cell is not part of those, then it is simply "unclustered" and no caching can happen on such cells. Of course, on structured meshes all or nearly all cells are parts of a very small number of clusters.

[GrahamBenHarper]

That's right. I'll add a couple things on here. It doesn't have to be just clustering. Really any sort of structure/similarity measuring or detection works, even AI/ML-based approaches 🙂 You just have to look at the reference to physical mapping's Jacobian. It looks like you use CellSimilarity::translation to identify this, since the Jacobians are exactly equal in the case of pure translations.

Clustering (depending on which clustering algorithm you choose) doesn't necessarily guarantee you'll have perfect translations. The greedy approach we've also followed constructively builds a list based on an upfront tolerance $\varepsilon$ in case errors are a major concern. Loop over all mesh cells. Loop over everything in the compressed list that starts empty. If you find a match on the compressed list that satisfies something like $\int_{\hat{T}} \|J_i - J_j\|_F<\varepsilon$, write down that lookup index. Otherwise, append the current mesh cell candidate to the compressed list. It's like an $O(nm)$ algorithm where you have $n$ mesh cells and the final amount of unique shapes is $m$. There are ways to improve that, but it's not terrible if you know upfront that $m$ will be small. And that's right, all of this can be done upfront as soon as you're done modifying the mesh but before you start any computations with your FEValues.

The patch preconditioner version (this computes an $\ell_1$ difference in Ifpack2) can be seen in:
https://github.com/trilinos/Trilinos/blob/58b8b785f4db5aa7a7fc960b87639287b9eb13ac/packages/ifpack2/src/Ifpack2_DatabaseSchwarz_def.hpp#L496-L552

There are plenty of other little optimizations to be made in there, but I want to switch away from those Teuchos::SerialDenseMatrix objects to Kokkos first 😱

Also, as a final note, I think deal.ii also has some support for generating extruded meshes. I don't know if you do this already, but it crossed my mind that every time you generate an extruded mesh, you'll get CellSimilarity::translation for free for the majority of the mesh as long as you number your mesh cells to index over the extruded direction first.