From ce299ebfe1ca196fc64c96d86fccf45438dd5ea2 Mon Sep 17 00:00:00 2001
From: Emily Boudreaux <emily@boudreauxmail.com>
Date: Sun, 23 Feb 2025 14:11:50 -0500
Subject: [PATCH] feat(probe): added getRadius and start of ray cast solution

---
 src/probe/private/probe.cpp | 52 +++++++++++++++++++++++++++++++++++++
 src/probe/public/probe.h    |  6 +++++
 2 files changed, 58 insertions(+)

diff --git a/src/probe/private/probe.cpp b/src/probe/private/probe.cpp
index 6d26f66..921610c 100644
--- a/src/probe/private/probe.cpp
+++ b/src/probe/private/probe.cpp
@@ -3,10 +3,12 @@
 #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 "mfem.hpp"
 
@@ -39,6 +41,56 @@ void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh,
   }
 }
 
+double getMeshRadius(mfem::Mesh& mesh) {
+  int numVertices = mesh.GetNV(); // Number of vertices
+  double maxRadius = 0.0;
+
+  for (int i = 0; i < numVertices; i++) {
+      double* vertex; 
+      vertex = mesh.GetVertex(i); // Get vertex coordinates
+
+      // 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;
+      }
+  }
+  return maxRadius;
+}
+
+std::vector<double> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh,
+                                   const std::vector<double>& rayDirection,
+                                   int numSamples) {
+  LogManager& logManager = LogManager::getInstance();
+  quill::Logger* logger = logManager.getLogger("polyTest");
+  std::vector<double> samples;
+  samples.reserve(numSamples);
+  double x, y, z, r;
+  double radius = getMeshRadius(mesh);
+  mfem::Vector rayOrigin(3); rayOrigin = 0.0;
+  mfem::DenseMatrix rayPoints(3, numSamples);
+
+  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;
+    LOG_INFO(logger, "Ray point {}: ({}, {}, {})", i, x, y, z);
+  }
+
+  mfem::Array<int> elementIds;
+  mfem::Array<mfem::IntegrationPoint> ips;
+  mesh.FindPoints(rayPoints, elementIds, ips);
+  return samples;
+}
+
 LogManager::LogManager() {
   Config& config = Config::getInstance();
   quill::Backend::start();
diff --git a/src/probe/public/probe.h b/src/probe/public/probe.h
index 1ef7e09..abd19ca 100644
--- a/src/probe/public/probe.h
+++ b/src/probe/public/probe.h
@@ -33,6 +33,12 @@ namespace Probe {
    */
   void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh,
                 const std::string& windowTitle = "solution");
+  
+  double getMeshRadius(mfem::Mesh& mesh);
+
+  std::vector<double> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh,
+                                    const std::vector<double>& rayDirection,
+                                    int numSamples);
 
   /**
    * @brief Class to manage logging operations.