diff --git a/remhos.cpp b/remhos.cpp
index 51e5413..1bda644 100644
--- a/remhos.cpp
+++ b/remhos.cpp
@@ -574,12 +574,8 @@ int main(int argc, char *argv[])
    VectorGridFunctionCoefficient v_sub_coef;
    Vector x0_sub;
 
-   if (order > 1)
+   if (use_subcell_RD)
    {
-      // The mesh corresponding to Bezier subcells of order p is constructed.
-      // NOTE: The mesh is assumed to consist of quads or hexes.
-      MFEM_VERIFY(order > 1, "This code should not be entered for order = 1.");
-
       // Get a uniformly refined mesh.
       const int btype = BasisType::ClosedUniform;
       subcell_mesh = new ParMesh(ParMesh::MakeRefined(pmesh, order, btype));
@@ -1037,11 +1033,13 @@ int main(int argc, char *argv[])
          x0.HostReadWrite(); v_sub_gf.HostReadWrite();
          x.HostReadWrite();
          add(x0, t, v_gf, x);
-         x0_sub.HostReadWrite(); v_sub_gf.HostReadWrite();
-         MFEM_VERIFY(xsub != NULL,
-                     "xsub == NULL/This code should not be entered for order = 1.");
-         xsub->HostReadWrite();
-         add(x0_sub, t, v_sub_gf, *xsub);
+         if (use_subcell_RD)
+         {
+            x0_sub.HostReadWrite(); v_sub_gf.HostReadWrite();
+            MFEM_VERIFY(xsub != NULL, "Subcell mesh not defined!");
+            xsub->HostReadWrite();
+            add(x0_sub, t, v_sub_gf, *xsub);
+         }
       }
 
       if (problem_num != 6 && problem_num != 7 && problem_num != 8)
@@ -1220,7 +1218,7 @@ int main(int argc, char *argv[])
    delete lom.pk;
    delete dc;
 
-   if (order > 1)
+   if (use_subcell_RD)
    {
       delete subcell_mesh;
       delete fec_sub;
diff --git a/remhos_tools.cpp b/remhos_tools.cpp
index a59247c..ddcc126 100644
--- a/remhos_tools.cpp
+++ b/remhos_tools.cpp
@@ -353,7 +353,8 @@ void SmoothnessIndicator::ComputeFromSparsity(const SparseMatrix &K,
 DofInfo::DofInfo(ParFiniteElementSpace &pfes_sltn, int btype)
    : bounds_type(btype),
      pmesh(pfes_sltn.GetParMesh()), pfes(pfes_sltn),
-     fec_bounds(pfes.GetOrder(0), pmesh->Dimension(), BasisType::GaussLobatto),
+     fec_bounds(std::max(pfes.GetOrder(0), 1),
+                pmesh->Dimension(), BasisType::GaussLobatto),
      pfes_bounds(pmesh, &fec_bounds),
      x_min(&pfes_bounds), x_max(&pfes_bounds)
 {