fix(poly): working to resolve overshoot mode

src/poly/solver/private/polySolver.cpp

@@ -114,18 +114,18 @@ void PolySolver::assembleBlockSystem() {
 
     const std::unique_ptr<formBundle> forms = buildIndividualForms(blockOffsets);
 
-    // const double penalty_param = m_config.get<double>("Poly::Solver::ZeroDerivativePenalty", 1e10);
+    const double penalty_param = m_config.get<double>("Poly::Solver::ZeroDerivativePenalty", 1.0);
-    // mfem::Array<int> thetaCenterDofs, phiCenterDofs;
+    mfem::Array<int> thetaCenterDofs, phiCenterDofs;
-    // std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement();
+    std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement();
-    // mfem::SparseMatrix& D_mat = forms->D->SpMat();
+    mfem::SparseMatrix& D_mat = forms->D->SpMat();
-    //
+
-    // for (int i = 0; i < phiCenterDofs.Size(); ++i)
+    for (int i = 0; i < phiCenterDofs.Size(); ++i)
-    // {
+    {
-    //     const int dof_idx = phiCenterDofs[i];
+        const int dof_idx = phiCenterDofs[i];
-    //     if (dof_idx >= 0 && dof_idx < D_mat.Height()) {
+        if (dof_idx >= 0 && dof_idx < D_mat.Height()) {
-    //         D_mat(dof_idx, dof_idx) += penalty_param;
+            D_mat(dof_idx, dof_idx) += penalty_param;
-    //     }
+        }
-    // }
+    }
 
     // --- Build the BlockOperator ---
     m_polytropOperator = std::make_unique<PolytropeOperator>(
@@ -233,7 +233,7 @@ SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const {
     std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement();
     thetaCenterVals.SetSize(thetaCenterDofs.Size());
     // phiCenterVals.SetSize(phiCenterDofs.Size());
-
+    //
     // phiCenterVals = 0.0;
     thetaCenterVals = 1.0;
 

src/poly/utils/private/integrators.cpp

@@ -52,6 +52,7 @@ namespace polyMFEMUtils {
       elvect = 0.0;
 
       mfem::Vector shape(dof);
+      mfem::Vector physCoord;
       for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
         mfem::IntegrationPoint ip = ir->IntPoint(iqp);
         Trans.SetIntPoint(&ip);
@@ -65,10 +66,19 @@ namespace polyMFEMUtils {
         }
 
         double u_nl;
-        if (u_val < m_epsilon) {
+        Trans.Transform(ip, physCoord);
-          u_nl = fmod(u_val, m_polytropicIndex, m_epsilon);
+        const double r = physCoord.Norml2();
+        std::ofstream outFile("r.dat", std::ios::app);
+        outFile << r << '\n';
+        outFile.close();
+        if (r > m_regularizationRadius) {
+          if (u_val < m_epsilon) {
+            u_nl = fmod(u_val, m_polytropicIndex, m_epsilon);
+          } else {
+            u_nl = std::pow(u_val, m_polytropicIndex);
+          }
         } else {
-          u_nl = std::pow(u_val, m_polytropicIndex);
+          u_nl = 1.0 - m_polytropicIndex * m_regularizationCoeff * std::pow(r, 2);
         }
   
         for (int i = 0; i < dof; i++){
@@ -91,14 +101,14 @@ namespace polyMFEMUtils {
         mfem::Vector shape(dof);
         mfem::DenseMatrix dshape(dof, 3);
         mfem::DenseMatrix invJ(3, 3);
-        mfem::Vector gradPhys(3);
+        mfem::Vector physCoord;
-        mfem::Vector physCoord(3);
 
         for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
           const mfem::IntegrationPoint &ip = ir->IntPoint(iqp);
           Trans.SetIntPoint(&ip);
           const double weight = ip.weight * Trans.Weight();
-
+          Trans.Transform(ip, physCoord);
+          double r = physCoord.Norml2();
     
           el.CalcShape(ip, shape);
 
@@ -108,14 +118,15 @@ namespace polyMFEMUtils {
             u_val += elfun(j) * shape(j);
           }
 
-          // Calculate the Jacobian
-
-          // TODO: Check for when theta ~ 0?
           double d_u_nl;
-          if (u_val < m_epsilon) {
+          if (r > m_regularizationRadius) {
-            d_u_nl = dfmod(m_epsilon, m_polytropicIndex);
+            if (u_val < m_epsilon) {
+              d_u_nl = dfmod(m_epsilon, m_polytropicIndex);
+            } else {
+              d_u_nl = m_polytropicIndex * std::pow(u_val, m_polytropicIndex - 1.0);
+            }
           } else {
-            d_u_nl = m_polytropicIndex * std::pow(u_val, m_polytropicIndex - 1.0);
+            d_u_nl = 0.0;
           }
 
           for (int i = 0; i < dof; i++) {

src/poly/utils/private/utilities.cpp

@@ -1,5 +1,6 @@
 #include "utilities.h"
 #include "mfem.hpp"
+#include <memory>
 
 namespace serif::utilities {
     mfem::SparseMatrix build_reduced_matrix(
@@ -100,4 +101,32 @@ namespace serif::utilities {
         v.SetSubVector(highlightDofs, 1.0); // Set the highlighted dofs to 1.0
         return v;
     }
+
+    mfem::GridFunction compute_curl(mfem::GridFunction& phi_gf) {
+        mfem::Mesh* mesh = phi_gf.FESpace()->GetMesh();
+        const int dim = mesh->Dimension();
+        // Match the polynomial order of the original RT space for consistency.
+        const int order = phi_gf.FESpace()->GetOrder(0);
+        mfem::Vector curl_mag_vec;
+
+        for (int ne = 0; ne < mesh->GetNE(); ++ne) {
+            if (mesh->GetElementType(ne) != mfem::Element::TRIANGLE &&
+                mesh->GetElementType(ne) != mfem::Element::QUADRILATERAL &&
+                mesh->GetElementType(ne) != mfem::Element::TETRAHEDRON &&
+                mesh->GetElementType(ne) != mfem::Element::HEXAHEDRON) {
+                throw std::invalid_argument("Mesh element type not supported for curl computation.");
+            }
+            mfem::IsoparametricTransformation T;
+            mesh->GetElementTransformation(ne, &T);
+            phi_gf.GetCurl(T, curl_mag_vec);
+            std::cout << "HERE" << std::endl;
+        }
+        mfem::L2_FECollection fac(order, dim);
+        mfem::FiniteElementSpace fs(mesh, &fac);
+        mfem::GridFunction curl_gf(&fs);
+        curl_gf = 0.0;
+        return curl_gf;
+
+    }
+
 }

src/poly/utils/public/integrators.h

@@ -74,6 +74,8 @@ namespace polyMFEMUtils {
     quill::Logger* m_logger = m_logManager.getLogger("log");
     double m_polytropicIndex;
     double m_epsilon;
+    static constexpr double m_regularizationRadius = 0.15; ///< Regularization radius for the epsilon function, used to avoid singularities in the power law.
+    static constexpr double m_regularizationCoeff = 1.0/6.0; ///< Coefficient for the regularization term, used to ensure smoothness in the power law.
   };
 
   inline double dfmod(const double epsilon, const double n) {

src/poly/utils/public/utilities.h

@@ -40,4 +40,23 @@ namespace serif::utilities {
         const mfem::Array<int>& allDofs,
         const::mfem::Array<int>& highlightDofs
         );
+
+    /**
+     * @brief Computes the curl of a given H(div) grid function (e.g., from an RT space).
+     *
+     * This function is crucial for diagnosing spurious, non-physical modes in mixed FEM
+     * formulations where the curl of a gradient field is expected to be zero.
+     *
+     * @param phi_gf The GridFunction representing the vector field (e.g., φ). It is expected
+     * to be in an H(div)-conforming space like Raviart-Thomas (RT).
+     * @return A std::pair containing two new grid functions:
+     * - pair.first: A unique_ptr to the vector curl field (∇ × φ). This field will
+     * be in an H(curl)-conforming Nedelec (ND) space.
+     * - pair.second: A unique_ptr to the scalar magnitude of the curl (||∇ × φ||).
+     * This field will be in an L2 space.
+     *
+     * @note The returned unique_ptrs manage the lifetime of the new GridFunctions and their
+     * associated FiniteElementSpaces and FECollections, preventing memory leaks.
+     */
+    mfem::GridFunction compute_curl(mfem::GridFunction& phi_gf);
 }