From 39e5117a245a0a2e465dab1c0200bc2fdd6d4768 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 1 Jul 2026 11:14:12 -0400 Subject: [PATCH] feat(python): added python bindings --- build-config/meson.build | 18 +- build-config/pybind/meson.build | 3 + build-config/python/meson.build | 5 + build-python/meson.build | 43 ++ meson.build | 8 +- meson_options.txt | 3 +- pyproject.toml | 25 + src/include/stroid/IO/mesh.h | 45 +- src/include/stroid/config/config.h | 10 + src/include/stroid/exceptions/exceptions.h | 3 + src/include/stroid/exceptions/stroid_error.h | 25 + src/include/stroid/meson.build | 4 +- src/include/stroid/refinement/uniform.h | 7 + src/include/stroid/{stroid.h.in => stroid.h} | 63 +-- src/include/stroid/utils/mesh_stats.h | 174 ++++++ src/include/stroid/utils/mesh_utils.h | 8 + src/include/stroid/utils/types.h | 65 +++ src/include/stroid/version.h.in | 46 ++ src/lib/IO/mesh.cpp | 505 +++++++++++++++++- src/lib/refinement/uniform.cpp | 39 ++ src/lib/topology/mapping.cpp | 188 +++++-- src/lib/topology/optimize.cpp | 9 +- src/lib/utils/mesh_stats.cpp | 485 +++++++++++++++++ src/lib/utils/mesh_utils.cpp | 44 +- src/meson.build | 76 ++- src/python/IO/bindings.cpp | 77 +++ src/python/IO/bindings.h | 5 + src/python/bindings.cpp | 36 ++ src/python/config/bindings.cpp | 202 +++++++ src/python/config/bindings.h | 5 + src/python/exceptions/bindings.cpp | 14 + src/python/exceptions/bindings.h | 5 + src/python/refinement/bindings.cpp | 11 + src/python/refinement/bindings.h | 5 + src/python/stroid/__init__.py | 45 ++ src/python/utils/bindings.cpp | 187 +++++++ src/python/utils/bindings.h | 5 + .../packagefiles/pybind11/LICENSE.build | 19 + subprojects/packagefiles/pybind11/meson.build | 8 + subprojects/pybind11.wrap | 8 + 40 files changed, 2434 insertions(+), 99 deletions(-) create mode 100644 build-config/pybind/meson.build create mode 100644 build-config/python/meson.build create mode 100644 build-python/meson.build create mode 100644 pyproject.toml create mode 100644 src/include/stroid/exceptions/exceptions.h create mode 100644 src/include/stroid/exceptions/stroid_error.h create mode 100644 src/include/stroid/refinement/uniform.h rename src/include/stroid/{stroid.h.in => stroid.h} (60%) create mode 100644 src/include/stroid/utils/mesh_stats.h create mode 100644 src/include/stroid/utils/types.h create mode 100644 src/include/stroid/version.h.in create mode 100644 src/lib/refinement/uniform.cpp create mode 100644 src/lib/utils/mesh_stats.cpp create mode 100644 src/python/IO/bindings.cpp create mode 100644 src/python/IO/bindings.h create mode 100644 src/python/bindings.cpp create mode 100644 src/python/config/bindings.cpp create mode 100644 src/python/config/bindings.h create mode 100644 src/python/exceptions/bindings.cpp create mode 100644 src/python/exceptions/bindings.h create mode 100644 src/python/refinement/bindings.cpp create mode 100644 src/python/refinement/bindings.h create mode 100644 src/python/stroid/__init__.py create mode 100644 src/python/utils/bindings.cpp create mode 100644 src/python/utils/bindings.h create mode 100644 subprojects/packagefiles/pybind11/LICENSE.build create mode 100644 subprojects/packagefiles/pybind11/meson.build create mode 100644 subprojects/pybind11.wrap diff --git a/build-config/meson.build b/build-config/meson.build index ea4e282..fc72fdd 100644 --- a/build-config/meson.build +++ b/build-config/meson.build @@ -1,4 +1,20 @@ subdir('mfem') subdir('libconfig') subdir('CLI11') -subdir('magic_enum') \ No newline at end of file +subdir('magic_enum') + +if get_option('build_python') + subdir('python') + subdir('pybind') +endif + +if get_option('build_python') + stroid_pkg_dir = py_installation.get_install_dir() / 'stroid' + stroid_includedir = stroid_pkg_dir / 'include' + stroid_libdir = stroid_pkg_dir / 'lib' + stroid_pcdir = stroid_libdir / 'pkgconfig' +else + stroid_includedir = get_option('includedir') + stroid_libdir = get_option('libdir') + stroid_pcdir = get_option('libdir') / 'pkgconfig' +endif \ No newline at end of file diff --git a/build-config/pybind/meson.build b/build-config/pybind/meson.build new file mode 100644 index 0000000..4929711 --- /dev/null +++ b/build-config/pybind/meson.build @@ -0,0 +1,3 @@ +pybind11_proj = subproject('pybind11') +pybind11_dep = pybind11_proj.get_variable('pybind11_dep') +python3_dep = dependency('python3') \ No newline at end of file diff --git a/build-config/python/meson.build b/build-config/python/meson.build new file mode 100644 index 0000000..b8e91e4 --- /dev/null +++ b/build-config/python/meson.build @@ -0,0 +1,5 @@ +py_installation = import('python').find_installation('python3', pure: false) + +py_dep = py_installation.dependency() +py_module_prefix = '' +py_module_suffix = 'so' \ No newline at end of file diff --git a/build-python/meson.build b/build-python/meson.build new file mode 100644 index 0000000..8f915a0 --- /dev/null +++ b/build-python/meson.build @@ -0,0 +1,43 @@ +if get_option('build_python') + message('Building Python bindings...') + + stroid_py_deps = [ + py_dep, + pybind11_dep, + stroid_dep + ] + + if host_machine.system() == 'darwin' + stroid_ext_rpath = '@loader_path/lib' + else + stroid_ext_rpath = '$ORIGIN/lib' + endif + + py_sources = [ + meson.project_source_root() + '/src/python/bindings.cpp', + meson.project_source_root() + '/src/python/config/bindings.cpp', + meson.project_source_root() + '/src/python/exceptions/bindings.cpp', + meson.project_source_root() + '/src/python/IO/bindings.cpp', + meson.project_source_root() + '/src/python/refinement/bindings.cpp', + meson.project_source_root() + '/src/python/utils/bindings.cpp', + ] + + py_mod = py_installation.extension_module( + '_stroid', + sources: py_sources, + dependencies: stroid_py_deps, + install: true, + link_args: stroid_ext_rpath_args, + build_rpath: stroid_ext_rpath, + install_rpath: stroid_ext_rpath, + subdir: 'stroid', + ) + + py_installation.install_sources( + meson.project_source_root() + '/src/python/stroid/__init__.py', + subdir: 'stroid', + ) + +else + message('Python bindings disabled') +endif \ No newline at end of file diff --git a/meson.build b/meson.build index bc13b2c..a9c7a03 100644 --- a/meson.build +++ b/meson.build @@ -1,4 +1,4 @@ -project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.3.0', default_options : ['cpp_std=c++23']) +project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.4.0', default_options : ['cpp_std=c++23']) subdir('build-check') @@ -13,6 +13,10 @@ if get_option('build_tools') subdir('tools') endif +if get_option('build_python') + subdir('build-python') +endif + if get_option('pkg_config') pkg = import('pkgconfig') pkg.generate( @@ -20,7 +24,7 @@ if get_option('pkg_config') description: 'Stroid multi-block curvilinear mesh generation library', version: meson.project_version(), libraries: [ - stroid_lib + libstroid ], subdirs: ['stroid'], filebase: 'stroid', diff --git a/meson_options.txt b/meson_options.txt index 2e318cc..8b3ab8b 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -1,3 +1,4 @@ option('pkg_config', type: 'boolean', value: false, description: 'generate pkg-config file for stroid') option('build_tests', type: 'boolean', value: true, description: 'compile subproject tests') -option('build_tools', type: 'boolean', value: true, description: 'compile stroid command line tools') \ No newline at end of file +option('build_tools', type: 'boolean', value: true, description: 'compile stroid command line tools') +option('build_python', type: 'boolean', value: true, description: 'compile stroid python bindings') \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..45706da --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,25 @@ +[build-system] +requires = ["meson-python>=0.19.0", "meson>=1.9.1", "pybind11==3.0.0", "fourdst==0.10.6"] +build-backend = "mesonpy" + +[project] +name = "stroid" +dynamic = ["version"] +description = "O-grid mesh generation with multiple domains" +readme = "README.md" +license = { file = "LICENSE.txt" } + +authors = [ + {name = "Emily M. Boudreaux", email = "emily@boudreauxmail.com"}, +] +maintainers = [ + {name = "Emily M. Boudreaux", email = "emily@boudreauxmail.com"} +] + +[tool.meson-python.args] +setup = [ + '-Dbuild_tools=false', + '-Dbuild_tests=false', + '-Dpkg_config=false' +] +install = ['--skip-subprojects'] \ No newline at end of file diff --git a/src/include/stroid/IO/mesh.h b/src/include/stroid/IO/mesh.h index 2c1c6cc..0fce639 100644 --- a/src/include/stroid/IO/mesh.h +++ b/src/include/stroid/IO/mesh.h @@ -1,7 +1,12 @@ #pragma once #include +#include +#include + #include "mfem.hpp" +#include "stroid/utils/types.h" + namespace stroid::IO { /** * @brief Visualization modes for GLVis display. @@ -15,18 +20,43 @@ namespace stroid::IO { BOUNDARY_ELEMENT_ID }; + void SaveStroidMesh(const StroidMesh& mesh, const std::string& filename, const std::string& comment=""); + /** * @brief Save a mesh to MFEM's native `.mesh` format. * @param mesh Mesh to serialize. * @param filename Output path (including extension). */ void SaveMesh(const mfem::Mesh& mesh, const std::string& filename); + + /** + * @brief Overload of SaveMesh which accepts a StroidMesh type and will internally unpack it + * @param mesh StroidMesh to serialize. + * @param filename Path to save to + * + * @note This function is a utility wrapper to save a StroidMesh object in MFEM's native .mesh format. Data other than the mesh pointer + * in StroidMesh **will not be saved** (e.g. the reference mesh, the number of refinement levels, etc..). If you need to serialize an + * entire StroidMesh then please use the stroid::IO::SaveStroidMesh function + */ + void SaveMesh(const stroid::StroidMesh& mesh, const std::string& filename); /** * @brief Save a mesh as a ParaView VTU dataset. * @param mesh Mesh to export. * @param exportName Output base name (ParaView will add extensions). */ void SaveVTU(mfem::Mesh& mesh, const std::string& exportName); + + /** + * @brief Overload of SaveVTU which accepts a StroidMesh type and will internally unpack it + * @param mesh StroidMesh to serialize. + * @param filename Path to save to + * + * @note This function is a utility wrapper to save a StroidMesh object in MFEM's native .mesh format. Data other than the mesh pointer + * in StroidMesh **will not be saved** (e.g. the reference mesh, the number of refinement levels, etc..). If you need to serialize an + * entire StroidMesh then please use the stroid::IO::SaveStroidVTU function + */ + void SaveVTU(const stroid::StroidMesh& mesh, const std::string& exportName); + /** * @brief Stream a mesh to a running GLVis server for interactive viewing. * @param mesh Mesh to display. @@ -36,9 +66,22 @@ namespace stroid::IO { * @param visport GLVis server port. */ void ViewMesh(mfem::Mesh &mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport); + + void ViewMesh(const stroid::StroidMesh& mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport); + /** * @brief Visualize boundary face valence (1=surface, 2=internal). * @param mesh Mesh whose boundary faces are inspected. */ - void VisualizeFaceValence(mfem::Mesh& mesh); + void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport); + + void VisualizeFaceValence(const stroid::StroidMesh& mesh, const std::string &vishost, int visport); + + std::expected ParseStroidMesh(std::istream& is); + std::expected LoadStroidMesh(const std::string& filename); + +#ifdef MFEM_USE_MPI + std::expected ParseStroidMesh(std::istream& is, MPI_Comm comm); + std::expected LoadStroidMesh(const std::string& filename, MPI_Comm comm); +#endif } \ No newline at end of file diff --git a/src/include/stroid/config/config.h b/src/include/stroid/config/config.h index f9227f0..918245a 100644 --- a/src/include/stroid/config/config.h +++ b/src/include/stroid/config/config.h @@ -1,6 +1,8 @@ #pragma once #include +#include +#include namespace stroid::config { @@ -76,6 +78,13 @@ namespace stroid::config { */ std::optional core_steepness = 1.0; + /** + * @brief Continuity order for the core-envelope transition (0 = discontinuous, 1 = C1, 2 = C2). + * @section toml + * - [main].continuity_order + */ + std::optional continuity_order = 2; + /** * @brief Boundary attribute id for stellar surface * @section toml @@ -112,5 +121,6 @@ namespace stroid::config { std::optional vacuum_id = 3; std::optional optimization_methods = OptimizationMethods{true, true}; + }; } diff --git a/src/include/stroid/exceptions/exceptions.h b/src/include/stroid/exceptions/exceptions.h new file mode 100644 index 0000000..da53107 --- /dev/null +++ b/src/include/stroid/exceptions/exceptions.h @@ -0,0 +1,3 @@ +#pragma once + +#include "stroid/exceptions/stroid_error.h" \ No newline at end of file diff --git a/src/include/stroid/exceptions/stroid_error.h b/src/include/stroid/exceptions/stroid_error.h new file mode 100644 index 0000000..dd395d1 --- /dev/null +++ b/src/include/stroid/exceptions/stroid_error.h @@ -0,0 +1,25 @@ +#pragma once +#include +#include + +namespace stroid::exceptions { + class StroidError : public std::exception { + public: + explicit StroidError(std::string message) : m_msg(std::move(message)) {} + const char* what() const noexcept override { return m_msg.c_str(); } + private: + std::string m_msg; + }; + + class StroidContinuityError : public StroidError { + using StroidError::StroidError; + }; + + class StroidMeshError : public StroidError { + using StroidError::StroidError; + }; + + class StroidMissingReferenceMesh : public StroidMeshError { + using StroidMeshError::StroidMeshError; + }; +} diff --git a/src/include/stroid/meson.build b/src/include/stroid/meson.build index d285070..9d07260 100644 --- a/src/include/stroid/meson.build +++ b/src/include/stroid/meson.build @@ -21,8 +21,8 @@ config.set('STROID_VERSION_PATCH', ver_parts[2]) config.set('STROID_VERSION_TAG', ver_parts[3]) configure_file( - input : 'stroid.h.in', - output : 'stroid.h', + input : 'version.h.in', + output : 'version.h', configuration : config , install: true, install_dir: get_option('includedir') / 'stroid' diff --git a/src/include/stroid/refinement/uniform.h b/src/include/stroid/refinement/uniform.h new file mode 100644 index 0000000..6b028dd --- /dev/null +++ b/src/include/stroid/refinement/uniform.h @@ -0,0 +1,7 @@ +#pragma once + +#include "stroid/utils/types.h" + +namespace stroid::refinement { + void UniformRefinement(StroidMesh& mesh, size_t levels); +} \ No newline at end of file diff --git a/src/include/stroid/stroid.h.in b/src/include/stroid/stroid.h similarity index 60% rename from src/include/stroid/stroid.h.in rename to src/include/stroid/stroid.h index f221d63..e77aff0 100644 --- a/src/include/stroid/stroid.h.in +++ b/src/include/stroid/stroid.h @@ -7,6 +7,9 @@ #include "stroid/topology/optimize.h" #include "stroid/utils/mesh_utils.h" #include "stroid/IO/mesh.h" +#include "stroid/utils/types.h" +#include "stroid/refinement/uniform.h" +#include "stroid/utils/mesh_stats.h" /** * @namespace stroid @@ -45,46 +48,36 @@ * @endcode */ namespace stroid { - /** - * @brief Version helpers for the stroid library. - */ - struct version { - static constexpr int major = @STROID_VERSION_MAJOR@; - static constexpr int minor = @STROID_VERSION_MINOR@; - static constexpr int patch = @STROID_VERSION_PATCH@; - static constexpr const char* tag = "@STROID_VERSION_TAG@"; + inline StroidMesh GenerateMesh(const fourdst::config::Config& cfg) { + StroidMesh sm; + sm.config = *cfg; + auto reference = stroid::topology::BuildSkeleton(cfg); + stroid::topology::Finalize(*reference, cfg); + sm.refinement_levels = cfg->refinement_levels.value_or(0); - static std::string toString() { - std::string versionStr = std::to_string(major) + "." + - std::to_string(minor) + "." + - std::to_string(patch); - if (std::string(tag) != "") { - versionStr += "-" + std::string(tag); - } - return versionStr; + sm.reference_mesh = std::move(reference); + sm.mesh = utils::BuildProjected(*sm.reference_mesh, cfg); + if (cfg->optimization_methods.has_value() && cfg->optimization_methods.value().tmop.has_value() && cfg->optimization_methods.value().tmop.value()) { + stroid::topology::ApplyTMOP(*sm.mesh, cfg); } + return sm; + } + inline StroidMesh GenerateMesh(const stroid::config::MeshConfig& config) { + fourdst::config::Config cfg; + auto Mutator = [&config](config::MeshConfig& orig) { + orig = config; + }; - friend std::ostream& operator<<(std::ostream& os, const version&) { - os << toString(); - return os; - } - }; + cfg.mutate(Mutator); + return GenerateMesh(cfg); + } + inline StroidMesh GenerateMesh(const std::string& filename) { + fourdst::config::Config config; + config.load(filename); + return GenerateMesh(config); + } } -/** - * @namespace std - * @brief Standard library extensions used by stroid. - * - * Provides a `std::formatter` specialization for `stroid::version` so it can - * be used with `std::format` and related APIs. - */ -// Overload format struct -template <> -struct std::formatter : std::formatter { - auto format(const stroid::version& v, auto& ctx) { - return std::formatter::format(stroid::version::toString(), ctx); - } -}; /** * @namespace stroid::config diff --git a/src/include/stroid/utils/mesh_stats.h b/src/include/stroid/utils/mesh_stats.h new file mode 100644 index 0000000..32c5f45 --- /dev/null +++ b/src/include/stroid/utils/mesh_stats.h @@ -0,0 +1,174 @@ +#pragma once + +#include "mfem.hpp" +#include "stroid/utils/types.h" +#include "stroid/config/config.h" + +#include +#include +#include +#include + +namespace stroid::stats { + enum class MeshStatFeatures : uint32_t { + NONE = 0u, + RADIUS = 1u << 0, + AXES = 1u << 1, + ELLIPTICITY = 1u << 2, + BOWING = 1u << 3, + CONFORMITY = 1u << 4, + JACOBIAN = 1u << 5, + VOLUME_AREA = 1u << 6, + ELEMENT_COUNT = 1u << 7, + MESH_SIZE = 1u << 8, + OUTER_BOUNDS = 1u << 9, + CENTROID = 1u << 10, + CONFIG_META = 1u << 11, + BOUNDING_BOX = 1u << 12, + }; + + constexpr MeshStatFeatures operator|(MeshStatFeatures lhs, MeshStatFeatures rhs) { + return static_cast(static_cast(lhs) | static_cast(rhs)); + } + + constexpr MeshStatFeatures operator&(MeshStatFeatures lhs, MeshStatFeatures rhs) { + return static_cast(static_cast(lhs) & static_cast(rhs)); + } + + constexpr bool has_feature(MeshStatFeatures feature, MeshStatFeatures set) { + return (static_cast(set) & static_cast(feature)) != 0u; + } + + inline constexpr MeshStatFeatures MESH_STAT_DEFAULT = + MeshStatFeatures::RADIUS | MeshStatFeatures::AXES | MeshStatFeatures::ELLIPTICITY | + MeshStatFeatures::CONFORMITY | MeshStatFeatures::CONFIG_META; + + inline constexpr auto MESH_STAT_ALL = static_cast(0xFFFFFFFFu); + + struct RadiusStats { + double min = 0, max = 0, mean = 0, stddev = 0; + long n_samples = 0; + }; + + struct AxisStats { + double semi_major = 0; + double semi_minor = 0; + }; + + struct EllipticityStats { + double flattening = 0; + double polar_equatorial = 1; + double radius_uniformity = 1; + }; + + struct BowingStats { + double max_inward = 0; + double max_outward = 0; + double rms = 0; + double worst_at_radius = 0; + }; + + struct ConformityStats { + bool conforming = true; + long n_nonconforming_faces = 0; + }; + + struct JacobianStats { + double detJ_min; + double detJ_max; + long n_flipped; + double min_detJ_ratio; + double worst_ratio_at_radius; + double detJ_min_at_radius; + long n_elements; + }; + struct VolumeAreaStats { + double stellar_volume = 0, surface_area = 0; + double analytic_volume = 0, analytic_area = 0; + }; + + struct ElementCounts { + long total = 0, core = 0, envelope = 0, vacuum = 0, other = 0; + long n_vertices = 0; + }; + + struct MeshSizeStats { + double h_min = 0, h_max = 0, h_mean = 0, h_stddev = 0; + }; + + struct OuterBoundsStats { + double min = 0, max = 0, mean = 0; + long n_samples = 0; + }; + + struct CentroidStats { + double x = 0, y = 0, z = 0, offset = 0; + }; + + struct ConfigMeta { + double r_core = 0, r_star = 0, flattening = 0, r_infinity = 0; + int geom_order = 0; + size_t refinement_levels = 0; + bool has_external_domain = true; + }; + + struct BoundingBox { + double xMin = 0, xMax = 0, yMin = 0, yMax = 0, zMin = 0, zMax = 0; + bool valid = false; + + [[nodiscard]] double dx() const {return xMax - xMin;} + [[nodiscard]] double dy() const {return yMax - yMin;} + [[nodiscard]] double dz() const {return zMax - zMin;} + [[nodiscard]] double diag() const { + const double a = dx(), b = dy(), c = dz(); + return std::sqrt(a*a + b*b + c*c); + } + }; + + struct BoundingBoxStats { + BoundingBox core; + BoundingBox star; + BoundingBox vacuum; + }; + + struct MeshStats { + MeshStatFeatures computed = MeshStatFeatures::NONE; + std::optional radius; + std::optional axes; + std::optional ellipticity; + std::optional bowing; + std::optional conformity; + std::optional jacobian; + std::optional jacobian_stellar; + std::optional jacobian_vacuum; + std::optional volume; + std::optional element_counts; + std::optional mesh_size; + std::optional outer_bounds; + std::optional centroid; + std::optional config_meta; + std::optional bounding_box; + + std::vector warnings; + std::vector errors; + }; + + MeshStats ComputeMeshStats(const StroidMesh& sm, MeshStatFeatures features = MESH_STAT_DEFAULT, int sample_order = -1); + + std::string to_string(const MeshStats& s); + + inline std::ostream& operator<<(std::ostream& os, const MeshStats& s) { + return os << to_string(s); + } +} + +template <> +struct std::formatter { + static constexpr auto parse(const std::format_parse_context& ctx) { + return ctx.begin(); + } + + static auto format(const stroid::stats::MeshStats &s, std::format_context& ctx) { + return std::format_to(ctx.out(), "{}", stroid::stats::to_string(s)); + } +}; diff --git a/src/include/stroid/utils/mesh_utils.h b/src/include/stroid/utils/mesh_utils.h index 6de6c53..9d3fd82 100644 --- a/src/include/stroid/utils/mesh_utils.h +++ b/src/include/stroid/utils/mesh_utils.h @@ -2,6 +2,9 @@ #include "mfem.hpp" +#include "stroid/config/config.h" +#include "fourdst/config/config.h" + namespace stroid::utils { /** * @brief Mark elements with negative Jacobian determinant. @@ -15,4 +18,9 @@ namespace stroid::utils { * @param mesh Mesh to scan and update in-place. */ void MarkFlippedBoundaryElements(mfem::Mesh& mesh); + + void ExportJacobianRadialProfile(mfem::Mesh& mesh, const std::string& filename); + + std::unique_ptr BuildProjected(const mfem::Mesh& reference, const fourdst::config::Config& cfg); + } \ No newline at end of file diff --git a/src/include/stroid/utils/types.h b/src/include/stroid/utils/types.h new file mode 100644 index 0000000..2b95332 --- /dev/null +++ b/src/include/stroid/utils/types.h @@ -0,0 +1,65 @@ +#pragma once + +#include "mfem.hpp" + +#include "stroid/config/config.h" + +#include +#include +#include +#include +#include + +namespace stroid { + enum class MFEM_MESH_TYPE { + SERIAL, + PARALLEL + }; + + struct StroidMesh { + MFEM_MESH_TYPE type; + std::unique_ptr mesh; + std::unique_ptr reference_mesh; + config::MeshConfig config; + size_t refinement_levels; + + [[nodiscard]] std::expected as_mesh() const { + if (type == MFEM_MESH_TYPE::SERIAL) { + return mesh.get(); + } + return std::unexpected{"Mesh is not serial. Try calling as_par_mesh()"}; + } + + [[nodiscard]] std::expected ref_as_mesh() const { + if (type == MFEM_MESH_TYPE::SERIAL) { + return reference_mesh.get(); + } + return std::unexpected{"Reference mesh is not serial. Try calling as_par_mesh()"}; + } + + [[nodiscard]] std::expected>, std::string> mesh_stats(bool use_ref_mesh = false) const { + if (type != MFEM_MESH_TYPE::SERIAL) { + return std::unexpected{"Mesh is not serial. Mesh stats currently only supports serial meshes."}; + } + + mfem::Mesh* umesh; + if (use_ref_mesh) { + umesh = reference_mesh.get(); + } else { + umesh = mesh.get(); + } + std::unordered_map> mesh_stats; + mesh_stats.emplace("num_elements", umesh->GetNE()); + mesh_stats.emplace("num_vertices", umesh->GetNV()); + mesh_stats.emplace("num_edges", umesh->GetNEdges()); + mesh_stats.emplace("num_faces", umesh->GetNFaces()); + mesh_stats.emplace("num_boundary_elements", umesh->GetNBE()); + mesh_stats.emplace("max_bdr_attribute_id", umesh->bdr_attributes.Max()); + mesh_stats.emplace("min_bdr_attribute_id", umesh->bdr_attributes.Min()); + mesh_stats.emplace("max_element_attribute_id", umesh->attributes.Max()); + mesh_stats.emplace("min_element_attribute_id", umesh->attributes.Min()); + + return mesh_stats; + } + }; +} diff --git a/src/include/stroid/version.h.in b/src/include/stroid/version.h.in new file mode 100644 index 0000000..a023545 --- /dev/null +++ b/src/include/stroid/version.h.in @@ -0,0 +1,46 @@ +#pragma once + +#include +#include + +namespace stroid { + /** + * @brief Version helpers for the stroid library. + */ + struct version { + static constexpr int major = @STROID_VERSION_MAJOR@; + static constexpr int minor = @STROID_VERSION_MINOR@; + static constexpr int patch = @STROID_VERSION_PATCH@; + static constexpr const char* tag = "@STROID_VERSION_TAG@"; + + static std::string toString() { + std::string versionStr = std::to_string(major) + "." + + std::to_string(minor) + "." + + std::to_string(patch); + if (std::string(tag) != "") { + versionStr += "-" + std::string(tag); + } + return versionStr; + } + + friend std::ostream& operator<<(std::ostream& os, const version&) { + os << toString(); + return os; + } + }; +} + +/** + * @namespace std + * @brief Standard library extensions used by stroid. + * + * Provides a `std::formatter` specialization for `stroid::version` so it can + * be used with `std::format` and related APIs. + */ +// Overload format struct +template <> +struct std::formatter : std::formatter { + auto format(const stroid::version& v, auto& ctx) { + return std::formatter::format(stroid::version::toString(), ctx); + } +}; diff --git a/src/lib/IO/mesh.cpp b/src/lib/IO/mesh.cpp index 434aaa6..333e660 100644 --- a/src/lib/IO/mesh.cpp +++ b/src/lib/IO/mesh.cpp @@ -2,12 +2,428 @@ #include "stroid/config/config.h" #include "stroid/IO/mesh.h" +#include + +#include "stroid/version.h" + #include #include #include +#include +#include namespace stroid::IO { + namespace { + std::string format_header(const StroidMesh& mesh, const std::string& comment) { + auto now = std::chrono::system_clock::now(); + version v; + + std::stringstream vs; + vs << v; + + std::string header = std::format(R"(# STROID MESH +# NOTE: STROID MESH IS A THIN WRAPPER AROUND MFEM's NATIVE MESH FORMAT +# STRUCTURE: +# - Type : Serial or Parallel (S for Serial, P for Parallel) +# - mesh : the primary computational domain which can be of n order and be h-refined +# - reference mesh : a reference, linear order mesh, used to ensure that the primary mesh remains well formed +# - config : The configuration options initially used to generate the mesh +# - refinement-levels : the total number of refinement levels the primary mesh has been subjected too +# NOTE: EACH BLOCK OF DATA IS STORED BETWEEN "BEGIN BLOCK \n ... \nEND BLOCK +# PARSING THE UNDERLYING MFEM NATIVE MESH FORMAT CAN BE DONE WITH MFEM'S STREAM READER +# IF YOU EXTRACT THE RAW CONTENTS BETWEEN THOSE LINES +BEGIN BLOCK HEADER + MESH_TYPE:{} + REFINEMENT_LEVELS:{} + DATE_CREATED:{:%Y-%m-%d} + COMMENT:{} + STROID_VERSION:{} +END BLOCK HEADER)", + mesh.type == MFEM_MESH_TYPE::PARALLEL ? "P" : "S", + mesh.refinement_levels, + now, + comment, + vs.str(), + mesh.refinement_levels + ); + return header; + } + + std::string format_primary_mesh(const StroidMesh& mesh) { + std::stringstream ss; + ss.precision(8); + mesh.mesh->Print(ss); + + std::string pmesh = std::format("BEGIN BLOCK PMESH\n{}END BLOCK PMESH", ss.str()); + + return pmesh; + } + + template + std::string format_opt(const std::optional opt, T default_val) { + if (opt.has_value()) { + return std::format("{}", opt.value()); + } + return std::format("{}", default_val); + + } + + std::string format_reference_mesh(const StroidMesh& mesh) { + std::stringstream ss; + ss.precision(8); + mesh.reference_mesh->Print(ss); + + std::string rmesh = std::format("BEGIN BLOCK RMESH\n{}END BLOCK RMESH", ss.str()); + return rmesh; + } + + std::string format_config(const StroidMesh& mesh) { + config::MeshConfig d; + + config::OptimizationMethods d_opt = d.optimization_methods.value_or(config::OptimizationMethods{false, true}); + config::OptimizationMethods m_opt = mesh.config.optimization_methods.value_or(d_opt); + + std::string config_str = std::format(R"(BEGIN BLOCK CONFIG +# refiniment_levels: Initial number of levels of refinmenet, note the value in the header may be more up to date +# std::optional +# default: 4 +refinement_levels:{} + +# order: Polynomial / geometric order to use when constructing the mesh +# std::optional +# default: 3 +order:{} + +# include_external_domain: Whether or not to include the external domain in the mesh generally used for applying boundary conditions at infinity +# std::optional +# default: true +include_external_domain:{} + +# r_core: the radius of the stellar core region (in reference space) +# std::optional +# default: 0.25 +r_core:{} + +# r_star: the radius of the stellar surface (in reference space) +# std::optional +# default: 1.0 +r_star:{} + +# flattening: the flattening of the star (in reference space) where 0 is spherical and >0 is oblate. Note that this parameter is not equivalent to solving for the structure of a rotating model +# std::optional +# default: 0.0 +flattening:{} + +# r_infinity: the radius of the outer boundary of the mesh (in reference space) +# std::optional +# default: 6.0 +r_infinity:{} + +# r_instability: the radius inside which computations of geometry are skipped to avoid a core singularity +# std::optional +# default: 1e-14 +r_instability:{} + +# core_steepness: Controls the rate of transition of the core-to-envelope transition +# std::optional +# default: 1.0 +core_steepness:{} + +# continuity_order: order of continuity to force from teh core-envelope transition (0 = discontinuous, 1=C1 continuity, etc...) +# std::optional +# default: 2 +continuity_order:{} + +# surface_bdr_id: the boundary id to tag the stellar surface boundary elements as +# std::optional +# default: 1 +surface_bdr_id:{} + +# inf_bdr_id: the boundary id to tag the outer boundary elements as +# std::optional +# default: 2 +inf_bdr_id:{} + +# core_id: the material attribute to tag elements in the core region as +# std::optional +# default 1 +core_id:{} + +# envelope_id: the material attribute to tag elements in the envelope as +# std::optional +# default 2 +envelope_id:{} + +# vacuum_id: the material attribute to tag elements in the vacuum region as +# std::optional +# default 3 +vacuum_id:{} + +# optimization_method: struct for storing which optimization methods are being used +# includes tmop and smoothstep booleans +optimization_methods-tmop:{} +optimization_methods-smoothstep:{} +END BLOCK CONFIG)", + format_opt(mesh.config.refinement_levels, d.refinement_levels.value()), + format_opt(mesh.config.order, d.order.value()), + format_opt(mesh.config.include_external_domain, d.include_external_domain.value()), + format_opt(mesh.config.r_core, d.r_core.value()), + format_opt(mesh.config.r_star, d.r_star.value()), + format_opt(mesh.config.flattening, d.flattening.value()), + format_opt(mesh.config.r_infinity, d.r_infinity.value()), + format_opt(mesh.config.r_instability, d.r_instability.value()), + format_opt(mesh.config.core_steepness, d.core_steepness.value()), + format_opt(mesh.config.continuity_order, d.continuity_order.value()), + format_opt(mesh.config.surface_bdr_id, d.surface_bdr_id.value()), + format_opt(mesh.config.inf_bdr_id, d.inf_bdr_id.value()), + format_opt(mesh.config.core_id, d.core_id.value()), + format_opt(mesh.config.envelope_id, d.envelope_id.value()), + format_opt(mesh.config.vacuum_id, d.vacuum_id.value()), + m_opt.tmop.value_or(false), + m_opt.smoothstep.value_or(true)); + + return config_str; + } + } + + namespace { + constexpr std::string_view BEGIN_PREFIX = "BEGIN BLOCK "; + constexpr std::string_view END_PREFIX = "END BLOCK "; + + std::string_view trim(std::string_view s) { + const auto b = s.find_first_not_of(" \t\r\n"); + if (b == std::string_view::npos) return {}; + const auto e = s.find_last_not_of(" \t\r\n"); + return s.substr(b, e - b + 1); + } + + + std::expected parse_bool(std::string_view v) { + std::string s(trim(v)); + std::ranges::transform(s, s.begin(), + [](const unsigned char c) { return static_cast(std::tolower(c)); }); + if (s == "true" || s == "1") return true; + if (s == "false" || s == "0") return false; + return std::unexpected(std::format("invalid bool value '{}'", v)); + } + + template + std::expected parse_int(std::string_view v) { + const std::string_view s = trim(v); + T out{}; + if (const auto res = std::from_chars(s.data(), s.data() + s.size(), out); res.ec != std::errc{} || res.ptr != s.data() + s.size()) + return std::unexpected(std::format("invalid integer value '{}'", v)); + return out; + } + + std::expected parse_double(std::string_view v) { + const std::string_view s = trim(v); + double out{}; + if (const auto [ptr, ec] = std::from_chars(s.data(), s.data() + s.size(), out); ec != std::errc{} || ptr != s.data() + s.size()) + return std::unexpected(std::format("invalid floating-point value '{}'", v)); + return out; + } + + std::expected, std::string> extract_blocks(std::istream& is) { + std::map blocks; + std::string line; + std::string current; + std::string buffer; + bool in_block = false; + + while (std::getline(is, line)) { + const std::string_view t = trim(line); + + if (!in_block) { + if (t.starts_with(BEGIN_PREFIX)) { + current = std::string(trim(t.substr(BEGIN_PREFIX.size()))); + if (current.empty()) + return std::unexpected("found 'BEGIN BLOCK' with no block name"); + if (blocks.contains(current)) + return std::unexpected(std::format("duplicate block '{}'", current)); + buffer.clear(); + in_block = true; + } + } else { + if (t.starts_with(END_PREFIX)) { + if (const std::string end_name(trim(t.substr(END_PREFIX.size()))); end_name != current) + return std::unexpected(std::format( + "mismatched block markers: opened '{}' but closed '{}'", + current, end_name)); + blocks.emplace(std::move(current), std::move(buffer)); + current.clear(); + buffer.clear(); + in_block = false; + } else { + std::string_view raw = line; + if (!raw.empty() && raw.back() == '\r') raw.remove_suffix(1); + buffer.append(raw); + buffer.push_back('\n'); + } + } + } + + if (in_block) + return std::unexpected(std::format("unterminated block '{}' (missing END BLOCK)", current)); + return blocks; + } + + std::expected parse_header(const std::string& content, StroidMesh& out) { + std::istringstream iss(content); + std::string line; + std::optional type; + std::optional ref_levels; + + while (std::getline(iss, line)) { + const std::string_view t = trim(line); + if (t.empty() || t.starts_with('#')) continue; + + const auto colon = t.find(':'); + if (colon == std::string_view::npos) continue; + + const std::string_view key = trim(t.substr(0, colon)); + const std::string_view val = trim(t.substr(colon + 1)); + + if (key == "MESH_TYPE") { + if (val == "P") type = MFEM_MESH_TYPE::PARALLEL; + else if (val == "S") type = MFEM_MESH_TYPE::SERIAL; + else return std::unexpected(std::format("unknown MESH_TYPE '{}'", val)); + } else if (key == "REFINEMENT_LEVELS") { + auto r = parse_int(val); + if (!r) return std::unexpected("REFINEMENT_LEVELS: " + r.error()); + ref_levels = *r; + } + } + + if (!type) return std::unexpected("HEADER block missing MESH_TYPE"); + out.type = *type; + out.refinement_levels = ref_levels.value_or(0); + return {}; + } + + std::expected parse_config(const std::string& content) { + config::MeshConfig cfg; + config::OptimizationMethods opt = + cfg.optimization_methods.value_or(config::OptimizationMethods{}); + + using Handler = std::function(std::string_view)>; + + auto as_int = [](std::optional* f) { return [f](const std::string_view v) -> std::expected { auto r = parse_int(v); if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; }; + auto as_size = [](std::optional* f) { return [f](const std::string_view v) -> std::expected { auto r = parse_int(v); if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; }; + auto as_double = [](std::optional* f) { return [f](const std::string_view v) -> std::expected { auto r = parse_double(v); if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; }; + auto as_bool = [](std::optional* f) { return [f](const std::string_view v) -> std::expected { auto r = parse_bool(v); if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; }; + + const std::unordered_map handlers = { + {"refinement_levels", as_int(&cfg.refinement_levels)}, + {"order", as_int(&cfg.order)}, + {"include_external_domain", as_bool(&cfg.include_external_domain)}, + {"r_core", as_double(&cfg.r_core)}, + {"r_star", as_double(&cfg.r_star)}, + {"flattening", as_double(&cfg.flattening)}, + {"r_infinity", as_double(&cfg.r_infinity)}, + {"r_instability", as_double(&cfg.r_instability)}, + {"core_steepness", as_double(&cfg.core_steepness)}, + {"continuity_order", as_size(&cfg.continuity_order)}, + {"surface_bdr_id", as_size(&cfg.surface_bdr_id)}, + {"inf_bdr_id", as_size(&cfg.inf_bdr_id)}, + {"core_id", as_size(&cfg.core_id)}, + {"envelope_id", as_size(&cfg.envelope_id)}, + {"vacuum_id", as_size(&cfg.vacuum_id)}, + {"optimization_methods-tmop", as_bool(&opt.tmop)}, + {"optimization_methods-smoothstep", as_bool(&opt.smoothstep)}, + }; + + std::istringstream iss(content); + std::string line; + while (std::getline(iss, line)) { + const std::string_view t = trim(line); + if (t.empty() || t.starts_with('#')) continue; + + const auto colon = t.find(':'); + if (colon == std::string_view::npos) continue; + + const std::string_view key = trim(t.substr(0, colon)); + const std::string_view val = trim(t.substr(colon + 1)); + + const auto it = handlers.find(key); + if (it == handlers.end()) continue; + if (auto r = it->second(val); !r) + return std::unexpected(std::format("{}: {}", key, r.error())); + } + + cfg.optimization_methods = opt; + return cfg; + } + + std::expected, std::string> load_serial_mesh(const std::string& raw) { + if (trim(raw).empty()) return std::unexpected("empty mesh block"); + std::istringstream iss(raw); + try { + return std::make_unique(iss); + } catch (const std::exception& e) { + return std::unexpected(std::string("MFEM failed to parse mesh: ") + e.what()); + } + } + + struct ParsedMeta { + StroidMesh mesh; + std::string pmesh_raw; + std::string rmesh_raw; + }; + + std::expected parse_metadata(std::istream& is) { + auto blocks = extract_blocks(is); + if (!blocks) return std::unexpected(blocks.error()); + + auto need = [&](std::string_view name) -> std::expected { + const auto it = blocks->find(std::string(name)); + if (it == blocks->end()) + return std::unexpected(std::format("missing required block '{}'", name)); + return it->second; + }; + + ParsedMeta pm{}; + + const auto header = need("HEADER"); + if (!header) return std::unexpected(header.error()); + if (auto r = parse_header(*header, pm.mesh); !r) return std::unexpected(r.error()); + + const auto config = need("CONFIG"); + if (!config) return std::unexpected(config.error()); + auto cfg = parse_config(*config); + if (!cfg) return std::unexpected("CONFIG block -> " + cfg.error()); + pm.mesh.config = std::move(*cfg); + + const auto pmesh = need("PMESH"); + if (!pmesh) return std::unexpected(pmesh.error()); + pm.pmesh_raw = *pmesh; + + const auto rmesh = need("RMESH"); + if (!rmesh) return std::unexpected(rmesh.error()); + pm.rmesh_raw = *rmesh; + + return pm; + } + + } + + void SaveStroidMesh(const StroidMesh &mesh, const std::string &filename, const std::string &comment) { + std::ofstream ofs(filename); + + // First Write a header with some information + std::string header = format_header(mesh, comment); + + std::string pmesh = format_primary_mesh(mesh); + std::string rmesh = format_reference_mesh(mesh); + + std::string config = format_config(mesh); + + ofs << header << "\n"; + ofs << pmesh << "\n"; + ofs << rmesh << "\n"; + ofs << config << "\n"; + } void SaveMesh(const mfem::Mesh& mesh, const std::string& filename) { std::ofstream ofs(filename); @@ -15,6 +431,10 @@ namespace stroid::IO { mesh.Print(ofs); } + void SaveMesh(const stroid::StroidMesh &mesh, const std::string &filename) { + SaveMesh(*mesh.mesh, filename); + } + void SaveVTU(mfem::Mesh &mesh, const std::string &exportName) { mfem::ParaViewDataCollection pd(exportName, &mesh); pd.SetDataFormat(mfem::VTKFormat::BINARY); @@ -22,6 +442,10 @@ namespace stroid::IO { pd.Save(); } + void SaveVTU(const stroid::StroidMesh &mesh, const std::string &exportName) { + SaveVTU(*mesh.mesh, exportName); + } + void ViewMesh(mfem::Mesh &mesh, const std::string& title, const VISUALIZATION_MODE mode, const std::string &vishost, int visport) { mfem::socketstream sol_sock(vishost.c_str(), visport); if (!sol_sock.is_open()) { @@ -61,7 +485,12 @@ namespace stroid::IO { sol_sock << "keys iMj\n"; sol_sock << std::flush; } - void VisualizeFaceValence(mfem::Mesh& mesh) { + + void ViewMesh(const stroid::StroidMesh &mesh, const std::string &title, VISUALIZATION_MODE mode, const std::string &vishost, int visport) { + ViewMesh(*mesh.mesh, title, mode, vishost, visport); + } + + void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport) { mfem::L2_FECollection fec(0, 3); mfem::FiniteElementSpace fes(&mesh, &fec); mfem::GridFunction valence_gf(&fes); @@ -78,13 +507,81 @@ namespace stroid::IO { } // View in GLVis - char vishost[] = "localhost"; - int visport = 19916; - mfem::socketstream sol_sock(vishost, visport); + mfem::socketstream sol_sock(vishost.c_str(), visport); if (sol_sock.is_open()) { sol_sock << "solution\n" << mesh << valence_gf; sol_sock << "window_title 'Boundary Valence: 1=Surface, 2=Internal'\n"; sol_sock << "keys am\n" << std::flush; } } + + void VisualizeFaceValence(const stroid::StroidMesh &mesh, const std::string &vishost, int visport) { + VisualizeFaceValence(*mesh.mesh, vishost, visport); + } + + std::expected ParseStroidMesh(std::istream& is) { + auto pm = parse_metadata(is); + if (!pm) return std::unexpected(pm.error()); + + if (pm->mesh.type != MFEM_MESH_TYPE::SERIAL) { + return std::unexpected( + "parsed a PARALLEL StroidMesh, but ParseStroidMesh(std::istream&) can only " + "reconstruct serial meshes; use the MPI-aware overload " + "ParseStroidMesh(std::istream&, MPI_Comm) (requires MFEM_USE_MPI)"); + } + + auto m = load_serial_mesh(pm->pmesh_raw); + if (!m) return std::unexpected("PMESH -> " + m.error()); + auto rm = load_serial_mesh(pm->rmesh_raw); + if (!rm) return std::unexpected("RMESH -> " + rm.error()); + + pm->mesh.mesh = std::move(*m); + pm->mesh.reference_mesh = std::move(*rm); + return std::move(pm->mesh); + } + + std::expected LoadStroidMesh(const std::string& filename) { + std::ifstream ifs(filename); + if (!ifs.is_open()) + return std::unexpected(std::format("could not open file '{}'", filename)); + return ParseStroidMesh(ifs); + } + +#ifdef MFEM_USE_MPI + std::expected ParseStroidMesh(std::istream& is, MPI_Comm comm) { + auto pm = parse_metadata(is); + if (!pm) return std::unexpected(pm.error()); + + auto build = [&](const std::string& raw) + -> std::expected, std::string> { + if (trim(raw).empty()) return std::unexpected("empty mesh block"); + std::istringstream iss(raw); + try { + if (pm->mesh.type == MFEM_MESH_TYPE::PARALLEL) + return std::unique_ptr(new mfem::ParMesh(comm, iss)); + return std::make_unique(iss); + } catch (const std::exception& e) { + return std::unexpected(std::string("MFEM failed to parse mesh: ") + e.what()); + } + }; + + auto m = build(pm->pmesh_raw); + if (!m) return std::unexpected("PMESH -> " + m.error()); + auto rm = build(pm->rmesh_raw); + if (!rm) return std::unexpected("RMESH -> " + rm.error()); + + pm->mesh.mesh = std::move(*m); + pm->mesh.reference_mesh = std::move(*rm); + return std::move(pm->mesh); + } + + std::expected LoadStroidMesh(const std::string& filename, MPI_Comm comm) { + std::ifstream ifs(filename); + if (!ifs.is_open()) + return std::unexpected(std::format("could not open file '{}'", filename)); + return ParseStroidMesh(ifs, comm); + } +#endif // MFEM_USE_MPI + + } \ No newline at end of file diff --git a/src/lib/refinement/uniform.cpp b/src/lib/refinement/uniform.cpp new file mode 100644 index 0000000..2c51640 --- /dev/null +++ b/src/lib/refinement/uniform.cpp @@ -0,0 +1,39 @@ +#include "mfem.hpp" + +#include "stroid/refinement/uniform.h" +#include "stroid/utils/types.h" +#include "stroid/utils/mesh_utils.h" +#include "stroid/exceptions/exceptions.h" +#include "stroid/topology/topology.h" +#include "stroid/topology/optimize.h" + +namespace stroid::refinement { + void UniformRefinement(StroidMesh &mesh, const size_t levels) { + if (!mesh.reference_mesh) { + throw exceptions::StroidMissingReferenceMesh("UniformRefinement requires a reference mesh to be present in the StroidMesh object. This should be present by construction and the fact that is is missing represents a bug. Please report this to the stroid developers on GitHub or by email at emily.boudreaux@dartmouth.edu"); + } + + if (levels == 0) { + return; + } + + for (size_t i = 0; i < levels; i++) { + mesh.reference_mesh->UniformRefinement(); + } + + mesh.refinement_levels += levels; + + fourdst::config::Config cfg; + auto Mutator = [&mesh](config::MeshConfig& orig) { + orig = mesh.config; + }; + + cfg.mutate(Mutator); + + mesh.mesh = utils::BuildProjected(*mesh.reference_mesh, cfg); + + + + topology::OptimizeMesh(*mesh.mesh, cfg); + } +} diff --git a/src/lib/topology/mapping.cpp b/src/lib/topology/mapping.cpp index c0b42fd..87cda49 100644 --- a/src/lib/topology/mapping.cpp +++ b/src/lib/topology/mapping.cpp @@ -1,6 +1,64 @@ #include "stroid/topology/mapping.h" +#include "stroid/exceptions/exceptions.h" #include #include +#include +#include +#include +#include + +namespace { + template + consteval int nCr() { + if constexpr (k > n) { + return 0; + } else { + if constexpr (constexpr int kk = (k * 2 > n) ? (n - k) : k; kk == 0) { + return 1; + } else { + int result = n; + + for (int i = 2; i <= kk; ++i) { + result *= (n - i + 1); + result /= i; + } + + return result; + } + } + } + + template + double GeneralizedSmoothstep(const double x) { + if (x <= 0.0) return 0.0; + if (x >= 1.0) return 1.0; + + double sum = 0.0; + + auto compute_term = [&](std::integral_constant) { + return nCr() * std::pow(1.0 - x, k); + }; + + auto unroller = [&](std::index_sequence) { + return (compute_term(std::integral_constant{}) + ...); + }; + + sum = unroller(std::make_index_sequence{}); + return sum * std::pow(x, n + 1); + } + + template + constexpr auto make_smoothstep_dispatch_table(std::index_sequence) { + return std::array{ + &GeneralizedSmoothstep... + }; + } + + constexpr int MAX_SMOOTHSTEP_ORDER = 10; + constexpr auto smoothstep_dispatch = make_smoothstep_dispatch_table( + std::make_index_sequence{} + ); +} namespace stroid::topology { void ApplyEquiangular(mfem::Vector &pos) { @@ -33,51 +91,43 @@ namespace stroid::topology { } void TransformPoint(mfem::Vector &pos, const fourdst::config::Config &config, int attribute_id) { - double l_inf = 0.0; - for (int i = 0; i < pos.Size(); ++i) { - l_inf = std::max(l_inf, std::abs(pos(i))); - } + double X = pos(0); + double Y = pos(1); + double Z = pos(2); - if (l_inf < config->r_instability) return; + double maxAbs = std::max({std::abs(X), std::abs(Y), std::abs(Z)}); + if (maxAbs < 1e-14) return; - // Gnomonic projection - const double r_log = pos.Norml2(); - mfem::Vector unit_dir = pos; - unit_dir /= r_log; + double cx = X / maxAbs; + double cy = Y / maxAbs; + double cz = Z / maxAbs; - ApplyEquiangular(unit_dir); - unit_dir /= unit_dir.Norml2(); // Re-normalize + double sx = cx * std::sqrt(1.0 - cy*cy/2.0 - cz*cz/2.0 + cy*cy*cz*cz/3.0); + double sy = cy * std::sqrt(1.0 - cx*cx/2.0 - cz*cz/2.0 + cx*cx*cz*cz/3.0); + double sz = cz * std::sqrt(1.0 - cx*cx/2.0 - cy*cy/2.0 + cx*cx*cy*cy/3.0); - if (l_inf <= config->r_core) { - const double t = l_inf / config->r_core.value(); - double alpha = std::pow(t, config->core_steepness.value()); + mfem::Vector unit_dir(3); + unit_dir(0) = sx; + unit_dir(1) = sy; + unit_dir(2) = sz; - // Smoothstep function to apply C1 continuity - alpha = alpha * alpha * (3.0 - 2.0 * alpha); + if (maxAbs <= config->r_core.value()) { + double nx = X / config->r_core.value(); + double ny = Y / config->r_core.value(); + double nz = Z / config->r_core.value(); - mfem::Vector pos_cartesian = pos; - mfem::Vector pos_spherical = unit_dir; + pos(0) = nx * std::sqrt(1.0 - ny*ny/2.0 - nz*nz/2.0 + ny*ny*nz*nz/3.0); + pos(1) = ny * std::sqrt(1.0 - nx*nx/2.0 - nz*nz/2.0 + nx*nx*nz*nz/3.0); + pos(2) = nz * std::sqrt(1.0 - nx*nx/2.0 - ny*ny/2.0 + nx*nx*ny*ny/3.0); - pos_spherical *= l_inf; - bool run_smoothstep = false; - - - if (config->optimization_methods.has_value() && config->optimization_methods.value().smoothstep.has_value() && config->optimization_methods.value().smoothstep.value()) { - run_smoothstep = true; - } - - - if (run_smoothstep) { - for (int d = 0; d < pos.Size(); ++d) { - pos(d) = (1.0 - alpha) * pos_cartesian(d) + alpha * pos_spherical(d); - } - } + pos *= config->r_core.value(); ApplySpheroidal(pos, config); return; } - if (l_inf <= config->r_star) { - const double xi = (l_inf - config->r_core.value()) / (config->r_star.value() - config->r_core.value()); + + if (maxAbs <= config->r_star.value()) { + const double xi = (maxAbs - config->r_core.value()) / (config->r_star.value() - config->r_core.value()); const double r_phys = config->r_core.value() + xi * (config->r_star.value() - config->r_core.value()); pos = unit_dir; @@ -86,8 +136,74 @@ namespace stroid::topology { ApplySpheroidal(pos, config); } else { pos = unit_dir; - pos *= l_inf; + pos *= maxAbs; ApplySpheroidal(pos, config); } - }} + } + + // void TransformPoint(mfem::Vector &pos, const fourdst::config::Config &config, int attribute_id) { + // double l_inf = 0.0; + // for (int i = 0; i < pos.Size(); ++i) { + // l_inf = std::max(l_inf, std::abs(pos(i))); + // } + // + // if (l_inf < config->r_instability) return; + // + // // Gnomonic projection + // const double r_log = pos.Norml2(); + // mfem::Vector unit_dir = pos; + // unit_dir /= r_log; + // + // ApplyEquiangular(unit_dir); + // unit_dir /= unit_dir.Norml2(); // Re-normalize + // + // if (l_inf <= config->r_core) { + // const double t = l_inf / config->r_core.value(); + // double alpha = std::pow(t, config->core_steepness.value()); + // const size_t order = config->continuity_order.value_or(2); + // if (order < 1 || order > MAX_SMOOTHSTEP_ORDER) { + // const std::string err_msg = std::format("Invalid continuity order: {}. Continuity order must be between (inclusive) 1 and {}. To push to higher orders you must update MAX_SMOOTHSTEP_ORDER in src/lib/topology/mapping.cpp and recompile.", order, MAX_SMOOTHSTEP_ORDER); + // throw exceptions::StroidContinuityError(err_msg); + // } + // + // alpha = smoothstep_dispatch[order - 1](alpha); // We use this funky method as it keeps smoothstep calculation largely offloaded to compile time rather than run-time + // + // mfem::Vector pos_cartesian = pos; + // mfem::Vector pos_spherical = unit_dir; + // + // pos_spherical *= l_inf; + // bool run_smoothstep = false; + // + // + // if (config->optimization_methods.has_value() && config->optimization_methods.value().smoothstep.has_value() && config->optimization_methods.value().smoothstep.value()) { + // run_smoothstep = true; + // } + // + // + // if (run_smoothstep) { + // for (int d = 0; d < pos.Size(); ++d) { + // pos(d) = (1.0 - alpha) * pos_cartesian(d) + alpha * pos_spherical(d); + // } + // } + // + // ApplySpheroidal(pos, config); + // return; + // } + // + // if (l_inf <= config->r_star) { + // const double xi = (l_inf - config->r_core.value()) / (config->r_star.value() - config->r_core.value()); + // const double r_phys = config->r_core.value() + xi * (config->r_star.value() - config->r_core.value()); + // + // pos = unit_dir; + // pos *= r_phys; + // + // ApplySpheroidal(pos, config); + // } else { + // pos = unit_dir; + // pos *= l_inf; + // + // ApplySpheroidal(pos, config); + // } + // } +} diff --git a/src/lib/topology/optimize.cpp b/src/lib/topology/optimize.cpp index 8d5a66b..648376b 100644 --- a/src/lib/topology/optimize.cpp +++ b/src/lib/topology/optimize.cpp @@ -172,7 +172,8 @@ class TMOPProgressBar : public mfem::IterativeSolverMonitor { fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list); mfem::TMOP_QualityMetric* metric = new mfem::TMOP_Metric_302(); - mfem::TargetConstructor* target_c = new mfem::TargetConstructor(mfem::TargetConstructor::IDEAL_SHAPE_UNIT_SIZE); + mfem::TargetConstructor* target_c = new mfem::TargetConstructor(mfem::TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE); + target_c->SetNodes(*mesh.GetNodes()); mfem::TMOP_Integrator* tmop_integrator = new mfem::TMOP_Integrator(metric, target_c); mfem::NonlinearForm a(fes); @@ -185,10 +186,10 @@ class TMOPProgressBar : public mfem::IterativeSolverMonitor { b = 0.0; mfem::MINRESSolver minres; - minres.SetMaxIter(500); + minres.SetMaxIter(750); minres.SetRelTol(1e-5); minres.SetAbsTol(0.0); - minres.SetPrintLevel(0); + minres.SetPrintLevel(-1); mfem::DSmoother jacobi(1, 1.0, 1); jacobi.SetPositiveDiagonal(true); @@ -206,7 +207,7 @@ class TMOPProgressBar : public mfem::IterativeSolverMonitor { } } - constexpr double newton_rtol = 1e-4; + constexpr double newton_rtol = 1e-8; mfem::TMOPNewtonSolver newton(ir, 0); newton.SetPreconditioner(minres); newton.SetOperator(a); diff --git a/src/lib/utils/mesh_stats.cpp b/src/lib/utils/mesh_stats.cpp new file mode 100644 index 0000000..3c90b0a --- /dev/null +++ b/src/lib/utils/mesh_stats.cpp @@ -0,0 +1,485 @@ +#include "stroid/utils/mesh_stats.h" + +#include +#include +#include +#include +#include + +namespace stroid::stats { + namespace { + double SpheroidRadius(const double ux, const double uy, const double uz, + const double r_star, const double flattening) { + const double a = r_star; + const double c = r_star * (1.0 - flattening); + const double inv = (ux*ux + uy*uy) / (a*a) + (uz*uz) / (c*c); + return (inv > 0.0) ? 1.0 / std::sqrt(inv) : 0.0; + } + } + + MeshStats ComputeMeshStats(const StroidMesh& sm, MeshStatFeatures features, int sample_order) { + MeshStats out; + out.computed = features; + + auto mesh_or = sm.as_mesh(); + if (!mesh_or) { + out.errors.push_back(mesh_or.error()); + return out; + } + mfem::Mesh* mesh = *mesh_or; + if (!mesh) { + out.errors.push_back("StroidMesh has no stored computational mesh to compute stats against"); + return out; // BUGFIX: was falling through to a null deref + } + + const auto& cfg = sm.config; + const double r_star = cfg.r_star.value_or(1.0); + const double flattening = cfg.flattening.value_or(0.0); + const int surf_bdr = static_cast(cfg.surface_bdr_id.value_or(-99)); + const int inf_bdr = static_cast(cfg.inf_bdr_id.value_or(-99)); + const int core_id = static_cast(cfg.core_id.value_or(-99)); + const int env_id = static_cast(cfg.envelope_id.value_or(-99)); + const int vac_id = static_cast(cfg.vacuum_id.value_or(-99)); + + int geom_order = -99; + if (mesh->GetNodes()) { + geom_order = mesh->GetNodes()->FESpace()->GetMaxElementOrder(); + } + const int sorder = (sample_order > 0) ? sample_order : (2 * geom_order + 4); + + if (has_feature(features, MeshStatFeatures::CONFIG_META)) { + ConfigMeta meta; + meta.r_core = cfg.r_core.value_or(-99.99); + meta.r_star = cfg.r_star.value_or(-99.99); + meta.r_infinity = cfg.r_infinity.value_or(-99.99); + meta.flattening = flattening; + meta.geom_order = geom_order; + meta.refinement_levels = sm.refinement_levels; + meta.has_external_domain = cfg.include_external_domain.value_or(false); + out.config_meta = meta; + } + + if (has_feature(features, MeshStatFeatures::CONFORMITY)) { + ConformityStats conformity; + conformity.conforming = mesh->Conforming(); + conformity.n_nonconforming_faces = conformity.conforming ? 0 : -99; // TODO: count + out.conformity = conformity; + } + + // ============================ SURFACE PASS ============================ + const bool needs_surface = + has_feature(features, MeshStatFeatures::RADIUS) || + has_feature(features, MeshStatFeatures::AXES) || + has_feature(features, MeshStatFeatures::ELLIPTICITY) || + has_feature(features, MeshStatFeatures::BOWING) || + has_feature(features, MeshStatFeatures::VOLUME_AREA); + + if (needs_surface) { + double r_min = std::numeric_limits::max(); + double r_max = std::numeric_limits::lowest(); + double r_sum = 0.0, r_sum_sq = 0.0; + double a_eq = 0.0, c_pol = 0.0; + double bow_in = 0.0, bow_out = 0.0, bow_sum_sq = 0.0, bow_worst_r = 0.0; + double bow_worst_mag = 0.0; + double area = 0.0; + long n_samples = 0; + mfem::Vector phys; + + for (int b = 0; b < mesh->GetNBE(); ++b) { + if (mesh->GetBdrAttribute(b) != surf_bdr) continue; + mfem::ElementTransformation* T = mesh->GetBdrElementTransformation(b); + const mfem::IntegrationRule& ir = mfem::IntRules.Get(T->GetGeometryType(), sorder); + + for (int q = 0; q < ir.GetNPoints(); ++q) { + const mfem::IntegrationPoint& ip = ir.IntPoint(q); + T->SetIntPoint(&ip); + T->Transform(ip, phys); + + const double x = phys(0), y = phys(1), z = phys(2); + const double r = phys.Norml2(); + const double rho = std::sqrt(x*x + y*y); + const double w = T->Weight() * ip.weight; + + r_min = std::min(r_min, r); + r_max = std::max(r_max, r); + r_sum += r; r_sum_sq += r*r; + a_eq = std::max(a_eq, rho); + c_pol = std::max(c_pol, std::abs(z)); + area += w; + ++n_samples; + + if (has_feature(features, MeshStatFeatures::BOWING) && r > 1e-14) { + const double rt = SpheroidRadius(x/r, y/r, z/r, r_star, flattening); + const double dev = r - rt; + bow_in = std::min(bow_in, dev); + bow_out = std::max(bow_out, dev); + bow_sum_sq += dev * dev; + if (std::abs(dev) > bow_worst_mag) { + bow_worst_mag = std::abs(dev); + bow_worst_r = r; + } + } + } + } + + if (n_samples == 0) { + out.warnings.push_back("No samples were collected from the surface boundary. Check that the " + "surface boundary ID is correct. Lacking a surface pass prevents the " + "following from reporting accurate results: " + "[RADIUS, AXES, ELLIPTICITY, BOWING, VOLUME_AREA]"); + } else { + if (has_feature(features, MeshStatFeatures::RADIUS)) { + const double mean = r_sum / n_samples; + const double var = std::max(0.0, r_sum_sq / n_samples - mean * mean); + out.radius = RadiusStats{ + .min = r_min, .max = r_max, .mean = mean, + .stddev = std::sqrt(var), .n_samples = n_samples + }; + } + if (has_feature(features, MeshStatFeatures::AXES)) { + out.axes = AxisStats{ .semi_major = a_eq, .semi_minor = c_pol }; + } + if (has_feature(features, MeshStatFeatures::ELLIPTICITY)) { + out.ellipticity = EllipticityStats{ + .flattening = (a_eq > 0) ? (a_eq - c_pol) / a_eq : 0.0, + .polar_equatorial = (a_eq > 0) ? c_pol / a_eq : 1.0, + .radius_uniformity = (r_max > 0) ? r_min / r_max : 1.0, + }; + } + if (has_feature(features, MeshStatFeatures::BOWING)) { + out.bowing = BowingStats{ + .max_inward = bow_in, .max_outward = bow_out, + .rms = std::sqrt(bow_sum_sq / n_samples), .worst_at_radius = bow_worst_r, + }; + } + if (has_feature(features, MeshStatFeatures::VOLUME_AREA)) { + const double a = r_star, c = r_star * (1.0 - flattening); + out.volume = VolumeAreaStats{ + .surface_area = area, + .analytic_area = (flattening == 0) ? 4.0 * M_PI * a * a : -99.99 + }; + } + } + } + + // ============================ VOLUME PASS ============================ + const bool needs_volume = + has_feature(features, MeshStatFeatures::VOLUME_AREA) || + has_feature(features, MeshStatFeatures::JACOBIAN) || + has_feature(features, MeshStatFeatures::ELEMENT_COUNT) || + has_feature(features, MeshStatFeatures::MESH_SIZE) || + has_feature(features, MeshStatFeatures::CENTROID) || + has_feature(features, MeshStatFeatures::BOUNDING_BOX); // BUGFIX: was a dangling ';' + + if (needs_volume) { + const bool need_bbox = has_feature(features, MeshStatFeatures::BOUNDING_BOX); + const bool need_jac = has_feature(features, MeshStatFeatures::JACOBIAN); + + // Bounding-box accumulators (seeded inverted so empty regions stay invalid). + double cxmin=+std::numeric_limits::max(), cxmax=-std::numeric_limits::max(); + double cymin=cxmin, cymax=cxmax, czmin=cxmin, czmax=cxmax; // core + double sxmin=cxmin, sxmax=cxmax, symin=cxmin, symax=cxmax, szmin=cxmin, szmax=cxmax; // stellar + double vxmin=cxmin, vxmax=cxmax, vymin=cxmin, vymax=cxmax, vzmin=cxmin, vzmax=cxmax; // vacuum + long n_core_box=0, n_stel_box=0, n_vac_box=0; + mfem::Vector bphys; + + // Per-region Jacobian accumulators. + struct JacAccum { + double detJ_min = std::numeric_limits::max(); + double detJ_max = -std::numeric_limits::max(); + double min_ratio = 1.0; + long n_flipped = 0; + long n_elem = 0; + double worst_ratio_r = -1.0; + double min_detJ_r = -1.0; + }; + JacAccum all_acc, stel_acc, vac_acc; + + auto jac_update = [](JacAccum& a, double dmin, double dmax, bool flip, double r) { + ++a.n_elem; + if (dmin < a.detJ_min) { a.detJ_min = dmin; a.min_detJ_r = r; } + if (dmax > a.detJ_max) a.detJ_max = dmax; + if (dmax > 1e-30) { + const double ratio = dmin / dmax; + if (ratio < a.min_ratio) { a.min_ratio = ratio; a.worst_ratio_r = r; } + } + if (flip) ++a.n_flipped; + }; + + double vol = 0.0, cx = 0.0, cy = 0.0, cz = 0.0; + double h_min = std::numeric_limits::max(); + double h_max = -std::numeric_limits::max(); + double h_sum = 0.0, h_sum_sq = 0.0; + long n_core = 0, n_env = 0, n_vac = 0, n_other = 0; + mfem::Vector phys; + + for (int e = 0; e < mesh->GetNE(); ++e) { + const int attr = mesh->GetAttribute(e); + if (attr == core_id) ++n_core; + else if (attr == env_id) ++n_env; + else if (attr == vac_id) ++n_vac; + else ++n_other; + + const bool stellar = (attr == core_id || attr == env_id); + + if (has_feature(features, MeshStatFeatures::MESH_SIZE) && stellar) { + const double h = mesh->GetElementSize(e); + h_min = std::min(h_min, h); + h_max = std::max(h_max, h); + h_sum += h; h_sum_sq += h*h; + } + + mfem::ElementTransformation* T = mesh->GetElementTransformation(e); + const mfem::IntegrationRule& ir = mfem::IntRules.Get(T->GetGeometryType(), sorder); + + double e_detmin = std::numeric_limits::max(); + double e_detmax = -std::numeric_limits::max(); + bool e_flip = false; + + for (int q = 0; q < ir.GetNPoints(); ++q) { + const mfem::IntegrationPoint& ip = ir.IntPoint(q); + T->SetIntPoint(&ip); + + if (need_bbox) { + T->Transform(ip, bphys); + const double X = bphys(0), Y = bphys(1), Z = bphys(2); + if (attr == core_id) { + cxmin=std::min(cxmin,X); cxmax=std::max(cxmax,X); + cymin=std::min(cymin,Y); cymax=std::max(cymax,Y); + czmin=std::min(czmin,Z); czmax=std::max(czmax,Z); + ++n_core_box; + } + if (stellar) { // stellar = core U envelope + sxmin=std::min(sxmin,X); sxmax=std::max(sxmax,X); + symin=std::min(symin,Y); symax=std::max(symax,Y); + szmin=std::min(szmin,Z); szmax=std::max(szmax,Z); + ++n_stel_box; + } + if (attr == vac_id) { + vxmin=std::min(vxmin,X); vxmax=std::max(vxmax,X); + vymin=std::min(vymin,Y); vymax=std::max(vymax,Y); + vzmin=std::min(vzmin,Z); vzmax=std::max(vzmax,Z); + ++n_vac_box; + } + } + + const double dJ = T->Jacobian().Det(); + e_detmin = std::min(e_detmin, dJ); + e_detmax = std::max(e_detmax, dJ); + if (dJ < 0.0) e_flip = true; + + if (stellar && (has_feature(features, MeshStatFeatures::VOLUME_AREA) || + has_feature(features, MeshStatFeatures::CENTROID))) { + const double w = std::abs(dJ) * ip.weight; + vol += w; + if (has_feature(features, MeshStatFeatures::CENTROID)) { + T->Transform(ip, phys); + cx += w*phys(0); cy += w*phys(1); cz += w*phys(2); + } + } + } + + if (need_jac) { + // Representative element radius (center) for locating the worst element. + const mfem::IntegrationPoint& cip = mfem::Geometries.GetCenter(T->GetGeometryType()); + T->SetIntPoint(&cip); + mfem::Vector cpt; + T->Transform(cip, cpt); + const double er = cpt.Norml2(); + + jac_update(all_acc, e_detmin, e_detmax, e_flip, er); + if (stellar) jac_update(stel_acc, e_detmin, e_detmax, e_flip, er); + else if (attr == vac_id) jac_update(vac_acc, e_detmin, e_detmax, e_flip, er); + } + } + + if (has_feature(features, MeshStatFeatures::ELEMENT_COUNT)) { + out.element_counts = ElementCounts{ + .total = n_core + n_env + n_vac + n_other, + .core = n_core, .envelope = n_env, .vacuum = n_vac, .other = n_other, + .n_vertices = mesh->GetNV(), + }; + } + + if (need_jac) { + auto finalize = [](const JacAccum& a) { + JacobianStats j; + j.n_elements = a.n_elem; + j.detJ_min = (a.n_elem > 0) ? a.detJ_min : 0.0; + j.detJ_max = (a.n_elem > 0) ? a.detJ_max : 0.0; + j.min_detJ_ratio = a.min_ratio; + j.n_flipped = a.n_flipped; + j.worst_ratio_at_radius = a.worst_ratio_r; + j.detJ_min_at_radius = a.min_detJ_r; + return j; + }; + out.jacobian = finalize(all_acc); + if (stel_acc.n_elem > 0) out.jacobian_stellar = finalize(stel_acc); + if (vac_acc.n_elem > 0) out.jacobian_vacuum = finalize(vac_acc); + } + + if (has_feature(features, MeshStatFeatures::MESH_SIZE)) { + const long ns = n_core + n_env; + const double mean = (ns > 0) ? h_sum / ns : 0.0; + const double var = (ns > 0) ? std::max(0.0, h_sum_sq / ns - mean * mean) : 0.0; + out.mesh_size = MeshSizeStats{ + .h_min = h_min, .h_max = h_max, .h_mean = mean, .h_stddev = std::sqrt(var), + }; + } + + if (has_feature(features, MeshStatFeatures::VOLUME_AREA)) { + if (!out.volume) out.volume.emplace(); // BUGFIX: surface pass may not have created it + out.volume->stellar_volume = vol; + const double a = r_star, c = r_star * (1.0 - flattening); + out.volume->analytic_volume = (4.0 / 3.0) * M_PI * a * a * c; + } + + if (has_feature(features, MeshStatFeatures::CENTROID)) { + if (vol <= 0) { + out.warnings.push_back("Stellar volume is zero or negative, cannot compute centroid."); + } else { + const double ccx = cx / vol, ccy = cy / vol, ccz = cz / vol; // BUGFIX: normalize + out.centroid = CentroidStats{ + .x = ccx, .y = ccy, .z = ccz, + .offset = std::sqrt(ccx*ccx + ccy*ccy + ccz*ccz) + }; + } + } + + if (need_bbox) { + BoundingBoxStats bb; + auto fill = [](BoundingBox& box, long n, + double xmn,double xmx,double ymn,double ymx,double zmn,double zmx) { + if (n > 0) { + box.valid = true; + box.xMin=xmn; box.xMax=xmx; + box.yMin=ymn; box.yMax=ymx; + box.zMin=zmn; box.zMax=zmx; + } + }; + fill(bb.core, n_core_box, cxmin,cxmax,cymin,cymax,czmin,czmax); + fill(bb.star, n_stel_box, sxmin,sxmax,symin,symax,szmin,szmax); + fill(bb.vacuum, n_vac_box, vxmin,vxmax,vymin,vymax,vzmin,vzmax); + out.bounding_box = bb; + } + } + + // ============================ OUTER BOUND PASS ============================ + if (has_feature(features, MeshStatFeatures::OUTER_BOUNDS)) { + double r_min = std::numeric_limits::max(); + double r_max = std::numeric_limits::lowest(); + double r_sum = 0.0; + long n_samples = 0; + mfem::Vector phys; + + for (int b = 0; b < mesh->GetNBE(); ++b) { + if (mesh->GetBdrAttribute(b) != inf_bdr) continue; + mfem::ElementTransformation* T = mesh->GetBdrElementTransformation(b); + const mfem::IntegrationRule& ir = mfem::IntRules.Get(T->GetGeometryType(), sorder); + for (int q = 0; q < ir.GetNPoints(); ++q) { + T->SetIntPoint(&ir.IntPoint(q)); + T->Transform(ir.IntPoint(q), phys); + const double r = phys.Norml2(); + r_min = std::min(r_min, r); r_max = std::max(r_max, r); + r_sum += r; ++n_samples; + } + } + + if (n_samples == 0) { + out.warnings.push_back("No samples found on the outer boundary, cannot compute outer bounds."); + } else { + out.outer_bounds = OuterBoundsStats{ + .min = r_min, .max = r_max, .mean = r_sum / n_samples, .n_samples = n_samples + }; + } + } + + return out; + } + + std::string to_string(const MeshStats& s) { + std::string o = "MeshStats:\n"; + auto line = [&](const std::string& l){ o += " =>" + l + "\n"; }; + + if (s.config_meta) { + const auto& m = *s.config_meta; + line(std::format( + "config: r_core={:0.4f}, r_star={:0.4f}, r_inf={:0.4f}, flattening={:0.4f}, " + "geometric order={}, refinement levels={}", + m.r_core, m.r_star, m.r_infinity, m.flattening, m.geom_order, m.refinement_levels)); + } + if (s.radius) { + const auto& r = *s.radius; + line(std::format("radius: min={:.6f} max={:.6f} mean={:.6f} std={:.3E} (n={})", + r.min, r.max, r.mean, r.stddev, r.n_samples)); + } + if (s.axes) { + line(std::format("axes: semi_major(eq)={:.6f} semi_minor(pol)={:.6f}", + s.axes->semi_major, s.axes->semi_minor)); + } + if (s.ellipticity) { + const auto& e = *s.ellipticity; + line(std::format("ellipticity: flattening={:.5f} c/a={:.5f} r_min/r_max={:.5f}", + e.flattening, e.polar_equatorial, e.radius_uniformity)); + } + if (s.bowing) { + const auto& b = *s.bowing; + line(std::format("bowing: max_inward={:.3E} max_outward={:.3E} rms={:.3E}", + b.max_inward, b.max_outward, b.rms)); + } + if (s.conformity) { + line(std::format("conforming: {}", s.conformity->conforming)); + } + + auto jac_line = [&](const std::string& label, const JacobianStats& j) { + line(std::format( + "jacobian[{}]: detJ=[{:.3E},{:.3E}] min_ratio={:.3E} (@r={:.4f}) " + "min_detJ@r={:.4f} flipped={} n={}", + label, j.detJ_min, j.detJ_max, j.min_detJ_ratio, j.worst_ratio_at_radius, + j.detJ_min_at_radius, j.n_flipped, j.n_elements)); + }; + if (s.jacobian) jac_line("all", *s.jacobian); + if (s.jacobian_stellar) jac_line("stellar", *s.jacobian_stellar); + if (s.jacobian_vacuum) jac_line("vacuum", *s.jacobian_vacuum); + + if (s.volume) { + const auto& v = *s.volume; + line(std::format("volume={:.6f} (analytic {:.6f}) area={:.6f}", + v.stellar_volume, v.analytic_volume, v.surface_area)); + } + if (s.element_counts) { + const auto& c = *s.element_counts; + line(std::format("elements: total={} core={} env={} vac={} other={} NV={}", + c.total, c.core, c.envelope, c.vacuum, c.other, c.n_vertices)); + } + if (s.mesh_size) { + line(std::format("h: min={:.4E} max={:.4E} mean={:.4E} std={:.4E}", + s.mesh_size->h_min, s.mesh_size->h_max, s.mesh_size->h_mean, s.mesh_size->h_stddev)); + } + if (s.bounding_box) { + const auto& bb = *s.bounding_box; + auto bline = [&](const char* nm, const BoundingBox& b){ + if (b.valid) + line(std::format("bbox[{}]: x[{:.4f},{:.4f}] y[{:.4f},{:.4f}] z[{:.4f},{:.4f}]", + nm, b.xMin,b.xMax, b.yMin,b.yMax, b.zMin,b.zMax)); + else + line(std::format("bbox[{}]: ", nm)); + }; + bline("core", bb.core); + bline("star", bb.star); + bline("vacuum", bb.vacuum); + } + if (s.outer_bounds) { + line(std::format("outer: min={:.4f} max={:.4f} mean={:.4f}", + s.outer_bounds->min, s.outer_bounds->max, s.outer_bounds->mean)); + } + if (s.centroid) { + line(std::format("centroid: x={:.6f} y={:.6f} z={:.6f} offset={:.6f}", + s.centroid->x, s.centroid->y, s.centroid->z, s.centroid->offset)); + } + + for (const auto& w : s.warnings) line("WARNING: " + w); + for (const auto& e : s.errors) line(std::format("ERROR: {}", e)); + return o; + } +} \ No newline at end of file diff --git a/src/lib/utils/mesh_utils.cpp b/src/lib/utils/mesh_utils.cpp index 9c691f2..f2faf9d 100644 --- a/src/lib/utils/mesh_utils.cpp +++ b/src/lib/utils/mesh_utils.cpp @@ -2,6 +2,8 @@ #include "mfem.hpp" #include +#include "stroid/topology/curvilinear.h" + namespace stroid::utils { void MarkFlippedElements(mfem::Mesh& mesh) { for (int i = 0; i < mesh.GetNE(); i++) { @@ -51,4 +53,44 @@ namespace stroid::utils { } } } -} \ No newline at end of file + + void ExportJacobianRadialProfile(mfem::Mesh& mesh, const std::string& filename) { + std::ofstream ofs(filename); + + if (!ofs.good()) { + throw std::runtime_error(std::format("Stroid: Could not open file {} for writing Jacobian radial profile", filename)); + } + + ofs << "Radius,DetJ,Attribute,ElementID\n"; + ofs.precision(10); + + const int sample_order = 2 * mesh.GetNodes()->FESpace()->GetMaxElementOrder() + 2; + for (int i = 0; i < mesh.GetNE(); ++i) { + mfem::ElementTransformation *T = mesh.GetElementTransformation(i); + const int attr = mesh.GetAttribute(i); + + const mfem::IntegrationRule &ir = mfem::IntRules.Get(T->GetGeometryType(), sample_order); + + for (int j = 0; j < ir.GetNPoints(); ++j) { + T->SetIntPoint(&ir.IntPoint(j)); + + mfem::Vector pos; + T->Transform(ir.IntPoint(j), pos); + + const double r = pos.Norml2(); + const double detJ = T->Jacobian().Det(); + + ofs << r << "," << detJ << "," << attr << ',' << i << "\n"; + } + } + ofs.close(); + std::println("Jacobian radial profile exported to {}", filename); + } + + std::unique_ptr BuildProjected(const mfem::Mesh& reference, const fourdst::config::Config& cfg) { + auto projected = std::make_unique(reference); + topology::PromoteToHighOrder(*projected, cfg); + topology::ProjectMesh(*projected, cfg); + return projected; + } +} diff --git a/src/meson.build b/src/meson.build index e50ec57..99d9587 100644 --- a/src/meson.build +++ b/src/meson.build @@ -13,21 +13,73 @@ stroid_sources = files( 'lib/topology/optimize.cpp', 'lib/IO/mesh.cpp', 'lib/utils/mesh_utils.cpp', + 'lib/utils/mesh_stats.cpp', + 'lib/refinement/uniform.cpp', ) -stroid_lib = static_library( - 'libstroid', - stroid_sources, - include_directories: stroid_include_files, - dependencies: dependencies, - install: true -) +if get_option('build_python') + if host_machine.system() == 'darwin' + stroid_lib_rpath_args = [ + '-Wl,-rpath,@loader_path', + '-Wl,-rpath,@loader_path/../../fourdst/lib', + '-Wl,-rpath,@loader_path/../../fourdst/lib/vendor', + ] + stroid_lib_rpath = '' -stroid_dep = declare_dependency( - link_with: stroid_lib, - include_directories: stroid_include_files, - dependencies: dependencies -) + stroid_ext_rpath_args = [ + '-Wl,-rpath,@loader_path/lib', + '-Wl,-rpath,@loader_path/../fourdst/lib', + '-Wl,-rpath,@loader_path/../fourdst/lib/vendor', + ] + stroid_ext_rpath = '' + else + stroid_lib_rpath_args = [] + stroid_lib_rpath = '$ORIGIN:' + '$ORIGIN/../../fourdst/lib:' + '$ORIGIN/../../fourdst/lib/vendor' + + stroid_ext_rpath_args = [] + stroid_ext_rpath = '$ORIGIN/lib:' + '$ORIGIN/../fourdst/lib:' + '$ORIGIN/../fourdst/lib/vendor' + endif + + libstroid = static_library( + 'libstroid', + stroid_sources, + include_directories: stroid_include_files, + dependencies: dependencies, + install_dir: stroid_libdir, + link_args: stroid_lib_rpath_args, + build_rpath: stroid_lib_rpath, + install_rpath: stroid_lib_rpath, + ) +else + libstroid = static_library( + 'libstroid', + stroid_sources, + include_directories: stroid_include_files, + dependencies: dependencies, + install: true, + ) +endif + + +if get_option('build_python') + stroid_iface_dep = declare_dependency( + dependencies: dependencies + ).partial_dependency(compile_args: true, includes: true) + + stroid_dep = declare_dependency( + link_with: libstroid, + include_directories: stroid_include_files, + dependencies: [stroid_iface_dep] + ) +else + stroid_dep = declare_dependency( + link_with: libstroid, + include_directories: stroid_include_files, + dependencies: dependencies + ) +endif + +meson.override_dependency('stroid', stroid_dep) install_subdir( 'include/stroid', diff --git a/src/python/IO/bindings.cpp b/src/python/IO/bindings.cpp new file mode 100644 index 0000000..11108bb --- /dev/null +++ b/src/python/IO/bindings.cpp @@ -0,0 +1,77 @@ +#include +#include +#include "bindings.h" + +#include "stroid/IO/mesh.h" + +namespace py = pybind11; + +void register_io_bindings(pybind11::module_ &m) { + py::enum_(m, "VISUALIZATION_MODE") + .value("NONE", stroid::IO::VISUALIZATION_MODE::NONE) + .value("ELEMENT_ID", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID) + .value("BOUNDARY_ELEMENT_ID", stroid::IO::VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID) + .export_values(); + + m.def( + "SaveStroidMesh", + &stroid::IO::SaveStroidMesh, + py::arg("mesh"), + py::arg("filename"), + py::arg("comment")="", + "Save a Stroid mesh to a file." + ); + m.def( + "SaveMesh", + py::overload_cast(&stroid::IO::SaveMesh), + py::arg("mesh"), + py::arg("filename") + ); + m.def( + "SaveVTU", + py::overload_cast(&stroid::IO::SaveVTU), + py::arg("mesh"), + py::arg("filename") + ); + m.def( + "ViewMesh", + py::overload_cast(&stroid::IO::ViewMesh), + py::arg("mesh"), + py::arg("title")="", + py::arg("mode")=stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, + py::arg("host")="localhost", + py::arg("port")=19916 + ); + + m.def( + "VisualizeFaceValence", + py::overload_cast(&stroid::IO::VisualizeFaceValence), + py::arg("mesh"), + py::arg("host")="localhost", + py::arg("port")=19916 + ); + + m.def( + "ParseStroidMesh", + [](const std::string& buf) { + std::stringstream ss; + ss << buf; + auto r = stroid::IO::ParseStroidMesh(ss); + if (!r.has_value()) { + throw std::runtime_error("Parsing failed: " + r.error()); + } + return std::move(r.value()); + } + ); + + m.def( + "LoadStroidMesh", + [](const std::string& filename) { + auto r = stroid::IO::LoadStroidMesh(filename); + if (!r.has_value()) { + throw std::runtime_error("Loading " + filename + " failed: " + r.error()); + } + return std::move(r.value()); + } + ); +} diff --git a/src/python/IO/bindings.h b/src/python/IO/bindings.h new file mode 100644 index 0000000..e0f42d7 --- /dev/null +++ b/src/python/IO/bindings.h @@ -0,0 +1,5 @@ +#pragma once + +#include + +void register_io_bindings(pybind11::module_& m); diff --git a/src/python/bindings.cpp b/src/python/bindings.cpp new file mode 100644 index 0000000..35706c7 --- /dev/null +++ b/src/python/bindings.cpp @@ -0,0 +1,36 @@ +#include +#include + +#include "config/bindings.h" +#include "exceptions/bindings.h" +#include "IO/bindings.h" +#include "refinement/bindings.h" +#include "utils/bindings.h" + +#include "stroid/exceptions/stroid_error.h" + +#include "stroid/stroid.h" +#include "stroid/version.h" + +PYBIND11_MODULE(_stroid, m) { + m.doc() = "Python bindings for stroid library."; + + register_utils_bindings(m); + + auto exceptionsMod = m.def_submodule("exceptions", "Exceptions Bindings"); + register_exceptions_bindings(exceptionsMod); + + auto configMod = m.def_submodule("config", "Config Bindings"); + register_config_bindings(configMod); + + auto IOMod = m.def_submodule("IO", "IO Bindings"); + register_io_bindings(IOMod); + + auto refinementMod = m.def_submodule("refinement", "Refinement Bindings"); + register_refinement_bindings(refinementMod); + + m.def("GenerateMesh", pybind11::overload_cast(&stroid::GenerateMesh), "Generate a mesh from a MeshConfig object."); + m.def("GenerateMesh", pybind11::overload_cast(&stroid::GenerateMesh), "Generate a mesh from a config file path."); + + +} diff --git a/src/python/config/bindings.cpp b/src/python/config/bindings.cpp new file mode 100644 index 0000000..ee0ca8f --- /dev/null +++ b/src/python/config/bindings.cpp @@ -0,0 +1,202 @@ +#include +#include +#include "bindings.h" + +#include "stroid/config/config.h" + +namespace py = pybind11; + +void register_config_bindings(pybind11::module_& m) { + py::class_(m, "OptimizationMethods") + .def(py::init([](bool tmop, bool smoothstep) { + return stroid::config::OptimizationMethods{tmop, smoothstep}; + }), py::arg("tmop") = false, py::arg("smoothstep") = true) + .def_property("tmop", + [](const stroid::config::OptimizationMethods& self) { + return self.tmop; + }, + [](stroid::config::OptimizationMethods& self, bool value) { + self.tmop = value; + } + ) + .def_property("smoothstep", + [](const stroid::config::OptimizationMethods& self) { + return self.smoothstep; + }, + [](stroid::config::OptimizationMethods& self, bool value) { + self.smoothstep = value; + } + ); + + py::class_(m, "MeshConfig") + .def(py::init([](py::kwargs kwargs) { + int ref_level = 4, order = 3; + size_t continuity_order = 2, surface_bdr_id = 1, inf_bdr_id = 2, core_id = 1, envelope_id = 2, vacuum_id=3; + bool include_external_domain = true; + double r_core = 0.25, r_star = 1.0, flattening = 0.0, r_inf = 6.0, r_instability = 1e-14, core_steepness = 1.0; + stroid::config::OptimizationMethods opt_method{.tmop = false, .smoothstep = true}; + + return stroid::config::MeshConfig{ + .refinement_levels = kwargs.contains("refinement_levels") ? kwargs["refinement_levels"].cast() : ref_level, + .order = kwargs.contains("order") ? kwargs["order"].cast() : order, + .include_external_domain = kwargs.contains("include_external_domain") ? kwargs["include_external_domain"].cast() : include_external_domain, + .r_core = kwargs.contains("r_core") ? kwargs["r_core"].cast() : r_core, + .r_star = kwargs.contains("r_star") ? kwargs["r_star"].cast() : r_star, + .flattening = kwargs.contains("flattening") ? kwargs["flattening"].cast() : flattening, + .r_infinity = kwargs.contains("r_infinity") ? kwargs["r_infinity"].cast() : r_inf, + .r_instability = kwargs.contains("r_instability") ? kwargs["r_instability"].cast() : r_instability, + .core_steepness = kwargs.contains("core_steepness") ? kwargs["core_steepness"].cast() : core_steepness, + .continuity_order = kwargs.contains("continuity_order") ? kwargs["continuity_order"].cast() : continuity_order, + .surface_bdr_id = kwargs.contains("surface_bdr_id") ? kwargs["surface_bdr_id"].cast() : surface_bdr_id, + .inf_bdr_id = kwargs.contains("inf_bdr_id") ? kwargs["inf_bdr_id"].cast() : inf_bdr_id, + .core_id = kwargs.contains("core_id") ? kwargs["core_id"].cast() : core_id, + .envelope_id = kwargs.contains("envelope_id") ? kwargs["envelope_id"].cast() : envelope_id, + .vacuum_id = kwargs.contains("vacuum_id") ? kwargs["vacuum_id"].cast() : vacuum_id, + .optimization_methods = kwargs.contains("optimization_methods") ? kwargs["optimization_methods"].cast() : opt_method + }; + })) + .def_property( + "refinement_levels", + [](const stroid::config::MeshConfig& self) { + return self.refinement_levels; + }, + [](stroid::config::MeshConfig& self, int value) { + self.refinement_levels = value; + } + ) + .def_property( + "order", + [](const stroid::config::MeshConfig& self) { + return self.order; + }, + [](stroid::config::MeshConfig& self, int value) { + self.order = value; + } + ) + .def_property( + "include_external_domain", + [](const stroid::config::MeshConfig& self) { + return self.include_external_domain; + }, + [](stroid::config::MeshConfig& self, bool value) { + self.include_external_domain = value; + } + ) + .def_property( + "r_core", + [](const stroid::config::MeshConfig& self) { + return self.r_core; + }, + [](stroid::config::MeshConfig& self, int value) { + self.order = value; + } + ) + .def_property( + "r_star", + [](const stroid::config::MeshConfig& self) { + return self.r_star; + }, + [](stroid::config::MeshConfig& self, double value) { + self.r_star = value; + } + ) + .def_property( + "flattening", + [](const stroid::config::MeshConfig& self) { + return self.flattening; + }, + [](stroid::config::MeshConfig& self, double value) { + self.flattening = value; + } + ) + .def_property( + "r_infinity", + [](const stroid::config::MeshConfig& self) { + return self.r_infinity; + }, + [](stroid::config::MeshConfig& self, double value) { + self.r_infinity = value; + } + ) + .def_property( + "r_instability", + [](const stroid::config::MeshConfig& self) { + return self.r_instability; + }, + [](stroid::config::MeshConfig& self, double value) { + self.r_instability = value; + } + ) + .def_property( + "core_steepness", + [](const stroid::config::MeshConfig& self) { + return self.core_steepness; + }, + [](stroid::config::MeshConfig& self, double value) { + self.core_steepness = value; + } + ) + .def_property( + "continuity_order", + [](const stroid::config::MeshConfig& self) { + return self.continuity_order; + }, + [](stroid::config::MeshConfig& self, size_t value) { + self.continuity_order = value; + } + ) + .def_property( + "surface_bdr_id", + [](const stroid::config::MeshConfig& self) { + return self.surface_bdr_id; + }, + [](stroid::config::MeshConfig& self, size_t value) { + self.surface_bdr_id = value; + } + ) + .def_property( + "inf_bdr_id", + [](const stroid::config::MeshConfig& self) { + return self.inf_bdr_id; + }, + [](stroid::config::MeshConfig& self, size_t value) { + self.inf_bdr_id = value; + } + ) + .def_property( + "core_id", + [](const stroid::config::MeshConfig& self) { + return self.core_id; + }, + [](stroid::config::MeshConfig& self, size_t value) { + self.core_id = value; + } + ) + .def_property( + "envelope_id", + [](const stroid::config::MeshConfig& self) { + return self.envelope_id; + }, + [](stroid::config::MeshConfig& self, size_t value) { + self.envelope_id = value; + } + ) + .def_property( + "vacuum_id", + [](const stroid::config::MeshConfig& self) { + return self.vacuum_id; + }, + [](stroid::config::MeshConfig& self, size_t value) { + self.vacuum_id = value; + } + ) + .def_property( + "optimization_methods", + [](const stroid::config::MeshConfig& self) { + return self.optimization_methods; + }, + [](stroid::config::MeshConfig& self, stroid::config::OptimizationMethods value) { + self.optimization_methods = value; + } + ); +} \ No newline at end of file diff --git a/src/python/config/bindings.h b/src/python/config/bindings.h new file mode 100644 index 0000000..9d40686 --- /dev/null +++ b/src/python/config/bindings.h @@ -0,0 +1,5 @@ +#pragma once + +#include + +void register_config_bindings(pybind11::module_& m); diff --git a/src/python/exceptions/bindings.cpp b/src/python/exceptions/bindings.cpp new file mode 100644 index 0000000..0998d18 --- /dev/null +++ b/src/python/exceptions/bindings.cpp @@ -0,0 +1,14 @@ +#include + +#include "bindings.h" +#include "stroid/exceptions/exceptions.h" + +namespace py = pybind11; + + +void register_exceptions_bindings(py::module_& m) { + py::register_exception(m, "StroidError"); + py::register_exception(m, "StroidContinuityError", m.attr("StroidError")); + py::register_exception(m, "StroidMeshError", m.attr("StroidError")); + py::register_exception(m, "StroidMissingReferenceMesh", m.attr("StroidMeshError")); +} \ No newline at end of file diff --git a/src/python/exceptions/bindings.h b/src/python/exceptions/bindings.h new file mode 100644 index 0000000..6dc1f0d --- /dev/null +++ b/src/python/exceptions/bindings.h @@ -0,0 +1,5 @@ +#pragma once + +#include + +void register_exceptions_bindings(pybind11::module_& m); \ No newline at end of file diff --git a/src/python/refinement/bindings.cpp b/src/python/refinement/bindings.cpp new file mode 100644 index 0000000..72a0e48 --- /dev/null +++ b/src/python/refinement/bindings.cpp @@ -0,0 +1,11 @@ +#include +#include +#include "bindings.h" + +#include "stroid/refinement/uniform.h" + +namespace py = pybind11; + +void register_refinement_bindings(pybind11::module_ &m) { + m.def("UniformRefinement", &stroid::refinement::UniformRefinement, py::arg("mesh"), py::arg("levels"), "Perform uniform refinement without breaking the higher order structure"); +} diff --git a/src/python/refinement/bindings.h b/src/python/refinement/bindings.h new file mode 100644 index 0000000..ee129bf --- /dev/null +++ b/src/python/refinement/bindings.h @@ -0,0 +1,5 @@ +#pragma once + +#include + +void register_refinement_bindings(pybind11::module_& m); diff --git a/src/python/stroid/__init__.py b/src/python/stroid/__init__.py new file mode 100644 index 0000000..266673d --- /dev/null +++ b/src/python/stroid/__init__.py @@ -0,0 +1,45 @@ +import io +import sys + +from ._stroid import * + +from ._stroid import config +from ._stroid import exceptions +from ._stroid import IO +from ._stroid import refinement +from ._stroid import stats +from ._stroid import GenerateMesh +from ._stroid import StroidMesh + +sys.modules['stroid.config'] = config +sys.modules['stroid.exceptions'] = exceptions +sys.modules['stroid.IO'] = IO +sys.modules["stroid.refinement"] = refinement +sys.modules["stroid.stats"] = stats + + +__all__ = ['config', 'exceptions', 'IO', 'refinement', 'stats', 'GenerateMesh', 'StroidMesh'] + +import importlib.metadata + + +try: + _meta = importlib.metadata.metadata('stroid') + __version__ = _meta['Version'] + __license__ = _meta['License'] + __description__ = _meta['Summary'] + __author__ = 'Emily M. Boudreaux' + __url__ = 'https://github.com/4D-STAR/stroid' +except importlib.metadata.PackageNotFoundError : + __version__ = 'unknown - Package not installed' + __license__ = 'GNU General Public License v3.0' + __email__ = 'emily.boudreaux@dartmouth.edu' + __url__ = 'https://github.com/4D-STAR/stroid' + + +import os +from pathlib import Path +from typing import List + +_PACKAGE_DIR = Path(__file__).resolve().parent + diff --git a/src/python/utils/bindings.cpp b/src/python/utils/bindings.cpp new file mode 100644 index 0000000..af47d51 --- /dev/null +++ b/src/python/utils/bindings.cpp @@ -0,0 +1,187 @@ +#include +#include +#include "bindings.h" + +#include "stroid/utils/types.h" +#include "stroid/utils/mesh_stats.h" +#include "stroid/utils/mesh_utils.h" + +namespace py = pybind11; + +void register_stats_bindings(pybind11::module_ &m) { + auto statsMod = m.def_submodule("stats", "Stats Bindings"); + py::enum_(statsMod, "MeshStatFeatures", py::arithmetic()) + .value("NONE", stroid::stats::MeshStatFeatures::NONE) + .value("RADIUS", stroid::stats::MeshStatFeatures::RADIUS) + .value("AXES", stroid::stats::MeshStatFeatures::AXES) + .value("ELLIPTICITY", stroid::stats::MeshStatFeatures::ELLIPTICITY) + .value("BOWING", stroid::stats::MeshStatFeatures::BOWING) + .value("CONFORMITY", stroid::stats::MeshStatFeatures::CONFORMITY) + .value("JACOBIAN", stroid::stats::MeshStatFeatures::JACOBIAN) + .value("VOLUME_AREA", stroid::stats::MeshStatFeatures::VOLUME_AREA) + .value("ELEMENT_COUNT", stroid::stats::MeshStatFeatures::ELEMENT_COUNT) + .value("MESH_SIZE", stroid::stats::MeshStatFeatures::MESH_SIZE) + .value("OUTER_BOUNDS", stroid::stats::MeshStatFeatures::OUTER_BOUNDS) + .value("CENTROID", stroid::stats::MeshStatFeatures::CENTROID) + .value("CONFIG_META", stroid::stats::MeshStatFeatures::CONFIG_META) + .value("BOUNDING_BOX", stroid::stats::MeshStatFeatures::BOUNDING_BOX) + .export_values(); + + py::class_(statsMod, "RadiusStats") + .def_readonly("min", &stroid::stats::RadiusStats::min) + .def_readonly("max", &stroid::stats::RadiusStats::max) + .def_readonly("mean", &stroid::stats::RadiusStats::mean) + .def_readonly("stddev", &stroid::stats::RadiusStats::stddev) + .def_readonly("n_samples", &stroid::stats::RadiusStats::n_samples); + + py::class_(statsMod, "AxisStats") + .def_readonly("semi_major", &stroid::stats::AxisStats::semi_major) + .def_readonly("semi_minor", &stroid::stats::AxisStats::semi_minor); + + py::class_(statsMod, "EllipticityStats") + .def_readonly("flattening", &stroid::stats::EllipticityStats::flattening) + .def_readonly("polar_equatorial", &stroid::stats::EllipticityStats::polar_equatorial) + .def_readonly("radius_uniformity", &stroid::stats::EllipticityStats::radius_uniformity); + + py::class_(statsMod, "BowingStats") + .def_readonly("max_inward", &stroid::stats::BowingStats::max_inward) + .def_readonly("max_outward", &stroid::stats::BowingStats::max_outward) + .def_readonly("rms", &stroid::stats::BowingStats::rms) + .def_readonly("worst_at_radius", &stroid::stats::BowingStats::worst_at_radius); + + py::class_(statsMod, "ConformityStats") + .def_readonly("conforming", &stroid::stats::ConformityStats::conforming) + .def_readonly("n_nonconforming_faces", &stroid::stats::ConformityStats::n_nonconforming_faces); + + py::class_(statsMod, "JacobianStats") + .def_readonly("detJ_min", &stroid::stats::JacobianStats::detJ_min) + .def_readonly("detJ_max", &stroid::stats::JacobianStats::detJ_max) + .def_readonly("n_flipped", &stroid::stats::JacobianStats::n_flipped) + .def_readonly("min_detJ_ratio", &stroid::stats::JacobianStats::min_detJ_ratio) + .def_readonly("worst_ratio_at_radius", &stroid::stats::JacobianStats::worst_ratio_at_radius) + .def_readonly("detJ_min_at_radius", &stroid::stats::JacobianStats::detJ_min_at_radius) + .def_readonly("n_elements", &stroid::stats::JacobianStats::n_elements); + + py::class_(statsMod, "VolumeAreaStats") + .def_readonly("stellar_volume", &stroid::stats::VolumeAreaStats::stellar_volume) + .def_readonly("surface_area", &stroid::stats::VolumeAreaStats::surface_area) + .def_readonly("analytic_volume", &stroid::stats::VolumeAreaStats::analytic_volume) + .def_readonly("analytic_area", &stroid::stats::VolumeAreaStats::analytic_area); + + py::class_(statsMod, "ElementCounts") + .def_readonly("total", &stroid::stats::ElementCounts::total) + .def_readonly("core", &stroid::stats::ElementCounts::core) + .def_readonly("envelope", &stroid::stats::ElementCounts::envelope) + .def_readonly("vacuum", &stroid::stats::ElementCounts::vacuum) + .def_readonly("other", &stroid::stats::ElementCounts::other) + .def_readonly("n_vertices", &stroid::stats::ElementCounts::n_vertices); + + py::class_(statsMod, "MeshSizeStats") + .def_readonly("h_min", &stroid::stats::MeshSizeStats::h_min) + .def_readonly("h_max", &stroid::stats::MeshSizeStats::h_max) + .def_readonly("h_mean", &stroid::stats::MeshSizeStats::h_mean) + .def_readonly("h_stddev", &stroid::stats::MeshSizeStats::h_stddev); + + py::class_(statsMod, "OuterBoundsStats") + .def_readonly("min", &stroid::stats::OuterBoundsStats::min) + .def_readonly("max", &stroid::stats::OuterBoundsStats::max) + .def_readonly("mean", &stroid::stats::OuterBoundsStats::mean) + .def_readonly("n_samples", &stroid::stats::OuterBoundsStats::n_samples); + + py::class_(statsMod, "CentroidStats") + .def_readonly("x", &stroid::stats::CentroidStats::x) + .def_readonly("y", &stroid::stats::CentroidStats::y) + .def_readonly("z", &stroid::stats::CentroidStats::z) + .def_readonly("offset", &stroid::stats::CentroidStats::offset); + + py::class_(statsMod, "ConfigMeta") + .def_readonly("r_core", &stroid::stats::ConfigMeta::r_core) + .def_readonly("r_star", &stroid::stats::ConfigMeta::r_star) + .def_readonly("flattening", &stroid::stats::ConfigMeta::flattening) + .def_readonly("r_infinity", &stroid::stats::ConfigMeta::r_infinity) + .def_readonly("geom_order", &stroid::stats::ConfigMeta::geom_order) + .def_readonly("refinement_levels", &stroid::stats::ConfigMeta::refinement_levels) + .def_readonly("has_external_domain", &stroid::stats::ConfigMeta::has_external_domain); + + py::class_(statsMod, "BoundingBox") + .def_readonly("xMin", &stroid::stats::BoundingBox::xMin) + .def_readonly("xMax", &stroid::stats::BoundingBox::xMax) + .def_readonly("yMin", &stroid::stats::BoundingBox::yMin) + .def_readonly("yMax", &stroid::stats::BoundingBox::yMax) + .def_readonly("zMin", &stroid::stats::BoundingBox::zMin) + .def_readonly("zMax", &stroid::stats::BoundingBox::zMax) + .def_readonly("valid", &stroid::stats::BoundingBox::valid) + .def("dx", &stroid::stats::BoundingBox::dx) + .def("dy", &stroid::stats::BoundingBox::dy) + .def("dz", &stroid::stats::BoundingBox::dz) + .def("diag", &stroid::stats::BoundingBox::diag); + + py::class_(statsMod, "BoundingBoxStats") + .def_readonly("core", &stroid::stats::BoundingBoxStats::core) + .def_readonly("star", &stroid::stats::BoundingBoxStats::star) + .def_readonly("vacuum", &stroid::stats::BoundingBoxStats::vacuum); + + py::class_(statsMod, "MeshStats") + .def_readonly("computed", &stroid::stats::MeshStats::computed) + .def_readonly("radius", &stroid::stats::MeshStats::radius) + .def_readonly("axes", &stroid::stats::MeshStats::axes) + .def_readonly("ellipticity", &stroid::stats::MeshStats::ellipticity) + .def_readonly("bowing", &stroid::stats::MeshStats::bowing) + .def_readonly("conformity", &stroid::stats::MeshStats::conformity) + .def_readonly("jacobian", &stroid::stats::MeshStats::jacobian) + .def_readonly("jacobian_stellar", &stroid::stats::MeshStats::jacobian_stellar) + .def_readonly("jacobian_vacuum", &stroid::stats::MeshStats::jacobian_vacuum) + .def_readonly("volume", &stroid::stats::MeshStats::volume) + .def_readonly("element_counts", &stroid::stats::MeshStats::element_counts) + .def_readonly("mesh_size", &stroid::stats::MeshStats::mesh_size) + .def_readonly("outer_bounds", &stroid::stats::MeshStats::outer_bounds) + .def_readonly("centroid", &stroid::stats::MeshStats::centroid) + .def_readonly("config_meta", &stroid::stats::MeshStats::config_meta) + .def_readonly("bounding_box", &stroid::stats::MeshStats::bounding_box) + .def_readonly("warnings", &stroid::stats::MeshStats::warnings) + .def_readonly("errors", &stroid::stats::MeshStats::errors) + .def("__repr__", [](const stroid::stats::MeshStats& self) { + return stroid::stats::to_string(self); + }); + + statsMod.attr("MESH_STAT_DEFAULT") = stroid::stats::MESH_STAT_DEFAULT; + statsMod.attr("MESH_STAT_ALL") = stroid::stats::MESH_STAT_ALL; + + statsMod.def( + "ComputeMeshStats", + &stroid::stats::ComputeMeshStats, + py::arg("mesh"), + py::arg("features") = stroid::stats::MESH_STAT_DEFAULT, + py::arg("sample_order")=-1 + ); +} + +void register_type_bindings(py::module_ &m) { + py::enum_(m, "MFEM_MESH_TYPE") + .value("SERIAL", stroid::MFEM_MESH_TYPE::SERIAL) + .value("PARALLEL", stroid::MFEM_MESH_TYPE::PARALLEL) + .export_values(); + + py::class_(m, "StroidMesh") + .def_property_readonly("type", [](const stroid::StroidMesh& self) { + return (self.type == stroid::MFEM_MESH_TYPE::SERIAL) ? "SERIAL" : "PARALLEL"; + }) + .def_readonly("config", &stroid::StroidMesh::config) + .def_readonly("refinement_levels", &stroid::StroidMesh::refinement_levels) + .def("has_mesh", [](const stroid::StroidMesh& self) { + return self.mesh != nullptr; + }) + .def("has_rmesh", [](const stroid::StroidMesh& self) { + return self.reference_mesh != nullptr; + }) + .def("mesh_stats", &stroid::StroidMesh::mesh_stats) + .def("__repr__", [](const stroid::StroidMesh& self) { + return std::format("", (self.type == stroid::MFEM_MESH_TYPE::SERIAL) ? "SERIAL" : "PARALLEL", self.mesh->GetNE(), self.mesh->GetNV()); + }); +} + +void register_utils_bindings(pybind11::module_ &m) { + register_type_bindings(m); + register_stats_bindings(m); +} + diff --git a/src/python/utils/bindings.h b/src/python/utils/bindings.h new file mode 100644 index 0000000..250b767 --- /dev/null +++ b/src/python/utils/bindings.h @@ -0,0 +1,5 @@ +#pragma once + +#include + +void register_utils_bindings(pybind11::module_& m); diff --git a/subprojects/packagefiles/pybind11/LICENSE.build b/subprojects/packagefiles/pybind11/LICENSE.build new file mode 100644 index 0000000..4c99270 --- /dev/null +++ b/subprojects/packagefiles/pybind11/LICENSE.build @@ -0,0 +1,19 @@ +Copyright (c) 2021 The Meson development team + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. \ No newline at end of file diff --git a/subprojects/packagefiles/pybind11/meson.build b/subprojects/packagefiles/pybind11/meson.build new file mode 100644 index 0000000..0cc9de8 --- /dev/null +++ b/subprojects/packagefiles/pybind11/meson.build @@ -0,0 +1,8 @@ +project('pybind11', 'cpp', + version : 'v3.0.0', + license : 'BSD-3-Clause') + +pybind11_incdir = include_directories('include') + +pybind11_dep = declare_dependency( + include_directories : pybind11_incdir) \ No newline at end of file diff --git a/subprojects/pybind11.wrap b/subprojects/pybind11.wrap new file mode 100644 index 0000000..6906926 --- /dev/null +++ b/subprojects/pybind11.wrap @@ -0,0 +1,8 @@ +[wrap-git] +url = https://github.com/pybind/pybind11.git +revision = v3.0.0 +depth = 1 +patch_directory = pybind11 + +[provide] +pybind11 = pybind11_dep \ No newline at end of file