CeedElemRestriction for H(curl)

[sebastiangrimberg]

High-order H(curl) spaces (p > 1) in 3D involve pairs of element face dofs which must be treated specially in the L-vector to E-vector transformation. These transformations are more complicated than the simple sign flip accounted for in #889. In general, the pairs of dofs are transfomed by 2x2 blocks with entries {-1, 0, 1} to write the transformed E-vector dofs as linear combinations of the L-vector dofs. Thus, if each pair of face dofs are arranged to be consecutive in the element dof ordering, we can represent the transformation required as a tridiagonal matrix with entries {-1, 0, 1} applied to the E-vector. Note this is a (more expensive) generalization of the simple orientation flip case where the transformation matrix is just diagonal with entries {-1, 1}. Note that an alternative approach is to apply some global ordering of the nodes of the mesh as a preprocessing step after which a simple signed restriction suffices (like in https://epubs.siam.org/doi/10.1137/08073901X). However this strategy has its own shortcomings in particular when doing things like AMR.

More information can be found in this paper: https://dl.acm.org/doi/pdf/10.1145/3524456 (section 4.1.4, for example). See also mfem/mfem#1046 and https://github.com/mfem/mfem/blob/master/fem/doftrans.hpp for examples of how this transformation is handled in MFEM.

This PR adds CeedElemRestrictionCreateCurlOriented, which is similar to the CeedElemRestrictionCreateOriented constructor introduced in #889. Currently only implemented for CPU backends, with the main change coming in CeedElemRestrictionApply_Ref_Core.

Also adds some missing methods which should make blocked and opt backends work for oriented restrictions.

TODO

• Open questions:
  • The tridiagonal matrix-vector product in CeedElemRestrictionApply_Ref_Core is kind of ugly and likely slow. Simplifying this might make the eventual GPU backends for oriented and curl-oriented element restriction easier.
  • The storage of curl_orients could be optimized because we know the entries are only {-1, 0, 1}. So we could pack the tridiagonal matrix into an array of char, for example, where three pairs of bits starting from the least significant represent T[i-1,i], T[i,i], and T[i+1,i].
• Add unit tests
• Update Python/Julia/Fortran/Rust APIs

[nbeams]

Hi @sebastiangrimberg , I was hoping you could help me understand what exactly you meant by

pairs of element face dofs which must be treated specially in the L-vector to E-vector transformation.

I was thinking you meant the pairs of face dofs which are located at the same place (the two tangential components at that location on the face). But what about orders higher than 2? I read through the Scroggs et al. paper you linked above, and based on my understanding of that paper, wouldn't having 3 or more edge dofs on an edge mean that if we need to flip the edge, we may need to swap dofs which are not necessarily consecutively ordered? Did you mean that dofs on edges would also be paired in ordering based on how they would be swapped if the edge were flipped, rather than going "in order" along the edge? But also, what about rotations/reflections of a 3rd-or-higher order face? What could the transformation matrix look like there?

(In my head I was picturing elements like figure 3.10 in http://www.logg.org/anders/pub/papers/KirbyLoggEtAl2012a.pdf,
then thinking about how the dofs would be ordered for 3rd order or higher)

[sebastiangrimberg]

Hi @nbeams, I'm sorry for the confusing explanation and let me see if I can make my thinking more clear in another go.

In the context of Figure 3.10 in the reference you provided, you can see for order > 3 we have multiple co-located face dofs for each triangular face, but they still come in pairs. The transformation matrix for each face will include a 2x2 block for each pair of face dofs. These transformations only swap each pair of face dofs, there are no couplings between pairs so that's why I say consecutive ordering in each pair will lead the entire transformation matrix for the element to appear with 2x2 blocks. For the edge degrees of freedom, each can be oriented with a simple sign flip.

Maybe specifically in the case of Figure 3.10, each triangular face of the tetrahedon has 15 dofs (9 edge dofs, 6 face dofs). If we order the element dofs as edge dofs followed by face dofs, I believe the transformation matrix would contain a 9x9 block which is diagonal with entries {-1, +1}, followed by 3, 2x2 blocks with entries {-1, 0, 1} orienting each pair of face dofs. The transformation for the entire element would combine the transformation for all faces.

To see this in code, this method performs the transformation for triangular elements in MFEM. You can see that for a triangle element of order p, we have p*(p-1) face dofs which are treated in pairs, multiplied by 2x2 dense matrices in the transformation depending on the orientation type of the face (which can take on one of 6 values depending on the ordering of the nodes).

[nbeams]

Thanks for the detailed response.

you can see for order > 3 we have multiple co-located face dofs for each triangular face, but they still come in pairs

Yes, I understand that part. My confusion is in lining this up with what is presented in the Scroggs paper, unless you do not mean for this PR to be an implementation of exactly that method, and/or I do not understand the paper fully. E.g., if we have 4 dofs on an edge, numbered 1-2-3-4, and we flip that edge, wouldn't that swap dofs 1 and 4 according to the paper? Even though they both just flip in direction. Or if we "reflected" the face for a high-order element, then wouldn't some pairs be swapped with each other (but of course they still remain "partnered" in their original pairs).

The linked code for MFEM does indeed seem to only be considering the face dofs, and skipping the edge dofs (everything starts after [number of edges] * nedofs).

[sebastiangrimberg]

Ah, I think I understand now. Indeed I don't mean to literally implement the exact matrix M from the paper as it is generally dense in their implementation (as you point out). But what I'm recognizing is that you can separate into two parts, the permutation of numbering (which is certainly not tridiagonal) followed but the orientation changes and rotations affecting pairs of face dofs (which is tridiagonal). Thus the same effective M as provided by a library like Basix should be compatible with this interface in libCEED by providing the correct permutations and tridiagonal transformation matrix.

Also, you are right that in MFEM the orientation changes are handled outside of this DofTransformation class, which was added relatively recently. Hence the reason why the code I pointed to is only touching face dofs. I do understand this is slightly confusing but has to do with the order in which things were added to the library I suppose. In the case where you have a nice tensor product element mesh with no triangular faces, you no longer need the general face dof transformations and can make do with only +/- sign flips of your dofs in the element restriction.

[nbeams]

Ah, thanks. That makes sense. I agree that a combination of offsets and a tridiagonal curl_orients could be used to describe the necessary transformations, provided the user had the information available in the correct format.

In some ways, the interface of libCEED makes this both "easier" and "harder" -- with libCEED (unlike a library like MFEM, which will create the mesh for you from high-level commands), it's the user's responsibility to give us correct information about the mesh and element transformations. But on the other hand, that means we need to make sure we're very clear about exactly what we expect the user to provide and how... I personally think the documentation for CurlOriented element restrictions would need to be improved before we add it to the interface, so users are clear on precisely what information they need to be able to get out of their mesh in order to use it correctly.

It seems this format would be straightforward to use with MFEM integration, but what about, e.g., our PETSc integration? (That is not a question for you, specifically, but more a general discussion question for anyone)

[sebastiangrimberg]

Here is an example of usage, from a draft MFEM PR (which has been tested):

• PR: libCEED interface: Partial and matrix-free assembly support for H(div) and H(curl) vector FE spaces mfem/mfem#3626
• Restriction construction: https://github.com/mfem/mfem/pull/3626/files#diff-66c54777b9e833630ecf5b28affce11129cb846ae6b0f6e6b263489ad55d0e34R281

[sebastiangrimberg]

I'll point you to the specific wording when I'm back at the desk in the morning. The convention is captured in tests/junit.py but that's not easy to parse at first glance

I got it, see 39c9572 and that fixed the failing OCCA tests. After some more GPU fixes in 52aff1d all the correctness tests are passing.

[sebastiangrimberg]

Once I get the go ahead and this is approved for merging, I will squash WIP/revert commits and force push if that sounds like a good plan.

[jedbrown]

Your squash plan sounds good. I think we also need

• add entry to doc/sphinx/source/releasenotes.md
• add your name to AUTHORS

[sebastiangrimberg]

See #1265 for the rebased + squashed branch, to be merged.

[jedbrown]

Thanks! What do you think of the claimed code coverage? https://app.codecov.io/gh/CEED/libCEED/pull/1168/blob/backends/ref/ceed-ref-restriction.c

There are sizable loop/indexing blocks that are purportedly not tested. Do you think the tests cover those or would you agree they are not covered? (If it's a codecov issue, we can fix later, but those are pretty big blocks to actually be untested/unreachable.)

[sebastiangrimberg]

I think this is right. Let me see if I can improve this by playing around with the unit tests a bit.

[sebastiangrimberg]

There are sizable loop/indexing blocks that are purportedly not tested. Do you think the tests cover those or would you agree they are not covered? (If it's a codecov issue, we can fix later, but those are pretty big blocks to actually be untested/unreachable.)

I think ae4e5d7 should offer a lot of coverage improvement for those loops, though I have to admit I'm not fully sure how to check. Looking at https://app.codecov.io/gh/CEED/libCEED/pull/1168/blob/backends/ref/ceed-ref-restriction.c doesn't show any data. No rush here, maybe it just takes time to populate? Thanks for flagging this!

[jedbrown]

Coverage looks off in this PR, but it looks reasonable (and much improved) for #1265: https://app.codecov.io/gh/CEED/libCEED/pull/1265

What do you think of the remaining red blocks? I'm not too concerned about coverage of the Fortran interface, but there's still a couple functions and loop nests

[sebastiangrimberg]

Coverage looks off in this PR, but it looks reasonable (and much improved) for #1265: https://app.codecov.io/gh/CEED/libCEED/pull/1265

What do you think of the remaining red blocks? I'm not too concerned about coverage of the Fortran interface, but there's still a couple functions and loop nests

Nice, I just added a new test t571 which should hopefully fix the remaining large untested code paths for curl-oriented element restrictions (fingers crossed).

[jeremylt]

I'm not worried about the missed lines in interface/ceed-elemrestriction.c and interface/ceed-fortran.c but the ones in the opt and blocked backends are strange.

[sebastiangrimberg]

Closed in favor of #1265.