diff --git a/src/poly/solver/meson.build b/src/poly/solver/meson.build index 62e5448..641151b 100644 --- a/src/poly/solver/meson.build +++ b/src/poly/solver/meson.build @@ -10,7 +10,7 @@ libPolySolver = static_library('polySolver', polySolver_sources, include_directories : include_directories('./public'), cpp_args: ['-fvisibility=default'], - dependencies: [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep, warning_control_dep, probe_dep, quill_dep, config_dep], + dependencies: [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep, macros_dep, probe_dep, quill_dep, config_dep], install: true ) @@ -18,5 +18,5 @@ polysolver_dep = declare_dependency( include_directories : include_directories('./public'), link_with : libPolySolver, sources : polySolver_sources, - dependencies : [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep, warning_control_dep, probe_dep, quill_dep, config_dep] + dependencies : [mfem_dep, meshio_dep, polycoeff_dep, polyutils_dep, macros_dep, probe_dep, quill_dep, config_dep] ) \ No newline at end of file diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 75446ef..c5f87b9 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -20,6 +20,7 @@ // *********************************************************************** */ #include "mfem.hpp" +#include #include #include @@ -102,6 +103,13 @@ void PolySolver::assembleNonlinearForm() { auto nonlinearIntegrator = std::make_unique(*nonlinearSourceCoeff, n); compositeIntegrator->add_integrator(nonlinearIntegrator.release()); + // Add the contraint term \gamma(\nabla \theta(0)\cdot\nabla v(0))^{2} + double gamma = config.get("Poly:Solver:Constraint:Gamma", 1e2); + auto constraintIntegrator = std::make_unique( + new polyMFEMUtils::ConstraintIntegrator(gamma, &mesh) + ); + compositeIntegrator->add_integrator(constraintIntegrator.release()); + nonlinearForm->AddDomainIntegrator(compositeIntegrator.release()); } @@ -110,7 +118,7 @@ void PolySolver::solve(){ mfem::FunctionCoefficient initCoeff ( [this](const mfem::Vector &x) { double r = x.Norml2(); - double theta = laneEmden::thetaSerieseExpansion(r, n, 12); + double theta = laneEmden::thetaSerieseExpansion(r, n, 10); return theta; // double radius = Probe::getMeshRadius(mesh); // double u = 1/radius; diff --git a/src/poly/utils/meson.build b/src/poly/utils/meson.build index adb5d1c..f7096fc 100644 --- a/src/poly/utils/meson.build +++ b/src/poly/utils/meson.build @@ -12,7 +12,7 @@ libpolyutils = static_library('polyutils', polyutils_sources, include_directories : include_directories('./public'), cpp_args: ['-fvisibility=default'], - dependencies: [mfem_dep, warning_control_dep, probe_dep, quill_dep, config_dep], + dependencies: [mfem_dep, macros_dep, probe_dep, quill_dep, config_dep], install: true ) @@ -20,5 +20,5 @@ polyutils_dep = declare_dependency( include_directories : include_directories('./public'), link_with : libpolyutils, sources : polyutils_sources, - dependencies : [mfem_dep, warning_control_dep, probe_dep, quill_dep, config_dep] + dependencies : [mfem_dep, macros_dep, probe_dep, quill_dep, config_dep] ) diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp index 438e72b..6c710b0 100644 --- a/src/poly/utils/private/polyMFEMUtils.cpp +++ b/src/poly/utils/private/polyMFEMUtils.cpp @@ -27,7 +27,8 @@ #include "quill/LogMacros.h" #include "polyMFEMUtils.h" -#include "probe.h" +#include "debug.h" + namespace polyMFEMUtils { NonlinearPowerIntegrator::NonlinearPowerIntegrator( @@ -204,6 +205,9 @@ namespace polyMFEMUtils { m_originCoordinateMatrix(1, 0) = 0.0; m_originCoordinateMatrix(2, 0) = 0.0; + m_originCoordinateVector.SetSize(3); + m_originCoordinateVector = 0.0; + m_mesh->FindPoints(m_originCoordinateMatrix, m_originElementIDs, m_originIntegrationPoints); if (m_originElementIDs.Size() == 0) { @@ -213,19 +217,29 @@ namespace polyMFEMUtils { LOG_INFO(m_logger, "The origin point is found in the mesh."); + // NOTE (EMB, March 2025): This function as it is currently written will break if the mesh is refined after being passed to the constructor + // This may or may not be an issue (it does seem unlikley that the mesh would be refined after being passed to the constructor) + // But if something mysteriously breaks in the future this is may be a good place to start looking mfem::Table* VETable = m_mesh->GetVertexToElementTable(); const int nVertices = VETable->Size(); LOG_INFO(m_logger, "The number of vertices in the mesh is {}", nVertices); std::vector originVertexIds; - mfem::Array row; - for (int i = 0; i < nVertices; i++) { - VETable->GetRow(i, row); - if (row[1] == m_originElementIDs[0]) { - originVertexIds.push_back(row[0]); + mfem::Array connectedElements; + + // -- Get all vertices connected to the origin element -- + for (int vertexID = 0; vertexID < nVertices; vertexID++) { + VETable->GetRow(vertexID, connectedElements); + for (int j = 0; j < connectedElements.Size(); j++) { + if (connectedElements[j] == m_originElementIDs[0]) { + originVertexIds.push_back(vertexID); + } } + } double minDistanceToOrigin = std::numeric_limits::max(); int minDistanceVertexId = -1; + + // -- Get the vertex closest to the origin ID -- for (const auto &vertexId : originVertexIds) { mfem::Vector vertex; const double* vcoord = m_mesh->GetVertex(vertexId); @@ -239,29 +253,68 @@ namespace polyMFEMUtils { } - if (minDistanceVertexId == -1 || minDistanceToOrigin > 1e-10) { LOG_ERROR(m_logger, "The origin vertex is not found in the mesh."); throw std::runtime_error("The origin vertex is not found in the mesh."); } // -- Find all elements connected to the origin vertex by looping through the VE table + VETable->GetRow(minDistanceVertexId, m_originConnectedElementIds); - // NOTE (EMB, March 2025): This function as it is currently written will break if the mesh is refined after being passed to the constructor - // This may or may not be an issue (it does seem unlikley that the mesh would be refined after being passed to the constructor) - // But if something mysteriously breaks in the future this is may be a good place to start looking - for (int i = 0; i < nVertices; i++) { - VETable->GetRow(i, row); - if (row[0] == minDistanceVertexId) { - m_originConnectedElementIds.push_back(row[1]); - } - } - - LOG_INFO(m_logger, "Found {} elements connected to the origin vertex.", m_originConnectedElementIds.size()); + LOG_INFO(m_logger, "Found {} elements connected to the origin vertex.", m_originConnectedElementIds.Size()); } 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) { + elmat = 0.0; + return; + } + + // -- Compute the derivitives using MFEM's build in shape routines -- + + int numDoF = el.GetDof(); + int dim = m_mesh->Dimension(); + + // -- Map the origin in physical space to the reference space of the element -- + mfem::Vector originReferenceCoordinate(dim); + mfem::IntegrationPoint originIntegrationPoint; + Trans.TransformBack(m_originCoordinateVector, originIntegrationPoint); + + // -- Compute the derivitives of the shape functions at the origin -- + mfem::DenseMatrix dshape(numDoF, dim); + el.CalcDShape(originIntegrationPoint, dshape); + + // -- Transform derivitives from reference space to physical space using the inverse of the Jacobian -- + mfem::DenseMatrix invJac(dim, dim); + invJac = Trans.InverseJacobian(); + + mfem::DenseMatrix dshapePhysical(numDoF, dim); + dshapePhysical = 0.0; + + for (int dofID = 0; dofID < numDoF; dofID++) { + for (int dimID = 0; dimID < dim; dimID++) { + for (int i = 0; i < dim; i++) { + dshapePhysical(dofID, dimID) += dshape(dofID, i) * invJac(i, dimID); + } + } + } + + // -- Assemble the element matrix contribution = gamma * (grad(phi_i) dot grad(phi_j)) -- + elmat.SetSize(numDoF); + elmat = 0.0; + for (int i = 0; i < numDoF; i++) { + for (int j = 0; j < numDoF; j++) { + double dotProduct = 0.0; + for (int dimID = 0; dimID < dim; dimID++) { + dotProduct += dshapePhysical(i, dimID) * dshapePhysical(j, dimID); + } + elmat(i, j) += m_gamma * dotProduct; + } + } } diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h index ce155e0..ca176b9 100644 --- a/src/poly/utils/public/polyMFEMUtils.h +++ b/src/poly/utils/public/polyMFEMUtils.h @@ -173,8 +173,9 @@ namespace polyMFEMUtils { mfem::Array m_originElementIDs; mfem::Array m_originIntegrationPoints; mfem::DenseMatrix m_originCoordinateMatrix; + mfem::Vector m_originCoordinateVector; mfem::Mesh* m_mesh; - std::vector m_originConnectedElementIds; + mfem::Array m_originConnectedElementIds; public: ConstraintIntegrator(const double gamma, mfem::Mesh* mesh); ~ConstraintIntegrator() = default;