From e5864ca31e2e04633f5c39c5d56a314fa23aebcc Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 25 Apr 2025 11:22:13 -0400 Subject: [PATCH] feat(probe): added glvisview helped method to probe namespace --- src/probe/meson.build | 23 +++++- src/probe/private/probe.cpp | 153 +++++++++++++++++++++++++++++++++++- src/probe/public/probe.h | 26 +++++- 3 files changed, 194 insertions(+), 8 deletions(-) diff --git a/src/probe/meson.build b/src/probe/meson.build index 653c5eb..cf6accf 100644 --- a/src/probe/meson.build +++ b/src/probe/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 +# +# *********************************************************************** # # Define the library probe_sources = files( 'private/probe.cpp', @@ -10,7 +30,8 @@ probe_headers = files( dependencies = [ config_dep, mfem_dep, - quill_dep + quill_dep, + macros_dep ] # Define the liblogger library so it can be linked against by other parts of the build system diff --git a/src/probe/private/probe.cpp b/src/probe/private/probe.cpp index 8ae8f11..78d512a 100644 --- a/src/probe/private/probe.cpp +++ b/src/probe/private/probe.cpp @@ -2,7 +2,7 @@ // // Copyright (C) 2025 -- The 4D-STAR Collaboration // File Author: Emily Boudreaux -// Last Modified: February 23, 2025 +// Last Modified: March 18, 2025 // // 4DSSE is free software; you can use it and/or modify // it under the terms and restrictions the GNU General Library Public @@ -23,16 +23,24 @@ #include "quill/Logger.h" #include "quill/sinks/ConsoleSink.h" #include "quill/sinks/FileSink.h" +#include "quill/LogMacros.h" + #include #include #include #include +#include +#include +#include +#include #include "mfem.hpp" #include "config.h" #include "probe.h" +#include "warning_control.h" + namespace Probe { @@ -47,18 +55,155 @@ void wait(int seconds) { } void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh, - const std::string& windowTitle) { + const std::string& windowTitle, const std::string& keyset) { Config& config = Config::getInstance(); + quill::Logger* logger = LogManager::getInstance().getLogger("log"); + std::string usedKeyset; if (config.get("Probe:GLVis:Visualization", true)) { + LOG_INFO(logger, "Visualizing solution using GLVis..."); + LOG_INFO(logger, "Window title: {}", windowTitle); + if (keyset == "") { + usedKeyset = config.get("Probe:GLVis:DefaultKeyset", ""); + } else { + usedKeyset = keyset; + } + LOG_INFO(logger, "Keyset: {}", usedKeyset); std::string vishost = config.get("Probe:GLVis:Host", "localhost"); - int visport = config.get("Probe:GLVis:Port", 19916); // Changed default port + int visport = config.get("Probe:GLVis:Port", 19916); + std::cout << "GLVis visualization enabled. Opening GLVis window... " << std::endl; + std::cout << "Using host: " << vishost << " and port: " << visport << std::endl; mfem::socketstream sol_sock(vishost.c_str(), visport); sol_sock.precision(8); sol_sock << "solution\n" << mesh << u - << "window_title '" << windowTitle << "'\n" << std::flush; // Added title + << "window_title '" << windowTitle << + "'\n" << "keys " << usedKeyset << "\n"; + sol_sock << std::flush; } } +void glVisView(mfem::Vector &vec, mfem::FiniteElementSpace &fes, const std::string &windowTitle, const std::string& keyset) { + mfem::GridFunction gf(&fes); + + DEPRECATION_WARNING_OFF + gf.SetData(vec); + DEPRECATION_WARNING_ON + glVisView(gf, *fes.GetMesh(), windowTitle, keyset); +} + +double getMeshRadius(mfem::Mesh& mesh) { + mesh.EnsureNodes(); + const mfem::GridFunction *nodes = mesh.GetNodes(); + double *node_data = nodes->GetData(); // THIS IS KEY + + int data_size = nodes->Size(); + double max_radius = 0.0; + + for (int i = 0; i < data_size; ++i) { + if (node_data[i] > max_radius) { + max_radius = node_data[i]; + } + } + return max_radius; + +} + +std::pair, std::vector> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh, + const std::vector& rayDirection, + int numSamples, std::string filename) { + Config& config = Config::getInstance(); + Probe::LogManager& logManager = Probe::LogManager::getInstance(); + quill::Logger* logger = logManager.getLogger("log"); + LOG_INFO(logger, "Getting ray solution..."); + // Check if the directory to write to exists + // If it does not exist and MakeDir is true create it + // Otherwise throw an exception + bool makeDir = config.get("Probe:GetRaySolution:MakeDir", true); + std::filesystem::path path = filename; + + if (makeDir) { + std::filesystem::path dir = path.parent_path(); + if (!std::filesystem::exists(dir)) { + LOG_INFO(logger, "Creating directory {}", dir.string()); + std::filesystem::create_directories(dir); + } + } else { + if (!std::filesystem::exists(path.parent_path())) { + throw(std::runtime_error("Directory " + path.parent_path().string() + " does not exist")); + } + } + + std::vector samples; + samples.reserve(numSamples); + double x, y, z, r, sampleValue; + double radius = getMeshRadius(mesh); + mfem::Vector rayOrigin(3); rayOrigin = 0.0; + mfem::DenseMatrix rayPoints(3, numSamples); + std::vector radialPoints; + radialPoints.reserve(numSamples); + mfem::ElementTransformation* Trans = nullptr; + mfem::Vector physicalCoords; + + for (int i = 0; i < numSamples; i++) { + r = i * radius / numSamples; + // Let rayDirection = (theta, phi) that the ray will be cast too + x = r * std::sin(rayDirection[0]) * std::cos(rayDirection[1]); + y = r * std::sin(rayDirection[0]) * std::sin(rayDirection[1]); + z = r * std::cos(rayDirection[0]); + rayPoints(0, i) = x; + rayPoints(1, i) = y; + rayPoints(2, i) = z; + } + + mfem::Array elementIds; + mfem::Array ips; + mesh.FindPoints(rayPoints, elementIds, ips); + for (int i = 0; i < elementIds.Size(); i++) { + Trans = mesh.GetElementTransformation(elementIds[i]); + Trans->Transform(ips[i], physicalCoords); + double r = std::sqrt(physicalCoords[0] * physicalCoords[0] + + physicalCoords[1] * physicalCoords[1] + + physicalCoords[2] * physicalCoords[2]); + radialPoints.push_back(r); + + int elementId = elementIds[i]; + mfem::IntegrationPoint ip = ips[i]; + if (elementId >= 0) { // Check if the point was found in an element + sampleValue = u.GetValue(elementId, ip); + LOG_DEBUG(logger, "Probe::getRaySolution() : Ray point {} found in element {} with r={:0.2f} and theta={:0.2f}", i, elementId, r, sampleValue); + samples.push_back(sampleValue); + } else { // If the point was not found in an element + samples.push_back(0.0); + } + } + std::pair, std::vector> samplesPair(radialPoints, samples); + + if (!filename.empty()) { + std::ofstream file(filename); + if (file.is_open()) { + file << "r,u\n"; + for (int i = 0; i < numSamples; i++) { + file << samplesPair.first[i] << "," << samplesPair.second[i] << "\n"; + } + file << std::endl; + file.close(); + } else { + throw(std::runtime_error("Could not open file " + filename)); + } + } + + return samplesPair; +} + +std::pair, std::vector> getRaySolution(mfem::Vector &vec, mfem::FiniteElementSpace &fes, + const std::vector& rayDirection, + int numSamples, std::string filename) { + mfem::GridFunction gf(&fes); + DEPRECATION_WARNING_OFF + gf.SetData(vec); + DEPRECATION_WARNING_ON + return getRaySolution(gf, *fes.GetMesh(), rayDirection, numSamples, filename); +} + LogManager::LogManager() { Config& config = Config::getInstance(); quill::Backend::start(); diff --git a/src/probe/public/probe.h b/src/probe/public/probe.h index ce863d0..f0d2709 100644 --- a/src/probe/public/probe.h +++ b/src/probe/public/probe.h @@ -2,7 +2,7 @@ // // Copyright (C) 2025 -- The 4D-STAR Collaboration // File Author: Emily Boudreaux -// Last Modified: February 23, 2025 +// Last Modified: April 03, 2025 // // 4DSSE is free software; you can use it and/or modify // it under the terms and restrictions the GNU General Library Public @@ -25,7 +25,7 @@ #include #include #include -#include +#include #include "mfem.hpp" #include "quill/Logger.h" @@ -50,9 +50,29 @@ namespace Probe { * @param u The GridFunction to visualize. * @param mesh The mesh associated with the GridFunction. * @param windowTitle The title of the visualization window. + * @param keyset The keyset to use for visualization. */ void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh, - const std::string& windowTitle = "solution"); + const std::string& windowTitle = "grid function", const std::string& keyset=""); + + /** + * @brief Visualize a vector using GLVis. + * @param vec The vector to visualize. + * @param mesh The mesh associated with the vector. + * @param windowTitle The title of the visualization window. + * @param keyset The keyset to use for visualization. + */ + void glVisView(mfem::Vector &vec, mfem::FiniteElementSpace &fes, + const std::string &windowTitle = "vector", const std::string& keyset=""); + + double getMeshRadius(mfem::Mesh& mesh); + + std::pair, std::vector> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh, + const std::vector& rayDirection, int numSamples, std::string filename=""); + + std::pair, std::vector> getRaySolution(mfem::Vector &vec, mfem::FiniteElementSpace &fes, + const std::vector& rayDirection, int numSamples, std::string filename=""); + /** * @brief Class to manage logging operations.