diff --git a/assets/dynamic/mesh/core.msh b/assets/dynamic/mesh/core.msh new file mode 100644 index 0000000..bda83cf Binary files /dev/null and b/assets/dynamic/mesh/core.msh differ diff --git a/assets/dynamic/mesh/core_hires.msh b/assets/dynamic/mesh/core_hires.msh new file mode 100644 index 0000000..9967986 Binary files /dev/null and b/assets/dynamic/mesh/core_hires.msh differ diff --git a/assets/dynamic/mesh/core_lowres.msh b/assets/dynamic/mesh/core_lowres.msh new file mode 100644 index 0000000..f0dbfdb Binary files /dev/null and b/assets/dynamic/mesh/core_lowres.msh differ diff --git a/assets/dynamic/mesh/core_midres.msh b/assets/dynamic/mesh/core_midres.msh new file mode 100644 index 0000000..4b62907 Binary files /dev/null and b/assets/dynamic/mesh/core_midres.msh differ diff --git a/assets/dynamic/mesh/sphere.msh b/assets/dynamic/mesh/sphere.msh index 416c655..1c68e68 100644 Binary files a/assets/dynamic/mesh/sphere.msh and b/assets/dynamic/mesh/sphere.msh differ diff --git a/src/meshIO/private/meshIO.cpp b/src/meshIO/private/meshIO.cpp index a894a14..66d028f 100644 --- a/src/meshIO/private/meshIO.cpp +++ b/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 #include #include +#include #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_; diff --git a/src/meshIO/public/meshIO.h b/src/meshIO/public/meshIO.h index 1a2c4db..f96cf63 100644 --- a/src/meshIO/public/meshIO.h +++ b/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. diff --git a/src/probe/meson.build b/src/probe/meson.build index 653c5eb..cf6accf 100644 --- a/src/probe/meson.build +++ b/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 diff --git a/src/probe/private/probe.cpp b/src/probe/private/probe.cpp index 8ae8f11..78d512a 100644 --- a/src/probe/private/probe.cpp +++ b/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 #include #include #include +#include +#include +#include +#include #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("Probe:GLVis:Visualization", true)) { + LOG_INFO(logger, "Visualizing solution using GLVis..."); + LOG_INFO(logger, "Window title: {}", windowTitle); + if (keyset == "") { + usedKeyset = config.get("Probe:GLVis:DefaultKeyset", ""); + } else { + usedKeyset = keyset; + } + LOG_INFO(logger, "Keyset: {}", usedKeyset); std::string vishost = config.get("Probe:GLVis:Host", "localhost"); - int visport = config.get("Probe:GLVis:Port", 19916); // Changed default port + int visport = config.get("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> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh, + const std::vector& 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("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 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; + // 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 elementIds; + 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(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> 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() { Config& config = Config::getInstance(); quill::Backend::start(); diff --git a/src/probe/public/probe.h b/src/probe/public/probe.h index ce863d0..f0d2709 100644 --- a/src/probe/public/probe.h +++ b/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 #include #include -#include +#include #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> 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. diff --git a/utils/meshGeneration/generateMesh.py b/utils/meshGeneration/generateMesh.py index eb76bb2..b758f43 100644 --- a/utils/meshGeneration/generateMesh.py +++ b/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') \ No newline at end of file + meshData = generate_spherical_mesh(args.radius, args.coreMeshRes, args.outerMeshRes, args.k) + write_mfem_mesh(meshData, args.output) diff --git a/utils/meshView/meshView.cpp b/utils/meshView/meshView.cpp new file mode 100644 index 0000000..0190c9e --- /dev/null +++ b/utils/meshView/meshView.cpp @@ -0,0 +1,28 @@ +#include "mfem.hpp" +#include +#include +#include + +#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); +} \ No newline at end of file diff --git a/utils/meshView/meson.build b/utils/meshView/meson.build new file mode 100644 index 0000000..5815967 --- /dev/null +++ b/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 +) \ No newline at end of file