feat(poly): find all connected elements to central vertex

In order to constrain the central slope we find all the elements connected to the central vertex. The slope will be approximated over these using the finite difference method

meson.build

@@ -1,3 +1,23 @@
+# ***********************************************************************
+#
+#   Copyright (C) 2025 -- The 4D-STAR Collaboration
+#   File Author: Emily Boudreaux
+#   Last Modified: March 19, 2025
+#
+#   4DSSE is free software; you can use it and/or modify
+#   it under the terms and restrictions the GNU General Library Public
+#   License version 3 (GPLv3) as published by the Free Software Foundation.
+#
+#   4DSSE is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+#   See the GNU Library General Public License for more details.
+#
+#   You should have received a copy of the GNU Library General Public License
+#   along with this software; if not, write to the Free Software
+#   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# *********************************************************************** #
 project('4DSSE', 'cpp', version: '0.0.1a', default_options: ['cpp_std=c++23'], meson_version: '>=1.6.0')
 
 # Add default visibility for all C++ targets

src/poly/utils/private/polyMFEMUtils.cpp

@@ -2,7 +2,7 @@
 //
 //   Copyright (C) 2025 -- The 4D-STAR Collaboration
 //   File Author: Emily Boudreaux
-//   Last Modified: March 18, 2025
+//   Last Modified: March 19, 2025
 //
 //   4DSSE is free software; you can use it and/or modify
 //   it under the terms and restrictions the GNU General Library Public
@@ -21,6 +21,10 @@
 #include "mfem.hpp"
 #include <cmath>
 #include <vector>
+#include <limits>
+#include <stdexcept>
+
+#include "quill/LogMacros.h"
 
 #include "polyMFEMUtils.h"
 #include "probe.h"
@@ -191,12 +195,69 @@ namespace polyMFEMUtils {
     }
   }
 
-  ConstraintIntegrator::ConstraintIntegrator(const double gamma): m_gamma(gamma) {
+  // TODO: break this up into smaller functions
+  // TODO: think of a more efficient way to find connected elements
+  ConstraintIntegrator::ConstraintIntegrator(const double gamma, mfem::Mesh* mesh): m_gamma(gamma), m_mesh(mesh) {
     LOG_INFO(m_logger, "Initializing Constraint Integrator...");
     m_originCoordinateMatrix.SetSize(3, 1);
     m_originCoordinateMatrix(0, 0) = 0.0;
     m_originCoordinateMatrix(1, 0) = 0.0;
     m_originCoordinateMatrix(2, 0) = 0.0;
+
+    m_mesh->FindPoints(m_originCoordinateMatrix, m_originElementIDs, m_originIntegrationPoints);
+
+    if (m_originElementIDs.Size() == 0) {
+      LOG_ERROR(m_logger, "The origin point is not found in the mesh.");
+      throw std::runtime_error("The origin point is not found in the mesh.");
+    }
+
+    LOG_INFO(m_logger, "The origin point is found in the mesh.");
+
+    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<int> originVertexIds;
+    mfem::Array<int> row;
+    for (int i = 0; i < nVertices; i++) {
+      VETable->GetRow(i, row);
+      if (row[1] == m_originElementIDs[0]) {
+        originVertexIds.push_back(row[0]);
+      }
+    }
+    double minDistanceToOrigin = std::numeric_limits<double>::max();
+    int minDistanceVertexId = -1;
+    for (const auto &vertexId : originVertexIds) {
+      mfem::Vector vertex;
+      const double* vcoord = m_mesh->GetVertex(vertexId);
+      // --- Note if this is run with a 2D or 1D mesh this may lead to a segfault ---
+      // TODO: Add a check for the dimension of the mesh
+      double distance = vcoord[0]*vcoord[0] + vcoord[1]*vcoord[1] + vcoord[2]*vcoord[2];
+      if (distance < minDistanceToOrigin) {
+        minDistanceToOrigin = distance;
+        minDistanceVertexId = vertexId;
+      }
+
+
+    }
+
+    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
+
+    // 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());
+
   }
 
   void ConstraintIntegrator::AssembleElementMatrix(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, mfem::DenseMatrix &elmat) {

src/poly/utils/public/polyMFEMUtils.h

@@ -2,7 +2,7 @@
 //
 //   Copyright (C) 2025 -- The 4D-STAR Collaboration
 //   File Author: Emily Boudreaux
-//   Last Modified: February 14, 2025
+//   Last Modified: March 19, 2025
 //
 //   4DSSE is free software; you can use it and/or modify
 //   it under the terms and restrictions the GNU General Library Public
@@ -23,11 +23,9 @@
 
 #include "mfem.hpp"
 #include <string>
-#include <array>
 #include <vector>
 #include "config.h"
 #include "probe.h"
-#include "quill/LogMacros.h"
 
 
 
@@ -175,8 +173,10 @@ namespace polyMFEMUtils {
     mfem::Array<int> m_originElementIDs;
     mfem::Array<mfem::IntegrationPoint> m_originIntegrationPoints;
     mfem::DenseMatrix m_originCoordinateMatrix;
+    mfem::Mesh* m_mesh;
+    std::vector<int> m_originConnectedElementIds;
   public:
-    ConstraintIntegrator(const double gamma);
+    ConstraintIntegrator(const double gamma, mfem::Mesh* mesh);
     ~ConstraintIntegrator() = default;
 
     void AssembleElementMatrix(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, mfem::DenseMatrix &elmat) override;