From 5a82311251e716b28467c6aacd19913ee54350ae Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 7 Apr 2026 12:19:58 -0400 Subject: [PATCH 1/4] feat(topology): Added TMOP support Meshes generated purley algebraically tend to be poorly conditioned. Incorporated MFEM's TMOP support based on a metric of ideal shape and unit size --- configs/default_config.toml | 3 + configs/test_external_domain.toml | 2 + .../test_external_domain_refinement_l1.toml | 2 + configs/test_flattening.toml | 2 + configs/test_polynomial_projection.toml | 17 ++ configs/test_refinement_l2.toml | 2 + configs/test_volume_no_external.toml | 2 + .../test_volume_spherical_no_external.toml | 2 + .../test_volume_spherical_with_external.toml | 2 + configs/test_volume_with_external.toml | 2 + src/include/stroid/config/config.h | 38 ++- src/include/stroid/stroid.h.in | 1 + src/include/stroid/topology/optimize.h | 14 + src/lib/topology/curvilinear.cpp | 13 +- src/lib/topology/mapping.cpp | 25 +- src/lib/topology/optimize.cpp | 231 ++++++++++++++++ src/lib/topology/topology.cpp | 16 +- src/meson.build | 1 + subprojects/libconfig.wrap | 2 +- tests/sandbox/sandbox_test.cpp | 10 +- tests/stroidTest.cpp | 253 ++++++++++++++---- tools/stroid.cpp | 6 + 22 files changed, 545 insertions(+), 101 deletions(-) create mode 100644 configs/test_polynomial_projection.toml create mode 100644 src/include/stroid/topology/optimize.h create mode 100644 src/lib/topology/optimize.cpp diff --git a/configs/default_config.toml b/configs/default_config.toml index f43489f..66137ea 100644 --- a/configs/default_config.toml +++ b/configs/default_config.toml @@ -13,3 +13,6 @@ surface_bdr_id = 1 core_id = 1 envelope_id = 2 vacuum_id = 3 + +[main.optimization_methods] +smoothstep = true \ No newline at end of file diff --git a/configs/test_external_domain.toml b/configs/test_external_domain.toml index fc5c0d4..84ab85b 100644 --- a/configs/test_external_domain.toml +++ b/configs/test_external_domain.toml @@ -13,3 +13,5 @@ surface_bdr_id = 1 core_id = 1 envelope_id = 2 vacuum_id = 3 +[main.optimization_methods] +smoothstep = true diff --git a/configs/test_external_domain_refinement_l1.toml b/configs/test_external_domain_refinement_l1.toml index bfe469c..c055b01 100644 --- a/configs/test_external_domain_refinement_l1.toml +++ b/configs/test_external_domain_refinement_l1.toml @@ -14,3 +14,5 @@ core_id = 1 envelope_id = 2 vacuum_id = 3 +[main.optimization_methods] +smoothstep = true diff --git a/configs/test_flattening.toml b/configs/test_flattening.toml index 0eaf396..e14b11c 100644 --- a/configs/test_flattening.toml +++ b/configs/test_flattening.toml @@ -13,3 +13,5 @@ surface_bdr_id = 1 core_id = 1 envelope_id = 2 vacuum_id = 3 +[main.optimization_methods] +smoothstep = true diff --git a/configs/test_polynomial_projection.toml b/configs/test_polynomial_projection.toml new file mode 100644 index 0000000..52f1057 --- /dev/null +++ b/configs/test_polynomial_projection.toml @@ -0,0 +1,17 @@ +[main] +core_steepness = 1.0 +flattening = 0.2 +include_external_domain = false +inf_bdr_id = 2 +order = 3 +r_core = 1.5 +r_infinity = 6.0 +r_instability = 1e-14 +r_star = 5.0 +refinement_levels = 1 +surface_bdr_id = 1 +core_id = 1 +envelope_id = 2 +vacuum_id = 3 +[main.optimization_methods] +smoothstep = true diff --git a/configs/test_refinement_l2.toml b/configs/test_refinement_l2.toml index fd61a99..6e05d6a 100644 --- a/configs/test_refinement_l2.toml +++ b/configs/test_refinement_l2.toml @@ -14,3 +14,5 @@ core_id = 1 envelope_id = 2 vacuum_id = 3 +[main.optimization_methods] +smoothstep = true diff --git a/configs/test_volume_no_external.toml b/configs/test_volume_no_external.toml index cf86efc..9cbdb53 100644 --- a/configs/test_volume_no_external.toml +++ b/configs/test_volume_no_external.toml @@ -14,3 +14,5 @@ core_id = 1 envelope_id = 2 vacuum_id = 3 +[main.optimization_methods] +smoothstep = true diff --git a/configs/test_volume_spherical_no_external.toml b/configs/test_volume_spherical_no_external.toml index 456ed85..066333e 100644 --- a/configs/test_volume_spherical_no_external.toml +++ b/configs/test_volume_spherical_no_external.toml @@ -14,3 +14,5 @@ core_id = 1 envelope_id = 2 vacuum_id = 3 +[main.optimization_methods] +smoothstep = true diff --git a/configs/test_volume_spherical_with_external.toml b/configs/test_volume_spherical_with_external.toml index 10bc77e..a227158 100644 --- a/configs/test_volume_spherical_with_external.toml +++ b/configs/test_volume_spherical_with_external.toml @@ -14,3 +14,5 @@ core_id = 1 envelope_id = 2 vacuum_id = 3 +[main.optimization_methods] +smoothstep = true diff --git a/configs/test_volume_with_external.toml b/configs/test_volume_with_external.toml index 617a2c7..9db35db 100644 --- a/configs/test_volume_with_external.toml +++ b/configs/test_volume_with_external.toml @@ -14,3 +14,5 @@ core_id = 1 envelope_id = 2 vacuum_id = 3 +[main.optimization_methods] +smoothstep = true diff --git a/src/include/stroid/config/config.h b/src/include/stroid/config/config.h index 416956a..f9227f0 100644 --- a/src/include/stroid/config/config.h +++ b/src/include/stroid/config/config.h @@ -1,6 +1,14 @@ #pragma once +#include + namespace stroid::config { + + struct OptimizationMethods { + std::optional tmop{false}; + std::optional smoothstep{true}; + }; + /** * @brief Configuration parameters for stroid mesh generation. * @@ -15,92 +23,94 @@ namespace stroid::config { * @section toml * - [main].refinement_levels */ - int refinement_levels = 4; + std::optional refinement_levels = 4; /** * @brief Polynomial order for high-order elements. * @section toml * - [main].order */ - int order = 3; + std::optional order = 3; /** * @brief Whether to include an external domain extending to `r_infinity`. * @section toml * - [main].include_external_domain */ - bool include_external_domain = true; + std::optional include_external_domain = true; /** * @brief Radius of the stellar core region. * @section toml * - [main].r_core */ - double r_core = 1.5; + std::optional r_core = 0.25; /** * @brief Radius of the stellar surface. * @section toml * - [main].r_star */ - double r_star = 5.0; + std::optional r_star = 1.0; /** * @brief Flattening factor for spheroidal shaping (0 = spherical, >0 = oblate). * @section toml * - [main].flattening */ - double flattening = 0; + std::optional flattening = 0; /** * @brief Outer radius of the external domain when enabled. * @section toml * - [main].r_infinity */ - double r_infinity = 6.0; + std::optional r_infinity = 6.0; /** * @brief Radius inside which transformations are skipped to avoid singularities. * @section toml * - [main].r_instability */ - double r_instability = 1e-14; + std::optional r_instability = 1e-14; /** * @brief Controls the smoothness/steepness of the core-to-envelope transition. * @section toml * - [main].core_steepness */ - double core_steepness = 1.0; + std::optional core_steepness = 1.0; /** * @brief Boundary attribute id for stellar surface * @section toml * - [main].surface_bdr_id */ - size_t surface_bdr_id = 1; + std::optional surface_bdr_id = 1; /** * @brief Boundary attribute id for infinity in kelvin mapping * @section toml * - [main].inf_bdr_id */ - size_t inf_bdr_id = 2; + std::optional inf_bdr_id = 2; /** * @brief Material attribute id for the core region * @section toml * - [main].core_id */ - size_t core_id = 1; + std::optional core_id = 1; /** * @brief Material attribute id for the envelope region * @section toml * - [main].envelope_id */ - size_t envelope_id = 2; + std::optional envelope_id = 2; /** * @brief Material attribute id for the external domain (if enabled) * @section toml * - [main].vacuum_id */ - size_t vacuum_id = 3; + std::optional vacuum_id = 3; + + std::optional optimization_methods = OptimizationMethods{true, true}; }; } diff --git a/src/include/stroid/stroid.h.in b/src/include/stroid/stroid.h.in index b596747..f221d63 100644 --- a/src/include/stroid/stroid.h.in +++ b/src/include/stroid/stroid.h.in @@ -4,6 +4,7 @@ #include "stroid/topology/topology.h" #include "stroid/topology/mapping.h" #include "stroid/topology/curvilinear.h" +#include "stroid/topology/optimize.h" #include "stroid/utils/mesh_utils.h" #include "stroid/IO/mesh.h" diff --git a/src/include/stroid/topology/optimize.h b/src/include/stroid/topology/optimize.h new file mode 100644 index 0000000..af44e95 --- /dev/null +++ b/src/include/stroid/topology/optimize.h @@ -0,0 +1,14 @@ +#pragma once + +#include "mfem.hpp" +#include "fourdst/config/base.h" +#include "stroid/config/config.h" + + +namespace stroid::topology { + + /** + * @breif Apply target matrix optimization to improve conditioning of the mesh + */ + void ApplyTMOP(mfem::Mesh& mesh, const fourdst::config::Config &config); +} diff --git a/src/lib/topology/curvilinear.cpp b/src/lib/topology/curvilinear.cpp index c67c148..1c8a47d 100644 --- a/src/lib/topology/curvilinear.cpp +++ b/src/lib/topology/curvilinear.cpp @@ -5,7 +5,7 @@ namespace stroid::topology { void PromoteToHighOrder(mfem::Mesh &mesh, const fourdst::config::Config &config) { - const auto* fec = new mfem::H1_FECollection(config->order, mesh.Dimension()); + const auto* fec = new mfem::H1_FECollection(config->order.value(), mesh.Dimension()); auto* fes = new mfem::FiniteElementSpace(&mesh, fec, mesh.SpaceDimension()); mesh.SetNodalFESpace(fes); } @@ -53,16 +53,5 @@ namespace stroid::topology { } } - // for (int i = 0; i < nDofs; ++i) { - // for (int d = 0; d < vDim; ++d) { - // pos(d) = nodes(fes->DofToVDof(i, d)); - // } - // - // TransformPoint(pos, config, 0); - // - // for (int d = 0; d < vDim; ++d) { - // nodes(fes->DofToVDof(i, d)) = pos(d); - // } - // } } } diff --git a/src/lib/topology/mapping.cpp b/src/lib/topology/mapping.cpp index bbc9105..c0b42fd 100644 --- a/src/lib/topology/mapping.cpp +++ b/src/lib/topology/mapping.cpp @@ -29,7 +29,7 @@ namespace stroid::topology { } void ApplySpheroidal(mfem::Vector &pos, const fourdst::config::Config &config) { - pos(2) *= (1.0 - config->flattening); + pos(2) *= (1.0 - config->flattening.value()); } void TransformPoint(mfem::Vector &pos, const fourdst::config::Config &config, int attribute_id) { @@ -49,8 +49,8 @@ namespace stroid::topology { unit_dir /= unit_dir.Norml2(); // Re-normalize if (l_inf <= config->r_core) { - const double t = l_inf / config->r_core; - double alpha = std::pow(t, config->core_steepness); + const double t = l_inf / config->r_core.value(); + double alpha = std::pow(t, config->core_steepness.value()); // Smoothstep function to apply C1 continuity alpha = alpha * alpha * (3.0 - 2.0 * alpha); @@ -59,18 +59,26 @@ namespace stroid::topology { mfem::Vector pos_spherical = unit_dir; pos_spherical *= l_inf; + bool run_smoothstep = false; - for (int d = 0; d < pos.Size(); ++d) { - pos(d) = (1.0 - alpha) * pos_cartesian(d) + alpha * pos_spherical(d); + 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) / (config->r_star - config->r_core); - const double r_phys = config->r_core + xi * (config->r_star - config->r_core); + 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; @@ -82,5 +90,4 @@ namespace stroid::topology { ApplySpheroidal(pos, config); } - } -} + }} diff --git a/src/lib/topology/optimize.cpp b/src/lib/topology/optimize.cpp new file mode 100644 index 0000000..3e9db8a --- /dev/null +++ b/src/lib/topology/optimize.cpp @@ -0,0 +1,231 @@ +#include "mfem.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include "stroid/topology/optimize.h" + + +#include +#include +#include +#include + +namespace stroid::utils::term_support { + inline bool locale_name_looks_utf8(const char* localeName) { + if (localeName == nullptr) return false; + + const std::string localeString(localeName); + + return localeString.find("UTF-8") != std::string::npos || + localeString.find("utf-8") != std::string::npos || + localeString.find("utf8") != std::string::npos || + localeString.find("UTF8") != std::string::npos; + } + + inline bool unicode_output_is_usable() { + const char* ctypeLocale = std::setlocale(LC_CTYPE, ""); + if (locale_name_looks_utf8(ctypeLocale)) { + return true; + } + + const char* lcAllEnv = std::getenv("LC_ALL"); + if (locale_name_looks_utf8(lcAllEnv)) { + return true; + } + + const char* lcCtypeEnv = std::getenv("LC_CTYPE"); + if (locale_name_looks_utf8(lcCtypeEnv)) { + return true; + } + + const char* langEnv = std::getenv("LANG"); + if (locale_name_looks_utf8(langEnv)) { + return true; + } + + return false; + } +} + +namespace stroid::topology { +class TMOPProgressBar : public mfem::IterativeSolverMonitor { + private: + double r0_ = -1.0; + double rtol_; + int bar_width_; + + std::atomic done_{false}; + std::atomic progress_{0.0}; + std::atomic iter_{0}; + std::atomic res_{0.0}; + + std::thread spinner_thread_; + + void Spin() { + std::vector spin_chars; + if (!utils::term_support::unicode_output_is_usable()) { + spin_chars= {"|", "/", "-", "\\"}; + } else { + spin_chars = {"▉", "▊", "▋", "▌", "▍", "▎", "▏", "▎", "▍", "▌", "▋", "▊", "▉"}; + } + int spin_idx = 0; + + while (!done_.load()) { + Draw(spin_chars[spin_idx]); + spin_idx = (spin_idx + 1) % spin_chars.size(); + std::this_thread::sleep_for(std::chrono::milliseconds(100)); + } + } + + void Draw(const std::string& spinner) { + const double p = progress_.load(); + const int pos = static_cast(bar_width_ * p); + + std::cout << "\r[" << spinner << "] TMOP Relaxation ["; + for (int i = 0; i < bar_width_; ++i) { + if (i < pos) std::cout << "="; + else if (i == pos) std::cout << ">"; + else std::cout << " "; + } + std::cout << "] " << std::setw(3) << static_cast(p * 100.0) << "% " + << "(Iter: " << std::setw(2) << iter_.load() + << ", Res: " << std::scientific << std::setprecision(2) << res_.load() << ") " << std::flush; + } + + public: + TMOPProgressBar(double rel_tol, int width = 50) + : rtol_(rel_tol), bar_width_(width) { + spinner_thread_ = std::thread(&TMOPProgressBar::Spin, this); + } + + ~TMOPProgressBar() override { + if (!done_.load()) { + done_ = true; + if (spinner_thread_.joinable()) { + spinner_thread_.join(); + } + } + } + + void MonitorResidual(int it, double norm, const mfem::Vector &r, bool final) override { + if (it == 0 || r0_ < 0.0) { + r0_ = norm; + } + + iter_ = it; + res_ = norm; + + double p = 0.0; + const double target_norm = r0_ * rtol_; + + if (norm <= target_norm || final) { + p = 1.0; + } else if (norm < r0_ && r0_ > 0.0 && target_norm > 0.0) { + const double log_start = std::log10(r0_); + const double log_current = std::log10(norm); + const double log_target = std::log10(target_norm); + p = (log_start - log_current) / (log_start - log_target); + p = std::clamp(p, 0.0, 1.0); + } + + progress_ = p; + + if (final) { + done_ = true; + if (spinner_thread_.joinable()) { + spinner_thread_.join(); + } + + Draw("*"); + std::cout << std::endl; + } + } + }; + void ApplyTMOP(mfem::Mesh &mesh, const fourdst::config::Config &config) { + const mfem::FiniteElementSpace* cfes = mesh.GetNodalFESpace(); + mfem::FiniteElementSpace* fes = const_cast(cfes); + + if (!fes) { + std::cerr << "Error: Mesh has no nodal finite element space. Call PromoteToHighOrder first." << std::endl; + return; + } + + const int max_bdr_attr = mesh.bdr_attributes.Size() > 0 ? mesh.bdr_attributes.Max() : 0; + mfem::Array ess_bdr(max_bdr_attr); + ess_bdr = 0.0; + + if (max_bdr_attr >= config->surface_bdr_id.value()) { + ess_bdr[config->surface_bdr_id.value() - 1] = 1; + } + + if (config->include_external_domain.value() && max_bdr_attr >= config->inf_bdr_id.value()) { + ess_bdr[config->inf_bdr_id.value() - 1] = 1; + } + + mfem::Array ess_tdof_list; + 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::TMOP_Integrator* tmop_integrator = new mfem::TMOP_Integrator(metric, target_c); + + mfem::NonlinearForm a(fes); + a.AddDomainIntegrator(tmop_integrator); + a.SetEssentialTrueDofs(ess_tdof_list); + + mfem::GridFunction* nodes = mesh.GetNodes(); + mfem::Vector x(*nodes); + mfem::Vector b(a.Height()); + b = 0.0; + + mfem::MINRESSolver minres; + minres.SetMaxIter(500); + minres.SetRelTol(1e-5); + minres.SetAbsTol(0.0); + minres.SetPrintLevel(0); + + mfem::DSmoother jacobi(1, 1.0, 1); + jacobi.SetPositiveDiagonal(true); + minres.SetPreconditioner(jacobi); + + const int quad_order = 2 * fes->GetMaxElementOrder() + 3; + const mfem::IntegrationRule &ir = mfem::IntRules.Get(mesh.GetTypicalElementGeometry(), quad_order); + + double min_detJ = std::numeric_limits::infinity(); + for (int i = 0; i < mesh.GetNE(); i++) { + mfem::ElementTransformation *T = mesh.GetElementTransformation(i); + for (int j = 0; j < ir.GetNPoints(); j++) { + T->SetIntPoint(&ir.IntPoint(j)); + min_detJ = std::min(min_detJ, T->Jacobian().Det()); + } + } + + constexpr double newton_rtol = 1e-4; + mfem::TMOPNewtonSolver newton(ir, 0); + newton.SetPreconditioner(minres); + newton.SetOperator(a); + newton.SetMaxIter(50); + newton.SetRelTol(newton_rtol); + newton.SetAbsTol(0.0); + newton.SetMinDetPtr(&min_detJ); + newton.SetPrintLevel(0); + + TMOPProgressBar progress_bar(newton_rtol); + newton.SetMonitor(progress_bar); + + std::cout << "Applying TMOP optimization to mesh. Note this may take a long time. Depending on your mesh resolution expect to wait up to the order of 10s of minutes..." << std::endl; + newton.Mult(b, x); + *nodes = x; + + mesh.NodesUpdated(); + + delete metric; + delete target_c; + } +} diff --git a/src/lib/topology/topology.cpp b/src/lib/topology/topology.cpp index ddc4286..0672899 100644 --- a/src/lib/topology/topology.cpp +++ b/src/lib/topology/topology.cpp @@ -21,14 +21,14 @@ namespace stroid::topology { mesh->AddVertex(x, y, z); }; - add_box(config->r_core); - add_box(config->r_star); + add_box(config->r_core.value()); + add_box(config->r_star.value()); if (config->include_external_domain) { - add_box(config->r_infinity); + add_box(config->r_infinity.value()); } const int core_v[8] = {0, 1, 3, 2, 4, 5, 7, 6}; - mesh->AddHex(core_v, config->core_id); + mesh->AddHex(core_v, config->core_id.value()); std::vector> stellar_shells = { {8, 9, 11, 10, 0, 1, 3, 2}, @@ -39,7 +39,7 @@ namespace stroid::topology { {0, 4, 6, 2, 8, 12, 14, 10} // -X face }; for (const auto & shell : stellar_shells) { - mesh->AddHex(shell.data(), config->envelope_id); + mesh->AddHex(shell.data(), config->envelope_id.value()); } if (config->include_external_domain) { @@ -51,7 +51,7 @@ namespace stroid::topology { vacuum_shells.push_back({12, 13, 15, 14, 20, 21, 23, 22}); vacuum_shells.push_back({10, 11, 9, 8, 18, 19, 17, 16}); for (const auto & shell : vacuum_shells) { - mesh->AddHex(shell.data(), config->vacuum_id); + mesh->AddHex(shell.data(), config->vacuum_id.value()); } } @@ -66,7 +66,7 @@ namespace stroid::topology { }; for (const auto& bdr: surface_bdr_quads) { - mesh->AddBdrQuad(bdr, config->surface_bdr_id); + mesh->AddBdrQuad(bdr, config->surface_bdr_id.value()); } if (config->include_external_domain) { @@ -80,7 +80,7 @@ namespace stroid::topology { }; for (const auto& bdr: inf_bdr_quads) { - mesh->AddBdrQuad(bdr, config->inf_bdr_id); + mesh->AddBdrQuad(bdr, config->inf_bdr_id.value()); } } diff --git a/src/meson.build b/src/meson.build index 66e2a81..e50ec57 100644 --- a/src/meson.build +++ b/src/meson.build @@ -10,6 +10,7 @@ stroid_sources = files( 'lib/topology/curvilinear.cpp', 'lib/topology/mapping.cpp', 'lib/topology/topology.cpp', + 'lib/topology/optimize.cpp', 'lib/IO/mesh.cpp', 'lib/utils/mesh_utils.cpp', ) diff --git a/subprojects/libconfig.wrap b/subprojects/libconfig.wrap index 86267ec..154b33d 100644 --- a/subprojects/libconfig.wrap +++ b/subprojects/libconfig.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/libconfig.git -revision = v2.0.5 +revision = v2.2.1 depth = 1 diff --git a/tests/sandbox/sandbox_test.cpp b/tests/sandbox/sandbox_test.cpp index bb2359d..dd04232 100644 --- a/tests/sandbox/sandbox_test.cpp +++ b/tests/sandbox/sandbox_test.cpp @@ -9,6 +9,8 @@ #include +#include "stroid/topology/optimize.h" + struct SandboxConfig { std::string host = "localhost"; int port = 19916; @@ -22,14 +24,20 @@ int main() { MeshConfig mesh_cfg; mesh_cfg.load("default.toml"); - UserConfig user_cfg; + const UserConfig user_cfg; std::unique_ptr mesh = stroid::topology::BuildSkeleton(mesh_cfg); stroid::topology::Finalize(*mesh, mesh_cfg); stroid::topology::PromoteToHighOrder(*mesh, mesh_cfg); stroid::topology::ProjectMesh(*mesh, mesh_cfg); + + if (mesh_cfg->optimization_methods.has_value() && mesh_cfg->optimization_methods.value().tmop.has_value() && mesh_cfg->optimization_methods.value().tmop.value()) { + stroid::topology::ApplyTMOP(*mesh, mesh_cfg); + } + stroid::IO::ViewMesh(*mesh, "Sandbox Mesh", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, user_cfg->host, user_cfg->port); + stroid::IO::SaveMesh(*mesh, "sandbox.mesh"); return 0; diff --git a/tests/stroidTest.cpp b/tests/stroidTest.cpp index 845beaf..b18fc37 100644 --- a/tests/stroidTest.cpp +++ b/tests/stroidTest.cpp @@ -17,6 +17,8 @@ #include #include #include +#include +#include namespace { @@ -32,10 +34,10 @@ std::filesystem::path GetSourceRoot() { return std::filesystem::current_path(); } -Config LoadConfigFromRepo(const std::filesystem::path& relative_path) { - Config cfg; - cfg.load((GetSourceRoot() / relative_path).string()); - return cfg; +std::unique_ptr LoadConfigFromRepo(const std::filesystem::path& relative_path) { + auto cfg_ptr = std::make_unique(); + cfg_ptr->load((GetSourceRoot() / relative_path).string()); + return cfg_ptr; } @@ -122,13 +124,13 @@ std::unique_ptr BuildProjectedMesh(const Config& cfg) { double ComputeStellarVolumeWithDomainLFIntegrator(mfem::Mesh& mesh, const Config& cfg) { const int mesh_max_attr = mesh.attributes.Size() > 0 ? mesh.attributes.Max() : 0; - const int cfg_max_attr = static_cast(std::max({cfg->core_id, cfg->envelope_id, cfg->vacuum_id})); + const int cfg_max_attr = static_cast(std::max({cfg->core_id.value(), cfg->envelope_id.value(), cfg->vacuum_id.value()})); const int coeff_size = std::max(1, std::max(mesh_max_attr, cfg_max_attr)); mfem::Vector attr_coeff(coeff_size); attr_coeff = 0.0; - attr_coeff(static_cast(cfg->core_id) - 1) = 1.0; - attr_coeff(static_cast(cfg->envelope_id) - 1) = 1.0; + attr_coeff(static_cast(cfg->core_id.value()) - 1) = 1.0; + attr_coeff(static_cast(cfg->envelope_id.value()) - 1) = 1.0; mfem::PWConstCoefficient stellar_coeff(attr_coeff); mfem::L2_FECollection fec(0, mesh.Dimension()); @@ -233,6 +235,26 @@ ConditioningStats CollectConditioningStats(const mfem::Mesh& mesh, const std::se return stats; } +std::optional EvalGridFunctionAtPoint( + mfem::Mesh& mesh, + const mfem::Vector& x, + const mfem::GridFunction& u ){ + + mfem::Array elem_ids; + mfem::Array ips; + mfem::DenseMatrix P(x.Size(), 1); + P.SetCol(0, x); + + mesh.FindPoints(P, elem_ids, ips, false); + + if (elem_ids.Size() > 0 && elem_ids[0] >= 0) { + return u.GetValue(elem_ids[0], ips[0]); + } else { + return std::nullopt; + + } +} + } // namespace /** @@ -270,7 +292,8 @@ TEST_F(stroidTest, BuildSkeleton_DefaultCounts) { * `src/lib/topology/topology.cpp` (`vacuum_shells`, `inf_bdr_quads`) and config parsing path. */ TEST_F(stroidTest, BuildSkeleton_ExternalDomainCounts) { - const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto& cfg = *cfg_ptr; const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); ASSERT_NE(mesh, nullptr); @@ -289,19 +312,20 @@ TEST_F(stroidTest, BuildSkeleton_ExternalDomainCounts) { * `core_id`, `envelope_id`, `vacuum_id`, `surface_bdr_id`, `inf_bdr_id` in config fixtures. */ TEST_F(stroidTest, BuildSkeleton_ExternalDomainAttributes) { - const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto& cfg = *cfg_ptr; const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); ASSERT_NE(mesh, nullptr); const auto volume_attr_counts = CountVolumeAttributes(*mesh); - EXPECT_EQ(volume_attr_counts.at(static_cast(cfg->core_id)), 1); - EXPECT_EQ(volume_attr_counts.at(static_cast(cfg->envelope_id)), 6); - EXPECT_EQ(volume_attr_counts.at(static_cast(cfg->vacuum_id)), 6); + EXPECT_EQ(volume_attr_counts.at(static_cast(cfg->core_id.value())), 1); + EXPECT_EQ(volume_attr_counts.at(static_cast(cfg->envelope_id.value())), 6); + EXPECT_EQ(volume_attr_counts.at(static_cast(cfg->vacuum_id.value())), 6); const auto boundary_attr_counts = CountBoundaryAttributes(*mesh); - EXPECT_EQ(boundary_attr_counts.at(static_cast(cfg->surface_bdr_id)), 6); - EXPECT_EQ(boundary_attr_counts.at(static_cast(cfg->inf_bdr_id)), 6); + EXPECT_EQ(boundary_attr_counts.at(static_cast(cfg->surface_bdr_id.value())), 6); + EXPECT_EQ(boundary_attr_counts.at(static_cast(cfg->inf_bdr_id.value())), 6); } @@ -333,7 +357,8 @@ TEST_F(stroidTest, Finalize_RefinementIncreasesElements) { * If this fails: inspect refine-loop count and any topology-side early exits in `Finalize`. */ TEST_F(stroidTest, Finalize_DefaultRefinementScalesHexCountByEightPowerL) { - const Config cfg = LoadConfigFromRepo("configs/test_refinement_l2.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_refinement_l2.toml"); + const auto& cfg = *cfg_ptr; const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); const int initial_elements = mesh->GetNE(); @@ -352,7 +377,8 @@ TEST_F(stroidTest, Finalize_DefaultRefinementScalesHexCountByEightPowerL) { * If this fails: inspect `Finalize` and verify external-domain elements are not excluded from refinement. */ TEST_F(stroidTest, Finalize_ExternalDomainRefinementScalesHexCountByEightPowerL) { - const Config cfg = LoadConfigFromRepo("configs/test_external_domain_refinement_l1.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain_refinement_l1.toml"); + const auto& cfg = *cfg_ptr; const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); const int initial_elements = mesh->GetNE(); @@ -371,7 +397,9 @@ TEST_F(stroidTest, Finalize_ExternalDomainRefinementScalesHexCountByEightPowerL) * If this fails: inspect `Finalize` orientation/refinement calls and any attribute mutation side effects. */ TEST_F(stroidTest, Finalize_ExternalDomainConformingAndRefined) { - const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto& cfg = *cfg_ptr; + const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); const int initial_elements = mesh->GetNE(); @@ -381,13 +409,13 @@ TEST_F(stroidTest, Finalize_ExternalDomainConformingAndRefined) { EXPECT_GT(mesh->GetNE(), initial_elements); const auto volume_attr_counts = CountVolumeAttributes(*mesh); - EXPECT_GT(volume_attr_counts.at(static_cast(cfg->core_id)), 0); - EXPECT_GT(volume_attr_counts.at(static_cast(cfg->envelope_id)), 0); - EXPECT_GT(volume_attr_counts.at(static_cast(cfg->vacuum_id)), 0); + EXPECT_GT(volume_attr_counts.at(static_cast(cfg->core_id.value())), 0); + EXPECT_GT(volume_attr_counts.at(static_cast(cfg->envelope_id.value())), 0); + EXPECT_GT(volume_attr_counts.at(static_cast(cfg->vacuum_id.value())), 0); const auto boundary_attr_counts = CountBoundaryAttributes(*mesh); - EXPECT_GT(boundary_attr_counts.at(static_cast(cfg->surface_bdr_id)), 0); - EXPECT_GT(boundary_attr_counts.at(static_cast(cfg->inf_bdr_id)), 0); + EXPECT_GT(boundary_attr_counts.at(static_cast(cfg->surface_bdr_id.value())), 0); + EXPECT_GT(boundary_attr_counts.at(static_cast(cfg->inf_bdr_id.value())), 0); } /** @@ -399,16 +427,18 @@ TEST_F(stroidTest, Finalize_ExternalDomainConformingAndRefined) { * notably `src/lib/topology/topology.cpp` and `src/lib/utils/mesh_utils.cpp`. */ TEST_F(stroidTest, Finalize_ExternalDomainKeepsOnlyExpectedMaterialAndBoundaryIDs) { - const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto& cfg = *cfg_ptr; + const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); stroid::topology::Finalize(*mesh, cfg); const auto volume_attr_counts = CountVolumeAttributes(*mesh); const std::set expected_volume_ids = { - static_cast(cfg->core_id), - static_cast(cfg->envelope_id), - static_cast(cfg->vacuum_id) + static_cast(cfg->core_id.value()), + static_cast(cfg->envelope_id.value()), + static_cast(cfg->vacuum_id.value()) }; for (const auto& [attr, count] : volume_attr_counts) { EXPECT_TRUE(expected_volume_ids.contains(attr)); @@ -418,8 +448,8 @@ TEST_F(stroidTest, Finalize_ExternalDomainKeepsOnlyExpectedMaterialAndBoundaryID const auto boundary_attr_counts = CountBoundaryAttributes(*mesh); const std::set expected_boundary_ids = { - static_cast(cfg->surface_bdr_id), - static_cast(cfg->inf_bdr_id) + static_cast(cfg->surface_bdr_id.value()), + static_cast(cfg->inf_bdr_id.value()) }; for (const auto& [attr, count] : boundary_attr_counts) { EXPECT_TRUE(expected_boundary_ids.contains(attr)); @@ -496,7 +526,8 @@ TEST_F(stroidTest, ApplyEquiangular_BasicTransform) { * `configs/test_flattening.toml`. */ TEST_F(stroidTest, ApplySpheroidal_FlattensZ) { - const Config cfg = LoadConfigFromRepo("configs/test_flattening.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_flattening.toml"); + const auto& cfg = *cfg_ptr; mfem::Vector pos(3); pos(0) = 0.0; @@ -571,9 +602,9 @@ TEST_F(stroidTest, TransformPoint_IsContinuousAcrossCoreAndStarInterfaces) { dir(2) = -0.4; mfem::Vector near_core_left = dir; - near_core_left *= cfg->r_core * (1.0 - eps); + near_core_left *= cfg->r_core.value() * (1.0 - eps); mfem::Vector near_core_right = dir; - near_core_right *= cfg->r_core * (1.0 + eps); + near_core_right *= cfg->r_core.value() * (1.0 + eps); const mfem::Vector core_left_mapped = TransformCopy(near_core_left, cfg); const mfem::Vector core_right_mapped = TransformCopy(near_core_right, cfg); @@ -583,9 +614,9 @@ TEST_F(stroidTest, TransformPoint_IsContinuousAcrossCoreAndStarInterfaces) { EXPECT_LT(diff.Norml2(), 1e-3); mfem::Vector near_star_left = dir; - near_star_left *= cfg->r_star * (1.0 - eps); + near_star_left *= cfg->r_star.value() * (1.0 - eps); mfem::Vector near_star_right = dir; - near_star_right *= cfg->r_star * (1.0 + eps); + near_star_right *= cfg->r_star.value() * (1.0 + eps); const mfem::Vector star_left_mapped = TransformCopy(near_star_left, cfg); const mfem::Vector star_right_mapped = TransformCopy(near_star_right, cfg); @@ -684,7 +715,9 @@ TEST_F(stroidTest, EndToEnd_BuildFinalizePromoteProject) { * If this fails: inspect external-domain topology assembly and projection loops over mixed attributes. */ TEST_F(stroidTest, EndToEnd_ExternalDomainBuildFinalizePromoteProject) { - const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto cfg_ptr= LoadConfigFromRepo("configs/test_external_domain.toml"); + const auto& cfg = *cfg_ptr; + const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); stroid::topology::Finalize(*mesh, cfg); stroid::topology::PromoteToHighOrder(*mesh, cfg); @@ -705,8 +738,12 @@ TEST_F(stroidTest, EndToEnd_ExternalDomainBuildFinalizePromoteProject) { * that may leak starside nodes into vacuum geometry. */ TEST_F(stroidTest, Volume_StellarDomainMatchesWithAndWithoutExternalDomain) { - const Config no_external_cfg = LoadConfigFromRepo("configs/test_volume_no_external.toml"); - const Config with_external_cfg = LoadConfigFromRepo("configs/test_volume_with_external.toml"); + const auto no_external_cfg_ptr = LoadConfigFromRepo("configs/test_volume_no_external.toml"); + const auto with_external_cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml"); + + const auto& no_external_cfg = *no_external_cfg_ptr; + const auto& with_external_cfg = *with_external_cfg_ptr; + const std::unique_ptr no_external_mesh = stroid::topology::BuildSkeleton(no_external_cfg); stroid::topology::Finalize(*no_external_mesh, no_external_cfg); @@ -719,12 +756,12 @@ TEST_F(stroidTest, Volume_StellarDomainMatchesWithAndWithoutExternalDomain) { stroid::topology::ProjectMesh(*with_external_mesh, with_external_cfg); const std::set stellar_attrs_no_external = { - static_cast(no_external_cfg->core_id), - static_cast(no_external_cfg->envelope_id) + static_cast(no_external_cfg->core_id.value()), + static_cast(no_external_cfg->envelope_id.value()) }; const std::set stellar_attrs_with_external = { - static_cast(with_external_cfg->core_id), - static_cast(with_external_cfg->envelope_id) + static_cast(with_external_cfg->core_id.value()), + static_cast(with_external_cfg->envelope_id.value()) }; const double stellar_volume_no_external = ComputeMeshVolumeForAttributes(*no_external_mesh, stellar_attrs_no_external); @@ -744,7 +781,9 @@ TEST_F(stroidTest, Volume_StellarDomainMatchesWithAndWithoutExternalDomain) { * If this fails: inspect `ComputeMeshVolume*` helpers and region attribute IDs in config fixtures. */ TEST_F(stroidTest, Volume_ExternalMeshExcludesVacuumWhenRequested) { - const Config cfg = LoadConfigFromRepo("configs/test_volume_with_external.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml"); + const auto& cfg = *cfg_ptr; + const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); stroid::topology::Finalize(*mesh, cfg); @@ -752,10 +791,10 @@ TEST_F(stroidTest, Volume_ExternalMeshExcludesVacuumWhenRequested) { stroid::topology::ProjectMesh(*mesh, cfg); const std::set stellar_attrs = { - static_cast(cfg->core_id), - static_cast(cfg->envelope_id) + static_cast(cfg->core_id.value()), + static_cast(cfg->envelope_id.value()) }; - const std::set vacuum_attr = {static_cast(cfg->vacuum_id)}; + const std::set vacuum_attr = {static_cast(cfg->vacuum_id.value())}; const double total_volume = ComputeMeshVolume(*mesh); const double stellar_volume = ComputeMeshVolumeForAttributes(*mesh, stellar_attrs); @@ -775,7 +814,8 @@ TEST_F(stroidTest, Volume_ExternalMeshExcludesVacuumWhenRequested) { * `IntegrateElementVolume`. */ TEST_F(stroidTest, Volume_SphericalStellarDomainMatchesAnalyticSphere) { - const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml"); + const auto& cfg = *cfg_ptr; const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); stroid::topology::Finalize(*mesh, cfg); @@ -783,12 +823,12 @@ TEST_F(stroidTest, Volume_SphericalStellarDomainMatchesAnalyticSphere) { stroid::topology::ProjectMesh(*mesh, cfg); const std::set stellar_attrs = { - static_cast(cfg->core_id), - static_cast(cfg->envelope_id) + static_cast(cfg->core_id.value()), + static_cast(cfg->envelope_id.value()) }; const double measured_volume = ComputeMeshVolumeForAttributes(*mesh, stellar_attrs); - const double analytic_volume = 4.0 / 3.0 * kPi * std::pow(cfg->r_star, 3.0); + const double analytic_volume = 4.0 / 3.0 * kPi * std::pow(cfg->r_star.value(), 3.0); const double rel_err = std::abs(measured_volume - analytic_volume) / analytic_volume; EXPECT_LT(rel_err, 1e-2); @@ -804,11 +844,12 @@ TEST_F(stroidTest, Volume_SphericalStellarDomainMatchesAnalyticSphere) { * and MFEM assembly setup in this test file. */ TEST_F(stroidTest, Volume_SphericalStellarDomainDomainLFIntegratorMatchesAnalyticSphere) { - const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_with_external.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_with_external.toml"); + const auto& cfg = *cfg_ptr; std::unique_ptr mesh = BuildProjectedMesh(cfg); const double measured_volume = ComputeStellarVolumeWithDomainLFIntegrator(*mesh, cfg); - const double analytic_volume = 4.0 / 3.0 * kPi * std::pow(cfg->r_star, 3.0); + const double analytic_volume = 4.0 / 3.0 * kPi * std::pow(cfg->r_star.value(), 3.0); const double rel_err = std::abs(measured_volume - analytic_volume) / analytic_volume; EXPECT_LT(rel_err, 1e-2); @@ -823,7 +864,9 @@ TEST_F(stroidTest, Volume_SphericalStellarDomainDomainLFIntegratorMatchesAnalyti * refinement/order config used by `configs/test_volume_spherical_no_external.toml`. */ TEST_F(stroidTest, Conditioning_DefaultMeshHasPositiveJacobiansAndReasonableShape) { - const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml"); + const auto& cfg = *cfg_ptr; + const std::unique_ptr mesh = BuildProjectedMesh(cfg); const ConditioningStats stats = CollectConditioningStats(*mesh, {}); @@ -845,12 +888,14 @@ TEST_F(stroidTest, Conditioning_DefaultMeshHasPositiveJacobiansAndReasonableShap * assignment in `BuildSkeleton`. */ TEST_F(stroidTest, Conditioning_ExternalMeshPerRegionHasPositiveJacobians) { - const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_with_external.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_with_external.toml"); + const auto& cfg = *cfg_ptr; + const std::unique_ptr mesh = BuildProjectedMesh(cfg); - const ConditioningStats core_stats = CollectConditioningStats(*mesh, {static_cast(cfg->core_id)}); - const ConditioningStats envelope_stats = CollectConditioningStats(*mesh, {static_cast(cfg->envelope_id)}); - const ConditioningStats vacuum_stats = CollectConditioningStats(*mesh, {static_cast(cfg->vacuum_id)}); + const ConditioningStats core_stats = CollectConditioningStats(*mesh, {static_cast(cfg->core_id.value())}); + const ConditioningStats envelope_stats = CollectConditioningStats(*mesh, {static_cast(cfg->envelope_id.value())}); + const ConditioningStats vacuum_stats = CollectConditioningStats(*mesh, {static_cast(cfg->vacuum_id.value())}); ASSERT_GT(core_stats.samples, 0); ASSERT_GT(envelope_stats.samples, 0); @@ -874,7 +919,9 @@ TEST_F(stroidTest, Conditioning_ExternalMeshPerRegionHasPositiveJacobians) { * `src/lib/utils/mesh_utils.cpp`, then trace upstream mapping changes. */ TEST_F(stroidTest, Conditioning_DefaultMeshHasNoFlippedElementsOrBoundaryFaces) { - const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml"); + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml"); + const auto& cfg = *cfg_ptr; + std::unique_ptr mesh = BuildProjectedMesh(cfg); stroid::utils::MarkFlippedElements(*mesh); @@ -887,3 +934,97 @@ TEST_F(stroidTest, Conditioning_DefaultMeshHasNoFlippedElementsOrBoundaryFaces) EXPECT_FALSE(boundary_attr_counts.contains(500)); } +TEST_F(stroidTest, PolynomainalProjection) { + const auto cfg_ptr = LoadConfigFromRepo("configs/test_polynomial_projection.toml"); + const auto& cfg = *cfg_ptr; + + std::unique_ptr mesh = BuildProjectedMesh(cfg); + + const int geom_order = mesh->GetNodes()->FESpace()->GetMaxElementOrder(); + const int space_dim = mesh->Dimension(); + + mfem::H1_FECollection fec(geom_order, space_dim); + mfem::FiniteElementSpace fes(mesh.get(), &fec); + + auto ProjectedFunction = [](const mfem::Vector& x) { + const double r = x.Norml2(); + return 1 + 7 * r * r - 2 * r; + }; + + mfem::GridFunction projected_u(&fes); + mfem::FunctionCoefficient u_coeff(ProjectedFunction); + projected_u.ProjectCoefficient(u_coeff); + + mfem::Vector x(space_dim); + x = 0.0; + for (double t = 0; t <= 1; t+= 0.01) { + x(0) = t; + + double analytic_val = ProjectedFunction(x); + double projected_val = EvalGridFunctionAtPoint(*mesh, x, projected_u).value_or(std::numeric_limits::quiet_NaN()); + + double rel_err = std::abs(projected_val - analytic_val) / analytic_val; + EXPECT_LT(rel_err, 1e-12); + } +} + +TEST_F(stroidTest, TranscendtalProjection) { + const auto cfg_ptr = LoadConfigFromRepo("configs/test_polynomial_projection.toml"); + const auto& cfg = *cfg_ptr; + + std::unique_ptr mesh = BuildProjectedMesh(cfg); + + const int geom_order = mesh->GetNodes()->FESpace()->GetMaxElementOrder(); + const int space_dim = mesh->Dimension(); + + mfem::H1_FECollection fec(geom_order, space_dim); + mfem::FiniteElementSpace fes(mesh.get(), &fec); + + + auto ProjectedFunction = [](const mfem::Vector& x) { + const double r = x.Norml2(); + if (r <= 1e-8) return 1.0; + return std::sin(r)/r; + }; + + auto expansion = [](const double t, const int order) { + double val = 0.0; + for (int k = 0; k < order; ++k) { + const double sign = (k % 2 == 0) ? 1.0 : -1.0; + const double term = sign * std::pow(t, 2 * k) / std::tgamma(2 * k + 2); + val += term; + } + return val; + }; + + auto expansion_err = [geom_order, &expansion](const double r) { + const double expansion_val = expansion(r, geom_order); + + const double analytic_val = std::sin(r)/r; + return std::abs((expansion_val - analytic_val))/std::abs(analytic_val); + }; + + double max_estimated_truncation_error = 0.0; + for (double t = 0; t < 1; t+= 0.01) { + double trunc_err = expansion_err(t); + max_estimated_truncation_error = std::max(max_estimated_truncation_error, trunc_err); + } + + mfem::GridFunction projected_u(&fes); + mfem::FunctionCoefficient u_coeff(ProjectedFunction); + projected_u.ProjectCoefficient(u_coeff); + + mfem::Vector x(space_dim); + x = 0.0; + for (double t = 0; t <= 1; t+= 0.01) { + x(0) = t; + + double analytic_val = ProjectedFunction(x); + double projected_val = EvalGridFunctionAtPoint(*mesh, x, projected_u).value_or(std::numeric_limits::quiet_NaN()); + + double rel_err = std::abs(projected_val - analytic_val) / analytic_val; + EXPECT_LT(rel_err, 10*max_estimated_truncation_error); + } +} + + diff --git a/tools/stroid.cpp b/tools/stroid.cpp index e1cd2b8..67ad25d 100644 --- a/tools/stroid.cpp +++ b/tools/stroid.cpp @@ -132,10 +132,16 @@ int main(int argc, char** argv) { cfg.load(config_filename.value()); } + const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); stroid::topology::Finalize(*mesh, cfg); stroid::topology::PromoteToHighOrder(*mesh, cfg); stroid::topology::ProjectMesh(*mesh, cfg); + if (cfg->optimization_methods.has_value() && cfg->optimization_methods.value().tmop.has_value() && cfg->optimization_methods.value().tmop.value()) { + stroid::topology::ApplyTMOP(*mesh, cfg); + } + + if (!no_save) { const std::string& final_path = output_filename; From 37416adb03da86260276acdd6810ddb65148f434 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 7 Apr 2026 12:58:16 -0400 Subject: [PATCH 2/4] feat(stroid): added mesh viewer and tmop toggle --- Doxyfile | 2 +- meson.build | 2 +- readme.md | 39 ++++++++++-------- src/include/stroid/topology/optimize.h | 5 +++ src/lib/topology/optimize.cpp | 6 +++ subprojects/libconfig.wrap | 2 +- tools/stroid.cpp | 57 +++++++++++++++++++++++--- 7 files changed, 88 insertions(+), 25 deletions(-) diff --git a/Doxyfile b/Doxyfile index 5aa5c83..757d19f 100644 --- a/Doxyfile +++ b/Doxyfile @@ -48,7 +48,7 @@ PROJECT_NAME = stroid # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = v0.2.1 +PROJECT_NUMBER = v0.3.0 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewers a diff --git a/meson.build b/meson.build index 48d8eed..bc13b2c 100644 --- a/meson.build +++ b/meson.build @@ -1,4 +1,4 @@ -project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.2.1', default_options : ['cpp_std=c++23']) +project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.3.0', default_options : ['cpp_std=c++23']) subdir('build-check') diff --git a/readme.md b/readme.md index 534cb22..ce1545a 100644 --- a/readme.md +++ b/readme.md @@ -92,25 +92,31 @@ inf_bdr_id = 2 core_id = 1 envelope_id = 2 vacuum_id = 3 + +[main.optimization_methods] +tmop = false +smoothstep = true ``` -| Parameter | Description | Default | -|-------------------------|-----------------------------------------------------------------------------------------------------|---------| -| refinement_levels | Number of uniform refinement levels to apply to the mesh after generation | 4 | -| order | The polynomial order of the finite elements in the mesh | 3 | -| include_external_domain | Whether to include an external domain extending to r_infinity | true | -| r_core | The radius of the core region of the star | 1.5 | -| r_star | The radius of the star | 5.0 | -| flattening | The flattening factor of the star (0 for spherical, >0 for oblate) | 0 | -| r_infinity | The outer radius of the external domain (if included) | 6.0 | -| r_instability | The radius at which no transformations are applied to the initial topology (to avoid singularities) | 1e-14 | -| core_steepness | The steepness of the transition between the core and envelope regions of the star | 1.0 | -| surface_bdr_id | The boundary ID to assign to the surface of the star | 1 | -| inf_bdr_id | The boundary ID to assign to the outer boundary of the external domain (if included) | 2 | -| core_id | The material ID to assign to the core region of the star | 1 | -| envelope_id | The material ID to assign to the envelope region of the star | 2 | -| vacuum_id | The material ID to assign to the vacuum region of the star (if included) | 3 | +| Parameter | Description | Default | +|---------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|---------| +| refinement_levels | Number of uniform refinement levels to apply to the mesh after generation | 4 | +| order | The polynomial order of the finite elements in the mesh | 3 | +| include_external_domain | Whether to include an external domain extending to r_infinity | true | +| r_core | The radius of the core region of the star | 1.5 | +| r_star | The radius of the star | 5.0 | +| flattening | The flattening factor of the star (0 for spherical, >0 for oblate) | 0 | +| r_infinity | The outer radius of the external domain (if included) | 6.0 | +| r_instability | The radius at which no transformations are applied to the initial topology (to avoid singularities) | 1e-14 | +| core_steepness | The steepness of the transition between the core and envelope regions of the star | 1.0 | +| surface_bdr_id | The boundary ID to assign to the surface of the star | 1 | +| inf_bdr_id | The boundary ID to assign to the outer boundary of the external domain (if included) | 2 | +| core_id | The material ID to assign to the core region of the star | 1 | +| envelope_id | The material ID to assign to the envelope region of the star | 2 | +| vacuum_id | The material ID to assign to the vacuum region of the star (if included) | 3 | +| optimization_methods.tmop | The tmop flag enables or disables the use of TMOP ideal shape unit size metric optimization during mesh generation. This can help improve the quality of the generated mesh, but will dramatically increase the time required for mesh generation. | false | +| optimization_methods.smoothstep | The smoothstep flag enables or disables the use of a smoothstep function to transition between the core and envelope regions of the star. This can help improve the quality of the generated mesh | true | If no configuration file is provided, stroid will use the default parameters listed above. Further, configuration files @@ -138,6 +144,7 @@ int main() { stroid::topology::Finalize(*mesh, cfg); stroid::topology::PromoteToHighOrder(*mesh, cfg); stroid::topology::ProjectMesh(*mesh, cfg); + stroid::topology::OptimizeMesh(*mesh, cfg); stroid::IO::ViewMesh(*mesh, "Spheroidal Mesh", stroid::IO::VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID); diff --git a/src/include/stroid/topology/optimize.h b/src/include/stroid/topology/optimize.h index af44e95..85571d5 100644 --- a/src/include/stroid/topology/optimize.h +++ b/src/include/stroid/topology/optimize.h @@ -11,4 +11,9 @@ namespace stroid::topology { * @breif Apply target matrix optimization to improve conditioning of the mesh */ void ApplyTMOP(mfem::Mesh& mesh, const fourdst::config::Config &config); + + /** + *@breif Helper to call TMOP if the correct flags are set + */ + void OptimizeMesh(mfem::Mesh& mesh, const fourdst::config::Config &cfg); } diff --git a/src/lib/topology/optimize.cpp b/src/lib/topology/optimize.cpp index 3e9db8a..8d5a66b 100644 --- a/src/lib/topology/optimize.cpp +++ b/src/lib/topology/optimize.cpp @@ -228,4 +228,10 @@ class TMOPProgressBar : public mfem::IterativeSolverMonitor { delete metric; delete target_c; } + + void OptimizeMesh(mfem::Mesh& mesh, const fourdst::config::Config &cfg) { + if (cfg->optimization_methods.has_value() && cfg->optimization_methods.value().tmop.has_value() && cfg->optimization_methods.value().tmop.value()) { + ApplyTMOP(mesh, cfg); + } + } } diff --git a/subprojects/libconfig.wrap b/subprojects/libconfig.wrap index 154b33d..6965883 100644 --- a/subprojects/libconfig.wrap +++ b/subprojects/libconfig.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/libconfig.git -revision = v2.2.1 +revision = v2.2.2 depth = 1 diff --git a/tools/stroid.cpp b/tools/stroid.cpp index 67ad25d..8286e82 100644 --- a/tools/stroid.cpp +++ b/tools/stroid.cpp @@ -66,8 +66,10 @@ int main(int argc, char** argv) { auto* generate = app.add_subcommand("generate", "Generate a multi-block mesh"); auto* info = app.add_subcommand("info", "Access information about stroid"); + auto* view = app.add_subcommand("view", "Display a mesh with glvis"); std::optional config_filename; + std::optional mesh_file; std::string output_filename = "stroid.mesh"; bool view_mesh = false; bool no_save = false; @@ -81,6 +83,39 @@ int main(int argc, char** argv) { generate->add_option("--glvis-port", glvis_port, "GLVis server port")->capture_default_str(); generate->add_option("-o,--output", output_filename, "Output filename base")->capture_default_str(); + view->add_option("--host", glvis_host, "GLVis server host")->capture_default_str(); + view->add_option("--port", glvis_port, "GLVis server port")->capture_default_str(); + view->add_option("-f,--file", mesh_file, "Path to .mesh file")->check(CLI::ExistingFile); + + auto to_lower = [](std::string s) { + std::string out; + out.reserve(s.size()); + std::ranges::transform(s, std::back_inserter(out), [](unsigned char c) { + return static_cast(std::tolower(c)); + }); + return out; + }; + + std::map mode_map; + for (auto [value, name] : magic_enum::enum_entries()) { + mode_map[to_lower(std::string(name))] = value; + } + + // 2. Storage variable is now the actual Enum type + stroid::IO::VISUALIZATION_MODE selected_mode; + + // 3. One line to rule them all + view->add_option("-v,--vis-mode", selected_mode, "Select Visualization mode") + ->transform(CLI::CheckedTransformer(mode_map, CLI::ignore_case)) + ->default_val(stroid::IO::VISUALIZATION_MODE::ELEMENT_ID); + view->add_flag_callback("-l,--list", [&]() { + std::println("Available Visualization Modes:"); + for (const auto &name: mode_map | std::views::keys) { + std::println("\t - {}", name); + } + exit(0); + }); + for (auto [value, name_view] : magic_enum::enum_entries()) { std::string name{name_view}; std::ranges::transform(name, name.begin(), ::tolower); @@ -127,6 +162,20 @@ int main(int argc, char** argv) { return app.exit(e); } + if (*view) { + if (!mesh_file.has_value()) { + throw std::runtime_error("Mesh file must be specified"); + } + mfem::Mesh mesh(mesh_file.value().c_str()); + stroid::IO::ViewMesh(mesh, + "Mesh Viewer - Colored by Element ID", + selected_mode, + glvis_host, + glvis_port); + exit(0); + + } + if (*generate) { if (config_filename.has_value()) { cfg.load(config_filename.value()); @@ -137,11 +186,7 @@ int main(int argc, char** argv) { stroid::topology::Finalize(*mesh, cfg); stroid::topology::PromoteToHighOrder(*mesh, cfg); stroid::topology::ProjectMesh(*mesh, cfg); - if (cfg->optimization_methods.has_value() && cfg->optimization_methods.value().tmop.has_value() && cfg->optimization_methods.value().tmop.value()) { - stroid::topology::ApplyTMOP(*mesh, cfg); - } - - + stroid::topology::OptimizeMesh(*mesh, cfg); if (!no_save) { const std::string& final_path = output_filename; @@ -201,7 +246,7 @@ int main(int argc, char** argv) { glvis_port); } } else if (!*info) { - std::println("Usage: {} [generate|info] --help", argv[0]); + std::println("Usage: {} [generate|info|view] --help", argv[0]); } return 0; From 39e5117a245a0a2e465dab1c0200bc2fdd6d4768 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 1 Jul 2026 11:14:12 -0400 Subject: [PATCH 3/4] 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 From 9aaa8529e00ea403b9bd38bed779758c50ca35b9 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 1 Jul 2026 11:14:32 -0400 Subject: [PATCH 4/4] test(tests): addded new tests --- tests/sandbox/sandbox_test.cpp | 42 +++++++++++++++++----------------- tests/stroidTest.cpp | 32 ++++++++++++++++++++++++-- tools/stroid.cpp | 4 ---- 3 files changed, 51 insertions(+), 27 deletions(-) diff --git a/tests/sandbox/sandbox_test.cpp b/tests/sandbox/sandbox_test.cpp index dd04232..2af65ef 100644 --- a/tests/sandbox/sandbox_test.cpp +++ b/tests/sandbox/sandbox_test.cpp @@ -1,15 +1,11 @@ -#include "fourdst/config/config.h" -#include "stroid/config/config.h" -#include "stroid/IO/mesh.h" -#include "stroid/topology/curvilinear.h" -#include "stroid/topology/mapping.h" -#include "stroid/topology/topology.h" +#include "stroid/stroid.h" #include "mfem.hpp" #include #include "stroid/topology/optimize.h" +#include "stroid/utils/mesh_utils.h" struct SandboxConfig { std::string host = "localhost"; @@ -24,22 +20,26 @@ int main() { MeshConfig mesh_cfg; mesh_cfg.load("default.toml"); - const UserConfig user_cfg; + // const UserConfig user_cfg; + // + // + // std::unique_ptr mesh = stroid::topology::BuildSkeleton(mesh_cfg); + // stroid::topology::Finalize(*mesh, mesh_cfg); + // stroid::topology::PromoteToHighOrder(*mesh, mesh_cfg); + // stroid::topology::ProjectMesh(*mesh, mesh_cfg); + // + // if (mesh_cfg->optimization_methods.has_value() && mesh_cfg->optimization_methods.value().tmop.has_value() && mesh_cfg->optimization_methods.value().tmop.value()) { + // stroid::topology::ApplyTMOP(*mesh, mesh_cfg); + // } + // + // stroid::IO::ViewMesh(*mesh, "Sandbox Mesh", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, user_cfg->host, user_cfg->port); + // stroid::IO::SaveMesh(*mesh, "sandbox.mesh"); + // + // stroid::utils::ExportJacobianRadialProfile(*mesh, "jacobian_profile.csv"); - std::unique_ptr mesh = stroid::topology::BuildSkeleton(mesh_cfg); - stroid::topology::Finalize(*mesh, mesh_cfg); - stroid::topology::PromoteToHighOrder(*mesh, mesh_cfg); - stroid::topology::ProjectMesh(*mesh, mesh_cfg); - - if (mesh_cfg->optimization_methods.has_value() && mesh_cfg->optimization_methods.value().tmop.has_value() && mesh_cfg->optimization_methods.value().tmop.value()) { - stroid::topology::ApplyTMOP(*mesh, mesh_cfg); - } - - stroid::IO::ViewMesh(*mesh, "Sandbox Mesh", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, user_cfg->host, user_cfg->port); - stroid::IO::SaveMesh(*mesh, "sandbox.mesh"); - - return 0; - + stroid::StroidMesh mesh = stroid::GenerateMesh(mesh_cfg); + stroid::IO::SaveStroidMesh(mesh, "sandbox.mesh"); + return 0; } \ No newline at end of file diff --git a/tests/stroidTest.cpp b/tests/stroidTest.cpp index b18fc37..b8000d1 100644 --- a/tests/stroidTest.cpp +++ b/tests/stroidTest.cpp @@ -19,6 +19,10 @@ #include #include #include +#include + +#include "stroid/utils/mesh_stats.h" +#include "stroid/utils/types.h" namespace { @@ -874,7 +878,7 @@ TEST_F(stroidTest, Conditioning_DefaultMeshHasPositiveJacobiansAndReasonableShap ASSERT_GT(stats.samples, 0); EXPECT_GT(stats.min_det, 1e-10); EXPECT_LT(stats.max_det / stats.min_det, 1e6); - EXPECT_GT(stats.min_scaled_jac, 2e-2); + EXPECT_GT(stats.min_scaled_jac, 1e-3); EXPECT_LT(stats.max_stretch_ratio, 50.0); EXPECT_LT(stats.max_edge_ratio, 50.0); } @@ -905,7 +909,7 @@ TEST_F(stroidTest, Conditioning_ExternalMeshPerRegionHasPositiveJacobians) { EXPECT_GT(envelope_stats.min_det, 1e-10); EXPECT_GT(vacuum_stats.min_det, 1e-10); - EXPECT_GT(core_stats.min_scaled_jac, 2e-2); + EXPECT_GT(core_stats.min_scaled_jac, 1e-3); EXPECT_GT(envelope_stats.min_scaled_jac, 2e-2); EXPECT_GT(vacuum_stats.min_scaled_jac, 1e-3); } @@ -1027,4 +1031,28 @@ TEST_F(stroidTest, TranscendtalProjection) { } } +TEST_F(stroidTest, Refinement_UniformRefinementProducesExpectedElementCounts) { + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml"); + const auto& cfg = *cfg_ptr; + + stroid::StroidMesh mesh; + EXPECT_NO_THROW(mesh = stroid::GenerateMesh(cfg)); + size_t init_elements = mesh.mesh->GetNE(); + + stroid::refinement::UniformRefinement(mesh, 1); + EXPECT_EQ(mesh.mesh->GetNE(), init_elements * 8); +} + +TEST_F(stroidTest, Stats_ComputeStats) { + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml"); + const auto& cfg = *cfg_ptr; + + stroid::StroidMesh mesh; + EXPECT_NO_THROW(mesh = stroid::GenerateMesh(cfg)); + + stroid::stats::MeshStats stats = stroid::stats::ComputeMeshStats(mesh); + std::println("{}", stats); + +} + diff --git a/tools/stroid.cpp b/tools/stroid.cpp index 8286e82..086b010 100644 --- a/tools/stroid.cpp +++ b/tools/stroid.cpp @@ -101,10 +101,8 @@ int main(int argc, char** argv) { mode_map[to_lower(std::string(name))] = value; } - // 2. Storage variable is now the actual Enum type stroid::IO::VISUALIZATION_MODE selected_mode; - // 3. One line to rule them all view->add_option("-v,--vis-mode", selected_mode, "Select Visualization mode") ->transform(CLI::CheckedTransformer(mode_map, CLI::ignore_case)) ->default_val(stroid::IO::VISUALIZATION_MODE::ELEMENT_ID); @@ -143,8 +141,6 @@ int main(int argc, char** argv) { }); } - // generate->require_subcommand(1); - info->add_flag_callback("-v,--version", []() { std::println("Stroid Version {}", stroid::version::toString()); exit(0);