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;