diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp
index 94e0881..5f9576e 100644
--- a/src/poly/solver/private/polySolver.cpp
+++ b/src/poly/solver/private/polySolver.cpp
@@ -70,37 +70,40 @@ namespace laneEmden {
 }
 
 
-PolySolver::PolySolver(const double n, const double order) {
-
-    // --- Check the polytropic index ---
-    if (n > 4.99 || n < 0.0) {
-        LOG_ERROR(m_logger, "The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is {}", m_polytropicIndex);
-        throw std::runtime_error("The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is " + std::to_string(m_polytropicIndex));
-    }
-
-    m_polytropicIndex = n;
-    m_feOrder = order;
-
-    // Load and rescale the mesh
-    const ResourceManager& rm = ResourceManager::getInstance();
-    const Resource& genericResource = rm.getResource("mesh:polySphere");
-    const auto &meshIO = std::get<std::unique_ptr<MeshIO>>(genericResource);
-    meshIO->LinearRescale(polycoeff::x1(m_polytropicIndex));
-    m_mesh = std::make_unique<mfem::Mesh>(meshIO->GetMesh());
+PolySolver::PolySolver(mfem::Mesh& mesh, const double n, const double order)
+    : m_config(Config::getInstance()),
+      m_logManager(Probe::LogManager::getInstance()),
+      m_logger(m_logManager.getLogger("log")),
+      m_polytropicIndex(n),
+      m_feOrder(order),
+      m_mesh(mesh) {
 
     // Use feOrder - 1 for the RT space to satisfy Brezzi-Babuska condition
     // for the H1 and RT [H(div)] spaces
-    m_fecH1 = std::make_unique<mfem::H1_FECollection>(m_feOrder, m_mesh->SpaceDimension());
-    m_fecRT = std::make_unique<mfem::RT_FECollection>(m_feOrder - 1, m_mesh->SpaceDimension());
+    m_fecH1 = std::make_unique<mfem::H1_FECollection>(m_feOrder, m_mesh.SpaceDimension());
+    m_fecRT = std::make_unique<mfem::RT_FECollection>(m_feOrder - 1, m_mesh.SpaceDimension());
 
-    m_feTheta = std::make_unique<mfem::FiniteElementSpace>(m_mesh.get(), m_fecH1.get());
-    m_fePhi = std::make_unique<mfem::FiniteElementSpace>(m_mesh.get(), m_fecRT.get());
+    m_feTheta = std::make_unique<mfem::FiniteElementSpace>(&m_mesh, m_fecH1.get());
+    m_fePhi = std::make_unique<mfem::FiniteElementSpace>(&m_mesh, m_fecRT.get());
 
     m_theta = std::make_unique<mfem::GridFunction>(m_feTheta.get());
     m_phi = std::make_unique<mfem::GridFunction>(m_fePhi.get());
 
     assembleBlockSystem();
+}
 
+PolySolver::PolySolver(const double n, const double order)
+    : PolySolver(prepareMesh(n), n, order){}
+
+mfem::Mesh& PolySolver::prepareMesh(const double n) {
+    if (n > 4.99 || n < 0.0) {
+        throw std::runtime_error("The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is " + std::to_string(n));
+    }
+    const ResourceManager& rm = ResourceManager::getInstance();
+    const Resource& genericResource = rm.getResource("mesh:polySphere");
+    const auto &meshIO = std::get<std::unique_ptr<MeshIO>>(genericResource);
+    meshIO->LinearRescale(polycoeff::x1(n));
+    return meshIO->GetMesh();
 }
 
 PolySolver::~PolySolver() = default;
@@ -139,7 +142,7 @@ std::unique_ptr<formBundle> PolySolver::buildIndividualForms(const mfem::Array<i
     );
 
     // --- -M negation -> M ---
-    mfem::Vector negOneVec(m_mesh->SpaceDimension());
+    mfem::Vector negOneVec(m_mesh.SpaceDimension());
     negOneVec = -1.0;
 
     m_negationCoeff = std::make_unique<mfem::VectorConstantCoefficient>(negOneVec);
@@ -214,7 +217,7 @@ SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const {
 
     thetaCenterVals = 1.0;
 
-    mfem::Array<int> ess_brd(m_mesh->bdr_attributes.Max());
+    mfem::Array<int> ess_brd(m_mesh.bdr_attributes.Max());
     ess_brd = 1;
 
     mfem::Array<double> thetaSurfaceVals, phiSurfaceVals;
@@ -240,14 +243,14 @@ SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const {
 std::pair<mfem::Array<int>, mfem::Array<int>> PolySolver::findCenterElement() const {
     mfem::Array<int> thetaCenterDofs;
     mfem::Array<int> phiCenterDofs;
-    mfem::DenseMatrix centerPoint(m_mesh->SpaceDimension(), 1);
+    mfem::DenseMatrix centerPoint(m_mesh.SpaceDimension(), 1);
     centerPoint(0, 0) = 0.0;
     centerPoint(1, 0) = 0.0;
     centerPoint(2, 0) = 0.0;
 
     mfem::Array<int> elementIDs;
     mfem::Array<mfem::IntegrationPoint> ips;
-    m_mesh->FindPoints(centerPoint, elementIDs, ips);
+    m_mesh.FindPoints(centerPoint, elementIDs, ips);
     mfem::Array<int> tempDofs;
     for (int i = 0; i < elementIDs.Size(); i++) {
         m_feTheta->GetElementDofs(elementIDs[i], tempDofs);
@@ -263,7 +266,7 @@ void PolySolver::setInitialGuess() const {
     mfem::FunctionCoefficient thetaInitGuess (
         [this](const mfem::Vector &x) {
             const double r = x.Norml2();
-            const double radius = Probe::getMeshRadius(*m_mesh);
+            // const double radius = Probe::getMeshRadius(*m_mesh);
             // const double u = 1/radius;
 
             // return (-1.0/radius) * r + 1;
@@ -272,18 +275,18 @@ void PolySolver::setInitialGuess() const {
         }
     );
 
-    mfem::VectorFunctionCoefficient phiSurfaceVectors (m_mesh->SpaceDimension(),
+    mfem::VectorFunctionCoefficient phiSurfaceVectors (m_mesh.SpaceDimension(),
         [this](const mfem::Vector &x, mfem::Vector &y) {
             const double r = x.Norml2();
             mfem::Vector xh(x);
             xh /= r; // Normalize the vector
-            y.SetSize(m_mesh->SpaceDimension());
+            y.SetSize(m_mesh.SpaceDimension());
             y = xh;
             y *= polycoeff::thetaSurfaceFlux(m_polytropicIndex);
         }
     );
     // We want to apply specific boundary conditions to the surface
-    mfem::Array<int> ess_brd(m_mesh->bdr_attributes.Max());
+    mfem::Array<int> ess_brd(m_mesh.bdr_attributes.Max());
     ess_brd = 1;
 
     // θ = 0 at surface
@@ -294,11 +297,19 @@ void PolySolver::setInitialGuess() const {
 
     mfem::GradientGridFunctionCoefficient phiInitGuess (m_theta.get());
     m_phi->ProjectCoefficient(phiInitGuess);
+
+    // Note that this will not result in perfect boundary conditions
+    // because it must maintain H(div) continuity, this is
+    // why inhomogenous boundary conditions enforcement is needed for φ
+    // This manifests in PolytropeOperator::Mult where we do not
+    // just zero out the essential dof elements in the residuals vector
+    // for φ; rather, we need to set this to something which will push the
+    // solver towards a more consistent answer (x_φ - target)
     m_phi->ProjectBdrCoefficientNormal(phiSurfaceVectors, ess_brd);
 
     if (m_config.get<bool>("Poly:Solver:ViewInitialGuess", false)) {
-        Probe::glVisView(*m_theta, *m_mesh, "θ init");
-        Probe::glVisView(*m_phi, *m_mesh, "φ init");
+        Probe::glVisView(*m_theta, m_mesh, "θ init");
+        Probe::glVisView(*m_phi, m_mesh, "φ init");
     }
 
 }
diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h
index b463973..6448c9a 100644
--- a/src/poly/solver/public/polySolver.h
+++ b/src/poly/solver/public/polySolver.h
@@ -28,6 +28,7 @@
 #include "4DSTARTypes.h"
 #include "operator.h"
 #include "config.h"
+#include "meshIO.h"
 #include "probe.h"
 #include "quill/Logger.h"
 
@@ -60,16 +61,16 @@ public: // Public methods
 
     double getN() const { return m_polytropicIndex; }
     double getOrder() const { return m_feOrder; }
-    mfem::Mesh* getMesh() const { return m_mesh.get(); }
+    mfem::Mesh& getMesh() const { return m_mesh; }
     mfem::GridFunction& getSolution() const { return *m_theta; }
 
 private: // Private Attributes
-    Config& m_config = Config::getInstance();
-    Probe::LogManager& m_logManager = Probe::LogManager::getInstance();
-    quill::Logger* m_logger = m_logManager.getLogger("log");
+    Config& m_config;
+    Probe::LogManager& m_logManager;
+    quill::Logger* m_logger;
     double m_polytropicIndex, m_feOrder;
-    
-    std::unique_ptr<mfem::Mesh> m_mesh;
+
+    mfem::Mesh& m_mesh;
     std::unique_ptr<mfem::H1_FECollection> m_fecH1;
     std::unique_ptr<mfem::RT_FECollection> m_fecRT;
 
@@ -87,6 +88,11 @@ private: // Private Attributes
 
 
 private: // Private methods
+    PolySolver(mfem::Mesh& mesh, double n, double order);
+
+    static mfem::Mesh& prepareMesh(double n);
+
+
     void assembleBlockSystem();
 
     /**