diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 50c6cc9..767618c 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -114,18 +114,18 @@ void PolySolver::assembleBlockSystem() { const std::unique_ptr forms = buildIndividualForms(blockOffsets); - // const double penalty_param = m_config.get("Poly::Solver::ZeroDerivativePenalty", 1e10); - // mfem::Array thetaCenterDofs, phiCenterDofs; - // std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement(); - // mfem::SparseMatrix& D_mat = forms->D->SpMat(); - // - // for (int i = 0; i < phiCenterDofs.Size(); ++i) - // { - // const int dof_idx = phiCenterDofs[i]; - // if (dof_idx >= 0 && dof_idx < D_mat.Height()) { - // D_mat(dof_idx, dof_idx) += penalty_param; - // } - // } + const double penalty_param = m_config.get("Poly::Solver::ZeroDerivativePenalty", 1.0); + mfem::Array thetaCenterDofs, phiCenterDofs; + std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement(); + mfem::SparseMatrix& D_mat = forms->D->SpMat(); + + for (int i = 0; i < phiCenterDofs.Size(); ++i) + { + const int dof_idx = phiCenterDofs[i]; + if (dof_idx >= 0 && dof_idx < D_mat.Height()) { + D_mat(dof_idx, dof_idx) += penalty_param; + } + } // --- Build the BlockOperator --- m_polytropOperator = std::make_unique( @@ -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; diff --git a/src/poly/utils/private/integrators.cpp b/src/poly/utils/private/integrators.cpp index 056bc0c..b2f122e 100644 --- a/src/poly/utils/private/integrators.cpp +++ b/src/poly/utils/private/integrators.cpp @@ -50,8 +50,9 @@ namespace polyMFEMUtils { int dof = el.GetDof(); elvect.SetSize(dof); 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) { - u_nl = fmod(u_val, m_polytropicIndex, m_epsilon); + Trans.Transform(ip, physCoord); + 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(3); + mfem::Vector physCoord; 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) { - d_u_nl = dfmod(m_epsilon, m_polytropicIndex); + if (r > m_regularizationRadius) { + 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++) { diff --git a/src/poly/utils/private/utilities.cpp b/src/poly/utils/private/utilities.cpp index f0ffe08..74ae125 100644 --- a/src/poly/utils/private/utilities.cpp +++ b/src/poly/utils/private/utilities.cpp @@ -1,5 +1,6 @@ #include "utilities.h" #include "mfem.hpp" +#include 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; + + } + } diff --git a/src/poly/utils/public/integrators.h b/src/poly/utils/public/integrators.h index fe34bf2..9437665 100644 --- a/src/poly/utils/public/integrators.h +++ b/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) { diff --git a/src/poly/utils/public/utilities.h b/src/poly/utils/public/utilities.h index a1a1e42..a20cac3 100644 --- a/src/poly/utils/public/utilities.h +++ b/src/poly/utils/public/utilities.h @@ -40,4 +40,23 @@ namespace serif::utilities { const mfem::Array& allDofs, const::mfem::Array& 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); }