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

src/include/stroid/config/config.h

@@ -1,6 +1,14 @@
 #pragma once
 
+#include <cstdint>
+
 namespace stroid::config {
+
+    struct OptimizationMethods {
+        std::optional<bool> tmop{false};
+        std::optional<bool> 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<int> refinement_levels = 4;
         /**
          * @brief Polynomial order for high-order elements.
          * @section toml
          * - [main].order
          */
-        int order = 3;
+        std::optional<int> 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<bool> include_external_domain = true;
 
         /**
          * @brief Radius of the stellar core region.
          * @section toml
          * - [main].r_core
          */
-        double r_core = 1.5;
+        std::optional<double> r_core = 0.25;
         /**
          * @brief Radius of the stellar surface.
          * @section toml
          * - [main].r_star
          */
-        double r_star = 5.0;
+        std::optional<double> r_star = 1.0;
         /**
          * @brief Flattening factor for spheroidal shaping (0 = spherical, >0 = oblate).
          * @section toml
          * - [main].flattening
          */
-        double flattening = 0;
+        std::optional<double> flattening = 0;
 
         /**
          * @brief Outer radius of the external domain when enabled.
          * @section toml
          * - [main].r_infinity
          */
-        double r_infinity = 6.0;
+        std::optional<double> 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<double> 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<double> 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<size_t> 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<size_t> inf_bdr_id = 2;
 
         /**
          * @brief Material attribute id for the core region
          * @section toml
          * - [main].core_id
          */
-        size_t core_id = 1;
+        std::optional<size_t> core_id = 1;
 
         /**
          * @brief Material attribute id for the envelope region
          * @section toml
          * - [main].envelope_id
          */
-        size_t envelope_id = 2;
+        std::optional<size_t> 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<size_t> vacuum_id = 3;
+
+        std::optional<OptimizationMethods> optimization_methods = OptimizationMethods{true, true};
     };
 }

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"
 

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::MeshConfig> &config);
+}

src/lib/topology/curvilinear.cpp

@@ -5,7 +5,7 @@
 
 namespace stroid::topology {
     void PromoteToHighOrder(mfem::Mesh &mesh, const fourdst::config::Config<config::MeshConfig> &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);
-        //     }
-        // }
     }
 }

src/lib/topology/mapping.cpp

@@ -29,7 +29,7 @@ namespace stroid::topology {
     }
 
     void ApplySpheroidal(mfem::Vector &pos, const fourdst::config::Config<config::MeshConfig> &config) {
-        pos(2) *= (1.0 - config->flattening);
+        pos(2) *= (1.0 - config->flattening.value());
     }
 
     void TransformPoint(mfem::Vector &pos, const fourdst::config::Config<config::MeshConfig> &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);
         }
-    }
-}
+    }}

src/lib/topology/optimize.cpp

@@ -0,0 +1,231 @@
+#include "mfem.hpp"
+
+#include <thread>
+#include <atomic>
+#include <chrono>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <algorithm>
+
+#include "stroid/topology/optimize.h"
+
+
+#include <clocale>
+#include <cstdlib>
+#include <cstring>
+#include <string>
+
+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<bool> done_{false};
+        std::atomic<double> progress_{0.0};
+        std::atomic<int> iter_{0};
+        std::atomic<double> res_{0.0};
+
+        std::thread spinner_thread_;
+
+        void Spin() {
+            std::vector<std::string> 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<int>(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<int>(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::MeshConfig> &config) {
+        const mfem::FiniteElementSpace* cfes = mesh.GetNodalFESpace();
+        mfem::FiniteElementSpace* fes = const_cast<mfem::FiniteElementSpace*>(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<int> 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<int> 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<double>::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;
+    }
+}

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<std::array<int, 8>> 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());
             }
         }
 

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',
 )