diff --git a/configs/test_external_domain_refinement_l1.toml b/configs/test_external_domain_refinement_l1.toml new file mode 100644 index 0000000..37bad7b --- /dev/null +++ b/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 + diff --git a/configs/test_refinement_l2.toml b/configs/test_refinement_l2.toml new file mode 100644 index 0000000..f270e7d --- /dev/null +++ b/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 + diff --git a/configs/test_volume_no_external.toml b/configs/test_volume_no_external.toml new file mode 100644 index 0000000..cf86efc --- /dev/null +++ b/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 + diff --git a/configs/test_volume_spherical_no_external.toml b/configs/test_volume_spherical_no_external.toml new file mode 100644 index 0000000..456ed85 --- /dev/null +++ b/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 + diff --git a/configs/test_volume_spherical_with_external.toml b/configs/test_volume_spherical_with_external.toml new file mode 100644 index 0000000..10bc77e --- /dev/null +++ b/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 + diff --git a/configs/test_volume_with_external.toml b/configs/test_volume_with_external.toml new file mode 100644 index 0000000..617a2c7 --- /dev/null +++ b/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 + diff --git a/src/include/stroid/config/config.h b/src/include/stroid/config/config.h index 083a331..416956a 100644 --- a/src/include/stroid/config/config.h +++ b/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; }; } diff --git a/src/lib/topology/curvilinear.cpp b/src/lib/topology/curvilinear.cpp index 3a98151..ae0150c 100644 --- a/src/lib/topology/curvilinear.cpp +++ b/src/lib/topology/curvilinear.cpp @@ -3,6 +3,7 @@ #include #include +#include namespace stroid::topology { void PromoteToHighOrder(mfem::Mesh &mesh, const fourdst::config::Config &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 processed(nDofs, false); + mfem::Array 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); + // } + // } } } diff --git a/src/lib/topology/mapping.cpp b/src/lib/topology/mapping.cpp index fba882a..bbc9105 100644 --- a/src/lib/topology/mapping.cpp +++ b/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) { - 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, 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); } } diff --git a/src/lib/topology/topology.cpp b/src/lib/topology/topology.cpp index db0cff8..ddc4286 100644 --- a/src/lib/topology/topology.cpp +++ b/src/lib/topology/topology.cpp @@ -10,7 +10,7 @@ namespace stroid::topology { std::unique_ptr BuildSkeleton(const fourdst::config::Config & 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(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> 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> 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; diff --git a/src/lib/utils/mesh_utils.cpp b/src/lib/utils/mesh_utils.cpp index 7c94f53..9c691f2 100644 --- a/src/lib/utils/mesh_utils.cpp +++ b/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); } } \ No newline at end of file diff --git a/tests/sandbox/meson.build b/tests/sandbox/meson.build new file mode 100644 index 0000000..e69de29 diff --git a/tests/sandbox/sandbox_test.cpp b/tests/sandbox/sandbox_test.cpp new file mode 100644 index 0000000..e69de29