diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp
index 16c659b..f4f23c2 100644
--- a/src/poly/utils/private/operator.cpp
+++ b/src/poly/utils/private/operator.cpp
@@ -54,24 +54,34 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
         return;
     }
 
-    m_M -> EliminateTestEssentialBC(m_theta_ess_tdofs.first);
-    m_M -> EliminateTrialEssentialBC(m_phi_ess_tdofs.first);
+    m_Mmat = std::make_unique<mfem::SparseMatrix>(m_M->SpMat());
+    m_Qmat = std::make_unique<mfem::SparseMatrix>(m_Q->SpMat());
+    m_Dmat = std::make_unique<mfem::SparseMatrix>(m_D->SpMat());
 
-    m_Q -> EliminateTestEssentialBC(m_phi_ess_tdofs.first);
-    m_Q -> EliminateTrialEssentialBC(m_theta_ess_tdofs.first);
+    // Remove the essential dofs from the constant matrices
+    for (const auto& dof: m_theta_ess_tdofs.first) {
+        m_Mmat->EliminateRow(dof);
+        m_Qmat->EliminateCol(dof);
+    }
 
-    m_D -> EliminateEssentialBC(m_phi_ess_tdofs.first);
+    for (const auto& dof: m_phi_ess_tdofs.first) {
+        m_Mmat->EliminateCol(dof);
+        m_Qmat->EliminateRow(dof);
+        m_Dmat->EliminateRowCol(dof);
+    }
 
-    m_negM_op = std::make_unique<mfem::ScaledOperator>(m_M.get(), -1.0);
-    m_negQ_op = std::make_unique<mfem::ScaledOperator>(m_Q.get(), -1.0);
+    m_negM_mat = std::make_unique<mfem::ScaledOperator>(m_Mmat.get(), -1.0);
+    m_negQ_mat = std::make_unique<mfem::ScaledOperator>(m_Qmat.get(), -1.0);
 
-    m_schurCompliment = std::make_unique<SchurCompliment>(*m_Q, *m_D, *m_M);
+    m_schurCompliment = std::make_unique<SchurCompliment>(*m_Qmat, *m_Dmat, *m_Mmat);
 
     // Set up the constant parts of the jacobian now
+
+    // We use the assembled matrices with their boundary conditions already removed for the jacobian
     m_jacobian = std::make_unique<mfem::BlockOperator>(m_blockOffsets);
-    m_jacobian->SetBlock(0, 1, m_negM_op.get()); //<- -M (constant)
-    m_jacobian->SetBlock(1, 0, m_negQ_op.get()); //<- -Q (constant)
-    m_jacobian->SetBlock(1, 1, m_D.get());       //<-  D (constant)
+    m_jacobian->SetBlock(0, 1, m_negM_mat.get()); //<- -M (constant)
+    m_jacobian->SetBlock(1, 0, m_negQ_mat.get()); //<- -Q (constant)
+    m_jacobian->SetBlock(1, 1, m_Dmat.get());     //<-  D (constant)
 
     m_invSchurCompliment = std::make_unique<GMRESInverter>(*m_schurCompliment);
 
@@ -240,9 +250,9 @@ void GMRESInverter::Mult(const mfem::Vector &x, mfem::Vector &y) const {
 
 
 SchurCompliment::SchurCompliment(
-    const mfem::MixedBilinearForm &QOp,
-    const mfem::BilinearForm &DOp,
-    const mfem::MixedBilinearForm &MOp,
+    const mfem::SparseMatrix &QOp,
+    const mfem::SparseMatrix &DOp,
+    const mfem::SparseMatrix &MOp,
     const mfem::Solver &GradInvOp) :
 mfem::Operator(DOp.Height(), DOp.Width())
 {
@@ -252,9 +262,9 @@ mfem::Operator(DOp.Height(), DOp.Width())
 }
 
 SchurCompliment::SchurCompliment(
-    const mfem::MixedBilinearForm &QOp,
-    const mfem::BilinearForm &DOp,
-    const mfem::MixedBilinearForm &MOp) :
+    const mfem::SparseMatrix &QOp,
+    const mfem::SparseMatrix &DOp,
+    const mfem::SparseMatrix &MOp) :
 mfem::Operator(DOp.Height(), DOp.Width())
 {
     updateConstantTerms(QOp, DOp, MOp);
@@ -262,7 +272,7 @@ mfem::Operator(DOp.Height(), DOp.Width())
     m_nTheta = m_MOp->Height();
 }
 
-void SchurCompliment::SetOperator(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp, const mfem::Solver &GradInvOp) {
+void SchurCompliment::SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp) {
     updateConstantTerms(QOp, DOp, MOp);
     updateInverseNonlinearJacobian(GradInvOp);
 }
@@ -271,7 +281,7 @@ void SchurCompliment::updateInverseNonlinearJacobian(const mfem::Solver &gradInv
     m_GradInvOp = &gradInv;
 }
 
-void SchurCompliment::updateConstantTerms(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp) {
+void SchurCompliment::updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp) {
     m_QOp = &QOp;
     m_DOp = &DOp;
     m_MOp = &MOp;
diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h
index 1cffe05..051845f 100644
--- a/src/poly/utils/public/operator.h
+++ b/src/poly/utils/public/operator.h
@@ -28,22 +28,23 @@
 
 class SchurCompliment final : public mfem::Operator {
 public:
-    SchurCompliment(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp, const mfem::Solver &GradInvOp);
-    SchurCompliment(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp);
+    SchurCompliment(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp);
+    SchurCompliment(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp);
     ~SchurCompliment() override = default;
     void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
-    void SetOperator(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp, const mfem::Solver &GradInvOp);
+    void SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &
+                     GradInvOp);
     void updateInverseNonlinearJacobian(const mfem::Solver &gradInv);
 
 private:
-    void updateConstantTerms(const mfem::MixedBilinearForm &QOp, const mfem::BilinearForm &DOp, const mfem::MixedBilinearForm &MOp);
+    void updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp);
 
 private:
 
     // Note that these are not owned by this class
-    const mfem::MixedBilinearForm* m_QOp = nullptr;
-    const mfem::BilinearForm* m_DOp = nullptr;
-    const mfem::MixedBilinearForm* m_MOp = nullptr;
+    const mfem::SparseMatrix* m_QOp = nullptr;
+    const mfem::SparseMatrix* m_DOp = nullptr;
+    const mfem::SparseMatrix* m_MOp = nullptr;
     const mfem::Solver* m_GradInvOp = nullptr;
     int m_nPhi = 0;
     int m_nTheta = 0;
@@ -101,32 +102,33 @@ private:
     std::unique_ptr<mfem::BilinearForm> m_D;
     std::unique_ptr<mfem::NonlinearForm> m_f;
 
+    // These are used to store the matrix representations
+    // for the bi-linear forms. This might seem counterintuitive
+    // for a matrix-free approach. However, these will be computed
+    // regardless due to MFEM's implementation of these operators.
+    // Further since these do not change it is not a performance issue.
+
+    // We need to store these separately because the jacobian and preconditioner
+    // must be computed from the boundary aware operators (which will be these
+    // matrices) whereas the residuals must be computed from the raw, physical,
+    // operators
+    std::unique_ptr<mfem::SparseMatrix> m_Mmat;
+    std::unique_ptr<mfem::SparseMatrix> m_Qmat;
+    std::unique_ptr<mfem::SparseMatrix> m_Dmat;
+
     const mfem::Array<int> m_blockOffsets;
 
     SSE::MFEMArrayPair m_theta_ess_tdofs;
     SSE::MFEMArrayPair m_phi_ess_tdofs;
 
-    // std::unique_ptr<mfem::SparseMatrix> m_Mmat, m_Qmat, m_Dmat;
-    std::unique_ptr<mfem::ScaledOperator> m_negM_op;
-    std::unique_ptr<mfem::ScaledOperator> m_negQ_op;
+    std::unique_ptr<mfem::ScaledOperator> m_negM_mat;
+    std::unique_ptr<mfem::ScaledOperator> m_negQ_mat;
     mutable std::unique_ptr<mfem::BlockOperator> m_jacobian;
 
-    // mutable std::unique_ptr<mfem::SparseMatrix> m_invSchurCompliment;
-    // mutable std::unique_ptr<mfem::SparseMatrix> m_invNonlinearJacobian;
     mutable std::unique_ptr<SchurCompliment> m_schurCompliment;
     mutable std::unique_ptr<GMRESInverter> m_invSchurCompliment;
     mutable std::unique_ptr<mfem::Solver> m_invNonlinearJacobian;
 
-    /*
-     * The schur preconditioner has the form
-     *
-     *  ⎡ḟ(θ)^-1 0 ⎤
-     * ⎣ 0     S^-1 ⎦
-     *
-     * Where S is the Schur compliment of the system
-     *
-     */
-
     mutable std::unique_ptr<mfem::BlockDiagonalPreconditioner> m_schurPreconditioner;
 
     bool m_isFinalized = false;