diff --git a/fem/fe/fe_base.cpp b/fem/fe/fe_base.cpp
index 63da7e4428..67d13aad87 100644
--- a/fem/fe/fe_base.cpp
+++ b/fem/fe/fe_base.cpp
@@ -2535,11 +2535,16 @@ const DofToQuad &NodalTensorFiniteElement::GetDofToQuad(
    DofToQuad::Mode mode,
    const ElementDofOrdering ordering) const
 {
+   MFEM_VERIFY(!(mode == DofToQuad::Mode::TENSOR &&
+                 ordering == ElementDofOrdering::NATIVE),
+               "Invalide combination of DofToQuad::Mode and ElementDofOrdering.");
    for (int i = 0; i < dof2quad_array.Size(); i++)
    {
       const DofToQuad &d2q = *dof2quad_array[i];
       if (d2q.IntRule == &ir && d2q.mode == mode && d2q.ordering == ordering) { return d2q; }
    }
+   //First create the DofToQuad map for ElementDofOrdering::NATIVE.
+   //Either return d2q, or create the FULL and LEXICOGRAPHIC map and return.
    auto &d2q = GetDofToQuad(ir, mode);
    if (mode == DofToQuad::Mode::FULL &&
        ordering == ElementDofOrdering::LEXICOGRAPHIC)
@@ -2548,17 +2553,84 @@ const DofToQuad &NodalTensorFiniteElement::GetDofToQuad(
       auto *d2q_new = new DofToQuad(d2q);
       d2q_new->ordering = ElementDofOrdering::LEXICOGRAPHIC;
       const int nqpt = ir.GetNPoints();
-      for (int i = 0; i < nqpt; i++)
+
+
+      if (range_type == SCALAR)
       {
-         for (int d = 0; d < dim; d++)
+         for (int i = 0; i < nqpt; i++)
          {
             for (int j = 0; j < dof; j++)
             {
-               d2q_new->G[i+nqpt*(d+dim*j)] = d2q_new->Gt[j+dof*(i+nqpt*d)] = d2q.G[i+nqpt*
-                                                                                    (d+dim*dof_map[j])];
+               d2q_new->B[i+nqpt*j] = d2q_new->Bt[j+dof*i] = d2q.B[i+nqpt*dof_map[j]];
             }
          }
       }
+      else if (range_type == VECTOR)
+      {
+         for (int i = 0; i < nqpt; i++)
+         {
+            for (int d = 0; d < dim; d++)
+            {
+               for (int j = 0; j < dof; j++)
+               {
+                  d2q_new->B[i+nqpt*(d+dim*j)] = d2q_new->Bt[j+dof*(i+nqpt*d)] = d2q.B[i+nqpt*
+                                                                                       (d+dim*dof_map[j])];
+               }
+            }
+         }
+      }
+      else
+      {
+         // Skip B and Bt for unknown range type
+      }
+      switch (deriv_type)
+      {
+         case GRAD:
+         {
+            for (int i = 0; i < nqpt; i++)
+            {
+               for (int d = 0; d < dim; d++)
+               {
+                  for (int j = 0; j < dof; j++)
+                  {
+                     d2q_new->G[i+nqpt*(d+dim*j)] = d2q_new->Gt[j+dof*(i+nqpt*d)] = d2q.G[i+nqpt*
+                                                                                          (d+dim*dof_map[j])];
+                  }
+               }
+            }
+            break;
+         }
+         case DIV:
+         {
+            for (int i = 0; i < nqpt; i++)
+            {
+               for (int j = 0; j < dof; j++)
+               {
+                  d2q_new->G[i+nqpt*j] = d2q_new->Gt[j+dof*i] = d2q.G[i+nqpt*dof_map[j]];
+               }
+            }
+            break;
+         }
+         case CURL:
+         {
+            for (int i = 0; i < nqpt; i++)
+            {
+               for (int d = 0; d < cdim; d++)
+               {
+                  for (int j = 0; j < dof; j++)
+                  {
+                     d2q_new->G[i+nqpt*(d+cdim*j)] = d2q_new->Gt[j+dof*(i+nqpt*d)] = d2q.G[i+nqpt*
+                                                                                           (d+cdim*dof_map[j])];
+                  }
+               }
+            }
+            break;
+         }
+         case NONE:
+         default:
+            // Skip G and Gt for unknown derivative type
+            break;
+      }
       dof2quad_array.Append(d2q_new);
       return *d2q_new;
    }
diff --git a/miniapps/solvers/block_fespace_operator.hpp b/miniapps/solvers/block_fespace_operator.hpp
index 87b05126fd..260f319c2b 100644
--- a/miniapps/solvers/block_fespace_operator.hpp
+++ b/miniapps/solvers/block_fespace_operator.hpp
@@ -46,6 +46,9 @@ public:
    const Operator* GetProlongation () const override;
    const Operator* GetRestriction () const override;
    void Mult(const Vector &x, Vector &y) const override {A.Mult(x,y);};
+   /// @brief Wraps BlockOperator::SetBlock. Eventually would like this class to inherit
+   /// from BlockOperator instead, but can't easily due to ownership of offset data
+   /// in BlockOperator being by reference.
    void SetBlock( int   iRow,
                   int   iCol,
                   Operator *  op,