Feature request: physical shape-function gradients for RAVIART_THOMAS (H(div)) elements

[joe61vette]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Summary

Please implement HDivFETransformation::map_dphi in libMesh so that the physical
shape-function gradients (dphi) of RAVIART_THOMAS elements are computed via the
contravariant-Piola gradient mapping. Today the method is a no-op, so dphi is identically zero
for H(div) elements in both volume and face reinit.

Why it's needed

For my task, the gap is the viscous operator. A consistent primal Symmetric
Interior Penalty (SIP) discretization of ∇·(μ∇u) requires the physical gradient of the RT
test/trial functions (∇w, ∇u). Because map_dphi is a no-op, any kernel that contracts an RT
gradient contributes exactly zero — the SIP operator is inert.

At present, I have a working approximation (an orthogonal two-point flux on the reconstructed
velocity, the DG-vector analog of FVOrthogonalDiffusion, but it captures only the face-normal
gradient and is not optimal-order convergent on the simplex meshes RT requires. A real
map_dphi would replace it with a consistent, convergent SIP operator.

What is required

• Implement the contravariant-Piola gradient mapping in
  HDivFETransformation::map_dphi (currently only issues a libmesh_warning and returns).
  The reference-space derivatives (shape_deriv) already exist; what is missing is the mapping
  to physical space, including the Jacobian and its spatial derivatives.
• Populate dphi (and the component arrays) consistently in both volume and face reinit,
  matching how map_phi already maps the values.
• A minimal verification: on a curved/affine simplex, confirm the mapped dphi reproduces the
  analytic gradient of a known RT field (e.g. a patch test / MMS), rather than zero.

References

• libMesh: HDivFETransformation::map_dphi (the no-op) and HDivFETransformation::map_phi
  (the working value mapping) — fe_transformation_base / hdiv_fe_transformation.
• Existing precedent: the covariant Piola map for NEDELEC_ONE (H(curl)) implements its own
  transformations; H(div) has the value map but not the gradient map.

Thanks,
Joe (& Claude-Code)

[GiudGiud]

@lindsayad @roystgnr what are your thoughts on this? who implemented it for NEDELEC?

@joe61vette have you tried having claude give it a go? it has a clear test for checking the implementation, this "harness" is in my experience what makes it reliably successful

[lindsayad]

Joe's method seems to exactly correspond to https://link.springer.com/article/10.1007/s10915-006-9107-7 as well as I can tell. A symmetric interior penalty discontinuous Galerkin method with the H(div) velocity space achieved through selection of the RT element whereas my IP-HDG method implementation obtains an H(div) velocity space through hybridization. So the gradient of the H(div) functions are always integrated over element interiors and never over faces; jumps in the tangential components of the velocity are handled in the typical symmetric interior penalty DG way.

[lindsayad]

I would be happy to add map_dphi in libmesh for the H(div) transformation unless you've already started on it @joe61vette. Then I'd add yet another example to our NS FE suite for flow-transport compatible discretizations with again the transported scalar not tied to the FE space of the pressure Lagrange multiplier since the velocity space is exactly/pointwise/strongly divergence free.

@joe61vette do you think you'd ever contribute your porous flow finite element implementation to MOOSE?

[lindsayad]

I've opened libMesh/libmesh#4524 and assigned myself, but I'll assign to you @joe61vette if you've already begun work on the issue

[joe61vette]

Thanks Alex. It would be great if you could provide that capability. I get nervous working in MOOSE (I still owe you ADDGVectorKernel) and libmesh is really beyond my comprehension.

Also, once this app is fully operational and tested, I'll be able to make it available to you for possible inclusion in the NS module.

[lindsayad]

FYI @joe61vette: libMesh/libmesh#4529