diff --git a/include/solid_mechanics/mf_nh_operator.h b/include/solid_mechanics/mf_nh_operator.h
index a0395d8..cccfc3e 100644
--- a/include/solid_mechanics/mf_nh_operator.h
+++ b/include/solid_mechanics/mf_nh_operator.h
@@ -871,10 +871,18 @@ NeoHookOperator<dim, Number>::do_operation_on_cell(FECellIntegrator &phi) const
                       {
                         VectorizedArrayType sum =
                           inv_jac[e][0] *
+#if DEAL_II_GTE(9, 6, 0)
                           ref_grads[(d * n_q_points + q) * dim + 0];
+#else
+                          ref_grads[(d * dim + 0) * n_q_points + q];
+#endif
                         for (unsigned int f = 1; f < dim; ++f)
                           sum += inv_jac[e][f] *
+#if DEAL_II_GTE(9, 6, 0)
                                  ref_grads[(d * n_q_points + q) * dim + f];
+#else
+                                 ref_grads[(d * dim + f) * n_q_points + q];
+#endif
                         F[d][e] = sum;
 
                         // since we already have the inverse Jacobian,
@@ -883,10 +891,18 @@ NeoHookOperator<dim, Number>::do_operation_on_cell(FECellIntegrator &phi) const
                         // same otherwise)
                         VectorizedArrayType sum2 =
                           inv_jac[e][0] *
+#if DEAL_II_GTE(9, 6, 0)
                           x_grads[(d * n_q_points + q) * dim + 0];
+#else
+                          x_grads[(d * dim + 0) * n_q_points + q];
+#endif
                         for (unsigned int f = 1; f < dim; ++f)
                           sum2 += inv_jac[e][f] *
+#if DEAL_II_GTE(9, 6, 0)
                                   x_grads[(d * n_q_points + q) * dim + f];
+#else
+                                  x_grads[(d * dim + f) * n_q_points + q];
+#endif
                         grad_Nx_v[d][e] = sum2;
                       }
                     F[d][d] += one;