feat(probe): added glvisview helped method to probe namespace

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

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 <stdexcept>              
 #include <string>
 #include <iostream>
 #include <chrono>
+#include <cmath>
+#include <vector>
+#include <utility>
+#include <filesystem>
 
 #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<bool>("Probe:GLVis:Visualization", true)) {
+    LOG_INFO(logger, "Visualizing solution using GLVis...");
+    LOG_INFO(logger, "Window title: {}", windowTitle);
+    if (keyset == "") {
+      usedKeyset = config.get<std::string>("Probe:GLVis:DefaultKeyset", "");
+    } else {
+      usedKeyset = keyset;
+    }
+    LOG_INFO(logger, "Keyset: {}", usedKeyset);
     std::string vishost = config.get<std::string>("Probe:GLVis:Host", "localhost");
-    int visport = config.get<int>("Probe:GLVis:Port", 19916); // Changed default port
+    int visport = config.get<int>("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<double>, std::vector<double>> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh,
+                                   const std::vector<double>& 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<bool>("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<double> 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<double> 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<int> elementIds;
+  mfem::Array<mfem::IntegrationPoint> 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<double>, std::vector<double>> 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<double>, std::vector<double>> getRaySolution(mfem::Vector &vec, mfem::FiniteElementSpace &fes,
+                                   const std::vector<double>& 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();

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 <string>
 #include <map>
 #include <vector>
-#include <iostream>
+#include <utility>
 
 #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<double>, std::vector<double>> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh,
+    const std::vector<double>& rayDirection, int numSamples, std::string filename="");
+
+  std::pair<std::vector<double>, std::vector<double>> getRaySolution(mfem::Vector &vec, mfem::FiniteElementSpace &fes,
+    const std::vector<double>& rayDirection, int numSamples, std::string filename="");
+
 
   /**
    * @brief Class to manage logging operations.