feat(probe): added getRadius and start of ray cast solution

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();