Hello Christopher,
Thanks for your reply. My application is a mixed higher-order space-time formulation of a parabolic PDE which needs a higher order H(div)-conforming function space, which was exactly our reason for choosing NGSolve.
I was unable to get a "clean" solution working (using a TensorProductFESpace) in Python, and I was unable to build NGSolve from source so I could work in C++, so I resorted to building my necessary matrices in space and time separately, and build the space-time system matrix using Kronecker products. In building one of my space matrices, I hit the bug mentioned above.
Indeed, one mesh is a refinement of the other, so my hope was to somehow get the prolongation operator from coarser to finer mesh and be on my merry way. It seems that this method is not exposed in the Python API (only a ProlongateInline() method, which I never got to work). I was able to work around this (in 2D at least) by creating a number of GridFunctions with setting `nested=True` and unit vectors as data, refining the mesh, Update()'ing the GridFunctions, and reading the resulting vectors, but this only seemed to work for lowest-order spaces.
I would love to get into the nitty gritty of NGSolve using C++, however I was unable to build it from source on both of my machines (one OSX 10.14, one OSX 10.12). If I manage to fix this problem, my next NGSolve project will work with C++ under the hood
Sorry for the lengthy reply!