diff --git a/Doxyfile b/Doxyfile index abcc940..21fb150 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.5.0 +PROJECT_NUMBER = v0.6.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_NC.png b/assets/imgs/ExampleMesh_NC.png new file mode 100644 index 0000000..647a85a Binary files /dev/null and b/assets/imgs/ExampleMesh_NC.png differ diff --git a/assets/imgs/ExampleMesh_multi-block.png b/assets/imgs/ExampleMesh_multi-block.png index 59d523c..b43f65a 100644 Binary files a/assets/imgs/ExampleMesh_multi-block.png and b/assets/imgs/ExampleMesh_multi-block.png differ diff --git a/assets/imgs/ExampleMesh_spherified.png b/assets/imgs/ExampleMesh_spherified.png index d07c0b4..aae633d 100644 Binary files a/assets/imgs/ExampleMesh_spherified.png and b/assets/imgs/ExampleMesh_spherified.png differ diff --git a/configs/nonconforming_vacuum.toml b/configs/nonconforming_vacuum.toml new file mode 100644 index 0000000..0a3f610 --- /dev/null +++ b/configs/nonconforming_vacuum.toml @@ -0,0 +1,17 @@ +[main] +# Absolute minimum depths from the initial block topology. +refinement_levels = 4 +vacuum_refinement_levels = 2 +# Omit to inherit refinement_levels at the vacuum outer boundary. +# vacuum_outer_refinement_levels = 4 +order = 3 +include_external_domain = true +core_mapping = "multi_block" +r_core = 0.25 +r_star = 1.0 +r_infinity = 6.0 +flattening = 0.0 + +[main.optimization_methods] +tmop = false +smoothstep = true diff --git a/meson.build b/meson.build index 0c15932..45580a9 100644 --- a/meson.build +++ b/meson.build @@ -1,4 +1,4 @@ -project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.5.0', default_options : ['cpp_std=c++23']) +project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.6.0', default_options : ['cpp_std=c++23']) subdir('build-check') diff --git a/readme.md b/readme.md index 58955fb..8455ef5 100644 --- a/readme.md +++ b/readme.md @@ -103,11 +103,11 @@ smoothstep = true | Parameter | Description | Default | |---------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|---------------| -| refinement_levels | Number of uniform refinement levels to apply to the mesh after generation | 4 | +| refinement_levels | Stellar minimum depth, or uniform depth when vacuum overrides are omitted | 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 | +| r_core | The radius of the core region of the star | 0.25 | +| r_star | The radius of the star | 1.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 | @@ -123,7 +123,11 @@ smoothstep = true 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. +need only include parameters that differ from the defaults. For compatibility with older TOML files, +an omitted `core_mapping` uses `"spherified"`, and omitted TMOP controls leave optimization disabled. +Set `core_mapping = "multi_block"` explicitly to use the conditioned mapping in a TOML file. +Default-constructed C++ and Python `MeshConfig` objects select `"multi_block"`; other omitted TOML +geometry fields use the defaults from `MeshConfig`. ### Conditioned core mapping @@ -148,32 +152,98 @@ build/tools/geometry_quality_experiment --orders 4 --refinements 2 \ --contraction-probe --probe-order 3 --output core_comparison.csv ``` +### Nonconforming vacuum refinement + +Stroid can keep the star and both ends of the vacuum well resolved while using +coarser elements in the vacuum interior. Refinement is isotropic: each refinement +splits a hexahedron into eight children. Note however that only one geometric polynomial `order` applies +to every region. + +```toml +[main] +refinement_levels = 4 +vacuum_refinement_levels = 2 +# Optional: omitted outer depth inherits refinement_levels (4 here). +# vacuum_outer_refinement_levels = 4 +order = 3 +include_external_domain = true +core_mapping = "multi_block" + +[main.optimization_methods] +tmop = false +smoothstep = true +``` + +`configs/nonconforming_vacuum.toml` provides a complete example. The three depth +settings are absolute minimum targets measured from the initial block topology: + +| Setting | Applies to | Default | +|----------------------------------|------------------------------------------|-----------------------------| +| `refinement_levels` | Core and envelope | `4` | +| `vacuum_refinement_levels` | Vacuum interior | Inherit `refinement_levels` | +| `vacuum_outer_refinement_levels` | Cells touching the vacuum outer boundary | Inherit `refinement_levels` | + +Omitting both vacuum overrides preserves uniform generation. Supplying either +activates the local refinement policy and requires `include_external_domain = true`. +All levels must be nonnegative integers. + +Stroid enforces that vacuum cells touching the stellar surface match the stellar face subdivision. That is to say that +the inner boundary of the vacuum region is conforming to the outer boundary of the stellar region. Further, the +outer-boundary cells receive the outer target, and automatic grading limits neighboring refinement depths to one level. +This two layer approach is intended to allow for refinement when using compactification maps. + +```python +import stroid + +cfg = stroid.config.MeshConfig( + refinement_levels=4, + vacuum_refinement_levels=2, + vacuum_outer_refinement_levels=None, # Inherit stellar depth. + order=3, + core_mapping="multi_block", + optimization_methods=stroid.config.OptimizationMethods(tmop=False), +) +mesh = stroid.GenerateMesh(cfg) +features = stroid.stats.MESH_STAT_DEFAULT | stroid.stats.MeshStatFeatures.ELEMENT_COUNT +stats = stroid.stats.ComputeMeshStats(mesh, features) +print(stats.element_counts.vacuum) +print(stats.refinement.vacuum.min_depth, stats.refinement.vacuum.max_depth) +print(stats.refinement.geometry_dofs, stats.refinement.geometry_true_dofs) +print(stats.conformity.conforming, stats.conformity.n_nonconforming_faces) + +stroid.IO.SaveStroidMesh(mesh, "graded.stroid") +restored = stroid.IO.LoadStroidMesh("graded.stroid") +stroid.refinement.UniformRefinement(restored, 1) +``` + + +The `UniformRefinement(mesh, n)` function adds `n` levels to every current leaf while preserving the +existing grading, and rebuilds the geometry and exterior coordinate. Note that this means that a non-conforming mesh +that has been Uniformly refined will still be non-conforming, but the refinement will be applied to all leaves. + +#### Viewing curved meshes in GLVis + +It is important to note --- and potentially confusing if not understood --- that GLVis approximates curved faces with +flat triangles. At a hanging interface, the same subdivision count on a coarse face and its finer neighbors samples the +curved surface at different locations. This can produce apparent gaps even when the finite-element face transformations +agree. These gaps are not indications that the mesh itself is non-conforming; rather, they are a visualization artifact. + ### C++ Interface Stroid can be used as a library in C++ projects. After installation, include the stroid header and link against the stroid library. A basic example of using stroid in C++ is shown below (note that you will need a glvis instance running on localhost:19916 to visualize the mesh): ```c++ -#include -#include "mfem.hpp" - -#include "stroid/config/config.h" -#include "stroid/IO/mesh.h" -#include "stroid/topology/curvilinear.h" -#include "stroid/topology/topology.h" - -#include "fourdst/config/config.h" +#include "stroid/stroid.h" int main() { - const fourdst::config::Config cfg; + stroid::config::MeshConfig cfg; + cfg.refinement_levels = 4; + cfg.vacuum_refinement_levels = 2; + cfg.optimization_methods = stroid::config::OptimizationMethods{false, true}; - const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); - stroid::topology::Finalize(*mesh, cfg); - stroid::topology::PromoteToHighOrder(*mesh, cfg); - stroid::topology::ProjectMesh(*mesh, cfg); - stroid::topology::OptimizeMesh(*mesh, cfg); - - - stroid::IO::ViewMesh(*mesh, "Spheroidal Mesh", stroid::IO::VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID); + auto mesh = stroid::GenerateMesh(cfg); + stroid::IO::SaveStroidMesh(mesh, "graded.stroid"); + stroid::IO::ViewMesh(mesh, "Spheroidal Mesh", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, "localhost", 19916); } ``` @@ -184,8 +254,12 @@ An example mesh with the default configuration parameters is shown below (colora 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. +An example of a non-conforming mesh generated with stroid. Note that the gaps between elements are a visualization artifact +rather than true gaps within the mesh. +![Non Conforming Mesh](assets/imgs/ExampleMesh_NC.png) + +Note that both of these meshes are shown with 3 levels of refinement and polynomial order 3. Blue shows the core +domain, yellow shows the envelope domain, while purple shows the vacuum domain. ## Funding diff --git a/src/include/stroid/IO/mesh.h b/src/include/stroid/IO/mesh.h index 0fce639..fb76cbf 100644 --- a/src/include/stroid/IO/mesh.h +++ b/src/include/stroid/IO/mesh.h @@ -2,6 +2,7 @@ #include #include #include +#include #include "mfem.hpp" @@ -57,6 +58,14 @@ namespace stroid::IO { */ void SaveVTU(const stroid::StroidMesh& mesh, const std::string& exportName); + /** + * @brief Make a display-only mesh copy with matching face subdivisions. + * This is purely for visualization and should not be used for any science goals. + * @param mesh Source mesh whose geometry is to be displayed. + * @return Independently owned copy with no hanging faces. + */ + std::unique_ptr MakeConformingVisualizationMesh(const mfem::Mesh& mesh); + /** * @brief Stream a mesh to a running GLVis server for interactive viewing. * @param mesh Mesh to display. @@ -64,18 +73,27 @@ namespace stroid::IO { * @param mode Attribute visualization mode. * @param vishost GLVis server host. * @param visport GLVis server port. + * @param conforming_display Refine a display-only copy at hanging interfaces + * to prevent GLVis tessellation gaps. Set false to inspect the original + * element layout, which can show rendering gaps on curved interfaces. + * */ - void ViewMesh(mfem::Mesh &mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport); + void ViewMesh(mfem::Mesh &mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport, bool conforming_display=true); - void ViewMesh(const stroid::StroidMesh& mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport); + void ViewMesh(const stroid::StroidMesh& mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport, bool conforming_display=true); /** - * @brief Visualize boundary face valence (1=surface, 2=internal). + * @brief Color boundary-adjacent elements by face valence (1=surface, 2=internal). + * Untagged elements are zero; elements touching several tagged faces use + * the maximum valence. Values are computed before display subdivision. * @param mesh Mesh whose boundary faces are inspected. + * @param vishost GLVis server host. + * @param visport GLVis server port. + * @param conforming_display Use the same display-only refinement as ViewMesh. */ - void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport); + void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport, bool conforming_display=true); - void VisualizeFaceValence(const stroid::StroidMesh& mesh, const std::string &vishost, int visport); + void VisualizeFaceValence(const stroid::StroidMesh& mesh, const std::string &vishost, int visport, bool conforming_display=true); std::expected ParseStroidMesh(std::istream& is); std::expected LoadStroidMesh(const std::string& filename); @@ -84,4 +102,4 @@ namespace stroid::IO { std::expected ParseStroidMesh(std::istream& is, MPI_Comm comm); std::expected LoadStroidMesh(const std::string& filename, MPI_Comm comm); #endif -} \ No newline at end of file +} diff --git a/src/include/stroid/config/config.h b/src/include/stroid/config/config.h index 999576f..1f14ef0 100644 --- a/src/include/stroid/config/config.h +++ b/src/include/stroid/config/config.h @@ -19,15 +19,35 @@ namespace stroid::config { * These values are typically loaded via * `fourdst::config::Config` from a TOML file. * The README shows the expected TOML layout under the `[main]` table. - * Unspecified keys use the defaults defined here. + * Unspecified geometry keys use the defaults defined here. ResolveDefaults preserves + * the historical fallback for omitted mapping and optimization controls in TOML files. */ struct MeshConfig { /** - * @brief Number of uniform refinement passes applied after topology creation. + * @brief Stellar refinement depth, or uniform depth when vacuum overrides are absent. * @section toml * - [main].refinement_levels */ std::optional refinement_levels = 4; + /** + * @brief Minimum refinement depth in the vacuum interior; unset inherits `refinement_levels`. + * + * Setting either vacuum override enables local isotropic refinement. Vacuum cells at the + * stellar surface match the stellar refinement, and automatic one-level grading can + * raise the interior depth above this minimum. Requires an external domain. + * @section toml + * - [main].vacuum_refinement_levels + */ + std::optional vacuum_refinement_levels = std::nullopt; + /** + * @brief Minimum refinement depth at the vacuum outer boundary; unset inherits `refinement_levels`. + * + * Boundary-adjacent cells are refined automatically, with grading toward the vacuum + * interior. Geometry uses the same polynomial order in every region. + * @section toml + * - [main].vacuum_outer_refinement_levels + */ + std::optional vacuum_outer_refinement_levels = std::nullopt; /** * @brief Polynomial order for high-order elements. * @section toml @@ -139,6 +159,39 @@ namespace stroid::config { }; + /** + * @brief Fill omitted configuration values. + */ + inline MeshConfig ResolveDefaults(const MeshConfig& mesh_config) { + const MeshConfig defaults; + MeshConfig resolved = mesh_config; + auto resolve = [](auto& value, const auto& default_value) { + if (!value.has_value()) value = default_value; + }; + + resolve(resolved.refinement_levels, defaults.refinement_levels); + resolve(resolved.order, defaults.order); + resolve(resolved.include_external_domain, defaults.include_external_domain); + resolve(resolved.r_core, defaults.r_core); + resolve(resolved.r_star, defaults.r_star); + resolve(resolved.flattening, defaults.flattening); + resolve(resolved.r_infinity, defaults.r_infinity); + resolve(resolved.r_instability, defaults.r_instability); + resolve(resolved.core_steepness, defaults.core_steepness); + resolve(resolved.continuity_order, defaults.continuity_order); + resolve(resolved.surface_bdr_id, defaults.surface_bdr_id); + resolve(resolved.inf_bdr_id, defaults.inf_bdr_id); + resolve(resolved.core_id, defaults.core_id); + resolve(resolved.envelope_id, defaults.envelope_id); + resolve(resolved.vacuum_id, defaults.vacuum_id); + resolved.core_mapping = resolved.core_mapping.value_or("spherified"); + resolved.optimization_methods = resolved.optimization_methods.value_or(OptimizationMethods{}); + resolved.optimization_methods->tmop = resolved.optimization_methods->tmop.value_or(false); + resolved.optimization_methods->smoothstep = resolved.optimization_methods->smoothstep.value_or(true); + + return resolved; + } + inline std::string to_string(const MeshConfig &mesh_config) { auto opt_2_string = [](const OptimizationMethods& opt) { std::stringstream ss; @@ -160,6 +213,8 @@ namespace stroid::config { ss << "MeshConfig:\n"; ss << std::format(" refinement_levels: {}\n", mesh_config.refinement_levels.value_or(4)); + ss << std::format(" vacuum_refinement_levels: {}\n", mesh_config.vacuum_refinement_levels.has_value() ? std::to_string(*mesh_config.vacuum_refinement_levels) : "inherit"); + ss << std::format(" vacuum_outer_refinement_levels: {}\n", mesh_config.vacuum_outer_refinement_levels.has_value() ? std::to_string(*mesh_config.vacuum_outer_refinement_levels) : "inherit"); 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)); diff --git a/src/include/stroid/refinement/uniform.h b/src/include/stroid/refinement/uniform.h index 6b028dd..e01a5a5 100644 --- a/src/include/stroid/refinement/uniform.h +++ b/src/include/stroid/refinement/uniform.h @@ -3,5 +3,11 @@ #include "stroid/utils/types.h" namespace stroid::refinement { + /** + * @brief Refine every current leaf, preserving any existing nonconforming hierarchy. + * Rebuilds constrained geometry and the exterior coordinate from the refined + * reference mesh. The saved generation configuration is unchanged; the + * refinement counter and actual regional depths increase by @p levels. + */ void UniformRefinement(StroidMesh& mesh, size_t levels); -} \ No newline at end of file +} diff --git a/src/include/stroid/stroid.h b/src/include/stroid/stroid.h index 3864ccc..ce7be45 100644 --- a/src/include/stroid/stroid.h +++ b/src/include/stroid/stroid.h @@ -49,7 +49,11 @@ * @endcode */ namespace stroid { - inline StroidMesh GenerateMesh(const fourdst::config::Config& cfg) { + inline StroidMesh GenerateMesh(const fourdst::config::Config& input) { + fourdst::config::Config cfg; + cfg.mutate([&input](config::MeshConfig& value) { + value = config::ResolveDefaults(*input); + }); StroidMesh sm; sm.type = MFEM_MESH_TYPE::SERIAL; sm.config = *cfg; @@ -59,9 +63,7 @@ namespace stroid { sm.reference_mesh = std::move(reference); sm.mesh = utils::BuildProjected(*sm.reference_mesh, cfg); - 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); - } + stroid::topology::OptimizeMesh(*sm.mesh, cfg); sm.exterior_coordinate = stroid::topology::BuildExteriorCoordinate(*sm.mesh, *sm.reference_mesh, cfg); return sm; } diff --git a/src/include/stroid/topology/topology.h b/src/include/stroid/topology/topology.h index 3327893..a45253d 100644 --- a/src/include/stroid/topology/topology.h +++ b/src/include/stroid/topology/topology.h @@ -15,9 +15,11 @@ namespace stroid::topology { */ std::unique_ptr BuildSkeleton(const fourdst::config::Config & config); /** - * @brief Finalize topology, validate orientation, and apply uniform refinement. + * @brief Finalize topology, validate orientation, and apply the configured refinement policy. * @param mesh Mesh to finalize in-place. - * @param config Mesh configuration (uses `refinement_levels`). + * @param config Mesh configuration. Vacuum refinement overrides enable a + * balanced hierarchy with fine layers at both vacuum boundaries and a + * conforming stellar interface. Without overrides, refinement is uniform. */ void Finalize(mfem::Mesh& mesh, const fourdst::config::Config &config); } diff --git a/src/include/stroid/utils/mesh_stats.h b/src/include/stroid/utils/mesh_stats.h index 32c5f45..2dff3d0 100644 --- a/src/include/stroid/utils/mesh_stats.h +++ b/src/include/stroid/utils/mesh_stats.h @@ -25,6 +25,7 @@ namespace stroid::stats { CENTROID = 1u << 10, CONFIG_META = 1u << 11, BOUNDING_BOX = 1u << 12, + REFINEMENT = 1u << 13, }; constexpr MeshStatFeatures operator|(MeshStatFeatures lhs, MeshStatFeatures rhs) { @@ -41,7 +42,7 @@ namespace stroid::stats { inline constexpr MeshStatFeatures MESH_STAT_DEFAULT = MeshStatFeatures::RADIUS | MeshStatFeatures::AXES | MeshStatFeatures::ELLIPTICITY | - MeshStatFeatures::CONFORMITY | MeshStatFeatures::CONFIG_META; + MeshStatFeatures::CONFORMITY | MeshStatFeatures::CONFIG_META | MeshStatFeatures::REFINEMENT; inline constexpr auto MESH_STAT_ALL = static_cast(0xFFFFFFFFu); @@ -70,9 +71,27 @@ namespace stroid::stats { struct ConformityStats { bool conforming = true; + bool hierarchy_enabled = false; + // Fine patches are counted once; their coarse master faces are excluded. long n_nonconforming_faces = 0; }; + struct RegionRefinementStats { + // An absent region has both depths set to -1. + int min_depth = -1; + int max_depth = -1; + }; + + struct RefinementStats { + RegionRefinementStats all; + RegionRefinementStats core; + RegionRefinementStats envelope; + RegionRefinementStats vacuum; + // Scalar nodal counts, independent of the coordinate vector dimension. + long geometry_dofs = 0; + long geometry_true_dofs = 0; + }; + struct JacobianStats { double detJ_min; double detJ_max; @@ -138,6 +157,7 @@ namespace stroid::stats { std::optional ellipticity; std::optional bowing; std::optional conformity; + std::optional refinement; std::optional jacobian; std::optional jacobian_stellar; std::optional jacobian_vacuum; diff --git a/src/lib/IO/mesh.cpp b/src/lib/IO/mesh.cpp index 285f311..2fd871f 100644 --- a/src/lib/IO/mesh.cpp +++ b/src/lib/IO/mesh.cpp @@ -41,7 +41,7 @@ namespace stroid::IO { # - 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 +# - refinement-levels : stellar baseline depth plus subsequent uniform passes; local depths are stored in the mesh hierarchy # NOTE: EACH BLOCK OF DATA IS STORED BETWEEN "BEGIN BLOCK \n ... \nEND BLOCK # PARSING THE UNDERLYING MFEM NATIVE MESH FORMAT CAN BE DONE WITH MFEM'S STREAM READER # IF YOU EXTRACT THE RAW CONTENTS BETWEEN THOSE LINES @@ -102,6 +102,16 @@ END BLOCK HEADER)", # default: 4 refinement_levels:{} +# vacuum_refinement_levels: Minimum vacuum interior depth; inherit uses refinement_levels +# std::optional +# default: inherit +vacuum_refinement_levels:{} + +# vacuum_outer_refinement_levels: Minimum vacuum outer-boundary depth; inherit uses refinement_levels +# std::optional +# default: inherit +vacuum_outer_refinement_levels:{} + # order: Polynomial / geometric order to use when constructing the mesh # std::optional # default: 3 @@ -183,6 +193,8 @@ optimization_methods-smoothstep:{} core_mapping:{} END BLOCK CONFIG)", format_opt(mesh.config.refinement_levels, d.refinement_levels.value()), + mesh.config.vacuum_refinement_levels.has_value() ? std::to_string(*mesh.config.vacuum_refinement_levels) : "inherit", + mesh.config.vacuum_outer_refinement_levels.has_value() ? std::to_string(*mesh.config.vacuum_outer_refinement_levels) : "inherit", format_opt(mesh.config.order, d.order.value()), format_opt(mesh.config.include_external_domain, d.include_external_domain.value()), format_opt(mesh.config.r_core, d.r_core.value()), @@ -395,12 +407,26 @@ END BLOCK CONFIG)", std::expected parse_config(const std::string& content) { config::MeshConfig cfg; + // Files written before core_mapping was introduced use the single core cube. + cfg.core_mapping = "spherified"; config::OptimizationMethods opt = cfg.optimization_methods.value_or(config::OptimizationMethods{}); using Handler = std::function(std::string_view)>; auto as_int = [](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_optional_int = [](std::optional* f) { + return [f](const std::string_view v) -> std::expected { + if (trim(v) == "inherit") { + f->reset(); + return {}; + } + auto r = parse_int(v); + if (!r) return std::unexpected(r.error()); + *f = *r; + return {}; + }; + }; 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 {}; }; }; @@ -408,6 +434,8 @@ END BLOCK CONFIG)", const std::unordered_map handlers = { {"refinement_levels", as_int(&cfg.refinement_levels)}, + {"vacuum_refinement_levels", as_optional_int(&cfg.vacuum_refinement_levels)}, + {"vacuum_outer_refinement_levels", as_optional_int(&cfg.vacuum_outer_refinement_levels)}, {"order", as_int(&cfg.order)}, {"include_external_domain", as_bool(&cfg.include_external_domain)}, {"r_core", as_double(&cfg.r_core)}, @@ -617,6 +645,37 @@ END BLOCK CONFIG)", return pm; } + bool HasHangingFaces(const mfem::Mesh& mesh) { + for (int face = 0; face < mesh.GetNumFaces(); ++face) { + if (mesh.GetFaceInformation(face).IsNonconformingCoarse()) return true; + } + return false; + } + + void RefineVisualizationMesh(mfem::Mesh& mesh, mfem::GridFunction* field = nullptr) { + while (true) { + std::vector marked(static_cast(mesh.GetNE()), false); + for (int face = 0; face < mesh.GetNumFaces(); ++face) { + const auto info = mesh.GetFaceInformation(face); + if (info.IsNonconformingCoarse()) { + marked[static_cast(info.element[0].index)] = true; + } + } + + mfem::Array refinements; + for (int element = 0; element < mesh.GetNE(); ++element) { + if (marked[static_cast(element)]) refinements.Append(element); + } + if (refinements.Size() == 0) break; + + mesh.GeneralRefinement(refinements, 1, 0); + if (field) { + field->FESpace()->Update(); + field->Update(); + } + } + } + } void SaveStroidMesh(const StroidMesh &mesh, const std::string &filename, const std::string &comment) { @@ -659,7 +718,13 @@ END BLOCK CONFIG)", SaveVTU(*mesh.mesh, exportName); } - void ViewMesh(mfem::Mesh &mesh, const std::string& title, const VISUALIZATION_MODE mode, const std::string &vishost, int visport) { + std::unique_ptr MakeConformingVisualizationMesh(const mfem::Mesh& mesh) { + auto display_mesh = std::make_unique(mesh); + RefineVisualizationMesh(*display_mesh); + return display_mesh; + } + + void ViewMesh(mfem::Mesh &mesh, const std::string& title, const VISUALIZATION_MODE mode, const std::string &vishost, int visport, const bool conforming_display) { mfem::socketstream sol_sock(vishost.c_str(), visport); if (!sol_sock.is_open()) { std::cerr << "Unable to connect to GLVis server at " @@ -667,8 +732,13 @@ END BLOCK CONFIG)", return; } - mfem::L2_FECollection fec(0, mesh.Dimension()); - mfem::FiniteElementSpace fes(&mesh, &fec); + std::unique_ptr display_mesh; + if (conforming_display && HasHangingFaces(mesh)) { + display_mesh = std::make_unique(mesh); + } + mfem::Mesh& viewed_mesh = display_mesh ? *display_mesh : mesh; + mfem::L2_FECollection fec(0, viewed_mesh.Dimension()); + mfem::FiniteElementSpace fes(&viewed_mesh, &fec); mfem::GridFunction attr_gf(&fes); attr_gf = 0.0; @@ -692,44 +762,55 @@ END BLOCK CONFIG)", break; } - sol_sock.precision(8); - sol_sock << "solution\n" << mesh << attr_gf; - sol_sock << "window_title '" << title << "'\n"; + // Transfer source coloring so boundary-adjacent regions keep their + // original extent when visualization-only children are introduced. + if (display_mesh) RefineVisualizationMesh(*display_mesh, &attr_gf); + + sol_sock.precision(std::numeric_limits::max_digits10); + sol_sock << "solution\n" << viewed_mesh << attr_gf; + sol_sock << "window_title '" << title + << (display_mesh ? " (display subdivisions)" : "") << "'\n"; sol_sock << "keys iMj\n"; sol_sock << std::flush; } - void ViewMesh(const stroid::StroidMesh &mesh, const std::string &title, VISUALIZATION_MODE mode, const std::string &vishost, int visport) { - ViewMesh(*mesh.mesh, title, mode, vishost, visport); + void ViewMesh(const stroid::StroidMesh &mesh, const std::string &title, VISUALIZATION_MODE mode, const std::string &vishost, int visport, const bool conforming_display) { + ViewMesh(*mesh.mesh, title, mode, vishost, visport, conforming_display); } - void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport) { - mfem::L2_FECollection fec(0, 3); - mfem::FiniteElementSpace fes(&mesh, &fec); + void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport, const bool conforming_display) { + mfem::socketstream sol_sock(vishost.c_str(), visport); + if (!sol_sock.is_open()) return; + + std::unique_ptr display_mesh; + if (conforming_display && HasHangingFaces(mesh)) { + display_mesh = std::make_unique(mesh); + } + mfem::Mesh& viewed_mesh = display_mesh ? *display_mesh : mesh; + mfem::L2_FECollection fec(0, viewed_mesh.Dimension()); + mfem::FiniteElementSpace fes(&viewed_mesh, &fec); mfem::GridFunction valence_gf(&fes); + valence_gf = 0.0; for (int i = 0; i < mesh.GetNBE(); i++) { - int f, o; - mesh.GetBdrElementFace(i, &f, &o); - - int e1, e2; - mesh.GetFaceElements(f, &e1, &e2); - - int valence = (e2 >= 0) ? 2 : 1; - valence_gf(i) = static_cast(valence); + const int face = mesh.GetBdrElementFaceIndex(i); + const double valence = mesh.GetFaceInformation(face).IsInterior() ? 2.0 : 1.0; + int element, side; + mesh.GetBdrElementAdjacentElement(i, element, side); + valence_gf(element) = std::max(valence_gf(element), valence); } - // View in GLVis - mfem::socketstream sol_sock(vishost.c_str(), visport); - if (sol_sock.is_open()) { - sol_sock << "solution\n" << mesh << valence_gf; - sol_sock << "window_title 'Boundary Valence: 1=Surface, 2=Internal'\n"; - sol_sock << "keys am\n" << std::flush; - } + if (display_mesh) RefineVisualizationMesh(*display_mesh, &valence_gf); + + sol_sock.precision(std::numeric_limits::max_digits10); + sol_sock << "solution\n" << viewed_mesh << valence_gf; + sol_sock << "window_title 'Boundary Valence: 1=Surface, 2=Internal" + << (display_mesh ? " (display subdivisions)" : "") << "'\n"; + sol_sock << "keys am\n" << std::flush; } - void VisualizeFaceValence(const stroid::StroidMesh &mesh, const std::string &vishost, int visport) { - VisualizeFaceValence(*mesh.mesh, vishost, visport); + void VisualizeFaceValence(const stroid::StroidMesh &mesh, const std::string &vishost, int visport, const bool conforming_display) { + VisualizeFaceValence(*mesh.mesh, vishost, visport, conforming_display); } std::expected ParseStroidMesh(std::istream& is) { diff --git a/src/lib/refinement/uniform.cpp b/src/lib/refinement/uniform.cpp index 26b6d45..736c527 100644 --- a/src/lib/refinement/uniform.cpp +++ b/src/lib/refinement/uniform.cpp @@ -8,6 +8,8 @@ #include "stroid/topology/topology.h" #include "stroid/topology/optimize.h" +#include + namespace stroid::refinement { void UniformRefinement(StroidMesh &mesh, const size_t levels) { if (!mesh.reference_mesh) { @@ -19,14 +21,21 @@ namespace stroid::refinement { } 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"); + 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(); + if (levels > std::numeric_limits::max() - mesh.refinement_levels) { + throw std::overflow_error("Uniform refinement level count would overflow."); } - mesh.refinement_levels += levels; + StroidMesh refined; + refined.type = mesh.type; + refined.config = mesh.config; + refined.refinement_levels = mesh.refinement_levels + levels; + refined.reference_mesh = std::make_unique(*mesh.reference_mesh); + for (size_t i = 0; i < levels; i++) { + refined.reference_mesh->UniformRefinement(); + } fourdst::config::Config cfg; auto Mutator = [&mesh](config::MeshConfig& orig) { @@ -35,8 +44,13 @@ 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); + refined.mesh = utils::BuildProjected(*refined.reference_mesh, cfg); + topology::OptimizeMesh(*refined.mesh, cfg); + refined.exterior_coordinate = topology::BuildExteriorCoordinate(*refined.mesh, *refined.reference_mesh, cfg); + + mesh.mesh.swap(refined.mesh); + mesh.reference_mesh.swap(refined.reference_mesh); + mesh.exterior_coordinate.swap(refined.exterior_coordinate); + mesh.refinement_levels = refined.refinement_levels; } -} \ No newline at end of file +} diff --git a/src/lib/topology/curvilinear.cpp b/src/lib/topology/curvilinear.cpp index 65c55c2..e954aa2 100644 --- a/src/lib/topology/curvilinear.cpp +++ b/src/lib/topology/curvilinear.cpp @@ -26,7 +26,8 @@ namespace { } 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; + const double tolerance = std::min(1.0e-8, 1024.0 * std::numeric_limits::epsilon() * + std::max(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)); @@ -40,9 +41,7 @@ namespace { 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()); - auto* fes = new mfem::FiniteElementSpace(&mesh, fec, mesh.SpaceDimension()); - mesh.SetNodalFESpace(fes); + mesh.SetCurvature(config->order.value(), false, mesh.SpaceDimension(), mfem::Ordering::byNODES); } void ProjectMesh(mfem::Mesh &mesh, const fourdst::config::Config &config) { @@ -59,16 +58,15 @@ namespace stroid::topology { const int nElem = mesh.GetNE(); std::vector processed(nDofs, false); - mfem::Array vdofs; + mfem::Array dofs; mfem::Vector pos(vDim); for (int elemID = 0; elemID < nElem; ++elemID) { const int attrID = mesh.GetAttribute(elemID); - fes->GetElementVDofs(elemID, vdofs); + fes->GetElementDofs(elemID, dofs); - for (int dofID = 0; dofID < vdofs.Size(); ++dofID) { - const int vDof = vdofs[dofID]; - const int scalar_dof = (fes->GetOrdering() == mfem::Ordering::byNODES) ? vDof / vDim : vDof % nDofs; + for (int dofID = 0; dofID < dofs.Size(); ++dofID) { + const int scalar_dof = dofs[dofID] >= 0 ? dofs[dofID] : -1 - dofs[dofID]; if (processed[scalar_dof]) { continue; // Skip already processed dofs. This avoids doing multiple transformations of a node if it was already transformed by a neighbor @@ -88,6 +86,12 @@ namespace stroid::topology { } } + // A mapped hanging node must lie on the coarse polynomial face. Mapping + // every node independently does not preserve this geometric constraint. + mfem::Vector true_nodes; + nodes.GetTrueDofs(true_nodes); + nodes.SetFromTrueDofs(true_nodes); + mesh.NodesUpdated(); } std::unique_ptr BuildExteriorCoordinate( @@ -164,6 +168,10 @@ namespace stroid::topology { } } + mfem::Vector true_values; + field->values->GetTrueDofs(true_values); + field->values->SetFromTrueDofs(true_values); + 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); diff --git a/src/lib/topology/optimize.cpp b/src/lib/topology/optimize.cpp index 648376b..c28562c 100644 --- a/src/lib/topology/optimize.cpp +++ b/src/lib/topology/optimize.cpp @@ -181,7 +181,8 @@ class TMOPProgressBar : public mfem::IterativeSolverMonitor { a.SetEssentialTrueDofs(ess_tdof_list); mfem::GridFunction* nodes = mesh.GetNodes(); - mfem::Vector x(*nodes); + mfem::Vector x; + nodes->GetTrueDofs(x); mfem::Vector b(a.Height()); b = 0.0; @@ -222,7 +223,7 @@ class TMOPProgressBar : public mfem::IterativeSolverMonitor { std::cout << "Applying TMOP optimization to mesh. Note this may take a long time. Depending on your mesh resolution expect to wait up to the order of 10s of minutes..." << std::endl; newton.Mult(b, x); - *nodes = x; + nodes->SetFromTrueDofs(x); mesh.NodesUpdated(); diff --git a/src/lib/topology/topology.cpp b/src/lib/topology/topology.cpp index 63cbaa8..9132fcc 100644 --- a/src/lib/topology/topology.cpp +++ b/src/lib/topology/topology.cpp @@ -3,13 +3,116 @@ #include #include #include +#include +#include #include "stroid/config/config.h" #include "fourdst/config/config.h" +namespace { + void ValidateRefinement(const stroid::config::MeshConfig& config) { + if (config.refinement_levels.value_or(4) < 0 || + config.vacuum_refinement_levels.value_or(0) < 0 || + config.vacuum_outer_refinement_levels.value_or(0) < 0) { + throw std::invalid_argument("Refinement levels must be non-negative."); + } + if (config.order.value_or(3) < 1) { + throw std::invalid_argument("Geometry order must be at least one."); + } + + if (!config.vacuum_refinement_levels && !config.vacuum_outer_refinement_levels) return; + if (!config.include_external_domain.value_or(true)) { + throw std::invalid_argument("Vacuum refinement overrides require an external domain."); + } + + const double r_core = config.r_core.value_or(0.25); + const double r_star = config.r_star.value_or(1.0); + const double r_infinity = config.r_infinity.value_or(6.0); + if (!std::isfinite(r_core) || !std::isfinite(r_star) || !std::isfinite(r_infinity) || + r_core <= 0.0 || r_star <= r_core || r_infinity <= r_star) { + throw std::invalid_argument("Vacuum refinement requires finite radii with 0 < r_core < r_star < r_infinity."); + } + if (!std::isfinite(config.flattening.value_or(0.0)) || config.flattening.value_or(0.0) >= 1.0) { + throw std::invalid_argument("Vacuum refinement requires finite flattening < 1."); + } + + const auto core = config.core_id.value_or(1); + const auto envelope = config.envelope_id.value_or(2); + const auto vacuum = config.vacuum_id.value_or(3); + const auto surface = config.surface_bdr_id.value_or(1); + const auto outer = config.inf_bdr_id.value_or(2); + const auto valid_id = [](size_t id) { + return id > 0 && id <= static_cast(std::numeric_limits::max()); + }; + if (!valid_id(core) || !valid_id(envelope) || !valid_id(vacuum) || + !valid_id(surface) || !valid_id(outer) || core == envelope || + core == vacuum || envelope == vacuum || surface == outer) { + throw std::invalid_argument("Vacuum refinement requires distinct positive material IDs and distinct positive boundary IDs representable as int."); + } + } + + void RefineReference(mfem::Mesh& mesh, const stroid::config::MeshConfig& config) { + const int stellar_level = config.refinement_levels.value_or(4); + const int bulk_level = config.vacuum_refinement_levels.value_or(stellar_level); + const int outer_level = config.vacuum_outer_refinement_levels.value_or(stellar_level); + if (!config.include_external_domain.value_or(true) || + (bulk_level == stellar_level && outer_level == stellar_level)) { + for (int level = 0; level < stellar_level; ++level) mesh.UniformRefinement(); + return; + } + + mesh.EnsureNCMesh(); + const int vacuum_id = static_cast(config.vacuum_id.value_or(3)); + const double r_star = config.r_star.value_or(1.0); + const double r_infinity = config.r_infinity.value_or(6.0); + const double tolerance = 128.0 * std::numeric_limits::epsilon() * r_infinity; + + while (true) { + std::vector marked(static_cast(mesh.GetNE()), false); + for (int element = 0; element < mesh.GetNE(); ++element) { + int target = stellar_level; + if (mesh.GetAttribute(element) == vacuum_id) { + target = bulk_level; + double minimum = std::numeric_limits::infinity(); + double maximum = 0.0; + const mfem::Element* hex = mesh.GetElement(element); + for (int vertex = 0; vertex < hex->GetNVertices(); ++vertex) { + const double* position = mesh.GetVertex(hex->GetVertices()[vertex]); + const double radius = std::max({std::abs(position[0]), std::abs(position[1]), std::abs(position[2])}); + minimum = std::min(minimum, radius); + maximum = std::max(maximum, radius); + } + if (std::abs(minimum - r_star) <= tolerance) target = std::max(target, stellar_level); + if (std::abs(maximum - r_infinity) <= tolerance) target = std::max(target, outer_level); + } + marked[static_cast(element)] = mesh.ncmesh->GetElementDepth(element) < target; + } + + for (int face = 0; face < mesh.GetNumFaces(); ++face) { + const auto info = mesh.GetFaceInformation(face); + if (!info.IsNonconformingFine() || !info.IsLocal()) continue; + const int first = info.element[0].index; + const int second = info.element[1].index; + if ((mesh.GetAttribute(first) == vacuum_id) == (mesh.GetAttribute(second) == vacuum_id)) continue; + const int coarse = mesh.ncmesh->GetElementDepth(first) < mesh.ncmesh->GetElementDepth(second) ? first : second; + marked[static_cast(coarse)] = true; + } + + mfem::Array refinements; + for (int element = 0; element < mesh.GetNE(); ++element) { + if (marked[static_cast(element)]) refinements.Append(element); + } + if (refinements.Size() == 0) break; + mesh.GeneralRefinement(refinements, 1, 1); + } + mesh.CheckBdrElementOrientation(true); + } +} + namespace stroid::topology { std::unique_ptr BuildSkeleton(const fourdst::config::Config & config) { + ValidateRefinement(*config); 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); @@ -129,19 +232,12 @@ namespace stroid::topology { // ReSharper disable once CppUseInternalLinkage void Finalize(mfem::Mesh& mesh, const fourdst::config::Config &config) { + ValidateRefinement(*config); mesh.FinalizeTopology(); mesh.Finalize(); mesh.CheckElementOrientation(true); mesh.CheckBdrElementOrientation(true); - for (int i = 0; i < config->refinement_levels; ++i) { - mesh.UniformRefinement(); - } - - if (!mesh.Conforming()) { - std::cerr << "WARNING: Mesh has been detected to be non conforming!" << std::endl; - } - - + RefineReference(mesh, *config); } } diff --git a/src/lib/utils/mesh_stats.cpp b/src/lib/utils/mesh_stats.cpp index 3c90b0a..ed97102 100644 --- a/src/lib/utils/mesh_stats.cpp +++ b/src/lib/utils/mesh_stats.cpp @@ -61,11 +61,39 @@ namespace stroid::stats { if (has_feature(features, MeshStatFeatures::CONFORMITY)) { ConformityStats conformity; - conformity.conforming = mesh->Conforming(); - conformity.n_nonconforming_faces = conformity.conforming ? 0 : -99; // TODO: count + conformity.hierarchy_enabled = mesh->ncmesh != nullptr; + for (int f = 0; f < mesh->GetNFaces(); ++f) { + if (mesh->GetFaceInformation(f).IsNonconformingFine()) { + ++conformity.n_nonconforming_faces; + } + } + conformity.conforming = conformity.n_nonconforming_faces == 0; out.conformity = conformity; } + if (has_feature(features, MeshStatFeatures::REFINEMENT)) { + RefinementStats refinement; + auto update_depth = [](RegionRefinementStats& region, const int depth) { + if (region.min_depth < 0) region.min_depth = depth; + region.min_depth = std::min(region.min_depth, depth); + region.max_depth = std::max(region.max_depth, depth); + }; + for (int e = 0; e < mesh->GetNE(); ++e) { + const int depth = mesh->ncmesh ? mesh->ncmesh->GetElementDepth(e) : + static_cast(sm.refinement_levels); + update_depth(refinement.all, depth); + const int attr = mesh->GetAttribute(e); + if (attr == core_id) update_depth(refinement.core, depth); + else if (attr == env_id) update_depth(refinement.envelope, depth); + else if (attr == vac_id) update_depth(refinement.vacuum, depth); + } + if (const auto* fes = mesh->GetNodalFESpace()) { + refinement.geometry_dofs = fes->GetNDofs(); + refinement.geometry_true_dofs = fes->GetNConformingDofs(); + } + out.refinement = refinement; + } + // ============================ SURFACE PASS ============================ const bool needs_surface = has_feature(features, MeshStatFeatures::RADIUS) || @@ -428,7 +456,18 @@ namespace stroid::stats { b.max_inward, b.max_outward, b.rms)); } if (s.conformity) { - line(std::format("conforming: {}", s.conformity->conforming)); + const auto& c = *s.conformity; + line(std::format("conforming: {} (hierarchy={}, hanging face patches={})", + c.conforming, c.hierarchy_enabled, c.n_nonconforming_faces)); + } + if (s.refinement) { + const auto& r = *s.refinement; + line(std::format( + "refinement depth: all=[{},{}] core=[{},{}] env=[{},{}] vac=[{},{}]", + r.all.min_depth, r.all.max_depth, r.core.min_depth, r.core.max_depth, + r.envelope.min_depth, r.envelope.max_depth, r.vacuum.min_depth, r.vacuum.max_depth)); + line(std::format("scalar geometry DOFs: total={} true={}", + r.geometry_dofs, r.geometry_true_dofs)); } auto jac_line = [&](const std::string& label, const JacobianStats& j) { @@ -482,4 +521,4 @@ namespace stroid::stats { for (const auto& e : s.errors) line(std::format("ERROR: {}", e)); return o; } -} \ No newline at end of file +} diff --git a/src/python/IO/bindings.cpp b/src/python/IO/bindings.cpp index 11108bb..b31477c 100644 --- a/src/python/IO/bindings.cpp +++ b/src/python/IO/bindings.cpp @@ -35,20 +35,29 @@ void register_io_bindings(pybind11::module_ &m) { ); m.def( "ViewMesh", - py::overload_cast(&stroid::IO::ViewMesh), + py::overload_cast(&stroid::IO::ViewMesh), py::arg("mesh"), py::arg("title")="", py::arg("mode")=stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, py::arg("host")="localhost", - py::arg("port")=19916 + py::arg("port")=19916, + py::arg("conforming_display")=false, + "Display the mesh in GLVis. By default, subdivide a temporary copy to avoid " + "rendering gaps at curved hanging interfaces, preserving the source geometry " + "and coloring. Extra display edges do not change computational DOFs. Set " + "conforming_display=False to inspect the original element layout." ); m.def( "VisualizeFaceValence", - py::overload_cast(&stroid::IO::VisualizeFaceValence), + py::overload_cast(&stroid::IO::VisualizeFaceValence), py::arg("mesh"), py::arg("host")="localhost", - py::arg("port")=19916 + py::arg("port")=19916, + py::arg("conforming_display")=true, + "Display boundary-adjacent element valence: zero for untagged elements, one " + "for surface faces, and two for internal faces (maximum if several touch an " + "element). Values are preserved through optional display-only subdivision." ); m.def( diff --git a/src/python/config/bindings.cpp b/src/python/config/bindings.cpp index 8b1ba23..6b5b3da 100644 --- a/src/python/config/bindings.cpp +++ b/src/python/config/bindings.cpp @@ -38,6 +38,8 @@ void register_config_bindings(pybind11::module_& m) { return stroid::config::MeshConfig{ .refinement_levels = kwargs.contains("refinement_levels") ? kwargs["refinement_levels"].cast() : ref_level, + .vacuum_refinement_levels = kwargs.contains("vacuum_refinement_levels") ? kwargs["vacuum_refinement_levels"].cast>() : std::nullopt, + .vacuum_outer_refinement_levels = kwargs.contains("vacuum_outer_refinement_levels") ? kwargs["vacuum_outer_refinement_levels"].cast>() : std::nullopt, .order = kwargs.contains("order") ? kwargs["order"].cast() : order, .include_external_domain = kwargs.contains("include_external_domain") ? kwargs["include_external_domain"].cast() : include_external_domain, .r_core = kwargs.contains("r_core") ? kwargs["r_core"].cast() : r_core, @@ -53,7 +55,7 @@ void register_config_bindings(pybind11::module_& m) { .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, - .core_mapping = kwargs.contains("core_mapping") ? kwargs["core_mapping"].cast() : "spherified" + .core_mapping = kwargs.contains("core_mapping") ? kwargs["core_mapping"].cast() : "multi_block" }; })) .def_property( @@ -65,6 +67,26 @@ void register_config_bindings(pybind11::module_& m) { self.refinement_levels = value; } ) + .def_property( + "vacuum_refinement_levels", + [](const stroid::config::MeshConfig& self) { + return self.vacuum_refinement_levels; + }, + [](stroid::config::MeshConfig& self, std::optional value) { + self.vacuum_refinement_levels = value; + }, + "Minimum vacuum interior depth, or None to inherit refinement_levels. Automatic grading may refine further." + ) + .def_property( + "vacuum_outer_refinement_levels", + [](const stroid::config::MeshConfig& self) { + return self.vacuum_outer_refinement_levels; + }, + [](stroid::config::MeshConfig& self, std::optional value) { + self.vacuum_outer_refinement_levels = value; + }, + "Minimum vacuum outer-boundary depth, or None to inherit refinement_levels." + ) .def_property( "order", [](const stroid::config::MeshConfig& self) { @@ -88,8 +110,8 @@ void register_config_bindings(pybind11::module_& m) { [](const stroid::config::MeshConfig& self) { return self.r_core; }, - [](stroid::config::MeshConfig& self, int value) { - self.order = value; + [](stroid::config::MeshConfig& self, double value) { + self.r_core = value; } ) .def_property( diff --git a/src/python/refinement/bindings.cpp b/src/python/refinement/bindings.cpp index 72a0e48..aeac353 100644 --- a/src/python/refinement/bindings.cpp +++ b/src/python/refinement/bindings.cpp @@ -7,5 +7,5 @@ namespace py = pybind11; void register_refinement_bindings(pybind11::module_ &m) { - m.def("UniformRefinement", &stroid::refinement::UniformRefinement, py::arg("mesh"), py::arg("levels"), "Perform uniform refinement without breaking the higher order structure"); + m.def("UniformRefinement", &stroid::refinement::UniformRefinement, py::arg("mesh"), py::arg("levels"), "Refine every current leaf by the requested additional levels, preserving hanging-node constraints and rebuilding high-order geometry. Initial configuration targets remain unchanged."); } diff --git a/src/python/utils/bindings.cpp b/src/python/utils/bindings.cpp index af47d51..a5d5797 100644 --- a/src/python/utils/bindings.cpp +++ b/src/python/utils/bindings.cpp @@ -25,7 +25,14 @@ void register_stats_bindings(pybind11::module_ &m) { .value("CENTROID", stroid::stats::MeshStatFeatures::CENTROID) .value("CONFIG_META", stroid::stats::MeshStatFeatures::CONFIG_META) .value("BOUNDING_BOX", stroid::stats::MeshStatFeatures::BOUNDING_BOX) - .export_values(); + .value("REFINEMENT", stroid::stats::MeshStatFeatures::REFINEMENT) + .export_values() + .def("__or__", [](stroid::stats::MeshStatFeatures lhs, stroid::stats::MeshStatFeatures rhs) { + return lhs | rhs; + }, py::is_operator()) + .def("__and__", [](stroid::stats::MeshStatFeatures lhs, stroid::stats::MeshStatFeatures rhs) { + return lhs & rhs; + }, py::is_operator()); py::class_(statsMod, "RadiusStats") .def_readonly("min", &stroid::stats::RadiusStats::min) @@ -51,8 +58,21 @@ void register_stats_bindings(pybind11::module_ &m) { py::class_(statsMod, "ConformityStats") .def_readonly("conforming", &stroid::stats::ConformityStats::conforming) + .def_readonly("hierarchy_enabled", &stroid::stats::ConformityStats::hierarchy_enabled) .def_readonly("n_nonconforming_faces", &stroid::stats::ConformityStats::n_nonconforming_faces); + py::class_(statsMod, "RegionRefinementStats") + .def_readonly("min_depth", &stroid::stats::RegionRefinementStats::min_depth) + .def_readonly("max_depth", &stroid::stats::RegionRefinementStats::max_depth); + + py::class_(statsMod, "RefinementStats") + .def_readonly("all", &stroid::stats::RefinementStats::all) + .def_readonly("core", &stroid::stats::RefinementStats::core) + .def_readonly("envelope", &stroid::stats::RefinementStats::envelope) + .def_readonly("vacuum", &stroid::stats::RefinementStats::vacuum) + .def_readonly("geometry_dofs", &stroid::stats::RefinementStats::geometry_dofs) + .def_readonly("geometry_true_dofs", &stroid::stats::RefinementStats::geometry_true_dofs); + py::class_(statsMod, "JacobianStats") .def_readonly("detJ_min", &stroid::stats::JacobianStats::detJ_min) .def_readonly("detJ_max", &stroid::stats::JacobianStats::detJ_max) @@ -128,6 +148,7 @@ void register_stats_bindings(pybind11::module_ &m) { .def_readonly("ellipticity", &stroid::stats::MeshStats::ellipticity) .def_readonly("bowing", &stroid::stats::MeshStats::bowing) .def_readonly("conformity", &stroid::stats::MeshStats::conformity) + .def_readonly("refinement", &stroid::stats::MeshStats::refinement) .def_readonly("jacobian", &stroid::stats::MeshStats::jacobian) .def_readonly("jacobian_stellar", &stroid::stats::MeshStats::jacobian_stellar) .def_readonly("jacobian_vacuum", &stroid::stats::MeshStats::jacobian_vacuum) @@ -174,7 +195,13 @@ void register_type_bindings(py::module_ &m) { .def("has_rmesh", [](const stroid::StroidMesh& self) { return self.reference_mesh != nullptr; }) - .def("mesh_stats", &stroid::StroidMesh::mesh_stats) + .def("mesh_stats", [](const stroid::StroidMesh& self, bool use_ref_mesh) { + auto result = self.mesh_stats(use_ref_mesh); + if (!result.has_value()) { + throw std::runtime_error(result.error()); + } + return result.value(); + }, py::arg("use_ref_mesh") = false) .def("__repr__", [](const stroid::StroidMesh& self) { return std::format("", (self.type == stroid::MFEM_MESH_TYPE::SERIAL) ? "SERIAL" : "PARALLEL", self.mesh->GetNE(), self.mesh->GetNV()); }); @@ -184,4 +211,3 @@ void register_utils_bindings(pybind11::module_ &m) { register_type_bindings(m); register_stats_bindings(m); } - diff --git a/tests/meson.build b/tests/meson.build index 95ec39b..aaf4536 100644 --- a/tests/meson.build +++ b/tests/meson.build @@ -3,9 +3,10 @@ gtest_main = dependency('gtest_main', required: true) gtest_nomain_dep = dependency('gtest', main: false, required : true) threads_dep = dependency('threads') -# Test files for const test_sources = [ - 'stroidTest.cpp' + 'stroidTest.cpp', + 'nonconformingTest.cpp', + 'visualizationTest.cpp' ] foreach test_file : test_sources diff --git a/tests/nonconformingTest.cpp b/tests/nonconformingTest.cpp new file mode 100644 index 0000000..beadc4c --- /dev/null +++ b/tests/nonconformingTest.cpp @@ -0,0 +1,486 @@ +#include + +#include "stroid/stroid.h" + +#include +#include +#include +#include +#include +#include + +namespace { + constexpr double kPi = 3.14159265358979323846; + + stroid::config::MeshConfig Configuration(int order = 2, int stellar_level = 2, + int bulk_level = 0, int outer_level = 0) { + stroid::config::MeshConfig config; + config.order = order; + config.refinement_levels = stellar_level; + config.vacuum_refinement_levels = bulk_level; + config.vacuum_outer_refinement_levels = outer_level; + config.optimization_methods = stroid::config::OptimizationMethods{false, true}; + return config; + } + + std::map CountAttributes(const mfem::Mesh& mesh, bool boundary = false) { + std::map counts; + for (int element = 0; element < (boundary ? mesh.GetNBE() : mesh.GetNE()); ++element) { + ++counts[boundary ? mesh.GetBdrAttribute(element) : mesh.GetAttribute(element)]; + } + return counts; + } + + void ExpectConstrainedField(mfem::GridFunction& values) { + mfem::Vector independent; + values.GetTrueDofs(independent); + mfem::GridFunction reconstructed(values.FESpace()); + reconstructed.SetFromTrueDofs(independent); + reconstructed -= values; + EXPECT_LT(reconstructed.Normlinf(), 5.0e-12); + } + + void ExpectGeometryAndCoordinateTraces(stroid::StroidMesh& generated, bool require_hanging_faces = true) { + mfem::Mesh& mesh = *generated.mesh; + ASSERT_NE(generated.exterior_coordinate, nullptr); + mfem::GridFunction& coordinate = *generated.exterior_coordinate->values; + ASSERT_EQ(coordinate.FESpace()->GetMesh(), &mesh); + ExpectConstrainedField(*mesh.GetNodes()); + ExpectConstrainedField(coordinate); + const int vacuum = static_cast(generated.config.vacuum_id.value()); + int hanging_faces = 0; + int stellar_faces = 0; + double geometry_error = 0.0; + double coordinate_error = 0.0; + double stellar_trace_error = 0.0; + mfem::Vector first(3), second(3); + + for (int face = 0; face < mesh.GetNumFaces(); ++face) { + const auto info = mesh.GetFaceInformation(face); + if (!info.IsLocal()) continue; + auto* transformation = mesh.GetFaceElementTransformations(face); + ASSERT_NE(transformation->Elem1, nullptr); + ASSERT_NE(transformation->Elem2, nullptr); + const bool first_vacuum = transformation->Elem1->Attribute == vacuum; + const bool second_vacuum = transformation->Elem2->Attribute == vacuum; + const bool stellar_interface = first_vacuum != second_vacuum; + if (info.IsNonconformingFine()) { + ++hanging_faces; + EXPECT_EQ(first_vacuum, second_vacuum); + } + if (stellar_interface) { + ++stellar_faces; + EXPECT_TRUE(info.IsConforming()) << "Stellar interface face " << face; + } + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + mfem::IntegrationPoint point; + point.Set2(i / 3.0, j / 3.0); + transformation->SetAllIntPoints(&point); + const auto& first_point = transformation->Elem1->GetIntPoint(); + const auto& second_point = transformation->Elem2->GetIntPoint(); + transformation->Elem1->Transform(first_point, first); + transformation->Elem2->Transform(second_point, second); + first -= second; + geometry_error = std::max(geometry_error, first.Norml2()); + const double first_value = coordinate.GetValue(transformation->Elem1No, first_point); + const double second_value = coordinate.GetValue(transformation->Elem2No, second_point); + coordinate_error = std::max(coordinate_error, std::abs(first_value - second_value)); + if (stellar_interface) { + stellar_trace_error = std::max({stellar_trace_error, std::abs(first_value), std::abs(second_value)}); + } + } + } + } + if (require_hanging_faces) EXPECT_GT(hanging_faces, 0); + EXPECT_GT(stellar_faces, 0); + EXPECT_LT(geometry_error, 5.0e-12); + EXPECT_LT(coordinate_error, 5.0e-12); + EXPECT_LT(stellar_trace_error, 5.0e-12); + + int outer_faces = 0; + double outer_trace_error = 0.0; + for (int boundary = 0; boundary < mesh.GetNBE(); ++boundary) { + if (mesh.GetBdrAttribute(boundary) != static_cast(generated.config.inf_bdr_id.value())) continue; + ++outer_faces; + auto* transformation = mesh.GetBdrFaceTransformations(boundary); + const auto& quadrature = mfem::IntRules.Get(mfem::Geometry::SQUARE, 6); + for (int q = 0; q < quadrature.GetNPoints(); ++q) { + transformation->SetAllIntPoints(&quadrature.IntPoint(q)); + outer_trace_error = std::max(outer_trace_error, + std::abs(coordinate.GetValue(transformation->Elem1No, transformation->Elem1->GetIntPoint()) - 1.0)); + } + } + EXPECT_GT(outer_faces, 0); + EXPECT_LT(outer_trace_error, 5.0e-12); + + double minimum = 1.0; + double maximum = 0.0; + double interior_error = 0.0; + for (int element = 0; element < mesh.GetNE(); ++element) { + const auto& quadrature = mfem::IntRules.Get(mfem::Geometry::CUBE, 6); + for (int q = 0; q < quadrature.GetNPoints(); ++q) { + const double value = coordinate.GetValue(element, quadrature.IntPoint(q)); + ASSERT_TRUE(std::isfinite(value)); + if (mesh.GetAttribute(element) == vacuum) { + minimum = std::min(minimum, value); + maximum = std::max(maximum, value); + } else { + interior_error = std::max(interior_error, std::abs(value)); + } + } + } + EXPECT_GE(minimum, -5.0e-12); + EXPECT_LE(maximum, 1.0 + 5.0e-12); + EXPECT_LT(interior_error, 5.0e-12); + } + + void ExpectPositiveJacobians(mfem::Mesh& mesh, int excluded_attribute = -1) { + double minimum = std::numeric_limits::infinity(); + int minimum_element = -1; + for (int element = 0; element < mesh.GetNE(); ++element) { + if (mesh.GetAttribute(element) == excluded_attribute) continue; + auto* transformation = mesh.GetElementTransformation(element); + auto inspect = [&](const mfem::IntegrationPoint& point) { + transformation->SetIntPoint(&point); + const double determinant = transformation->Jacobian().Det(); + ASSERT_TRUE(std::isfinite(determinant)); + if (determinant < minimum) { + minimum = determinant; + minimum_element = element; + } + }; + const auto& quadrature = mfem::IntRules.Get(mfem::Geometry::CUBE, 2 * transformation->Order() + 2); + for (int q = 0; q < quadrature.GetNPoints(); ++q) inspect(quadrature.IntPoint(q)); + for (int i = 0; i <= 2; ++i) { + for (int j = 0; j <= 2; ++j) { + for (int k = 0; k <= 2; ++k) { + mfem::IntegrationPoint point; + point.Set3(i / 2.0, j / 2.0, k / 2.0); + inspect(point); + } + } + } + } + EXPECT_GT(minimum, 0.0) << "Element " << minimum_element; + } + + double StellarVolume(stroid::StroidMesh& generated) { + double volume = 0.0; + for (int element = 0; element < generated.mesh->GetNE(); ++element) { + if (generated.mesh->GetAttribute(element) == static_cast(generated.config.vacuum_id.value())) continue; + auto* transformation = generated.mesh->GetElementTransformation(element); + const auto& quadrature = mfem::IntRules.Get(mfem::Geometry::CUBE, 3 * transformation->Order() + 3); + for (int q = 0; q < quadrature.GetNPoints(); ++q) { + const auto& point = quadrature.IntPoint(q); + transformation->SetIntPoint(&point); + volume += point.weight * transformation->Jacobian().Det(); + } + } + return volume; + } + + double SurfaceRadiusError(stroid::StroidMesh& generated) { + double error = 0.0; + mfem::Vector physical(3); + for (int boundary = 0; boundary < generated.mesh->GetNBE(); ++boundary) { + if (generated.mesh->GetBdrAttribute(boundary) != static_cast(generated.config.surface_bdr_id.value())) continue; + auto* transformation = generated.mesh->GetBdrElementTransformation(boundary); + const auto& quadrature = mfem::IntRules.Get(mfem::Geometry::SQUARE, 8); + for (int q = 0; q < quadrature.GetNPoints(); ++q) { + transformation->Transform(quadrature.IntPoint(q), physical); + physical(2) /= 1.0 - generated.config.flattening.value(); + error = std::max(error, std::abs(physical.Norml2() - generated.config.r_star.value())); + } + } + return error; + } +} + +TEST(NonconformingMesh, UnspecifiedVacuumLevelsPreserveUniformGeneration) { + auto config = Configuration(2, 1); + config.vacuum_refinement_levels.reset(); + config.vacuum_outer_refinement_levels.reset(); + auto generated = stroid::GenerateMesh(config); + EXPECT_TRUE(generated.mesh->Conforming()); + EXPECT_EQ(generated.mesh->GetNE(), 19 * 8); + ExpectGeometryAndCoordinateTraces(generated, false); + + config.vacuum_refinement_levels = 1; + config.vacuum_outer_refinement_levels = 1; + auto explicit_levels = stroid::GenerateMesh(config); + EXPECT_EQ(explicit_levels.mesh->GetNE(), generated.mesh->GetNE()); + mfem::H1_FECollection collection(2, 3); + mfem::FiniteElementSpace space(explicit_levels.mesh.get(), &collection); + EXPECT_EQ(space.GetTrueVSize(), space.GetVSize()); + EXPECT_NEAR(StellarVolume(explicit_levels), StellarVolume(generated), 1.0e-12); +} + +TEST(NonconformingMesh, RejectsNegativeRefinementTargets) { + for (int field = 0; field < 3; ++field) { + auto config = Configuration(); + if (field == 0) config.refinement_levels = -1; + if (field == 1) config.vacuum_refinement_levels = -1; + if (field == 2) config.vacuum_outer_refinement_levels = -1; + EXPECT_THROW(stroid::GenerateMesh(config), std::invalid_argument); + } +} + +TEST(NonconformingMesh, DefaultOuterLevelProtectsBothSurfacesAndSavesBulkDofs) { + auto config = Configuration(2, 3); + config.vacuum_outer_refinement_levels.reset(); + auto generated = stroid::GenerateMesh(config); + ASSERT_NE(generated.reference_mesh->ncmesh, nullptr); + const auto attributes = CountAttributes(*generated.mesh); + EXPECT_EQ(attributes.at(1), 7 * 512); + EXPECT_EQ(attributes.at(2), 6 * 512); + EXPECT_LT(attributes.at(3), 6 * 512); + const auto boundaries = CountAttributes(*generated.mesh, true); + EXPECT_EQ(boundaries.at(1), 6 * 64); + EXPECT_EQ(boundaries.at(2), 6 * 64); + int vacuum_minimum = 3; + int vacuum_maximum = 0; + for (int element = 0; element < generated.mesh->GetNE(); ++element) { + ASSERT_EQ(generated.mesh->GetAttribute(element), generated.reference_mesh->GetAttribute(element)); + if (generated.mesh->GetAttribute(element) == 3) { + const int depth = generated.reference_mesh->ncmesh->GetElementDepth(element); + vacuum_minimum = std::min(vacuum_minimum, depth); + vacuum_maximum = std::max(vacuum_maximum, depth); + } + } + EXPECT_LT(vacuum_minimum, 3); + EXPECT_EQ(vacuum_maximum, 3); + for (int face = 0; face < generated.reference_mesh->GetNumFaces(); ++face) { + const auto information = generated.reference_mesh->GetFaceInformation(face); + if (!information.IsLocal()) continue; + const int first = generated.reference_mesh->ncmesh->GetElementDepth(information.element[0].index); + const int second = generated.reference_mesh->ncmesh->GetElementDepth(information.element[1].index); + EXPECT_LE(std::abs(first - second), 1); + } + mfem::H1_FECollection collection(2, 3); + mfem::FiniteElementSpace reduced(generated.mesh.get(), &collection); + EXPECT_LT(reduced.GetTrueVSize(), reduced.GetVSize()); + const auto stats = stroid::stats::ComputeMeshStats(generated, + stroid::stats::MeshStatFeatures::REFINEMENT | stroid::stats::MeshStatFeatures::CONFORMITY | + stroid::stats::MeshStatFeatures::ELEMENT_COUNT); + ASSERT_TRUE(stats.refinement.has_value()); + ASSERT_TRUE(stats.conformity.has_value()); + ASSERT_TRUE(stats.element_counts.has_value()); + EXPECT_EQ(stats.refinement->vacuum.min_depth, vacuum_minimum); + EXPECT_EQ(stats.refinement->vacuum.max_depth, vacuum_maximum); + EXPECT_EQ(stats.refinement->core.min_depth, 3); + EXPECT_EQ(stats.refinement->core.max_depth, 3); + EXPECT_EQ(stats.refinement->geometry_dofs, reduced.GetVSize()); + EXPECT_EQ(stats.refinement->geometry_true_dofs, reduced.GetTrueVSize()); + EXPECT_TRUE(stats.conformity->hierarchy_enabled); + EXPECT_FALSE(stats.conformity->conforming); + EXPECT_GT(stats.conformity->n_nonconforming_faces, 0); + EXPECT_EQ(stats.element_counts->vacuum, attributes.at(3)); + config.vacuum_refinement_levels.reset(); + auto uniform = stroid::GenerateMesh(config); + mfem::FiniteElementSpace full(uniform.mesh.get(), &collection); + EXPECT_LT(reduced.GetTrueVSize(), full.GetTrueVSize()); + EXPECT_NEAR(StellarVolume(generated), StellarVolume(uniform), 2.0e-12); + EXPECT_NEAR(SurfaceRadiusError(generated), SurfaceRadiusError(uniform), 2.0e-13); + ExpectGeometryAndCoordinateTraces(generated); + ExpectPositiveJacobians(*generated.mesh); +} + +TEST(NonconformingMesh, OuterTargetCanExceedStellarTarget) { + auto config = Configuration(2, 1, 0, 3); + auto generated = stroid::GenerateMesh(config); + EXPECT_EQ(CountAttributes(*generated.mesh, true).at(2), 6 * 64); + ExpectGeometryAndCoordinateTraces(generated); + ExpectPositiveJacobians(*generated.mesh); +} + +TEST(NonconformingMesh, InterfaceClosureRaisesCoarseStellarBoundaryToMatchVacuum) { + auto generated = stroid::GenerateMesh(Configuration(2, 0, 0, 3)); + ASSERT_NE(generated.reference_mesh->ncmesh, nullptr); + EXPECT_EQ(CountAttributes(*generated.mesh, true).at(2), 6 * 64); + int envelope_maximum = 0; + for (int element = 0; element < generated.reference_mesh->GetNE(); ++element) { + if (generated.reference_mesh->GetAttribute(element) == static_cast(generated.config.envelope_id.value())) { + envelope_maximum = std::max(envelope_maximum, generated.reference_mesh->ncmesh->GetElementDepth(element)); + } + } + EXPECT_GT(envelope_maximum, 0); + for (int face = 0; face < generated.reference_mesh->GetNumFaces(); ++face) { + const auto information = generated.reference_mesh->GetFaceInformation(face); + if (!information.IsLocal()) continue; + const int first = generated.reference_mesh->ncmesh->GetElementDepth(information.element[0].index); + const int second = generated.reference_mesh->ncmesh->GetElementDepth(information.element[1].index); + EXPECT_LE(std::abs(first - second), 1); + } + ExpectGeometryAndCoordinateTraces(generated); + ExpectPositiveJacobians(*generated.mesh); +} + +TEST(NonconformingMesh, RefinementAndExteriorCoordinateAreInvariantUnderSmallLengthScales) { + auto config = Configuration(2, 1, 0, 3); + auto reference = stroid::GenerateMesh(config); + constexpr double scale = 1.0e-15; + config.r_core = config.r_core.value() * scale; + config.r_star = config.r_star.value() * scale; + config.r_infinity = config.r_infinity.value() * scale; + auto scaled = stroid::GenerateMesh(config); + ASSERT_EQ(scaled.mesh->GetNE(), reference.mesh->GetNE()); + ASSERT_EQ(scaled.mesh->GetNodes()->Size(), reference.mesh->GetNodes()->Size()); + double coordinate_error = 0.0; + for (int dof = 0; dof < scaled.mesh->GetNodes()->Size(); ++dof) { + coordinate_error = std::max(coordinate_error, + std::abs((*scaled.mesh->GetNodes())(dof) / scale - (*reference.mesh->GetNodes())(dof))); + } + EXPECT_LT(coordinate_error, 5.0e-12); + ASSERT_EQ(scaled.exterior_coordinate->values->Size(), reference.exterior_coordinate->values->Size()); + mfem::Vector difference(*scaled.exterior_coordinate->values); + difference -= *reference.exterior_coordinate->values; + EXPECT_LT(difference.Normlinf(), 5.0e-12); + ExpectGeometryAndCoordinateTraces(scaled); + ExpectPositiveJacobians(*scaled.mesh); +} + +TEST(NonconformingMesh, CurvedGeometryAndScalarConstraintsAcrossOrdersAndMappings) { + for (const std::string mapping : {"multi_block", "spherified"}) { + for (const int order : {1, 2, 3, 4}) { + SCOPED_TRACE(mapping + " order=" + std::to_string(order)); + auto config = Configuration(order); + config.core_mapping = mapping; + config.flattening = 0.2; + config.core_id = 11; + config.envelope_id = 17; + config.vacuum_id = 23; + config.surface_bdr_id = 31; + config.inf_bdr_id = 37; + auto generated = stroid::GenerateMesh(config); + EXPECT_EQ(CountAttributes(*generated.mesh).size(), 3); + EXPECT_EQ(CountAttributes(*generated.mesh, true).size(), 2); + ExpectGeometryAndCoordinateTraces(generated); + // The legacy spherified core has known corner degeneracies. Its + // envelope and vacuum must still remain strictly oriented. + ExpectPositiveJacobians(*generated.mesh, mapping == "spherified" ? 11 : -1); + } + } +} + +TEST(NonconformingMesh, NoVacuumLeavesOnlyTheUniformStellarMesh) { + auto config = Configuration(2, 1, 0, 3); + config.include_external_domain = false; + EXPECT_THROW(stroid::GenerateMesh(config), std::invalid_argument); + config.vacuum_refinement_levels.reset(); + config.vacuum_outer_refinement_levels.reset(); + auto generated = stroid::GenerateMesh(config); + EXPECT_EQ(generated.mesh->GetNE(), 13 * 8); + EXPECT_EQ(generated.exterior_coordinate, nullptr); + EXPECT_EQ(CountAttributes(*generated.mesh).size(), 2); + ExpectPositiveJacobians(*generated.mesh); +} + +TEST(NonconformingMesh, SerializationAndSubsequentRefinementPreserveHierarchy) { + auto generated = stroid::GenerateMesh(Configuration()); + const auto path = std::filesystem::temp_directory_path() / "stroid_nonconforming_round_trip.smesh"; + stroid::IO::SaveStroidMesh(generated, path.string(), "Nonconforming hierarchy regression"); + auto result = stroid::IO::LoadStroidMesh(path.string()); + ASSERT_TRUE(result.has_value()) << result.error(); + auto loaded = std::move(*result); + ASSERT_NE(loaded.reference_mesh->ncmesh, nullptr); + ASSERT_NE(loaded.mesh->ncmesh, nullptr); + EXPECT_EQ(loaded.config.vacuum_refinement_levels, generated.config.vacuum_refinement_levels); + EXPECT_EQ(loaded.config.vacuum_outer_refinement_levels, generated.config.vacuum_outer_refinement_levels); + ASSERT_EQ(loaded.mesh->GetNE(), generated.mesh->GetNE()); + for (int element = 0; element < loaded.mesh->GetNE(); ++element) { + EXPECT_EQ(loaded.reference_mesh->ncmesh->GetElementDepth(element), + generated.reference_mesh->ncmesh->GetElementDepth(element)); + } + EXPECT_NEAR(StellarVolume(loaded), StellarVolume(generated), 2.0e-12); + ExpectGeometryAndCoordinateTraces(loaded); + stroid::refinement::UniformRefinement(loaded, 1); + EXPECT_EQ(loaded.mesh->GetNE(), generated.mesh->GetNE() * 8); + EXPECT_EQ(loaded.mesh->GetNE(), loaded.reference_mesh->GetNE()); + ExpectGeometryAndCoordinateTraces(loaded); + ExpectPositiveJacobians(*loaded.mesh); + std::error_code error; + std::filesystem::remove(path, error); + EXPECT_FALSE(error); +} + +TEST(NonconformingMesh, LinearPhysicalPatchSolveUsesIndependentDofs) { + auto generated = stroid::GenerateMesh(Configuration()); + mfem::H1_FECollection collection(2, 3); + mfem::FiniteElementSpace space(generated.mesh.get(), &collection); + ASSERT_LT(space.GetTrueVSize(), space.GetVSize()); + mfem::FunctionCoefficient exact([](const mfem::Vector& point) { + return 1.0 + 0.3 * point(0) - 0.2 * point(1) + 0.1 * point(2); + }); + mfem::Array boundary(generated.mesh->bdr_attributes.Max()); + boundary = 0; + boundary[static_cast(generated.config.inf_bdr_id.value()) - 1] = 1; + mfem::Array essential; + space.GetEssentialTrueDofs(boundary, essential); + mfem::GridFunction solution(&space); + solution = 0.0; + solution.ProjectBdrCoefficient(exact, boundary); + mfem::LinearForm rhs(&space); + rhs = 0.0; + mfem::ConstantCoefficient one(1.0); + mfem::BilinearForm form(&space); + const auto& quadrature = mfem::IntRules.Get(mfem::Geometry::CUBE, 10); + auto* diffusion = new mfem::DiffusionIntegrator(one); + diffusion->SetIntRule(&quadrature); + form.AddDomainIntegrator(diffusion); + form.Assemble(); + mfem::OperatorPtr system; + mfem::Vector independent, forcing; + form.FormLinearSystem(essential, solution, rhs, system, independent, forcing); + EXPECT_EQ(system->Height(), space.GetTrueVSize()); + mfem::GSSmoother preconditioner(static_cast(*system)); + mfem::CGSolver solver; + solver.SetOperator(*system); + solver.SetPreconditioner(preconditioner); + solver.SetRelTol(1.0e-13); + solver.SetAbsTol(1.0e-14); + solver.SetMaxIter(1500); + solver.SetPrintLevel(-1); + solver.Mult(forcing, independent); + ASSERT_TRUE(solver.GetConverged()); + mfem::Vector residual(forcing.Size()); + system->Mult(independent, residual); + residual -= forcing; + EXPECT_LT(residual.Norml2() / forcing.Norml2(), 2.0e-12); + form.RecoverFEMSolution(independent, rhs, solution); + EXPECT_LT(solution.ComputeL2Error(exact), 2.0e-8); + ExpectConstrainedField(solution); +} + +TEST(NonconformingMesh, StellarVolumeAndSurfaceShapeConverge) { + auto coarse_config = Configuration(2, 2); + coarse_config.vacuum_outer_refinement_levels.reset(); + auto coarse = stroid::GenerateMesh(coarse_config); + auto fine_config = coarse_config; + fine_config.refinement_levels = 3; + auto fine = stroid::GenerateMesh(fine_config); + const double exact_volume = 4.0 * kPi / 3.0; + const double coarse_volume_error = std::abs(StellarVolume(coarse) - exact_volume); + const double fine_volume_error = std::abs(StellarVolume(fine) - exact_volume); + EXPECT_GT(coarse_volume_error, 1.0e-10); + EXPECT_LT(fine_volume_error, 0.5 * coarse_volume_error); + EXPECT_LT(SurfaceRadiusError(fine), 0.5 * SurfaceRadiusError(coarse)); +} + +TEST(NonconformingMesh, TMOPPreservesHangingConstraintsAndBoundaryTraces) { + auto config = Configuration(1); + auto initial = stroid::GenerateMesh(config); + config.optimization_methods = stroid::config::OptimizationMethods{true, true}; + auto generated = stroid::GenerateMesh(config); + mfem::Array marker(generated.mesh->bdr_attributes.Max()); + marker = 1; + mfem::Array essential; + generated.mesh->GetNodalFESpace()->GetEssentialTrueDofs(marker, essential); + mfem::Vector initial_nodes, optimized_nodes; + initial.mesh->GetNodes()->GetTrueDofs(initial_nodes); + generated.mesh->GetNodes()->GetTrueDofs(optimized_nodes); + ASSERT_EQ(initial_nodes.Size(), optimized_nodes.Size()); + for (const int dof : essential) EXPECT_NEAR(initial_nodes(dof), optimized_nodes(dof), 2.0e-13); + ExpectGeometryAndCoordinateTraces(generated); + ExpectPositiveJacobians(*generated.mesh); +} diff --git a/tests/stroidTest.cpp b/tests/stroidTest.cpp index c30b265..c79077f 100644 --- a/tests/stroidTest.cpp +++ b/tests/stroidTest.cpp @@ -358,11 +358,11 @@ std::optional EvalGridFunctionAtPoint( class stroidTest : public ::testing::Test {}; /** - * @brief Verifies the baseline block topology cardinalities in the no-vacuum case. + * @brief Verifies the default multi-block topology including the vacuum. * @details * Rationale: this is the fastest canary for accidental edits in block construction order, * vertex indexing, or boundary-face assembly. - * Method: build the default skeleton and assert exact counts (3D, 16 vertices, 7 hexes, 6 bdr quads). + * Method: build the default skeleton and assert exact counts (3D, 32 vertices, 19 hexes, 12 bdr quads). * If this fails: inspect `stroid::topology::BuildSkeleton` in `src/lib/topology/topology.cpp`, * especially `add_box`, `stellar_shells`, and `surface_bdr_quads`, plus ID defaults in * `src/include/stroid/config/config.h`. @@ -373,8 +373,8 @@ TEST_F(stroidTest, BuildSkeleton_DefaultCounts) { ASSERT_NE(mesh, nullptr); EXPECT_EQ(mesh->Dimension(), 3); - EXPECT_EQ(mesh->GetNV(), 24); - EXPECT_EQ(mesh->GetNE(), 13); + EXPECT_EQ(mesh->GetNV(), 32); + EXPECT_EQ(mesh->GetNE(), 19); EXPECT_EQ(mesh->GetNBE(), 12); } @@ -1355,8 +1355,8 @@ void ExpectCoreFaceContinuity(mfem::Mesh& mesh, int coreAttribute) { } // namespace -TEST_F(stroidTest, MultiBlockCore_TopologyCountsAndAttributesAreOptIn) { - EXPECT_EQ(stroid::config::MeshConfig{}.core_mapping.value(), "spherified"); +TEST_F(stroidTest, MultiBlockCore_DefaultTopologyCountsAndAttributes) { + EXPECT_EQ(stroid::config::MeshConfig{}.core_mapping.value(), "multi_block"); for (const bool external : {false, true}) { SCOPED_TRACE(external); auto cfg = MultiBlockConfiguration(2, 0, external); diff --git a/tests/visualizationTest.cpp b/tests/visualizationTest.cpp new file mode 100644 index 0000000..c758d5d --- /dev/null +++ b/tests/visualizationTest.cpp @@ -0,0 +1,182 @@ +#include + +#include "stroid/stroid.h" + +#include +#include +#include +#include + +namespace { + std::string Serialize(const mfem::Mesh& mesh) { + std::ostringstream stream; + stream.precision(std::numeric_limits::max_digits10); + mesh.Print(stream); + return stream.str(); + } + + int HangingFaces(const mfem::Mesh& mesh) { + int count = 0; + for (int face = 0; face < mesh.GetNumFaces(); ++face) { + count += mesh.GetFaceInformation(face).IsNonconformingFine(); + } + return count; + } + + void TessellatedPoint(mfem::ElementTransformation& transformation, + const mfem::IntegrationPoint& point, int subdivisions, + mfem::Vector& value) { + const double coordinates[] = {point.x, point.y, point.z}; + int cell[3]; + double local[3]; + for (int d = 0; d < 3; ++d) { + const double scaled = coordinates[d] * subdivisions; + cell[d] = std::clamp(static_cast(std::floor(scaled)), 0, subdivisions - 1); + local[d] = scaled - cell[d]; + } + value = 0.0; + mfem::Vector corner(3); + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < 2; ++j) { + for (int k = 0; k < 2; ++k) { + mfem::IntegrationPoint sample; + sample.Set3(static_cast(cell[0] + i) / subdivisions, + static_cast(cell[1] + j) / subdivisions, + static_cast(cell[2] + k) / subdivisions); + transformation.Transform(sample, corner); + const double weight = (i ? local[0] : 1.0 - local[0]) * + (j ? local[1] : 1.0 - local[1]) * + (k ? local[2] : 1.0 - local[2]); + value.Add(weight, corner); + } + } + } + } + + double FaceGap(mfem::Mesh& mesh, int subdivisions = 0) { + double maximum = 0.0; + mfem::Vector first(3), second(3); + const auto& quadrature = mfem::IntRules.Get(mfem::Geometry::SQUARE, 6); + for (int face = 0; face < mesh.GetNumFaces(); ++face) { + if (!mesh.GetFaceInformation(face).IsLocal()) continue; + auto* transformation = mesh.GetFaceElementTransformations(face); + for (int q = 0; q < quadrature.GetNPoints(); ++q) { + transformation->SetAllIntPoints(&quadrature.IntPoint(q)); + const auto first_point = transformation->Elem1->GetIntPoint(); + const auto second_point = transformation->Elem2->GetIntPoint(); + if (subdivisions > 0) { + TessellatedPoint(*transformation->Elem1, first_point, subdivisions, first); + TessellatedPoint(*transformation->Elem2, second_point, subdivisions, second); + } else { + transformation->Elem1->Transform(first_point, first); + transformation->Elem2->Transform(second_point, second); + } + first -= second; + maximum = std::max(maximum, first.Norml2()); + } + } + return maximum; + } + + stroid::StroidMesh CurvedVacuumMesh() { + stroid::config::MeshConfig config; + config.order = 3; + config.refinement_levels = 1; + config.vacuum_refinement_levels = 0; + config.vacuum_outer_refinement_levels = 2; + config.flattening = 0.15; + config.optimization_methods = stroid::config::OptimizationMethods{false, true}; + return stroid::GenerateMesh(config); + } +} + +TEST(Visualization, ConformingDisplayEliminatesCurvedTessellationGaps) { + auto generated = CurvedVacuumMesh(); + auto& source = *generated.mesh; + ASSERT_GT(HangingFaces(source), 0); + EXPECT_LT(FaceGap(source), 5.0e-12); + // Reproduce visible cracks even though the finite-element traces coincide. + EXPECT_GT(FaceGap(source, 2), 1.0e-4); + const auto original = Serialize(source); + const auto original_reference = Serialize(*generated.reference_mesh); + + auto display = stroid::IO::MakeConformingVisualizationMesh(source); + ASSERT_NE(display, nullptr); + EXPECT_GT(display->GetNE(), source.GetNE()); + EXPECT_EQ(HangingFaces(*display), 0); + EXPECT_EQ(display->GetNodalFESpace()->GetNDofs(), + display->GetNodalFESpace()->GetTrueVSize() / display->SpaceDimension()); + EXPECT_LT(FaceGap(*display), 5.0e-12); + for (int subdivisions : {1, 2, 3, 4}) { + EXPECT_LT(FaceGap(*display, subdivisions), 5.0e-12) << subdivisions; + } + + std::istringstream stream(Serialize(*display)); + mfem::Mesh reloaded(stream, 1, 1, true); + EXPECT_EQ(HangingFaces(reloaded), 0); + EXPECT_LT(FaceGap(reloaded, 2), 5.0e-12); + EXPECT_EQ(Serialize(source), original); + EXPECT_EQ(Serialize(*generated.reference_mesh), original_reference); +} + +TEST(Visualization, DisplayRefinementRestrictsExistingGeometryAndAttributes) { + auto generated = CurvedVacuumMesh(); + auto display = stroid::IO::MakeConformingVisualizationMesh(*generated.mesh); + auto reference = stroid::IO::MakeConformingVisualizationMesh(*generated.reference_mesh); + ASSERT_EQ(display->GetNE(), reference->GetNE()); + + mfem::DenseMatrix centers(3, reference->GetNE()); + mfem::Vector point(3), expected(3), actual(3); + for (int element = 0; element < reference->GetNE(); ++element) { + reference->GetElementCenter(element, point); + centers.SetCol(element, point); + } + mfem::Array parents; + mfem::Array parent_points; + ASSERT_EQ(generated.reference_mesh->FindPoints(centers, parents, parent_points, false), + reference->GetNE()); + double maximum_error = 0.0; + for (int element = 0; element < display->GetNE(); ++element) { + const int parent = parents[element]; + ASSERT_GE(parent, 0); + EXPECT_EQ(display->GetAttribute(element), generated.mesh->GetAttribute(parent)); + mfem::InverseElementTransformation inverse( + generated.reference_mesh->GetElementTransformation(parent)); + auto* logical = reference->GetElementTransformation(element); + auto* physical = display->GetElementTransformation(element); + auto* original = generated.mesh->GetElementTransformation(parent); + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + mfem::IntegrationPoint sample, parent_sample; + sample.Set3(i / 2.0, j / 2.0, k / 2.0); + logical->Transform(sample, point); + ASSERT_EQ(inverse.Transform(point, parent_sample), + mfem::InverseElementTransformation::Inside); + original->Transform(parent_sample, expected); + physical->Transform(sample, actual); + actual -= expected; + maximum_error = std::max(maximum_error, actual.Norml2()); + } + } + } + } + EXPECT_LT(maximum_error, 5.0e-12); +} + +TEST(Visualization, AlreadyMatchingFacesNeedNoAdditionalElements) { + stroid::config::MeshConfig config; + config.order = 3; + config.refinement_levels = 1; + config.optimization_methods = stroid::config::OptimizationMethods{false, true}; + for (bool hierarchy : {false, true}) { + config.vacuum_refinement_levels = hierarchy ? std::optional(1) : std::nullopt; + auto generated = stroid::GenerateMesh(config); + if (hierarchy) generated.mesh->EnsureNCMesh(); + const auto original = Serialize(*generated.mesh); + auto display = stroid::IO::MakeConformingVisualizationMesh(*generated.mesh); + EXPECT_EQ(display->GetNE(), generated.mesh->GetNE()); + EXPECT_EQ(Serialize(*display), original); + EXPECT_EQ(Serialize(*generated.mesh), original); + } +} diff --git a/tools/meson.build b/tools/meson.build index 33d5d0c..f98b6ff 100644 --- a/tools/meson.build +++ b/tools/meson.build @@ -1,4 +1,5 @@ -executable('stroid', 'stroid.cpp', dependencies: [stroid_dep, cli11_dep, magic_enum_dep], install: true) +stroid_cli = 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) + +executable('vacuum_refinement_experiment', 'vacuum_refinement_experiment.cpp', dependencies: [stroid_dep, cli11_dep], build_by_default: false, install: false) diff --git a/tools/stroid.cpp b/tools/stroid.cpp index 086b010..56d127a 100644 --- a/tools/stroid.cpp +++ b/tools/stroid.cpp @@ -7,6 +7,7 @@ #include #include #include +#include // ReSharper disable once CppUnusedIncludeDirective #include "mfem.hpp" @@ -72,12 +73,14 @@ int main(int argc, char** argv) { std::optional mesh_file; std::string output_filename = "stroid.mesh"; bool view_mesh = false; + bool original_elements = false; bool no_save = false; std::string glvis_host = "localhost"; int glvis_port = 19916; generate->add_option("-c,--config", config_filename, "Path to configuration file")->check(CLI::ExistingFile); generate->add_flag("-v,--view", view_mesh, "View the generated mesh using GLVis"); + generate->add_flag("--original-elements", original_elements, "Display original cells; curved hanging faces may show GLVis tessellation gaps"); generate->add_flag("-n,--no-save", no_save, "Do not save the generated mesh to a file"); generate->add_option("--glvis-host", glvis_host, "GLVis server host")->capture_default_str(); generate->add_option("--glvis-port", glvis_port, "GLVis server port")->capture_default_str(); @@ -85,6 +88,7 @@ int main(int argc, char** argv) { view->add_option("--host", glvis_host, "GLVis server host")->capture_default_str(); view->add_option("--port", glvis_port, "GLVis server port")->capture_default_str(); + view->add_flag("--original-elements", original_elements, "Display original cells; curved hanging faces may show GLVis tessellation gaps"); view->add_option("-f,--file", mesh_file, "Path to .mesh file")->check(CLI::ExistingFile); auto to_lower = [](std::string s) { @@ -167,7 +171,8 @@ int main(int argc, char** argv) { "Mesh Viewer - Colored by Element ID", selected_mode, glvis_host, - glvis_port); + glvis_port, + !original_elements); exit(0); } @@ -178,11 +183,8 @@ int main(int argc, char** argv) { } - const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); - stroid::topology::Finalize(*mesh, cfg); - stroid::topology::PromoteToHighOrder(*mesh, cfg); - stroid::topology::ProjectMesh(*mesh, cfg); - stroid::topology::OptimizeMesh(*mesh, cfg); + auto generated = stroid::GenerateMesh(cfg); + mfem::Mesh* mesh = generated.mesh.get(); if (!no_save) { const std::string& final_path = output_filename; @@ -193,6 +195,7 @@ int main(int argc, char** argv) { std::cerr << "WARNING! Saving to MFEM format without the standard '.mesh' extension. File will be called " << final_path << std::endl; } std::ofstream ofs(final_path); + ofs.precision(std::numeric_limits::max_digits10); mesh->Print(ofs); break; } @@ -201,12 +204,19 @@ int main(int argc, char** argv) { std::cerr << "WARNING! Saving to VTU format without the standard '.vtu' extension. File will be called " << final_path << std::endl; } std::ofstream ofs(final_path); + ofs.precision(std::numeric_limits::max_digits10); + // MFEM's stream overload writes an open Piece, allowing callers + // to append fields. Supply the enclosing document for a mesh export. + ofs << "\n\n"; mesh->PrintVTU(ofs, out_cfg.vtu.ref, out_cfg.vtu.format, out_cfg.vtu.high_order_output, out_cfg.vtu.compression_level, out_cfg.vtu.bdr_elements); + ofs << "\n\n\n"; break; } case MESH_FORMATS::VTK: { @@ -239,11 +249,12 @@ int main(int argc, char** argv) { "Spheroidal Mesh - Colored by Element ID", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, glvis_host, - glvis_port); + glvis_port, + !original_elements); } } else if (!*info) { std::println("Usage: {} [generate|info|view] --help", argv[0]); } return 0; -} \ No newline at end of file +} diff --git a/tools/vacuum_refinement_experiment.cpp b/tools/vacuum_refinement_experiment.cpp new file mode 100644 index 0000000..bd6dafe --- /dev/null +++ b/tools/vacuum_refinement_experiment.cpp @@ -0,0 +1,337 @@ +#include "stroid/stroid.h" +#include "CLI/CLI.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + constexpr double kPi = 3.14159265358979323846; + + class ScopedOutputRedirect { + std::streambuf* original; + public: + ScopedOutputRedirect() : original(std::cout.rdbuf(std::cerr.rdbuf())) {} + ~ScopedOutputRedirect() { std::cout.rdbuf(original); } + }; + + struct GeometryResults { + int stellar_elements = 0; + int vacuum_elements = 0; + int hanging_faces = 0; + size_t samples = 0; + size_t nonpositive_samples = 0; + double min_signed_det = std::numeric_limits::infinity(); + double max_condition = 0.0; + double max_face_mismatch = 0.0; + double max_coordinate_mismatch = 0.0; + double coordinate_constraint_residual = 0.0; + double stellar_volume = 0.0; + double stellar_volume_error = 0.0; + double surface_radius_error = 0.0; + double outer_radius_error = 0.0; + }; + + GeometryResults InspectGeometry(stroid::StroidMesh& generated, int grid_points) { + GeometryResults result; + auto& mesh = *generated.mesh; + const int vacuum = static_cast(generated.config.vacuum_id.value()); + const auto& coordinate = *generated.exterior_coordinate->values; + mfem::Vector independent; + coordinate.GetTrueDofs(independent); + mfem::GridFunction reconstructed(generated.exterior_coordinate->space.get()); + reconstructed.SetFromTrueDofs(independent); + reconstructed -= coordinate; + result.coordinate_constraint_residual = reconstructed.Normlinf(); + + for (int element = 0; element < mesh.GetNE(); ++element) { + const bool is_vacuum = mesh.GetAttribute(element) == vacuum; + is_vacuum ? ++result.vacuum_elements : ++result.stellar_elements; + auto* transformation = mesh.GetElementTransformation(element); + auto inspect = [&](const mfem::IntegrationPoint& point) { + transformation->SetIntPoint(&point); + const auto& jacobian = transformation->Jacobian(); + const double determinant = jacobian.Det(); + const double smallest = jacobian.CalcSingularvalue(2); + const double largest = jacobian.CalcSingularvalue(0); + ++result.samples; + if (!(determinant > 0.0) || !std::isfinite(determinant)) ++result.nonpositive_samples; + result.min_signed_det = std::min(result.min_signed_det, determinant); + result.max_condition = std::max(result.max_condition, + smallest > 0.0 ? largest / smallest : std::numeric_limits::infinity()); + }; + const auto& quadrature = mfem::IntRules.Get(mfem::Geometry::CUBE, 3 * transformation->Order() + 3); + for (int q = 0; q < quadrature.GetNPoints(); ++q) { + const auto& point = quadrature.IntPoint(q); + inspect(point); + if (!is_vacuum) result.stellar_volume += point.weight * transformation->Jacobian().Det(); + } + 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)); + inspect(point); + } + } + } + } + const double stellar_radius = generated.config.r_star.value(); + const double flattening = generated.config.flattening.value(); + const double analytic_volume = 4.0 * kPi * std::pow(stellar_radius, 3) * (1.0 - flattening) / 3.0; + result.stellar_volume_error = std::abs(result.stellar_volume - analytic_volume); + + mfem::Vector first(3), second(3); + for (int face = 0; face < mesh.GetNumFaces(); ++face) { + const auto information = mesh.GetFaceInformation(face); + if (!information.IsLocal()) continue; + if (information.IsNonconformingFine()) ++result.hanging_faces; + auto* transformation = mesh.GetFaceElementTransformations(face); + const bool stellar_interface = (transformation->Elem1->Attribute == vacuum) != + (transformation->Elem2->Attribute == vacuum); + if (stellar_interface && !information.IsConforming()) { + throw std::runtime_error("The stellar-vacuum interface is nonconforming."); + } + for (int i = 0; i < grid_points; ++i) { + for (int j = 0; j < grid_points; ++j) { + mfem::IntegrationPoint point; + point.Set2(static_cast(i) / (grid_points - 1), + static_cast(j) / (grid_points - 1)); + transformation->SetAllIntPoints(&point); + const auto& first_point = transformation->Elem1->GetIntPoint(); + const auto& second_point = transformation->Elem2->GetIntPoint(); + transformation->Elem1->Transform(first_point, first); + transformation->Elem2->Transform(second_point, second); + first -= second; + result.max_face_mismatch = std::max(result.max_face_mismatch, first.Norml2()); + const double difference = coordinate.GetValue(transformation->Elem1No, first_point) - + coordinate.GetValue(transformation->Elem2No, second_point); + result.max_coordinate_mismatch = std::max(result.max_coordinate_mismatch, std::abs(difference)); + } + } + } + + const auto& surface_quadrature = mfem::IntRules.Get(mfem::Geometry::SQUARE, 10); + for (int boundary = 0; boundary < mesh.GetNBE(); ++boundary) { + const bool surface = mesh.GetBdrAttribute(boundary) == static_cast(generated.config.surface_bdr_id.value()); + const bool outer = mesh.GetBdrAttribute(boundary) == static_cast(generated.config.inf_bdr_id.value()); + if (!surface && !outer) continue; + const double radius = surface ? stellar_radius : generated.config.r_infinity.value(); + auto* transformation = mesh.GetBdrElementTransformation(boundary); + for (int q = 0; q < surface_quadrature.GetNPoints(); ++q) { + transformation->Transform(surface_quadrature.IntPoint(q), first); + first(2) /= 1.0 - flattening; + double& error = surface ? result.surface_radius_error : result.outer_radius_error; + error = std::max(error, std::abs(first.Norml2() - radius)); + } + } + return result; + } + + struct SolveResults { + int dofs = 0; + int true_dofs = 0; + int iterations = 0; + bool converged = false; + double relative_residual = 0.0; + double l2_error = 0.0; + double stellar_l2_error = 0.0; + double vacuum_l2_error = 0.0; + double h1_error = 0.0; + }; + + SolveResults SolveManufacturedProblem(stroid::StroidMesh& generated, int solution_order) { + auto& mesh = *generated.mesh; + const double outer_sixth_power = std::pow(generated.config.r_infinity.value(), 6); + mfem::FunctionCoefficient exact([outer_sixth_power](const mfem::Vector& point) { + const double radius_squared = point * point; + return std::exp(-radius_squared) + 0.1 * std::pow(radius_squared, 3) / outer_sixth_power; + }); + mfem::VectorFunctionCoefficient gradient(3, [outer_sixth_power](const mfem::Vector& point, mfem::Vector& value) { + const double radius_squared = point * point; + value = point; + value *= -2.0 * std::exp(-radius_squared) + 0.6 * radius_squared * radius_squared / outer_sixth_power; + }); + mfem::FunctionCoefficient forcing([outer_sixth_power](const mfem::Vector& point) { + const double radius_squared = point * point; + return (6.0 - 4.0 * radius_squared) * std::exp(-radius_squared) + - 4.2 * radius_squared * radius_squared / outer_sixth_power; + }); + mfem::H1_FECollection collection(solution_order, 3); + mfem::FiniteElementSpace space(&mesh, &collection); + SolveResults result; + result.dofs = space.GetVSize(); + result.true_dofs = space.GetTrueVSize(); + mfem::Array boundary(mesh.bdr_attributes.Max()); + boundary = 0; + boundary[static_cast(generated.config.inf_bdr_id.value()) - 1] = 1; + mfem::Array essential; + space.GetEssentialTrueDofs(boundary, essential); + mfem::GridFunction solution(&space); + solution = 0.0; + solution.ProjectBdrCoefficient(exact, boundary); + const int quadrature_order = 2 * solution_order + 3 * generated.config.order.value() + 4; + const auto& quadrature = mfem::IntRules.Get(mfem::Geometry::CUBE, quadrature_order); + mfem::LinearForm rhs(&space); + auto* load = new mfem::DomainLFIntegrator(forcing); + load->SetIntRule(&quadrature); + rhs.AddDomainIntegrator(load); + rhs.Assemble(); + mfem::ConstantCoefficient one(1.0); + mfem::BilinearForm form(&space); + auto* diffusion = new mfem::DiffusionIntegrator(one); + diffusion->SetIntRule(&quadrature); + form.AddDomainIntegrator(diffusion); + form.Assemble(); + mfem::OperatorPtr system; + mfem::Vector independent, system_rhs; + form.FormLinearSystem(essential, solution, rhs, system, independent, system_rhs); + mfem::GSSmoother preconditioner(static_cast(*system)); + mfem::CGSolver solver; + solver.SetOperator(*system); + solver.SetPreconditioner(preconditioner); + solver.SetRelTol(1.0e-11); + solver.SetAbsTol(1.0e-14); + solver.SetMaxIter(2000); + solver.SetPrintLevel(-1); + solver.Mult(system_rhs, independent); + result.converged = solver.GetConverged(); + result.iterations = solver.GetNumIterations(); + mfem::Vector residual(system_rhs.Size()); + system->Mult(independent, residual); + residual -= system_rhs; + result.relative_residual = residual.Norml2() / system_rhs.Norml2(); + form.RecoverFEMSolution(independent, rhs, solution); + const mfem::IntegrationRule* rules[mfem::Geometry::NumGeom]{}; + rules[mfem::Geometry::CUBE] = &quadrature; + result.l2_error = solution.ComputeL2Error(exact, rules); + result.h1_error = solution.ComputeH1Error(&exact, &gradient, rules); + for (int element = 0; element < mesh.GetNE(); ++element) { + auto* transformation = mesh.GetElementTransformation(element); + double element_error = 0.0; + for (int q = 0; q < quadrature.GetNPoints(); ++q) { + const auto& point = quadrature.IntPoint(q); + transformation->SetIntPoint(&point); + const double difference = solution.GetValue(element, point) - exact.Eval(*transformation, point); + element_error += point.weight * transformation->Weight() * difference * difference; + } + if (mesh.GetAttribute(element) == static_cast(generated.config.vacuum_id.value())) { + result.vacuum_l2_error += element_error; + } else { + result.stellar_l2_error += element_error; + } + } + result.stellar_l2_error = std::sqrt(result.stellar_l2_error); + result.vacuum_l2_error = std::sqrt(result.vacuum_l2_error); + return result; + } +} + +int main(int argc, char* argv[]) { + std::vector orders{1, 2}; + std::vector refinements{2, 3}; + std::vector bulk_levels{0}; + int outer_level = -1; + int solution_order = 1; + int grid_points = 3; + double flattening = 0.0; + double infinity_radius = 6.0; + std::string output_path; + CLI::App app{"Compare uniform and graded vacuum meshes, geometry, and an H1 manufactured Poisson solution; TMOP is disabled."}; + app.add_option("--orders", orders, "Geometry orders, comma separated")->delimiter(',')->check(CLI::Range(1, 6)); + app.add_option("--refinements", refinements, "Stellar refinement levels, comma separated")->delimiter(',')->check(CLI::Range(0, 4)); + app.add_option("--bulk-levels", bulk_levels, "Vacuum background target levels, comma separated")->delimiter(',')->check(CLI::Range(0, 4)); + app.add_option("--outer-level", outer_level, "Outer vacuum target; -1 inherits stellar level")->check(CLI::Range(-1, 4)); + app.add_option("--solution-order", solution_order, "H1 polynomial order for the manufactured solve")->check(CLI::Range(1, 4)); + app.add_option("--grid-points", grid_points, "Closed tensor sample grid per coordinate, in addition to quadrature")->check(CLI::Range(2, 9)); + app.add_option("--flattening", flattening, "Spheroidal flattening")->check(CLI::Range(0.0, 0.9)); + app.add_option("--infinity-radius", infinity_radius, "Finite geometric outer radius; stellar radius is one"); + app.add_option("--output", output_path, "New CSV output file; default stdout"); + try { + app.parse(argc, argv); + } catch (const CLI::ParseError& error) { + return app.exit(error); + } + try { + if (!std::isfinite(infinity_radius) || infinity_radius <= 1.0) { + throw std::invalid_argument("Require a finite infinity-radius greater than one."); + } + 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 << "policy,geometry_order,stellar_level,bulk_target,outer_target,solution_order,flattening,r_infinity,elements,stellar_elements,vacuum_elements,vacuum_min_depth,vacuum_max_depth,geometry_dofs,geometry_true_dofs,solution_dofs,solution_true_dofs,hanging_faces,samples,nonpositive_samples,min_signed_det,max_condition,max_face_mismatch,max_coordinate_mismatch,coordinate_constraint_residual,stellar_volume,stellar_volume_error,surface_radius_error,outer_radius_error,l2_error,stellar_l2_error,vacuum_l2_error,h1_error,solver_converged,solver_iterations,relative_residual\n"; + std::cerr << "Evaluating actual FE geometry using signed Jacobians at quadrature points and a closed grid.\n" + "Manufactured problem: -Delta u=f, u=exp(-r^2)+0.1*r^6/R^6, exact Dirichlet data only at the outer boundary.\n" + "This finite-domain Poisson diagnostic measures discretization error; it does not implement compactified physics.\n"; + bool verified = true; + for (const int order : orders) { + for (const int refinement : refinements) { + std::vector cases{-1}; + cases.insert(cases.end(), bulk_levels.begin(), bulk_levels.end()); + for (const int bulk : cases) { + const bool uniform = bulk < 0; + const std::string policy = uniform ? "uniform" : "graded"; + stroid::config::MeshConfig config; + config.order = order; + config.refinement_levels = refinement; + config.flattening = flattening; + config.r_infinity = infinity_radius; + config.optimization_methods = stroid::config::OptimizationMethods{false, true}; + if (!uniform) { + config.vacuum_refinement_levels = bulk; + if (outer_level >= 0) config.vacuum_outer_refinement_levels = outer_level; + } + std::cerr << "Inspecting " << policy << ", geometry order " << order << ", stellar level " << refinement + << ", bulk " << (uniform ? refinement : bulk) << '\n'; + stroid::StroidMesh generated; + { + ScopedOutputRedirect redirect; + generated = stroid::GenerateMesh(config); + } + const auto geometry = InspectGeometry(generated, grid_points); + const auto stats = stroid::stats::ComputeMeshStats(generated, stroid::stats::MeshStatFeatures::REFINEMENT); + const auto solve = SolveManufacturedProblem(generated, solution_order); + verified = verified && geometry.nonpositive_samples == 0 && geometry.max_face_mismatch < 1.0e-10 + && geometry.max_coordinate_mismatch < 1.0e-11 && geometry.coordinate_constraint_residual < 1.0e-11 + && solve.converged && solve.relative_residual < 1.0e-9 && std::isfinite(solve.h1_error); + const int outer = uniform ? refinement : outer_level < 0 ? refinement : outer_level; + output << policy << ',' << order << ',' << refinement << ',' << (uniform ? refinement : bulk) << ',' << outer + << ',' << solution_order << ',' << flattening << ',' << infinity_radius << ',' << generated.mesh->GetNE() + << ',' << geometry.stellar_elements << ',' << geometry.vacuum_elements + << ',' << stats.refinement->vacuum.min_depth << ',' << stats.refinement->vacuum.max_depth + << ',' << stats.refinement->geometry_dofs << ',' << stats.refinement->geometry_true_dofs + << ',' << solve.dofs << ',' << solve.true_dofs << ',' << geometry.hanging_faces + << ',' << geometry.samples << ',' << geometry.nonpositive_samples << ',' << geometry.min_signed_det + << ',' << geometry.max_condition << ',' << geometry.max_face_mismatch << ',' << geometry.max_coordinate_mismatch + << ',' << geometry.coordinate_constraint_residual << ',' << geometry.stellar_volume << ',' << geometry.stellar_volume_error + << ',' << geometry.surface_radius_error << ',' << geometry.outer_radius_error << ',' << solve.l2_error + << ',' << solve.stellar_l2_error << ',' << solve.vacuum_l2_error << ',' << solve.h1_error + << ',' << solve.converged << ',' << solve.iterations << ',' << solve.relative_residual << '\n'; + output.flush(); + if (!output) throw std::runtime_error("Failed to write experiment output."); + } + } + } + if (!verified) { + std::cerr << "Numerical verification failed: inspect Jacobian, trace, or solver columns.\n"; + return 1; + } + } catch (const std::exception& error) { + std::cerr << "Vacuum refinement experiment failed: " << error.what() << '\n'; + return 1; + } + return 0; +}