From 77d8cc8e8606c2e3ab97d5417e426f0eebb4a821 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 3 Mar 2025 09:55:24 -0500 Subject: [PATCH] feat(probe): default glvis keysets and vector version of glvisView glVisView function now accepts a keyset to send and has an overloaded version which takes a vector and finite element space instead of just a grid function and mesh --- src/probe/private/probe.cpp | 93 ++++++++++++++++++++++++++----------- src/probe/public/probe.h | 16 +++++-- 2 files changed, 79 insertions(+), 30 deletions(-) diff --git a/src/probe/private/probe.cpp b/src/probe/private/probe.cpp index 4f38045..b5c0d0f 100644 --- a/src/probe/private/probe.cpp +++ b/src/probe/private/probe.cpp @@ -4,11 +4,14 @@ #include "quill/sinks/ConsoleSink.h" #include "quill/sinks/FileSink.h" #include "quill/LogMacros.h" + #include #include #include #include #include +#include +#include #include "mfem.hpp" @@ -31,7 +34,7 @@ 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(); if (config.get("Probe:GLVis:Visualization", true)) { std::string vishost = config.get("Probe:GLVis:Host", "localhost"); @@ -41,50 +44,56 @@ void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh, 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"; // Added title + if (!keyset.empty()) { + sol_sock << "keys " << keyset << '\n'; + } + sol_sock << std::flush; } } -void glVisView(mfem::Vector &vec, mfem::FiniteElementSpace &fes, const std::string &windowTitle) { +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); + glVisView(gf, *fes.GetMesh(), windowTitle, keyset); } double getMeshRadius(mfem::Mesh& mesh) { - int numVertices = mesh.GetNV(); // Number of vertices - double maxRadius = 0.0; + mesh.EnsureNodes(); + const mfem::GridFunction *nodes = mesh.GetNodes(); + double *node_data = nodes->GetData(); // THIS IS KEY - for (int i = 0; i < numVertices; i++) { - double* vertex; - vertex = mesh.GetVertex(i); // Get vertex coordinates + int data_size = nodes->Size(); + double max_radius = 0.0; - // Compute the Euclidean distance from the origin - double radius = std::sqrt(vertex[0] * vertex[0] + - vertex[1] * vertex[1] + - vertex[2] * vertex[2]); - - if (radius > maxRadius) { - maxRadius = radius; - } + for (int i = 0; i < data_size; ++i) { + if (node_data[i] > max_radius) { + max_radius = node_data[i]; + } } - return maxRadius; + return max_radius; + } -std::vector getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh, +std::pair, std::vector> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh, const std::vector& rayDirection, - int numSamples) { - LogManager& logManager = LogManager::getInstance(); - quill::Logger* logger = logManager.getLogger("polyTest"); + int numSamples, std::string filename) { + Probe::LogManager& logManager = Probe::LogManager::getInstance(); + quill::Logger* logger = logManager.getLogger("log"); + LOG_INFO(logger, "Getting ray solution..."); 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; @@ -101,18 +110,50 @@ std::vector getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh, 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(i, ip); - LOG_INFO(logger, "Sample {}: Value = {:0.2E}", i, sampleValue); + sampleValue = u.GetValue(elementId, ip); + LOG_INFO(logger, "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); } } - return samples; + 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() { diff --git a/src/probe/public/probe.h b/src/probe/public/probe.h index af178a0..dfe81b8 100644 --- a/src/probe/public/probe.h +++ b/src/probe/public/probe.h @@ -6,6 +6,7 @@ #include #include #include +#include #include "mfem.hpp" #include "quill/Logger.h" @@ -30,23 +31,30 @@ 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 = "grid function"); + 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 &windowTitle = "vector", const std::string& keyset=""); + double getMeshRadius(mfem::Mesh& mesh); - std::vector getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh, - const std::vector& rayDirection, int numSamples); + 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.