diff --git a/Doxyfile b/Doxyfile index 757d19f..abcc940 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.3.0 +PROJECT_NUMBER = v0.5.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/assets/imgs/ExampleMesh_multi-block.png b/assets/imgs/ExampleMesh_multi-block.png new file mode 100644 index 0000000..59d523c Binary files /dev/null and b/assets/imgs/ExampleMesh_multi-block.png differ diff --git a/assets/imgs/ExampleMesh_spherified.png b/assets/imgs/ExampleMesh_spherified.png new file mode 100644 index 0000000..d07c0b4 Binary files /dev/null and b/assets/imgs/ExampleMesh_spherified.png differ diff --git a/build-config/mfem/meson.build b/build-config/mfem/meson.build index 4a6c646..6476093 100644 --- a/build-config/mfem/meson.build +++ b/build-config/mfem/meson.build @@ -1,16 +1,22 @@ -cmake = import('cmake') -mfem_cmake_options = cmake.subproject_options() -mfem_cmake_options.add_cmake_defines({ - 'MFEM_ENABLE_EXAMPLES': 'OFF', - 'MFEM_ENABLE_TESTING': 'OFF', - 'MFEM_ENABLE_MINIAPPS': 'OFF', - 'MFEM_USE_BENCMARK': 'OFF', - 'BUILD_SHARED_LIBS': 'OFF', - 'BUILD_STATIC_LIBS': 'ON', -}) -mfem_cmake_options.set_install(true) +mfem_dep = dependency('mfem', required : false) -mfem_sp = cmake.subproject( - 'mfem', - options: mfem_cmake_options) -mfem_dep = mfem_sp.dependency('mfem') +if not mfem_dep.found() + cmake = import('cmake') + mfem_cmake_options = cmake.subproject_options() + mfem_cmake_options.add_cmake_defines({ + 'MFEM_ENABLE_EXAMPLES': 'OFF', + 'MFEM_ENABLE_TESTING': 'OFF', + 'MFEM_ENABLE_MINIAPPS': 'OFF', + 'MFEM_USE_BENCMARK': 'OFF', + 'BUILD_SHARED_LIBS': 'OFF', + 'BUILD_STATIC_LIBS': 'ON', + }) + mfem_cmake_options.set_install(true) + + mfem_sp = cmake.subproject( + 'mfem', + options: mfem_cmake_options) + mfem_dep = mfem_sp.dependency('mfem') +else + message('Using system-installed MFEM library') +endif diff --git a/configs/conditioned_core.toml b/configs/conditioned_core.toml new file mode 100644 index 0000000..8434f89 --- /dev/null +++ b/configs/conditioned_core.toml @@ -0,0 +1,21 @@ +[main] +core_mapping = "multi_block" +refinement_levels = 2 +order = 3 +include_external_domain = true +r_core = 0.25 +r_star = 1.0 +r_infinity = 6.0 +flattening = 0.0 +r_instability = 1e-14 +core_steepness = 1.0 +continuity_order = 2 +surface_bdr_id = 1 +inf_bdr_id = 2 +core_id = 1 +envelope_id = 2 +vacuum_id = 3 + +[main.optimization_methods] +tmop = false +smoothstep = true diff --git a/meson.build b/meson.build index a9c7a03..0c15932 100644 --- a/meson.build +++ b/meson.build @@ -1,4 +1,4 @@ -project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.4.0', default_options : ['cpp_std=c++23']) +project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.5.0', default_options : ['cpp_std=c++23']) subdir('build-check') @@ -28,6 +28,7 @@ if get_option('pkg_config') ], subdirs: ['stroid'], filebase: 'stroid', - install_dir: join_paths(get_option('libdir'), 'pkgconfig') + install_dir: join_paths(get_option('libdir'), 'pkgconfig'), + requires: ['fourdst_config'] ) endif \ No newline at end of file diff --git a/readme.md b/readme.md index ce1545a..58955fb 100644 --- a/readme.md +++ b/readme.md @@ -92,6 +92,8 @@ inf_bdr_id = 2 core_id = 1 envelope_id = 2 vacuum_id = 3 +core_mapping = "multi_block" + [main.optimization_methods] tmop = false @@ -99,29 +101,53 @@ 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 | -| 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 | +| 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 | + | core_mapping | The core mapping strategy to use for the mesh generation. Options are "spherified" (legacy) or "multi_block" (conditioned). The multi_block strategy is strongly preferred for its improved condition number. | "multi_block" | If no configuration file is provided, stroid will use the default parameters listed above. Further, configuration files need only include parameters that differ from the defaults, any parameters not specified will use the default values. +### Conditioned core mapping + +There are two core mapping strategies, spherified and multi_block. Generally multi_block should be strongly preferred. The +`core_mapping = "multi_block"` strategy avoids the radial rank loss at the eight corners of the spherified core +block. It uses a Cartesian center plus six transition blocks inside the core. The inner cube has circumscribed radius +`r_core / 2`; its six faces connect linearly to the existing spherical `r_core` interface. If enabled, spheroidal flattening is +applied afterwards. + +```python +cfg = stroid.config.MeshConfig(core_mapping="multi_block", refinement_levels=2) +cfg.optimization_methods = stroid.config.OptimizationMethods(tmop=False) +mesh = stroid.GenerateMesh(cfg) +``` + +The optional, non-installed `geometry_quality_experiment` target may be used to measure the actual high-order geometry +at quadrature points, vertices, edges, and near-corner probes. You may build and run it explicitly: + +```bash +meson compile -C build geometry_quality_experiment +build/tools/geometry_quality_experiment --orders 4 --refinements 2 \ + --contraction-probe --probe-order 3 --output core_comparison.csv +``` + ### C++ Interface Stroid can be used as a library in C++ projects. After installation, include the stroid header and link against the stroid library. @@ -153,7 +179,14 @@ int main() { ## Example Meshes An example mesh with the default configuration parameters is shown below (coloration indicates attribute IDs of different regions): -![Example Mesh](assets/imgs/ExampleMesh.png) +![Example Mesh](assets/imgs/ExampleMesh_multi-block.png) + +The legacy spherified core mapping strategy is shown below as well +![Example Spheried Mesh](assets/imgs/ExampleMesh_spherified.png) + +Note that both of these meshes are shown with 3 levels of refinement and polynomial order 3. Blue shows the stellar +domain while purple shows the vacuum domain. + ## Funding Stroid is developed as part of the 4D-STAR project. @@ -161,4 +194,4 @@ Stroid is developed as part of the 4D-STAR project. 4D-STAR is funded by European Research Council (ERC) under the Horizon Europe programme (Synergy Grant agreement No. 101071505: 4D-STAR) Work for this project is funded by the European Union. Views and opinions expressed are however those of the author(s) -only and do not necessarily reflect those of the European Union or the European Research Council. \ No newline at end of file +only and do not necessarily reflect those of the European Union or the European Research Council. diff --git a/src/include/stroid.zip b/src/include/stroid.zip new file mode 100644 index 0000000..eddd818 Binary files /dev/null and b/src/include/stroid.zip differ diff --git a/src/include/stroid/config/config.h b/src/include/stroid/config/config.h index 918245a..999576f 100644 --- a/src/include/stroid/config/config.h +++ b/src/include/stroid/config/config.h @@ -3,6 +3,8 @@ #include #include #include +#include +#include namespace stroid::config { @@ -122,5 +124,60 @@ namespace stroid::config { std::optional optimization_methods = OptimizationMethods{true, true}; + /** + * @brief Core mapping strategy: legacy "spherified" or conditioned "multi_block". + * + * spherified generates a either two or three inscribed cubes then projects them into spheres. + * multi_block generates a multi-block topology with a single core block and six envelope blocks, then projects the core block into a sphere and the envelope blocks into a spheroid. + * + * multi_block is strongly preferred for its ~1000x improved condition number, Spherified is only provided for legacy compatibility. + * + * @section toml + * - [main].core_mapping + */ + std::optional core_mapping = "multi_block"; + }; + + inline std::string to_string(const MeshConfig &mesh_config) { + auto opt_2_string = [](const OptimizationMethods& opt) { + std::stringstream ss; + ss << ""; + return ss.str(); + }; + + std::stringstream ss; + + OptimizationMethods opt = mesh_config.optimization_methods.value_or(OptimizationMethods{false, true}); + std::string opt_string = opt_2_string(opt); + + ss << "MeshConfig:\n"; + ss << std::format(" refinement_levels: {}\n", mesh_config.refinement_levels.value_or(4)); + ss << std::format(" order: {}\n", mesh_config.order.value_or(3)); + ss << std::format(" include_external_domain: {}\n", mesh_config.include_external_domain.value_or(true)); + ss << std::format(" r_core: {}\n", mesh_config.r_core.value_or(0.25)); + ss << std::format(" r_star: {}\n", mesh_config.r_star.value_or(1.0)); + ss << std::format(" flattening: {}\n", mesh_config.flattening.value_or(0.0)); + ss << std::format(" r_infinity: {}\n", mesh_config.r_infinity.value_or(6.0)); + ss << std::format(" r_instability: {}\n", mesh_config.r_instability.value_or(1e-14)); + ss << std::format(" core_steepness: {}\n", mesh_config.core_steepness.value_or(1.0)); + ss << std::format(" continuity_order: {}\n", mesh_config.continuity_order.value_or(2)); + ss << std::format(" surface_bdr_id: {}\n", mesh_config.surface_bdr_id.value_or(1)); + ss << std::format(" inf_bdr_id: {}\n", mesh_config.inf_bdr_id.value_or(2)); + ss << std::format(" core_id: {}\n", mesh_config.core_id.value_or(1)); + ss << std::format(" envelope_id: {}\n", mesh_config.envelope_id.value_or(2)); + ss << std::format(" vacuum_id: {}\n", mesh_config.vacuum_id.value_or(3)); + ss << std::format(" optimization_methods: {}\n", opt_string); + ss << std::format(" core_mapping: {}\n", mesh_config.core_mapping.value_or("spherified")); + + return ss.str(); + + } } diff --git a/src/include/stroid/stroid.h b/src/include/stroid/stroid.h index c9eadef..3864ccc 100644 --- a/src/include/stroid/stroid.h +++ b/src/include/stroid/stroid.h @@ -51,6 +51,7 @@ namespace stroid { inline StroidMesh GenerateMesh(const fourdst::config::Config& cfg) { StroidMesh sm; + sm.type = MFEM_MESH_TYPE::SERIAL; sm.config = *cfg; auto reference = stroid::topology::BuildSkeleton(cfg); stroid::topology::Finalize(*reference, cfg); @@ -61,6 +62,7 @@ namespace stroid { 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); } + sm.exterior_coordinate = stroid::topology::BuildExteriorCoordinate(*sm.mesh, *sm.reference_mesh, cfg); return sm; } inline StroidMesh GenerateMesh(const stroid::config::MeshConfig& config) { diff --git a/src/include/stroid/topology/curvilinear.h b/src/include/stroid/topology/curvilinear.h index 3b8a7da..16adf83 100644 --- a/src/include/stroid/topology/curvilinear.h +++ b/src/include/stroid/topology/curvilinear.h @@ -3,6 +3,7 @@ #include "mfem.hpp" #include "stroid/config/config.h" #include "fourdst/config/config.h" +#include "stroid/utils/types.h" namespace stroid::topology { /** @@ -18,4 +19,17 @@ namespace stroid::topology { * @param config Mesh configuration (uses radii, flattening, and mapping parameters). */ void ProjectMesh(mfem::Mesh& mesh, const fourdst::config::Config &config); + + /** + * @brief Build a scalar grid function representing the compactification coordinate for a mesh. This ranges from 0-1 with 0 at the stellar surface and 1 at the compactified infinity. + * @param mesh Reference to the underlying serial MFEM mesh which has been promoted to high-order and projected into the curvilinear domain. + * @param reference_mesh reference to the underlying serial which has not been promoted to high-order or projected into the curvilinear domain. This is used to compute the compactification coordinate. + * @param config Config file + * @return Unique pointer to a scalar mesh field representing the compactification coordinate. + */ + std::unique_ptr BuildExteriorCoordinate( + mfem::Mesh& mesh, + mfem::Mesh& reference_mesh, + const fourdst::config::Config& config + ); } \ No newline at end of file diff --git a/src/include/stroid/topology/mapping.h b/src/include/stroid/topology/mapping.h index e1f53b6..3761a32 100644 --- a/src/include/stroid/topology/mapping.h +++ b/src/include/stroid/topology/mapping.h @@ -28,8 +28,23 @@ namespace stroid::topology { /** * @brief Map a point from the initial block topology to the curvilinear domain. * @param pos Position vector updated in-place. - * @param config Mesh configuration (uses radii, flattening, instability radius, and core steepness). + * @param config Mesh configuration (uses radii, flattening, and `core_mapping`). + * The `multi_block` strategy requires the matching skeleton from BuildSkeleton; + * changing only the mapping on a legacy core element is not supported. * @param attribute_id Element attribute ID (currently unused). */ void TransformPoint(mfem::Vector& pos, const fourdst::config::Config &config, int attribute_id); -} \ No newline at end of file + + /** + * @brief Compute the compactification coordinate for a point in the curvilinear domain. This ranges from 0-1 with 0 at the stellar surface and 1 at the compactified infinity. + * @param logical_position Logical position of the point in the curvilinear domain. + * @param attribute Element attribute ID (used to determine the exterior coordinate). + * @param config Mesh configuration (uses radii and flattening). + * @return Compactification coordinate ranging from 0 (stellar surface) to 1 (compactified infinity). + */ + double ComputeExteriorCoordinate( + const mfem::Vector& logical_position, + int attribute, + const fourdst::config::Config& config + ); +} diff --git a/src/include/stroid/topology/topology.h b/src/include/stroid/topology/topology.h index c63c887..3327893 100644 --- a/src/include/stroid/topology/topology.h +++ b/src/include/stroid/topology/topology.h @@ -8,7 +8,9 @@ namespace stroid::topology { /** * @brief Build the initial multi-block mesh topology for the star model. - * @param config Mesh configuration (uses radii and domain flags). + * @param config Mesh configuration (uses radii, domain flags, and `core_mapping`). + * The legacy `spherified` core uses one block; `multi_block` uses an + * inner Cartesian block and six core transition blocks. * @return Newly allocated mesh skeleton (not yet refined or curved). */ std::unique_ptr BuildSkeleton(const fourdst::config::Config & config); diff --git a/src/include/stroid/utils/types.h b/src/include/stroid/utils/types.h index 2b95332..412a6d2 100644 --- a/src/include/stroid/utils/types.h +++ b/src/include/stroid/utils/types.h @@ -16,10 +16,16 @@ namespace stroid { PARALLEL }; + struct ScalarMeshField { + std::unique_ptr space; + std::unique_ptr values; + }; + struct StroidMesh { MFEM_MESH_TYPE type; std::unique_ptr mesh; std::unique_ptr reference_mesh; + std::unique_ptr exterior_coordinate; config::MeshConfig config; size_t refinement_levels; @@ -61,5 +67,17 @@ namespace stroid { return mesh_stats; } + + std::unique_ptr clone() const { + std::unique_ptr new_mesh; + new_mesh->type = type; + new_mesh->mesh = std::make_unique(*mesh); + new_mesh->reference_mesh = std::make_unique(*reference_mesh); + new_mesh->config = config; + new_mesh->refinement_levels = refinement_levels; + + return new_mesh; + } + }; } diff --git a/src/lib/IO/mesh.cpp b/src/lib/IO/mesh.cpp index 670c517..285f311 100644 --- a/src/lib/IO/mesh.cpp +++ b/src/lib/IO/mesh.cpp @@ -1,12 +1,16 @@ #include "mfem.hpp" #include "stroid/config/config.h" #include "stroid/IO/mesh.h" +#include "stroid/topology/curvilinear.h" +#include #include +#include #include "stroid/version.h" #include +#include #include #include #include @@ -17,6 +21,7 @@ #include #include #include +#include namespace stroid::IO { @@ -34,6 +39,7 @@ namespace stroid::IO { # - 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 +# - exterior coordinate : a scalar material coordinate which is zero at the stellar surface and one at infinity # - 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 @@ -58,7 +64,7 @@ END BLOCK HEADER)", std::string format_primary_mesh(const StroidMesh& mesh) { std::stringstream ss; - ss.precision(8); + ss.precision(std::numeric_limits::max_digits10); mesh.mesh->Print(ss); std::string pmesh = std::format("BEGIN BLOCK PMESH\n{}END BLOCK PMESH", ss.str()); @@ -77,7 +83,7 @@ END BLOCK HEADER)", std::string format_reference_mesh(const StroidMesh& mesh) { std::stringstream ss; - ss.precision(8); + ss.precision(std::numeric_limits::max_digits10); mesh.reference_mesh->Print(ss); std::string rmesh = std::format("BEGIN BLOCK RMESH\n{}END BLOCK RMESH", ss.str()); @@ -170,6 +176,11 @@ vacuum_id:{} # includes tmop and smoothstep booleans optimization_methods-tmop:{} optimization_methods-smoothstep:{} + +# core_mapping: Core mapping strategy, either spherified or multi_block +# std::optional +# default: spherified +core_mapping:{} END BLOCK CONFIG)", format_opt(mesh.config.refinement_levels, d.refinement_levels.value()), format_opt(mesh.config.order, d.order.value()), @@ -187,10 +198,53 @@ END BLOCK CONFIG)", 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)); + m_opt.smoothstep.value_or(true), + format_opt(mesh.config.core_mapping, d.core_mapping.value())); return config_str; } + + std::string format_exterior_coordinate(const StroidMesh& mesh) { + const bool include_external_domain = mesh.config.include_external_domain.value_or(true); + + if (!include_external_domain) { + if (mesh.exterior_coordinate) throw std::runtime_error("A mesh without an external domain cannot contain an exterior-coordinate field."); + return "BEGIN BLOCK EXTERIOR_COORDINATE\nPRESENT:false\nEND BLOCK EXTERIOR_COORDINATE"; + } + + if (!mesh.exterior_coordinate || !mesh.exterior_coordinate->space || !mesh.exterior_coordinate->values) { + throw std::runtime_error("A mesh with an external domain must contain a complete exterior-coordinate field before it can be saved."); + } + if (mesh.exterior_coordinate->space->GetMesh() != mesh.mesh.get()) { + throw std::runtime_error("The exterior-coordinate finite-element space is attached to the wrong mesh."); + } + if (mesh.exterior_coordinate->values->FESpace() != mesh.exterior_coordinate->space.get()) { + throw std::runtime_error("The exterior-coordinate grid function is attached to the wrong finite-element space."); + } + + const int scalar_dofs = mesh.exterior_coordinate->space->GetNDofs(); + if (mesh.exterior_coordinate->values->Size() != scalar_dofs) { + throw std::runtime_error("The exterior-coordinate grid function has an invalid size."); + } + + std::stringstream ss; + ss << std::setprecision(std::numeric_limits::max_digits10); + ss << "BEGIN BLOCK EXTERIOR_COORDINATE\n"; + ss << "PRESENT:true\n"; + ss << "NDOFS:" << scalar_dofs << '\n'; + ss << "VALUES:\n"; + + for (int dof = 0; dof < scalar_dofs; ++dof) { + const double coordinate = (*mesh.exterior_coordinate->values)(dof); + if (!std::isfinite(coordinate) || coordinate < 0.0 || coordinate > 1.0) { + throw std::runtime_error(std::format("Exterior-coordinate DOF {} has invalid value {}.", dof, coordinate)); + } + ss << coordinate << '\n'; + } + + ss << "END BLOCK EXTERIOR_COORDINATE"; + return ss.str(); + } } namespace { @@ -350,6 +404,7 @@ END BLOCK CONFIG)", 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 {}; }; }; + auto as_string = [](std::optional* f) { return [f](const std::string_view v) -> std::expected { *f = std::string(v); return {}; }; }; const std::unordered_map handlers = { {"refinement_levels", as_int(&cfg.refinement_levels)}, @@ -369,6 +424,7 @@ END BLOCK CONFIG)", {"vacuum_id", as_size(&cfg.vacuum_id)}, {"optimization_methods-tmop", as_bool(&opt.tmop)}, {"optimization_methods-smoothstep", as_bool(&opt.smoothstep)}, + {"core_mapping", as_string(&cfg.core_mapping)}, }; std::istringstream iss(content); @@ -407,8 +463,122 @@ END BLOCK CONFIG)", StroidMesh mesh; std::string pmesh_raw; std::string rmesh_raw; + std::optional exterior_coordinate_raw; }; + struct ParsedExteriorCoordinate { + bool present{false}; + int scalar_dofs{0}; + std::vector values; + }; + + std::expected parse_exterior_coordinate(const std::string& content) { + ParsedExteriorCoordinate parsed; + std::optional present; + std::optional scalar_dofs; + bool reading_values = false; + std::istringstream iss(content); + std::string line; + + while (std::getline(iss, line)) { + const std::string_view value = trim(line); + if (value.empty() || value.starts_with('#')) continue; + + if (reading_values) { + auto coordinate = parse_double(value); + if (!coordinate) return std::unexpected("EXTERIOR_COORDINATE value -> " + coordinate.error()); + parsed.values.push_back(*coordinate); + continue; + } + + const auto colon = value.find(':'); + if (colon == std::string_view::npos) return std::unexpected(std::format("invalid EXTERIOR_COORDINATE line '{}'.", value)); + + const std::string_view key = trim(value.substr(0, colon)); + const std::string_view field_value = trim(value.substr(colon + 1)); + + if (key == "PRESENT") { + auto result = parse_bool(field_value); + if (!result) return std::unexpected("EXTERIOR_COORDINATE PRESENT -> " + result.error()); + present = *result; + } else if (key == "NDOFS") { + auto result = parse_int(field_value); + if (!result) return std::unexpected("EXTERIOR_COORDINATE NDOFS -> " + result.error()); + scalar_dofs = *result; + } else if (key == "VALUES") { + if (!field_value.empty()) return std::unexpected("EXTERIOR_COORDINATE VALUES must not contain an inline value."); + reading_values = true; + } else { + return std::unexpected(std::format("unknown EXTERIOR_COORDINATE key '{}'.", key)); + } + } + + if (!present.has_value()) return std::unexpected("EXTERIOR_COORDINATE block is missing PRESENT."); + parsed.present = *present; + + if (!parsed.present) { + if (scalar_dofs.has_value() || !parsed.values.empty()) return std::unexpected("An absent exterior coordinate cannot contain NDOFS or VALUES."); + return parsed; + } + + if (!scalar_dofs.has_value() || *scalar_dofs < 0) return std::unexpected("EXTERIOR_COORDINATE block has an invalid or missing NDOFS."); + if (static_cast(parsed.values.size()) != *scalar_dofs) { + return std::unexpected(std::format("EXTERIOR_COORDINATE expected {} values but found {}.", *scalar_dofs, parsed.values.size())); + } + + parsed.scalar_dofs = *scalar_dofs; + return parsed; + } + + std::expected restore_exterior_coordinate(StroidMesh& mesh, const std::optional& raw) { + fourdst::config::Config config; + config.mutate([&mesh](config::MeshConfig& value) { value = mesh.config; }); + + try { + mesh.exterior_coordinate = topology::BuildExteriorCoordinate(*mesh.mesh, *mesh.reference_mesh, config); + } catch (const std::exception& exception) { + return std::unexpected(std::string("failed to reconstruct exterior coordinate: ") + exception.what()); + } + + if (!raw.has_value()) return {}; + + auto parsed = parse_exterior_coordinate(*raw); + if (!parsed) return std::unexpected(parsed.error()); + + const bool include_external_domain = mesh.config.include_external_domain.value_or(true); + if (!parsed->present) { + if (include_external_domain) return std::unexpected("EXTERIOR_COORDINATE is absent even though the mesh includes an external domain."); + if (mesh.exterior_coordinate) return std::unexpected("An exterior-coordinate field was reconstructed for a mesh without an external domain."); + return {}; + } + + if (!include_external_domain) return std::unexpected("EXTERIOR_COORDINATE is present for a mesh without an external domain."); + if (!mesh.exterior_coordinate || !mesh.exterior_coordinate->space || !mesh.exterior_coordinate->values) { + return std::unexpected("Unable to allocate the exterior-coordinate field while loading the mesh."); + } + if (parsed->scalar_dofs != mesh.exterior_coordinate->space->GetNDofs()) { + return std::unexpected(std::format("EXTERIOR_COORDINATE contains {} DOFs but the reconstructed space has {}.", parsed->scalar_dofs, mesh.exterior_coordinate->space->GetNDofs())); + } + + constexpr double consistency_tolerance = 1.0e-12; + + for (int dof = 0; dof < parsed->scalar_dofs; ++dof) { + const double stored_coordinate = parsed->values[static_cast(dof)]; + const double reconstructed_coordinate = (*mesh.exterior_coordinate->values)(dof); + + if (!std::isfinite(stored_coordinate) || stored_coordinate < 0.0 || stored_coordinate > 1.0) { + return std::unexpected(std::format("EXTERIOR_COORDINATE DOF {} has invalid stored value {}.", dof, stored_coordinate)); + } + if (std::abs(stored_coordinate - reconstructed_coordinate) > consistency_tolerance) { + return std::unexpected(std::format("EXTERIOR_COORDINATE DOF {} is inconsistent with the reference mesh: stored value {}, reconstructed value {}.", dof, stored_coordinate, reconstructed_coordinate)); + } + + (*mesh.exterior_coordinate->values)(dof) = stored_coordinate; + } + + return {}; + } + std::expected parse_metadata(std::istream& is) { auto blocks = extract_blocks(is); if (!blocks) return std::unexpected(blocks.error()); @@ -440,6 +610,10 @@ END BLOCK CONFIG)", if (!rmesh) return std::unexpected(rmesh.error()); pm.rmesh_raw = *rmesh; + if (const auto exterior_coordinate = blocks->find("EXTERIOR_COORDINATE"); exterior_coordinate != blocks->end()) { + pm.exterior_coordinate_raw = exterior_coordinate->second; + } + return pm; } @@ -455,16 +629,18 @@ END BLOCK CONFIG)", std::string rmesh = format_reference_mesh(mesh); std::string config = format_config(mesh); + std::string exterior_coordinate = format_exterior_coordinate(mesh); ofs << header << "\n"; ofs << pmesh << "\n"; ofs << rmesh << "\n"; ofs << config << "\n"; + ofs << exterior_coordinate << "\n"; } void SaveMesh(const mfem::Mesh& mesh, const std::string& filename) { std::ofstream ofs(filename); - ofs.precision(8); + ofs.precision(std::numeric_limits::max_digits10); mesh.Print(ofs); } @@ -574,6 +750,7 @@ END BLOCK CONFIG)", pm->mesh.mesh = std::move(*m); pm->mesh.reference_mesh = std::move(*rm); + if (auto result = restore_exterior_coordinate(pm->mesh, pm->exterior_coordinate_raw); !result) return std::unexpected(result.error()); return std::move(pm->mesh); } @@ -609,6 +786,7 @@ END BLOCK CONFIG)", pm->mesh.mesh = std::move(*m); pm->mesh.reference_mesh = std::move(*rm); + if (auto result = restore_exterior_coordinate(pm->mesh, pm->exterior_coordinate_raw); !result) return std::unexpected(result.error()); return std::move(pm->mesh); } @@ -621,4 +799,4 @@ END BLOCK CONFIG)", #endif // MFEM_USE_MPI -} \ No newline at end of file +} diff --git a/src/lib/refinement/uniform.cpp b/src/lib/refinement/uniform.cpp index 2c51640..26b6d45 100644 --- a/src/lib/refinement/uniform.cpp +++ b/src/lib/refinement/uniform.cpp @@ -4,6 +4,7 @@ #include "stroid/utils/types.h" #include "stroid/utils/mesh_utils.h" #include "stroid/exceptions/exceptions.h" +#include "stroid/topology/curvilinear.h" #include "stroid/topology/topology.h" #include "stroid/topology/optimize.h" @@ -17,6 +18,10 @@ namespace stroid::refinement { return; } + if (!mesh.mesh) { + throw exceptions::StroidMissingReferenceMesh("UniformRefinement requires a primary mesh to be present in the StroidMesh object. This should be present by construction and the fact that it is missing represents a bug. Please report this to the stroid developers on GitHub or by email at emily.boudreaux@dartmouth.edu"); + } + mesh.exterior_coordinate.reset(); for (size_t i = 0; i < levels; i++) { mesh.reference_mesh->UniformRefinement(); } @@ -31,9 +36,7 @@ namespace stroid::refinement { cfg.mutate(Mutator); mesh.mesh = utils::BuildProjected(*mesh.reference_mesh, cfg); - - - topology::OptimizeMesh(*mesh.mesh, cfg); + mesh.exterior_coordinate = topology::BuildExteriorCoordinate(*mesh.mesh, *mesh.reference_mesh, cfg); } -} +} \ No newline at end of file diff --git a/src/lib/topology/curvilinear.cpp b/src/lib/topology/curvilinear.cpp index 1c8a47d..65c55c2 100644 --- a/src/lib/topology/curvilinear.cpp +++ b/src/lib/topology/curvilinear.cpp @@ -3,6 +3,41 @@ #include +namespace { + double compute_exterior_coordinate( + const mfem::Vector& logical_position, + const int attribute, + const fourdst::config::Config& config + ) { + if (!config->include_external_domain.value_or(true) || attribute != static_cast(config->vacuum_id.value_or(3))) return 0.0; + + const double r_star = config->r_star.value_or(1.0); + const double r_infinity = config->r_infinity.value_or(6.0); + const double radial_extent = r_infinity - r_star; + + if (!std::isfinite(r_star) || !std::isfinite(r_infinity) || r_star <= 0.0 || radial_extent <= 0.0) { + throw std::invalid_argument("Exterior-coordinate construction requires finite radii with 0 < r_star < r_infinity."); + } + + double logical_radius = 0.0; + for (int d = 0; d < logical_position.Size(); ++d) { + if (!std::isfinite(logical_position(d))) throw std::runtime_error("Reference mesh produced a non-finite logical position."); + logical_radius = std::max(logical_radius, std::abs(logical_position(d))); + } + + double coordinate = (logical_radius - r_star) / radial_extent; + const double tolerance = 1024.0 * std::numeric_limits::epsilon() * std::max({1.0, std::abs(r_star), std::abs(r_infinity)}) / radial_extent; + + if (coordinate < -tolerance || coordinate > 1.0 + tolerance) { + throw std::runtime_error(std::format("Logical exterior coordinate {} lies outside [0, 1].", coordinate)); + } + + if (std::abs(coordinate) <= tolerance) coordinate = 0.0; + if (std::abs(coordinate - 1.0) <= tolerance) coordinate = 1.0; + return coordinate; + } +} + namespace stroid::topology { void PromoteToHighOrder(mfem::Mesh &mesh, const fourdst::config::Config &config) { const auto* fec = new mfem::H1_FECollection(config->order.value(), mesh.Dimension()); @@ -54,4 +89,89 @@ namespace stroid::topology { } } + + std::unique_ptr BuildExteriorCoordinate( + mfem::Mesh& mesh, + mfem::Mesh& reference_mesh, + const fourdst::config::Config& config + ) { + if (!config->include_external_domain.value_or(true)) return nullptr; + if (mesh.Dimension() != reference_mesh.Dimension() || mesh.SpaceDimension() != reference_mesh.SpaceDimension()) { + throw std::invalid_argument("Primary and reference meshes must have matching dimensions when constructing the exterior coordinate."); + } + if (mesh.GetNE() != reference_mesh.GetNE()) { + throw std::invalid_argument("Primary and reference meshes must have the same number of elements when constructing the exterior coordinate."); + } + if (mesh.GetNodalFESpace() == nullptr) { + throw std::invalid_argument("Exterior-coordinate construction requires a primary mesh with a nodal finite-element space."); + } + + for (int element_id = 0; element_id < mesh.GetNE(); ++element_id) { + if (mesh.GetElementGeometry(element_id) != reference_mesh.GetElementGeometry(element_id)) { + throw std::invalid_argument(std::format("Primary and reference element {} have different geometries.", element_id)); + } + if (mesh.GetAttribute(element_id) != reference_mesh.GetAttribute(element_id)) { + throw std::invalid_argument(std::format("Primary and reference element {} have different attributes.", element_id)); + } + } + + auto field = std::make_unique(); + const mfem::FiniteElementCollection* collection = mesh.GetNodalFESpace()->FEColl(); + field->space = std::make_unique(&mesh, collection); + field->values = std::make_unique(field->space.get()); + *field->values = 0.0; + + const int scalar_dofs = field->space->GetNDofs(); + std::vector processed(static_cast(scalar_dofs), false); + mfem::Array element_dofs; + mfem::Vector logical_position(reference_mesh.SpaceDimension()); + + const double consistency_tolerance = 4096.0 * std::numeric_limits::epsilon(); + + for (int element_id = 0; element_id < mesh.GetNE(); ++element_id) { + const mfem::FiniteElement& element = *field->space->GetFE(element_id); + const mfem::IntegrationRule& nodes = element.GetNodes(); + mfem::ElementTransformation* reference_transformation = reference_mesh.GetElementTransformation(element_id); + + if (reference_transformation == nullptr) throw std::runtime_error(std::format("Reference element {} has no element transformation.", element_id)); + + field->space->GetElementDofs(element_id, element_dofs); + if (nodes.GetNPoints() != element_dofs.Size()) { + throw std::runtime_error(std::format("Element {} has {} nodal points but {} scalar DOFs.", element_id, nodes.GetNPoints(), element_dofs.Size())); + } + + for (int local_dof = 0; local_dof < element_dofs.Size(); ++local_dof) { + const int encoded_dof = element_dofs[local_dof]; + const int global_dof = encoded_dof >= 0 ? encoded_dof : -1 - encoded_dof; + + if (global_dof < 0 || global_dof >= scalar_dofs) { + throw std::runtime_error(std::format("Element {} references invalid scalar DOF {}.", element_id, global_dof)); + } + + reference_transformation->Transform(nodes.IntPoint(local_dof), logical_position); + const double coordinate = compute_exterior_coordinate(logical_position, mesh.GetAttribute(element_id), config); + + if (processed[static_cast(global_dof)]) { + const double existing_coordinate = (*field->values)(global_dof); + if (std::abs(existing_coordinate - coordinate) > consistency_tolerance) { + throw std::runtime_error(std::format("Exterior coordinate is inconsistent at shared scalar DOF {}: existing value {}, new value {} from element {}.", global_dof, existing_coordinate, coordinate, element_id)); + } + continue; + } + + (*field->values)(global_dof) = coordinate; + processed[static_cast(global_dof)] = true; + } + } + + for (int dof = 0; dof < scalar_dofs; ++dof) { + if (!processed[static_cast(dof)]) throw std::runtime_error(std::format("Exterior-coordinate scalar DOF {} was not assigned.", dof)); + const double coordinate = (*field->values)(dof); + if (!std::isfinite(coordinate) || coordinate < 0.0 || coordinate > 1.0) { + throw std::runtime_error(std::format("Exterior-coordinate scalar DOF {} has invalid value {}.", dof, coordinate)); + } + } + + return field; + } } diff --git a/src/lib/topology/mapping.cpp b/src/lib/topology/mapping.cpp index 87cda49..0b683cc 100644 --- a/src/lib/topology/mapping.cpp +++ b/src/lib/topology/mapping.cpp @@ -96,7 +96,14 @@ namespace stroid::topology { double Z = pos(2); double maxAbs = std::max({std::abs(X), std::abs(Y), std::abs(Z)}); - if (maxAbs < 1e-14) return; + const bool multi_block = config->core_mapping.value_or("spherified") == "multi_block"; + const double inner_radius = config->r_core.value() / 2.0; + if (multi_block && maxAbs <= inner_radius) { + pos /= std::sqrt(3.0); + ApplySpheroidal(pos, config); + return; + } + if (!multi_block && maxAbs < 1e-14) return; double cx = X / maxAbs; double cy = Y / maxAbs; @@ -112,6 +119,15 @@ namespace stroid::topology { unit_dir(2) = sz; if (maxAbs <= config->r_core.value()) { + if (multi_block) { + const double t = (maxAbs - inner_radius) / inner_radius; + const double inner_scale = inner_radius / std::sqrt(3.0); + pos(0) = (1.0 - t) * inner_scale * cx + t * config->r_core.value() * sx; + pos(1) = (1.0 - t) * inner_scale * cy + t * config->r_core.value() * sy; + pos(2) = (1.0 - t) * inner_scale * cz + t * config->r_core.value() * sz; + ApplySpheroidal(pos, config); + return; + } double nx = X / config->r_core.value(); double ny = Y / config->r_core.value(); double nz = Z / config->r_core.value(); @@ -142,68 +158,31 @@ 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))); - // } - // - // 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); - // } - // } + double ComputeExteriorCoordinate( + const mfem::Vector& logical_position, + const int attribute, + const fourdst::config::Config& config + ) { + if (!config->include_external_domain.value() || attribute != static_cast(config->vacuum_id.value())) return 0.0; + + const double logical_radius = std::max({ + std::abs(logical_position(0)), + std::abs(logical_position(1)), + std::abs(logical_position(2)) + }); + + const double r_star = config->r_star.value(); + const double r_infinity = config->r_infinity.value(); + const double coordinate = (logical_radius - r_star) / (r_infinity - r_star); + constexpr double tolerance = 64.0 * std::numeric_limits::epsilon(); + + if (coordinate < -tolerance || coordinate > 1.0 + tolerance) { + throw std::runtime_error("Logical exterior coordinate lies outside [0, 1]."); + } + + if (std::abs(coordinate) <= tolerance) return 0.0; + if (std::abs(coordinate - 1.0) <= tolerance) return 1.0; + return coordinate; + } + } diff --git a/src/lib/topology/topology.cpp b/src/lib/topology/topology.cpp index 0672899..63cbaa8 100644 --- a/src/lib/topology/topology.cpp +++ b/src/lib/topology/topology.cpp @@ -1,6 +1,8 @@ #include "mfem.hpp" #include #include +#include +#include #include "stroid/config/config.h" #include "fourdst/config/config.h" @@ -8,9 +10,31 @@ 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 = config->include_external_domain ? 12 : 6; + const std::string core_mapping = config->core_mapping.value_or("spherified"); + if (core_mapping != "spherified" && core_mapping != "multi_block") { + throw std::invalid_argument("Unknown core_mapping: " + core_mapping); + } + + const bool multi_block = core_mapping == "multi_block"; + const bool include_external_domain = config->include_external_domain.value_or(true); + if (multi_block) { + const double r_core = config->r_core.value(); + const double r_star = config->r_star.value(); + const double r_infinity = config->r_infinity.value_or(6.0); + const double flattening = config->flattening.value(); + if (!std::isfinite(r_core) || !std::isfinite(r_star) || r_core <= 0.0 || r_star <= r_core || + (include_external_domain && (!std::isfinite(r_infinity) || r_infinity <= r_star))) { + throw std::invalid_argument("multi_block requires 0 < r_core < r_star < r_infinity (when external)."); + } + if (!std::isfinite(flattening) || flattening >= 1.0) { + throw std::invalid_argument("multi_block requires finite flattening < 1."); + } + } + + const int offset = multi_block ? 8 : 0; + int nVert = (include_external_domain ? 24 : 16) + offset; + int nElem = (include_external_domain ? 13 : 7) + (multi_block ? 6 : 0); + int nBev = include_external_domain ? 12 : 6; auto mesh = std::make_unique(3, nVert, nElem, nBev, 3); @@ -21,9 +45,12 @@ namespace stroid::topology { mesh->AddVertex(x, y, z); }; + if (multi_block) { + add_box(config->r_core.value() / 2.0); + } add_box(config->r_core.value()); add_box(config->r_star.value()); - if (config->include_external_domain) { + if (include_external_domain) { add_box(config->r_infinity.value()); } @@ -38,11 +65,18 @@ namespace stroid::topology { {1, 3, 7, 5, 9, 11, 15, 13}, // +X face {0, 4, 6, 2, 8, 12, 14, 10} // -X face }; + if (multi_block) { + for (const auto & shell : stellar_shells) { + mesh->AddHex(shell.data(), config->core_id.value()); + } + } for (const auto & shell : stellar_shells) { - mesh->AddHex(shell.data(), config->envelope_id.value()); + auto vertices = shell; + for (auto & vertex : vertices) vertex += offset; + mesh->AddHex(vertices.data(), config->envelope_id.value()); } - if (config->include_external_domain) { + if (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}); @@ -51,7 +85,9 @@ 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.value()); + auto vertices = shell; + for (auto & vertex : vertices) vertex += offset; + mesh->AddHex(vertices.data(), config->vacuum_id.value()); } } @@ -66,10 +102,12 @@ namespace stroid::topology { }; for (const auto& bdr: surface_bdr_quads) { - mesh->AddBdrQuad(bdr, config->surface_bdr_id.value()); + int vertices[4]; + for (int i = 0; i < 4; ++i) vertices[i] = bdr[i] + offset; + mesh->AddBdrQuad(vertices, config->surface_bdr_id.value()); } - if (config->include_external_domain) { + if (include_external_domain) { const int inf_bdr_quads[6][4] = { {16, 17, 21, 20}, {17, 19, 23, 21}, @@ -80,13 +118,16 @@ namespace stroid::topology { }; for (const auto& bdr: inf_bdr_quads) { - mesh->AddBdrQuad(bdr, config->inf_bdr_id.value()); + int vertices[4]; + for (int i = 0; i < 4; ++i) vertices[i] = bdr[i] + offset; + mesh->AddBdrQuad(vertices, config->inf_bdr_id.value()); } } return mesh; } + // ReSharper disable once CppUseInternalLinkage void Finalize(mfem::Mesh& mesh, const fourdst::config::Config &config) { mesh.FinalizeTopology(); mesh.Finalize(); diff --git a/src/meson.build b/src/meson.build index 99d9587..3b93bae 100644 --- a/src/meson.build +++ b/src/meson.build @@ -41,7 +41,7 @@ if get_option('build_python') endif libstroid = static_library( - 'libstroid', + 'stroid', stroid_sources, include_directories: stroid_include_files, dependencies: dependencies, @@ -52,11 +52,11 @@ if get_option('build_python') ) else libstroid = static_library( - 'libstroid', + 'stroid', stroid_sources, include_directories: stroid_include_files, dependencies: dependencies, - install: true, + install: true ) endif @@ -81,7 +81,9 @@ endif meson.override_dependency('stroid', stroid_dep) +message('stroid include dir: ' + stroid_includedir) install_subdir( 'include/stroid', - install_dir: get_option('includedir') / 'stroid' + install_dir: stroid_includedir, + exclude_files: ['version.h.in'] ) \ No newline at end of file diff --git a/src/python/config/bindings.cpp b/src/python/config/bindings.cpp index ee0ca8f..8b1ba23 100644 --- a/src/python/config/bindings.cpp +++ b/src/python/config/bindings.cpp @@ -52,7 +52,8 @@ void register_config_bindings(pybind11::module_& m) { .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 + .optimization_methods = kwargs.contains("optimization_methods") ? kwargs["optimization_methods"].cast() : opt_method, + .core_mapping = kwargs.contains("core_mapping") ? kwargs["core_mapping"].cast() : "spherified" }; })) .def_property( @@ -198,5 +199,20 @@ void register_config_bindings(pybind11::module_& m) { [](stroid::config::MeshConfig& self, stroid::config::OptimizationMethods value) { self.optimization_methods = value; } - ); -} \ No newline at end of file + ) + .def_property( + "core_mapping", + [](const stroid::config::MeshConfig& self) { + return self.core_mapping; + }, + [](stroid::config::MeshConfig& self, const std::string& value) { + if (value != "spherified" && value != "multi_block") { + throw std::invalid_argument("Invalid core_mapping value. Must be 'spherified' or 'multi_block'."); + } + self.core_mapping = value; + } + ) + .def("__repr__", [](const stroid::config::MeshConfig& self) { + return stroid::config::to_string(self); + }); +} diff --git a/tests/python/core_mapping_smoke.py b/tests/python/core_mapping_smoke.py new file mode 100644 index 0000000..f9d5022 --- /dev/null +++ b/tests/python/core_mapping_smoke.py @@ -0,0 +1,126 @@ +"""Focused core-mapping smoke checks for a built or staged STROID Python module. + +Example: + python core_mapping_smoke.py --module _stroid --module-dir build/build-python + python core_mapping_smoke.py --module-dir /path/to/staged/site-packages +""" + +import argparse +import importlib +import json +from pathlib import Path +import sys +import tempfile + + +def element_counts(stroid, mesh): + result = stroid.stats.ComputeMeshStats( + mesh, stroid.stats.MeshStatFeatures.ELEMENT_COUNT + ) + assert not result.errors, result.errors + assert result.element_counts is not None + return result.element_counts + + +def run(stroid): + default = stroid.config.MeshConfig() + assert default.core_mapping == "spherified" + default.core_mapping = "multi_block" + assert default.core_mapping == "multi_block" + assert "core_mapping: multi_block" in repr(default) + + summaries = [] + with tempfile.TemporaryDirectory(prefix="stroid-python-smoke-") as output: + output_path = Path(output) + config_path = output_path / "multi_block.toml" + config_path.write_text( + '[main]\ncore_mapping = "multi_block"\nrefinement_levels = 0\n' + 'order = 3\ninclude_external_domain = false\n' + 'r_core = 0.25\nr_star = 1.0\nr_infinity = 6.0\n' + 'flattening = 0.0\nr_instability = 1e-14\n' + 'core_steepness = 1.0\ncontinuity_order = 2\n' + 'surface_bdr_id = 1\ninf_bdr_id = 2\n' + 'core_id = 1\nenvelope_id = 2\nvacuum_id = 3\n' + '[main.optimization_methods]\ntmop = false\nsmoothstep = true\n' + ) + configured_mesh = stroid.GenerateMesh(str(config_path)) + assert configured_mesh.config.core_mapping == "multi_block" + assert element_counts(stroid, configured_mesh).total == 13 + + for mapping in ("spherified", "multi_block"): + for external in (False, True): + config = stroid.config.MeshConfig( + core_mapping=mapping, + refinement_levels=0, + order=3, + include_external_domain=external, + optimization_methods=stroid.config.OptimizationMethods( + tmop=False, smoothstep=True + ), + ) + mesh = stroid.GenerateMesh(config) + assert mesh.has_mesh() and mesh.has_rmesh() + counts = element_counts(stroid, mesh) + expected_core = 7 if mapping == "multi_block" else 1 + expected_total = expected_core + 6 + (6 if external else 0) + assert counts.total == expected_total + assert counts.core == expected_core + assert counts.envelope == 6 + assert counts.vacuum == (6 if external else 0) + + path = output_path / f"{mapping}-{external}.smesh" + stroid.IO.SaveStroidMesh(mesh, str(path), "Python core-mapping smoke") + loaded = stroid.IO.LoadStroidMesh(str(path)) + assert loaded.config.core_mapping == mapping + assert element_counts(stroid, loaded).total == expected_total + + stroid.refinement.UniformRefinement(loaded, 1) + assert loaded.config.core_mapping == mapping + assert loaded.refinement_levels == 1 + refined_counts = element_counts(stroid, loaded) + assert refined_counts.total == expected_total * 8 + assert refined_counts.core == expected_core * 8 + + if mapping == "spherified": + legacy = "\n".join( + line for line in path.read_text().splitlines() + if not line.startswith("core_mapping:") + ) + legacy_mesh = stroid.IO.ParseStroidMesh(legacy) + assert legacy_mesh.config.core_mapping == "spherified" + assert element_counts(stroid, legacy_mesh).total == expected_total + + summaries.append({ + "mapping": mapping, + "external": external, + "initial_elements": expected_total, + "refined_elements": refined_counts.total, + }) + + invalid = stroid.config.MeshConfig( + core_mapping="unknown", refinement_levels=0 + ) + try: + stroid.GenerateMesh(invalid) + except (ValueError, RuntimeError): + pass + else: + raise AssertionError("Unsupported core_mapping was accepted") + + return summaries + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--module", default="stroid") + parser.add_argument("--module-dir", type=Path) + args = parser.parse_args() + if args.module_dir is not None: + sys.path.insert(0, str(args.module_dir.resolve())) + stroid = importlib.import_module(args.module) + summaries = run(stroid) + print(json.dumps({"module": stroid.__file__, "cases": summaries}, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/tests/stroidTest.cpp b/tests/stroidTest.cpp index b8000d1..c30b265 100644 --- a/tests/stroidTest.cpp +++ b/tests/stroidTest.cpp @@ -7,10 +7,13 @@ #include "stroid/topology/mapping.h" #include "stroid/topology/topology.h" #include "stroid/utils/mesh_utils.h" +#include "stroid/stroid.h" #include #include +#include #include +#include #include #include #include @@ -126,6 +129,94 @@ std::unique_ptr BuildProjectedMesh(const Config& cfg) { return mesh; } +void ExpectExteriorCoordinateRange(stroid::StroidMesh& stroid_mesh) { + ASSERT_NE(stroid_mesh.mesh, nullptr); + ASSERT_NE(stroid_mesh.exterior_coordinate, nullptr); + ASSERT_NE(stroid_mesh.exterior_coordinate->space, nullptr); + ASSERT_NE(stroid_mesh.exterior_coordinate->values, nullptr); + ASSERT_EQ(stroid_mesh.exterior_coordinate->space->GetMesh(), stroid_mesh.mesh.get()); + ASSERT_EQ(stroid_mesh.exterior_coordinate->values->FESpace(), stroid_mesh.exterior_coordinate->space.get()); + + mfem::Mesh& mesh = *stroid_mesh.mesh; + mfem::GridFunction& coordinate = *stroid_mesh.exterior_coordinate->values; + const int vacuum_attribute = static_cast(stroid_mesh.config.vacuum_id.value()); + bool sampled_vacuum = false; + + for (int element_id = 0; element_id < mesh.GetNE(); ++element_id) { + const mfem::FiniteElement& element = *stroid_mesh.exterior_coordinate->space->GetFE(element_id); + const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(element.GetGeomType(), 2 * element.GetOrder() + 4); + + for (int q = 0; q < integration_rule.GetNPoints(); ++q) { + const double value = coordinate.GetValue(element_id, integration_rule.IntPoint(q)); + EXPECT_TRUE(std::isfinite(value)); + + if (mesh.GetAttribute(element_id) == vacuum_attribute) { + sampled_vacuum = true; + EXPECT_GE(value, -1.0e-12); + EXPECT_LE(value, 1.0 + 1.0e-12); + } else { + EXPECT_NEAR(value, 0.0, 1.0e-12); + } + } + } + + EXPECT_TRUE(sampled_vacuum); +} + +void ExpectExteriorCoordinateBoundaryTraces(stroid::StroidMesh& stroid_mesh) { + ASSERT_NE(stroid_mesh.mesh, nullptr); + ASSERT_NE(stroid_mesh.exterior_coordinate, nullptr); + ASSERT_NE(stroid_mesh.exterior_coordinate->values, nullptr); + + mfem::Mesh& mesh = *stroid_mesh.mesh; + mfem::GridFunction& coordinate = *stroid_mesh.exterior_coordinate->values; + const int vacuum_attribute = static_cast(stroid_mesh.config.vacuum_id.value()); + const int infinity_boundary = static_cast(stroid_mesh.config.inf_bdr_id.value()); + int stellar_vacuum_faces = 0; + int infinity_faces = 0; + + for (int face_id = 0; face_id < mesh.GetNumFaces(); ++face_id) { + mfem::FaceElementTransformations* transformation = mesh.GetFaceElementTransformations(face_id); + if (transformation == nullptr || transformation->Elem1 == nullptr || transformation->Elem2 == nullptr) continue; + + const bool element_1_vacuum = transformation->Elem1->Attribute == vacuum_attribute; + const bool element_2_vacuum = transformation->Elem2->Attribute == vacuum_attribute; + if (element_1_vacuum == element_2_vacuum) continue; + + ++stellar_vacuum_faces; + const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 6); + + for (int q = 0; q < integration_rule.GetNPoints(); ++q) { + const mfem::IntegrationPoint& face_point = integration_rule.IntPoint(q); + transformation->SetAllIntPoints(&face_point); + + const int vacuum_element = element_1_vacuum ? transformation->Elem1No : transformation->Elem2No; + const mfem::IntegrationPoint& vacuum_point = element_1_vacuum ? transformation->Elem1->GetIntPoint() : transformation->Elem2->GetIntPoint(); + EXPECT_NEAR(coordinate.GetValue(vacuum_element, vacuum_point), 0.0, 1.0e-12); + } + } + + for (int boundary_element = 0; boundary_element < mesh.GetNBE(); ++boundary_element) { + if (mesh.GetBdrAttribute(boundary_element) != infinity_boundary) continue; + + mfem::FaceElementTransformations* transformation = mesh.GetBdrFaceTransformations(boundary_element); + ASSERT_NE(transformation, nullptr); + ASSERT_NE(transformation->Elem1, nullptr); + ++infinity_faces; + + const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 6); + + for (int q = 0; q < integration_rule.GetNPoints(); ++q) { + const mfem::IntegrationPoint& face_point = integration_rule.IntPoint(q); + transformation->SetAllIntPoints(&face_point); + EXPECT_NEAR(coordinate.GetValue(transformation->Elem1No, transformation->Elem1->GetIntPoint()), 1.0, 1.0e-12); + } + } + + EXPECT_GT(stellar_vacuum_faces, 0); + EXPECT_GT(infinity_faces, 0); +} + 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.value(), cfg->envelope_id.value(), cfg->vacuum_id.value()})); @@ -1043,6 +1134,134 @@ TEST_F(stroidTest, Refinement_UniformRefinementProducesExpectedElementCounts) { EXPECT_EQ(mesh.mesh->GetNE(), init_elements * 8); } +TEST_F(stroidTest, ExteriorCoordinate_HasValidRangeAndExactBoundaryTraces) { + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml"); + const auto& cfg = *cfg_ptr; + + stroid::StroidMesh mesh; + ASSERT_NO_THROW(mesh = stroid::GenerateMesh(cfg)); + + ExpectExteriorCoordinateRange(mesh); + ExpectExteriorCoordinateBoundaryTraces(mesh); +} + +TEST_F(stroidTest, ExteriorCoordinate_IsRebuiltAfterUniformRefinement) { + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml"); + const auto& cfg = *cfg_ptr; + + stroid::StroidMesh mesh; + ASSERT_NO_THROW(mesh = stroid::GenerateMesh(cfg)); + ASSERT_NE(mesh.exterior_coordinate, nullptr); + + const int initial_elements = mesh.mesh->GetNE(); + const int initial_coordinate_dofs = mesh.exterior_coordinate->space->GetNDofs(); + + ASSERT_NO_THROW(stroid::refinement::UniformRefinement(mesh, 1)); + ASSERT_NE(mesh.exterior_coordinate, nullptr); + EXPECT_EQ(mesh.mesh->GetNE(), initial_elements * 8); + EXPECT_GT(mesh.exterior_coordinate->space->GetNDofs(), initial_coordinate_dofs); + + ExpectExteriorCoordinateRange(mesh); + ExpectExteriorCoordinateBoundaryTraces(mesh); +} + +TEST_F(stroidTest, ExteriorCoordinate_SurvivesSaveAndLoad) { + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml"); + const auto& cfg = *cfg_ptr; + + stroid::StroidMesh original; + ASSERT_NO_THROW(original = stroid::GenerateMesh(cfg)); + ASSERT_NE(original.exterior_coordinate, nullptr); + + const std::filesystem::path path = std::filesystem::temp_directory_path() / "stroid_exterior_coordinate_round_trip.smesh"; + ASSERT_NO_THROW(stroid::IO::SaveStroidMesh(original, path.string(), "Exterior-coordinate round-trip test")); + + auto loaded_result = stroid::IO::LoadStroidMesh(path.string()); + if (!loaded_result.has_value()) FAIL() << loaded_result.error(); + stroid::StroidMesh loaded = std::move(*loaded_result); + + ASSERT_NE(loaded.exterior_coordinate, nullptr); + ASSERT_EQ(loaded.exterior_coordinate->space->GetNDofs(), original.exterior_coordinate->space->GetNDofs()); + ASSERT_EQ(loaded.exterior_coordinate->values->Size(), original.exterior_coordinate->values->Size()); + + for (int dof = 0; dof < original.exterior_coordinate->values->Size(); ++dof) { + EXPECT_DOUBLE_EQ((*loaded.exterior_coordinate->values)(dof), (*original.exterior_coordinate->values)(dof)); + } + + ExpectExteriorCoordinateRange(loaded); + ExpectExteriorCoordinateBoundaryTraces(loaded); + + std::error_code error; + std::filesystem::remove(path, error); + EXPECT_FALSE(error); +} + +TEST_F(stroidTest, ExteriorCoordinate_IsAbsentWithoutExternalDomainAcrossSaveAndLoad) { + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml"); + const auto& cfg = *cfg_ptr; + + stroid::StroidMesh original; + ASSERT_NO_THROW(original = stroid::GenerateMesh(cfg)); + EXPECT_EQ(original.exterior_coordinate, nullptr); + + const std::filesystem::path path = std::filesystem::temp_directory_path() / "stroid_no_exterior_coordinate_round_trip.smesh"; + ASSERT_NO_THROW(stroid::IO::SaveStroidMesh(original, path.string(), "No-exterior-coordinate round-trip test")); + + auto loaded_result = stroid::IO::LoadStroidMesh(path.string()); + if (!loaded_result.has_value()) FAIL() << loaded_result.error(); + EXPECT_EQ(loaded_result->exterior_coordinate, nullptr); + + std::error_code error; + std::filesystem::remove(path, error); + EXPECT_FALSE(error); +} + +TEST_F(stroidTest, ExteriorCoordinate_IsReconstructedWhenLoadingLegacyFiles) { + const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml"); + const auto& cfg = *cfg_ptr; + + stroid::StroidMesh original; + ASSERT_NO_THROW(original = stroid::GenerateMesh(cfg)); + ASSERT_NE(original.exterior_coordinate, nullptr); + + const std::filesystem::path path = std::filesystem::temp_directory_path() / "stroid_legacy_exterior_coordinate.smesh"; + ASSERT_NO_THROW(stroid::IO::SaveStroidMesh(original, path.string(), "Legacy exterior-coordinate reconstruction test")); + + std::ifstream input(path); + ASSERT_TRUE(input.is_open()); + std::string contents((std::istreambuf_iterator(input)), std::istreambuf_iterator()); + + constexpr std::string_view begin_marker = "BEGIN BLOCK EXTERIOR_COORDINATE"; + constexpr std::string_view end_marker = "END BLOCK EXTERIOR_COORDINATE"; + const size_t begin = contents.find(begin_marker); + const size_t end_begin = contents.find(end_marker); + ASSERT_NE(begin, std::string::npos); + ASSERT_NE(end_begin, std::string::npos); + + size_t end = end_begin + end_marker.size(); + if (end < contents.size() && contents[end] == '\n') ++end; + contents.erase(begin, end - begin); + + std::istringstream legacy_stream(contents); + auto loaded_result = stroid::IO::ParseStroidMesh(legacy_stream); + if (!loaded_result.has_value()) FAIL() << loaded_result.error(); + + stroid::StroidMesh loaded = std::move(*loaded_result); + ASSERT_NE(loaded.exterior_coordinate, nullptr); + ASSERT_EQ(loaded.exterior_coordinate->values->Size(), original.exterior_coordinate->values->Size()); + + for (int dof = 0; dof < original.exterior_coordinate->values->Size(); ++dof) { + EXPECT_NEAR((*loaded.exterior_coordinate->values)(dof), (*original.exterior_coordinate->values)(dof), 1.0e-12); + } + + ExpectExteriorCoordinateRange(loaded); + ExpectExteriorCoordinateBoundaryTraces(loaded); + + std::error_code error; + std::filesystem::remove(path, error); + EXPECT_FALSE(error); +} + TEST_F(stroidTest, Stats_ComputeStats) { const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml"); const auto& cfg = *cfg_ptr; @@ -1055,4 +1274,360 @@ TEST_F(stroidTest, Stats_ComputeStats) { } +namespace { +std::unique_ptr MultiBlockConfiguration(int order, int refinement, bool external, double flattening = 0.0) { + auto cfg = std::make_unique(); + cfg->mutate([&](stroid::config::MeshConfig& value) { + value.core_mapping = "multi_block"; + value.order = order; + value.refinement_levels = refinement; + value.include_external_domain = external; + value.flattening = flattening; + value.optimization_methods = stroid::config::OptimizationMethods{false, false}; + }); + return cfg; +} + +// Unlike CollectConditioningStats, this uses the signed determinant, actual +// singular values, and a closed sample grid including vertices/edges/faces. +// Column-length ratios and open quadrature points miss the old core-corner defect. +void ExpectClosedGridCoreConditioning(mfem::Mesh& mesh, int coreAttribute, double maximumCondition = 10.0) { + int coreElements = 0; + double largestCondition = 0.0; + double smallestDeterminant = std::numeric_limits::infinity(); + for (int element = 0; element < mesh.GetNE(); ++element) { + if (mesh.GetAttribute(element) != coreAttribute) continue; + ++coreElements; + auto* transformation = mesh.GetElementTransformation(element); + ASSERT_EQ(transformation->GetGeometryType(), mfem::Geometry::CUBE); + for (double x : {0.0, 0.01, 0.5, 0.99, 1.0}) { + for (double y : {0.0, 0.01, 0.5, 0.99, 1.0}) { + for (double z : {0.0, 0.01, 0.5, 0.99, 1.0}) { + mfem::IntegrationPoint point; + point.Set3(x, y, z); + transformation->SetIntPoint(&point); + const auto& jacobian = transformation->Jacobian(); + const double determinant = jacobian.Det(); + const double minimumSingular = jacobian.CalcSingularvalue(2); + const double maximumSingular = jacobian.CalcSingularvalue(0); + ASSERT_TRUE(std::isfinite(determinant)); + ASSERT_GT(determinant, 0.0) << "element=" << element << " point=" << x << ',' << y << ',' << z; + ASSERT_TRUE(std::isfinite(minimumSingular)); + ASSERT_GT(minimumSingular, 0.0) << "element=" << element; + const double condition = maximumSingular / minimumSingular; + ASSERT_TRUE(std::isfinite(condition)); + ASSERT_LT(condition, maximumCondition) + << "element=" << element << " point=" << x << ',' << y << ',' << z; + smallestDeterminant = std::min(smallestDeterminant, determinant); + largestCondition = std::max(largestCondition, condition); + } + } + } + } + EXPECT_GT(coreElements, 0); + EXPECT_GT(smallestDeterminant, 0.0); + EXPECT_LT(largestCondition, maximumCondition); +} + +void ExpectCoreFaceContinuity(mfem::Mesh& mesh, int coreAttribute) { + int faces = 0; + mfem::Vector left(3), right(3); + for (int face = 0; face < mesh.GetNumFaces(); ++face) { + auto* transformation = mesh.GetFaceElementTransformations(face); + if (transformation == nullptr || transformation->Elem1 == nullptr || transformation->Elem2 == nullptr) continue; + if (transformation->Elem1->Attribute != coreAttribute && transformation->Elem2->Attribute != coreAttribute) continue; + ++faces; + for (double x : {0.0, 0.25, 0.5, 0.75, 1.0}) { + for (double y : {0.0, 0.25, 0.5, 0.75, 1.0}) { + mfem::IntegrationPoint point; + point.Set2(x, y); + transformation->SetAllIntPoints(&point); + transformation->Elem1->Transform(transformation->Elem1->GetIntPoint(), left); + transformation->Elem2->Transform(transformation->Elem2->GetIntPoint(), right); + left -= right; + EXPECT_LT(left.Norml2(), 2.0e-12) << "face=" << face; + } + } + } + EXPECT_GT(faces, 0); +} + +} // namespace + +TEST_F(stroidTest, MultiBlockCore_TopologyCountsAndAttributesAreOptIn) { + EXPECT_EQ(stroid::config::MeshConfig{}.core_mapping.value(), "spherified"); + for (const bool external : {false, true}) { + SCOPED_TRACE(external); + auto cfg = MultiBlockConfiguration(2, 0, external); + auto mesh = stroid::topology::BuildSkeleton(*cfg); + ASSERT_NE(mesh, nullptr); + EXPECT_EQ(mesh->GetNV(), external ? 32 : 24); + EXPECT_EQ(mesh->GetNE(), external ? 19 : 13); + EXPECT_EQ(mesh->GetNBE(), external ? 12 : 6); + const auto volumes = CountVolumeAttributes(*mesh); + EXPECT_EQ(volumes.at(1), 7); + EXPECT_EQ(volumes.at(2), 6); + EXPECT_EQ(volumes.contains(3), external); + if (external) EXPECT_EQ(volumes.at(3), 6); + const auto boundaries = CountBoundaryAttributes(*mesh); + EXPECT_EQ(boundaries.at(1), 6); + EXPECT_EQ(boundaries.contains(2), external); + if (external) EXPECT_EQ(boundaries.at(2), 6); + + cfg->mutate([](stroid::config::MeshConfig& value) { value.core_mapping = "spherified"; }); + auto legacy = stroid::topology::BuildSkeleton(*cfg); + EXPECT_EQ(legacy->GetNE(), external ? 13 : 7); + EXPECT_EQ(CountVolumeAttributes(*legacy).at(1), 1); + } +} + +TEST_F(stroidTest, MultiBlockCore_RejectsUnknownMappingAndInvalidGeometryConfiguration) { + auto cfg = MultiBlockConfiguration(2, 0, true); + cfg->mutate([](stroid::config::MeshConfig& value) { value.core_mapping = "not_a_core_mapping"; }); + EXPECT_THROW(stroid::topology::BuildSkeleton(*cfg), std::invalid_argument); + cfg = MultiBlockConfiguration(2, 0, true); + cfg->mutate([](stroid::config::MeshConfig& value) { value.r_core = value.r_star; }); + EXPECT_THROW(stroid::topology::BuildSkeleton(*cfg), std::invalid_argument); + cfg = MultiBlockConfiguration(2, 0, true); + cfg->mutate([](stroid::config::MeshConfig& value) { value.r_infinity = value.r_star; }); + EXPECT_THROW(stroid::topology::BuildSkeleton(*cfg), std::invalid_argument); + cfg = MultiBlockConfiguration(2, 0, true); + cfg->mutate([](stroid::config::MeshConfig& value) { value.flattening = 1.0; }); + EXPECT_THROW(stroid::topology::BuildSkeleton(*cfg), std::invalid_argument); +} + +TEST_F(stroidTest, MultiBlockCore_MapHasAffineInnerCubeAndContinuousSphericalInterface) { + auto cfg = MultiBlockConfiguration(4, 0, true); + const double radius = (*cfg)->r_core.value(); + for (int axis = 0; axis < 3; ++axis) { + for (double sign : {-1.0, 1.0}) { + for (double a : {-1.0, -0.4, 0.0, 0.6, 1.0}) { + for (double b : {-1.0, -0.3, 0.0, 0.7, 1.0}) { + mfem::Vector direction(3); + direction(axis) = sign; + direction((axis + 1) % 3) = a; + direction((axis + 2) % 3) = b; + mfem::Vector inner(direction); + inner *= radius / 2.0; + mfem::Vector expected(inner); + expected /= std::sqrt(3.0); + mfem::Vector mapped = TransformCopy(inner, *cfg, 1); + mapped -= expected; + EXPECT_LT(mapped.Norml2(), 2.0e-14); + for (double interfaceRadius : {radius / 2.0, radius}) { + mfem::Vector inside(direction), outside(direction); + inside *= interfaceRadius * (1.0 - 1.0e-8); + outside *= interfaceRadius * (1.0 + 1.0e-8); + mapped = TransformCopy(inside, *cfg, 1); + mapped -= TransformCopy(outside, *cfg, interfaceRadius == radius ? 2 : 1); + EXPECT_LT(mapped.Norml2(), 1.0e-7 * radius); + } + mfem::Vector coreInterface(direction); + coreInterface *= radius; + EXPECT_NEAR(TransformCopy(coreInterface, *cfg, 1).Norml2(), radius, 2.0e-14); + } + } + } + } + auto mesh = stroid::GenerateMesh(*cfg); + ASSERT_NE(mesh.mesh, nullptr); + ExpectCoreFaceContinuity(*mesh.mesh, 1); +} + +TEST_F(stroidTest, MultiBlockCore_ClosedGridSignedJacobiansAndSvdAcrossOrdersAndRefinements) { + for (int order = 1; order <= 6; ++order) { + for (int refinement = 0; refinement <= 2; ++refinement) { + SCOPED_TRACE("order=" + std::to_string(order) + " refinement=" + std::to_string(refinement)); + auto cfg = MultiBlockConfiguration(order, refinement, false); + auto mesh = stroid::GenerateMesh(*cfg); + ASSERT_NE(mesh.mesh, nullptr); + const int factor = 1 << (3 * refinement); + EXPECT_EQ(mesh.mesh->GetNE(), 13 * factor); + EXPECT_EQ(CountVolumeAttributes(*mesh.mesh).at(1), 7 * factor); + ExpectClosedGridCoreConditioning(*mesh.mesh, 1); + } + } +} + +TEST_F(stroidTest, MultiBlockCore_MapIsScaleInvariantBelowLegacyRadiusCutoff) { + constexpr double scale = 1.0e-15; + auto reference = MultiBlockConfiguration(2, 0, true); + reference->mutate([](stroid::config::MeshConfig& value) { value.r_infinity = 5.0; }); + auto scaled = MultiBlockConfiguration(2, 0, true); + scaled->mutate([](stroid::config::MeshConfig& value) { + value.r_core = 2.5e-16; + value.r_star = 1.0e-15; + value.r_infinity = 5.0e-15; + }); + const std::array, 10> points{{ + {{0.0, 0.0, 0.0}}, + {{0.05, -0.04, 0.1}}, + {{0.125, 0.08, -0.02}}, + {{0.18, -0.09, 0.12}}, + {{-0.2, -0.2, -0.2}}, + {{0.25, 0.12, -0.2}}, + {{0.6, -0.2, 0.4}}, + {{1.0, 0.7, -0.3}}, + {{3.0, -1.3, 0.4}}, + {{-5.0, 2.1, -1.0}} + }}; + for (const auto& coordinates : points) { + mfem::Vector point(3); + for (int component = 0; component < 3; ++component) point(component) = coordinates[component]; + const double logicalRadius = std::max({std::abs(point(0)), std::abs(point(1)), std::abs(point(2))}); + const int attribute = logicalRadius <= 0.25 ? 1 : logicalRadius <= 1.0 ? 2 : 3; + const auto expected = TransformCopy(point, *reference, attribute); + point *= scale; + auto actual = TransformCopy(point, *scaled, attribute); + actual /= scale; + for (int component = 0; component < 3; ++component) { + EXPECT_NEAR(actual(component), expected(component), 2.0e-13) + << "logical radius=" << logicalRadius << " component=" << component; + } + } +} + +TEST_F(stroidTest, MultiBlockCore_FlatteningCustomIdsAndExteriorCoordinateRemainConsistent) { + auto cfg = MultiBlockConfiguration(3, 1, true, 0.2); + cfg->mutate([](stroid::config::MeshConfig& value) { + value.core_id = 11; + value.envelope_id = 17; + value.vacuum_id = 23; + value.surface_bdr_id = 31; + value.inf_bdr_id = 37; + }); + auto mesh = stroid::GenerateMesh(*cfg); + ASSERT_NE(mesh.mesh, nullptr); + const auto volume = CountVolumeAttributes(*mesh.mesh); + EXPECT_EQ(volume.at(11), 7 * 8); + EXPECT_EQ(volume.at(17), 6 * 8); + EXPECT_EQ(volume.at(23), 6 * 8); + const auto boundary = CountBoundaryAttributes(*mesh.mesh); + EXPECT_EQ(boundary.at(31), 6 * 4); + EXPECT_EQ(boundary.at(37), 6 * 4); + ExpectClosedGridCoreConditioning(*mesh.mesh, 11); + ExpectCoreFaceContinuity(*mesh.mesh, 11); + ExpectExteriorCoordinateRange(mesh); + ExpectExteriorCoordinateBoundaryTraces(mesh); + mfem::Vector point(3); + point(0) = 0.25; + point(1) = 0.25; + point(2) = 0.25; + auto mapped = TransformCopy(point, *cfg, 11); + mapped(2) /= 0.8; + EXPECT_NEAR(mapped.Norml2(), 0.25, 2.0e-14); +} + +TEST_F(stroidTest, MultiBlockCore_OuterMappingAndSignedStellarVolumeMatchLegacy) { + auto cfg = MultiBlockConfiguration(3, 1, true); + auto legacyCfg = MultiBlockConfiguration(3, 1, true); + legacyCfg->mutate([](stroid::config::MeshConfig& value) { value.core_mapping = "spherified"; }); + const double coreRadius = (*cfg)->r_core.value(); + const double stellarRadius = (*cfg)->r_star.value(); + const double infinityRadius = (*cfg)->r_infinity.value(); + for (int axis = 0; axis < 3; ++axis) { + for (double sign : {-1.0, 1.0}) { + for (double a : {-1.0, -0.3, 0.0, 0.8, 1.0}) { + for (double b : {-1.0, 0.0, 0.4, 1.0}) { + mfem::Vector direction(3); + direction(axis) = sign; + direction((axis + 1) % 3) = a; + direction((axis + 2) % 3) = b; + for (double radius : {coreRadius, (coreRadius + stellarRadius) / 2.0, stellarRadius, + (stellarRadius + infinityRadius) / 2.0, infinityRadius}) { + mfem::Vector point(direction); + point *= radius; + const int attribute = radius <= stellarRadius ? 2 : 3; + auto difference = TransformCopy(point, *cfg, attribute); + difference -= TransformCopy(point, *legacyCfg, attribute); + EXPECT_LT(difference.Norml2(), 2.0e-14 * infinityRadius); + } + } + } + } + } + auto mesh = stroid::GenerateMesh(*cfg); + auto legacy = stroid::GenerateMesh(*legacyCfg); + const auto signedStellarVolume = [](mfem::Mesh& candidate) { + double volume = 0.0; + for (int element = 0; element < candidate.GetNE(); ++element) { + if (candidate.GetAttribute(element) == 3) continue; + auto* transformation = candidate.GetElementTransformation(element); + const auto& rule = mfem::IntRules.Get(transformation->GetGeometryType(), 3 * transformation->Order() + 2); + for (int q = 0; q < rule.GetNPoints(); ++q) { + const auto& point = rule.IntPoint(q); + transformation->SetIntPoint(&point); + volume += point.weight * transformation->Jacobian().Det(); + } + } + return volume; + }; + const double newVolume = signedStellarVolume(*mesh.mesh); + const double oldVolume = signedStellarVolume(*legacy.mesh); + EXPECT_GT(newVolume, 0.0); + EXPECT_NEAR(newVolume, oldVolume, 2.0e-11 * oldVolume); +} + +TEST_F(stroidTest, MultiBlockCore_SaveLoadConfigAndRefinementPreserveContracts) { + for (const bool external : {false, true}) { + SCOPED_TRACE(external); + auto cfg = MultiBlockConfiguration(3, 0, external); + auto original = stroid::GenerateMesh(*cfg); + EXPECT_EQ(original.type, stroid::MFEM_MESH_TYPE::SERIAL); + const auto path = std::filesystem::temp_directory_path() / + (external ? "stroid_multiblock_external_round_trip.smesh" : "stroid_multiblock_stellar_round_trip.smesh"); + stroid::IO::SaveStroidMesh(original, path.string(), "Multi-block core regression"); + auto result = stroid::IO::LoadStroidMesh(path.string()); + ASSERT_TRUE(result.has_value()) << result.error(); + auto loaded = std::move(*result); + EXPECT_EQ(loaded.type, stroid::MFEM_MESH_TYPE::SERIAL); + ASSERT_NE(loaded.mesh, nullptr); + ASSERT_NE(loaded.reference_mesh, nullptr); + EXPECT_EQ(loaded.config.core_mapping.value(), "multi_block"); + EXPECT_EQ(loaded.config.include_external_domain.value(), external); + EXPECT_EQ(loaded.mesh->GetNE(), original.mesh->GetNE()); + ASSERT_EQ(loaded.mesh->GetNodes()->Size(), original.mesh->GetNodes()->Size()); + for (int dof = 0; dof < original.mesh->GetNodes()->Size(); ++dof) { + EXPECT_NEAR((*loaded.mesh->GetNodes())(dof), (*original.mesh->GetNodes())(dof), 2.0e-14); + } + stroid::refinement::UniformRefinement(loaded, 1); + EXPECT_EQ(loaded.refinement_levels, 1); + EXPECT_EQ(loaded.mesh->GetNE(), original.mesh->GetNE() * 8); + EXPECT_EQ(CountVolumeAttributes(*loaded.mesh).at(1), 7 * 8); + ExpectClosedGridCoreConditioning(*loaded.mesh, 1); + ExpectCoreFaceContinuity(*loaded.mesh, 1); + if (external) { + ExpectExteriorCoordinateRange(loaded); + ExpectExteriorCoordinateBoundaryTraces(loaded); + } else { + EXPECT_EQ(loaded.exterior_coordinate, nullptr); + } + std::error_code error; + std::filesystem::remove(path, error); + EXPECT_FALSE(error); + } + + auto legacyCfg = MultiBlockConfiguration(2, 0, false); + legacyCfg->mutate([](stroid::config::MeshConfig& value) { value.core_mapping = "spherified"; }); + auto legacy = stroid::GenerateMesh(*legacyCfg); + EXPECT_EQ(legacy.type, stroid::MFEM_MESH_TYPE::SERIAL); + const auto path = std::filesystem::temp_directory_path() / "stroid_core_mapping_legacy_round_trip.smesh"; + stroid::IO::SaveStroidMesh(legacy, path.string(), "Legacy core mapping default regression"); + std::ifstream input(path); + std::string contents((std::istreambuf_iterator(input)), std::istreambuf_iterator()); + const auto marker = contents.find("\ncore_mapping:"); + ASSERT_NE(marker, std::string::npos); + const auto fieldStart = marker + 1; + const auto newline = contents.find('\n', fieldStart); + ASSERT_NE(newline, std::string::npos); + contents.erase(fieldStart, newline - fieldStart + 1); + std::istringstream legacyStream(contents); + auto restored = stroid::IO::ParseStroidMesh(legacyStream); + ASSERT_TRUE(restored.has_value()) << restored.error(); + EXPECT_EQ(restored->config.core_mapping.value(), "spherified"); + EXPECT_EQ(CountVolumeAttributes(*restored->mesh).at(1), 1); + std::error_code error; + std::filesystem::remove(path, error); + EXPECT_FALSE(error); +} diff --git a/tools/geometry_quality_experiment.cpp b/tools/geometry_quality_experiment.cpp new file mode 100644 index 0000000..f9948c6 --- /dev/null +++ b/tools/geometry_quality_experiment.cpp @@ -0,0 +1,348 @@ +#include "stroid/stroid.h" +#include "CLI/CLI.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + constexpr double infinity = std::numeric_limits::infinity(); + + struct SampleLocation { + int element = -1; + mfem::IntegrationPoint point; + std::string source = "none"; + }; + + struct ConditioningStats { + int elements = 0; + size_t samples = 0; + size_t nonpositive_samples = 0; + double min_det = infinity; + double min_sigma = infinity; + double max_condition = 0.0; + double min_scaled_jacobian = infinity; + double contraction_boundary = infinity; + SampleLocation det_location; + SampleLocation condition_location; + SampleLocation contraction_location; + }; + + class ScopedOutputRedirect { + std::streambuf* original; + public: + ScopedOutputRedirect() : original(std::cout.rdbuf(std::cerr.rdbuf())) {} + ~ScopedOutputRedirect() { std::cout.rdbuf(original); } + }; + + std::vector ClosedSamples(int grid_points) { + std::vector points; + for (int i = 0; i < grid_points; ++i) { + for (int j = 0; j < grid_points; ++j) { + for (int k = 0; k < grid_points; ++k) { + mfem::IntegrationPoint point; + point.Set3(static_cast(i) / (grid_points - 1), + static_cast(j) / (grid_points - 1), + static_cast(k) / (grid_points - 1)); + points.push_back(point); + } + } + } + for (const double offset : {0.005, 0.010885670927, 0.02}) { + for (int corner = 0; corner < 8; ++corner) { + mfem::IntegrationPoint point; + point.Set3((corner & 1) ? 1.0 - offset : offset, + (corner & 2) ? 1.0 - offset : offset, + (corner & 4) ? 1.0 - offset : offset); + points.push_back(point); + } + } + return points; + } + + double ColumnNorm(const mfem::DenseMatrix& matrix, int column) { + double norm_squared = 0.0; + for (int row = 0; row < 3; ++row) { + norm_squared += matrix(row, column) * matrix(row, column); + } + return std::sqrt(norm_squared); + } + + double ContractionBoundary(const mfem::DenseMatrix& jacobian, + const mfem::DenseMatrix& direction) { + std::array coefficients{}; + mfem::DenseMatrix mixed(3); + for (int mask = 0; mask < 8; ++mask) { + int degree = 0; + for (int column = 0; column < 3; ++column) { + const bool use_direction = (mask & (1 << column)) != 0; + degree += use_direction; + for (int row = 0; row < 3; ++row) { + mixed(row, column) = use_direction ? direction(row, column) : jacobian(row, column); + } + } + coefficients[degree] += mixed.Det(); + } + if (!(coefficients[0] > 0.0)) return 0.0; + + auto polynomial = [&](double x) { + return ((coefficients[3] * x + coefficients[2]) * x + coefficients[1]) * x + coefficients[0]; + }; + std::vector breaks{0.0, 1.0}; + auto add_break = [&](double x) { + if (std::isfinite(x) && x > 0.0 && x < 1.0) breaks.push_back(x); + }; + const double a = 3.0 * coefficients[3]; + const double b = 2.0 * coefficients[2]; + const double c = coefficients[1]; + if (a == 0.0) { + if (b != 0.0) add_break(-c / b); + } else { + const double discriminant = b * b - 4.0 * a * c; + if (discriminant >= 0.0) { + const double q = -0.5 * (b + std::copysign(std::sqrt(discriminant), b)); + if (q == 0.0) { + add_break(-b / (2.0 * a)); + } else { + add_break(q / a); + add_break(c / q); + } + } + } + std::sort(breaks.begin(), breaks.end()); + const double scale = std::abs(coefficients[0]) + std::abs(coefficients[1]) + + std::abs(coefficients[2]) + std::abs(coefficients[3]); + const double tolerance = 64.0 * std::numeric_limits::epsilon() * scale; + for (size_t i = 1; i < breaks.size(); ++i) { + double right = breaks[i]; + const double value = polynomial(right); + if (value > tolerance) continue; + if (std::abs(value) <= tolerance) return right; + double left = breaks[i - 1]; + for (int iteration = 0; iteration < 64; ++iteration) { + const double middle = 0.5 * (left + right); + if (polynomial(middle) > 0.0) left = middle; + else right = middle; + } + return right; + } + return infinity; + } + + std::unique_ptr BuildContractionProbe( + stroid::StroidMesh& mesh, mfem::FiniteElementSpace& space, double stellar_radius + ) { + auto values = std::make_unique(&space); + *values = 0.0; + std::vector processed(space.GetNDofs(), false); + mfem::Array dofs; + mfem::Vector physical(3), logical(3); + for (int element = 0; element < mesh.mesh->GetNE(); ++element) { + const auto& nodes = space.GetFE(element)->GetNodes(); + space.GetElementDofs(element, dofs); + auto* physical_transform = mesh.mesh->GetElementTransformation(element); + auto* logical_transform = mesh.reference_mesh->GetElementTransformation(element); + for (int local = 0; local < dofs.Size(); ++local) { + const int dof = dofs[local] >= 0 ? dofs[local] : -1 - dofs[local]; + if (processed[dof]) continue; + physical_transform->Transform(nodes.IntPoint(local), physical); + logical_transform->Transform(nodes.IntPoint(local), logical); + const double radius = physical.Norml2(); + const double logical_radius = std::max({std::abs(logical(0)), std::abs(logical(1)), std::abs(logical(2))}); + const double fraction = std::min(logical_radius / stellar_radius, 1.0); + for (int component = 0; component < 3; ++component) { + (*values)(space.DofToVDof(dof, component)) = radius > 0.0 + ? -stellar_radius * fraction * fraction * physical(component) / radius : 0.0; + } + processed[dof] = true; + } + } + return values; + } + + std::map InspectMesh(stroid::StroidMesh& mesh, int order, + int grid_points, bool contraction_probe, int probe_order) { + std::map result; + const auto closed_samples = ClosedSamples(grid_points); + mfem::H1_FECollection probe_collection(probe_order, 3); + mfem::FiniteElementSpace probe_space(mesh.mesh.get(), &probe_collection, 3); + std::unique_ptr probe; + if (contraction_probe) probe = BuildContractionProbe(mesh, probe_space, mesh.config.r_star.value()); + mfem::DenseMatrix probe_values, probe_shape, direction(3); + mfem::Array probe_dofs; + + for (int element = 0; element < mesh.mesh->GetNE(); ++element) { + const int attribute = mesh.mesh->GetAttribute(element); + ++result[attribute].elements; + ++result[0].elements; + auto* transform = mesh.mesh->GetElementTransformation(element); + const bool inspect_probe = contraction_probe && attribute != static_cast(mesh.config.vacuum_id.value()); + if (inspect_probe) { + probe_space.GetElementDofs(element, probe_dofs); + probe_values.SetSize(3, probe_dofs.Size()); + probe_shape.SetSize(probe_dofs.Size(), 3); + for (int local = 0; local < probe_dofs.Size(); ++local) { + const int dof = probe_dofs[local] >= 0 ? probe_dofs[local] : -1 - probe_dofs[local]; + for (int component = 0; component < 3; ++component) { + probe_values(component, local) = (*probe)(probe_space.DofToVDof(dof, component)); + } + } + } + auto inspect_point = [&](const mfem::IntegrationPoint& point, const std::string& source) { + transform->SetIntPoint(&point); + const mfem::DenseMatrix& jacobian = transform->Jacobian(); + const double determinant = jacobian.Det(); + const double sigma_min = jacobian.CalcSingularvalue(2); + const double sigma_max = jacobian.CalcSingularvalue(0); + const double condition = sigma_min > 0.0 ? sigma_max / sigma_min : infinity; + const double denominator = ColumnNorm(jacobian, 0) * ColumnNorm(jacobian, 1) * ColumnNorm(jacobian, 2); + const double scaled_jacobian = denominator > 0.0 ? determinant / denominator : 0.0; + double boundary = infinity; + if (inspect_probe) { + probe_space.GetFE(element)->CalcDShape(point, probe_shape); + mfem::Mult(probe_values, probe_shape, direction); + boundary = ContractionBoundary(jacobian, direction); + } + for (const int region : {0, attribute}) { + auto& stats = result[region]; + ++stats.samples; + if (!(determinant > 0.0)) ++stats.nonpositive_samples; + if (determinant < stats.min_det) { + stats.min_det = determinant; + stats.det_location = {element, point, source}; + } + stats.min_sigma = std::min(stats.min_sigma, sigma_min); + stats.min_scaled_jacobian = std::min(stats.min_scaled_jacobian, scaled_jacobian); + if (condition > stats.max_condition) { + stats.max_condition = condition; + stats.condition_location = {element, point, source}; + } + if (boundary < stats.contraction_boundary) { + stats.contraction_boundary = boundary; + stats.contraction_location = {element, point, source}; + } + } + }; + const auto& quadrature = mfem::IntRules.Get(transform->GetGeometryType(), 2 * order + 4); + for (int point = 0; point < quadrature.GetNPoints(); ++point) { + inspect_point(quadrature.IntPoint(point), "quadrature"); + } + for (const auto& point : closed_samples) inspect_point(point, "closed_grid_and_corner_probes"); + } + return result; + } + + void WriteLocation(std::ostream& output, const SampleLocation& location, stroid::StroidMesh& mesh) { + output << ',' << location.element << ',' << location.source; + if (location.element < 0) { + output << ",nan,nan,nan,nan,nan,nan,nan,nan,nan"; + return; + } + mfem::Vector physical(3), logical(3); + mesh.mesh->GetElementTransformation(location.element)->Transform(location.point, physical); + mesh.reference_mesh->GetElementTransformation(location.element)->Transform(location.point, logical); + output << ',' << location.point.x << ',' << location.point.y << ',' << location.point.z; + for (int component = 0; component < 3; ++component) output << ',' << physical(component); + for (int component = 0; component < 3; ++component) output << ',' << logical(component); + } +} + +int main(int argc, char** argv) { + std::vector orders{1, 2, 3, 4, 5, 6}; + std::vector refinements{0, 1, 2}; + std::vector mappings{"spherified", "multi_block"}; + std::string output_path; + int grid_points = 5; + int probe_order = 3; + double core_radius = 0.25; + double infinity_radius = 5.0; + double flattening = 0.0; + bool no_external = false; + bool contraction_probe = false; + CLI::App app{"Compare signed Jacobians and conditioning of the actual high-order STROID mesh; TMOP is disabled."}; + app.add_option("--orders", orders, "Geometry orders, comma separated")->delimiter(',')->check(CLI::Range(1, 8)); + app.add_option("--refinements", refinements, "Uniform refinement levels, comma separated")->delimiter(',')->check(CLI::Range(0, 3)); + app.add_option("--mappings", mappings, "Core mappings, comma separated")->delimiter(',')->check(CLI::IsMember({"spherified", "multi_block"})); + app.add_option("--grid-points", grid_points, "Closed tensor grid points per coordinate, plus near-corner probes")->check(CLI::Range(2, 15)); + app.add_option("--core-radius", core_radius, "Core radius; stellar radius is one"); + app.add_option("--infinity-radius", infinity_radius, "Outer reference radius"); + app.add_option("--flattening", flattening, "Spheroidal flattening"); + app.add_option("--output", output_path, "New CSV output file; defaults to stdout"); + app.add_flag("--no-external", no_external, "Omit exterior domain"); + app.add_flag("--contraction-probe", contraction_probe, "Inspect an interpolated unit logical-radius-squared radial contraction in stellar elements"); + app.add_option("--probe-order", probe_order, "H1 displacement order for the optional contraction probe")->check(CLI::Range(1, 8)); + try { + app.parse(argc, argv); + } catch (const CLI::ParseError& error) { + return app.exit(error); + } + + try { + if (!std::isfinite(core_radius) || core_radius <= 0.0 || core_radius >= 1.0 + || !std::isfinite(infinity_radius) || infinity_radius <= 1.0 + || !std::isfinite(flattening) || flattening < 0.0 || flattening >= 1.0) { + throw std::invalid_argument("Require 0 < core-radius < 1 < infinity-radius and 0 <= flattening < 1."); + } + std::ofstream file; + if (!output_path.empty()) { + if (std::filesystem::exists(output_path)) throw std::runtime_error("Refusing to overwrite existing output: " + output_path); + file.open(output_path); + if (!file) throw std::runtime_error("Could not open output: " + output_path); + } + std::ostream& output = output_path.empty() ? std::cout : file; + output << std::setprecision(17); + output << "mapping,order,refinement,r_core,r_star,r_infinity,flattening,external,grid_points,quadrature_order,probe_order,attribute,elements,samples,nonpositive_samples,min_signed_det,min_sigma,max_condition,min_scaled_jacobian,contraction_boundary_up_to_one"; + for (const std::string prefix : {"det", "condition", "contraction"}) { + output << ',' << prefix << "_element," << prefix << "_source," << prefix << "_xi," << prefix << "_eta," << prefix << "_zeta," << prefix << "_x," << prefix << "_y," << prefix << "_z," << prefix << "_logical_x," << prefix << "_logical_y," << prefix << "_logical_z"; + } + output << '\n'; + std::cerr << "Sampling actual FE geometry, not the analytical map. Attribute 0 aggregates all regions.\n" + "Signed determinants and scaled Jacobians retain orientation; inf boundary means no sampled root through alpha=1.\n" + "The optional contraction probe is a diagnostic field, not a Newton correction or a production exterior extension.\n"; + for (const auto& mapping : mappings) { + for (const int order : orders) { + for (const int refinement : refinements) { + stroid::config::MeshConfig config; + config.core_mapping = mapping; + config.order = order; + config.refinement_levels = refinement; + config.r_core = core_radius; + config.r_star = 1.0; + config.r_infinity = infinity_radius; + config.flattening = flattening; + config.include_external_domain = !no_external; + config.optimization_methods = stroid::config::OptimizationMethods{false, true}; + std::cerr << "Inspecting " << mapping << ", order " << order << ", refinement " << refinement << '\n'; + stroid::StroidMesh mesh; + { + ScopedOutputRedirect redirect; + mesh = stroid::GenerateMesh(config); + } + const auto regions = InspectMesh(mesh, order, grid_points, contraction_probe, probe_order); + for (const auto& [attribute, stats] : regions) { + output << mapping << ',' << order << ',' << refinement << ',' << core_radius << ",1," << infinity_radius << ',' << flattening << ',' << !no_external << ',' << grid_points << ',' << 2 * order + 4 << ',' << (contraction_probe ? probe_order : 0) << ',' << attribute << ',' << stats.elements << ',' << stats.samples << ',' << stats.nonpositive_samples << ',' << stats.min_det << ',' << stats.min_sigma << ',' << stats.max_condition << ',' << stats.min_scaled_jacobian << ',' << stats.contraction_boundary; + WriteLocation(output, stats.det_location, mesh); + WriteLocation(output, stats.condition_location, mesh); + WriteLocation(output, stats.contraction_location, mesh); + output << '\n'; + } + output.flush(); + if (!output) throw std::runtime_error("Failed to write experiment output."); + } + } + } + } catch (const std::exception& error) { + std::cerr << "Geometry quality experiment failed: " << error.what() << '\n'; + return 1; + } + return 0; +} diff --git a/tools/meson.build b/tools/meson.build index 379f7b9..33d5d0c 100644 --- a/tools/meson.build +++ b/tools/meson.build @@ -1 +1,4 @@ executable('stroid', 'stroid.cpp', dependencies: [stroid_dep, cli11_dep, magic_enum_dep], install: true) + +# Opt-in diagnostic driver; deliberately not part of the installed API/tools. +executable('geometry_quality_experiment', 'geometry_quality_experiment.cpp', dependencies: [stroid_dep, cli11_dep], build_by_default: false, install: false) diff --git a/tools/python/visualize_winding.py b/tools/python/visualize_winding.py deleted file mode 100644 index d722639..0000000 --- a/tools/python/visualize_winding.py +++ /dev/null @@ -1,46 +0,0 @@ -import matplotlib.pyplot as plt - -class Box: - def __init__(self, scale, vc="red", ec="blue", offset=0): - self.scale = scale - self.offset = offset - self.verticies = [[scale, -scale, -scale], [scale, scale, -scale], [-scale, -scale, -scale], [-scale, scale, -scale], [scale, -scale, scale], [scale, scale, scale], [-scale, -scale, scale], [-scale, scale, scale]] - self.edges = [[0, 1], [0, 2], [0, 4], [1, 3], [1, 5], [3, 7], [3, 2], [2, 6], [4, 5], [4, 6], [5, 7], [7, 6]] - self.edge_color = ec - self.vertex_color = vc - - def plot(self, ax, vlabel=True, elabel=False): - for vertID, vert in enumerate(self.verticies): - ax.scatter(vert[0], vert[1], vert[2], c=self.vertex_color) - if vlabel: - ax.text(vert[0], vert[1], vert[2], f"{self.offset + vertID}", fontsize=25) - for edge in self.edges: - ax.plot([self.verticies[edge[0]][0], self.verticies[edge[1]][0]], [self.verticies[edge[0]][1], self.verticies[edge[1]][1]], [self.verticies[edge[0]][2], self.verticies[edge[1]][2]], color=self.edge_color) - -class Wedge: - def __init__(self, A, B, ec="green"): - self.A = A - self.B = B - self.edge_color = ec - def plot(self, ax): - for vA, vB in zip(self.A.verticies, self.B.verticies): - ax.plot([vA[0], vB[0]], [vA[1], vB[1]], [vA[2], vB[2]], color=self.edge_color) - -def main(): - core = Box(0.5) - envelope = Box(2, offset=8) - star = Wedge(core, envelope) - infinity = Box(5, offset=16) - vacuum = Wedge(envelope, infinity) - - fig, ax = plt.subplots(1, 1, figsize=(10, 10), subplot_kw={"projection": "3d"}) - core.plot(ax) - envelope.plot(ax) - star.plot(ax) - infinity.plot(ax) - vacuum.plot(ax) - ax.view_init(30, 30) - plt.show() - -if __name__ == "__main__": - main() diff --git a/utils/build-wheels-macos_aarch64.sh b/utils/build-wheels-macos_aarch64.sh index a702a83..1453119 100755 --- a/utils/build-wheels-macos_aarch64.sh +++ b/utils/build-wheels-macos_aarch64.sh @@ -80,14 +80,10 @@ for PY_VERSION in "${PYTHON_VERSIONS[@]}"; do echo "➤ Repairing wheel with delocate" if [[ -n "${FOURDST_PIN}" ]]; then - FOURDST_LIB_PATH="$("$PY" -c 'import fourdst, os; print(os.pathsep.join(fourdst.get_lib_dirs()))')" - DELOCATE_DYLD_PATH="${FOURDST_LIB_PATH}:${DELOCATE_DYLD_PATH}" - DYLD_LIBRARY_PATH="${DELOCATE_DYLD_PATH}" \ delocate-wheel --require-archs arm64 \ -e composition -e logging -e const -e reflect_cpp \ -w "${FINAL_WHEEL_DIR}" -v "$CURRENT_WHEEL" else - DYLD_LIBRARY_PATH="${DELOCATE_DYLD_PATH}" \ delocate-wheel --require-archs arm64 \ -w "${FINAL_WHEEL_DIR}" -v "$CURRENT_WHEEL" fi @@ -109,4 +105,4 @@ done rm -rf "${TMPDIR}" rm -rf "${WHEEL_DIR}" -echo "✅ All builds complete. Artifacts in ${FINAL_WHEEL_DIR}" \ No newline at end of file +echo "All builds complete. Artifacts in ${FINAL_WHEEL_DIR}" \ No newline at end of file