Merge pull request #41 from tboudreaux/feature/meshGeneration

Feature/mesh generation

src/meshIO/private/meshIO.cpp

@@ -2,7 +2,7 @@
 //
 //   Copyright (C) 2025 -- The 4D-STAR Collaboration
 //   File Author: Emily Boudreaux
-//   Last Modified: February 16, 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
@@ -22,11 +22,12 @@
 #include <string>
 #include <iostream>
 #include <fstream>
+#include <stdexcept>
 
 #include "meshIO.h"
 
 
-MeshIO::MeshIO(const std::string mesh_file)
+MeshIO::MeshIO(const std::string &mesh_file, double scale_factor)
 {
     mesh_file_ = mesh_file;
     std::ifstream mesh_stream(mesh_file);
@@ -37,8 +38,9 @@ MeshIO::MeshIO(const std::string mesh_file)
     }
     else
     {
-        mesh_ = mfem::Mesh(mesh_stream, 1, 2);
+        mesh_ = mfem::Mesh(mesh_stream, 1, 10);
         loaded_ = true;
+        if (scale_factor != 1.0) { LinearRescale(scale_factor); }
     }
 }
 
@@ -46,6 +48,32 @@ MeshIO::~MeshIO()
 {
 }
 
+void MeshIO::LinearRescale(double scale_factor) {
+    if (!loaded_) {
+        throw std::runtime_error("Mesh not loaded before rescaling.");
+    }
+
+    if (scale_factor <= 0.0) {
+        throw std::invalid_argument("Scale factor must be positive.");
+    }
+
+    // Ensure there are nodes even for linear order meshes
+    mesh_.EnsureNodes();
+    const mfem::GridFunction *nodes = mesh_.GetNodes();
+    double *node_data = nodes->GetData();  // THIS IS KEY
+
+    int data_size = nodes->Size();
+
+    for (int i = 0; i < data_size; ++i) {
+        node_data[i] *= scale_factor;
+    }
+
+    // nodes->Update(); /updates the fespace
+    mesh_.NodesUpdated();
+ 
+
+}
+
 bool MeshIO::IsLoaded() const
 {
     return loaded_;

src/meshIO/public/meshIO.h

@@ -2,7 +2,7 @@
 //
 //   Copyright (C) 2025 -- The 4D-STAR Collaboration
 //   File Author: Emily Boudreaux
-//   Last Modified: February 16, 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
@@ -39,13 +39,19 @@ public:
      * @brief Constructor that initializes the MeshIO object with a mesh file.
      * @param mesh_file The name of the mesh file.
      */
-    MeshIO(const std::string mesh_file);
+    MeshIO(const std::string &mesh_file, double scale_factor = 1.0);
 
     /**
      * @brief Destructor for the MeshIO class.
      */
     ~MeshIO();
 
+    /**
+     * @brief Rescale the mesh by a linear factor.
+     * @param scale_factor The factor by which to scale the mesh.
+     */
+    void LinearRescale(double scale_factor);
+
     /**
      * @brief Get the mesh object.
      * @return Reference to the mesh object.

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.

utils/meshGeneration/generateMesh.py

@@ -1,15 +1,25 @@
 import pygmsh
 import meshio
+import numpy as np
+from math import sqrt
 
 import argparse
 
-def generate_spherical_mesh(radius=1, meshSize=0.1):
-    with pygmsh.geo.Geometry() as geo:
-        # Create a spherical (ball) geometry centered at (0,0,0)
-        sphere = geo.add_ball([0, 0, 0], radius)
-        # Generate a 2D mesh (i.e. surface mesh) of the sphere
-        myMesh = geo.generate_mesh(dim=3)
-    return myMesh
+def mesh_size_callback(dim, tag, x, y, z, lc, a, b, k):
+    r = np.sqrt(x**2 + y**2 + z**2)
+    size = a + (b-a) * (1-np.exp(-k*r))
+    return size
+
+def generate_spherical_mesh(radius=1.0, coreRes = 0.1, outRes = 0.1, k = 2):
+    with pygmsh.occ.Geometry() as geom:
+        originPoint = geom.add_point([0.0, 0.0, 0.0], mesh_size=coreRes)
+        geom.add_ball([0.0, 0.0, 0.0], 1.0)
+
+        geom.set_mesh_size_callback(lambda dim, tag, x, y, z, lc: mesh_size_callback(dim, tag, x, y, z, lc, coreRes, outRes, k))
+        mesh = geom.generate_mesh(verbose=True)
+
+    return mesh
+
 
 def write_mfem_mesh(meshData, filename):
     meshio.write(filename, meshData, file_format='gmsh22')
@@ -17,8 +27,11 @@ def write_mfem_mesh(meshData, filename):
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(description='Generate a spherical mesh')
     parser.add_argument('--radius', type=float, default=1.0, help='Radius of the sphere')
-    parser.add_argument('--meshSize', type=float, default=0.1, help='Mesh size')
+    parser.add_argument('--coreMeshRes', type=float, default=0.01, help='Core Mesh Resolution')
+    parser.add_argument('--outerMeshRes', type=float, default=0.1, help='Outer Mesh Resolution')
+    parser.add_argument('--k', type=float, default=2, help='Mesh size scaling factor')
+    parser.add_argument('--output', type=str, default='sphere.msh', help='Output file name')
     args = parser.parse_args()
 
-    meshData = generate_spherical_mesh(args.radius, args.meshSize)
-    write_mfem_mesh(meshData, 'sphere.msh')
+    meshData = generate_spherical_mesh(args.radius, args.coreMeshRes, args.outerMeshRes, args.k)
+    write_mfem_mesh(meshData, args.output)

utils/meshView/meshView.cpp

@@ -0,0 +1,28 @@
+#include "mfem.hpp"
+#include <string>
+#include <iostream>
+#include <cmath>
+
+#include "probe.h"
+#include "meshIO.h"
+#include "config.h"
+
+int main(int argv, char **argc) {
+    std::string CONFIG_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/testsConfig.yaml";
+    Config& config = Config::getInstance();
+    config.loadConfig(CONFIG_FILENAME);
+    // Read the mesh from the given mesh file
+    std::string meshFile = argc[1];
+    MeshIO meshIO(meshFile);
+    mfem::Mesh& mesh = meshIO.GetMesh();
+    mfem::H1_FECollection feCollection(1, mesh.SpaceDimension());
+    mfem::FiniteElementSpace feSpace(&mesh, &feCollection);
+    mfem::GridFunction u(&feSpace);
+    mfem::FunctionCoefficient initCoeff (
+        [&mesh](const mfem::Vector &x) {
+            return 1.0;
+        }
+    );
+    u.ProjectCoefficient(initCoeff);
+    Probe::glVisView(u, mesh, meshFile);
+}

utils/meshView/meson.build

@@ -0,0 +1,8 @@
+sources = files('meshView.cpp')
+
+meshView_exe = executable(
+    'meshView',
+    sources,
+    dependencies: [probe_dep, quill_dep, config_dep, mfem_dep, meshio_dep],
+    install_rpath: '@loader_path/../../src'  # Ensure runtime library path resolves correctly
+)