H(div)-Mixed FEM

[rezgarshakeri]

We would like to implement H(div)-Mixed FEM in libCEED. Based on the H(div) 2D example implementation here.

1. We need to be able to evaluate an inner product of the form (q, div(u) )
2. Since the dof of the element is defined as a normal on its edge/face (in 3D) thus element restriction operator contains +1 (when local and global dof are in same direction) and -1 (in opposite direction).

Note that the Nodal basis and its divergence are computed using the prime basis which we need to evaluate in the physical coordinate system.

@YohannDudouit @knepley @nocollier

[YohannDudouit]

I think H(div) should already be supported as a non-tensor basis (you might have to put each vectorial component in a different scalar field), similarly to other basis the interpolation operator B and gradient interpolation operator G should be provided by the user. Then, one can evaluate grad(u) using the standard interface, and then extract div(u) in the QFunction at each quadrature point. The Piola transformation should also be applied at the QFunction level.

Evaluating H(div) using the tensor structure of tensor Nedelec basis would require changes in libCEED as we currently assume the same basis and the same number of basis functions in each direction. If this is something you would be interested to implement I can provide you insights, but you should know that this will require to add new interface functions, and rewriting partially all the backends.

[jedbrown]

Thanks @YohannDudouit. For this project, we're interested in fairly low order, where tensor structure is not necessary, but we'd want to be able to extend to the tensor Nedelec elements.

What you're describing sounds like using an H^1 basis (perhaps with reduced continuity), versus using an H(div) basis (a bigger space than H^1). We'd like to support CEED_EVAL_DIV for efficiency and to keep the computed values bounded. But your suggestion to hack it like H^1 would be a short-term solution and provide a way to test.

[YohannDudouit]

I don't think it's an issue that non-tensor H(div) basis are bigger than Hˆ1, the libCEED interface is more or less agnostic of this. Whereas for the continuity, I thought it was enforced by the Piola transformation, which should be applied in the QFunction. Am I missing something?

Supporting CEED_EVAL_DIV for currently supported basis should be fairly straightforward for the cuda-shared and cuda-gen backends. I'm not sure about the other backends.

If supporting tensor Nedelec basis is something we want to do in the near future, we should probably organize some sort of technical meeting to see who would be available to work on the different aspects, brainstorm, and spread the work as it will have significant impact on most of the backends' code.

[jeremylt]

I think we accomplished this? Please re-open if there's more to do