To activate for instance for a P=3 simulation, use a Hex8 or Quad4 mesh and add a line under realms (e.g. below "use_edges")
-- "polynomial_order: 3"
Two output files will be written, a regular output with only the original Hex8/Quad4 nodes, and a "fine" output which outputs subelements of the promoted topology
Limited currently: only works with wall, inflow, and outflow bcs. No shifted gauss points.
Little effort to make the method P-scalable yet: setting P > 3 or 4 in 3D will make assembly time very large. The memory usage is also very, very large.