fix(poly): have a much closer to working polytrope, still small divergence near core

src/poly/solver/private/polySolver.cpp

@@ -63,9 +63,6 @@ namespace laneEmden {
     }
 }
 
-// TODO: Come back to this and think of a better way to get the mesh file
-const std::string SPHERICAL_MESH = std::string(getenv("MESON_SOURCE_ROOT")) + "/src/resources/mesh/core.msh";
-
 PolySolver::PolySolver(double n, double order, mfem::Mesh& mesh_)
 : logger(logManager.getLogger("log")),
   n(n),

src/poly/utils/private/polyMFEMUtils.cpp

@@ -23,6 +23,8 @@
 #include <vector>
 #include <limits>
 #include <stdexcept>
+#include <set>
+#include <unordered_map>
 
 #include "quill/LogMacros.h"
 
@@ -264,15 +266,18 @@ namespace polyMFEMUtils {
 
   }
 
+
   void ConstraintIntegrator::AssembleElementMatrix(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, mfem::DenseMatrix &elmat) {
     int elemID = Trans.ElementNo;
 
     // -- Check if the element is connected to the origin vertex --
-    int checkID = m_originConnectedElementIds.Find(elemID);
-    if (checkID == -1) {
+    bool isConnected = m_originConnectedElementIds.Find(elemID) != -1 ? true : false;
+
+    if (!isConnected) {
       elmat = 0.0;
       return;
     }
+    
 
     // -- Compute the derivitives using MFEM's build in shape routines --
 
@@ -312,7 +317,7 @@ namespace polyMFEMUtils {
         for (int dimID = 0; dimID < dim; dimID++) {
           dotProduct += dshapePhysical(i, dimID) * dshapePhysical(j, dimID);
         }
-        elmat(i, j) += m_gamma * dotProduct;
+        elmat(i, j) +=  m_gamma * dotProduct;
       }
     }
   }

src/poly/utils/public/polyMFEMUtils.h

@@ -26,6 +26,7 @@
 #include <vector>
 #include "config.h"
 #include "probe.h"
+#include <unordered_map>