From 0ec1b6e7514fa4f004c6469724d62f88f736a52a Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 19 Mar 2025 11:15:37 -0400 Subject: [PATCH] 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 | 20 ++++++++ src/poly/utils/private/polyMFEMUtils.cpp | 65 +++++++++++++++++++++++- src/poly/utils/public/polyMFEMUtils.h | 8 +-- 3 files changed, 87 insertions(+), 6 deletions(-) diff --git a/meson.build b/meson.build index cdf0c15..e54a0c2 100644 --- a/meson.build +++ b/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 diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp index 5314beb..438e72b 100644 --- a/src/poly/utils/private/polyMFEMUtils.cpp +++ b/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 #include +#include +#include + +#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 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]); + } + } + double minDistanceToOrigin = std::numeric_limits::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) { diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h index b38a342..ce155e0 100644 --- a/src/poly/utils/public/polyMFEMUtils.h +++ b/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 -#include #include #include "config.h" #include "probe.h" -#include "quill/LogMacros.h" @@ -175,8 +173,10 @@ namespace polyMFEMUtils { mfem::Array m_originElementIDs; mfem::Array m_originIntegrationPoints; mfem::DenseMatrix m_originCoordinateMatrix; + mfem::Mesh* m_mesh; + std::vector 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;