feat(src): Enabled Vacuume Region Generation

Stroid can now generate vacuume regions. These are confirming spheres marked with a seperate attribute that solver may then apply arbitrary mappings to. Further, the boundary of this region is also marked with a seperate boundary attrite allowing seperate boundary conditions to easily be applied.

configs/test_external_domain_refinement_l1.toml

@@ -0,0 +1,11 @@
+[main]
+core_steepness = 1.0
+flattening = 0.0
+include_external_domain = true
+order = 3
+r_core = 1.5
+r_infinity = 6.0
+r_instability = 1e-14
+r_star = 5.0
+refinement_levels = 1
+

configs/test_refinement_l2.toml

@@ -0,0 +1,11 @@
+[main]
+core_steepness = 1.0
+flattening = 0.0
+include_external_domain = false
+order = 3
+r_core = 1.5
+r_infinity = 6.0
+r_instability = 1e-14
+r_star = 5.0
+refinement_levels = 2
+

configs/test_volume_no_external.toml

@@ -0,0 +1,16 @@
+[main]
+refinement_levels = 2
+order = 3
+include_external_domain = false
+r_core = 1.5
+r_star = 5.0
+flattening = 0.08
+r_infinity = 6.0
+r_instability = 1e-14
+core_steepness = 1.0
+surface_bdr_id = 1
+inf_bdr_id = 2
+core_id = 1
+envelope_id = 2
+vacuum_id = 3
+

configs/test_volume_spherical_no_external.toml

@@ -0,0 +1,16 @@
+[main]
+refinement_levels = 2
+order = 3
+include_external_domain = false
+r_core = 1.5
+r_star = 5.0
+flattening = 0.0
+r_infinity = 6.0
+r_instability = 1e-14
+core_steepness = 1.0
+surface_bdr_id = 1
+inf_bdr_id = 2
+core_id = 1
+envelope_id = 2
+vacuum_id = 3
+

configs/test_volume_spherical_with_external.toml

@@ -0,0 +1,16 @@
+[main]
+refinement_levels = 2
+order = 3
+include_external_domain = true
+r_core = 1.5
+r_star = 5.0
+flattening = 0.0
+r_infinity = 6.0
+r_instability = 1e-14
+core_steepness = 1.0
+surface_bdr_id = 1
+inf_bdr_id = 2
+core_id = 1
+envelope_id = 2
+vacuum_id = 3
+

configs/test_volume_with_external.toml

@@ -0,0 +1,16 @@
+[main]
+refinement_levels = 2
+order = 3
+include_external_domain = true
+r_core = 1.5
+r_star = 5.0
+flattening = 0.08
+r_infinity = 6.0
+r_instability = 1e-14
+core_steepness = 1.0
+surface_bdr_id = 1
+inf_bdr_id = 2
+core_id = 1
+envelope_id = 2
+vacuum_id = 3
+

src/include/stroid/config/config.h

@@ -12,52 +12,95 @@ namespace stroid::config {
     struct MeshConfig {
         /**
          * @brief Number of uniform refinement passes applied after topology creation.
-         * @toml [main].refinement_levels
+         * @section toml
+         * - [main].refinement_levels
          */
         int refinement_levels = 4;
         /**
          * @brief Polynomial order for high-order elements.
-         * @toml [main].order
+         * @section toml
+         * - [main].order
          */
         int order = 3;
         /**
          * @brief Whether to include an external domain extending to `r_infinity`.
-         * @details Currently this flag does not affect mesh generation.
-         * @toml [main].include_external_domain
+         * @section toml
+         * - [main].include_external_domain
          */
-        bool include_external_domain = false;
+        bool include_external_domain = true;
 
         /**
          * @brief Radius of the stellar core region.
-         * @toml [main].r_core
+         * @section toml
+         * - [main].r_core
          */
         double r_core = 1.5;
         /**
          * @brief Radius of the stellar surface.
-         * @toml [main].r_star
+         * @section toml
+         * - [main].r_star
          */
         double r_star = 5.0;
         /**
          * @brief Flattening factor for spheroidal shaping (0 = spherical, >0 = oblate).
-         * @toml [main].flattening
+         * @section toml
+         * - [main].flattening
          */
         double flattening = 0;
 
         /**
          * @brief Outer radius of the external domain when enabled.
-         * @toml [main].r_infinity
+         * @section toml
+         * - [main].r_infinity
          */
         double r_infinity = 6.0;
 
         /**
          * @brief Radius inside which transformations are skipped to avoid singularities.
-         * @toml [main].r_instability
+         * @section toml
+         * - [main].r_instability
          */
         double r_instability = 1e-14;
         /**
          * @brief Controls the smoothness/steepness of the core-to-envelope transition.
-         * @toml [main].core_steepness
+         * @section toml
+         * - [main].core_steepness
          */
         double core_steepness = 1.0;
+
+        /**
+         * @brief Boundary attribute id for stellar surface
+         * @section toml
+         * - [main].surface_bdr_id
+         */
+        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;
+
+        /**
+         * @brief Material attribute id for the core region
+         * @section toml
+         * - [main].core_id
+         */
+        size_t core_id = 1;
+
+        /**
+         * @brief Material attribute id for the envelope region
+         * @section toml
+         * - [main].envelope_id
+         */
+        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;
     };
 }

src/lib/topology/curvilinear.cpp

@@ -3,6 +3,7 @@
 
 #include <iostream>
 #include <memory>
+#include <sys/proc.h>
 
 namespace stroid::topology {
     void PromoteToHighOrder(mfem::Mesh &mesh, const fourdst::config::Config<config::MeshConfig> &config) {
@@ -22,20 +23,48 @@ namespace stroid::topology {
 
         const int vDim = fes->GetVDim();
         const int nDofs = fes->GetNDofs();
+        const int nElem = mesh.GetNE();
 
+        std::vector<bool> processed(nDofs, false);
+        mfem::Array<int> vdofs;
         mfem::Vector pos(vDim);
 
-        for (int i = 0; i < nDofs; ++i) {
-            for (int d = 0; d < vDim; ++d) {
-                pos(d) = nodes(fes->DofToVDof(i, d));
-            }
+        for (int elemID = 0; elemID < nElem; ++elemID) {
+            const int attrID = mesh.GetAttribute(elemID);
+            fes->GetElementVDofs(elemID, vdofs);
 
-            TransformPoint(pos, config, 0);
+            for (int dofID = 0; dofID < vdofs.Size(); ++dofID) {
+                const int vDof = vdofs[dofID];
+                const int scalar_dof = (fes->GetOrdering() == mfem::Ordering::byNODES) ? vDof / vDim : vDof % nDofs;
 
-            for (int d = 0; d < vDim; ++d) {
-                nodes(fes->DofToVDof(i, d)) = pos(d);
+                if (processed[scalar_dof]) {
+                    continue; // Skip already processed dofs. This avoids doing multiple transformations of a node if it was already transformed by a neighbor
+                }
+
+                for (int d = 0; d < vDim; ++d) {
+                    pos(d) = nodes(fes->DofToVDof(scalar_dof, d));
+                }
+
+                TransformPoint(pos, config, attrID);
+
+                for (int d = 0; d < vDim; ++d) {
+                    nodes(fes->DofToVDof(scalar_dof, d)) = pos(d);
+                }
+
+                processed[scalar_dof] = true;
             }
         }
 
+        // 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

@@ -32,20 +32,6 @@ namespace stroid::topology {
         pos(2) *= (1.0 - config->flattening);
     }
 
-    void ApplyKelvin(mfem::Vector &pos, const fourdst::config::Config<config::MeshConfig> &config) {
-        const double r = pos.Norml2();
-        if (r <= config->r_star) {
-            return;
-        }
-
-        double xi = (r - config->r_star) / (config->r_infinity - config->r_star);
-        xi = std::min(0.999, std::max(0.0, xi)); // Clamp xi to [0, 0.999]
-
-        const double r_new = config->r_star + xi / (1.0 - xi);
-        const double scale = r_new / r;
-        pos *= scale;
-    }
-
     void TransformPoint(mfem::Vector &pos, const fourdst::config::Config<config::MeshConfig> &config, int attribute_id) {
         double l_inf = 0.0;
         for (int i = 0; i < pos.Size(); ++i) {
@@ -82,9 +68,6 @@ namespace stroid::topology {
             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);
@@ -97,7 +80,6 @@ namespace stroid::topology {
             pos = unit_dir;
             pos *= l_inf;
 
-            ApplyKelvin(pos, config);
             ApplySpheroidal(pos, config);
         }
     }

src/lib/topology/topology.cpp

@@ -10,7 +10,7 @@ namespace stroid::topology {
     std::unique_ptr<mfem::Mesh> BuildSkeleton(const fourdst::config::Config<config::MeshConfig> & config) {
         int nVert = config->include_external_domain ? 24 : 16;
         int nElem = config->include_external_domain ? 13 : 7;
-        int nBev  = 6;
+        int nBev  = config->include_external_domain ? 12 : 6;
 
         auto mesh = std::make_unique<mfem::Mesh>(3, nVert, nElem, nBev, 3);
 
@@ -28,9 +28,9 @@ namespace stroid::topology {
         }
 
         const int core_v[8] = {0, 1, 3, 2, 4, 5, 7, 6};
-        mesh->AddHex(core_v, 1);
+        mesh->AddHex(core_v, config->core_id);
 
-        int shells[6][8] = {
+        std::vector<std::array<int, 8>> stellar_shells = {
             {8, 9, 11, 10, 0, 1, 3, 2},
             {4, 5, 7, 6, 12, 13, 15, 14}, // +Z face
             {0, 1, 5, 4, 8, 9, 13, 12},   // -Y face
@@ -38,9 +38,25 @@ namespace stroid::topology {
             {1, 3, 7, 5, 9, 11, 15, 13},  // +X face
             {0, 4, 6, 2, 8, 12, 14, 10}   // -X face
         };
-        for (const auto & shell : shells) mesh->AddHex(shell, 2);
+        for (const auto & shell : stellar_shells) {
+            mesh->AddHex(shell.data(), config->envelope_id);
+        }
 
-        const int bdr_quads[6][4] = {
+        if (config->include_external_domain) {
+            std::vector<std::array<int, 8>> vacuum_shells;
+            vacuum_shells.push_back({8, 9, 13, 12, 16, 17, 21, 20});
+            vacuum_shells.push_back({9, 11, 15, 13, 17, 19, 23, 21});
+            vacuum_shells.push_back({11, 10, 14, 15, 19, 18, 22, 23});
+            vacuum_shells.push_back({10, 8, 12, 14, 18, 16, 20, 22});
+            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);
+            }
+        }
+
+
+        const int surface_bdr_quads[6][4] = {
             {12, 13, 15, 14},
             {13, 9, 11, 15},
             {9, 8, 10, 11},
@@ -49,8 +65,23 @@ namespace stroid::topology {
             {14, 15, 11, 10}
         };
 
-        for (const auto& bdr: bdr_quads) {
-            mesh->AddBdrQuad(bdr, 1);
+        for (const auto& bdr: surface_bdr_quads) {
+            mesh->AddBdrQuad(bdr, config->surface_bdr_id);
+        }
+
+        if (config->include_external_domain) {
+            const int inf_bdr_quads[6][4] = {
+                {16, 17, 21, 20},
+                {17, 19, 23, 21},
+                {19, 18, 22, 23},
+                {18, 16, 20, 22},
+                {18, 19, 17, 16},
+                {20, 21, 23, 22}
+            };
+
+            for (const auto& bdr: inf_bdr_quads) {
+                mesh->AddBdrQuad(bdr, config->inf_bdr_id);
+            }
         }
 
         return mesh;

src/lib/utils/mesh_utils.cpp

@@ -4,8 +4,6 @@
 
 namespace stroid::utils {
     void MarkFlippedElements(mfem::Mesh& mesh) {
-        size_t total_flipped = 0;
-        size_t total_elements = mesh.GetNE();
         for (int i = 0; i < mesh.GetNE(); i++) {
             mfem::ElementTransformation *T = mesh.GetElementTransformation(i);
 
@@ -22,17 +20,11 @@ namespace stroid::utils {
 
             if (is_flipped) {
                 mesh.SetAttribute(i, 999);
-                total_flipped++;
             }
-
-
         }
-        std::println("Marked {}/{} elements as flipped.", total_flipped, total_elements);
     }
 
     void MarkFlippedBoundaryElements(mfem::Mesh& mesh) {
-        size_t total_flipped = 0;
-        size_t total_boundary_elements = mesh.GetNBE();
         for (int i = 0; i < mesh.GetNBE(); i++) {
             mfem::ElementTransformation *T = mesh.GetBdrElementTransformation(i);
             const mfem::IntegrationRule &ir = mfem::IntRules.Get(T->GetGeometryType(), 2 * T->Order());
@@ -56,9 +48,7 @@ namespace stroid::utils {
 
             if (is_flipped) {
                 mesh.SetBdrAttribute(i, 500);
-                total_flipped++;
             }
         }
-        std::println("Marked {}/{} boundary elements as flipped.", total_flipped, total_boundary_elements);
     }
 }