From 5a82311251e716b28467c6aacd19913ee54350ae Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 7 Apr 2026 12:19:58 -0400 Subject: [PATCH] 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;